StRoot: StFilterMaker/StGammaFilterMaker.cxx Source File

00001 // STL
00002 #include <vector>
00003 
00004 // ROOT Libraries
00005 #include "TVector3.h"
00006 #include "TFile.h"
00007 #include "TTree.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 // StMuDst Libraries
00030 #include "StBFChain/StBFChain.h"
00031 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00032 
00033 // St_g2t Libraries
00034 #include "tables/St_g2t_event_Table.h"
00035 #include "tables/St_g2t_pythia_Table.h"
00036 #include "tables/St_g2t_vertex_Table.h"
00037 #include "tables/St_particle_Table.h"
00038 #include "StSpinPool/StJetSkimEvent/StPythiaEvent.h"
00039 
00040 // Class Header
00041 #include "StGammaFilterMaker.h"
00042 
00043 ClassImp(StGammaFilterMaker)
00044 
00045 
00046 //                      Constructor                       //
00048 StGammaFilterMaker::StGammaFilterMaker(const char *name): 
00049 StMaker(name), mPythiaFile(0), mPythiaTree(0), mPythiaEvent(0)
00050 {
00051 
00052     LOG_DEBUG << "StGammaFilterMaker()" << endm;
00053 
00054     mFirst = true;
00055     mPythiaName = "pythia.root";
00056 
00057     // 2009 Defaults
00058     mSeedEnergyThreshold = 5.0;
00059     mClusterEtThreshold = 6.55;
00060 
00061     // Zero counters
00062     mTotal = 0;
00063     mAccepted = 0;
00064     
00065     // Calorimeter Geometries
00066     
00067     mMomentum = new TVector3();
00068         
00069     // BEMC Geometry   
00070     mBemcGeom = StEmcGeom::instance("bemc");
00071     assert(mBemcGeom);
00072     
00073     mPosition = new StEmcPosition();
00074     assert(mPosition);
00075     
00076 }
00077 
00079 //                       Destructor                       //
00081 StGammaFilterMaker::~StGammaFilterMaker()
00082 {
00083     LOG_DEBUG << "~StGammaFilterMaker()" << endm;
00084     
00085     // Clean up
00086     mMomentum->Delete();
00087     delete mPosition;
00088     
00089 }
00090 
00092 //                      Maker Init                        //
00094 Int_t StGammaFilterMaker::Init()
00095 {
00096 
00097     LOG_DEBUG << "Init()" << endm;
00098 
00099     // Force the BFC to listen to the filter
00100     SetAttr(".Privilege", 1);
00101     
00102     LOG_INFO  << "Init() - Using gamma filter on the BEMC" << endm;
00103     
00104     LOG_INFO << "Init() : Seed energy threshold = " << mSeedEnergyThreshold << " GeV " << endm;
00105     LOG_INFO << "Init() : Cluster eT threshold = " << mClusterEtThreshold << " GeV " << endm;
00106     
00107     // Create pythia output file
00108     StBFChain *bfChain = dynamic_cast<StBFChain*>(GetMakerInheritsFrom("StBFChain"));
00109     if(bfChain)
00110     {
00111         TString name(bfChain->GetFileOut());
00112         name.ReplaceAll(".root",".pythia.root");
00113         mPythiaName = name;
00114     }    
00115 
00116     LOG_INFO << "Init() : Storing pythia event information in " << mPythiaName << endm;
00117     mPythiaEvent = new StPythiaEvent;
00118     mPythiaFile = new TFile(mPythiaName, "recreate");
00119     mPythiaTree = new TTree("pythiaTree", "Pythia Tree");
00120     mPythiaTree->Branch("PythiaEvents", "StPythiaEvent", &mPythiaEvent);
00121 
00122     return kStOK;
00123     
00124 }
00125 
00127 //                     Maker Clear                        //
00129 void StGammaFilterMaker::Clear(const Option_t*)
00130 {
00131     LOG_DEBUG << "Clear()" << endm;
00132 
00133     // Clear the StEmcRawHit arrays
00134     StMcCalorimeterHit **point1 = mBemcTowerHits;
00135     for(unsigned int i = 0; i < nBemcTowers + 1; ++i)
00136     {
00137         *point1 = NULL;
00138         ++point1;
00139     }
00140 
00141     mPythiaEvent->Clear();
00142     
00143 }
00144 
00146 //                 Set filter thresholds                  //
00148 void StGammaFilterMaker::setThresholds(double seed, double cluster)
00149 {
00150 
00151     LOG_DEBUG << "setThresholds()" << endm;
00152     
00153     if(seed > 0) mSeedEnergyThreshold = seed;
00154     if(cluster > 0) mClusterEtThreshold = cluster;
00155     
00156 }
00157 
00159 //  Calculate multiplier between energy deposited in the  //
00160 //    BEMC sampling calorimeter and the true particle     //
00161 //   energy, i.e. the inverse of the sampling fraction,   //  
00162 //          as a function of tower pseudorapidity         //
00163 //                                                        //
00164 //        LOW_EM compatible coeffiecients taken from      //
00165 //   StEmcSimulatorMaker/StEmcSimpleSimulator.cxx v1.16   //
00167 
00168 double StGammaFilterMaker::BEMCSamplingMultiplier(double eta)
00169 {
00170 
00171     double x = fabs(eta);
00172     double c0 = 14.365;
00173     double c1 = -0.512;
00174     double c2 = 0.668;
00175     
00176     return c0 + c1 * x + c2 * x * x;
00177 
00178 }
00180 //                      Apply filter                      //
00182 Int_t StGammaFilterMaker::Make()
00183 {
00184 
00185     LOG_DEBUG << "Make()" << endm;
00186     
00187     // Acquire StEvent from the chain
00188     StEvent *mEvent = dynamic_cast<StEvent*>(GetDataSet("StEvent"));
00189     if(!mEvent)
00190     {
00191         LOG_WARN << "Make() - StEvent not found!" << endm;
00192         return kStWarn;
00193     }
00194     
00195     // Acquire StMcEvent from the chain
00196     StMcEvent *mMcEvent = dynamic_cast<StMcEvent*>(GetDataSet("StMcEvent"));
00197     
00198     ++mTotal;
00199 
00200     // Store the pythia event record
00201     fStorePythia();
00202 
00203     // Store the first event to ensure that a MuDst is created
00204     // in the case of no events passing the filter
00205     if(mFirst)
00206     {
00207 
00208         LOG_INFO << "Make() : Storing first event without applying filter" << endm;   
00209         
00210         ++mAccepted;
00211         mFirst = false;
00212 
00213         return kStOK;
00214     }    
00215 
00216     // Check for functional StMcEvent
00217     if(!mMcEvent)
00218     {
00219         LOG_WARN << "Make() - Bad StMcEvent!" << endm;
00220         return kStWarn;
00221     }
00222 
00223     // Acquire vertex information from the geant record
00224     // and abort the event if the vertex is too extreme
00225     double zVertex = mMcEvent->primaryVertex()->position().z();
00226 
00227     // Look for clusters in the chosen calorimeter
00228     Int_t status = 0;
00229     
00230     // Pick up geant calorimeter depositions from StMcEvent
00231     StMcEmcHitCollection *mcEmcCollection = mMcEvent->bemcHitCollection();
00232     assert(mcEmcCollection);
00233     status = makeBEMC(mcEmcCollection, zVertex);
00234     
00235     if(status == kStSKIP)
00236     {
00237         LOG_INFO << "Make() : Aborting Event " << mEvent->id() << " due to gamma filter rejection" << endm;
00238     }
00239     
00240     if(status == kStOK)
00241     {
00242         ++mAccepted;
00243         LOG_INFO << "Make() : Event " << mEvent->id() << " accepted!" << endm;
00244     }
00245     
00246     return status;
00247  
00248 }
00249 
00251 //       Construct crude BEMC clusters of 3x3 towers,     // 
00252 //        returning kStOK if a satisifactory cluster      // 
00253 //             is found and kStSKIP otherwise             //
00255 
00256 Int_t StGammaFilterMaker::makeBEMC(StMcEmcHitCollection *mcEmcCollection, double zVertex)
00257 {
00258 
00259     LOG_DEBUG << "makeBEMC()" << endm;
00260 
00262     //   Instantiate variables  //
00264     
00265     // Tower identification
00266     int bemcId = 0;
00267     
00268     unsigned int m = 0;
00269     unsigned int e = 0;
00270     unsigned int s = 0;
00271     
00272     float eta = 0;
00273 
00274     // Tower position
00275     float x = 0;
00276     float y = 0;
00277     float z = 0;
00278     double r = 0;
00279     
00281     //    Fill mBemcTowerHits   //
00282     //       with hit data      //
00284     for(unsigned int m = 1; m < mcEmcCollection->numberOfModules() + 1; ++m)
00285     {
00286         
00287         const StMcEmcModuleHitCollection *module = mcEmcCollection->module(m);
00288         const vector<StMcCalorimeterHit*> hits = module->detectorHits();
00289         
00290         for(unsigned int l = 0; l < hits.size(); ++l)
00291         {
00292             mBemcGeom->getId(hits[l]->module(), hits[l]->eta(), hits[l]->sub(), bemcId);
00293             mBemcTowerHits[bemcId] = hits[l];
00294         }
00295         
00296     }
00297     
00299     //    Search for clusters   //
00301     double seedEnergy = 0;
00302     double clusterEnergy = 0;
00303     double clusterEt = 0;
00304     unsigned int neighborId = 0;
00305     StMcCalorimeterHit *neighborHit = NULL;
00306     
00307     for(unsigned int n = 1; n < nBemcTowers + 1; ++n)
00308     {
00309         
00310         StMcCalorimeterHit *hit = mBemcTowerHits[n];
00311         if(!hit) continue;
00312         
00313         mBemcGeom->getEta(n, eta);
00314         
00315         seedEnergy = hit->dE() * BEMCSamplingMultiplier(eta);
00316         
00317         // Start with a seed tower...
00318         if(seedEnergy > mSeedEnergyThreshold)
00319         {
00320         
00321             clusterEnergy = seedEnergy;
00322         
00323             m = hit->module();
00324             e = hit->eta();
00325             s = hit->sub();
00326         
00327             // and accumulate energy from the neighboring towers
00328             for(int i = -1; i < 2; ++i)
00329             {
00330                 for(int j = -1; j < 2; ++j)
00331                 {
00332                     
00333                     if( !i && !j ) 
00334                     {
00335                         LOG_DEBUG << "\t\tTower " << n << ", E = " << seedEnergy << " (Seed)" << endm;
00336                         continue;
00337                     }
00338                     
00339                     neighborId = mPosition->getNextId(1, m, (int)e, s, i, j);
00340                     if(neighborId < 1) continue;
00341                     if(neighborId > 4800) continue;
00342                     
00343                     neighborHit = mBemcTowerHits[neighborId];
00344                     
00345                     if(!neighborHit) continue;
00346                     mBemcGeom->getEta(neighborId, eta);
00347                     LOG_DEBUG << "\t\tTower " << neighborId << ", E = " << neighborHit->dE() * BEMCSamplingMultiplier(eta) << endm;
00348                     clusterEnergy += neighborHit->dE() * BEMCSamplingMultiplier(eta);
00349                     
00350                 }
00351                 
00352             }
00353             
00354             // Project cluster energy into transverse plane
00355             mBemcGeom->getXYZ((int)n, x, y, z);
00356             z -= zVertex;
00357             r = sqrt(x * x + y * y + z * z);
00358             x *= clusterEnergy / r;
00359             y *= clusterEnergy / r;
00360             z *= clusterEnergy / r;
00361             mMomentum->SetXYZ(x, y, z);
00362             clusterEt = mMomentum->Perp();
00363 
00364             LOG_DEBUG << "\tCluster Energy = " << clusterEnergy << " GeV" << endm;
00365             LOG_DEBUG << "\tCluster eT  (corrected for vertex) = " << clusterEt << " GeV" << endm;
00366             LOG_DEBUG << "\tCluster eta (corrected for vertex) = " << mMomentum->Eta() << endm;
00367             LOG_DEBUG << "\tCluster phi = " << mMomentum->Phi() << endm;
00368             LOG_DEBUG << "" << endm;
00369             
00370             // Apply filter
00371             if(clusterEt > mClusterEtThreshold)
00372             {
00373                 return kStOk;
00374             }
00375             
00376         } // Seed tower check
00377         
00378     } // Seed tower loop
00379     
00380     return kStSKIP;
00381     
00382 }
00383  
00385 //                 Display filter stats                   //
00387 Int_t StGammaFilterMaker::Finish()
00388 {
00389 
00390     LOG_DEBUG << "::Finish()" << endm;
00391 
00392     mPythiaFile->Write();
00393     mPythiaFile->Close();
00394     
00395     LOG_INFO << "Finish() : " << GetName() << " finishing with " << endm;
00396     LOG_INFO << "Finish() : \t" << mAccepted << " of " << mTotal << " events passing the filter" << endm;
00397     
00398     return kStOK;
00399 
00400 }
00401   
00403 //                Store Pythia Event Record               //
00405 void StGammaFilterMaker::fStorePythia()
00406 {
00407 
00408     TDataSet* geant = GetDataSet("geant");
00409     
00410     if(geant) 
00411     {
00412     
00413         TDataSetIter iter(geant);
00414     
00415         // Global event information
00416         St_g2t_event* g2tEvent = (St_g2t_event*)iter("g2t_event");
00417         if(g2tEvent)
00418         {
00419 
00420             g2t_event_st* eventTable = (g2t_event_st*)g2tEvent->GetTable();
00421             if(eventTable) 
00422             {
00423                 mPythiaEvent->setRunId(eventTable->n_run);
00424                 mPythiaEvent->setEventId(eventTable->n_event);
00425             }
00426     
00427         }
00428 
00429         // Pythia event information
00430         St_g2t_pythia* pythiaEvent = (St_g2t_pythia*)iter("g2t_pythia");
00431         if(pythiaEvent) 
00432         {
00433         
00434             g2t_pythia_st* pythiaTable = (g2t_pythia_st*)pythiaEvent->GetTable();
00435             if(pythiaTable) 
00436             {
00437 
00438                 mPythiaEvent->setProcessId(pythiaTable->subprocess_id);
00439                 mPythiaEvent->setS(pythiaTable->mand_s);
00440                 mPythiaEvent->setT(pythiaTable->mand_t);
00441                 mPythiaEvent->setU(pythiaTable->mand_u);
00442                 mPythiaEvent->setPt(pythiaTable->hard_p);
00443                 mPythiaEvent->setCosTheta(pythiaTable->cos_th);
00444                 mPythiaEvent->setX1(pythiaTable->bjor_1);
00445                 mPythiaEvent->setX2(pythiaTable->bjor_2);
00446            
00447             }
00448     
00449         }
00450 
00451         // Vertex information
00452         St_g2t_vertex* vertexEvent = (St_g2t_vertex*)iter("g2t_vertex");
00453         if(vertexEvent) 
00454         {
00455      
00456             g2t_vertex_st* vertexTable = (g2t_vertex_st*)vertexEvent->GetTable();
00457             if(vertexTable) mPythiaEvent->setVertex(TVector3(vertexTable[0].ge_x));
00458         }
00459 
00460         // Pythia particle record
00461         St_particle* particleEvent = (St_particle*)iter("particle");
00462         if(particleEvent) 
00463         {
00464 
00465             particle_st* particleTable = (particle_st*)particleEvent->GetTable();
00466             if(particleTable) 
00467             {
00468 
00469                 for(int i = 0; i < particleEvent->GetNRows(); ++i) 
00470                 {
00471                     mPythiaEvent->addParticle(TParticle(particleTable[i].idhep,                  // pdg
00472                                                         particleTable[i].isthep,                 // status
00473                                                         particleTable[i].jmohep[0],              // mother1
00474                                                         particleTable[i].jmohep[1],              // mother2
00475                                                         particleTable[i].jdahep[0],              // daughter1
00476                                                         particleTable[i].jdahep[1],              // daughter2
00477                                                         TLorentzVector(particleTable[i].phep),   // momentum and energy
00478                                                         TLorentzVector(particleTable[i].vhep))); // production vertex and time
00479                 }
00480 
00481            }
00482 
00483        }
00484 
00485     }
00486 
00487     mPythiaTree->Fill();
00488 
00489     return;
00490 
00491 }