00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00085
00086
00087
00088
00089
00092 #include "StPointCollection.h"
00093 #include "StPi0Candidate.h"
00094 #include "emc_def.h"
00095 #include <Stiostream.h>
00096 #include <TBrowser.h>
00097 #include "StEmcUtil/geometry/StEmcGeom.h"
00098 #include "StEmcUtil/projection/StEmcPosition.h"
00099 #include "StEvent.h"
00100 #include "StEventTypes.h"
00101 #include "StarCallf77.h"
00102 #include "TMath.h"
00103
00104
00105 #define assndx F77_NAME(assndx,ASSNDX)
00106 extern "C"
00107 {
00108 void type_of_call assndx ( Int_t &, Float_t *, Int_t &, Int_t &,
00109 Int_t &, Int_t *, Float_t &,Int_t*,Int_t &);
00110 }
00111
00112 ClassImp(StPointCollection)
00113
00114 const TString detname[] =
00115 {"Bemc", "Bsmde", "Bsmdp"
00116 };
00117
00118 StMatchVecClus matchlist_bemc_clus[Epc::nModule][Epc::nPhiBin];
00119 StMatchVecClus matchlist_bprs_clus[Epc::nModule][Epc::nPhiBin];
00120 StMatchVecClus matchlist_bsmde_clus[Epc::nModule][Epc::nPhiBin];
00121 StMatchVecClus matchlist_bsmdp_clus[Epc::nModule][Epc::nPhiBin];
00122
00123 FloatVector HitTrackEta;
00124 FloatVector HitTrackPhi;
00125 FloatVector HitTrackMom;
00126
00127
00128 StPointCollection::StPointCollection():TDataSet("Default")
00129 {
00130 SetTitle("EmcPoints");
00131 mPrint = kTRUE;
00132 mBField = 0.5;
00133 mNPoints =0;
00134 mNPointsReal=0;
00135 mPosition = new StEmcPosition();
00136 }
00137
00138 StPointCollection::StPointCollection(const Char_t *Name):TDataSet(Name)
00139 {
00140 SetTitle("EmcPoints");
00141 mPrint = kTRUE;
00142 mBField = 0.5;
00143 mNPoints =0;
00144 mNPointsReal=0;
00145 mPosition = new StEmcPosition();
00146 }
00147
00148 StPointCollection::~StPointCollection()
00149 {
00150 mPoints.Delete();
00151
00152 mNPoints =0;
00153 mNPointsReal=0;
00154 delete mPosition;
00155 }
00156
00157 void
00158 StPointCollection::Browse(TBrowser* b)
00159 {
00160 TDataSet::Browse(b);
00161 }
00162
00163 Int_t StPointCollection::makeEmcPoints(StEvent* event)
00164 {
00165 LOG_DEBUG <<"Finding EMC Points ..." << endm;
00166 if(!event)
00167 return 0;
00168
00169 StEmcCollection* emc = event->emcCollection();
00170 if(!emc)
00171 return 0;
00172
00173 StSPtrVecEmcPoint& points = emc->barrelPoints();
00174 if(points.size())
00175 points.clear();
00176
00177 StEmcClusterCollection* cluster[4] = {0,0,0,0};
00178
00179 for(Int_t i = 0;i<4;i++)
00180 {
00181 StDetectorId EmcId = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00182 ;
00183 StEmcDetector* EmcDet = emc->detector(EmcId);
00184 if(EmcDet)
00185 cluster[i] = EmcDet->cluster();
00186 }
00187
00188
00189
00190 findMatchedClusters(cluster[0],cluster[1],cluster[2],cluster[3]);
00191
00192
00193 LOG_DEBUG <<"Making points for non-matched clusters ..." <<endm;
00194 ClusterSort(cluster[0],cluster[1],cluster[2],cluster[3]);
00195
00196
00197 for(Int_t im=0;im<Epc::nModule;im++)
00198 {
00199 for(Int_t is=0;is<Epc::nPhiBin;is++)
00200 {
00201 if(matchlist_bemc_clus[im][is].size()>0)
00202 {
00203 matchClusters(matchlist_bemc_clus[im][is],
00204 matchlist_bprs_clus[im][is],
00205 matchlist_bsmde_clus[im][is],
00206 matchlist_bsmdp_clus[im][is]);
00207 }
00208 }
00209 }
00210
00211
00212 matchToTracks(event);
00213 return 1;
00214 }
00215
00216 Int_t StPointCollection::findMatchedClusters(StEmcClusterCollection* Bemccluster,
00217 StEmcClusterCollection *Bprscluster,
00218 StEmcClusterCollection *Bsmdecluster,
00219 StEmcClusterCollection *Bsmdpcluster)
00220 {
00221 LOG_DEBUG <<"Making points for already matched clusters ..." <<endm;
00222 if(Bemccluster)
00223 {
00224 Int_t Ncluster0=Bemccluster->numberOfClusters();
00225 if(Ncluster0>0)
00226 {
00227 const StSPtrVecEmcCluster& emcclusters= Bemccluster->clusters();
00228 for(UInt_t i=0;i<emcclusters.size();i++)
00229 {
00230 StEmcCluster *btow=(StEmcCluster*)emcclusters[i];
00231 UInt_t matchId = btow->GetUniqueID();
00232 if(matchId!=0)
00233 {
00234 StEmcCluster *bprs = NULL;
00235 StEmcCluster *bsmde = NULL;
00236 StEmcCluster *bsmdp = NULL;
00237 if(Bprscluster)
00238 {
00239 const StSPtrVecEmcCluster& clusters= Bprscluster->clusters();
00240 for(UInt_t i=0;i<clusters.size();i++)
00241 if( ((StEmcCluster*)clusters[i])->GetUniqueID() == matchId)
00242 {
00243 bprs = (StEmcCluster*)clusters[i];
00244 break;
00245 }
00246 }
00247 if(Bsmdecluster)
00248 {
00249 const StSPtrVecEmcCluster& clusters= Bsmdecluster->clusters();
00250 for(UInt_t i=0;i<clusters.size();i++)
00251 if( ((StEmcCluster*)clusters[i])->GetUniqueID() == matchId)
00252 {
00253 bsmde = (StEmcCluster*)clusters[i];
00254 break;
00255 }
00256 }
00257 if(Bsmdpcluster)
00258 {
00259 const StSPtrVecEmcCluster& clusters= Bsmdpcluster->clusters();
00260 for(UInt_t i=0;i<clusters.size();i++)
00261 if( ((StEmcCluster*)clusters[i])->GetUniqueID() == matchId)
00262 {
00263 bsmdp = (StEmcCluster*)clusters[i];
00264 break;
00265 }
00266 }
00267
00268 makePoint(btow,bprs,bsmde,bsmdp);
00269 }
00270 }
00271 }
00272 }
00273 return 0;
00274 }
00275 StEmcPoint* StPointCollection::makePoint(StEmcCluster* btow,StEmcCluster* bprs,StEmcCluster* bsmde,StEmcCluster* bsmdp,Float_t fraction)
00276 {
00277 if(!btow)
00278 return NULL;
00279
00280 LOG_DEBUG <<"Making point"<< endm;
00281
00282 Int_t Category = 0;
00283 if(btow)
00284 Category = Category | 1;
00285 if(bprs)
00286 Category = Category | 2;
00287 if(bsmde)
00288 Category = Category | 4;
00289 if(bsmdp)
00290 Category = Category | 8;
00291
00292 Float_t en[4] = {0,0,0,0};
00293 Float_t si[4] = {0,0,0,0};
00294
00295 Float_t eta = btow->eta();
00296 Float_t phi = btow->phi();
00297
00298
00299
00300
00301 Float_t energy = btow->energy()*fraction;
00302 if(fraction<0)
00303 energy = fabs(fraction);
00304 Float_t sigEta = btow->sigmaEta();
00305 Float_t sigPhi = btow->sigmaPhi();
00306 en[0] = btow->energy();
00307 si[0] = sqrt(eta*eta+phi*phi);
00308
00309 if(bprs)
00310 {
00311 en[1] = bprs->energy();
00312 si[1] = sqrt(bprs->sigmaEta()*bprs->sigmaEta()+bprs->sigmaPhi()*bprs->sigmaPhi());
00313 }
00314
00315 if(bsmde)
00316 {
00317 eta = bsmde->eta();
00318 sigEta = bsmde->sigmaEta();
00319 en[2] = bsmde->energy();
00320 si[2] = sqrt(bsmde->sigmaEta()*bsmde->sigmaEta()+bsmde->sigmaPhi()*bsmde->sigmaPhi());
00321 }
00322
00323 if(bsmdp)
00324 {
00325 phi = bsmdp->phi();
00326 sigPhi = bsmdp->sigmaPhi();
00327 en[3] = bsmdp->energy();
00328 si[3] = sqrt(bsmdp->sigmaEta()*bsmdp->sigmaEta()+bsmdp->sigmaPhi()*bsmdp->sigmaPhi());
00329 }
00330
00331 Float_t xp,yp,zp;
00332
00333
00334 xp=(StEpcCut::RAD_SMD_E())*cos(phi);
00335 yp=(StEpcCut::RAD_SMD_E())*sin(phi);
00336 zp=(StEpcCut::RAD_SMD_E())*sinh(eta);
00337 StThreeVectorF PointPosition(xp*centimeter, yp*centimeter, zp*centimeter);
00338
00339
00340 xp=(StEpcCut::RAD_SMD_E())*(cos(phi+sigPhi)-cos(phi-sigPhi))/2;
00341 yp=(StEpcCut::RAD_SMD_E())*(sin(phi+sigPhi)-sin(phi-sigPhi))/2;
00342 zp=(StEpcCut::RAD_SMD_E())*(sinh(eta+sigEta)-sinh(eta-sigEta))/2;
00343 StThreeVectorF ErrorPosition(xp*centimeter, yp*centimeter, zp*centimeter);
00344
00345 StThreeVectorF size(sigEta,sigPhi,0.0);
00346
00347
00348 Float_t ChiSquare = 0.0;
00349
00350
00351 StEmcPoint *point = new StEmcPoint();
00352 point->setQuality(Category);
00353 point->setPosition(PointPosition);
00354 point->setPositionError(ErrorPosition);
00355 point->setSize(size);
00356 point->setChiSquare(ChiSquare);
00357 point->setEnergy(energy);
00358 point->setDeltaEta(9999);
00359 point->setDeltaPhi(9999);
00360 for(Int_t i=0;i<4;i++)
00361 {
00362 StDetectorId id=static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00363 point->setEnergyInDetector(id,en[i]);
00364 point->setSizeAtDetector(id,si[i]);
00365 if(i==0 && btow)
00366 point->addCluster(id,btow);
00367 else
00368 point->addCluster(id,NULL);
00369 if(i==1 && bprs)
00370 point->addCluster(id,bprs);
00371 else
00372 point->addCluster(id,NULL);
00373 if(i==2 && bsmde)
00374 point->addCluster(id,bsmde);
00375 else
00376 point->addCluster(id,NULL);
00377 if(i==3 && bsmdp)
00378 point->addCluster(id,bsmdp);
00379 else
00380 point->addCluster(id,NULL);
00381 }
00382
00383 mPointsReal.Add(point);
00384 mNPointsReal++;
00385
00386 return point;
00387
00388 }
00389
00390 Int_t StPointCollection::matchClusters(const StMatchVecClus mvec,
00391 const StMatchVecClus prsvec,
00392 const StMatchVecClus evec,
00393 const StMatchVecClus pvec)
00394
00395 {
00396 Int_t na=evec.size();
00397 Int_t ma=pvec.size();
00398 Float_t smin;
00399 Int_t mode=1;
00400 Int_t Category=-1;
00401 Float_t EmcTot = 0;
00402 Float_t totAvg = 0;
00403 Int_t idw =Epc::nMaxNoOfClinBin;
00404 Int_t ida =Epc::nMaxNoOfClinBin;
00405 Float_t ep[Epc::nMaxNoOfClinBin][Epc::nMaxNoOfClinBin];
00406 Int_t iw[Epc::nMaxNoOfClinBin][Epc::nMaxNoOfClinBin];
00407 Int_t k[Epc::nMaxNoOfClinBin];
00408
00409
00410 for (Int_t iF=0;iF<Epc::nMaxNoOfClinBin;iF++)
00411 {
00412 k[iF]=0;
00413 for (Int_t iL=0;iL<Epc::nMaxNoOfClinBin;iL++)
00414 {
00415 ep[iF][iL]=0.0;
00416 iw[iF][iL]=0;
00417 }
00418 }
00419 if(evec.size()==0 && pvec.size()==0)
00420 {
00421 Category=0;
00422 na=mvec.size();
00423 ma=mvec.size();
00424 }
00425 if(evec.size()>0 && pvec.size()==0)
00426 {
00427 Category=1;
00428 na=mvec.size();
00429 ma=evec.size();
00430 }
00431 if(evec.size()==0 && pvec.size()>0)
00432 {
00433 Category=2;
00434 na=mvec.size();
00435 ma=pvec.size();
00436 }
00437 if(evec.size()>0 && pvec.size()>0)
00438 {
00439 Category=3;
00440 na=evec.size();
00441 ma=pvec.size();
00442 }
00443
00444
00445 if(mvec.size()>0)
00446 {
00447 for (UInt_t ims=0;ims<mvec.size();ims++)
00448 {
00449 StEmcCluster *cl0;
00450 cl0=(StEmcCluster*)mvec[ims];
00451 Float_t emen=cl0->energy();
00452 EmcTot+=emen;
00453 }
00454 }
00455
00456
00457 for(Int_t ie=0;ie<na;ie++)
00458 {
00459 StEmcCluster *cl1 = NULL;
00460 for(Int_t ip=0;ip<ma;ip++)
00461 {
00462 StEmcCluster *cl2 = NULL;
00463 switch (Category)
00464 {
00465 case 0:
00466 cl1 = (StEmcCluster*)mvec[ie];
00467 cl2 = (StEmcCluster*)mvec[ip];
00468 break;
00469 case 1:
00470 cl1 = (StEmcCluster*)mvec[ie];
00471 cl2 = (StEmcCluster*)evec[ip];
00472 break;
00473 case 2:
00474 cl1 = (StEmcCluster*)mvec[ie];
00475 cl2 = (StEmcCluster*)pvec[ip];
00476 break;
00477 case 3:
00478 cl1 = (StEmcCluster*)evec[ie];
00479 cl2 = (StEmcCluster*)pvec[ip];
00480 break;
00481 }
00482
00483 if(cl1 && cl2){
00484 Float_t diff=TMath::Abs((cl1->energy())-(cl2->energy()));
00485 Float_t summ= (cl1->energy())+(cl2->energy());
00486 ep[ip][ie]=diff/summ;
00487 }
00488 }
00489 }
00490
00491 assndx(mode,ep[0],na,ma,ida,k,smin,iw[0],idw);
00492
00493 int i1;
00494 switch (Category)
00495 {
00496 case 0:
00497 for(i1=0;i1<na;i1++)
00498 {
00499 if((k[i1]-1)>=0)
00500 {
00501 StEmcCluster *cl1;
00502 cl1 = (StEmcCluster*)mvec[i1];
00503 Float_t avg_en = cl1->energy();
00504 totAvg += avg_en;
00505 }
00506 }
00507 break;
00508 case 1:
00509 for(i1=0;i1<na;i1++)
00510 {
00511 if((k[i1]-1)>=0)
00512 {
00513 StEmcCluster *cl1;
00514 cl1 = (StEmcCluster*)mvec[i1];
00515 Float_t avg_en = cl1->energy();
00516 totAvg += avg_en;
00517 }
00518 }
00519 break;
00520 case 2:
00521 for(i1=0;i1<na;i1++)
00522 {
00523 if((k[i1]-1)>=0)
00524 {
00525 StEmcCluster *cl1;
00526 cl1 = (StEmcCluster*)mvec[i1];
00527 Float_t avg_en = cl1->energy();
00528 totAvg += avg_en;
00529 }
00530 }
00531 break;
00532 case 3:
00533 for(i1=0;i1<na;i1++)
00534 {
00535 if((k[i1]-1)>=0)
00536 {
00537 StEmcCluster *cl1;
00538 cl1 = (StEmcCluster*)evec[i1];
00539 StEmcCluster *cl2;
00540 cl2 = (StEmcCluster*)pvec[k[i1]-1];
00541 Float_t avg_en = 0.5*(cl1->energy()+cl2->energy());
00542 totAvg += avg_en;
00543 }
00544 }
00545 break;
00546 }
00547
00548 for(i1=0;i1<na;i1++)
00549 {
00550
00551 StEmcCluster *btow = NULL;
00552 StEmcCluster *bprs = NULL;
00553 StEmcCluster *bsmde = NULL;
00554 StEmcCluster *bsmdp = NULL;
00555 Float_t fraction = 1;
00556 Float_t eta, phi;
00557 eta = phi = -999.;
00558
00559 switch (Category)
00560 {
00561 case 0:
00562 if((k[i1]-1)>=0)
00563 {
00564 btow = (StEmcCluster*)mvec[i1];
00565 eta = btow->eta();
00566 phi = btow->phi();
00567 fraction = 1;
00568 }
00569 break;
00570 case 1:
00571 if((k[i1]-1)>=0)
00572 {
00573 btow = (StEmcCluster*)mvec[i1];
00574 bsmde = (StEmcCluster*)evec[k[i1]-1];
00575 eta = bsmde->eta();
00576 phi = btow->phi();
00577 fraction = EmcTot/totAvg;
00578 }
00579 break;
00580 case 2:
00581 if((k[i1]-1)>=0)
00582 {
00583 btow = (StEmcCluster*)mvec[i1];
00584 bsmdp = (StEmcCluster*)pvec[k[i1]-1];
00585 eta=btow->eta();
00586 phi=bsmdp->phi();
00587 fraction = EmcTot/totAvg;
00588 }
00589 break;
00590 case 3:
00591 if((k[i1]-1)>=0)
00592 {
00593 bsmde = (StEmcCluster*)evec[i1];
00594 bsmdp = (StEmcCluster*)pvec[k[i1]-1];
00595 fraction = -fabs(EmcTot*0.5*(bsmde->energy()+bsmdp->energy())/totAvg);
00596 eta = bsmde->eta();
00597 phi = bsmdp->phi();
00598 Float_t delta = 999999;
00599 for (UInt_t ims=0;ims<mvec.size();ims++)
00600 {
00601 StEmcCluster *cl0 = (StEmcCluster*)mvec[ims];
00602 Float_t de = sqrt((eta-cl0->eta())*(eta-cl0->eta()) +
00603 (phi-cl0->phi())*(phi-cl0->phi()));
00604 if(de<delta)
00605 {
00606 btow = cl0;
00607 delta = de;
00608 }
00609 }
00610 }
00611 break;
00612 }
00613
00614 if(prsvec.size()>0)
00615 {
00616 Float_t delta = 999999;
00617 for (UInt_t ims=0;ims<prsvec.size();ims++)
00618 {
00619 StEmcCluster *cl0 = (StEmcCluster*)prsvec[ims];
00620 Float_t de = sqrt((eta-cl0->eta())*(eta-cl0->eta()) +
00621 (phi-cl0->phi())*(phi-cl0->phi()));
00622 if(de<delta)
00623 {
00624 bprs = cl0;
00625 delta = de;
00626 }
00627 }
00628 }
00629 makePoint(btow,bprs,bsmde,bsmdp,fraction);
00630 }
00631 return 1;
00632 }
00633
00634 void StPointCollection::ClusterSort(StEmcClusterCollection* Bemccluster,
00635 StEmcClusterCollection* Bprscluster,
00636 StEmcClusterCollection* Bsmdecluster,
00637 StEmcClusterCollection* Bsmdpcluster)
00638 {
00639
00640
00641
00642
00643
00644
00645 const Int_t eta_shift_fix=1;
00646 LOG_DEBUG <<" I am inside PointCalc***"<<endm;
00647 for(Int_t i1=0;i1<Epc::nModule;i1++)
00648 {
00649 for(Int_t i2=0;i2<Epc::nPhiBin;i2++)
00650 {
00651 matchlist_bemc_clus[i1][i2].clear();
00652 matchlist_bprs_clus[i1][i2].clear();
00653 matchlist_bsmde_clus[i1][i2].clear();
00654 matchlist_bsmdp_clus[i1][i2].clear();
00655 }
00656 }
00657
00658 StEmcGeom* GeomIn = StEmcGeom::getEmcGeom("bemc");
00659
00660
00661 if(Bemccluster)
00662 {
00663 Int_t Ncluster0=Bemccluster->numberOfClusters();
00664 if(Ncluster0>0)
00665 {
00666 const StSPtrVecEmcCluster& emcclusters= Bemccluster->clusters();
00667 for(UInt_t i=0;i<emcclusters.size();i++)
00668 {
00669 StEmcCluster *cl1=(StEmcCluster*)emcclusters[i];
00670 LOG_DEBUG <<"BEMC cluster UniqueId = "<<cl1->GetUniqueID()<<endm;
00671 if(cl1->GetUniqueID()==0)
00672 {
00673 Float_t eta_emc=cl1->eta();
00674 Float_t phi_emc=cl1->phi();
00675
00676 Int_t ebin,pbin;
00677 Int_t emc_module=0;
00678 Int_t & imd =emc_module;
00679 Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
00680 if(testb==0)
00681 emc_module=imd;
00682 if(testb==0)
00683 {
00684 Int_t emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
00685
00686
00687 if(!eta_shift_fix && emc_phi_bin>0)
00688 {
00689 if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<0.2)
00690 {
00691 emc_phi_bin--;
00692 }
00693 }
00694 if(emc_phi_bin>9)
00695 {
00696 emc_phi_bin=9;
00697 }
00698
00699 matchlist_bemc_clus[emc_module-1][emc_phi_bin].push_back(cl1);
00700 }
00701 }
00702 }
00703 }
00704 }
00705
00706
00707 if(Bprscluster)
00708 {
00709 Int_t Ncluster1=Bprscluster->numberOfClusters();
00710 if(Ncluster1>0)
00711 {
00712 const StSPtrVecEmcCluster& emcclusters= Bprscluster->clusters();
00713 for(UInt_t i=0;i<emcclusters.size();i++)
00714 {
00715 StEmcCluster *cl2=(StEmcCluster*)emcclusters[i];
00716 if(cl2->GetUniqueID()==0)
00717 {
00718 Float_t eta_emc=cl2->eta();
00719 Float_t phi_emc=cl2->phi();
00720
00721 Int_t ebin,pbin;
00722 Int_t emc_module=0;
00723 Int_t emc_phi_bin=0;
00724 Int_t & imd =emc_module;
00725 Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
00726 if(testb==0)
00727 emc_module=imd;
00728 if(testb==0)
00729 {
00730 emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
00731
00732 if(!eta_shift_fix && emc_phi_bin>0)
00733 {
00734 if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<0.2)
00735 {
00736 emc_phi_bin--;
00737 }
00738 }
00739 if(emc_phi_bin>9)
00740 {
00741 emc_phi_bin=9;
00742 }
00743
00744 matchlist_bprs_clus[emc_module-1][emc_phi_bin].push_back(cl2);
00745 }
00746 }
00747 }
00748 }
00749 }
00750
00751
00752 if(Bsmdecluster)
00753 {
00754 Int_t Ncluster2=Bsmdecluster->numberOfClusters();
00755 if(Ncluster2>0)
00756 {
00757 const StSPtrVecEmcCluster& emcclusters= Bsmdecluster->clusters();
00758 for(UInt_t i=0;i<emcclusters.size();i++)
00759 {
00760 StEmcCluster *cl3=(StEmcCluster*)emcclusters[i];
00761 if(cl3->GetUniqueID()==0)
00762 {
00763 Float_t eta_emc=cl3->eta();
00764 Float_t phi_emc=cl3->phi();
00765
00766 Int_t ebin,pbin;
00767 Int_t emc_module=0;
00768 Int_t emc_phi_bin=0;
00769 Int_t & imd =emc_module;
00770 Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
00771 if(testb==0)
00772 emc_module=imd;
00773 if(testb==0)
00774 {
00775 emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
00776
00777 if(!eta_shift_fix && emc_phi_bin>0)
00778 {
00779 if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<.2)
00780 {
00781 emc_phi_bin--;
00782 }
00783 }
00784 if(emc_phi_bin>9)
00785 {
00786 emc_phi_bin=9;
00787 }
00788
00789 matchlist_bsmde_clus[emc_module-1][emc_phi_bin].push_back(cl3);
00790 }
00791 }
00792 }
00793 }
00794 }
00795
00796
00797 if(Bsmdpcluster)
00798 {
00799 Int_t Ncluster3=Bsmdpcluster->numberOfClusters();
00800 if(Ncluster3>0)
00801 {
00802 const StSPtrVecEmcCluster& emcclusters= Bsmdpcluster->clusters();
00803 for(UInt_t i=0;i<emcclusters.size();i++)
00804 {
00805 StEmcCluster *cl4=(StEmcCluster*)emcclusters[i];
00806 if(cl4->GetUniqueID()==0)
00807 {
00808 Float_t eta_emc=cl4->eta();
00809 Float_t phi_emc=cl4->phi();
00810
00811 Int_t ebin,pbin;
00812 Int_t emc_module=0;
00813 Int_t & imd =emc_module;
00814 Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
00815 if(testb==0)
00816 emc_module=imd;
00817 if(testb==0)
00818 {
00819 Int_t emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
00820
00821 if(!eta_shift_fix && emc_phi_bin>0)
00822 {
00823 if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<=0.01)
00824 {
00825 emc_phi_bin--;
00826 }
00827 }
00828 if(emc_phi_bin>9)
00829 {
00830 emc_phi_bin=9;
00831 }
00832
00833 matchlist_bsmdp_clus[emc_module-1][emc_phi_bin].push_back(cl4);
00834 }
00835 }
00836 }
00837 }
00838 }
00839 }
00840
00841 Int_t StPointCollection::matchToTracks(StEvent* event)
00842 {
00843 if(!event)
00844 return 0;
00845
00846 Float_t field = 0.5;
00847 StEventSummary *evtSummary = event->summary();
00848 if (evtSummary)
00849 field = evtSummary->magneticField()/10;
00850
00851 Int_t nR = NPointsReal();
00852 StEmcGeom* geom = StEmcGeom::instance("bemc");
00853 if(nR>0)
00854 {
00855 LOG_DEBUG << "Matching to tracks... NP = " << nR << endm;
00856 StSPtrVecTrackNode& tracks=event->trackNodes();
00857 Int_t nTracks = tracks.size();
00858 StThreeVectorD momentum,position;
00859 for(Int_t t=0;t<nTracks;t++)
00860 {
00861 StTrack *track = tracks[t]->track(0);
00862 if(track)
00863 {
00864 if(track->geometry())
00865 {
00866 Bool_t tok = mPosition->trackOnEmc(&position,&momentum,
00867 track,(Double_t)field,
00868 (Double_t)geom->Radius());
00869 if(tok)
00870 {
00871 Float_t eta = position.pseudoRapidity();
00872 Float_t phi = position.phi();
00873 if(fabs(eta)<1)
00874 {
00875 TIter next(PointsReal());
00876 StEmcPoint *cl;
00877
00878 for(Int_t i=0; i<nR; i++)
00879 {
00880 cl = (StEmcPoint*)next();
00881 if(cl)
00882 {
00883 StThreeVectorF pos = cl->position();
00884 Float_t etaP = pos.pseudoRapidity();
00885 Float_t phiP = pos.phi();
00886 Float_t D = sqrt(cl->deltaEta()*cl->deltaEta()+cl->deltaPhi()*cl->deltaPhi());
00887 Float_t d = sqrt((eta-etaP)*(eta-etaP)+(phi-phiP)*(phi-phiP));
00888 if(d<D)
00889 {
00890 cl->setDeltaEta(eta-etaP);
00891 cl->setDeltaPhi(phi-phiP);
00892 }
00893
00894 StThreeVectorF err = cl->positionError();
00895 Float_t etaE = err.pseudoRapidity();
00896 Float_t phiE = err.phi();
00897
00898 Float_t dPhi = fabs(phi-phiP);
00899 if (dPhi>TMath::Pi())
00900 dPhi=2*TMath::Pi()-dPhi;
00901
00902 if(fabs(eta-etaP)<fabs(etaE) && dPhi<fabs(phiE))
00903 {
00904 Int_t Category = cl->quality();
00905 Category = Category | 16;
00906 cl->setQuality(Category);
00907 cl->addTrack(track);
00908 }
00909 }
00910 }
00911 }
00912 }
00913 }
00914 }
00915 }
00916 }
00917
00918 return 0;
00919 }
