StRoot: StAssociationMaker/EMC/StEmcAssociationMaker.cxx Source File

00001 /********************************************************************************
00002  *
00003  * Author: Marcia Maria de Moura
00004  *
00005  * See README file for description
00006  *
00007  ********************************************************************************/
00008 
00009 #include "StEmcAssociationMaker.h"
00010 
00011 #include <Stiostream.h>
00012 #include <math.h>
00013 
00014 #include "StEmcUtil/geometry/StEmcGeom.h"
00015 #include "StMcEvent/StMcCalorimeterHit.hh"
00016 #include "StMcEventTypes.hh"
00017 #include "StMcEvent.hh"
00018 #include "StEventTypes.h"
00019 #include <StMessMgr.h>
00020 
00021 ClassImp(StEmcAssociation)
00022 ClassImp(StEmcClusterAssociation)
00023 ClassImp(StEmcPointAssociation)
00024 ClassImp(StEmcAssociationMaker)
00025 //------------------------------------------------------------------------------
00026 StEmcAssociation::StEmcAssociation(StMcTrack *t) 
00027 { mTrack = t; }
00028 StEmcAssociation::~StEmcAssociation() 
00029 { }
00030 //------------------------------------------------------------------------------
00031 StEmcClusterAssociation::StEmcClusterAssociation(StMcTrack* t, StEmcCluster* c, float ft,float fe):StEmcAssociation(t) 
00032 { mCluster=c;  mFTrack = ft; mFEmc = fe;}
00033 StEmcClusterAssociation::~StEmcClusterAssociation()
00034 { }
00035 //------------------------------------------------------------------------------
00036 StEmcPointAssociation::StEmcPointAssociation(StMcTrack* t, StEmcPoint* p, int at):StEmcAssociation(t) 
00037 { mPoint=p; mAssocType = at; }
00038 StEmcPointAssociation::~StEmcPointAssociation() 
00039 { }
00040 //------------------------------------------------------------------------------
00041 StEmcAssociationMaker::StEmcAssociationMaker(const char* name):StMaker(name)
00042 {
00043   for(Int_t i=0;i<NDETECTORS;i++)
00044   {
00045     mTrackCluster[i] = NULL;
00046     mClusterTrack[i] = NULL;
00047   }
00048   mTrackPoint = NULL;
00049   mPointTrack = NULL;
00050   
00051   mPrint = kTRUE;
00052 }
00053 //------------------------------------------------------------------------------
00054 StEmcAssociationMaker::~StEmcAssociationMaker()
00055 {
00056 }
00057 void StEmcAssociationMaker::Clear(const char* a)
00058 {
00059   if(mPrint) cout <<"Cleaning old stuff from EMC association\n";
00060   for(Int_t i=0;i<NDETECTORS;i++)
00061   {
00062     if(mPrint) cout <<"Cleaning for detector "<<i<<endl;
00063     if(mTrackCluster[i]) 
00064     { 
00065       for(multiEmcTrackClusterIter j=mTrackCluster[i]->begin(); j!=mTrackCluster[i]->end(); j++) 
00066       {
00067         SafeDelete((*j).second);
00068       }
00069       mTrackCluster[i]->clear(); 
00070       SafeDelete(mTrackCluster[i]); 
00071       mTrackCluster[i]=NULL; 
00072     }
00073     if(mClusterTrack[i]) 
00074     { 
00075       mClusterTrack[i]->clear(); 
00076       SafeDelete(mClusterTrack[i]); 
00077       mClusterTrack[i]=NULL;
00078     }
00079   }
00080   if (mPrint) cout <<"Cleaning points Association\n";
00081   if(mTrackPoint) 
00082   { 
00083     for(multiEmcTrackPointIter j=mTrackPoint->begin(); j!=mTrackPoint->end(); j++) 
00084     {
00085       SafeDelete((*j).second);
00086     }
00087     mTrackPoint->clear(); 
00088     SafeDelete(mTrackPoint); 
00089     mTrackPoint=NULL;
00090   }
00091   if(mPointTrack) 
00092   { 
00093     mPointTrack->clear(); 
00094     SafeDelete(mPointTrack); 
00095     mPointTrack=NULL;
00096   }
00097   return;
00098 }
00099 //------------------------------------------------------------------------------
00100 Int_t StEmcAssociationMaker::Init()
00101 {
00102   return StMaker::Init();
00103 }
00104 //------------------------------------------------------------------------------
00105 Int_t StEmcAssociationMaker::Make()
00106 {
00107   if(mPrint) gMessMgr->Info() << "StEmcAssociationMaker::Make()" << endm;
00108   
00109   // Getting McEvent object
00110   StMcEvent* mcEvent= (StMcEvent*) GetDataSet("StMcEvent");
00111   if (!mcEvent) return kStWarn;
00112   StSPtrVecMcTrack& tracks=mcEvent->tracks();
00113   if(tracks.size()==0) return kStWarn;
00114 
00115   // Getting Event object
00116   StEvent* event=(StEvent*) GetDataSet("StEvent");
00117   if (!event) return kStWarn;
00118   
00119   StEmcCollection* emcCollection=event->emcCollection();
00120   if(!emcCollection) return kStWarn;
00121 
00122   // Starting doing cluster association for each detector
00123   for (Int_t detnum=0; detnum<NDETECTORS; detnum++) // For detnum<4, only barrel EMC
00124   {
00125     if(mPrint) cout <<"Doing association for detector "<<detnum<<endl;
00126     StDetectorId detId=static_cast<StDetectorId>(detnum+kBarrelEmcTowerId);
00127     StEmcDetector* detector=emcCollection->detector(detId);
00128     StEmcClusterCollection* clusterColl=NULL;
00129     if (detector) clusterColl=detector->cluster();
00130     if (clusterColl)
00131     { 
00132       StSPtrVecEmcCluster& clusters=clusterColl->clusters();
00133 
00134       // Allocate cluster matching maps
00135       if(mPrint) cout<<"Track size = "<<tracks.size()<<"  cluster size = "<<clusters.size()<<endl;
00136       if(!mTrackCluster[detnum]) mTrackCluster[detnum]=new multiEmcTrackCluster;
00137       if(!mClusterTrack[detnum]) mClusterTrack[detnum]=new multiEmcClusterTrack;
00138       
00139       for (UInt_t i=0; i<tracks.size(); i++)
00140       {
00141         StPtrVecMcCalorimeterHit hits;
00142         switch (detnum)
00143         {
00144           case 0: { hits=tracks[i]->bemcHits(); break; }
00145           case 1: { hits=tracks[i]->bprsHits(); break; }
00146           case 2: { hits=tracks[i]->bsmdeHits(); break; }
00147           case 3: { hits=tracks[i]->bsmdpHits(); break; }
00148           default: { gMessMgr->Warning()<<"Detector does not exist"<<endm; break; }
00149         }
00150 
00151         if (hits.size()!=0);
00152         {
00153           // Calculating track energy by summing over hits energies
00154           Float_t energy=0;
00155           for (UInt_t i2=0; i2<hits.size(); i2++)
00156           {
00157             Float_t hitEnergy=dEToEnergy(hits[i2],detnum);
00158             energy+=hitEnergy;
00159           } 
00160           
00161           for (UInt_t j=0; j<clusters.size(); j++)
00162           {
00163             StPtrVecEmcRawHit& clHits=clusters[j]->hit();
00164             Float_t trackEnergyFraction=0;
00165             Float_t clusterEnergyFraction=0;
00166             Int_t assoc=0;
00167             for (UInt_t k=0; k<hits.size(); k++)
00168             {
00169               UInt_t module=hits[k]->module();
00170               UInt_t eta=hits[k]->eta();
00171               Int_t sub=hits[k]->sub();
00172               
00173               for (UInt_t l=0; l<clHits.size(); l++)
00174               {
00175                 UInt_t clModule=clHits[l]->module();
00176                 UInt_t clEta=clHits[l]->eta();
00177                 Int_t clSub=abs(clHits[l]->sub());
00178                 //Doing comparision between hit track and hit cluster
00179                 if (module==clModule && eta==clEta && sub==clSub)
00180                 {
00181                   assoc=1;
00182                   Float_t hitEnergy=dEToEnergy(hits[k],detnum);
00183                   trackEnergyFraction+=hitEnergy/energy;
00184                   Float_t clEnergy=clusters[j]->energy();
00185                   clusterEnergyFraction+=hitEnergy/clEnergy;
00186                 }
00187               }
00188             }
00189             if (assoc) {
00190               StEmcClusterAssociation *c=new StEmcClusterAssociation(tracks[i],clusters[j],trackEnergyFraction,clusterEnergyFraction);
00191               mTrackCluster[detnum]->insert(multiEmcTrackClusterValue(tracks[i],c));
00192 
00193               mClusterTrack[detnum]->insert(multiEmcClusterTrackValue(clusters[j],c));
00194             }             
00195           }
00196         }
00197       }
00198     }
00199   }
00200   // Finishing doing cluster association
00201   
00202   // Starting doing point association 
00203   if(mPrint) cout <<"Doing point association\n";
00204   StSPtrVecEmcPoint& barrelPoints=emcCollection->barrelPoints();
00205   if(barrelPoints.size()!=0)
00206   {
00207     // Allocate point matching maps
00208     if(!mTrackPoint) mTrackPoint = new multiEmcTrackPoint;
00209     if(!mPointTrack) mPointTrack = new multiEmcPointTrack;
00210     
00211     for (UInt_t i=0; i<tracks.size(); i++)
00212     {
00213       StPtrVecMcCalorimeterHit hits;
00214 
00215       for (UInt_t j=0; j<barrelPoints.size(); j++)
00216       {
00217         Int_t assoc=0;
00218         for(Int_t detnum=0;detnum<NDETECTORS;detnum++)
00219         {
00220           switch (detnum)
00221           {
00222             case 0: { hits=tracks[i]->bemcHits(); break; }
00223             case 1: { hits=tracks[i]->bprsHits(); break; }
00224             case 2: { hits=tracks[i]->bsmdeHits(); break; }
00225             case 3: { hits=tracks[i]->bsmdpHits(); break; }
00226             default: { gMessMgr->Warning()<<"Detector does not exist"<<endm; break; }
00227           }
00228           StDetectorId detId=static_cast<StDetectorId>(detnum+kBarrelEmcTowerId);
00229           StPtrVecEmcCluster& cluster=barrelPoints[j]->cluster(detId);
00230           for (UInt_t k=0; k<hits.size(); k++) if(hits[k])
00231           {
00232             UInt_t module=hits[k]->module();
00233             UInt_t eta=hits[k]->eta();
00234             Int_t sub=hits[k]->sub();
00235 
00236             for (UInt_t j2=0; j2<cluster.size(); j2++) if(cluster[j2])
00237             {
00238               StPtrVecEmcRawHit& clHit=cluster[j2]->hit();
00239               for (UInt_t l=0; l<clHit.size(); l++) if(clHit[l])
00240               {
00241                 UInt_t clModule=clHit[l]->module();
00242                 UInt_t clEta=clHit[l]->eta();
00243                 Int_t clSub=abs(clHit[l]->sub());
00244 
00245                 //Doing comparision between hit track and hit cluster from a point
00246                 if (module==clModule && eta==clEta && sub==clSub) { assoc+=(1<<(detnum)); goto nextDetector;}
00247               }
00248             }
00249           }
00250           nextDetector: continue;
00251         } 
00252         if (assoc) {
00253           StEmcPointAssociation *p= new StEmcPointAssociation(tracks[i],barrelPoints[j],assoc);
00254           mTrackPoint->insert(multiEmcTrackPointValue(tracks[i],p));
00255           mPointTrack->insert(multiEmcPointTrackValue(barrelPoints[j],p));
00256         }
00257       }
00258     }   
00259   }  
00260   // Finishing doing point association 
00261   return kStOK;
00262 }
00263 //------------------------------------------------------------------------------
00264 Int_t StEmcAssociationMaker::Finish()
00265 {
00266   // Finishing the maker
00267   gMessMgr->Info() << "StEmcAssociationMaker::Finish()" << endm;
00268   return kStOK;
00269 }
00270 //------------------------------------------------------------------------------
00271 Float_t StEmcAssociationMaker::dEToEnergy(StMcCalorimeterHit* hit, Int_t detnum)
00272 {
00273   // Get hit deposited energy and convert to energy for each EMC detector. 
00274   
00275   Float_t P0[]={14.69,559.7,0.1185e6,0.1260e6};
00276   Float_t P1[]={-0.1022,-109.9,-0.3292e5,-0.1395e5};
00277   Float_t P2[]={0.7484,-97.81,0.3113e5,0.1971e5};
00278 
00279   UInt_t m=hit->module();
00280   UInt_t e=hit->eta();
00281   Float_t dE=hit->dE();
00282   Float_t hitEnergy = 0;
00283 
00284   StEmcGeom* emcGeom=StEmcGeom::instance(detnum+1);
00285   Float_t Eta;
00286   emcGeom->getEta(m,e,Eta);
00287   Float_t x=fabs(Eta);
00288   Float_t sf=P0[detnum]+P1[detnum]*x+P2[detnum]*x*x;
00289   hitEnergy=dE*sf;
00290   if (hitEnergy<=0) hitEnergy=0;
00291   return hitEnergy;
00292 }
00293 //---------------------------------------------------------------------
00294 multiEmcTrackCluster* StEmcAssociationMaker::getTrackClusterMap(const char* detname)
00295 {
00296   Int_t det = getDetNum(detname);
00297   if(det>0) return mTrackCluster[det-1];
00298   return NULL;
00299 }     
00300 //---------------------------------------------------------------------
00301 multiEmcClusterTrack* StEmcAssociationMaker::getClusterTrackMap(const char* detname)
00302 {
00303   Int_t det = getDetNum(detname);
00304   if(det>0) return mClusterTrack[det-1];
00305   return NULL;
00306 }     
00307 //---------------------------------------------------------------------
00308 Int_t StEmcAssociationMaker::getDetNum(const char* detname)
00309 {
00310   const TString det[] = {"bemc","bprs","bsmde","bsmdp"};
00311   Int_t detnum = 0;
00312   for (Int_t i=0; i<NDETECTORS; i++)  
00313     if (!strcmp(detname,det[i])) detnum = i+1;
00314   return detnum;
00315 }     
00316 //---------------------------------------------------------------------
00317 void StEmcAssociationMaker::printMaps()
00318 {
00319   const TString det[] = {"bemc","bprs","bsmde","bsmdp"};
00320   for(int i=0;i<NDETECTORS;i++)
00321   {
00322     cout <<"-----------------------------------------------------------\n";
00323     cout <<"ASSOCIATION FOR DETECTOR "<<det[i].Data()<<endl;
00324     multiEmcTrackCluster* map = getTrackClusterMap(i+1);
00325     if(map) 
00326     { 
00327       cout <<"Track->Cluster Association Map\n";
00328       for(multiEmcTrackClusterIter j=map->begin(); j!=map->end(); j++) 
00329       {
00330         StMcTrack* track = (StMcTrack*)(*j).first;
00331         StEmcClusterAssociation* value = (StEmcClusterAssociation*)(*j).second;
00332         if(track && value)
00333         {
00334           StEmcCluster *c = (StEmcCluster*) value->getCluster();
00335           
00336           if(c) 
00337           {
00338             cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00339                  <<" Cl = "<<c<<" E = "<<c->energy()<<" eta = "<<c->eta()<<" phi = "<<c->phi()
00340                  <<" FrTr = "<<value->getFractionTrack()<<" FrCl = "<<value->getFractionCluster()<<endl;
00341           }
00342         }
00343       }      
00344     }
00345     cout <<endl;
00346     multiEmcClusterTrack* map1 = getClusterTrackMap(i+1);
00347     if(map1) 
00348     { 
00349       cout <<"Cluster->Track Association Map\n";
00350       for(multiEmcClusterTrackIter j=map1->begin(); j!=map1->end(); j++) 
00351       {
00352         StEmcCluster *c = (StEmcCluster*)(*j).first;
00353         StEmcClusterAssociation* value = (StEmcClusterAssociation*)(*j).second;
00354         if(c && value)
00355         {               
00356           StMcTrack* track = (StMcTrack*)value->getTrack();
00357           if(track) 
00358           {
00359             cout <<" Cl = "<<c<<" E = "<<c->energy()<<" eta = "<<c->eta()<<" phi = "<<c->phi()
00360                  <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00361                  <<" FrTr = "<<value->getFractionTrack()<<" FrCl = "<<value->getFractionCluster()<<endl;
00362           }
00363         }
00364       }      
00365     }
00366   }
00367   multiEmcTrackPoint *mapPoint = getTrackPointMap();
00368   if(mapPoint)
00369   {
00370     cout <<"Track->Point Association Map\n";
00371     for(multiEmcTrackPointIter j = mapPoint->begin(); j!=mapPoint->end(); j++)
00372     {
00373       StMcTrack *track = (StMcTrack*)(*j).first;
00374       StEmcPointAssociation *value = (StEmcPointAssociation*)(*j).second;
00375       if(value)
00376       {
00377         StEmcPoint *point = (StEmcPoint*)value->getPoint(); 
00378         if(track && point)
00379         {
00380           cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00381                <<" Point = "<<point<<" E = "<<point->energy()
00382                <<" Assoc. = "<<value->getAssociation();
00383           for(int i=1;i<=4;i++) cout <<" det "<<i<<" A = "<< value->getAssociation(i);
00384           cout <<endl;
00385         }
00386       }
00387     }
00388   }
00389   multiEmcPointTrack *mapPoint1 = getPointTrackMap();
00390   if(mapPoint1)
00391   {
00392     cout <<"Point->Track Association Map\n";
00393     for(multiEmcPointTrackIter j = mapPoint1->begin(); j!=mapPoint1->end(); j++)
00394     {
00395       StEmcPoint *point = (StEmcPoint*)(*j).first;
00396       StEmcPointAssociation *value = (StEmcPointAssociation*)(*j).second;
00397       if(value)
00398       {
00399         StMcTrack *track = (StMcTrack*)value->getTrack();
00400         if(track && point)
00401         {
00402           cout <<" McTrack = "<<track<<" GeantId = "<<track->geantId()<<" pt = "<<track->pt()<<" TReta = "<<track->pseudoRapidity()
00403                <<" Point = "<<point<<" E = "<<point->energy()
00404                <<" Assoc. = "<<value->getAssociation();
00405           for(int i=1;i<=4;i++) cout <<" det "<<i<<" A = "<< value->getAssociation(i);
00406           cout <<endl;
00407         }
00408       }
00409     }
00410   }
00411 }
00412 void StEmcAssociationMaker::printTracks()
00413 {
00414   StMcEvent* mcEvent=(StMcEvent*) GetDataSet("StMcEvent");
00415   if (!mcEvent) return;
00416   StSPtrVecMcTrack& tracks=mcEvent->tracks();
00417   if(tracks.size()==0) return;
00418 
00419   StMcVertex *v = mcEvent->primaryVertex();
00420   if(v)
00421   {
00422     cout<<"Primary Vertex (x,y,z) = "<<v->position().x()<<"  "<<v->position().y()<<"  "<<v->position().z()<<"  "<<endl;
00423   }
00424   for(UInt_t i = 0; i<tracks.size();i++)
00425   {
00426     cout <<"Track " <<tracks[i]<<" Geant Id = "<<tracks[i]->geantId()<<"  pT = "<< tracks[i]->pt()
00427          <<"  eta = "<<tracks[i]->pseudoRapidity()
00428          <<"  phi = "<<tracks[i]->momentum().phi()<<endl;
00429     cout <<"     Parent = "<<tracks[i]->parent()<<endl;
00430     v = tracks[i]->startVertex();
00431     if(v) cout <<"     Start vertex (x,y,z) = "<<v->position().x()<<"  "<<v->position().y()<<"  "<<v->position().z()<<"  "<<endl;  
00432     v = tracks[i]->stopVertex();
00433     if(v) cout <<"     Stop vertex (x,y,z)  = "<<v->position().x()<<"  "<<v->position().y()<<"  "<<v->position().z()<<"  "<<endl;  
00434   }
00435 }