StRoot: StFilterMaker/StEemcGammaFilterMaker.cxx Source File

00001 // $Id: StEemcGammaFilterMaker.cxx,v 1.3 2010/08/09 21:51:22 seluzhen Exp $
00002 
00003 // STL
00004 #include <vector>
00005 
00006 // ROOT Libraries
00007 #include "TVector3.h"
00008 
00009 // StEvent Libraries
00010 #include "StEvent.h"
00011 
00012 // StEvent Libraries
00013 #include "StEventTypes.h"
00014 
00015 // StMcEvent Libraries
00016 #include "StMcEvent/StMcCalorimeterHit.hh"
00017 #include "StMcEvent/StMcEmcModuleHitCollection.hh"
00018 #include "StMcEvent/StMcEmcHitCollection.hh"
00019 #include "StMcEvent/StMcEvent.hh"
00020 #include "StMcEvent/StMcVertex.hh"
00021 
00022 // StEmc Libraries
00023 #include "StEmcClusterCollection.h"
00024 #include "StEmcPoint.h"
00025 #include "StEmcUtil/geometry/StEmcGeom.h"
00026 #include "StEmcUtil/projection/StEmcPosition.h"
00027 #include "StEmcUtil/others/emcDetectorName.h"
00028 
00029 // Endcap Libraries
00030 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00031 #include "StEEmcSimulatorMaker/StEEmcFastMaker.h"
00032 
00033 // Class Header
00034 #include "StEemcGammaFilterMaker.h"
00035 
00036 #include "tables/St_eemcGammaFilterMakerParams_Table.h"
00037 
00038 ClassImp(StEemcGammaFilterMaker)
00039 
00040 
00041 //                      Constructor                       //
00043 StEemcGammaFilterMaker::StEemcGammaFilterMaker():StMaker("eemcGammaFilterMaker.1.00") 
00044 {
00045   
00046   LOG_DEBUG << "StEemcGammaFilterMaker()" << endm;
00047   
00048   // Set reasonable defaults
00049   
00050   mMaxVertex = 120.0;
00051   
00052   mSeedEnergyThreshold = 3.8; //2.8; //3.4;//2.8;//2 sigma
00053   mClusterEtThreshold = 5.0; //3.8; //4.5;// 4.2;
00054   mEemcSamplingFraction = StEEmcFastMaker::getSamplingFraction();
00055   
00056   
00057   mMomentum = new TVector3();
00058   
00060   //       Calorimeter Geometry       //
00062   
00063   //EEMC Geometry
00064   mEemcGeom = new   EEmcGeomSimple();
00065   
00066   // Zero counters
00067   mTotal = 0;
00068   mAccepted = 0;
00069   mVertexRejected = 0;
00070   mFilterRejected = 0;
00071   //     mFilterMode = 1;
00072   //     mUseDbParams = true;
00073   mFilterMode = 1;
00074   mUseDbParams = false;
00075   
00076   if(!mFilterMode)
00077   {
00078     LOG_INFO <<" StEemcGammaFilterMaker is running in test mode.";
00079     LOG_INFO <<" set mFilterMode to 1 to run in analysis mode."<<endm;
00080   }
00081   
00082 }
00083 
00085 //                       Destructor                       //
00087 StEemcGammaFilterMaker::~StEemcGammaFilterMaker()
00088 {
00089   LOG_DEBUG << "~StEemcGammaFilterMaker()" << endm;
00090   
00091   // Clean up
00092   delete mMomentum;
00093   delete mEemcGeom;
00094   
00095 }
00096 
00098 //                      Maker Init                        //
00100 Int_t StEemcGammaFilterMaker::Init()
00101 {
00102   
00103   LOG_DEBUG << "Init()" << endm;
00104   
00105   // Force the BFC to listen to the filter
00106   SetAttr(".Privilege", 1);
00107   
00108   // Display maker information
00109   
00110   if (this->getUseDbParams())
00111   { // get parameter values from the database
00112   LOG_INFO << "Getting parameters from the database" << endm;
00113   TDataSet *DB = this->GetInputDB("Calibrations/emc");
00114   if (DB)
00115   {
00116     St_eemcGammaFilterMakerParams* params = (St_eemcGammaFilterMakerParams*)DB->Find("eemcGammaFilterMakerParams");
00117     if (params)
00118     {
00119       eemcGammaFilterMakerParams_st *p = params->GetTable();
00120       if (p)
00121       {
00122         this->mEemcSamplingFraction = p->eemcSamplingFraction;
00123         this->mSeedEnergyThreshold = p->seedEnergyThreshold;
00124         this->mClusterEtThreshold = p->clusterEtThreshold;
00125         this->mMaxVertex = p->maxVertexZ;
00126         this->mFilterMode = p->filterMode;
00127         LOG_DEBUG << "Successfully read parameters from the database table" << endm;
00128       } 
00129       else
00130       {
00131         LOG_ERROR << "The database table has no entry!" << endl;
00132         return kStFATAL;
00133       }
00134     } 
00135     else
00136     {
00137       LOG_ERROR << "Cannot find the table in the database!" << endm;
00138       return kStFATAL;
00139     }
00140   } 
00141   else
00142   {
00143     LOG_ERROR << "Cannot read the database!" << endm;
00144     return kStFATAL;
00145   }
00146   }//done getting db parameters
00147   
00148   LOG_INFO << "StEEmcDbMaker::Init() : Using gamma filter on the EEMC" << endm;
00149   LOG_INFO << "StEEmcDbMaker::Init() : EEMC Sampling Fraction = " << mEemcSamplingFraction << endm;
00150   
00151   LOG_INFO << "StEEmcDbMaker::Init() : Seed energy threshold = " << mSeedEnergyThreshold << " GeV " << endm;
00152   LOG_INFO << "StEEmcDbMaker::Init() : Cluster eT threshold = " << mClusterEtThreshold << " GeV " << endm;
00153   LOG_INFO << "StEEmcDbMaker::Init() : Maximum vertex = +/- " << mMaxVertex << " cm" << endm;
00154   if (!mFilterMode)
00155     LOG_INFO << "StEEmcDbMaker::Init() : Running the TEST mode (accepting all events). Set mFilterMode=1 to actually reject events in BFC"<< endm;
00156   
00157   
00158   
00159   // If using the EEMC ensure that a sampling 
00160   // fraction has been entered by the user 
00161   if (mEemcSamplingFraction != StEEmcFastMaker::getSamplingFraction())
00162   {
00163     LOG_WARN << "StEEmcDbMaker::Init() : EEMC Sampling fraction (" << mEemcSamplingFraction << ") is different from what is in the Fast Simulator (" << StEEmcFastMaker::getSamplingFraction() <<")" << endm;
00164   }
00165   
00166   
00167   return kStOK;
00168   
00169 }
00170 
00172 //                     Maker Clear                        //
00174 void StEemcGammaFilterMaker::Clear(const Option_t*)
00175 {
00176   LOG_DEBUG << "Clear()" << endm;
00177   
00178   // Clear the StEmcRawHit arrays
00179   
00180   double *point2 = mEemcTowerHits;
00181   for(unsigned int i = 0; i < nEemcTowers; ++i)
00182   {
00183     *point2 = 0;
00184     ++point2;
00185   }
00186   
00187   
00188   
00189 }
00190 
00191 
00192 
00194 //       Set maximum vertex allowed by the filter         //
00196 void StEemcGammaFilterMaker::setMaxVertex(double maxVertex)
00197 {
00198   LOG_DEBUG << "setMaxVertex()" << endm;
00199   mMaxVertex = maxVertex;
00200 }
00201 
00203 //                 Set filter thresholds                  //
00205 void StEemcGammaFilterMaker::setThresholds(double seed, double cluster)
00206 {
00207   
00208   LOG_DEBUG << "setThresholds()" << endm;
00209   
00210   if(seed > 0) mSeedEnergyThreshold = seed;
00211   if(cluster > 0) mClusterEtThreshold = cluster;
00212   
00213 }
00214 
00216 //               Set EEMC sampling fraction               //
00218 void StEemcGammaFilterMaker::setEEMCSamplingFraction(double f)
00219 {
00220   mEemcSamplingFraction = f;
00221 }
00222 
00223 
00225 //                      Apply filter                      //
00227 Int_t StEemcGammaFilterMaker::Make()
00228 {
00229   
00230   LOG_DEBUG << "Make()" << endm;
00231   
00232   // Acquire StEvent from the chain
00233   StEvent *mEvent = static_cast<StEvent*>(GetDataSet("StEvent"));
00234   
00235   // Acquire StMcEvent from the chain
00236   StMcEvent *mMcEvent = static_cast<StMcEvent*>(GetDataSet("StMcEvent"));
00237   
00238   ++mTotal;
00239   
00240   // Store the first event to ensure that a MuDst is created
00241   // in the case of no events passing the filter
00242   //     if(mEvent->id() == 1)
00243   //     {
00244     //         LOG_INFO << "Make() : Storing first event without applying filter" << endm;   
00245     //         
00246     //         ++mAccepted;
00247     //         return kStOK;
00248     //     }
00249     
00250     // Check for functional StMcEvent
00251     if(!mMcEvent)
00252     {
00253       LOG_WARN << "Reject() : Bad StMcEvent!" << endm;
00254       return kStWarn;
00255     }
00256     
00257     // Acquire vertex information from the geant record
00258     // and abort the event if the vertex is too extreme
00259     double zVertex = mMcEvent->primaryVertex()->position().z();
00260     
00261     if(fabs(zVertex) > mMaxVertex)
00262     {
00263       LOG_INFO << "Reject() : Aborting Event " << mEvent->id()
00264       << " due to extreme geant vertex of " << zVertex << " cm " << endm;
00265       
00266       ++mVertexRejected;
00267       if(mFilterMode)
00268         return kStSKIP;
00269       else
00270         return kStOK;
00271       
00272     }
00273     
00274     // Look for clusters in the chosen calorimeter
00275     Int_t status = 0;
00276     
00277     // Pick up calorimeter information from StEvent
00278     if (mEvent) {
00279       StEmcCollection *emcCollection = mEvent->emcCollection();
00280       status = makeEEMC(emcCollection, zVertex);
00281     }
00282     
00283     if(status == kStSKIP)
00284     {
00285       ++mFilterRejected;
00286       LOG_INFO << "Make() : Aborting Event " << mEvent->id() << " due to gamma filter rejection" << endm;
00287       LOG_DEBUG << "Make() : Vertex = " << zVertex << " (cm)" << endm;
00288     }
00289     
00290     if(status == kStOK)
00291     {
00292       ++mAccepted;
00293       LOG_INFO << "Make() : Event " << mEvent->id() << " accepted!" << endm;
00294     }
00295     
00296     if (mFilterMode==0) status = kStOK;
00297     return status;
00298     
00299 }
00300 
00301 
00303 //       Construct crude EEMC clusters of 3x3 towers,     // 
00304 //        returning kStOK if a satisifactory cluster      // 
00305 //             is found and kStSKIP otherwise             //
00307 Int_t StEemcGammaFilterMaker::makeEEMC(StEmcCollection *emcCollection, double zVertex)
00308 {
00309   
00310   LOG_DEBUG << "makeEEMC()" << endm;
00311   
00313   //   Instantiate variables  //
00315   
00316   unsigned int nSub = 5;
00317   int eta = 0;
00318   int phi = 0;
00319   
00320   double clusterEnergy = 0;
00321   double clusterEt = 0;
00322   double seedEnergy = 0;
00323   
00324   int neighborEta = 0;
00325   int neighborPhi = 0;
00326   
00327   // Access EEMC information
00328   StEmcDetector* eemc = emcCollection ? emcCollection->detector(kEndcapEmcTowerId) : 0;  
00329   if(!eemc) return kStSKIP;
00330   
00332   //    Fill mEemcTowerHits   //
00333   //       with hit data      //
00335   for(unsigned int m = 1; m < eemc->numberOfModules() + 1; ++m)
00336   {
00337     
00338     StEmcModule *module = eemc->module(m);
00339     if (module) {
00340       StSPtrVecEmcRawHit &rawHits = module->hits();
00341       
00342       for(unsigned int l = 0 ; l < rawHits.size(); ++l)
00343       {
00344         
00345         StEmcRawHit *tempHit = rawHits[l];
00346         if (tempHit) {
00347           eta = tempHit->eta() - 1;
00348           phi = (m - 1) * nSub + tempHit->sub() - 1;
00349           mEemcTowerHits[eta * nEemcPhiBins + phi] = tempHit->energy() / mEemcSamplingFraction;
00350         }
00351       }
00352     }
00353     
00354   }
00355   
00357   //    Search for clusters   //
00359   
00360   for(unsigned int e = 0; e < nEemcEtaBins; ++e)
00361   {
00362     
00363     for(unsigned int p = 0; p < nEemcPhiBins; ++p)
00364     {
00365       
00366       clusterEnergy = 0;
00367       
00368       seedEnergy = mEemcTowerHits[e * nEemcPhiBins + p];
00369       
00370       // Start with a seed tower...
00371       if(seedEnergy > mSeedEnergyThreshold)
00372       {
00373         
00374         clusterEnergy = seedEnergy;
00375         
00376         // and accumulate energy from the neighboring towers
00377         for(int i = -1; i < 2; ++i)
00378         {
00379           
00380           for(int j = -1; j < 2; ++j)
00381           {
00382             
00383             if( !i && !j ) 
00384             {
00385               LOG_DEBUG << "\t\tTower (eta, phi) = (" << e << ", " << p 
00386               << "), E = " << seedEnergy << " (Seed)" << endm;
00387               continue;
00388             }
00389             
00390             neighborEta = e + i;
00391             if(neighborEta < 0) continue;
00392             if(neighborEta > int(nEemcEtaBins - 1)) continue;
00393             
00394             neighborPhi = (p + j) % nEemcPhiBins;
00395             if(neighborPhi == - 1) neighborPhi = nEemcPhiBins - 1;
00396             LOG_DEBUG << "\t\tTower (eta, phi) = (" << neighborEta << ", " << neighborPhi 
00397             << "), E = " << mEemcTowerHits[neighborEta * nEemcPhiBins + neighborPhi] << endm;
00398             clusterEnergy += mEemcTowerHits[neighborEta * nEemcPhiBins + neighborPhi];
00399             
00400           }
00401           
00402         }
00403         
00404         // Project cluster energy into transverse plane
00405         if (mMomentum && mEemcGeom) {
00406           *mMomentum = mEemcGeom->getTowerCenter(p / nSub, p % nSub, e);
00407           mMomentum->SetZ(mMomentum->Z() - zVertex);
00408           mMomentum->SetX(mMomentum->X() * clusterEnergy / mMomentum->Mag() );
00409           mMomentum->SetY(mMomentum->Y() * clusterEnergy / mMomentum->Mag() );
00410           mMomentum->SetZ(mMomentum->Z() * clusterEnergy / mMomentum->Mag() );
00411           clusterEt = mMomentum->Perp();
00412           
00413           LOG_DEBUG << "\tCluster Energy = " << clusterEnergy << " GeV" << endm;
00414           LOG_DEBUG << "\tCluster eT  (corrected for vertex) = " << clusterEt << " GeV" << endm;
00415           LOG_DEBUG << "\tCluster eta (corrected for vertex) = " << mMomentum->Eta() << endm;
00416           LOG_DEBUG << "\tCluster phi = " << mMomentum->Phi() << endm;
00417           LOG_DEBUG << "" << endm;
00418         }
00419         
00420         // Apply filter
00421         if(clusterEt > mClusterEtThreshold)
00422         {
00423           return kStOK;
00424         }
00425         
00426         
00427       } // Seed tower check
00428       
00429     } // Seed tower phi loop
00430     
00431   } // Seed tower eta loop
00432   
00433   return kStSKIP;
00434   
00435 }
00436 
00438 //                 Display filter stats                   //
00440 Int_t StEemcGammaFilterMaker::Finish()
00441 {
00442   LOG_DEBUG << "::Finish()" << endm;
00443   
00444   LOG_INFO << "Finish() : " << GetName() << " finishing with " << endm;
00445   LOG_INFO << "Finish() : \t" << mAccepted << " of " << mTotal << " events passing the filter" << endm;
00446   LOG_INFO << "Finish() : \t" << mVertexRejected << " rejected for bad vertex" << endm;
00447   LOG_INFO << "Finish() : \t" << mFilterRejected << " rejected for no clusters" << endm;
00448   
00449   return kStOK;
00450 }
00451 
00452 // $Log: StEemcGammaFilterMaker.cxx,v $
00453 // Revision 1.3  2010/08/09 21:51:22  seluzhen
00454 // updated comment field
00455 //