StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StGammaCandidateMaker.cxx
1 #include "TVector2.h"
2 
3 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
4 #include "StEEmcPool/StEEmcClusterMaker/StEEmcGenericClusterMaker.h"
5 
6 #include "SystemOfUnits.h"
7 #include "StEmcUtil/geometry/StEmcGeom.h"
8 #include "StEmcUtil/projection/StEmcPosition.h"
9 
10 #include "StBarrelEmcCluster.h"
11 #include "StBarrelEmcClusterMaker.h"
12 #include "StGammaFitter.h"
13 #include "StGammaEvent.h"
14 #include "StGammaEventMaker.h"
15 #include "StGammaRawMaker.h"
16 
17 #include "StGammaCandidateMaker.h"
18 
19 ClassImp(StGammaCandidate);
20 
22 // Constructor //
24 StGammaCandidateMaker::StGammaCandidateMaker(const char *name): StMaker(name)
25 {
26 
27  mUseBemc = false;
28  mUseEemc = false;
29 
30  mStrictBemcStatus = true;
31 
32  mCompressLevel = 2;
33 
34  mMinimumEt = 0.0;
35  mRadius = 0.7;
36  mBsmdRange = 0.05237;
37  mEsmdRange = 20.0;
38 
39 }
40 
42 // Destructor //
44 StGammaCandidateMaker::~StGammaCandidateMaker()
45 {}
46 
48 // Maker Init //
50 Int_t StGammaCandidateMaker::Init()
51 {
52  return StMaker::Init();
53 }
54 
56 // Maker Clear //
58 void StGammaCandidateMaker::Clear( Option_t *opts )
59 {
60  mId = 0;
61 }
62 
64 // Maker Make //
66 Int_t StGammaCandidateMaker::Make()
67 {
68 
69  if(mUseBemc) MakeBarrel();
70  if(mUseEemc) MakeEndcap();
71 
72  Compress();
73 
74  return kStOK;
75 
76 }
77 
78 
80 // Return the vertex corrected momentum //
81 // for the given EEMC cluster //
83 TVector3 momentum(StEEmcCluster cluster, TVector3 vertex)
84 {
85 
86  TVector3 d = cluster.momentum().Unit();
87 
88  d.SetMag( kEEmcZSMD / d.CosTheta() );
89  d = d - vertex;
90  d = d.Unit();
91  d *= cluster.energy();
92 
93  return d;
94 
95 }
96 
98 // Compile BEMC Clusters //
100 Int_t StGammaCandidateMaker::MakeBarrel()
101 {
102 
103  // Retrieve all those gamma makers from the chain
104  StGammaEventMaker *mGammaEventMaker = dynamic_cast<StGammaEventMaker*>(GetMakerInheritsFrom("StGammaEventMaker"));
105  if(!mGammaEventMaker)
106  {
107  LOG_WARN << "MakeBarrel() - No StGammaEventMaker found!" << endm;
108  return kStWarn;
109  }
110 
111  StGammaEvent *mGammaEvent = mGammaEventMaker->event();
112  if(!mGammaEvent)
113  {
114  LOG_WARN << "MakeBarrel() - No StGammaEvent found!" << endm;
115  return kStWarn;
116  }
117 
118  StGammaRawMaker *mGammaRawMaker = dynamic_cast<StGammaRawMaker*>(GetMakerInheritsFrom("StGammaRawMaker"));
119  if(!mGammaRawMaker)
120  {
121  LOG_WARN << "MakeBarrel() - No StGammaRawMaker found!" << endm;
122  return kStWarn;
123  }
124 
125  StBarrelEmcClusterMaker *mBemcClusterMaker = dynamic_cast<StBarrelEmcClusterMaker*>(GetMakerInheritsFrom("StBarrelEmcClusterMaker"));
126  if(!mBemcClusterMaker)
127  {
128  LOG_WARN << "MakeBarrel() - No StBarrelEmcClusterMaker found!" << endm;
129  return kStWarn;
130  }
131 
132  // Create a StEmcGeom instance for track projections
133  StEmcGeom* geom = StEmcGeom::instance("bemc");
134 
135  StEmcGeom *smdEtaGeom = new StEmcGeom("bsmde");
136  StEmcGeom *smdPhiGeom = new StEmcGeom("bsmdp");
137 
138  StEmcPosition emcPosition;
139 
140  // Loop over BEMC clusters
141  for(unsigned int i = 0; i < mBemcClusterMaker->clusters().size(); ++i)
142  {
143 
144  StBarrelEmcCluster* cluster = mBemcClusterMaker->clusters()[i];
145 
146  // Require minimum Et
147  if(cluster->momentum().Pt() < mMinimumEt) continue;
148 
149  // Require that seed towers have non fail status
150  if(cluster->seed()->fail) continue;
151 
152  // If strict status is observed, require that all
153  // neighboring towers also have a non fail status
154  if(mStrictBemcStatus)
155  {
156 
157  bool fail = false;
158 
159  for(int deta = -1; deta <= 1; ++deta)
160  {
161  for(int dphi = -1; dphi <= 1; ++dphi)
162  {
163  if(StGammaTower *tower = cluster->tower(deta, dphi))
164  {
165  if(tower->fail) fail = true;
166  }
167  }
168  }
169 
170  if(fail) continue;
171 
172  }
173 
174  // Create gamma candidate
175  StGammaCandidate *candidate = mGammaEvent->newCandidate();
176  candidate->SetTowerClusterId(i);
177  candidate->SetId(nextId());
178  candidate->SetDetectorId(StGammaCandidate::kBEmc);
179  candidate->SetEnergy(cluster->energy());
180  candidate->SetPosition(cluster->position());
181  candidate->SetMomentum(cluster->momentum());
182 
183  // Start with the seed tower
184  StGammaTower *seed = cluster->seed();
185  candidate->SetSeedEnergy(seed->energy);
186  candidate->SetTowerId(seed->id);
187 
188  // Loop over neighboring towers in the cluster
189  float energy = 0;
190 
191  for(int deta = -1; deta <= 1; ++deta)
192  {
193 
194  for(int dphi = -1; dphi <= 1; ++dphi)
195  {
196 
197  if(StGammaTower *tower = cluster->tower(deta, dphi))
198  {
199 
200  // Skip any failing towers
201  if(tower->fail) continue;
202 
203  // Associate tower with the candidate
204  candidate->addMyTower(tower);
205  tower->candidates.Add(candidate);
206 
207  // Associate preshower with the candidate
208  if(StGammaTower *preshower = mGammaRawMaker->tower(tower->id, kBEmcPres))
209  {
210  candidate->addMyPreshower1(preshower);
211  preshower->candidates.Add(candidate);
212  energy += preshower->energy;
213  }
214 
215  // Associate tracks with the candidate
216  for(int k = 0; k < mGammaEvent->numberOfTracks(); ++k)
217  {
218 
219  StGammaTrack* track = mGammaEvent->track(k);
220  if(!track) continue;
221 
222  TVector3 position = track->positionAtRadius(geom->Radius());
223  if(position != TVector3(0, 0, 0))
224  {
225 
226  int id;
227  if(geom->getId(position.Phi(), position.Eta(), id) == 0 && id == static_cast<int>(tower->id))
228  {
229  candidate->addMyTrack(track);
230  track->candidates.Add(candidate);
231  }
232 
233  }
234 
235  }
236 
237  } // if(tower)
238 
239  } // dphi
240 
241  } // deta
242 
243 
244  // Set candidate preshower sum
245  candidate->SetBPrsEnergy(energy);
246 
247  // Associate tracks falling within a cone of mRadius around the candidate
248  for (int k = 0; k < mGammaEvent->numberOfTracks(); ++k)
249  {
250 
251  if (StGammaTrack* track = mGammaEvent->track(k))
252  {
253 
254  float deta = cluster->momentum().Eta() - track->eta();
255  float dphi = cluster->momentum().Phi() - track->phi();
256  dphi = TVector2::Phi_mpi_pi(dphi);
257  float r = hypot(deta, dphi);
258 
259  if(r <= mRadius)
260  {
261  candidate->addTrack(track);
262  if(!track->candidates.FindObject(candidate)) track->candidates.Add(candidate);
263  }
264 
265  }
266 
267  }
268 
269  // Associate towers falling within a cone of radius mRadius around the candidate
270  for(int k = 0; k < mGammaEvent->numberOfTowers(); ++k)
271  {
272 
273  StGammaTower* tower = mGammaEvent->tower(k);
274  if(tower && tower->layer == kBEmcTower)
275  {
276 
277  float deta = cluster->momentum().Eta() - tower->eta;
278  float dphi = cluster->momentum().Phi() - tower->phi;
279  dphi = TVector2::Phi_mpi_pi(dphi);
280  float r = hypot(deta, dphi);
281 
282  if(r <= mRadius)
283  {
284  candidate->addTower(tower);
285  if(!tower->candidates.FindObject(candidate)) tower->candidates.Add(candidate);
286 
287  // Associate tracks projecting to the tower with the tower
288  for(int t = 0; t < mGammaEvent->numberOfTracks(); ++t)
289  {
290 
291  StGammaTrack* track = mGammaEvent->track(t);
292  if(!track) continue;
293 
294  TVector3 position = track->positionAtRadius(geom->Radius());
295  if(position != TVector3(0, 0, 0))
296  {
297 
298  int id;
299  if(geom->getId(position.Phi(), position.Eta(), id) == 0 && id == static_cast<int>(tower->id))
300  {
301  if(!tower->tracks.FindObject(track)) tower->tracks.Add(track);
302  }
303 
304  }
305 
306  }
307 
308  }
309 
310  }
311 
312  }
313 
314  // Associate preshower hits falling within a cone of radius mRadius around the candidate
315  for(int k = 0; k < mGammaEvent->numberOfPreshower1(); ++k)
316  {
317 
318  StGammaTower* preshower = mGammaEvent->preshower1(k);
319  if(preshower && preshower->layer == kBEmcPres)
320  {
321 
322  float deta = cluster->momentum().Eta() - preshower->eta;
323  float dphi = cluster->momentum().Phi() - preshower->phi;
324  dphi = TVector2::Phi_mpi_pi(dphi);
325  float r = hypot(deta, dphi);
326 
327  if (r <= mRadius)
328  {
329  candidate->addPreshower1(preshower);
330  if(!preshower->candidates.FindObject(candidate)) preshower->candidates.Add(candidate);
331  }
332 
333  }
334 
335  }
336 
337  // Associate BSMDE and BSMDP strips falling within mSmdRange of the candidate
338  float smdEtaEnergy = 0;
339  float smdPhiEnergy = 0;
340 
341  int centralId = 0;
342  int module = 0;
343  int eta = 0;
344  int sub = 0;
345 
346  smdEtaGeom->getId(cluster->position().Phi(), cluster->position().Eta(), centralId);
347 
348  int sector = 0; // Dummy variable in the BEMC
349 
350  for(int dPhiStrip = -1; dPhiStrip <= 1; ++dPhiStrip)
351  {
352 
353  for(int dEtaStrip = -20; dEtaStrip <= 20; ++dEtaStrip)
354  {
355 
356  int id = emcPosition.getNextId(3, centralId, dEtaStrip, dPhiStrip);
357  if(id == 0) continue;
358 
359  smdEtaGeom->getBin(id, module, eta, sub);
360 
361  // Find vertex corrected position of the strip
362  float x, y, z;
363  smdEtaGeom->getXYZ(id, x, y, z);
364  TVector3 vEta(x, y, z);
365 
366  float deta = cluster->position().Eta() - vEta.Eta();
367  float dphi = cluster->position().Phi() - vEta.Phi();
368  dphi = TVector2::Phi_mpi_pi(dphi);
369  float r = hypot(deta, dphi);
370 
371  // If the strip is within range then associate the strip
372  if(r <= mBsmdRange)
373  {
374 
375  // Use strip is already created in StGammaRawMaker
376  if(StGammaStrip* strip = mGammaRawMaker->strip(sector, kBEmcSmdEta, id))
377  {
378 
379  candidate->addSmdEta(strip);
380  strip->candidates.Add(candidate);
381  smdEtaEnergy += strip->energy;
382 
383  }
384  // Otherwise create a new, empty strip
385  else
386  {
387 
388  float eta;
389 
390  smdEtaGeom->getEta(id, eta);
391 
392  strip = mGammaEvent->newStrip();
393 
394  strip->index = id;
395  strip->sector = module;
396  strip->plane = kBEmcSmdEta;
397  // strip->stat Filled in StGammaRawMaker::AddEtaStrip()
398  // strip->fail Filled in StGammaRawMaker::AddEtaStrip()
399  strip->energy = 0;
400  strip->position = 230.705 * sinh(eta);
401 
402  double offset = fabs(eta) < 0.5 ? 0.73 : 0.94;
403 
404  strip->left = strip->position - offset;
405  strip->right = strip->position + offset;
406 
407  mGammaRawMaker->AddEtaStrip(strip);
408 
409  candidate->addSmdEta(strip);
410  strip->candidates.Add(candidate);
411 
412  } // if exist
413 
414  } // if range
415 
416  } // dEtaStrips
417 
418  } // dPhiStrips
419 
420  centralId = 0;
421  smdPhiGeom->getId(cluster->position().Phi(), cluster->position().Eta(), centralId);
422 
423  for(int dEtaStrip = -1; dEtaStrip <= 1; ++dEtaStrip)
424  {
425 
426  for(int dPhiStrip = -20; dPhiStrip <= 20; ++dPhiStrip)
427  {
428 
429  int id = emcPosition.getNextId(4, centralId, dEtaStrip, dPhiStrip);
430  if(id == 0) continue;
431 
432  smdEtaGeom->getBin(id, module, eta, sub);
433 
434  // Find vertex corrected position of the strip
435  float x, y, z;
436  smdPhiGeom->getXYZ(id, x, y, z);
437  TVector3 vPhi(x, y, z);
438 
439  float deta = cluster->position().Eta() - vPhi.Eta();
440  float dphi = cluster->position().Phi() - vPhi.Phi();
441  dphi = TVector2::Phi_mpi_pi(dphi);
442  float r = hypot(deta, dphi);
443 
444  // If the strip is within range then associate the strip
445  if(r < mBsmdRange)
446  {
447 
448  // Use strip is already created in StGammaRawMaker
449  if(StGammaStrip* strip = mGammaRawMaker->strip(sector, kBEmcSmdPhi, id))
450  {
451 
452  candidate->addSmdPhi(strip);
453  strip->candidates.Add(candidate);
454  smdPhiEnergy += strip->energy;
455 
456  }
457  // Otherwise create a new, empty strip
458  else
459  {
460 
461  float phi;
462 
463  smdPhiGeom->getPhi(id, phi);
464 
465  strip = mGammaEvent->newStrip();
466 
467  strip->index = id;
468  strip->sector = module;
469  strip->plane = kBEmcSmdPhi;
470  //strip->stat Filled in StGammaRawMaker::AddPhiStrip()
471  //strip->fail Filled in StGammaRawMaker::AddPhiStrip()
472  strip->energy = 0;
473  strip->position = phi;
474 
475  double offset = 0.00293;
476 
477  strip->left = phi - offset;
478  strip->right = phi + offset;
479 
480  mGammaRawMaker->AddPhiStrip(strip);
481 
482  candidate->addSmdPhi(strip);
483  strip->candidates.Add(candidate);
484 
485  } // if exist
486 
487  } // if range
488 
489  } // dEtaStrips
490 
491  } // dPhiStrips
492 
493  candidate->SetSmdEtaEnergy(smdEtaEnergy);
494  candidate->SetSmdPhiEnergy(smdPhiEnergy);
495 
496  } // Clusters
497 
498  return kStOK;
499 
500 }
501 
503 // Compile EEMC Clusters //
505 Int_t StGammaCandidateMaker::MakeEndcap()
506 {
507 
508  // Retrieve all those gamma makers from the chain
509  StGammaEventMaker *mGammaEventMaker = dynamic_cast<StGammaEventMaker*>(GetMakerInheritsFrom("StGammaEventMaker"));
510  if(!mGammaEventMaker)
511  {
512  LOG_WARN << "MakeEndcap() - No StGammaEventMaker found!" << endm;
513  return kStWarn;
514  }
515 
516  StGammaEvent *mGammaEvent = mGammaEventMaker->event();
517  if(!mGammaEvent)
518  {
519  LOG_WARN << "MakeEndcap() - No StGammaEvent found!" << endm;
520  return kStWarn;
521  }
522 
523  StGammaRawMaker *mGammaRawMaker = dynamic_cast<StGammaRawMaker*>(GetMakerInheritsFrom("StGammaRawMaker"));
524  if(!mGammaRawMaker)
525  {
526  LOG_WARN << "MakeEndcap() - No StGammaRawMaker found!" << endm;
527  return kStWarn;
528  }
529 
530  StEEmcGenericClusterMaker *mEEclusters = dynamic_cast<StEEmcGenericClusterMaker*>(GetMakerInheritsFrom("StEEmcGenericClusterMaker"));
531  if(!mEEclusters)
532  {
533  LOG_DEBUG << "MakeEndcap() - No StEEmcGenericClusterMaker found!" << endm;
534  return kStWarn;
535  }
536 
537  StEEmcA2EMaker *mEEanalysis = dynamic_cast<StEEmcA2EMaker*>(GetMakerInheritsFrom("StEEmcA2EMaker"));
538  if(!mEEanalysis)
539  {
540  LOG_DEBUG << "MakeEndcap() - No StEEmcA2EMaker found!" << endm;
541  return kStWarn;
542  }
543 
544 
545  // Loop over each sector,
546  for(Int_t sector = 0; sector < kEEmcNumSectors; sector++)
547  {
548 
549  // then the clusters in each sector
550  StEEmcClusterVec_t clusters = mEEclusters->clusters(sector, 0);
551  for(UInt_t i = 0; i < clusters.size(); i++)
552  {
553 
554  // Correct the cluster momentum for vertex
555  StEEmcCluster cluster = clusters[i];
556  TVector3 position = getEEmcClusterPosition(cluster);
557  if(position == TVector3(0,0,0)) position = cluster.position();
558  TVector3 pcluster = position - mGammaEvent->vertex();
559  pcluster.SetMag(cluster.energy());
560 
561  // Require minimum transverse energy
562  Float_t ET = pcluster.Perp();
563  if(ET < mMinimumEt) continue;
564 
565  // Creat a new candidate
566  StGammaCandidate *mCandidate = mGammaEvent->newCandidate();
567 
568  mCandidate->SetId( nextId() );
569  mCandidate->SetTowerClusterId( cluster.key() );
570  mCandidate->SetDetectorId( StGammaCandidate::kEEmc );
571 
572  // Set candidate kinematics
573  mCandidate->SetMomentum( pcluster );
574  mCandidate->SetPosition( cluster.position() );
575  mCandidate->SetEnergy( cluster.energy() );
576 
577  // Set seed information
578  StEEmcTower seed = cluster.tower(0);
579  mCandidate->SetSeedEnergy( seed.energy() );
580  mCandidate->SetTowerId( seed.index() );
581 
582 
583  // Associate towers, preshowers, and postshowers with the candidate
584  // while summing up the total energy response in each
585  Float_t epre1 = 0.;
586  Float_t epre2 = 0.;
587  Float_t epost = 0.;
588 
589  // Loop over neighboring towers
590  for(Int_t j = 0; j < cluster.numberOfTowers(); j++)
591  {
592 
593  // Grab the jth neighboring tower
594  StEEmcTower tower = cluster.tower(j);
595  Int_t index = tower.index();
596 
597  // and its pre/postshowers
598  StEEmcTower pre1 = mEEanalysis->tower(index, 1);
599  StEEmcTower pre2 = mEEanalysis->tower(index, 2);
600  StEEmcTower post = mEEanalysis->tower(index, 3);
601 
602  // Increment the energy sums
603  if( !pre1.fail() ) epre1 += pre1.energy();
604  if( !pre2.fail() ) epre2 += pre2.energy();
605  if( !post.fail() ) epost += post.energy();
606 
607  // Remember, if the respective detector element did not have
608  // sufficient energy deposited in it then these pointers
609  // will be null
610  StGammaTower *gtower = mGammaRawMaker->tower(index, kEEmcTower);
611  StGammaTower *gpre1 = mGammaRawMaker->tower(index, kEEmcPre1);
612  StGammaTower *gpre2 = mGammaRawMaker->tower(index, kEEmcPre2);
613  StGammaTower *gpost = mGammaRawMaker->tower(index, kEEmcPost);
614 
615 
616  // Associate clustered towers with the candidate
617  if(gtower)
618  {
619  mCandidate->addMyTower(gtower);
620  gtower->candidates.Add(mCandidate);
621  }
622 
623  // Associate clustered preshower1 tiles with the candidate
624  if(gpre1 && !pre1.fail() && pre1.energy() > 0.0)
625  {
626  mCandidate->addMyPreshower1(gpre1);
627  gpre1->candidates.Add(mCandidate);
628  }
629 
630  // Associate clustered preshower2 tiles with the candidate
631  if(gpre2 && !pre2.fail() && pre2.energy() > 0.0)
632  {
633  mCandidate->addMyPreshower2(gpre2);
634  gpre2->candidates.Add(mCandidate);
635  }
636 
637  // Associate clustered postshower tiles with the candidate
638  if(gpost && !post.fail() && post.energy() > 0.0)
639  {
640  mCandidate->addMyPostshower(gpost);
641  gpost->candidates.Add(mCandidate);
642  }
643 
644  // Associate tracks extrapolating to the tower to the cluster
645  if(gtower)
646  {
647 
648  // Loop over all tracks in the event
649  for(int k = 0; k < mGammaEvent->numberOfTracks(); ++k)
650  {
651 
652  StGammaTrack* track = mGammaEvent->track(k);
653 
654  // Skip misbehaving or null tracks
655  if (!track || track->pz() < 0) continue;
656 
657  // Extrapolate track to the SMD distance
658  EEmcGeomSimple& geom = EEmcGeomSimple::Instance();
659  TVector3 position = track->positionAtZ(geom.getZ1());
660 
661  if(position != TVector3(0, 0, 0))
662  {
663 
664  // If the tower to which the track extrapolates is equal
665  // to the current tower then associate track with cluster
666  int sector, subsector, etabin;
667 
668  if (geom.getTower(position, sector, subsector, etabin))
669  {
670 
671  bool match = gtower->sector() == sector;
672  match &= gtower->subsector() == subsector;
673  match &= gtower->etabin() == etabin;
674 
675  if(match)
676  {
677  mCandidate->addMyTrack(track);
678  track->candidates.Add(mCandidate);
679  }
680 
681  }
682 
683  }
684 
685 
686  } // Tracks
687 
688  } // if(gtower)
689 
690  } // Cluster towers
691 
692 
693  // Set preshower and postshower energy sums
694  mCandidate->SetPre1Energy(epre1);
695  mCandidate->SetPre2Energy(epre2);
696  mCandidate->SetPostEnergy(epost);
697 
698  // Now associate tracks, towers, etc, lying within a cone
699  // of mRadium around the candidate to the candidate
700  // Only these objects will be stored in the StGammaEvent
701 
702  Float_t candidateEta = pcluster.Eta();
703  Float_t candidatePhi = pcluster.Phi();
704 
705  // Check all tracks
706  for(Int_t i = 0; i < mGammaEvent->numberOfTracks(); i++)
707  {
708 
709  StGammaTrack *track = mGammaEvent->track(i);
710  if (!track) continue;
711 
712  Float_t trackEta = track->eta();
713  Float_t trackPhi = track->phi();
714  Float_t dEta = trackEta - candidateEta;
715  Float_t dPhi = TVector2::Phi_mpi_pi(trackPhi - candidatePhi);
716 
717  Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
718 
719  if(r <= mRadius)
720  {
721  mCandidate->addTrack(track);
722  if(!track->candidates.FindObject(mCandidate)) track->candidates.Add(mCandidate);
723  }
724 
725  }
726 
727  // Check all towers
728  for(Int_t i = 0; i < mGammaEvent->numberOfTowers(); i++)
729  {
730 
731  StGammaTower *tower = mGammaEvent->tower(i);
732  if(!tower) continue;
733 
734  Float_t towerEta = tower->eta;
735  Float_t towerPhi = tower->phi;
736  Float_t dEta = towerEta - candidateEta;
737  Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
738 
739  Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
740 
741  if(r <= mRadius)
742  {
743  mCandidate->addTower(tower);
744  if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
745  }
746 
747  }
748 
749  // Check all preshowers
750  for(Int_t i = 0; i < mGammaEvent->numberOfPreshower1(); i++)
751  {
752 
753  StGammaTower *tower = mGammaEvent->preshower1(i);
754  if(!tower) continue;
755 
756  Float_t towerEta = tower->eta;
757  Float_t towerPhi = tower->phi;
758  Float_t dEta = towerEta - candidateEta;
759  Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
760 
761  Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
762 
763  if(r <= mRadius)
764  {
765  mCandidate->addPreshower1(tower);
766  if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
767  }
768 
769  }
770 
771 
772  for(Int_t i = 0; i < mGammaEvent->numberOfPreshower2(); i++)
773  {
774 
775  StGammaTower *tower = mGammaEvent->preshower2(i);
776  if(!tower) continue;
777 
778  Float_t towerEta = tower->eta;
779  Float_t towerPhi = tower->phi;
780  Float_t dEta = towerEta - candidateEta;
781  Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
782 
783  Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
784 
785  if(r <= mRadius)
786  {
787  mCandidate->addPreshower2(tower);
788  if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
789  }
790 
791  }
792 
793  // Check all postshowers
794  for(Int_t i = 0; i < mGammaEvent->numberOfPostshower(); i++)
795  {
796 
797  StGammaTower *tower = mGammaEvent->postshower(i);
798  if (!tower) continue;
799 
800  Float_t towerEta = tower->eta;
801  Float_t towerPhi = tower->phi;
802  Float_t dEta = towerEta - candidateEta;
803  Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
804 
805  Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
806 
807  if(r <= mRadius)
808  {
809  mCandidate->addPostshower(tower);
810  if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
811  }
812 
813  }
814 
815  // Now associate SMD strips within mSMDRange of the cluster, being
816  // care in how one defines the center from which mSMDRange is defined
817 
818  Float_t umin[kEEmcNumSectors], vmin[kEEmcNumSectors];
819  Float_t umax[kEEmcNumSectors], vmax[kEEmcNumSectors];
820  Float_t umid[kEEmcNumSectors], vmid[kEEmcNumSectors];
821  Int_t ntow[kEEmcNumSectors];
822 
823  // First, find the geometric center of the tower cluster in each sector
824  for(Int_t ii = 0; ii < 12; ii++)
825  {
826  umin[ii] = 0.0;
827  vmin[ii] = 0.0;
828  umax[ii] = 0.0;
829  vmax[ii] = 0.0;
830  umid[ii] = -1.0;
831  vmid[ii] = -1.0;
832  ntow[ii] = 0;
833  }
834 
835  EEmcSmdMap *eemap = EEmcSmdMap::instance();
836  for(Int_t itow = 0; itow < cluster.numberOfTowers(); itow++)
837  {
838 
839  StEEmcTower tower = cluster.tower(itow);
840  Int_t U, V;
841  eemap->getMiddleU( tower.sector(), tower.subsector(), tower.etabin(), U );
842  eemap->getMiddleV( tower.sector(), tower.subsector(), tower.etabin(), V );
843  ntow[ tower.sector() ]++;
844  umid[ tower.sector() ]+=0.5+(Float_t)U;// from middle of strip
845  vmid[ tower.sector() ]+=0.5+(Float_t)V;
846 
847  }
848 
849 
850  for(Int_t isec = 0; isec < 12; isec++)
851  {
852 
853  if(ntow[isec])
854  {
855  umid[isec] /= ntow[isec];
856  vmid[isec] /= ntow[isec];
857  umin[isec] = TMath::Max(umid[isec] - mEsmdRange * 2.0, 0. );
858  vmin[isec] = TMath::Max(vmid[isec] - mEsmdRange * 2.0, 0. );
859  umax[isec] = TMath::Min(umid[isec] + mEsmdRange * 2.0, 287. );
860  vmax[isec] = TMath::Min(vmid[isec] + mEsmdRange * 2.0, 287. );
861  }
862 
863  }
864 
865  // Now associate SMD strips within mSMDRange of the center
866  for(Int_t isec = 0; isec < 12; isec++)
867  {
868 
869  if(!ntow[isec]) continue;
870 
871  // Start with the U plane
872  for(Int_t i = (Int_t)umin[isec]; i < (Int_t)umax[isec]; i++)
873  {
874 
875  StGammaStrip *strip = mGammaRawMaker->strip(isec, 0, i);
876  if(strip)
877  {
878  strip->position = i;
879  mCandidate->addSmdu(strip);
880  strip->candidates.Add(mCandidate);
881  }
882 
883  }
884 
885  // Finish with the V plane
886  for(Int_t i = (Int_t)vmin[isec]; i < (Int_t)vmax[isec]; i++)
887  {
888 
889  StGammaStrip *strip = mGammaRawMaker->strip(isec,1,i);
890  if(strip)
891  {
892  strip->position = i;
893  mCandidate->addSmdv(strip);
894  strip->candidates.Add(mCandidate);
895  }
896 
897  }
898 
899  } // Sectors (SMD)
900 
901  // Sum SMDU energy
902  float smduEnergy = 0;
903 
904  for(int i = 0; i < mCandidate->numberOfSmdu(); ++i)
905  {
906  StGammaStrip* strip = mCandidate->smdu(i);
907  smduEnergy += strip->energy;
908  }
909 
910  mCandidate->SetSmduEnergy(smduEnergy);
911 
912  // Sum SMDV energy
913  float smdvEnergy = 0;
914  for (int i = 0; i < mCandidate->numberOfSmdv(); ++i)
915  {
916  StGammaStrip* strip = mCandidate->smdv(i);
917  smdvEnergy += strip->energy;
918  }
919 
920  mCandidate->SetSmdvEnergy(smdvEnergy);
921 
922  } // Clusters
923 
924  } // Sectors (Clusters)
925 
926  // Run the gamma fitter
927  for(int i = 0; i < mGammaEvent->numberOfCandidates(); ++i)
928  {
929 
930  StGammaCandidate* candidate = mGammaEvent->candidate(i);
931  if (candidate->detectorId() == StGammaCandidate::kEEmc)
932  {
933  StGammaFitterResult ufit;
934  StGammaFitterResult vfit;
935  int ustatus = StGammaFitter::instance()->fit(candidate, &ufit, 0);
936  if(ustatus == 0) candidate->SetSmdFit(ufit,0);
937  int vstatus = StGammaFitter::instance()->fit(candidate, &vfit, 1);
938  if(vstatus == 0) candidate->SetSmdFit(vfit,1);
939  }
940 
941  }
942 
943  return kStOK;
944 
945 }
946 
948 // Return the position of an EEMC Cluster //
950 TVector3 StGammaCandidateMaker::getEEmcClusterPosition(const StEEmcCluster& cluster)
951 {
952 
953  // Get gamma raw maker
954  StGammaRawMaker *mGammaRawMaker = dynamic_cast<StGammaRawMaker*>(GetMakerInheritsFrom("StGammaRawMaker"));
955  if(!mGammaRawMaker)
956  {
957  LOG_WARN << "MakeEndcap() - No StGammaRawMaker found!" << endm;
958  return TVector3(-999,-999,-999);
959  }
960 
961  // Get cluster seed tower
962  StEEmcTower seed = cluster.tower(0);
963 
964  // Get ranges of ESMD strips under seed tower
965  int sector = seed.sector();
966  int subsector = seed.subsector();
967  int etabin = seed.etabin();
968  int xmin[2], xmax[2];
969 
970  EEmcSmdMap::instance()->getRangeU(sector, subsector, etabin, xmin[0], xmax[0]);
971  EEmcSmdMap::instance()->getRangeV(sector, subsector, etabin, xmin[1], xmax[1]);
972 
973  // Find max strips within ranges
974  StGammaStrip* maxStrips[2] = {0, 0};
975  for(int plane = 0; plane < 2; ++plane)
976  {
977 
978  for(int i = xmin[plane]; i <= xmax[plane]; ++i)
979  {
980 
981  StGammaStrip* strip = mGammaRawMaker->strip(sector, plane, i);
982  if(strip)
983  {
984  if(!maxStrips[plane] || strip->energy > maxStrips[plane]->energy)
985  {
986  maxStrips[plane] = strip;
987  }
988  }
989 
990  }
991 
992  }
993 
994  // Get intersection of max strips and check that it lies within
995  // the fiducial volume of the seed tower.
996  if(maxStrips[0] && maxStrips[1])
997  {
998 
999  TVector3 position = EEmcSmdGeom::instance()->getIntersection(sector, maxStrips[0]->index, maxStrips[1]->index);
1000  int sec, sub, eta;
1001 
1002  bool success = position.z() != -999;
1003  success &= EEmcGeomSimple::Instance().getTower(position, sec, sub, eta);
1004  success &= sector == sec;
1005  success &= subsector == sub;
1006  success &= etabin == eta;
1007 
1008  if(success) return position;
1009 
1010  }
1011 
1012  // Failed to calculate cluster position using ESMD strips
1013  return TVector3(0, 0, 0);
1014 
1015 }
1016 
1018 // Clean up the StGammaEvent, dropping detector //
1019 // objects not associated with a candidate //
1021 Int_t StGammaCandidateMaker::Compress()
1022 {
1023 
1024  // No compression
1025  if(mCompressLevel == 0) return kStOk;
1026 
1027  // Retrieve the StGammaEvent via StGammaEventMaker
1028  StGammaEventMaker *mGammaEventMaker = dynamic_cast<StGammaEventMaker*>(GetMakerInheritsFrom("StGammaEventMaker"));
1029  if(!mGammaEventMaker)
1030  {
1031  LOG_WARN << "Compress() - No StGammaEventMaker found!" << endm;
1032  return kStWarn;
1033  }
1034 
1035  StGammaEvent *mGammaEvent = mGammaEventMaker->event();
1036  if(!mGammaEvent)
1037  {
1038  LOG_WARN << "Compress() - No StGammaEvent found!" << endm;
1039  return kStWarn;
1040  }
1041 
1042  Compress<StGammaStrip>(mGammaEvent->mStrips);
1043 
1044  // Compress SMD strips only
1045  if(mCompressLevel == 1) return kStOk;
1046 
1047  // Drop strips, tracks, towers without candidates
1048  Compress<StGammaTrack>(mGammaEvent->mTracks);
1049  Compress<StGammaTower>(mGammaEvent->mTowers);
1050  Compress<StGammaTower>(mGammaEvent->mPreshower1);
1051  Compress<StGammaTower>(mGammaEvent->mPreshower2);
1052  Compress<StGammaTower>(mGammaEvent->mPostshower);
1053 
1054  // Compress all
1055  if(mCompressLevel == 2) return kStOk;
1056 
1057  LOG_WARN << Form("Unknown compression level (%d)", mCompressLevel) << endm;
1058 
1059  return kStWarn;
1060 
1061 }
1062 
1064 // Template function for dropping objects in //
1065 // the TClonesArray without an association //
1066 // to a StGammaCandidate //
1068 template<class T>
1069 void StGammaCandidateMaker::Compress(TClonesArray* clones)
1070 {
1071  TIter next(clones);
1072  while (T* x = (T*)next()) if (x->candidates.IsEmpty()) clones->Remove(x);
1073  clones->Compress();
1074 }
StEEmcCluster::momentum
TVector3 momentum() const
Definition: StEEmcCluster.h:69
StEEmcCluster::tower
StEEmcTower & tower(Int_t t)
Get the specified tower within the cluster.
Definition: StEEmcCluster.h:79
EEmcGeomSimple::getZ1
Float_t getZ1() const
gets lower Z edge of EEMC (preshower)
Definition: EEmcGeomSimple.h:73
StEEmcA2EMaker
EEmc ADC –&gt; energy maker.
Definition: StEEmcA2EMaker.h:19
EEmcGeomSimple::Instance
static EEmcGeomSimple & Instance()
returns a reference to a static instance of EEmcGeomSimple
Definition: EEmcGeomSimple.h:178
StGammaCandidateMaker::Clear
void Clear(Option_t *opts="")
User defined functions.
Definition: StGammaCandidateMaker.cxx:58
StGammaEvent::tower
StGammaTower * tower(Int_t i) const
Return ith track.
Definition: StGammaEvent.h:66
StGammaEvent::preshower1
StGammaTower * preshower1(Int_t i) const
Return ith tower.
Definition: StGammaEvent.h:67
StGammaEvent::numberOfPostshower
Int_t numberOfPostshower() const
Return number of pre2.
Definition: StGammaEvent.h:61
StEEmcCluster::position
TVector3 position() const
Definition: StEEmcCluster.h:73
StEmcPosition::getNextId
Int_t getNextId(Int_t det, Int_t m, Int_t e, Int_t s, Int_t nEta, Int_t nPhi) const
Return neighbor id (works for all detectors 1=bemc, 2=bprs, 3=bsmde, 4=bsmdp)
Definition: StEmcPosition.cxx:384
StGammaEvent::numberOfPreshower1
Int_t numberOfPreshower1() const
Return number of towers.
Definition: StGammaEvent.h:59
StGammaEvent::numberOfCandidates
Int_t numberOfCandidates() const
Return number of strips.
Definition: StGammaEvent.h:63
StEEmcGenericClusterMaker::clusters
StEEmcClusterVec_t & clusters(Int_t sec, Int_t layer)
Return a vector of tower clusters.
Definition: StEEmcGenericClusterMaker.h:129
StGammaFitter::instance
static StGammaFitter * instance()
Access to single instance of this singleton class.
Definition: StGammaFitter.cxx:67
StEEmcTower::etabin
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
Definition: StEEmcTower.h:45
StEEmcElement::fail
void fail(unsigned f)
Set a fail bit for this element.
Definition: StEEmcElement.h:25
StGammaEvent::vertex
TVector3 & vertex()
Returns muDst file from which event originated.
Definition: StGammaEvent.h:79
StEEmcTower::subsector
Int_t subsector() const
Returns subsector of this tower, pre- or postshower element.
Definition: StEEmcTower.h:43
StGammaEvent::newCandidate
StGammaCandidate * newCandidate()
Add a new SMD strip.
Definition: StGammaEvent.cxx:153
StEEmcTower::index
void index(Int_t i)
Definition: StEEmcTower.cxx:76
StEEmcTower
Base class for representing tower, preshower and postshower elements.
Definition: StEEmcTower.h:11
StGammaEvent::newStrip
StGammaStrip * newStrip()
Add a new postshower element.
Definition: StGammaEvent.cxx:146
StEEmcTower::sector
Int_t sector() const
Returns sector of this tower, pre- or postshower element.
Definition: StEEmcTower.h:41
StGammaFitter::fit
int fit(StGammaCandidate *candidate, StGammaFitterResult *fits, Int_t plane=0)
Fit transverse SMD profile to predetermined peak in u- and v-plane.
Definition: StGammaFitter.cxx:73
StGammaEvent::track
StGammaTrack * track(Int_t i) const
Return number of candidates.
Definition: StGammaEvent.h:65
StGammaTower::subsector
Int_t subsector()
Definition: StGammaTower.cxx:22
StEEmcCluster::numberOfTowers
Int_t numberOfTowers() const
Get the number of towers in cluster.
Definition: StEEmcCluster.h:76
kStWarn
Definition: Stypes.h:42
EEmcGeomSimple
EEMC simple geometry.
Definition: EEmcGeomSimple.h:22
kStOK
Definition: Stypes.h:40
StGammaTrack::phi
Float_t phi() const
eta at vertex
Definition: StGammaTrack.h:65
StEEmcCluster::energy
Float_t energy() const
Get energy of this cluster.
Definition: StEEmcCluster.h:62
EEmcSmdGeom::getIntersection
TVector3 getIntersection(Int_t iSec, Int_t iUStrip, Int_t iVStrip, const TVector3 &vertex) const
Definition: EEmcSmdGeom.cxx:628
StGammaTrack::candidates
TRefArray candidates
Returns (0,0,0) if failed.
Definition: StGammaTrack.h:49
StEEmcCluster
A base class for describing clusters of EEMC towers.
Definition: StEEmcCluster.h:50
StGammaEvent::numberOfPreshower2
Int_t numberOfPreshower2() const
Return number of pre1.
Definition: StGammaEvent.h:60
StEEmcA2EMaker::tower
StEEmcTower & tower(Int_t index, Int_t layer=0)
Definition: StEEmcA2EMaker.h:69
StGammaEvent::preshower2
StGammaTower * preshower2(Int_t i) const
Return ith pre1.
Definition: StGammaEvent.h:68
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
StGammaTrack::pz
Float_t pz() const
pt at vertex
Definition: StGammaTrack.h:63
StGammaTrack::positionAtRadius
TVector3 positionAtRadius(Double_t radius) const
Returns outer helix (last measured point)
Definition: StGammaTrack.cxx:48
EEmcGeomSimple::getTower
bool getTower(const TVector3 &r, int &sec, int &sub, int &etabin, Float_t &dphi, Float_t &deta) const
Definition: EEmcGeomSimple.cxx:159
StGammaTrack::positionAtZ
TVector3 positionAtZ(Double_t z) const
Returns (0,0,0) if failed.
Definition: StGammaTrack.cxx:85
StGammaEvent::postshower
StGammaTower * postshower(Int_t i) const
Return ith pre2.
Definition: StGammaEvent.h:69
StEEmcElement::energy
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Definition: StEEmcElement.h:21
StGammaEvent::numberOfTowers
Int_t numberOfTowers() const
Return number of tracks.
Definition: StGammaEvent.h:58
StGammaEvent::candidate
StGammaCandidate * candidate(Int_t i) const
Return ith strip.
Definition: StGammaEvent.h:71
kStOk
Definition: Stypes.h:41
StGammaTrack::eta
Float_t eta() const
pz at vertex
Definition: StGammaTrack.h:64
Generated by   doxygen 1.8.5