pi0ClusterAna/src/pi0ClusterAna.cc

0001 #include "pi0ClusterAna.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 #include <TTree.h>
0019
0020 //for emc clusters
0021 #include <calobase/RawCluster.h>
0022 #include <calobase/RawClusterContainer.h>
0023 #include <calobase/RawClusterUtility.h>
0024 #include <calobase/RawTowerGeomContainer.h>
0025 #include <calobase/RawTower.h>
0026 #include <calobase/RawTowerContainer.h>
0027
0028 //for vetex information
0029 #include <g4vertex/GlobalVertex.h>
0030 #include <g4vertex/GlobalVertexMap.h>
0031
0032 //tracking
0033 #include <trackbase_historic/SvtxTrack.h>
0034 #include <trackbase_historic/SvtxTrackMap.h>
0035 #include <trackbase_historic/SvtxVertex.h>
0036 #include <trackbase_historic/SvtxVertexMap.h>
0037
0038 //truth information
0039 #include <g4main/PHG4TruthInfoContainer.h>
0040 #include <g4main/PHG4Particle.h>
0041 #include <g4main/PHG4VtxPoint.h>
0042 #pragma GCC diagnostic push
0043 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0044 #include <HepMC/GenEvent.h>
0045 #include <HepMC/GenParticle.h>
0046 #include <HepMC/GenVertex.h>
0047 #include <HepMC/IteratorRange.h>
0048 #include <HepMC/SimpleVector.h>
0049 #include <HepMC/GenParticle.h>
0050 #pragma GCC diagnostic pop
0051 #include <CLHEP/Geometry/Point3D.h>
0052
0053
0054 //____________________________________________________________________________..
0055 pi0ClusterAna::pi0ClusterAna(const std::string &name, const std::string &outName = "pi0ClusterAnaOut"):
0056 SubsysReco(name)
0057   , clusters_Towers(nullptr)
0058   , truth_photon(nullptr)
0059   , truth_pi0(nullptr)
0060   //, caloevalstack(NULL)
0061   , m_eta_center()
0062   , m_phi_center()
0063   , m_tower_energy()
0064 , m_cluster_eta()
0065   , m_cluster_phi()
0066 , m_cluster_e()
0067   , m_cluster_chi2()
0068   , m_cluster_prob()
0069 , m_cluster_nTowers()
0070   , m_asym()
0071   , m_deltaR()
0072   , m_lead_E()
0073 , m_sublead_E()
0074 , m_lead_phi()
0075 , m_lead_eta()
0076   , m_sublead_phi()
0077   , m_sublead_eta()
0078   , m_pi0_E()
0079 , m_pi0_eta()
0080 , m_pi0_phi()
0081   , Outfile(outName)
0082
0083 {
0084
0085
0086   std::cout << "pi0ClusterAna::pi0ClusterAna(const std::string &name) Calling ctor" << std::endl;
0087 }
0088
0089 //____________________________________________________________________________..
0090 pi0ClusterAna::~pi0ClusterAna()
0091 {
0092   std::cout << "pi0ClusterAna::~pi0ClusterAna() Calling dtor" << std::endl;
0093 }
0094
0095 //____________________________________________________________________________..
0096 int pi0ClusterAna::Init(PHCompositeNode *topNode)
0097 {
0098
0099   clusters_Towers = new TTree("TreeClusterTower","Tree for cluster and tower information");
0100   clusters_Towers -> Branch("towerEtaCEnter",&m_eta_center);
0101   clusters_Towers -> Branch("towerPhiCenter",& m_phi_center);
0102   clusters_Towers -> Branch("towerEnergy",&m_tower_energy);
0103   clusters_Towers -> Branch("clusterEta",&m_cluster_eta);
0104   clusters_Towers -> Branch("clusterPhi", &m_cluster_phi);
0105   clusters_Towers -> Branch("clusterE",&m_cluster_e);
0106   clusters_Towers -> Branch("clusterChi2",&m_cluster_chi2);
0107   clusters_Towers -> Branch("clusterProb",&m_cluster_prob);
0108   clusters_Towers -> Branch("clusterNTowers",&m_cluster_nTowers);
0109
0110   truth_photon = new TTree("TreeTruthPhoton", "Tree for truth photons");
0111   truth_photon -> Branch("pairAsym",&m_asym);
0112   truth_photon -> Branch("pairDeltaR",&m_deltaR);
0113   truth_photon -> Branch("leadPhotonPhi",&m_lead_phi);
0114   truth_photon -> Branch("leadPhotonEta",&m_lead_eta);
0115   truth_photon -> Branch("subleadPhotonPhi",&m_sublead_phi);
0116   truth_photon -> Branch("subleadPhotonEta",&m_sublead_eta);
0117   truth_photon -> Branch("leadPhotonE", &m_lead_E);
0118   truth_photon -> Branch("subLeadPhotonE", &m_sublead_E);
0119
0120   truth_pi0 = new TTree("TreeTruthPi0", "Tree for Truth pi0");
0121   truth_pi0 -> Branch("pi0E",&m_pi0_E);
0122   truth_pi0 -> Branch("pi0_phi",&m_pi0_phi);
0123   truth_pi0 -> Branch("pi0_eta",&m_pi0_eta);
0124
0125   //so that the histos actually get written out
0126   Fun4AllServer *se = Fun4AllServer::instance();
0127   se -> Print("NODETREE");
0128   hm = new Fun4AllHistoManager("MYHISTOS");
0129
0130   se -> registerHistoManager(hm);
0131
0132   se -> registerHisto(truth_pi0 -> GetName(),truth_pi0);
0133   se -> registerHisto(truth_photon -> GetName(), truth_photon);
0134   se -> registerHisto(clusters_Towers -> GetName(), clusters_Towers);
0135
0136   out = new TFile(Outfile.c_str(),"RECREATE");
0137
0138   std::cout << "pi0ClusterAna::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 //____________________________________________________________________________..
0143 int pi0ClusterAna::InitRun(PHCompositeNode *topNode)
0144 {
0145   std::cout << "pi0ClusterAna::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0146
0147   return Fun4AllReturnCodes::EVENT_OK;
0148 }
0149
0150 //____________________________________________________________________________..
0151 int pi0ClusterAna::process_event(PHCompositeNode *topNode)
0152 {
0153
0154   //Information on clusters
0155   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0156   //RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0157   if(!clusterContainer)
0158     {
0159       std::cout << PHWHERE << "pi0Efficiency::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0160
0161       return 0;
0162     }
0163
0164
0165   //Vertex information
0166   GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
0167   if (!vtxContainer)
0168     {
0169       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;
0170       assert(vtxContainer);  // force quit
0171
0172       return 0;
0173     }
0174
0175   if (vtxContainer->empty())
0176     {
0177       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;
0178       return 0;
0179     }
0180
0181   //More vertex information
0182   GlobalVertex *vtx = vtxContainer->begin()->second;
0183   if(!vtx)
0184     {
0185
0186       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find vtx from vtxContainer"  << std::endl;
0187       return Fun4AllReturnCodes::ABORTEVENT;
0188     }
0189
0190   //Tower geometry node for location information
0191   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0192   if (!towergeom)
0193     {
0194       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node TOWERGEOM_CEMC"  << std::endl;
0195       return Fun4AllReturnCodes::ABORTEVENT;
0196     }
0197
0198   //Raw tower information
0199   RawTowerContainer *towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC");
0200   if(!towerContainer)
0201     {
0202       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node TOWER_CALIB_CEMC"  << std::endl;
0203       return Fun4AllReturnCodes::ABORTEVENT;
0204     }
0205
0206   //truth particle information
0207   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0208   if(!truthinfo)
0209     {
0210       std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node G4TruthInfo"  << std::endl;
0211       return Fun4AllReturnCodes::ABORTEVENT;
0212     }
0213
0214
0215   float pi0Eta = -99999;
0216   float photonEtaMax = 1.1;
0217   float mesonEtaMax = 0.3;
0218   //Now we go and find our truth pi0s and just take its eta coordinate for now
0219   PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange();
0220   PHG4TruthInfoContainer::ConstIterator truthIter;
0221
0222   PHG4Particle *truthPar = NULL;
0223
0224   for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
0225     {
0226       truthPar = truthIter->second;
0227
0228       if(truthPar -> get_pid() == 111 && truthPar -> get_parent_id() == 0)
0229     {
0230       pi0Eta = getEta(truthPar);
0231       if(abs(pi0Eta) >= mesonEtaMax)
0232         {
0233           return 0;
0234         }
0235     }
0236     }
0237
0238   int firstphotonflag = 0;
0239   int firstfirstphotonflag = 0;
0240   int secondphotonflag = 0;
0241
0242   //int secondsecondphotonflag = 0;
0243
0244   PHG4TruthInfoContainer::Range truthRangeDecay1 = truthinfo->GetSecondaryParticleRange();
0245   PHG4TruthInfoContainer::ConstIterator truthIterDecay1;
0246
0247
0248   TLorentzVector photon1;
0249   TLorentzVector photon2;
0250   int nParticles = 0;
0251
0252   //loop over all our decay photons.
0253   //Make sure they fall within the desired acceptance
0254   //Toss Dalitz decays
0255   //Retain photon kinematics to determine lead photon for eta binning
0256   for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
0257     {
0258       PHG4Particle *decay = truthIterDecay1 -> second;
0259
0260       int dumtruthpid = decay -> get_pid();
0261       int dumparentid = decay -> get_parent_id();
0262
0263       //if the parent is the pi0 and it's a photon and we haven't marked one yet
0264       if(dumparentid == 1 && dumtruthpid == 22 && !firstphotonflag)
0265     {
0266       if(abs(getEta(decay)) > photonEtaMax)
0267         {
0268           return 0;
0269         }
0270       photon1.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e());
0271
0272       firstphotonflag = 1;
0273     }
0274
0275       if(dumparentid == 1 && dumtruthpid == 22 && firstphotonflag && firstfirstphotonflag)
0276     {
0277       if(abs(getEta(decay)) > photonEtaMax)
0278         {
0279           return 0;
0280         }
0281       photon2.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e()) ;
0282
0283       secondphotonflag = 1;
0284     }
0285
0286       //Need this flag to make it skip the first photon slot
0287       if(firstphotonflag) firstfirstphotonflag = 1;
0288       if(dumparentid == 1)nParticles ++;
0289     }
0290
0291   if((!firstphotonflag || !secondphotonflag) && nParticles > 1) //Dalitz
0292     {
0293       return 0;
0294     }
0295   else if((!firstphotonflag || !secondphotonflag) && nParticles == 1) //One photon falls outside simulation acceptance
0296     {
0297       return 0;
0298     }
0299   else if((!firstphotonflag || !secondphotonflag) && nParticles == 0) //Both photons fall outside simulation acceptance
0300     {
0301       return 0;
0302     }
0303
0304
0305   TLorentzVector pi0;
0306   pi0.SetPxPyPzE(truthPar -> get_px(), truthPar -> get_py(), truthPar -> get_pz(), truthPar -> get_e());
0307   TLorentzVector leadPho, subLeadPho;
0308   if(photon1.Energy() >= photon2.Energy())
0309     {
0310       leadPho = photon1;
0311       subLeadPho = photon2;
0312     }
0313   else
0314     {
0315       leadPho = photon2;
0316       subLeadPho = photon1;
0317     }
0318
0319   float asym = abs(photon1.Energy() - photon2.Energy())/(photon1.Energy() + photon2.Energy());
0320
0321   m_asym.push_back(asym);
0322
0323   float deltaR = photon1.DeltaR(photon2);
0324
0325   m_deltaR.push_back(deltaR);
0326
0327   m_lead_phi.push_back(leadPho.Phi());
0328   m_sublead_phi.push_back(subLeadPho.Phi());
0329
0330   m_lead_eta.push_back(leadPho.PseudoRapidity());
0331   m_sublead_eta.push_back(subLeadPho.PseudoRapidity());
0332
0333   m_lead_E.push_back(leadPho.Energy());
0334   m_sublead_E.push_back(subLeadPho.Energy());
0335
0336   m_pi0_E.push_back(pi0.Energy());
0337   m_pi0_phi.push_back(pi0.Phi());
0338   m_pi0_eta.push_back(pi0.PseudoRapidity());
0339
0340   //grab all the towers and fill their energies.
0341   RawTowerContainer::ConstRange tower_range = towerContainer -> getTowers();
0342   for(RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter!= tower_range.second; tower_iter++)
0343     {
0344       int phibin = tower_iter -> second -> get_binphi();
0345       int etabin = tower_iter -> second -> get_bineta();
0346       double phi = towergeom -> get_phicenter(phibin);
0347       double eta = towergeom -> get_etacenter(etabin);
0348       double energy = tower_iter -> second -> get_energy();
0349
0350       m_phi_center.push_back(phi);
0351       m_eta_center.push_back(eta);
0352       m_tower_energy.push_back(energy);
0353     }
0354
0355   //This is how we iterate over items in the container.
0356   RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
0357   RawClusterContainer::ConstIterator clusterIter;
0358
0359   for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0360     {
0361       RawCluster *recoCluster = clusterIter -> second;
0362
0363       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0364       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0365       float clusE = E_vec_cluster.mag();
0366       float clus_eta = E_vec_cluster.pseudoRapidity();
0367       float clus_phi = E_vec_cluster.phi();
0368
0369       m_cluster_eta.push_back(clus_eta);
0370       m_cluster_phi.push_back(clus_phi);
0371       m_cluster_e.push_back(clusE);
0372
0373       m_cluster_chi2.push_back(recoCluster -> get_chi2());
0374       m_cluster_prob.push_back(recoCluster -> get_prob());
0375       m_cluster_nTowers.push_back(recoCluster -> getNTowers());
0376     }
0377
0378
0379   clusters_Towers -> Fill();
0380   truth_photon -> Fill();
0381   truth_pi0 -> Fill();
0382
0383
0384   return Fun4AllReturnCodes::EVENT_OK;
0385 }
0386
0387 //____________________________________________________________________________..
0388 int pi0ClusterAna::ResetEvent(PHCompositeNode *topNode)
0389 {
0390   //std::cout << "pi0ClusterAna::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0391
0392   m_eta_center.clear();
0393   m_phi_center.clear();
0394   m_tower_energy.clear();
0395   m_cluster_eta.clear();
0396   m_cluster_phi.clear();
0397   m_cluster_e.clear();
0398   m_cluster_chi2.clear();
0399   m_cluster_prob.clear();
0400   m_cluster_nTowers.clear();
0401   m_asym.clear();
0402   m_deltaR.clear();
0403   m_lead_E.clear();
0404   m_sublead_E.clear();
0405   m_lead_phi.clear();
0406   m_lead_eta.clear();
0407   m_sublead_phi.clear();
0408   m_sublead_eta.clear();
0409   m_pi0_E.clear();
0410   m_pi0_eta.clear();
0411   m_pi0_phi.clear();
0412   return Fun4AllReturnCodes::EVENT_OK;
0413 }
0414
0415 //____________________________________________________________________________..
0416 int pi0ClusterAna::EndRun(const int runnumber)
0417 {
0418   std::cout << "pi0ClusterAna::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0419   return Fun4AllReturnCodes::EVENT_OK;
0420 }
0421
0422 //____________________________________________________________________________..
0423 int pi0ClusterAna::End(PHCompositeNode *topNode)
0424 {
0425   std::cout << "pi0ClusterAna::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0426   out -> cd();
0427
0428   truth_pi0 -> Write();
0429   truth_photon -> Write();
0430   clusters_Towers -> Write();
0431
0432   out -> Close();
0433   delete out;
0434
0435   hm -> dumpHistos(Outfile.c_str(),"UPDATE");
0436
0437
0438   return Fun4AllReturnCodes::EVENT_OK;
0439 }
0440
0441 //____________________________________________________________________________..
0442 int pi0ClusterAna::Reset(PHCompositeNode *topNode)
0443 {
0444   std::cout << "pi0ClusterAna::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0445   return Fun4AllReturnCodes::EVENT_OK;
0446 }
0447
0448 //____________________________________________________________________________..
0449 void pi0ClusterAna::Print(const std::string &what) const
0450 {
0451   std::cout << "pi0ClusterAna::Print(const std::string &what) const Printing info for " << what << std::endl;
0452 }
0453 //____________________________________________________________________________..
0454 float pi0ClusterAna::getEta(PHG4Particle *particle)
0455 {
0456   float px = particle -> get_px();
0457   float py = particle -> get_py();
0458   float pz = particle -> get_pz();
0459   float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
0460
0461   return 0.5*log((p+pz)/(p-pz));
0462 }