StRoot: StEEmcPool/StEEmcAssociationMaker/StEEmcAssociationMaker.cxx Source File

00001 /********************************************************************************
00002  *
00003  * Author: Marcia Maria de Moura
00004  *
00005  * See README file for description
00006  *
00007  ********************************************************************************/
00008 
00009 #include "StEEmcAssociationMaker.h"
00010 
00011 #include <Stiostream.h>
00012 #include <math.h>
00013 
00014 #include "StEmcUtil/geometry/StEmcGeom.h"
00015 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
00016 #include "StMcEventTypes.hh"
00017 #include "StMcEvent.hh"
00018 #include "StEventTypes.h"
00019 #include <StMessMgr.h>
00020 
00021 const TString eemcDetname[] = {"etow","eprs","esmdu","esmdv"};
00022 
00023 ClassImp(StEEmcAssociation)
00024 ClassImp(StEEmcClusterAssociation)
00025 ClassImp(StEEmcPointAssociation)
00026 ClassImp(StEEmcAssociationMaker)
00027 //------------------------------------------------------------------------------
00028 StEEmcAssociation::StEEmcAssociation(StMcTrack *t) 
00029 { mTrack = t; }
00030 StEEmcAssociation::~StEEmcAssociation() 
00031 { }
00032 //------------------------------------------------------------------------------
00033 StEEmcClusterAssociation::StEEmcClusterAssociation(StMcTrack* t, StEmcCluster* c, float ft,float fe):StEEmcAssociation(t) 
00034 { mCluster=c;  mFTrack = ft; mFEmc = fe;}
00035 StEEmcClusterAssociation::~StEEmcClusterAssociation()
00036 { }
00037 //------------------------------------------------------------------------------
00038 StEEmcPointAssociation::StEEmcPointAssociation(StMcTrack* t, StEmcPoint* p, int at):StEEmcAssociation(t) 
00039 { mPoint=p; mAssocType = at; }
00040 StEEmcPointAssociation::~StEEmcPointAssociation() 
00041 { }
00042 //------------------------------------------------------------------------------
00043 StEEmcAssociationMaker::StEEmcAssociationMaker(const char* name):StMaker(name)
00044 {
00045   for(Int_t i=0;i<NDETECTORS;i++)
00046   {
00047     mTrackCluster[i] = NULL;
00048     mClusterTrack[i] = NULL;
00049   }
00050   mTrackPoint = NULL;
00051   mPointTrack = NULL;
00052   
00053   mPrint = kTRUE;
00054 }
00055 //------------------------------------------------------------------------------
00056 StEEmcAssociationMaker::~StEEmcAssociationMaker()
00057 {
00058 }
00059 void StEEmcAssociationMaker::Clear(const char* a)
00060 {
00061   if(mPrint) cout <<"Cleaning old stuff from EEMC association\n";
00062   for(Int_t i=0;i<NDETECTORS;i++)
00063   {
00064     if(mPrint) cout <<"Cleaning for detector "<<i<<endl;
00065     if(mTrackCluster[i]) 
00066     { 
00067       for(multiEEmcTrackClusterIter j=mTrackCluster[i]->begin(); j!=mTrackCluster[i]->end(); j++) 
00068       {
00069         SafeDelete((*j).second);
00070       }
00071       mTrackCluster[i]->clear(); 
00072       SafeDelete(mTrackCluster[i]); 
00073       mTrackCluster[i]=NULL; 
00074     }
00075     if(mClusterTrack[i]) 
00076     { 
00077       mClusterTrack[i]->clear(); 
00078       SafeDelete(mClusterTrack[i]); 
00079       mClusterTrack[i]=NULL;
00080     }
00081   }
00082   if (mPrint) cout <<"Cleaning points Association\n";
00083   if(mTrackPoint) 
00084   { 
00085     for(multiEEmcTrackPointIter j=mTrackPoint->begin(); j!=mTrackPoint->end(); j++) 
00086     {
00087       SafeDelete((*j).second);
00088     }
00089     mTrackPoint->clear(); 
00090     SafeDelete(mTrackPoint); 
00091     mTrackPoint=NULL;
00092   }
00093   if(mPointTrack) 
00094   { 
00095     mPointTrack->clear(); 
00096     SafeDelete(mPointTrack); 
00097     mPointTrack=NULL;
00098   }
00099   return;
00100 }
00101 //------------------------------------------------------------------------------
00102 Int_t StEEmcAssociationMaker::Init()
00103 {
00104   return StMaker::Init();
00105 }
00106 //------------------------------------------------------------------------------
00107 Int_t StEEmcAssociationMaker::Make()
00108 {
00109   if(mPrint) gMessMgr->Info() << "StEEmcAssociationMaker::Make()" << endm;
00110   
00111   // Getting McEvent object
00112   StMcEvent* mcEvent= (StMcEvent*) GetDataSet("StMcEvent");
00113   if (!mcEvent) return kStWarn;
00114   StSPtrVecMcTrack& tracks=mcEvent->tracks();
00115   if(tracks.size()==0) return kStWarn;
00116 
00117   // Getting Event object
00118   StEvent* event=(StEvent*) GetDataSet("StEvent");
00119   if (!event) return kStWarn;
00120   
00121   if(Debug()==2) if(mPrint) printHits(event);
00122 
00123   StEmcCollection* emcCollection=event->emcCollection();
00124   if(!emcCollection) return kStWarn;
00125 
00126   // Starting doing cluster association for each detector
00127   for (Int_t detnum=0; detnum<NDETECTORS; detnum++) // For detnum<4, only endcap EMC
00128   {
00129     if(mPrint) cout <<"Doing association for detector "<<detnum<<endl;
00130     StDetectorId detId=static_cast<StDetectorId>(detnum+kEndcapEmcTowerId);
00131     StEmcDetector* detector=emcCollection->detector(detId);
00132     StEmcClusterCollection* clusterColl=NULL;
00133     if (detector) clusterColl=detector->cluster();
00134     if (clusterColl)
00135     { 
00136       StSPtrVecEmcCluster& clusters=clusterColl->clusters();
00137 
00138       // Allocate cluster matching maps
00139       if(mPrint) cout<<"Track size = "<<tracks.size()<<"  cluster size = "<<clusters.size()<<endl;
00140       if(!mTrackCluster[detnum]) mTrackCluster[detnum]=new multiEEmcTrackCluster;
00141       if(!mClusterTrack[detnum]) mClusterTrack[detnum]=new multiEEmcClusterTrack;
00142       
00143       for (UInt_t i=0; i<tracks.size(); i++)
00144       {
00145         StPtrVecMcCalorimeterHit hits;
00146         switch (detnum)
00147         {
00148           case 0: { hits=tracks[i]->eemcHits(); break; }
00149           case 1: { hits=tracks[i]->eprsHits(); break; }
00150           case 2: { hits=tracks[i]->esmduHits(); break; }
00151           case 3: { hits=tracks[i]->esmdvHits(); break; }
00152           default: { gMessMgr->Warning()<<"Detector does not exist"<<endm; break; }
00153         }
00154 
00155         if (hits.size()!=0);
00156         {
00157           // Calculating track energy by summing over hits energies
00158           Float_t energy=0;
00159           for (UInt_t i2=0; i2<hits.size(); i2++)
00160           {
00161             Float_t hitEnergy=hits[i2]->dE();
00162             energy+=hitEnergy;
00163           } 
00164           
00165           for (UInt_t j=0; j<clusters.size(); j++)
00166           {
00167             StPtrVecEmcRawHit& clHits=clusters[j]->hit();
00168             Float_t trackEnergyFraction=0;
00169             Float_t clusterEnergyFraction=0;
00170             Int_t assoc=0;
00171             for (UInt_t k=0; k<hits.size(); k++)
00172             {
00173               UInt_t module=hits[k]->module();
00174               UInt_t eta=hits[k]->eta();
00175               Int_t sub=hits[k]->sub();
00176               
00177               for (UInt_t l=0; l<clHits.size(); l++)
00178               {
00179                 UInt_t clModule=clHits[l]->module();
00180                 UInt_t clEta=clHits[l]->eta();
00181                 Int_t clSub=abs(clHits[l]->sub());
00182 
00183                 //compare between hit track and hit cluster
00184                 if (module==clModule && eta==clEta  
00185                     && ((detnum > 1)? 1 : sub==clSub))  
00186                 {
00187                   assoc=1;
00188                   if(Debug()) if(mPrint) { 
00189                     cout << eemcDetname[detnum] 
00190                        << " cluster = " << j  
00191                        << ", m = " <<module <<" " << clModule
00192                        << ", e = " << eta << " " << clEta
00193                        << ", s = " << sub << " " << clSub << endl; 
00194                   }
00195                   Float_t hitEnergy=hits[k]->dE();
00196                   trackEnergyFraction+=hitEnergy/energy;
00197                   Float_t clEnergy=clusters[j]->energy();
00198                   clusterEnergyFraction+=hitEnergy/clEnergy;
00199                 }
00200               }
00201             }
00202             if (assoc) {
00203               StEEmcClusterAssociation *c=new StEEmcClusterAssociation(tracks[i],clusters[j],trackEnergyFraction,clusterEnergyFraction);
00204               mTrackCluster[detnum]->insert(multiEEmcTrackClusterValue(tracks[i],c));
00205 
00206               mClusterTrack[detnum]->insert(multiEEmcClusterTrackValue(clusters[j],c));
00207             }             
00208           }
00209         }
00210       }
00211     }
00212   }
00213   // Finishing doing cluster association
00214   
00215   // Starting doing point association 
00216   if(mPrint) cout <<"Doing point association\n";
00217   StSPtrVecEmcPoint& endcapPoints=emcCollection->endcapPoints();
00218   cout <<"points.size() = " << endcapPoints.size() << endl;
00219   for (UInt_t j=0; j<endcapPoints.size(); j++) {
00220       int nCl[4];      
00221       for(Int_t detnum=0;detnum<NDETECTORS;detnum++){
00222         StDetectorId detId=static_cast<StDetectorId>(detnum+kEndcapEmcTowerId);
00223         StPtrVecEmcCluster& cluster=endcapPoints[j]->cluster(detId);
00224         nCl[detnum] = cluster.size();
00225       }
00226       if(nCl[0]!=0 && nCl[1]!=0 && nCl[2]!=1 && nCl[3]!=1 &&
00227          endcapPoints[j]->nTracks()!=0) endcapPoints[j]->print();
00228   }
00229   if(endcapPoints.size()!=0)
00230   {
00231     // Allocate point matching maps
00232     if(!mTrackPoint) mTrackPoint = new multiEEmcTrackPoint;
00233     if(!mPointTrack) mPointTrack = new multiEEmcPointTrack;
00234     
00235     for (UInt_t i=0; i<tracks.size(); i++)
00236     {
00237       StPtrVecMcCalorimeterHit hits;
00238 
00239       for (UInt_t j=0; j<endcapPoints.size(); j++)
00240       {
00241         Int_t assoc=0;
00242         for(Int_t detnum=0;detnum<NDETECTORS;detnum++)
00243         {
00244           switch (detnum)
00245           {
00246             case 0: { hits=tracks[i]->eemcHits(); break; }
00247             case 1: { hits=tracks[i]->eprsHits(); break; }
00248             case 2: { hits=tracks[i]->esmduHits(); break; }
00249             case 3: { hits=tracks[i]->esmdvHits(); break; }
00250             default: { gMessMgr->Warning()<<"Detector does not exist"<<endm; break; }
00251           }
00252           StDetectorId detId=static_cast<StDetectorId>(detnum+kEndcapEmcTowerId);
00253           StPtrVecEmcCluster& cluster=endcapPoints[j]->cluster(detId);
00254           for (UInt_t k=0; k<hits.size(); k++) if(hits[k])
00255           {
00256             UInt_t module=hits[k]->module();
00257             UInt_t eta=hits[k]->eta();
00258             Int_t sub=hits[k]->sub();
00259 
00260             for (UInt_t j2=0; j2<cluster.size(); j2++) if(cluster[j2])
00261             {
00262               StPtrVecEmcRawHit& clHit=cluster[j2]->hit();
00263               for (UInt_t l=0; l<clHit.size(); l++) if(clHit[l])
00264               {
00265                 UInt_t clModule=clHit[l]->module();
00266                 UInt_t clEta=clHit[l]->eta();
00267                 Int_t clSub=abs(clHit[l]->sub());
00268 
00269                 // compare between hit track and hit cluster from a point
00270                 if (module==clModule && eta==clEta 
00271                       && ((detnum > 1)? 1 : sub==clSub)) 
00272                 { 
00273                     assoc+=(1<<(detnum)); 
00274                     if(Debug()) if(mPrint) {
00275                       cout <<"assoc = " << assoc << " " 
00276                        << eemcDetname[detnum] 
00277                        << ", m = " <<module <<" " << clModule
00278                        << ", e = " << eta << " " << clEta
00279                        << ", s = " << sub << " " << clSub << endl; 
00280                     }
00281                     goto nextDetector;
00282                 }
00283               }
00284             }
00285           }
00286           nextDetector: continue;
00287         } 
00288         if (assoc) {
00289           StEEmcPointAssociation *p= new StEEmcPointAssociation(tracks[i],endcapPoints[j],assoc);
00290           mTrackPoint->insert(multiEEmcTrackPointValue(tracks[i],p));
00291           mPointTrack->insert(multiEEmcPointTrackValue(endcapPoints[j],p));
00292         }
00293       }
00294     }   
00295   }  
00296   if(Debug()) if(mPrint) printMaps();
00297   // Finishing doing point association 
00298   return kStOK;
00299 }
00300 //------------------------------------------------------------------------------
00301 Int_t StEEmcAssociationMaker::Finish()
00302 {
00303   // Finishing the maker
00304   gMessMgr->Info() << "StEEmcAssociationMaker::Finish()" << endm;
00305   return kStOK;
00306 }
00307 //---------------------------------------------------------------------
00308 multiEEmcTrackCluster* StEEmcAssociationMaker::getTrackClusterMap(const char* detname)
00309 {
00310   Int_t det = getDetNum(detname);
00311   if(det>0) return mTrackCluster[det-1];
00312   return NULL;
00313 }     
00314 //---------------------------------------------------------------------
00315 multiEEmcClusterTrack* StEEmcAssociationMaker::getClusterTrackMap(const char* detname)
00316 {
00317   Int_t det = getDetNum(detname);
00318   if(det>0) return mClusterTrack[det-1];
00319   return NULL;
00320 }     
00321 //---------------------------------------------------------------------
00322 Int_t StEEmcAssociationMaker::getDetNum(const char* detname)
00323 {
00324   Int_t detnum = 0;
00325   for (Int_t i=0; i<NDETECTORS; i++)  
00326     if (!strcmp(detname,eemcDetname[i])) detnum = i+1;
00327   return detnum;
00328 }     
00329 //---------------------------------------------------------------------
00330 void StEEmcAssociationMaker::printMaps()
00331 {
00332   for(int i=0;i<NDETECTORS;i++)
00333   {
00334     cout <<"-----------------------------------------------------------\n";
00335     cout <<"ASSOCIATION FOR DETECTOR "<<eemcDetname[i].Data()<<endl;
00336     multiEEmcTrackCluster* map = getTrackClusterMap(i+1);
00337     if(map) 
00338     { 
00339       cout <<"Track->Cluster Association Map\n";
00340       for(multiEEmcTrackClusterIter j=map->begin(); j!=map->end(); j++) 
00341       {
00342         StMcTrack* track = (StMcTrack*)(*j).first;
00343         StEEmcClusterAssociation* value = (StEEmcClusterAssociation*)(*j).second;
00344         if(track && value)
00345         {
00346           StEmcCluster *c = (StEmcCluster*) value->getCluster();
00347           
00348           if(c) 
00349           {
00350             cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00351                  <<" Cl = "<<c<<" E = "<<c->energy()<<" eta = "<<c->eta()<<" phi = "<<c->phi()
00352                  <<" FrTr = "<<value->getFractionTrack()<<" FrCl = "<<value->getFractionCluster()<<endl;
00353           }
00354         }
00355       }      
00356     }
00357     
00358     multiEEmcClusterTrack* map1 = getClusterTrackMap(i+1);
00359     if(map1) 
00360     { 
00361       cout <<"Cluster->Track Association Map\n";
00362       for(multiEEmcClusterTrackIter j=map1->begin(); j!=map1->end(); j++) 
00363       {
00364         StEmcCluster *c = (StEmcCluster*)(*j).first;
00365         StEEmcClusterAssociation* value = (StEEmcClusterAssociation*)(*j).second;
00366         if(c && value)
00367         {               
00368           StMcTrack* track = (StMcTrack*)value->getTrack();
00369           if(track) 
00370           {
00371             cout <<" Cl = "<<c<<" E = "<<c->energy()<<" eta = "<<c->eta()<<" phi = "<<c->phi()
00372                  <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00373                  <<" FrTr = "<<value->getFractionTrack()<<" FrCl = "<<value->getFractionCluster()<<endl;
00374           }
00375         }
00376       }      
00377     }
00378   }
00379   multiEEmcTrackPoint *mapPoint = getTrackPointMap();
00380   if(mapPoint)
00381   {
00382     cout <<"Track->Point Association Map\n";
00383     for(multiEEmcTrackPointIter j = mapPoint->begin(); j!=mapPoint->end(); j++)
00384     {
00385       StMcTrack *track = (StMcTrack*)(*j).first;
00386       StEEmcPointAssociation *value = (StEEmcPointAssociation*)(*j).second;
00387       if(value)
00388       {
00389         StEmcPoint *point = (StEmcPoint*)value->getPoint(); 
00390         if(track && point)
00391         {
00392           cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00393                <<" Point = "<<point<<" E = "<<point->energy()
00394                <<" Assoc. = "<<value->getAssociation();
00395           for(int i=1;i<=4;i++) cout <<" det "<<i<<" A = "<< value->getAssociation(i);
00396           cout <<endl;
00397         }
00398       }
00399     }
00400   }
00401   multiEEmcPointTrack *mapPoint1 = getPointTrackMap();
00402   if(mapPoint1)
00403   {
00404     cout <<"Point->Track Association Map\n";
00405     for(multiEEmcPointTrackIter j = mapPoint1->begin(); j!=mapPoint1->end(); j++)
00406     {
00407       StEmcPoint *point = (StEmcPoint*)(*j).first;
00408       StEEmcPointAssociation *value = (StEEmcPointAssociation*)(*j).second;
00409       if(value)
00410       {
00411         StMcTrack *track = (StMcTrack*)value->getTrack();
00412         if(track && point)
00413         {
00414           cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00415                <<" Point = "<<point<<" E = "<<point->energy()
00416                <<" Assoc. = "<<value->getAssociation();
00417           for(int i=1;i<=4;i++) cout <<" det "<<i<<" A = "<< value->getAssociation(i);
00418           cout <<endl;
00419         }
00420       }
00421     }
00422   }
00423 }
00424 //-------------------------------------------------------------------
00425 void StEEmcAssociationMaker::printHits(StEvent *event)
00426 {
00427   StEmcCollection* emccol=(StEmcCollection*)event->emcCollection();
00428 
00429     gMessMgr->Info()<<"****** EEMC Association printHits() *******" <<endm;
00430   for(Int_t i=0; i<NDETECTORS; i++)
00431   {
00432     StDetectorId id = static_cast<StDetectorId>(i+kEndcapEmcTowerId);
00433     StEmcDetector* detector=emccol->detector(id);
00434     gMessMgr->Info()<<"****************** hits in detector "
00435                                  << eemcDetname[i].Data()<<endm;
00436     if(detector) for(int j=1;j<=kEEmcNumSectors;j++)
00437     {
00438       StEmcModule* sector = detector->module(j);
00439       StSPtrVecEmcRawHit& rawHit=sector->hits();
00440       if(rawHit.size()>0)
00441           gMessMgr->Info()<<"Number of hits for module "
00442                                       <<j<<" = "<<rawHit.size()<<endm;
00443       for(UInt_t k=0;k<rawHit.size();k++)
00444         if(i < 2) // for Tower and Prs
00445           gMessMgr->Info()<<"Hit number = "<<k
00446                      <<"  sector = " << rawHit[k]->module()
00447                      <<"  eta = "<<rawHit[k]->eta()
00448                      <<"  sub = "<< rawHit[k]->sub()
00449                      <<"  adc = "<< rawHit[k]->adc()
00450                      <<"  energy = "<<rawHit[k]->energy()<<endm;
00451         else
00452           gMessMgr->Info()<<"Hit number = "<<k
00453                      <<"  sector = " << rawHit[k]->module()
00454                      <<"  strip = "<<rawHit[k]->eta()
00455                      <<"  adc = "<< rawHit[k]->adc()
00456                      <<"  energy = "<<rawHit[k]->energy()<<endm;
00457     }
00458   }
00459 }
00460 
00462 //
00463 // $Id: StEEmcAssociationMaker.cxx,v 1.2 2005/09/29 14:58:10 fisyak Exp $
00464 // $Log: StEEmcAssociationMaker.cxx,v $
00465 // Revision 1.2  2005/09/29 14:58:10  fisyak
00466 // Persistent StMcEvent
00467 //
00468 // Revision 1.1.1.1  2005/05/31 18:54:13  wzhang
00469 // First version
00470 //
00471 //