StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcAssociationMaker.cxx
1 /********************************************************************************
2  *
3  * Author: Marcia Maria de Moura
4  *
5  * See README file for description
6  *
7  ********************************************************************************/
8 
9 #include "StEEmcAssociationMaker.h"
10 
11 #include <Stiostream.h>
12 #include <math.h>
13 
14 #include "StEmcUtil/geometry/StEmcGeom.h"
15 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
16 #include "StMcEventTypes.hh"
17 #include "StMcEvent.hh"
18 #include "StEventTypes.h"
19 #include <StMessMgr.h>
20 
21 const TString eemcDetname[] = {"etow","eprs","esmdu","esmdv"};
22 
23 ClassImp(StEEmcAssociation)
24 ClassImp(StEEmcClusterAssociation)
25 ClassImp(StEEmcPointAssociation)
26 ClassImp(StEEmcAssociationMaker)
27 //------------------------------------------------------------------------------
28 StEEmcAssociation::StEEmcAssociation(StMcTrack *t)
29 { mTrack = t; }
30 StEEmcAssociation::~StEEmcAssociation()
31 { }
32 //------------------------------------------------------------------------------
33 StEEmcClusterAssociation::StEEmcClusterAssociation(StMcTrack* t, StEmcCluster* c, float ft,float fe):StEEmcAssociation(t)
34 { mCluster=c; mFTrack = ft; mFEmc = fe;}
35 StEEmcClusterAssociation::~StEEmcClusterAssociation()
36 { }
37 //------------------------------------------------------------------------------
38 StEEmcPointAssociation::StEEmcPointAssociation(StMcTrack* t, StEmcPoint* p, int at):StEEmcAssociation(t)
39 { mPoint=p; mAssocType = at; }
40 StEEmcPointAssociation::~StEEmcPointAssociation()
41 { }
42 //------------------------------------------------------------------------------
43 StEEmcAssociationMaker::StEEmcAssociationMaker(const char* name):StMaker(name)
44 {
45  for(Int_t i=0;i<NDETECTORS;i++)
46  {
47  mTrackCluster[i] = NULL;
48  mClusterTrack[i] = NULL;
49  }
50  mTrackPoint = NULL;
51  mPointTrack = NULL;
52 
53  mPrint = kTRUE;
54 }
55 //------------------------------------------------------------------------------
56 StEEmcAssociationMaker::~StEEmcAssociationMaker()
57 {
58 }
59 void StEEmcAssociationMaker::Clear(const char* a)
60 {
61  if(mPrint) cout <<"Cleaning old stuff from EEMC association\n";
62  for(Int_t i=0;i<NDETECTORS;i++)
63  {
64  if(mPrint) cout <<"Cleaning for detector "<<i<<endl;
65  if(mTrackCluster[i])
66  {
67  for(multiEEmcTrackClusterIter j=mTrackCluster[i]->begin(); j!=mTrackCluster[i]->end(); j++)
68  {
69  SafeDelete((*j).second);
70  }
71  mTrackCluster[i]->clear();
72  SafeDelete(mTrackCluster[i]);
73  mTrackCluster[i]=NULL;
74  }
75  if(mClusterTrack[i])
76  {
77  mClusterTrack[i]->clear();
78  SafeDelete(mClusterTrack[i]);
79  mClusterTrack[i]=NULL;
80  }
81  }
82  if (mPrint) cout <<"Cleaning points Association\n";
83  if(mTrackPoint)
84  {
85  for(multiEEmcTrackPointIter j=mTrackPoint->begin(); j!=mTrackPoint->end(); j++)
86  {
87  SafeDelete((*j).second);
88  }
89  mTrackPoint->clear();
90  SafeDelete(mTrackPoint);
91  mTrackPoint=NULL;
92  }
93  if(mPointTrack)
94  {
95  mPointTrack->clear();
96  SafeDelete(mPointTrack);
97  mPointTrack=NULL;
98  }
99  return;
100 }
101 //------------------------------------------------------------------------------
102 Int_t StEEmcAssociationMaker::Init()
103 {
104  return StMaker::Init();
105 }
106 //------------------------------------------------------------------------------
107 Int_t StEEmcAssociationMaker::Make()
108 {
109  if(mPrint) gMessMgr->Info() << "StEEmcAssociationMaker::Make()" << endm;
110 
111  // Getting McEvent object
112  StMcEvent* mcEvent= (StMcEvent*) GetDataSet("StMcEvent");
113  if (!mcEvent) return kStWarn;
114  StSPtrVecMcTrack& tracks=mcEvent->tracks();
115  if(tracks.size()==0) return kStWarn;
116 
117  // Getting Event object
118  StEvent* event=(StEvent*) GetDataSet("StEvent");
119  if (!event) return kStWarn;
120 
121  if(Debug()==2) if(mPrint) printHits(event);
122 
123  StEmcCollection* emcCollection=event->emcCollection();
124  if(!emcCollection) return kStWarn;
125 
126  // Starting doing cluster association for each detector
127  for (Int_t detnum=0; detnum<NDETECTORS; detnum++) // For detnum<4, only endcap EMC
128  {
129  if(mPrint) cout <<"Doing association for detector "<<detnum<<endl;
130  StDetectorId detId=static_cast<StDetectorId>(detnum+kEndcapEmcTowerId);
131  StEmcDetector* detector=emcCollection->detector(detId);
132  StEmcClusterCollection* clusterColl=NULL;
133  if (detector) clusterColl=detector->cluster();
134  if (clusterColl)
135  {
136  StSPtrVecEmcCluster& clusters=clusterColl->clusters();
137 
138  // Allocate cluster matching maps
139  if(mPrint) cout<<"Track size = "<<tracks.size()<<" cluster size = "<<clusters.size()<<endl;
140  if(!mTrackCluster[detnum]) mTrackCluster[detnum]=new multiEEmcTrackCluster;
141  if(!mClusterTrack[detnum]) mClusterTrack[detnum]=new multiEEmcClusterTrack;
142 
143  for (UInt_t i=0; i<tracks.size(); i++)
144  {
145  StPtrVecMcCalorimeterHit hits;
146  switch (detnum)
147  {
148  case 0: { hits=tracks[i]->eemcHits(); break; }
149  case 1: { hits=tracks[i]->eprsHits(); break; }
150  case 2: { hits=tracks[i]->esmduHits(); break; }
151  case 3: { hits=tracks[i]->esmdvHits(); break; }
152  default: { gMessMgr->Warning()<<"Detector does not exist"<<endm; break; }
153  }
154 
155  if (hits.size()!=0);
156  {
157  // Calculating track energy by summing over hits energies
158  Float_t energy=0;
159  for (UInt_t i2=0; i2<hits.size(); i2++)
160  {
161  Float_t hitEnergy=hits[i2]->dE();
162  energy+=hitEnergy;
163  }
164 
165  for (UInt_t j=0; j<clusters.size(); j++)
166  {
167  StPtrVecEmcRawHit& clHits=clusters[j]->hit();
168  Float_t trackEnergyFraction=0;
169  Float_t clusterEnergyFraction=0;
170  Int_t assoc=0;
171  for (UInt_t k=0; k<hits.size(); k++)
172  {
173  UInt_t module=hits[k]->module();
174  UInt_t eta=hits[k]->eta();
175  Int_t sub=hits[k]->sub();
176 
177  for (UInt_t l=0; l<clHits.size(); l++)
178  {
179  UInt_t clModule=clHits[l]->module();
180  UInt_t clEta=clHits[l]->eta();
181  Int_t clSub=abs(clHits[l]->sub());
182 
183  //compare between hit track and hit cluster
184  if (module==clModule && eta==clEta
185  && ((detnum > 1)? 1 : sub==clSub))
186  {
187  assoc=1;
188  if(Debug()) if(mPrint) {
189  cout << eemcDetname[detnum]
190  << " cluster = " << j
191  << ", m = " <<module <<" " << clModule
192  << ", e = " << eta << " " << clEta
193  << ", s = " << sub << " " << clSub << endl;
194  }
195  Float_t hitEnergy=hits[k]->dE();
196  trackEnergyFraction+=hitEnergy/energy;
197  Float_t clEnergy=clusters[j]->energy();
198  clusterEnergyFraction+=hitEnergy/clEnergy;
199  }
200  }
201  }
202  if (assoc) {
203  StEEmcClusterAssociation *c=new StEEmcClusterAssociation(tracks[i],clusters[j],trackEnergyFraction,clusterEnergyFraction);
204  mTrackCluster[detnum]->insert(multiEEmcTrackClusterValue(tracks[i],c));
205 
206  mClusterTrack[detnum]->insert(multiEEmcClusterTrackValue(clusters[j],c));
207  }
208  }
209  }
210  }
211  }
212  }
213  // Finishing doing cluster association
214 
215  // Starting doing point association
216  if(mPrint) cout <<"Doing point association\n";
217  StSPtrVecEmcPoint& endcapPoints=emcCollection->endcapPoints();
218  cout <<"points.size() = " << endcapPoints.size() << endl;
219  for (UInt_t j=0; j<endcapPoints.size(); j++) {
220  int nCl[4];
221  for(Int_t detnum=0;detnum<NDETECTORS;detnum++){
222  StDetectorId detId=static_cast<StDetectorId>(detnum+kEndcapEmcTowerId);
223  StPtrVecEmcCluster& cluster=endcapPoints[j]->cluster(detId);
224  nCl[detnum] = cluster.size();
225  }
226  if(nCl[0]!=0 && nCl[1]!=0 && nCl[2]!=1 && nCl[3]!=1 &&
227  endcapPoints[j]->nTracks()!=0) endcapPoints[j]->print();
228  }
229  if(endcapPoints.size()!=0)
230  {
231  // Allocate point matching maps
232  if(!mTrackPoint) mTrackPoint = new multiEEmcTrackPoint;
233  if(!mPointTrack) mPointTrack = new multiEEmcPointTrack;
234 
235  for (UInt_t i=0; i<tracks.size(); i++)
236  {
237  StPtrVecMcCalorimeterHit hits;
238 
239  for (UInt_t j=0; j<endcapPoints.size(); j++)
240  {
241  Int_t assoc=0;
242  for(Int_t detnum=0;detnum<NDETECTORS;detnum++)
243  {
244  switch (detnum)
245  {
246  case 0: { hits=tracks[i]->eemcHits(); break; }
247  case 1: { hits=tracks[i]->eprsHits(); break; }
248  case 2: { hits=tracks[i]->esmduHits(); break; }
249  case 3: { hits=tracks[i]->esmdvHits(); break; }
250  default: { gMessMgr->Warning()<<"Detector does not exist"<<endm; break; }
251  }
252  StDetectorId detId=static_cast<StDetectorId>(detnum+kEndcapEmcTowerId);
253  StPtrVecEmcCluster& cluster=endcapPoints[j]->cluster(detId);
254  for (UInt_t k=0; k<hits.size(); k++) if(hits[k])
255  {
256  UInt_t module=hits[k]->module();
257  UInt_t eta=hits[k]->eta();
258  Int_t sub=hits[k]->sub();
259 
260  for (UInt_t j2=0; j2<cluster.size(); j2++) if(cluster[j2])
261  {
262  StPtrVecEmcRawHit& clHit=cluster[j2]->hit();
263  for (UInt_t l=0; l<clHit.size(); l++) if(clHit[l])
264  {
265  UInt_t clModule=clHit[l]->module();
266  UInt_t clEta=clHit[l]->eta();
267  Int_t clSub=abs(clHit[l]->sub());
268 
269  // compare between hit track and hit cluster from a point
270  if (module==clModule && eta==clEta
271  && ((detnum > 1)? 1 : sub==clSub))
272  {
273  assoc+=(1<<(detnum));
274  if(Debug()) if(mPrint) {
275  cout <<"assoc = " << assoc << " "
276  << eemcDetname[detnum]
277  << ", m = " <<module <<" " << clModule
278  << ", e = " << eta << " " << clEta
279  << ", s = " << sub << " " << clSub << endl;
280  }
281  goto nextDetector;
282  }
283  }
284  }
285  }
286  nextDetector: continue;
287  }
288  if (assoc) {
289  StEEmcPointAssociation *p= new StEEmcPointAssociation(tracks[i],endcapPoints[j],assoc);
290  mTrackPoint->insert(multiEEmcTrackPointValue(tracks[i],p));
291  mPointTrack->insert(multiEEmcPointTrackValue(endcapPoints[j],p));
292  }
293  }
294  }
295  }
296  if(Debug()) if(mPrint) printMaps();
297  // Finishing doing point association
298  return kStOK;
299 }
300 //------------------------------------------------------------------------------
301 Int_t StEEmcAssociationMaker::Finish()
302 {
303  // Finishing the maker
304  gMessMgr->Info() << "StEEmcAssociationMaker::Finish()" << endm;
305  return kStOK;
306 }
307 //---------------------------------------------------------------------
308 multiEEmcTrackCluster* StEEmcAssociationMaker::getTrackClusterMap(const char* detname)
309 {
310  Int_t det = getDetNum(detname);
311  if(det>0) return mTrackCluster[det-1];
312  return NULL;
313 }
314 //---------------------------------------------------------------------
315 multiEEmcClusterTrack* StEEmcAssociationMaker::getClusterTrackMap(const char* detname)
316 {
317  Int_t det = getDetNum(detname);
318  if(det>0) return mClusterTrack[det-1];
319  return NULL;
320 }
321 //---------------------------------------------------------------------
322 Int_t StEEmcAssociationMaker::getDetNum(const char* detname)
323 {
324  Int_t detnum = 0;
325  for (Int_t i=0; i<NDETECTORS; i++)
326  if (!strcmp(detname,eemcDetname[i])) detnum = i+1;
327  return detnum;
328 }
329 //---------------------------------------------------------------------
330 void StEEmcAssociationMaker::printMaps()
331 {
332  for(int i=0;i<NDETECTORS;i++)
333  {
334  cout <<"-----------------------------------------------------------\n";
335  cout <<"ASSOCIATION FOR DETECTOR "<<eemcDetname[i].Data()<<endl;
336  multiEEmcTrackCluster* map = getTrackClusterMap(i+1);
337  if(map)
338  {
339  cout <<"Track->Cluster Association Map\n";
340  for(multiEEmcTrackClusterIter j=map->begin(); j!=map->end(); j++)
341  {
342  StMcTrack* track = (StMcTrack*)(*j).first;
343  StEEmcClusterAssociation* value = (StEEmcClusterAssociation*)(*j).second;
344  if(track && value)
345  {
346  StEmcCluster *c = (StEmcCluster*) value->getCluster();
347 
348  if(c)
349  {
350  cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
351  <<" Cl = "<<c<<" E = "<<c->energy()<<" eta = "<<c->eta()<<" phi = "<<c->phi()
352  <<" FrTr = "<<value->getFractionTrack()<<" FrCl = "<<value->getFractionCluster()<<endl;
353  }
354  }
355  }
356  }
357 
358  multiEEmcClusterTrack* map1 = getClusterTrackMap(i+1);
359  if(map1)
360  {
361  cout <<"Cluster->Track Association Map\n";
362  for(multiEEmcClusterTrackIter j=map1->begin(); j!=map1->end(); j++)
363  {
364  StEmcCluster *c = (StEmcCluster*)(*j).first;
365  StEEmcClusterAssociation* value = (StEEmcClusterAssociation*)(*j).second;
366  if(c && value)
367  {
368  StMcTrack* track = (StMcTrack*)value->getTrack();
369  if(track)
370  {
371  cout <<" Cl = "<<c<<" E = "<<c->energy()<<" eta = "<<c->eta()<<" phi = "<<c->phi()
372  <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
373  <<" FrTr = "<<value->getFractionTrack()<<" FrCl = "<<value->getFractionCluster()<<endl;
374  }
375  }
376  }
377  }
378  }
379  multiEEmcTrackPoint *mapPoint = getTrackPointMap();
380  if(mapPoint)
381  {
382  cout <<"Track->Point Association Map\n";
383  for(multiEEmcTrackPointIter j = mapPoint->begin(); j!=mapPoint->end(); j++)
384  {
385  StMcTrack *track = (StMcTrack*)(*j).first;
386  StEEmcPointAssociation *value = (StEEmcPointAssociation*)(*j).second;
387  if(value)
388  {
389  StEmcPoint *point = (StEmcPoint*)value->getPoint();
390  if(track && point)
391  {
392  cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
393  <<" Point = "<<point<<" E = "<<point->energy()
394  <<" Assoc. = "<<value->getAssociation();
395  for(int i=1;i<=4;i++) cout <<" det "<<i<<" A = "<< value->getAssociation(i);
396  cout <<endl;
397  }
398  }
399  }
400  }
401  multiEEmcPointTrack *mapPoint1 = getPointTrackMap();
402  if(mapPoint1)
403  {
404  cout <<"Point->Track Association Map\n";
405  for(multiEEmcPointTrackIter j = mapPoint1->begin(); j!=mapPoint1->end(); j++)
406  {
407  StEmcPoint *point = (StEmcPoint*)(*j).first;
408  StEEmcPointAssociation *value = (StEEmcPointAssociation*)(*j).second;
409  if(value)
410  {
411  StMcTrack *track = (StMcTrack*)value->getTrack();
412  if(track && point)
413  {
414  cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
415  <<" Point = "<<point<<" E = "<<point->energy()
416  <<" Assoc. = "<<value->getAssociation();
417  for(int i=1;i<=4;i++) cout <<" det "<<i<<" A = "<< value->getAssociation(i);
418  cout <<endl;
419  }
420  }
421  }
422  }
423 }
424 //-------------------------------------------------------------------
425 void StEEmcAssociationMaker::printHits(StEvent *event)
426 {
427  StEmcCollection* emccol=(StEmcCollection*)event->emcCollection();
428 
429  gMessMgr->Info()<<"****** EEMC Association printHits() *******" <<endm;
430  for(Int_t i=0; i<NDETECTORS; i++)
431  {
432  StDetectorId id = static_cast<StDetectorId>(i+kEndcapEmcTowerId);
433  StEmcDetector* detector=emccol->detector(id);
434  gMessMgr->Info()<<"****************** hits in detector "
435  << eemcDetname[i].Data()<<endm;
436  if(detector) for(int j=1;j<=kEEmcNumSectors;j++)
437  {
438  StEmcModule* sector = detector->module(j);
439  StSPtrVecEmcRawHit& rawHit=sector->hits();
440  if(rawHit.size()>0)
441  gMessMgr->Info()<<"Number of hits for module "
442  <<j<<" = "<<rawHit.size()<<endm;
443  for(UInt_t k=0;k<rawHit.size();k++)
444  if(i < 2) // for Tower and Prs
445  gMessMgr->Info()<<"Hit number = "<<k
446  <<" sector = " << rawHit[k]->module()
447  <<" eta = "<<rawHit[k]->eta()
448  <<" sub = "<< rawHit[k]->sub()
449  <<" adc = "<< rawHit[k]->adc()
450  <<" energy = "<<rawHit[k]->energy()<<endm;
451  else
452  gMessMgr->Info()<<"Hit number = "<<k
453  <<" sector = " << rawHit[k]->module()
454  <<" strip = "<<rawHit[k]->eta()
455  <<" adc = "<< rawHit[k]->adc()
456  <<" energy = "<<rawHit[k]->energy()<<endm;
457  }
458  }
459 }
460 
462 //
463 // $Id: StEEmcAssociationMaker.cxx,v 1.2 2005/09/29 14:58:10 fisyak Exp $
464 // $Log: StEEmcAssociationMaker.cxx,v $
465 // Revision 1.2 2005/09/29 14:58:10 fisyak
466 // Persistent StMcEvent
467 //
468 // Revision 1.1.1.1 2005/05/31 18:54:13 wzhang
469 // First version
470 //
471 //
StEEmcAssociationMaker::getClusterTrackMap
multiEEmcClusterTrack * getClusterTrackMap(Int_t i)
returns multimap for association betwwen clusters and MC tracks
Definition: StEEmcAssociationMaker.h:218
StEEmcAssociationMaker::getTrackClusterMap
multiEEmcTrackCluster * getTrackClusterMap(Int_t i)
returns multimap for association betwwen MC tracks and clusters
Definition: StEEmcAssociationMaker.h:216
StEEmcAssociationMaker::printHits
void printHits(StEvent *)
Set print log.
Definition: StEEmcAssociationMaker.cxx:425
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StEEmcAssociationMaker::getDetNum
Int_t getDetNum(const char *)
returns detector number for each EMC sub detector
Definition: StEEmcAssociationMaker.cxx:322
StEEmcAssociationMaker::getPointTrackMap
multiEEmcPointTrack * getPointTrackMap()
returns multimap for association betwwen points and MC tracks
Definition: StEEmcAssociationMaker.h:221
StEEmcAssociationMaker::getTrackPointMap
multiEEmcTrackPoint * getTrackPointMap()
returns multimap for association betwwen MC tracks and points
Definition: StEEmcAssociationMaker.h:220
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
Generated by   doxygen 1.8.5