StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StPointCollection.cxx
1 //
2 // $id$
3 //
4 // $Log: StPointCollection.cxx,v $
5 // Revision 1.28 2007/07/25 16:53:20 kocolosk
6 // bugfix for the bugfix in 1.27
7 //
8 // Revision 1.27 2007/07/24 15:41:44 kocolosk
9 // bugFix from Oleksandr:
10 // http://www.star.bnl.gov/HyperNews-star/get/emc2/2444.html
11 //
12 // Revision 1.26 2007/01/22 19:13:50 kocolosk
13 // use STAR logger for all output
14 //
15 // Revision 1.25 2006/09/20 13:44:28 kocolosk
16 // fix autobuild warnings
17 //
18 // Revision 1.24 2006/03/24 19:31:15 suaide
19 // fixed break segmentation error.
20 //
21 // Revision 1.23 2005/05/23 12:35:14 suaide
22 // New Point maker code
23 //
24 // Revision 1.22 2004/08/13 13:08:01 suaide
25 // small fixed introduced by Marco to remove ineficiencies on the
26 // edges of phi bins
27 //
28 // Revision 1.20 2003/10/21 15:35:25 suaide
29 // fix a break segmentation introduced when a memory leak was fixed
30 //
31 // Revision 1.19 2003/10/12 02:56:51 perev
32 // LeakOff TClonesArray::Delete added
33 //
34 // Revision 1.18 2003/09/17 00:55:59 suaide
35 // Fixed bug. Chain was crashing because some data members were not initialized
36 //
37 // Revision 1.17 2003/09/02 17:58:03 perev
38 // gcc 3.2 updates + WarnOff
39 //
40 // Revision 1.16 2003/05/26 13:44:34 suaide
41 // added setPrint() method
42 //
43 // Revision 1.15 2003/01/23 04:03:21 jeromel
44 // Include fixed
45 //
46 // Revision 1.14 2001/12/03 22:24:28 pavlinov
47 // tuned for case of no tracks
48 //
49 // Revision 1.13 2001/12/01 02:44:50 pavlinov
50 // Cleanp for events with zero number of tracks
51 //
52 // Revision 1.12 2001/11/06 23:35:27 suaide
53 // fixed bug in the way we get magnetic field
54 //
55 // Revision 1.11 2001/10/24 13:55:05 suaide
56 // small bugs fixed
57 //
58 // Revision 1.10 2001/09/29 01:15:17 pavlinov
59 // Clean up for production
60 //
61 // Revision 1.9 2001/09/24 15:14:55 pavlinov
62 // No public constructor for StEmcGeom
63 //
64 // Revision 1.8 2001/08/18 22:13:49 subhasis
65 // phi-cluster attached to cat#3 points corrected
66 //
67 // Revision 1.7 2001/04/25 17:27:44 perev
68 // HPcorrs
69 //
70 // Revision 1.6 2001/04/24 23:06:29 subhasis
71 // clusters attached to Points, QA hists are made for all category separately
72 //
73 // Revision 1.3 2000/08/29 20:33:04 subhasis
74 // Modified to accept input from StEvent and writing output to StEvent for Emc
75 //
76 // Revision 1.1 2000/05/15 21:18:33 subhasis
77 // initial version
78 //
79 // Pi0 Candidate Finder Maker for EMC
80 //
81 //
82 // Author: Subhasis Chattopadhyay , Jan 2000
83 //
85 //
86 // StPointCollection
87 //
88 
89 //
92 #include "StPointCollection.h"
93 #include "StPi0Candidate.h"
94 #include "emc_def.h"
95 #include <Stiostream.h>
96 #include <TBrowser.h>
97 #include "StEmcUtil/geometry/StEmcGeom.h"
98 #include "StEmcUtil/projection/StEmcPosition.h"
99 #include "StEvent.h"
100 #include "StEventTypes.h"
101 #include "StarCallf77.h"
102 #include "TMath.h"
103 
104 // declaring cernlib routine (mathlib, H301) assndx to be used for matching.
105 #define assndx F77_NAME(assndx,ASSNDX)
106 extern "C"
107 {
108  void type_of_call assndx ( Int_t &, Float_t *, Int_t &, Int_t &,
109  Int_t &, Int_t *, Float_t &,Int_t*,Int_t &);
110 }
111 
112 ClassImp(StPointCollection)
113 
114 const TString detname[] =
115  {"Bemc", "Bsmde", "Bsmdp"
116  };
117 // Extern for sorted emc-smd
118 StMatchVecClus matchlist_bemc_clus[Epc::nModule][Epc::nPhiBin];
119 StMatchVecClus matchlist_bprs_clus[Epc::nModule][Epc::nPhiBin];
120 StMatchVecClus matchlist_bsmde_clus[Epc::nModule][Epc::nPhiBin];
121 StMatchVecClus matchlist_bsmdp_clus[Epc::nModule][Epc::nPhiBin];
122 
123 FloatVector HitTrackEta;
124 FloatVector HitTrackPhi;
125 FloatVector HitTrackMom;
126 
127 //_____________________________________________________________________________
128 StPointCollection::StPointCollection():TDataSet("Default")
129 {
130  SetTitle("EmcPoints");
131  mPrint = kTRUE;
132  mBField = 0.5;
133  mNPoints =0;
134  mNPointsReal=0;
135  mPosition = new StEmcPosition();
136 }
137 //_____________________________________________________________________________
138 StPointCollection::StPointCollection(const Char_t *Name):TDataSet(Name)
139 {
140  SetTitle("EmcPoints");
141  mPrint = kTRUE;
142  mBField = 0.5;
143  mNPoints =0;
144  mNPointsReal=0;
145  mPosition = new StEmcPosition();
146 }
147 //_____________________________________________________________________________
148 StPointCollection::~StPointCollection()
149 {
150  mPoints.Delete();
151  //mPointsReal.Delete(); // The objects saved here are owned by StEvent
152  mNPoints =0;
153  mNPointsReal=0;
154  delete mPosition;
155 }
156 
157 void
158 StPointCollection::Browse(TBrowser* b)
159 {
160  TDataSet::Browse(b);
161 }
162 //*************** FIND EMC POINTS **********************************
163 Int_t StPointCollection::makeEmcPoints(StEvent* event)
164 {
165  LOG_DEBUG <<"Finding EMC Points ..." << endm;
166  if(!event)
167  return 0;
168 
169  StEmcCollection* emc = event->emcCollection();
170  if(!emc)
171  return 0;
172 
173  StSPtrVecEmcPoint& points = emc->barrelPoints();
174  if(points.size())
175  points.clear();
176 
177  StEmcClusterCollection* cluster[4] = {0,0,0,0};
178 
179  for(Int_t i = 0;i<4;i++)
180  {
181  StDetectorId EmcId = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
182  ;
183  StEmcDetector* EmcDet = emc->detector(EmcId);
184  if(EmcDet)
185  cluster[i] = EmcDet->cluster();
186  }
187 
188  // first find points that were already matched in the cluster finder
189 
190  findMatchedClusters(cluster[0],cluster[1],cluster[2],cluster[3]);
191 
192  //Sort BEMC, SMDe, SMDp, PRS clusters according to location
193  LOG_DEBUG <<"Making points for non-matched clusters ..." <<endm;
194  ClusterSort(cluster[0],cluster[1],cluster[2],cluster[3]);
195 
196  // MATCHING BEMC detector clusters
197  for(Int_t im=0;im<Epc::nModule;im++)
198  {
199  for(Int_t is=0;is<Epc::nPhiBin;is++)
200  {
201  if(matchlist_bemc_clus[im][is].size()>0)
202  {
203  matchClusters(matchlist_bemc_clus[im][is],
204  matchlist_bprs_clus[im][is],
205  matchlist_bsmde_clus[im][is],
206  matchlist_bsmdp_clus[im][is]);
207  }
208  }
209  }
210 
211  // MATCHING points to tracks
212  matchToTracks(event);
213  return 1;
214 }
215 //*************** Make points for the clusters already matched ********
216 Int_t StPointCollection::findMatchedClusters(StEmcClusterCollection* Bemccluster,
217  StEmcClusterCollection *Bprscluster,
218  StEmcClusterCollection *Bsmdecluster,
219  StEmcClusterCollection *Bsmdpcluster)
220 {
221  LOG_DEBUG <<"Making points for already matched clusters ..." <<endm;
222  if(Bemccluster)
223  {
224  Int_t Ncluster0=Bemccluster->numberOfClusters();
225  if(Ncluster0>0)
226  {
227  const StSPtrVecEmcCluster& emcclusters= Bemccluster->clusters();
228  for(UInt_t i=0;i<emcclusters.size();i++)
229  {
230  StEmcCluster *btow=(StEmcCluster*)emcclusters[i];
231  UInt_t matchId = btow->GetUniqueID();
232  if(matchId!=0)
233  {
234  StEmcCluster *bprs = NULL;
235  StEmcCluster *bsmde = NULL;
236  StEmcCluster *bsmdp = NULL;
237  if(Bprscluster)
238  {
239  const StSPtrVecEmcCluster& clusters= Bprscluster->clusters();
240  for(UInt_t i=0;i<clusters.size();i++)
241  if( ((StEmcCluster*)clusters[i])->GetUniqueID() == matchId)
242  {
243  bprs = (StEmcCluster*)clusters[i];
244  break;
245  }
246  }
247  if(Bsmdecluster)
248  {
249  const StSPtrVecEmcCluster& clusters= Bsmdecluster->clusters();
250  for(UInt_t i=0;i<clusters.size();i++)
251  if( ((StEmcCluster*)clusters[i])->GetUniqueID() == matchId)
252  {
253  bsmde = (StEmcCluster*)clusters[i];
254  break;
255  }
256  }
257  if(Bsmdpcluster)
258  {
259  const StSPtrVecEmcCluster& clusters= Bsmdpcluster->clusters();
260  for(UInt_t i=0;i<clusters.size();i++)
261  if( ((StEmcCluster*)clusters[i])->GetUniqueID() == matchId)
262  {
263  bsmdp = (StEmcCluster*)clusters[i];
264  break;
265  }
266  }
267 
268  makePoint(btow,bprs,bsmde,bsmdp);
269  }
270  }
271  }
272  }
273  return 0;
274 }
275 StEmcPoint* StPointCollection::makePoint(StEmcCluster* btow,StEmcCluster* bprs,StEmcCluster* bsmde,StEmcCluster* bsmdp,Float_t fraction)
276 {
277  if(!btow)
278  return NULL;
279 
280  LOG_DEBUG <<"Making point"<< endm;
281 
282  Int_t Category = 0;
283  if(btow)
284  Category = Category | 1;
285  if(bprs)
286  Category = Category | 2;
287  if(bsmde)
288  Category = Category | 4;
289  if(bsmdp)
290  Category = Category | 8;
291 
292  Float_t en[4] = {0,0,0,0};
293  Float_t si[4] = {0,0,0,0};
294 
295  Float_t eta = btow->eta();
296  Float_t phi = btow->phi();
297 
298  // if fraction >0 it calculates the point energy as fraction*Energy from BTOW cluster
299  // if fraction <0 the point energy is taken as -fraction
300 
301  Float_t energy = btow->energy()*fraction;
302  if(fraction<0)
303  energy = fabs(fraction);
304  Float_t sigEta = btow->sigmaEta();
305  Float_t sigPhi = btow->sigmaPhi();
306  en[0] = btow->energy();
307  si[0] = sqrt(eta*eta+phi*phi);
308 
309  if(bprs)
310  {
311  en[1] = bprs->energy();
312  si[1] = sqrt(bprs->sigmaEta()*bprs->sigmaEta()+bprs->sigmaPhi()*bprs->sigmaPhi());
313  }
314 
315  if(bsmde)
316  {
317  eta = bsmde->eta();
318  sigEta = bsmde->sigmaEta();
319  en[2] = bsmde->energy();
320  si[2] = sqrt(bsmde->sigmaEta()*bsmde->sigmaEta()+bsmde->sigmaPhi()*bsmde->sigmaPhi());
321  }
322 
323  if(bsmdp)
324  {
325  phi = bsmdp->phi();
326  sigPhi = bsmdp->sigmaPhi();
327  en[3] = bsmdp->energy();
328  si[3] = sqrt(bsmdp->sigmaEta()*bsmdp->sigmaEta()+bsmdp->sigmaPhi()*bsmdp->sigmaPhi());
329  }
330 
331  Float_t xp,yp,zp;
332 
333  // Point position
334  xp=(StEpcCut::RAD_SMD_E())*cos(phi);
335  yp=(StEpcCut::RAD_SMD_E())*sin(phi);
336  zp=(StEpcCut::RAD_SMD_E())*sinh(eta);
337  StThreeVectorF PointPosition(xp*centimeter, yp*centimeter, zp*centimeter);
338 
339  //Error in location of Point
340  xp=(StEpcCut::RAD_SMD_E())*(cos(phi+sigPhi)-cos(phi-sigPhi))/2;
341  yp=(StEpcCut::RAD_SMD_E())*(sin(phi+sigPhi)-sin(phi-sigPhi))/2;
342  zp=(StEpcCut::RAD_SMD_E())*(sinh(eta+sigEta)-sinh(eta-sigEta))/2;
343  StThreeVectorF ErrorPosition(xp*centimeter, yp*centimeter, zp*centimeter);
344 
345  StThreeVectorF size(sigEta,sigPhi,0.0);
346 
347  // Chisquare
348  Float_t ChiSquare = 0.0;
349 
350 
351  StEmcPoint *point = new StEmcPoint();
352  point->setQuality(Category);
353  point->setPosition(PointPosition);
354  point->setPositionError(ErrorPosition);
355  point->setSize(size);
356  point->setChiSquare(ChiSquare);
357  point->setEnergy(energy);
358  point->setDeltaEta(9999);
359  point->setDeltaPhi(9999);
360  for(Int_t i=0;i<4;i++)
361  {
362  StDetectorId id=static_cast<StDetectorId>(i+kBarrelEmcTowerId);
363  point->setEnergyInDetector(id,en[i]);
364  point->setSizeAtDetector(id,si[i]);
365  if(i==0 && btow)
366  point->addCluster(id,btow);
367  else
368  point->addCluster(id,NULL);
369  if(i==1 && bprs)
370  point->addCluster(id,bprs);
371  else
372  point->addCluster(id,NULL);
373  if(i==2 && bsmde)
374  point->addCluster(id,bsmde);
375  else
376  point->addCluster(id,NULL);
377  if(i==3 && bsmdp)
378  point->addCluster(id,bsmdp);
379  else
380  point->addCluster(id,NULL);
381  }
382 
383  mPointsReal.Add(point);
384  mNPointsReal++;
385 
386  return point;
387 
388 }
389 //--------------------------------------------------------------------------
390 Int_t StPointCollection::matchClusters(const StMatchVecClus mvec,
391  const StMatchVecClus prsvec,
392  const StMatchVecClus evec,
393  const StMatchVecClus pvec)
394 
395 {
396  Int_t na=evec.size();
397  Int_t ma=pvec.size();
398  Float_t smin;
399  Int_t mode=1;
400  Int_t Category=-1; //this one *always* changes to 0,1,2,3 in the if-block below
401  Float_t EmcTot = 0;
402  Float_t totAvg = 0;
403  Int_t idw =Epc::nMaxNoOfClinBin;
404  Int_t ida =Epc::nMaxNoOfClinBin;
405  Float_t ep[Epc::nMaxNoOfClinBin][Epc::nMaxNoOfClinBin];
406  Int_t iw[Epc::nMaxNoOfClinBin][Epc::nMaxNoOfClinBin];
407  Int_t k[Epc::nMaxNoOfClinBin];
408 
409 
410  for (Int_t iF=0;iF<Epc::nMaxNoOfClinBin;iF++)
411  {
412  k[iF]=0;
413  for (Int_t iL=0;iL<Epc::nMaxNoOfClinBin;iL++)
414  {
415  ep[iF][iL]=0.0;
416  iw[iF][iL]=0;
417  }
418  }
419  if(evec.size()==0 && pvec.size()==0)
420  {
421  Category=0;
422  na=mvec.size();
423  ma=mvec.size();
424  }
425  if(evec.size()>0 && pvec.size()==0)
426  {
427  Category=1;
428  na=mvec.size();
429  ma=evec.size();
430  }
431  if(evec.size()==0 && pvec.size()>0)
432  {
433  Category=2;
434  na=mvec.size();
435  ma=pvec.size();
436  }
437  if(evec.size()>0 && pvec.size()>0)
438  {
439  Category=3;
440  na=evec.size();
441  ma=pvec.size();
442  }
443 
444  // getting total BTOW energy in the patch
445  if(mvec.size()>0)
446  {
447  for (UInt_t ims=0;ims<mvec.size();ims++)
448  {
449  StEmcCluster *cl0;
450  cl0=(StEmcCluster*)mvec[ims];
451  Float_t emen=cl0->energy();
452  EmcTot+=emen;
453  }
454  }
455 
456  //
457  for(Int_t ie=0;ie<na;ie++)
458  {
459  StEmcCluster *cl1 = NULL;
460  for(Int_t ip=0;ip<ma;ip++)
461  {
462  StEmcCluster *cl2 = NULL;
463  switch (Category)
464  {
465  case 0:
466  cl1 = (StEmcCluster*)mvec[ie];
467  cl2 = (StEmcCluster*)mvec[ip];
468  break;
469  case 1:
470  cl1 = (StEmcCluster*)mvec[ie];
471  cl2 = (StEmcCluster*)evec[ip];
472  break;
473  case 2:
474  cl1 = (StEmcCluster*)mvec[ie];
475  cl2 = (StEmcCluster*)pvec[ip];
476  break;
477  case 3:
478  cl1 = (StEmcCluster*)evec[ie];
479  cl2 = (StEmcCluster*)pvec[ip];
480  break;
481  }
482 
483  if(cl1 && cl2){
484  Float_t diff=TMath::Abs((cl1->energy())-(cl2->energy()));
485  Float_t summ= (cl1->energy())+(cl2->energy());
486  ep[ip][ie]=diff/summ;
487  }
488  }
489  }
490 
491  assndx(mode,ep[0],na,ma,ida,k,smin,iw[0],idw);
492 
493  int i1;
494  switch (Category)
495  {
496  case 0:
497  for(i1=0;i1<na;i1++)
498  {
499  if((k[i1]-1)>=0)
500  {
501  StEmcCluster *cl1;
502  cl1 = (StEmcCluster*)mvec[i1];
503  Float_t avg_en = cl1->energy();
504  totAvg += avg_en;
505  }
506  }
507  break;
508  case 1:
509  for(i1=0;i1<na;i1++)
510  {
511  if((k[i1]-1)>=0)
512  {
513  StEmcCluster *cl1;
514  cl1 = (StEmcCluster*)mvec[i1];
515  Float_t avg_en = cl1->energy();
516  totAvg += avg_en;
517  }
518  }
519  break;
520  case 2:
521  for(i1=0;i1<na;i1++)
522  {
523  if((k[i1]-1)>=0)
524  {
525  StEmcCluster *cl1;
526  cl1 = (StEmcCluster*)mvec[i1];
527  Float_t avg_en = cl1->energy();
528  totAvg += avg_en;
529  }
530  }
531  break;
532  case 3:
533  for(i1=0;i1<na;i1++)
534  {
535  if((k[i1]-1)>=0)
536  {
537  StEmcCluster *cl1;
538  cl1 = (StEmcCluster*)evec[i1];
539  StEmcCluster *cl2;
540  cl2 = (StEmcCluster*)pvec[k[i1]-1];
541  Float_t avg_en = 0.5*(cl1->energy()+cl2->energy());
542  totAvg += avg_en;
543  }
544  }
545  break;
546  }
547 
548  for(i1=0;i1<na;i1++)
549  {
550 
551  StEmcCluster *btow = NULL;
552  StEmcCluster *bprs = NULL;
553  StEmcCluster *bsmde = NULL;
554  StEmcCluster *bsmdp = NULL;
555  Float_t fraction = 1;
556  Float_t eta, phi;
557  eta = phi = -999.; //these should always be initialized, but it would be nice to prove it
558 
559  switch (Category)
560  {
561  case 0:
562  if((k[i1]-1)>=0)
563  {
564  btow = (StEmcCluster*)mvec[i1];
565  eta = btow->eta();
566  phi = btow->phi();
567  fraction = 1;
568  }
569  break;
570  case 1:
571  if((k[i1]-1)>=0)
572  {
573  btow = (StEmcCluster*)mvec[i1];
574  bsmde = (StEmcCluster*)evec[k[i1]-1];
575  eta = bsmde->eta();
576  phi = btow->phi();
577  fraction = EmcTot/totAvg;
578  }
579  break;
580  case 2:
581  if((k[i1]-1)>=0)
582  {
583  btow = (StEmcCluster*)mvec[i1];
584  bsmdp = (StEmcCluster*)pvec[k[i1]-1];
585  eta=btow->eta();
586  phi=bsmdp->phi();
587  fraction = EmcTot/totAvg;
588  }
589  break;
590  case 3:
591  if((k[i1]-1)>=0)
592  {
593  bsmde = (StEmcCluster*)evec[i1];
594  bsmdp = (StEmcCluster*)pvec[k[i1]-1];
595  fraction = -fabs(EmcTot*0.5*(bsmde->energy()+bsmdp->energy())/totAvg);
596  eta = bsmde->eta();
597  phi = bsmdp->phi();
598  Float_t delta = 999999;
599  for (UInt_t ims=0;ims<mvec.size();ims++)
600  {
601  StEmcCluster *cl0 = (StEmcCluster*)mvec[ims];
602  Float_t de = sqrt((eta-cl0->eta())*(eta-cl0->eta()) +
603  (phi-cl0->phi())*(phi-cl0->phi()));
604  if(de<delta)
605  {
606  btow = cl0;
607  delta = de;
608  }
609  }
610  }
611  break;
612  }
613 
614  if(prsvec.size()>0)
615  {
616  Float_t delta = 999999;
617  for (UInt_t ims=0;ims<prsvec.size();ims++)
618  {
619  StEmcCluster *cl0 = (StEmcCluster*)prsvec[ims];
620  Float_t de = sqrt((eta-cl0->eta())*(eta-cl0->eta()) +
621  (phi-cl0->phi())*(phi-cl0->phi()));
622  if(de<delta)
623  {
624  bprs = cl0;
625  delta = de;
626  }
627  }
628  }
629  makePoint(btow,bprs,bsmde,bsmdp,fraction);
630  }
631  return 1;
632 }
633 //-----------------------------------------------------------------------
634 void StPointCollection::ClusterSort(StEmcClusterCollection* Bemccluster,
635  StEmcClusterCollection* Bprscluster,
636  StEmcClusterCollection* Bsmdecluster,
637  StEmcClusterCollection* Bsmdpcluster)
638 {
639  // The TObject::GetUniqueID() method is used in the cluster
640  // finder to allow matching at clustering level. StEpcMaker
641  // will try to match only clusters which have UniqueID = 0
642  // clusters with UniqueID!=0 were matched at cluster finder
643  // StEpcMaker will, in these cases, use the UniqueID as plain
644  // matching information to create the corresponding points.
645  const Int_t eta_shift_fix=1;
646  LOG_DEBUG <<" I am inside PointCalc***"<<endm;
647  for(Int_t i1=0;i1<Epc::nModule;i1++)
648  {
649  for(Int_t i2=0;i2<Epc::nPhiBin;i2++)
650  {
651  matchlist_bemc_clus[i1][i2].clear();
652  matchlist_bprs_clus[i1][i2].clear();
653  matchlist_bsmde_clus[i1][i2].clear();
654  matchlist_bsmdp_clus[i1][i2].clear();
655  }
656  }
657 
658  StEmcGeom* GeomIn = StEmcGeom::getEmcGeom("bemc");
659 
660  //BEMC
661  if(Bemccluster)
662  {
663  Int_t Ncluster0=Bemccluster->numberOfClusters();
664  if(Ncluster0>0)
665  {
666  const StSPtrVecEmcCluster& emcclusters= Bemccluster->clusters();
667  for(UInt_t i=0;i<emcclusters.size();i++)
668  {
669  StEmcCluster *cl1=(StEmcCluster*)emcclusters[i];
670  LOG_DEBUG <<"BEMC cluster UniqueId = "<<cl1->GetUniqueID()<<endm;
671  if(cl1->GetUniqueID()==0)
672  {
673  Float_t eta_emc=cl1->eta();
674  Float_t phi_emc=cl1->phi();
675  //Get the module number
676  Int_t ebin,pbin;
677  Int_t emc_module=0;
678  Int_t & imd =emc_module;
679  Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
680  if(testb==0)
681  emc_module=imd;
682  if(testb==0)
683  {
684  Int_t emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
685  //keeping the cluster very close to phibin boundry to the previous bin
686  //
687  if(!eta_shift_fix && emc_phi_bin>0)
688  {
689  if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<0.2)
690  {
691  emc_phi_bin--;
692  }
693  }
694  if(emc_phi_bin>9)
695  {
696  emc_phi_bin=9;
697  }
698  //copy cl1 pointer to StMatchvec
699  matchlist_bemc_clus[emc_module-1][emc_phi_bin].push_back(cl1);
700  }
701  }
702  }
703  }
704  }
705 
706  //PRS
707  if(Bprscluster)
708  {
709  Int_t Ncluster1=Bprscluster->numberOfClusters();
710  if(Ncluster1>0)
711  {
712  const StSPtrVecEmcCluster& emcclusters= Bprscluster->clusters();
713  for(UInt_t i=0;i<emcclusters.size();i++)
714  {
715  StEmcCluster *cl2=(StEmcCluster*)emcclusters[i];
716  if(cl2->GetUniqueID()==0)
717  {
718  Float_t eta_emc=cl2->eta();
719  Float_t phi_emc=cl2->phi();
720  //Get the module number
721  Int_t ebin,pbin;
722  Int_t emc_module=0;
723  Int_t emc_phi_bin=0;
724  Int_t & imd =emc_module;
725  Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
726  if(testb==0)
727  emc_module=imd;
728  if(testb==0)
729  {
730  emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
731  //keeping the cluster very close to phibin boundry to the previous bin
732  if(!eta_shift_fix && emc_phi_bin>0)
733  {
734  if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<0.2)
735  {
736  emc_phi_bin--;
737  }
738  }
739  if(emc_phi_bin>9)
740  {
741  emc_phi_bin=9;
742  }
743  //copy cl1 pointer to StMatchvec
744  matchlist_bprs_clus[emc_module-1][emc_phi_bin].push_back(cl2);
745  }
746  }
747  }
748  }
749  }
750 
751  //BSMD_ETA
752  if(Bsmdecluster)
753  {
754  Int_t Ncluster2=Bsmdecluster->numberOfClusters();
755  if(Ncluster2>0)
756  {
757  const StSPtrVecEmcCluster& emcclusters= Bsmdecluster->clusters();
758  for(UInt_t i=0;i<emcclusters.size();i++)
759  {
760  StEmcCluster *cl3=(StEmcCluster*)emcclusters[i];
761  if(cl3->GetUniqueID()==0)
762  {
763  Float_t eta_emc=cl3->eta();
764  Float_t phi_emc=cl3->phi();
765  //Get the module number
766  Int_t ebin,pbin;
767  Int_t emc_module=0;
768  Int_t emc_phi_bin=0;
769  Int_t & imd =emc_module;
770  Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
771  if(testb==0)
772  emc_module=imd;
773  if(testb==0)
774  {
775  emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
776  //keeping the cluster very close to phibin boundry to the previous bin
777  if(!eta_shift_fix && emc_phi_bin>0)
778  {
779  if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<.2)
780  {
781  emc_phi_bin--;
782  }
783  }
784  if(emc_phi_bin>9)
785  {
786  emc_phi_bin=9;
787  }
788  //copy cl1 pointer to StMatchvec
789  matchlist_bsmde_clus[emc_module-1][emc_phi_bin].push_back(cl3);
790  }
791  }
792  }
793  }
794  }
795 
796  // BSMDP
797  if(Bsmdpcluster)
798  {
799  Int_t Ncluster3=Bsmdpcluster->numberOfClusters();
800  if(Ncluster3>0)
801  {
802  const StSPtrVecEmcCluster& emcclusters= Bsmdpcluster->clusters();
803  for(UInt_t i=0;i<emcclusters.size();i++)
804  {
805  StEmcCluster *cl4=(StEmcCluster*)emcclusters[i];
806  if(cl4->GetUniqueID()==0)
807  {
808  Float_t eta_emc=cl4->eta();
809  Float_t phi_emc=cl4->phi();
810  //Get the module number
811  Int_t ebin,pbin;
812  Int_t emc_module=0;
813  Int_t & imd =emc_module;
814  Int_t testb=GeomIn->getBin(phi_emc,eta_emc,imd,ebin,pbin);
815  if(testb==0)
816  emc_module=imd;
817  if(testb==0)
818  {
819  Int_t emc_phi_bin=Int_t(TMath::Abs(eta_emc*10));
820  //keeping the cluster very close to phibin boundry to the previous bin
821  if(!eta_shift_fix && emc_phi_bin>0)
822  {
823  if((TMath::Abs(eta_emc*10)-Float_t(emc_phi_bin))<=0.01)
824  {
825  emc_phi_bin--;
826  }
827  }
828  if(emc_phi_bin>9)
829  {
830  emc_phi_bin=9;
831  }
832  //copy cl1 pointer to StMatchvec
833  matchlist_bsmdp_clus[emc_module-1][emc_phi_bin].push_back(cl4);
834  }
835  }
836  }
837  }
838  }
839 }
840 //--------------------------------------------------------------------------
841 Int_t StPointCollection::matchToTracks(StEvent* event)
842 {
843  if(!event)
844  return 0;
845 
846  Float_t field = 0.5;
847  StEventSummary *evtSummary = event->summary();
848  if (evtSummary)
849  field = evtSummary->magneticField()/10;
850 
851  Int_t nR = NPointsReal();
852  StEmcGeom* geom = StEmcGeom::instance("bemc");
853  if(nR>0)
854  {
855  LOG_DEBUG << "Matching to tracks... NP = " << nR << endm;
856  StSPtrVecTrackNode& tracks=event->trackNodes();
857  Int_t nTracks = tracks.size();
858  StThreeVectorD momentum,position;
859  for(Int_t t=0;t<nTracks;t++)
860  {
861  StTrack *track = tracks[t]->track(0);
862  if(track)
863  {
864  if(track->geometry())
865  {
866  Bool_t tok = mPosition->trackOnEmc(&position,&momentum,
867  track,(Double_t)field,
868  (Double_t)geom->Radius());
869  if(tok)
870  {
871  Float_t eta = position.pseudoRapidity();
872  Float_t phi = position.phi();
873  if(fabs(eta)<1)
874  {
875  TIter next(PointsReal());
876  StEmcPoint *cl;
877 
878  for(Int_t i=0; i<nR; i++)
879  {
880  cl = (StEmcPoint*)next();
881  if(cl)
882  {
883  StThreeVectorF pos = cl->position();
884  Float_t etaP = pos.pseudoRapidity();
885  Float_t phiP = pos.phi();
886  Float_t D = sqrt(cl->deltaEta()*cl->deltaEta()+cl->deltaPhi()*cl->deltaPhi());
887  Float_t d = sqrt((eta-etaP)*(eta-etaP)+(phi-phiP)*(phi-phiP));
888  if(d<D)
889  {
890  cl->setDeltaEta(eta-etaP);
891  cl->setDeltaPhi(phi-phiP);
892  }
893 
894  StThreeVectorF err = cl->positionError();
895  Float_t etaE = err.pseudoRapidity();
896  Float_t phiE = err.phi();
897 
898  Float_t dPhi = fabs(phi-phiP);
899  if (dPhi>TMath::Pi())
900  dPhi=2*TMath::Pi()-dPhi;
901 
902  if(fabs(eta-etaP)<fabs(etaE) && dPhi<fabs(phiE))
903  {
904  Int_t Category = cl->quality();
905  Category = Category | 16;
906  cl->setQuality(Category);
907  cl->addTrack(track);
908  }
909  }
910  }
911  }
912  }
913  }
914  }
915  }
916  }
917 
918  return 0;
919 }
StEmcPosition::trackOnEmc
Bool_t trackOnEmc(StThreeVectorD *position, StThreeVectorD *momentum, const StTrack *const track, Double_t magField, Double_t emcRadius=225.405) const
Track projection utility.
Definition: StEmcPosition.cxx:270
StPointCollection::Browse
virtual void Browse(TBrowser *b)
Browse this dataset (called by TBrowser).
Definition: StPointCollection.cxx:158
TDataSet::Browse
virtual void Browse(TBrowser *b)
Browse this dataset (called by TBrowser).
Definition: TDataSet.cxx:297
StEmcGeom::getBin
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
Definition: StEmcGeom.h:321
Generated by   doxygen 1.8.5