StRoot: StEmcMixerMaker/StEmcMixerMaker.cxx Source File

00001 #include "StEmcMixerMaker.h"
00002 #include <Stiostream.h>
00003 #include <math.h>
00004 #include "SystemOfUnits.h"
00005 #include <stdlib.h>
00006 #include <string.h>
00007 #include "TMath.h"
00008 #include "TFile.h"
00009 #include "Stypes.h"
00010 #include "math_constants.h"
00011 #include "StEmcSimulatorMaker/StEmcSimulatorMaker.h"
00012 #include "StEventTypes.h"
00013 #include "StEvent.h"
00014 #include "StEmcUtil/geometry/StEmcGeom.h"
00015 #include "tables/St_emcStatus_Table.h"
00016 #include "tables/St_smdStatus_Table.h"
00017                                                   
00018 ClassImp(StEmcMixerMaker)
00019 
00020 //_____________________________________________________________________________
00021 StEmcMixerMaker::StEmcMixerMaker(const char *name):StMaker(name)
00022 {
00023   mAddHits = kTRUE;
00024   mClear = kTRUE;
00025   mUseDB = kTRUE;
00026   mEmbedAll = kTRUE;
00027   mFakeTrackEmbed = kFALSE;
00028   for(Int_t i=0;i<NDETECTORS;i++) mGeom[i]=StEmcGeom::instance(i+1);
00029 }
00030 //_____________________________________________________________________________
00031 StEmcMixerMaker::~StEmcMixerMaker()
00032 {
00033 }
00034 //_____________________________________________________________________________
00035 Int_t StEmcMixerMaker::Init()
00036 {
00037   m_hit_1 = new TH1F("old_hit","old_hit",100,0.,30.);
00038   m_hit_2 = new TH1F("new_hit","new_hit",100,0.,30.);
00039   m_edep_1 = new TH1F("EDEP1","edep1",100,0.,30.);
00040   m_edep_2 = new TH1F("EDEP2","edep2",100,0.,30.);
00041   return StMaker::Init();
00042 }  
00043 //_____________________________________________________________________________
00044 Int_t StEmcMixerMaker::Make()
00045 {
00046   if(!getEvents()) return kStWarn;
00047   
00048   if(mUseDB) getDB(); // get status tables from DB. If not, set all status to 1 (valid channel)
00049   else 
00050     for(Int_t i=0;i<NDETECTORS;i++) 
00051       for(Int_t j=0;j<(EMCCHANNELS+(MAXCHANNELS-EMCCHANNELS)*(i>1));j++) 
00052         mStatus[i][j] = 1;
00053   
00054   clearPoints(); // clear EMC points
00055   clearClusters();  // clear EMC clusters
00056 
00057   if(mAddHits){ if(addHits()!=kStOk) { LOG_WARN <<" error in addhits***"<<endm; return kStWarn; } }
00058   if(mFakeTrackEmbed) addTracks(); 
00059   return kStOK;
00060 }
00061 
00062 //_____________________________________________________________________________
00063 Int_t StEmcMixerMaker::Finish() 
00064 {
00065   return StMaker::Finish();
00066 }
00067 //--------------------------------------------------------------------
00072 Int_t StEmcMixerMaker::addHits()
00073 {   
00074   StEmcCollection* emccol1=(StEmcCollection*)mEvent1->emcCollection(); 
00075   StEmcCollection* emccol2=(StEmcCollection*)mEvent2->emcCollection();
00076   Int_t stat_h2_all[MAXCHANNELS];
00077   Int_t stat_h2[MAXCHANNELS];
00078 
00079 // Loop over all four sub detectors
00080   Float_t old_edep_tot=0.0;
00081   Float_t new_edep_tot=0.0;
00082   for(Int_t i=0; i<NDETECTORS; i++) 
00083   {
00084     StDetectorId id = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00085     StEmcDetector* detector1=emccol1->detector(id);
00086     StEmcDetector* detector2=emccol2->detector(id);
00087     if(!detector1){ LOG_WARN <<"detector1 not loaded"<<endm; }
00088         if(!detector2){ LOG_WARN <<"detector2 not loaded"<<endm; }
00089    
00090     Float_t edep1_tot=0;
00091     Float_t edep2_tot=0;
00092  
00093     if(detector1 && detector2) for(UInt_t j=1;j<=NMODULES;j++) // loop over all modules
00094     {
00095       StEmcModule* module1 = detector1->module(j);
00096       StEmcModule* module2 = detector2->module(j);
00097       StSPtrVecEmcRawHit& rawHit1=module1->hits();
00098       StSPtrVecEmcRawHit& rawHit2=module2->hits();
00099 
00100       for(UInt_t k1=0;k1<rawHit1.size();k1++) {edep1_tot+=rawHit1[k1]->energy();}
00101       for(UInt_t k1=0;k1<rawHit2.size();k1++) {edep2_tot+=rawHit2[k1]->energy();}
00102 
00103       for(UInt_t is=0;is<rawHit2.size();is++) {stat_h2_all[is]=0;stat_h2[is]=0;}
00104  
00105       for(UInt_t k1=0;k1<rawHit1.size();k1++)  
00106       {
00107         Int_t m1, e1, s1; 
00108         m1=(Int_t)rawHit1[k1]->module();
00109         e1=(Int_t)rawHit1[k1]->eta();
00110         s1=abs(rawHit1[k1]->sub());                                                   
00111         for(UInt_t is=0;is<rawHit2.size();is++) {stat_h2[is]=0;}
00112 
00113         Float_t edep_add=0.0;
00114         UInt_t adc_add=0; 
00115         for(UInt_t k2=0;k2<rawHit2.size();k2++) if(checkHit(i,rawHit2[k2]))
00116         {
00117           if(stat_h2[k2]==1)continue;
00118           Int_t m2, e2, s2;
00119           m2=(Int_t)rawHit2[k2]->module();
00120           e2=(Int_t)rawHit2[k2]->eta();
00121           s2=abs(rawHit2[k2]->sub());          
00122           if(m1==m2 && e1==e2 && s1==s2)
00123           {
00124             stat_h2[k2]=1;
00125             stat_h2_all[k2]=1;
00126             edep_add=rawHit2[k2]->energy();
00127             adc_add=rawHit2[k2]->adc();
00128           }
00129           if(stat_h2[k2]==1)break;
00130         }
00131         Float_t new_edep=rawHit1[k1]->energy()+edep_add;
00132         UInt_t new_adc=rawHit1[k1]->adc()+adc_add;
00133         Float_t oldE=rawHit1[k1]->energy();
00134         old_edep_tot+=rawHit1[k1]->energy();
00135         rawHit1[k1]->setAdc(new_adc);
00136         rawHit1[k1]->setEnergy(new_edep);
00137         new_edep_tot+=new_edep;
00138         UInt_t calib = rawHit1[k1]->calibrationType();
00139         while(calib>127) calib-=128;
00140         if(edep_add>0) {
00141           LOG_DEBUG <<"EMBEDDED HIT -> det = "<<i+1<<"  m = "<<m1<<"  e = "<<e1<<"  s = "<<s1
00142                            <<"  calib = "<<rawHit1[k1]->calibrationType()
00143                            <<"  new calib = "<<calib
00144                            <<"  oldE = "<<oldE<<" EADD = "<<edep_add
00145                            <<"  newE = "<<rawHit1[k1]->energy()<<endm;
00146                 }
00147         rawHit1[k1]->setCalibrationType(calib);
00148       }
00149       
00150       //Add remainig hits
00151       if(mEmbedAll) for(UInt_t k2=0;k2<rawHit2.size();k2++) if(checkHit(i,rawHit2[k2]))
00152       {
00153         if(stat_h2_all[k2]==0)
00154         { 
00155           StEmcRawHit* hit = new StEmcRawHit(id, rawHit2[k2]->module(),
00156                                              rawHit2[k2]->eta(), rawHit2[k2]->sub(),
00157                                              rawHit2[k2]->adc(), rawHit2[k2]->energy()); 
00158           detector1->addHit(hit); 
00159           new_edep_tot+=rawHit2[k2]->energy();
00160         }
00161       }   
00162     } 
00163     m_hit_1->Fill(old_edep_tot);
00164     m_hit_2->Fill(new_edep_tot);
00165     m_edep_1->Fill(edep1_tot);
00166     m_edep_2->Fill(edep2_tot);
00167   }
00168   return kStOK;
00169 }
00170 //-------------------------------------------------------------------
00174 void StEmcMixerMaker::clearPoints()
00175 {
00176   StEmcCollection* emccol=(StEmcCollection*)mEvent1->emcCollection();
00177   if(emccol)
00178   {
00179     StSPtrVecEmcPoint& pvec = emccol->barrelPoints();
00180     if(mClear)pvec.clear();  
00181   }
00182 }
00183 //-------------------------------------------------------------------
00187 void StEmcMixerMaker::clearClusters()
00188 {
00189   StEmcCollection* emccol=(StEmcCollection*)mEvent1->emcCollection();
00190   if(emccol) for(Int_t i=0; i<NDETECTORS; i++)
00191   {
00192     StDetectorId id = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00193     StEmcDetector* detector=emccol->detector(id);
00194     if(detector)
00195     {
00196       if(detector->cluster())
00197       {
00198         StSPtrVecEmcCluster& cluster=detector->cluster()->clusters();
00199         if(mClear)cluster.clear();  
00200       }
00201     }
00202   }
00203 }
00204 //-------------------------------------------------------------------
00212 Bool_t StEmcMixerMaker::getEvents()
00213 {
00214   // the default event should be the one where the data will be embedded
00215   // for this, it should be created before the second event
00216   // if there is no EMC collection in the event there is no reason for
00217   // embedding.
00218   mEvent1 = (StEvent*)GetInputDS("StEvent");
00219   if(!mEvent1) return kFALSE;
00220   if(!mEvent1->emcCollection()) return kFALSE; // should have EMC data
00221 
00222   // getting second event. First try StEmcSimulatorMaker
00223   // if there is no EMC simulator, try to get from a second StEvent in the memory
00224   // if the event source is the EMC simulator, StEmcMixerMaker owns the StEvent
00225   // that is created to hold the StEmcCollection from simulator
00226   // This maker adds the StEvent to .data so it will be properly deleted
00227   // during Clear() action
00228   StEmcSimulatorMaker *sim = (StEmcSimulatorMaker*)GetMaker("EmcSimulator");
00229   if(sim)
00230   {
00231     StEmcCollection *ecol = sim->getEmcCollection();
00232     if(ecol)
00233     {
00234       mEvent2 = new StEvent();
00235       mEvent2->setEmcCollection(ecol);
00236       AddData(mEvent2);
00237     } else { LOG_WARN <<"No second event to embed"<<endm; return kFALSE; }
00238   }
00239   else // no EMC simulator. Events come from another source
00240   {
00241     StMaker *m = GetMaker("embedIO");
00242         if(!m) { LOG_WARN <<"No embedIO maker"<<endm; return kFALSE; }
00243     mEvent2 = (StEvent*)m->GetInputDS("StEvent");
00244     if(!mEvent2) { LOG_WARN <<"No second event to embed"<<endm; return kFALSE; }
00245     if(!mEvent2->emcCollection()) { LOG_WARN <<"No second event to embed"<<endm; return kFALSE; }
00246   }
00247   
00248   LOG_DEBUG <<"Event 1 = "<<mEvent1<<"   Event 2 ="<<mEvent2<<endm;
00249   
00250   // StEvent pointers should be different.
00251   if(mEvent1==mEvent2) { LOG_DEBUG <<"Identical events"<<endm; return kFALSE; }
00252   return kTRUE;
00253 }
00254 //-------------------------------------------------------------------
00255 void StEmcMixerMaker::printHits(StEvent *event)
00256 {
00257   const TString detname[] = {"bemc", "bprs", "bsmde", "bsmdp","eemc", "eprs", "esmde", "esmdp"};
00258   StEmcCollection* emccol=(StEmcCollection*)event->emcCollection();
00259   
00260   for(Int_t i=0; i<NDETECTORS; i++)
00261   {  
00262     StDetectorId id = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00263     StEmcDetector* detector=emccol->detector(id);
00264     LOG_INFO <<"****************** hits in detector "<< detname[i].Data()<<endm;
00265     if(detector) for(UInt_t j=1;j<=NMODULES;j++) 
00266     {
00267       StEmcModule* module = detector->module(j);
00268       StSPtrVecEmcRawHit& rawHit=module->hits();
00269       if(rawHit.size()>0) { LOG_INFO <<"Number of hits for module "<<j<<" = "<<rawHit.size()<<endm; }
00270       for(UInt_t k=0;k<rawHit.size();k++) {
00271         LOG_INFO <<"Hit number = "<<k<<"  module = " << rawHit[k]->module()<<"  eta = "<<rawHit[k]->eta() << "  sub = "<< rawHit[k]->sub()<< "  adc = "<< rawHit[k]->adc() <<"  energy = "<<rawHit[k]->energy()<<endm;
00272           }
00273     }
00274   }
00275   
00276 }
00277 //-------------------------------------------------------------------
00282 void StEmcMixerMaker::getDB()
00283 {
00284   TDataSet *Db=GetDataBase("Calibrations/emc/y3bemc");
00285   Int_t valid[NDETECTORS]; for(Int_t i=0;i<NDETECTORS;i++) valid[i]=0;
00286   St_emcStatus* s0 = (St_emcStatus*)Db->Find("bemcStatus");
00287   if(s0)
00288   {
00289     emcStatus_st* st=s0->GetTable();
00290     for(Int_t i=0;i<EMCCHANNELS;i++)  mStatus[0][i] = st[0].Status[i];
00291   } else for(Int_t i=0;i<EMCCHANNELS;i++) mStatus[0][i] = 0;
00292    
00293   Db=GetDataBase("Calibrations/emc/y3bprs");
00294   St_emcStatus* s1 = (St_emcStatus*)Db->Find("bprsStatus");
00295   if(s1)
00296   {
00297     emcStatus_st* st=s1->GetTable();
00298     for(Int_t i=0;i<EMCCHANNELS;i++)  mStatus[1][i] = st[0].Status[i];
00299   } else for(Int_t i=0;i<EMCCHANNELS;i++) mStatus[1][i] = 0;
00300     
00301   Db=GetDataBase("Calibrations/emc/y3bsmde");
00302   St_smdStatus* s2 = (St_smdStatus*)Db->Find("bsmdeStatus");
00303   if(s2)
00304   {
00305     smdStatus_st* st=s2->GetTable();
00306     for(Int_t i=0;i<SMDCHANNELS;i++)  mStatus[2][i] = st[0].Status[i];
00307   } else for(Int_t i=0;i<SMDCHANNELS;i++) mStatus[2][i] = 0;
00308     
00309   Db=GetDataBase("Calibrations/emc/y3bsmdp");
00310   St_smdStatus* s3 = (St_smdStatus*)Db->Find("bsmdpStatus");
00311   if(s3)
00312   {
00313     smdStatus_st* st=s3->GetTable();
00314     for(Int_t i=0;i<SMDCHANNELS;i++)  mStatus[3][i] = st[0].Status[i];
00315   } else for(Int_t i=0;i<SMDCHANNELS;i++) mStatus[3][i] = 0;
00316   
00317   for(Int_t i=0;i<NDETECTORS;i++) 
00318     for(Int_t j=0;j<(EMCCHANNELS*(i<2)+SMDCHANNELS*(i>1));j++) 
00319       if(mStatus[i][j]==1) valid[i]++;
00320   
00321   LOG_DEBUG <<"Date = "<<GetDate()<<"  time = "<<GetTime()<<endm;
00322   for(Int_t i=0;i<NDETECTORS;i++) { LOG_DEBUG <<"Number of valid channels for detector "<<i<<" = "<<valid[i]<<endm; }
00323 }
00324 //-------------------------------------------------------------------
00331 Bool_t StEmcMixerMaker::checkHit(Int_t d,StEmcRawHit* h)
00332 {
00333   if(!h) return kFALSE;
00334   Int_t m = (Int_t)h->module();
00335   Int_t e = (Int_t)h->eta();
00336   Int_t s = abs(h->sub());     
00337   Int_t id;
00338   mGeom[d]->getId(m,e,s,id);
00339   if(id>0) if(mStatus[d][id-1]==1) return kTRUE;
00340   return kFALSE;     
00341 }
00342 //--------------------------------------------------------------------
00348 Int_t StEmcMixerMaker::addTracks()
00349 {
00350   // this just adds the tracks from one StEvent into the tracks of another StEvent.
00351   StSPtrVecTrackNode& Node1 = mEvent1->trackNodes();
00352   StSPtrVecTrackNode& Node2 = mEvent2->trackNodes();
00353   for (unsigned int j=0; j < Node2.size(); j++)
00354   {
00355     StGlobalTrack*   gTrack = (StGlobalTrack*)(Node2[j]->track(global));
00356     StPrimaryTrack*  pTrack = (StPrimaryTrack*)(Node2[j]->track(primary));
00357     StTrackNode *node = new StTrackNode();
00358     if(gTrack) node->addTrack(new StGlobalTrack(*gTrack));
00359     if(pTrack) node->addTrack(new StPrimaryTrack(*pTrack));
00360     Node1.push_back(node);
00361   } 
00362   return 0;
00363 }