pi0Efficiency/src/pi0Efficiency.cc

0001 #include "pi0Efficiency.h"
0002
0003 //Fun4all stuff
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <fun4all/Fun4AllHistoManager.h>
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009 #include <phool/phool.h>
0010 #include <ffaobjects/EventHeader.h>
0011
0012 //ROOT stuff
0013 #include <TH1F.h>
0014 #include <TH2F.h>
0015 #include <TH3F.h>
0016 #include <TFile.h>
0017 #include <TLorentzVector.h>
0018
0019 //for emc clusters
0020 #include <calobase/RawCluster.h>
0021 #include <calobase/RawClusterContainer.h>
0022 #include <calobase/RawClusterUtility.h>
0023 #include <calobase/RawTowerGeomContainer.h>
0024 #include <calobase/RawTower.h>
0025 #include <calobase/RawTowerContainer.h>
0026
0027 //caloEvalStack for cluster to truth matching
0028 #include <g4eval/CaloEvalStack.h>
0029 #include <g4eval/CaloRawClusterEval.h>
0030
0031 //for vetex information
0032 #include <g4vertex/GlobalVertex.h>
0033 #include <g4vertex/GlobalVertexMap.h>
0034
0035 //tracking
0036 #include <trackbase_historic/SvtxTrack.h>
0037 #include <trackbase_historic/SvtxTrackMap.h>
0038 #include <trackbase_historic/SvtxVertex.h>
0039 #include <trackbase_historic/SvtxVertexMap.h>
0040
0041 //truth information
0042 #include <g4main/PHG4TruthInfoContainer.h>
0043 #include <g4main/PHG4Particle.h>
0044 #include </gpfs/mnt/gpfs02/sphenix/user/ahodges/macros_git/coresoftware/generators/phhepmc/PHHepMCGenEvent.h>
0045 #include </gpfs/mnt/gpfs02/sphenix/user/ahodges/macros_git/coresoftware/generators/phhepmc/PHHepMCGenEventMap.h>
0046 #include <g4main/PHG4VtxPoint.h>
0047 #include <HepMC/GenEvent.h>
0048 #include <HepMC/GenParticle.h>
0049 #include <HepMC/GenVertex.h>
0050 #include <HepMC/IteratorRange.h>
0051 #include <HepMC/SimpleVector.h>
0052 #include <HepMC/GenParticle.h>
0053 #include <CLHEP/Geometry/Point3D.h>
0054
0055 //____________________________________________________________________________..
0056 pi0Efficiency::pi0Efficiency(const std::string &name, const std::string &outName):
0057 SubsysReco(name),
0058   OutFile(outName)
0059 {
0060   std::cout << "pi0Efficiency::pi0Efficiency(const std::string &name) Calling ctor" << std::endl;
0061
0062
0063   for(int j = 0; j < nEtaBins; j++) ePi0InvMassEcut[j] = 0;
0064
0065   photonE = 0;
0066
0067   for(int i = 0; i < nEtaBins; i++)prim_pi0_E[i] = 0;
0068
0069   clusterE = 0;
0070
0071   truthPi0EAsym = 0;
0072
0073   truthPi0EDeltaR = 0;
0074
0075   hMassRat = 0;
0076
0077   pi0EScale = 0;
0078 }
0079
0080 //____________________________________________________________________________..
0081 pi0Efficiency::~pi0Efficiency()
0082 {
0083   std::cout << "pi0Efficiency::~pi0Efficiency() Calling dtor" << std::endl;
0084 }
0085
0086 //____________________________________________________________________________..
0087 int pi0Efficiency::Init(PHCompositeNode *topNode)
0088 {
0089   std::cout << "pi0Efficiency::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0090
0091
0092
0093   for(int j = 0; j < nEtaBins; j++)
0094     {
0095       ePi0InvMassEcut[j] = new TH3F(Form("ePi0InvMassEcut_Eta%d",j),Form("Pi0 energy vs. Inv Mass vs. Min pho Energy Eta%d", j), 500,0,50,500,-0.1,1,200,0,20);
0096     }
0097
0098   clusterE = new TH1F("clusterE","Cluster Energy",200,0,20);
0099
0100   for(int i = 0; i < nEtaBins; i++)prim_pi0_E[i] = new TH1F(Form("primPi0E_Eta%d",i),"Primary pi0 Energy",200,0,20);
0101
0102   photonE = new TH1F("photonE","Decay Photon Energy",200,0,20);
0103
0104   pi0EScale = new TH2F("pi0EScale","Pi0 energy fraction",200,0,2,200,0,20);
0105
0106   truthPi0EDeltaR = new TH2F("truthPi0EDeltaR","truth pi0 energy versus decay opening angle",200,0,20,100,0,.5);
0107
0108   truthPi0EAsym = new TH2F("truthPi0EAsym","truth pi0 energy vs. decay energy asymmetry",200,0,20,200,0,1);
0109
0110   hMassRat = new TH1F("hMassRat","ratio of pi0 mass reco from truth photons to primary mass",200,0,2);
0111
0112   Fun4AllServer *se = Fun4AllServer::instance();
0113   se -> Print("NODETREE");
0114   hm = new Fun4AllHistoManager("MYHISTOS");
0115
0116   se -> registerHistoManager(hm);
0117
0118   se -> registerHisto(clusterE -> GetName(), clusterE);
0119
0120   for(int i = 0; i < nEtaBins; i++)se -> registerHisto(prim_pi0_E[i] -> GetName(), prim_pi0_E[i]);
0121
0122   se -> registerHisto(photonE -> GetName(), photonE);
0123
0124   for(int j = 0; j < nEtaBins; j++) se -> registerHisto(ePi0InvMassEcut[j] -> GetName(), ePi0InvMassEcut[j]);
0125
0126   se -> registerHisto(truthPi0EAsym->GetName(), truthPi0EAsym);
0127
0128   se -> registerHisto(hMassRat -> GetName(), hMassRat);
0129
0130   se -> registerHisto(truthPi0EDeltaR -> GetName(), truthPi0EDeltaR);
0131
0132   out = new TFile(OutFile.c_str(),"RECREATE");
0133
0134   return Fun4AllReturnCodes::EVENT_OK;
0135 }
0136
0137 //____________________________________________________________________________..
0138 int pi0Efficiency::InitRun(PHCompositeNode *topNode)
0139 {
0140   std::cout << "pi0Efficiency::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0141   return Fun4AllReturnCodes::EVENT_OK;
0142 }
0143
0144 //____________________________________________________________________________..
0145 int pi0Efficiency::process_event(PHCompositeNode *topNode)
0146 {
0147   std::cout << "pi0Efficiency::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0148
0149
0150   //Information on clusters
0151   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0152   //RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0153   if(!clusterContainer)
0154     {
0155       std::cout << PHWHERE << "pi0Efficiency::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
0156
0157       return 0;
0158     }
0159
0160
0161   //Vertex information
0162   GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
0163   if (!vtxContainer)
0164     {
0165       std::cout << PHWHERE << "pi0Efficiency::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0166       assert(vtxContainer);  // force quit
0167
0168       return 0;
0169     }
0170
0171   if (vtxContainer->empty())
0172     {
0173       std::cout << PHWHERE << "pi0Efficiency::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0174       return 0;
0175     }
0176
0177   //More vertex information
0178   GlobalVertex *vtx = vtxContainer->begin()->second;
0179   if(!vtx)
0180     {
0181
0182       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find vtx from vtxContainer"  << std::endl;
0183       return Fun4AllReturnCodes::ABORTEVENT;
0184     }
0185
0186   //Tower geometry node for location information
0187   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0188   if (!towergeom)
0189     {
0190       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node TOWERGEOM_CEMC"  << std::endl;
0191       return Fun4AllReturnCodes::ABORTEVENT;
0192     }
0193
0194   //Raw tower information
0195   RawTowerContainer *towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC");
0196   if(!towerContainer)
0197     {
0198       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node TOWER_CALIB_CEMC"  << std::endl;
0199       return Fun4AllReturnCodes::ABORTEVENT;
0200     }
0201
0202   //truth particle information
0203   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0204   if(!truthinfo)
0205     {
0206       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node G4TruthInfo"  << std::endl;
0207       return Fun4AllReturnCodes::ABORTEVENT;
0208     }
0209
0210   int firstphotonflag = 0;
0211   int firstfirstphotonflag = 0;
0212   int secondphotonflag = 0;
0213   int secondsecondphotonflag = 0;
0214
0215   PHG4TruthInfoContainer::Range truthRangeDecay1 = truthinfo->GetSecondaryParticleRange();
0216   PHG4TruthInfoContainer::ConstIterator truthIterDecay1;
0217
0218   float photonEtaMax = 1.1;
0219   float mesonEtaMax = 0.3;
0220   TLorentzVector photon1;
0221   TLorentzVector photon2;
0222
0223   for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
0224     {
0225
0226       PHG4Particle *decay = truthIterDecay1 -> second;
0227
0228       int dumtruthpid = decay -> get_pid();
0229       int dumparentid = decay -> get_parent_id();
0230
0231       //if the parent is the pi0 and it's a photon and we haven't marked one yet
0232       if(dumparentid == 1 && dumtruthpid == 22 && !firstphotonflag)
0233     {
0234       if(abs(getEta(decay)) > photonEtaMax)
0235         {
0236           std::cout << "Photon 1 outside acceptance " << std::endl;
0237           return 0; //decay product outside acceptance.
0238         }
0239       photon1.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e());
0240
0241       firstphotonflag = 1;
0242     }
0243
0244       if(dumparentid == 1 && dumtruthpid == 22 && firstphotonflag && firstfirstphotonflag)
0245     {
0246       if(abs(getEta(decay)) > photonEtaMax)
0247         {
0248           std::cout << "Photon 2 outside acceptance " << std::endl;
0249           return 0; //decay product outside acceptance
0250         }
0251       photon2.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e()) ;
0252
0253       secondphotonflag = 1;
0254     }
0255
0256       //deal with dalitz decays
0257       if(dumparentid == 1 && firstphotonflag && secondphotonflag && secondsecondphotonflag)
0258     {
0259       std::cout << "Dalitz decay, skipping event" << std::endl;
0260       return 0;
0261     }
0262
0263
0264
0265       //Need these extra flags, otherwise the dalitz
0266       //check will always have first and second photon
0267       //flags marked true, so this allows the dalitz check
0268       //to occur after marking the second truth photon
0269       //as a decay component
0270       if(firstphotonflag) firstfirstphotonflag = 1;
0271       if(secondphotonflag) secondsecondphotonflag = 1;
0272     }
0273   photonE -> Fill(photon1.Energy());
0274   photonE -> Fill(photon2.Energy());
0275   float asym = abs(photon1.Energy() - photon2.Energy())/(photon1.Energy() + photon2.Energy());
0276   float deltaR = photon1.DeltaR(photon2);
0277
0278   //Now we go and find all our truth pi0s
0279   PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange();
0280   PHG4TruthInfoContainer::ConstIterator truthIter;
0281   std::vector<int> mPi0Barcode;
0282   Double_t pi0Mass = 0;
0283   PHG4Particle *truthPar = NULL;
0284   TLorentzVector truePi0;
0285
0286   for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
0287     {
0288        truthPar = truthIter->second;
0289
0290       if(truthPar -> get_pid() == 111 && truthPar -> get_parent_id() == 0)//for single particle gun, this list has one particle, the
0291                                                                           //one of interest, but best be safe out of habit
0292     {
0293       if(getEta(truthPar) >= mesonEtaMax)
0294         {
0295           std::cout << "Parent outside allowed spatial limit" << std::endl;
0296           return 0;
0297         }
0298       int etaBin = -1;
0299       if(photon1.Energy() >= photon2.Energy())etaBin = getEtaBin(photon1.PseudoRapidity());
0300       else etaBin = getEtaBin(photon2.PseudoRapidity());
0301       if(etaBin >= 0)prim_pi0_E[etaBin] -> Fill(truthPar -> get_e());
0302       prim_pi0_E[nEtaBins-1] -> Fill(truthPar -> get_e());
0303
0304       mPi0Barcode.push_back(truthPar -> get_barcode());
0305       truePi0.SetPxPyPzE(truthPar -> get_px(), truthPar -> get_py(), truthPar -> get_pz(), truthPar -> get_e());
0306       pi0Mass = truePi0.M();//sanity check later to make sure when we reconstruct the pi0 from truth photons we get the pi0 mass back
0307
0308     }
0309
0310     }
0311
0312   truthPi0EDeltaR -> Fill(truePi0.Energy(), deltaR);
0313   truthPi0EAsym -> Fill(truePi0.Energy(), asym);
0314
0315   //iterate over items in the container.
0316   RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
0317   RawClusterContainer::ConstIterator clusterIter;
0318
0319   //remove noise with a 300MeV energy cut and store
0320   //good clusters in a vector for later.
0321   std::vector<RawCluster*> goodRecoCluster;
0322   float minE = 0.3;
0323   CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0324
0325   for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0326     {
0327       RawCluster *recoCluster = clusterIter -> second;
0328
0329       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0330       float clusE = E_vec_cluster.mag();
0331       clusterE -> Fill(clusE);
0332       if(clusE < minE) continue;
0333
0334       goodRecoCluster.push_back(recoCluster);
0335     }
0336
0337   for(int i = 0; i < (int)goodRecoCluster.size(); ++i)
0338     {
0339       //grab the first good cluster.
0340       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0341       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*goodRecoCluster[i], vertex);
0342
0343       TLorentzVector pho1, pho2, pi0;
0344       pho1.SetPtEtaPhiE(E_vec_cluster.perp(), E_vec_cluster.pseudoRapidity(), E_vec_cluster.getPhi(),  E_vec_cluster.mag());
0345
0346        for(int j = 0; j <(int)goodRecoCluster.size(); ++j)
0347     {
0348
0349       if(j <= i) continue; //make sure we don't get redundant pairs
0350       CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*goodRecoCluster[j], vertex);
0351
0352       pho2.SetPtEtaPhiE(E_vec_cluster2.perp(), E_vec_cluster2.pseudoRapidity(), E_vec_cluster2.getPhi(),  E_vec_cluster2.mag());
0353
0354       pi0 = pho1 + pho2;
0355
0356       pi0EScale -> Fill(pi0.Energy()/truePi0.Energy(), truePi0.Energy());//check the pi0 energy scale.
0357
0358
0359       //isolating spatial kinematics to make sure we're really measuring
0360       //sPHENIX's ability to reconstruct pi0's within its volume
0361       if(abs(pho1.PseudoRapidity()) < photonEtaMax && abs(pho2.PseudoRapidity()) < photonEtaMax && abs(pi0.PseudoRapidity()) < mesonEtaMax)
0362         {
0363           int etaBin = -1;
0364           //determine psuedorapidity bin from lead photon
0365           if(pho1.Energy() >= pho2.Energy()) etaBin = getEtaBin(pho1.PseudoRapidity());
0366           else etaBin = getEtaBin(pho2.PseudoRapidity());
0367           if(etaBin >= 0)ePi0InvMassEcut[etaBin] -> Fill(pi0.Energy(), pi0.M(),std::min(pho1.Energy(), pho2.Energy()));
0368           ePi0InvMassEcut[nEtaBins-1] -> Fill(pi0.Energy(), pi0.M(),std::min(pho1.Energy(), pho2.Energy()));
0369         }
0370     }
0371     }
0372
0373
0374
0375
0376
0377   //make sure we aren't messing up the reconstruction :)
0378   TLorentzVector pi0fromTruPhoton = photon1 + photon2;
0379   hMassRat -> Fill(pi0fromTruPhoton.M()/pi0Mass);
0380   //std::cout << "Mass diff " << massdiff << " from truth reco pi0: " << pi0fromTruPhoton.M() << " and truth pi0 mass: " << pi0Mass << std::endl;
0381
0382   /*This doesn't work on things that aren't primary particles.
0383   //and here's where we do the sneaky business of finding out where the photons went.
0384   if(!cluster1 || !cluster2)
0385     {
0386       if(!cluster1)
0387     {
0388
0389       unmatchedLocale -> Fill(photon1.PseudoRapidity(), photon1.Phi(), truthPar -> get_e());
0390       unmatchedE -> Fill(photon1.E());
0391     }
0392       if(!cluster2)
0393     {
0394
0395       unmatchedLocale -> Fill(photon2.PseudoRapidity(), photon2.Phi(), truthPar -> get_e());
0396       unmatchedE -> Fill(photon2.E());
0397     }
0398       ePi0InvMassEcut[1] -> Fill(truthPar -> get_e(), 0., std::min(photon1.Energy(), photon2.Energy()));
0399     }
0400   else
0401     {
0402
0403       if(cluster1 -> get_id() == cluster2 -> get_id()) return 0;
0404
0405       TLorentzVector cpi0, c1, c2;
0406
0407       CLHEP::Hep3Vector eVecCluster1 = RawClusterUtility::GetECoreVec(*cluster1, vertex);
0408       CLHEP::Hep3Vector eVecCluster2 = RawClusterUtility::GetECoreVec(*cluster2, vertex);
0409
0410       c1.SetPtEtaPhiE(eVecCluster1.perp(), eVecCluster1.pseudoRapidity(), eVecCluster1.getPhi(), eVecCluster1.mag());
0411       c2.SetPtEtaPhiE(eVecCluster2.perp(), eVecCluster2.pseudoRapidity(), eVecCluster2.getPhi(), eVecCluster2.mag());
0412
0413       cpi0 = c1 + c2;
0414
0415       ePi0InvMassEcut[1] -> Fill(cpi0.Energy(), cpi0.M(),std::min(c1.Energy(), c1.Energy()));
0416
0417       invMassEPhi -> Fill(cpi0.M(), truthPar -> get_e(), cpi0.Phi());
0418       invMassPhoPi0 -> Fill(cpi0.M(), cpi0.Energy()/truthPar -> get_e(), c1.Energy()/photon1.Energy());
0419       invMassPhoPi0 -> Fill(cpi0.M(), cpi0.Energy()/truthPar -> get_e(), c1.Energy()/photon2.Energy());
0420       }*/
0421   /*Supplanting with method developed by Joe Osborn above
0422   PHG4TruthInfoContainer::Range truthRangeDecay1 = truthinfo->GetSecondaryParticleRange();
0423   PHG4TruthInfoContainer::ConstIterator truthIterDecay1;
0424
0425
0426   //Collect truth decay photons
0427   std::vector<PHG4Particle*> pi0Pho;
0428   std::vector<int> pi0PhoBarcode;
0429   for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
0430     {
0431       PHG4Particle *truthDecay1 = truthIterDecay1 -> second;
0432       PHG4Particle *mother = truthinfo -> GetParticle(truthDecay1 -> get_parent_id());
0433       int mBarcode = mother -> get_barcode();
0434
0435       if(truthDecay1 -> get_parent_id() != 0 && truthDecay1 -> get_pid() == 22 && getEta(truthDecay1) <= etaMax && (std::find(mPi0Barcode.begin(), mPi0Barcode.end(),mBarcode) != mPi0Barcode.end())) //want decay photons w/ strict matching back to a pi0
0436     {
0437       pi0PhoBarcode.push_back(mBarcode);//this is for the actual pairing process later where we want to match two photons together.
0438       photonE -> Fill(truthDecay1 -> get_e());
0439       pi0Pho.push_back(truthDecay1);
0440
0441
0442     }
0443     }
0444
0445   //Re-pair decay photons to extract their kinematic information, match to clusters, and hopefully
0446   //see what causes inefficiencies. Apparently GEANT considers literally every particle
0447   //created after the initial generation a secondary particle, you need to do this pairing and then
0448   //find the pair with the invariant mass closes to the pi0's (~O(10-9) difference)
0449   float massdiff = 10;
0450
0451   for(int i = 0; i < (int)pi0Pho.size(); i++)
0452     {
0453       for(int j = 0; j < (int)pi0Pho.size(); j++)
0454     {
0455       if(j <= i) continue;
0456
0457       if(pi0PhoBarcode.at(i) != pi0PhoBarcode.at(j)) continue;
0458
0459       TLorentzVector tpi0, pho1, pho2;
0460       pho1.SetPxPyPzE(pi0Pho.at(i) -> get_px(), pi0Pho.at(i) -> get_py(), pi0Pho.at(i) -> get_pz(), pi0Pho.at(i) -> get_e());
0461       pho2.SetPxPyPzE(pi0Pho.at(j) -> get_px(), pi0Pho.at(j) -> get_py(), pi0Pho.at(j) -> get_pz(), pi0Pho.at(j) -> get_e());
0462
0463       tpi0 = pho1 + pho2;
0464       if(std::floor(100000*tpi0.M()) == std::floor(100000*pi0Mass))
0465         {
0466           massdiff = abs(tpi0.M() - pi0Mass);
0467           RawCluster *cluster1 = clustereval -> best_cluster_from(pi0Pho.at(i));
0468           RawCluster *cluster2 = clustereval -> best_cluster_from(pi0Pho.at(j));
0469           //.13497
0470           if(!cluster1 || !cluster2)
0471         {
0472           ePi0InvMassEcut[1] -> Fill(tpi0.Energy(), 0., std::min(pi0Pho.at(i) -> get_e(), pi0Pho.at(j) -> get_e()));
0473         }
0474           else
0475         {
0476           TLorentzVector cpi0, c1, c2;
0477
0478           CLHEP::Hep3Vector eVecCluster1 = RawClusterUtility::GetECoreVec(*cluster1, vertex);
0479           CLHEP::Hep3Vector eVecCluster2 = RawClusterUtility::GetECoreVec(*cluster2, vertex);
0480
0481           c1.SetPtEtaPhiE(eVecCluster1.perp(), eVecCluster1.pseudoRapidity(), eVecCluster1.getPhi(), eVecCluster1.mag());
0482           c2.SetPtEtaPhiE(eVecCluster2.perp(), eVecCluster2.pseudoRapidity(), eVecCluster2.getPhi(), eVecCluster2.mag());
0483
0484           cpi0 = c1 + c2;
0485
0486           ePi0InvMassEcut[1] -> Fill(cpi0.Energy(),cpi0.M(),std::min(c1.Energy(), c1.Energy()));
0487
0488           invMassEPhi -> Fill(cpi0.M(), tpi0.Energy(), cpi0.Phi());
0489           invMassPhoPi0 -> Fill(cpi0.M(), cpi0.Energy()/tpi0.Energy(), c1.Energy()/pho1.Energy());
0490           invMassPhoPi0 -> Fill(cpi0.M(), cpi0.Energy()/tpi0.Energy(), c1.Energy()/pho2.Energy());
0491
0492
0493         }
0494         }
0495     }
0496     }*/
0497
0498
0499
0500
0501   goodRecoCluster.clear();
0502   return Fun4AllReturnCodes::EVENT_OK;
0503 }
0504
0505 //____________________________________________________________________________..
0506 int pi0Efficiency::ResetEvent(PHCompositeNode *topNode)
0507 {
0508   return Fun4AllReturnCodes::EVENT_OK;
0509 }
0510
0511 //____________________________________________________________________________..
0512 int pi0Efficiency::EndRun(const int runnumber)
0513 {
0514   std::cout << "pi0Efficiency::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0515   return Fun4AllReturnCodes::EVENT_OK;
0516 }
0517
0518 //____________________________________________________________________________..
0519 int pi0Efficiency::End(PHCompositeNode *topNode)
0520 {
0521   std::cout << "pi0Efficiency::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0522
0523   out -> cd();
0524
0525
0526   for(int j = 0; j < nEtaBins; j++) ePi0InvMassEcut[j] -> Write();
0527
0528   clusterE -> Write();
0529   for(int i = 0; i < nEtaBins; i++)prim_pi0_E[i] -> Write();
0530   photonE -> Write();
0531   pi0EScale -> Write();
0532   ///unmatchedLocale -> Write();
0533   hMassRat -> Write();
0534   truthPi0EDeltaR -> Write();
0535   truthPi0EAsym -> Write();
0536   //unmatchedE -> Write();
0537
0538
0539
0540   return Fun4AllReturnCodes::EVENT_OK;
0541 }
0542
0543 //____________________________________________________________________________..
0544 int pi0Efficiency::Reset(PHCompositeNode *topNode)
0545 {
0546   std::cout << "pi0Efficiency::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0547   return Fun4AllReturnCodes::EVENT_OK;
0548 }
0549
0550 //____________________________________________________________________________..
0551 void pi0Efficiency::Print(const std::string &what) const
0552 {
0553   std::cout << "pi0Efficiency::Print(const std::string &what) const Printing info for " << what << std::endl;
0554 }
0555 //____________________________________________________________________________..
0556 float pi0Efficiency::getEta(PHG4Particle *particle)
0557 {
0558   float px = particle -> get_px();
0559   float py = particle -> get_py();
0560   float pz = particle -> get_pz();
0561   float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
0562
0563   return 0.5*log((p+pz)/(p-pz));
0564 }
0565 //____________________________________________________________________________..
0566 int pi0Efficiency::getEtaBin(float eta)
0567 {
0568   for(int i = 0; i < nEtaBins-1; i++)
0569     {
0570       if(abs(eta) < (i+1)/(float)(nEtaBins-1) * 1.1 && abs(eta) >= (i+1-1)/(float)(nEtaBins-1) * 1.1) return i;
0571     }
0572  return -1;
0573 }