sPhenix code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-08-05 08:14:51
0001 0002 //____________________________________________________________________________.. 0003 0004 //our base code 0005 #include "pythiaEMCalAna.h" 0006 0007 #include <fun4all/Fun4AllReturnCodes.h> 0008 0009 #include <phool/PHCompositeNode.h> 0010 0011 //Fun4all stuff 0012 #include <fun4all/Fun4AllReturnCodes.h> 0013 #include <fun4all/Fun4AllServer.h> 0014 #include <fun4all/Fun4AllHistoManager.h> 0015 #include <phool/PHCompositeNode.h> 0016 #include <phool/getClass.h> 0017 #include <phool/phool.h> 0018 #include <ffaobjects/EventHeader.h> 0019 0020 //ROOT stuff 0021 #include <TH1F.h> 0022 #include <TH2F.h> 0023 #include <TH3F.h> 0024 #include <TFile.h> 0025 #include <TLorentzVector.h> 0026 #include <TTree.h> 0027 0028 //for emc clusters 0029 #include <calobase/RawCluster.h> 0030 #include <calobase/RawClusterContainer.h> 0031 #include <calobase/RawClusterUtility.h> 0032 #include <calobase/RawTowerGeomContainer.h> 0033 #include <calobase/RawTower.h> 0034 #include <calobase/RawTowerContainer.h> 0035 #include <calobase/TowerInfo.h> 0036 #include <calobase/TowerInfoContainer.h> 0037 0038 // Minimum Bias 0039 #include <calotrigger/MinimumBiasInfo.h> 0040 0041 //caloEvalStack for cluster to truth matching 0042 #include <g4eval/CaloEvalStack.h> 0043 /* #include <g4eval/CaloRawClusterEval.h> */ 0044 0045 //for vertex information 0046 #include <globalvertex/GlobalVertex.h> 0047 #include <globalvertex/GlobalVertexMap.h> 0048 /* #include <g4vertex/GlobalVertex.h> */ 0049 /* #include <g4vertex/GlobalVertexMap.h> */ 0050 0051 //tracking 0052 #include <trackbase_historic/SvtxTrack.h> 0053 #include <trackbase_historic/SvtxTrackMap.h> 0054 /* #include <trackbase_historic/SvtxVertex.h> */ 0055 /* #include <trackbase_historic/SvtxVertexMap.h> */ 0056 0057 //truth information 0058 #include <g4main/PHG4TruthInfoContainer.h> 0059 #include <g4main/PHG4Particle.h> 0060 #include <g4main/PHG4Shower.h> 0061 #include <g4main/PHG4VtxPoint.h> 0062 #include <phhepmc/PHHepMCGenEvent.h> 0063 #include <phhepmc/PHHepMCGenEventMap.h> 0064 #pragma GCC diagnostic push 0065 #pragma GCC diagnostic ignored "-Wdeprecated-declarations" 0066 #include <HepMC/GenEvent.h> 0067 #include <HepMC/GenParticle.h> 0068 #include <HepMC/GenVertex.h> 0069 #include <HepMC/IteratorRange.h> 0070 #include <HepMC/SimpleVector.h> 0071 #include <HepMC/GenParticle.h> 0072 #pragma GCC diagnostic pop 0073 0074 // PID stuff to print info on how many of each truth particle there were 0075 // Mostly useful as a sanity check 0076 0077 bool IsBaryon(int pdg_id) { 0078 /* std::cout << Form("Entering IsBaryon(%d)\n", pdg_id); */ 0079 bool ret = false; 0080 // baryon ids are 4-digit numbers 0081 if (abs(pdg_id) > 999 && abs(pdg_id) < 10000) ret = true; 0082 /* std::cout << Form("In IsBaryon(%d): returning %s\n", pdg_id, ret ? "true": "false"); */ 0083 return ret; 0084 } 0085 0086 bool IsNeutralBaryon(int pdg_id) { 0087 /* std::cout << Form("Entering IsNeutralBaryon(%d)\n", pdg_id); */ 0088 if (!IsBaryon(pdg_id)) return false; 0089 // baryon ids are 4-digit numbers 0090 int base_id = abs(pdg_id) % 10000; 0091 // the last digit is the angular momentum state 0092 // the first three digits are the quark flavors 0093 int q1 = (int)(base_id/1000); 0094 int q2 = (int)(base_id/100 % 10); 0095 int q3 = (int)(base_id/10 % 10); 0096 // odd numbers 1,3,5 are d,s,b; even numbers 2,4,6 are u,c,t 0097 float charge = 0.; 0098 if (q1%2 == 0) charge += 2./3.; 0099 else charge -= 1./3.; 0100 if (q2%2 == 0) charge += 2./3.; 0101 else charge -= 1./3.; 0102 if (q3%2 == 0) charge += 2./3.; 0103 else charge -= 1./3.; 0104 /* std::cout << Form("In IsNeutralBaryon: base_id=%d, charge=%f\n", base_id, charge); */ 0105 if (abs(charge) < 0.05) return true; 0106 else return false; 0107 } 0108 0109 bool IsNeutralMeson(int pdg_id) { 0110 /* std::cout << Form("Entering IsNeutralMeson(%d)\n", pdg_id); */ 0111 // excited states are >3 digit numbers; first reduce them to the base quark content 0112 int base_id = pdg_id % 1000; 0113 // meson ids are 3-digit numbers 0114 // the last digit is the angular momentum state 0115 // the first two digits are the quark flavors 0116 int q1 = (int)(base_id/100); 0117 int q2 = (int)(base_id/10 % 10); 0118 float charge = 0.; 0119 if (q1%2 == 0) charge += 2./3.; 0120 else charge -= 1./3.; 0121 if (q2%2 == 0) charge -= 2./3.; 0122 else charge += 1./3.; 0123 /* std::cout << Form("In IsNeutralMeson: base_id=%d, charge=%f\n", base_id, charge); */ 0124 if (abs(charge) < 0.05) return true; 0125 else return false; 0126 } 0127 0128 bool IsNeutralHadron(int pdg_id) { 0129 /* std::cout << Form("\nEntering IsNeutralHadron(%d)\n", pdg_id); */ 0130 bool ret = (IsNeutralBaryon(pdg_id) || IsNeutralMeson(pdg_id)); 0131 /* std::cout << Form("In IsNeutralHadron: returning %s\n", ret ? "true" : "false"); */ 0132 return ret; 0133 } 0134 0135 bool IsChargedHadron(int pdg_id) { 0136 if (pdg_id < 0) return true; 0137 if (pdg_id < 99) return false; 0138 return (!IsNeutralHadron(pdg_id)); 0139 } 0140 0141 //____________________________________________________________________________.. 0142 pythiaEMCalAna::pythiaEMCalAna(const std::string &name, const std::string &oname, bool isMC, bool hasPythia): 0143 SubsysReco(name), 0144 clusters_Towers(nullptr), 0145 truth_particles(nullptr), 0146 m_truthInfo(nullptr), 0147 m_theEvent(nullptr), 0148 m_clusterContainer(nullptr), 0149 m_clusterEval(nullptr), 0150 m_tower_energy(0), 0151 m_eta_center(0), 0152 m_phi_center(0), 0153 m_cluster_ID(0), 0154 m_cluster_e(0), 0155 m_cluster_eta(0), 0156 m_cluster_phi(0), 0157 m_cluster_ecore(0), 0158 m_cluster_chi2(0), 0159 m_cluster_prob(0), 0160 m_cluster_nTowers(0), 0161 m_cluster_allTowersE(0), 0162 m_cluster_allTowersEta(0), 0163 m_cluster_allTowersPhi(0), 0164 m_cluster_nParticles(0), 0165 m_cluster_primaryParticlePid(0), 0166 m_cluster_allSecondaryPids(0), 0167 /* m_cluster_maxE_trackID(0), */ 0168 /* m_cluster_maxE_PID(0), */ 0169 /* m_cluster_all_trackIDs(0), */ 0170 m_truth_Parent_Barcode(0), 0171 m_truth_Parent_Pid(0), 0172 m_truth_Parent_M(0), 0173 m_truth_Parent_E(0), 0174 m_truth_Parent_Eta(0), 0175 m_truth_Parent_Phi(0), 0176 m_truth_Parent_Pt(0), 0177 m_truth_Parent_xF(0), 0178 m_truth_Parent_EndVtx_x(0), 0179 m_truth_Parent_EndVtx_y(0), 0180 m_truth_Parent_EndVtx_z(0), 0181 /* m_truth_Parent_EndVtx_t(0), */ 0182 m_truth_Parent_EndVtx_r(0), 0183 /* m_truth_Parent_TotalNDaughters(0), */ 0184 /* m_truth_Parent_AllDaughterPids(0), */ 0185 /* m_truth_Parent_AllDaughterEnergyFractions(0), */ 0186 m_truth_Parent_TotalNClusters(0), 0187 m_truth_Parent_AllClusterIDs(0), 0188 m_truth_Parent_AllClusterEnergyFractions(0), 0189 m_truth_Decay1_Barcode(0), 0190 m_truth_Decay1_Pid(0), 0191 m_truth_Decay1_M(0), 0192 m_truth_Decay1_E(0), 0193 m_truth_Decay1_Eta(0), 0194 m_truth_Decay1_Phi(0), 0195 m_truth_Decay1_Pt(0), 0196 m_truth_Decay1_xF(0), 0197 m_truth_Decay1_EndVtx_x(0), 0198 m_truth_Decay1_EndVtx_y(0), 0199 m_truth_Decay1_EndVtx_z(0), 0200 /* m_truth_Decay1_EndVtx_t(0), */ 0201 m_truth_Decay1_EndVtx_r(0), 0202 m_truth_Decay1_TotalNClusters(0), 0203 m_truth_Decay1_BestClusterID(0), 0204 m_truth_Decay1_BestClusterEfraction(0), 0205 /* m_truth_Decay1_AllClusterIDs(0), */ 0206 /* m_truth_Decay1_AllClusterEnergyFractions(0), */ 0207 m_truth_Decay2_Barcode(0), 0208 m_truth_Decay2_Pid(0), 0209 m_truth_Decay2_M(0), 0210 m_truth_Decay2_E(0), 0211 m_truth_Decay2_Eta(0), 0212 m_truth_Decay2_Phi(0), 0213 m_truth_Decay2_Pt(0), 0214 m_truth_Decay2_xF(0), 0215 m_truth_Decay2_EndVtx_x(0), 0216 m_truth_Decay2_EndVtx_y(0), 0217 m_truth_Decay2_EndVtx_z(0), 0218 /* m_truth_Decay2_EndVtx_t(0), */ 0219 m_truth_Decay2_EndVtx_r(0), 0220 m_truth_Decay2_TotalNClusters(0), 0221 m_truth_Decay2_BestClusterID(0), 0222 m_truth_Decay2_BestClusterEfraction(0), 0223 /* m_truth_Decay2_AllClusterIDs(0), */ 0224 /* m_truth_Decay2_AllClusterEnergyFractions(0), */ 0225 /* m_truth_Cluster1_Id(0), */ 0226 /* m_truth_Cluster1_E(0), */ 0227 /* m_truth_Cluster1_Ecore(0), */ 0228 /* m_truth_Cluster1_Eta(0), */ 0229 /* m_truth_Cluster1_Phi(0), */ 0230 /* m_truth_Cluster1_Pt(0), */ 0231 /* m_truth_Cluster1_xF(0), */ 0232 /* m_truth_Cluster1_Chi2(0), */ 0233 /* m_truth_Cluster1_Decay1EnergyFraction(0), */ 0234 /* m_truth_Cluster1_Decay2EnergyFraction(0), */ 0235 /* m_truth_Cluster2_Id(0), */ 0236 /* m_truth_Cluster2_E(0), */ 0237 /* m_truth_Cluster2_Ecore(0), */ 0238 /* m_truth_Cluster2_Eta(0), */ 0239 /* m_truth_Cluster2_Phi(0), */ 0240 /* m_truth_Cluster2_Pt(0), */ 0241 /* m_truth_Cluster2_xF(0), */ 0242 /* m_truth_Cluster2_Chi2(0), */ 0243 /* m_truth_Cluster2_Decay1EnergyFraction(0), */ 0244 /* m_truth_Cluster2_Decay2EnergyFraction(0), */ 0245 fout(NULL), 0246 outname(oname), 0247 getEvent(-9999), 0248 /* hasHIJING(isAuAu), */ 0249 isMonteCarlo(isMC), 0250 hasPYTHIA(hasPythia), 0251 n_events_total(0), 0252 n_events_minbias(0), 0253 n_events_with_vertex(0), 0254 n_events_with_good_vertex(0), 0255 n_events_positiveCaloE(0), 0256 n_primaries(0), 0257 n_primary_photons(0), 0258 n_direct_photons(0), 0259 n_leptons(0), 0260 n_neutral_hadrons(0), 0261 n_neutral_hadrons_geant(0), 0262 n_neutral_hadrons_pythia(0), 0263 n_charged_hadrons(0), 0264 n_charged_hadrons_geant(0), 0265 n_charged_hadrons_pythia(0), 0266 n_pythia_decayed_hadrons(0), 0267 allTrackIDs(), 0268 primaryBarcodes() 0269 { 0270 std::cout << "pythiaEMCalAna::pythiaEMCalAna(const std::string &name) Calling ctor" << std::endl; 0271 } 0272 //____________________________________________________________________________.. 0273 pythiaEMCalAna::~pythiaEMCalAna() 0274 { 0275 std::cout << "pythiaEMCalAna::~pythiaEMCalAna() Calling dtor" << std::endl; 0276 } 0277 0278 //____________________________________________________________________________.. 0279 int pythiaEMCalAna::Init(PHCompositeNode *topNode) 0280 { 0281 std::cout << "pythiaEMCalAna::Init(PHCompositeNode *topNode) Initializing" << std::endl; 0282 0283 fout = new TFile(outname.c_str(),"RECREATE"); 0284 0285 clusters_Towers = new TTree("TreeClusterTower","Tree for cluster and tower information"); 0286 clusters_Towers -> Branch("towerEnergy",&m_tower_energy); 0287 clusters_Towers -> Branch("towerEtaCenter",&m_eta_center); 0288 clusters_Towers -> Branch("towerPhiCenter",& m_phi_center); 0289 clusters_Towers -> Branch("clusterID",&m_cluster_ID); 0290 clusters_Towers -> Branch("clusterE",&m_cluster_e); 0291 clusters_Towers -> Branch("clusterEta",&m_cluster_eta); 0292 clusters_Towers -> Branch("clusterPhi", &m_cluster_phi); 0293 clusters_Towers -> Branch("clusterEcore",&m_cluster_ecore); 0294 clusters_Towers -> Branch("clusterChi2",&m_cluster_chi2); 0295 clusters_Towers -> Branch("clusterProb",&m_cluster_prob); 0296 clusters_Towers -> Branch("clusterNTowers",&m_cluster_nTowers); 0297 clusters_Towers -> Branch("clusterAllTowersE",&m_cluster_allTowersE); 0298 clusters_Towers -> Branch("clusterAllTowersEta",&m_cluster_allTowersEta); 0299 clusters_Towers -> Branch("clusterAllTowersPhi",&m_cluster_allTowersPhi); 0300 if (isMonteCarlo) { 0301 clusters_Towers -> Branch("clusterNParticles",&m_cluster_nParticles); 0302 clusters_Towers -> Branch("clusterPrimaryParticlePid",&m_cluster_primaryParticlePid); 0303 clusters_Towers -> Branch("clusterAllParticlePids",&m_cluster_allSecondaryPids, 32000, 0); 0304 /* clusters_Towers -> Branch("clusterMaxE_trackID",&m_cluster_maxE_trackID); */ 0305 /* clusters_Towers -> Branch("clusterMaxE_PID",&m_cluster_maxE_PID); */ 0306 /* clusters_Towers -> Branch("clusterAll_trackIDs",&m_cluster_all_trackIDs, 32000, 0); */ 0307 } 0308 0309 if (isMonteCarlo) { 0310 truth_particles = new TTree("TreeTruthParticles", "Tree for Truth Particles"); 0311 truth_particles->Branch("truth_Parent_Barcode", &m_truth_Parent_Barcode); 0312 truth_particles->Branch("truth_Parent_Pid", &m_truth_Parent_Pid); 0313 truth_particles->Branch("truth_Parent_M", &m_truth_Parent_M); 0314 truth_particles->Branch("truth_Parent_E", &m_truth_Parent_E); 0315 truth_particles->Branch("truth_Parent_Eta", &m_truth_Parent_Eta); 0316 truth_particles->Branch("truth_Parent_Phi", &m_truth_Parent_Phi); 0317 truth_particles->Branch("truth_Parent_Pt", &m_truth_Parent_Pt); 0318 truth_particles->Branch("truth_Parent_xF", &m_truth_Parent_xF); 0319 truth_particles->Branch("truth_Parent_EndVtx_x", &m_truth_Parent_EndVtx_x); 0320 truth_particles->Branch("truth_Parent_EndVtx_y", &m_truth_Parent_EndVtx_y); 0321 truth_particles->Branch("truth_Parent_EndVtx_z", &m_truth_Parent_EndVtx_z); 0322 /* truth_particles->Branch("truth_Parent_EndVtx_t", &m_truth_Parent_EndVtx_t); */ 0323 truth_particles->Branch("truth_Parent_EndVtx_r", &m_truth_Parent_EndVtx_r); 0324 truth_particles->Branch("truth_Parent_TotalNClusters", &m_truth_Parent_TotalNClusters); 0325 truth_particles->Branch("truth_Parent_AllClusterIDs", &m_truth_Parent_AllClusterIDs); 0326 truth_particles->Branch("truth_Parent_AllClusterEnergyFractions", &m_truth_Parent_AllClusterEnergyFractions); 0327 truth_particles->Branch("truth_Decay1_Barcode", &m_truth_Decay1_Barcode); 0328 truth_particles->Branch("truth_Decay1_Pid", &m_truth_Decay1_Pid); 0329 truth_particles->Branch("truth_Decay1_M", &m_truth_Decay1_M); 0330 truth_particles->Branch("truth_Decay1_E", &m_truth_Decay1_E); 0331 truth_particles->Branch("truth_Decay1_Eta", &m_truth_Decay1_Eta); 0332 truth_particles->Branch("truth_Decay1_Phi", &m_truth_Decay1_Phi); 0333 truth_particles->Branch("truth_Decay1_Pt", &m_truth_Decay1_Pt); 0334 truth_particles->Branch("truth_Decay1_xF", &m_truth_Decay1_xF); 0335 truth_particles->Branch("truth_Decay1_EndVtx_x", &m_truth_Decay1_EndVtx_x); 0336 truth_particles->Branch("truth_Decay1_EndVtx_y", &m_truth_Decay1_EndVtx_y); 0337 truth_particles->Branch("truth_Decay1_EndVtx_z", &m_truth_Decay1_EndVtx_z); 0338 /* truth_particles->Branch("truth_Decay1_EndVtx_t", &m_truth_Decay1_EndVtx_t); */ 0339 truth_particles->Branch("truth_Decay1_EndVtx_r", &m_truth_Decay1_EndVtx_r); 0340 truth_particles->Branch("truth_Decay1_TotalNClusters", &m_truth_Decay1_TotalNClusters); 0341 truth_particles->Branch("truth_Decay1_BestClusterID", &m_truth_Decay1_BestClusterID); 0342 truth_particles->Branch("truth_Decay1_BestClusterEfraction", &m_truth_Decay1_BestClusterEfraction); 0343 /* truth_particles->Branch("truth_Decay1_AllClusterIDs", &m_truth_Decay1_AllClusterIDs); */ 0344 /* truth_particles->Branch("truth_Decay1_AllClusterEnergyFractions", &m_truth_Decay1_AllClusterEnergyFractions); */ 0345 truth_particles->Branch("truth_Decay2_Barcode", &m_truth_Decay2_Barcode); 0346 truth_particles->Branch("truth_Decay2_Pid", &m_truth_Decay2_Pid); 0347 truth_particles->Branch("truth_Decay2_M", &m_truth_Decay2_M); 0348 truth_particles->Branch("truth_Decay2_E", &m_truth_Decay2_E); 0349 truth_particles->Branch("truth_Decay2_Eta", &m_truth_Decay2_Eta); 0350 truth_particles->Branch("truth_Decay2_Phi", &m_truth_Decay2_Phi); 0351 truth_particles->Branch("truth_Decay2_Pt", &m_truth_Decay2_Pt); 0352 truth_particles->Branch("truth_Decay2_xF", &m_truth_Decay2_xF); 0353 truth_particles->Branch("truth_Decay2_EndVtx_x", &m_truth_Decay2_EndVtx_x); 0354 truth_particles->Branch("truth_Decay2_EndVtx_y", &m_truth_Decay2_EndVtx_y); 0355 truth_particles->Branch("truth_Decay2_EndVtx_z", &m_truth_Decay2_EndVtx_z); 0356 /* truth_particles->Branch("truth_Decay2_EndVtx_t", &m_truth_Decay2_EndVtx_t); */ 0357 truth_particles->Branch("truth_Decay2_EndVtx_r", &m_truth_Decay2_EndVtx_r); 0358 truth_particles->Branch("truth_Decay2_TotalNClusters", &m_truth_Decay2_TotalNClusters); 0359 truth_particles->Branch("truth_Decay2_BestClusterID", &m_truth_Decay2_BestClusterID); 0360 truth_particles->Branch("truth_Decay2_BestClusterEfraction", &m_truth_Decay2_BestClusterEfraction); 0361 /* truth_particles->Branch("truth_Decay2_AllClusterIDs", &m_truth_Decay2_AllClusterIDs); */ 0362 /* truth_particles->Branch("truth_Decay2_AllClusterEnergyFractions", &m_truth_Decay2_AllClusterEnergyFractions); */ 0363 /* truth_particles->Branch("truth_Cluster1_Id", &m_truth_Cluster1_Id); */ 0364 /* truth_particles->Branch("truth_Cluster1_E", &m_truth_Cluster1_E); */ 0365 /* truth_particles->Branch("truth_Cluster1_Ecore", &m_truth_Cluster1_Ecore); */ 0366 /* truth_particles->Branch("truth_Cluster1_Eta", &m_truth_Cluster1_Eta); */ 0367 /* truth_particles->Branch("truth_Cluster1_Phi", &m_truth_Cluster1_Phi); */ 0368 /* truth_particles->Branch("truth_Cluster1_Pt", &m_truth_Cluster1_Pt); */ 0369 /* truth_particles->Branch("truth_Cluster1_xF", &m_truth_Cluster1_xF); */ 0370 /* truth_particles->Branch("truth_Cluster1_Chi2", &m_truth_Cluster1_Chi2); */ 0371 /* truth_particles->Branch("truth_Cluster1_Decay1EnergyFraction", &m_truth_Cluster1_Decay1EnergyFraction); */ 0372 /* truth_particles->Branch("truth_Cluster1_Decay2EnergyFraction", &m_truth_Cluster1_Decay2EnergyFraction); */ 0373 /* truth_particles->Branch("truth_Cluster2_Id", &m_truth_Cluster2_Id); */ 0374 /* truth_particles->Branch("truth_Cluster2_E", &m_truth_Cluster2_E); */ 0375 /* truth_particles->Branch("truth_Cluster2_Ecore", &m_truth_Cluster2_Ecore); */ 0376 /* truth_particles->Branch("truth_Cluster2_Eta", &m_truth_Cluster2_Eta); */ 0377 /* truth_particles->Branch("truth_Cluster2_Phi", &m_truth_Cluster2_Phi); */ 0378 /* truth_particles->Branch("truth_Cluster2_Pt", &m_truth_Cluster2_Pt); */ 0379 /* truth_particles->Branch("truth_Cluster2_xF", &m_truth_Cluster2_xF); */ 0380 /* truth_particles->Branch("truth_Cluster2_Chi2", &m_truth_Cluster2_Chi2); */ 0381 /* truth_particles->Branch("truth_Cluster2_Decay1EnergyFraction", &m_truth_Cluster2_Decay1EnergyFraction); */ 0382 /* truth_particles->Branch("truth_Cluster2_Decay2EnergyFraction", &m_truth_Cluster2_Decay2EnergyFraction); */ 0383 } 0384 0385 /* clusters_Towers->MakeClass("cluster"); */ 0386 /* truth_particles->MakeClass("truthParticles"); */ 0387 0388 return Fun4AllReturnCodes::EVENT_OK; 0389 } 0390 0391 //____________________________________________________________________________.. 0392 int pythiaEMCalAna::InitRun(PHCompositeNode *topNode) 0393 { 0394 std::cout << "pythiaEMCalAna::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl; 0395 return Fun4AllReturnCodes::EVENT_OK; 0396 } 0397 0398 //____________________________________________________________________________.. 0399 int pythiaEMCalAna::process_event(PHCompositeNode *topNode) 0400 { 0401 /* std::cout << "pythiaEMCalAna::process_event(PHCompositeNode *topNode) Processing event..." << std::endl; */ 0402 n_events_total++; 0403 0404 //Tower geometry node for location information 0405 RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC"); 0406 if (!towergeom) 0407 { 0408 std::cout << PHWHERE << "singlePhotonClusterAna::process_event Could not find node TOWERGEOM_CEMC" << std::endl; 0409 return Fun4AllReturnCodes::ABORTEVENT; 0410 } 0411 0412 //Raw tower information 0413 RawTowerContainer *towerContainer = nullptr; 0414 TowerInfoContainer *towerInfoContainer = nullptr; 0415 // again, node has different names in MC and RD 0416 if (isMonteCarlo) { 0417 towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC"); 0418 /* std::cout << "isMonteCarlo=" << isMonteCarlo << ", towerContainer=" << towerContainer << "\n"; */ 0419 /* else { */ 0420 /* std::cout << "Greg info: isMonteCarlo is false\n"; */ 0421 /* towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWERINFO_CALIB_CEMC"); */ 0422 /* std::cout << "Greg info: towerContainer is " << towerContainer << std::endl; */ 0423 /* } */ 0424 if(!towerContainer) { 0425 std::cout << PHWHERE << "pythiaEMCalAna::process_event Could not find node TOWER_CALIB_CEMC" << std::endl; 0426 return Fun4AllReturnCodes::ABORTEVENT; 0427 } 0428 } 0429 else { 0430 towerInfoContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC"); 0431 if (!towerInfoContainer) { 0432 std::cout << PHWHERE << "pythiaEMCalAna::process_event Could not find node TOWERINFO_CALIB_CEMC" << std::endl; 0433 return Fun4AllReturnCodes::ABORTEVENT; 0434 } 0435 } 0436 0437 0438 CaloEvalStack *caloevalstack = nullptr; 0439 /* CaloRawClusterEval *clustereval = nullptr; */ 0440 if (isMonteCarlo) { 0441 caloevalstack = new CaloEvalStack(topNode, "CEMC"); 0442 caloevalstack->next_event(topNode); 0443 m_clusterEval = caloevalstack->get_rawcluster_eval(); 0444 if(!m_clusterEval) 0445 { 0446 std::cout << PHWHERE << "pythiaEMCalAna::process_event cluster eval not available" << std::endl; 0447 return Fun4AllReturnCodes::ABORTEVENT; 0448 } 0449 } 0450 0451 //Information on clusters 0452 /* RawClusterContainer *clusterContainer = nullptr; */ 0453 // Name of node is different in MC and RD 0454 if (isMonteCarlo) { 0455 /* m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC"); */ 0456 m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC"); 0457 } 0458 else { 0459 m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTERINFO_CEMC"); 0460 /* m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTERINFO_POS_COR_CEMC"); */ 0461 } 0462 if(!m_clusterContainer) 0463 { 0464 if (isMonteCarlo) std::cout << PHWHERE << "pythiaEMCalAna::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl; 0465 else std::cout << PHWHERE << "pythiaEMCalAna::process_event - Fatal Error - CLUSTERINFO_CEMC node is missing. " << std::endl; 0466 return 0; 0467 } 0468 0469 // Truth information 0470 /* PHG4TruthInfoContainer *truthinfo = nullptr; */ 0471 if (isMonteCarlo) { 0472 m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo"); 0473 if(!m_truthInfo) { 0474 std::cout << PHWHERE << "pythiaEMCalAna::process_event Could not find node G4TruthInfo" << std::endl; 0475 return Fun4AllReturnCodes::ABORTEVENT; 0476 } 0477 } 0478 0479 // Minbias info 0480 if (!isMonteCarlo) { 0481 MinimumBiasInfo *minBiasInfo = findNode::getClass<MinimumBiasInfo>(topNode,"MinimumBiasInfo"); 0482 bool isMinBias = (minBiasInfo) ? minBiasInfo->isAuAuMinimumBias() : false; 0483 if (!isMinBias) { 0484 /* std::cout << PHWHERE << "pythiaEMCalAna::process_event Event is not minbias" << std::endl; */ 0485 return Fun4AllReturnCodes::ABORTEVENT; 0486 } 0487 } 0488 n_events_minbias++; 0489 0490 //Vertex information 0491 GlobalVertex* gVtx = nullptr; 0492 PHG4VtxPoint* mcVtx = nullptr; 0493 // Problem is MC has a PrimaryVtx but no GlobalVertex, while RD has the opposite 0494 if (isMonteCarlo) { 0495 PHG4TruthInfoContainer::VtxRange vtx_range = m_truthInfo->GetPrimaryVtxRange(); 0496 PHG4TruthInfoContainer::ConstVtxIterator vtxIter = vtx_range.first; 0497 mcVtx = vtxIter->second; 0498 n_events_with_vertex++; 0499 if (abs(mcVtx->get_z()) > 10.0) { 0500 /* std::cout << "Greg info: vertex z is " << mcVtx->get_z() << "\n"; */ 0501 return Fun4AllReturnCodes::ABORTEVENT; 0502 } 0503 n_events_with_good_vertex++; 0504 } 0505 else { 0506 GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap"); 0507 if (!vtxContainer) 0508 { 0509 std::cout << PHWHERE << "pythiaEMCalAna::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; 0510 assert(vtxContainer); // force quit 0511 return 0; 0512 } 0513 if (vtxContainer->empty()) 0514 { 0515 std::cout << PHWHERE << "pythiaEMCalAna::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; 0516 return Fun4AllReturnCodes::ABORTEVENT; 0517 } 0518 0519 //More vertex information 0520 else { 0521 gVtx = vtxContainer->begin()->second; 0522 if(!gVtx) 0523 { 0524 std::cout << PHWHERE << "pythiaEMCalAna::process_event Could not find vtx from vtxContainer" << std::endl; 0525 return Fun4AllReturnCodes::ABORTEVENT; 0526 } 0527 n_events_with_vertex++; 0528 // Require vertex to be within 10cm of 0 0529 if (abs(gVtx->get_z()) > 10.0) { 0530 return Fun4AllReturnCodes::ABORTEVENT; 0531 } 0532 n_events_with_good_vertex++; 0533 } 0534 } 0535 0536 /* PHG4TruthInfoContainer::Range truthRange = m_truthInfo->GetPrimaryParticleRange(); */ 0537 /* PHG4TruthInfoContainer::ConstIterator truthIter; */ 0538 /* //from the HepMC event log */ 0539 /* for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) */ 0540 /* { */ 0541 /* //PHG4Particle* part = m_truthInfo->GetParticle(truthIter->second->get_trkid()) */ 0542 0543 //For eventgen ancestory information 0544 PHHepMCGenEventMap *genEventMap = nullptr; 0545 if (isMonteCarlo) { 0546 genEventMap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap"); 0547 if(!genEventMap) 0548 { 0549 std::cout << PHWHERE << "pythiaEMCalAna::process_event Could not find PHHepMCGenEventMap" << std::endl; 0550 return Fun4AllReturnCodes::ABORTEVENT; 0551 } 0552 /* std::cout << "genEventMap is " << genEventMap << " with size " << genEventMap->size() << "; more details:\n"; */ 0553 /* genEventMap->identify(); */ 0554 } 0555 0556 //event level information of the above 0557 PHHepMCGenEvent *genEvent = nullptr; // = genEventMap -> get(getEvent); 0558 /* if(!genEvent) */ 0559 /* { */ 0560 /* std::cout << PHWHERE << "pythiaEMCalAna::process_event Could not find PHHepMCGenEvent" << std::endl; */ 0561 /* return Fun4AllReturnCodes::ABORTEVENT; */ 0562 /* } */ 0563 /* std::cout << "genEvent is " << genEvent << " with size " << genEvent->size() << "; more details:\n"; */ 0564 /* genEvent->identify(); */ 0565 0566 /* HepMC::GenEvent *theEvent = nullptr; // = genEvent -> getEvent(); */ 0567 0568 //grab all the towers and fill their energies. 0569 bool write_towers = true; 0570 if (write_towers) { 0571 if (isMonteCarlo) { 0572 RawTowerContainer::ConstRange tower_range = towerContainer -> getTowers(); 0573 for(RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter!= tower_range.second; tower_iter++) 0574 { 0575 int phibin = tower_iter -> second -> get_binphi(); 0576 int etabin = tower_iter -> second -> get_bineta(); 0577 double phi = towergeom -> get_phicenter(phibin); 0578 double eta = towergeom -> get_etacenter(etabin); 0579 double energy = tower_iter -> second -> get_energy(); 0580 0581 m_tower_energy.push_back(energy); 0582 m_eta_center.push_back(eta); 0583 m_phi_center.push_back(phi); 0584 } 0585 n_events_positiveCaloE++; 0586 } 0587 else { 0588 unsigned int tower_range = towerInfoContainer->size(); 0589 double totalCaloE = 0; 0590 int goodTowerCtr = 0; 0591 for(unsigned int iter = 0; iter < tower_range; iter++) 0592 { 0593 unsigned int towerkey = towerInfoContainer->encode_key(iter); 0594 unsigned int etabin = towerInfoContainer->getTowerEtaBin(towerkey); 0595 unsigned int phibin = towerInfoContainer->getTowerPhiBin(towerkey); 0596 double eta = towergeom -> get_etacenter(etabin); 0597 double phi = towergeom -> get_phicenter(phibin); 0598 0599 TowerInfo* tower = towerInfoContainer->get_tower_at_channel(iter); 0600 // check if tower is good 0601 if(!tower->get_isGood()) continue; 0602 double energy = tower->get_energy(); 0603 goodTowerCtr++; 0604 totalCaloE += energy; 0605 m_tower_energy.push_back(energy); 0606 m_eta_center.push_back(eta); 0607 m_phi_center.push_back(phi); 0608 } 0609 if (totalCaloE < 0) { 0610 return Fun4AllReturnCodes::ABORTEVENT; 0611 } 0612 n_events_positiveCaloE++; 0613 } 0614 } // done writing tower info 0615 0616 //grab all the clusters in the event 0617 bool write_clusters = true; 0618 bool apply_energy_cut = true; 0619 float clusterMinEnergyCut = 0.3; // in GeV 0620 if (write_clusters) { 0621 RawClusterContainer::ConstRange clusterEnd = m_clusterContainer -> getClusters(); 0622 RawClusterContainer::ConstIterator clusterIter; 0623 /* std::cout << "\n\nBeginning cluster loop\n"; */ 0624 for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++) 0625 { 0626 RawCluster *recoCluster = clusterIter -> second; 0627 CLHEP::Hep3Vector vertex; 0628 if (isMonteCarlo) vertex = CLHEP::Hep3Vector(mcVtx->get_x(), mcVtx->get_y(), mcVtx->get_z()); 0629 else { 0630 if (gVtx) vertex = CLHEP::Hep3Vector(gVtx->get_x(), gVtx->get_y(), gVtx->get_z()); 0631 else vertex = CLHEP::Hep3Vector(0,0,0); 0632 } 0633 /* vertex = CLHEP::Hep3Vector(0., 0., 0.); */ 0634 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex); 0635 float clusE = recoCluster->get_energy(); 0636 float clusEcore = recoCluster->get_ecore(); 0637 float clus_eta = E_vec_cluster.pseudoRapidity(); 0638 float clus_phi = E_vec_cluster.phi(); 0639 0640 if (apply_energy_cut) { 0641 if (clusE < clusterMinEnergyCut) continue; 0642 } 0643 0644 m_cluster_ID.push_back(recoCluster->get_id()); 0645 m_cluster_e.push_back(clusE); 0646 m_cluster_eta.push_back(clus_eta); 0647 m_cluster_phi.push_back(clus_phi); 0648 m_cluster_ecore.push_back(clusEcore); 0649 m_cluster_chi2.push_back(recoCluster -> get_chi2()); 0650 m_cluster_prob.push_back(recoCluster -> get_prob()); 0651 m_cluster_nTowers.push_back(recoCluster -> getNTowers()); 0652 // get all the towers in this cluster 0653 std::vector<float> allE; 0654 std::vector<float> allEta; 0655 std::vector<float> allPhi; 0656 RawCluster::TowerConstRange tower_range = recoCluster->get_towers(); 0657 RawCluster::TowerConstIterator tower_iter; 0658 for (tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++) { 0659 double E = tower_iter->second; 0660 allE.push_back(E); 0661 unsigned int towerkey = tower_iter->first; 0662 /* unsigned int etabin = towerInfoContainer->getTowerEtaBin(towerkey); */ 0663 unsigned int etabin = RawTowerDefs::decode_index1(towerkey); 0664 /* unsigned int phibin = towerInfoContainer->getTowerPhiBin(towerkey); */ 0665 unsigned int phibin = RawTowerDefs::decode_index2(towerkey); 0666 double eta = towergeom -> get_etacenter(etabin); 0667 double phi = towergeom -> get_phicenter(phibin); 0668 allEta.push_back(eta); 0669 allPhi.push_back(phi); 0670 } 0671 m_cluster_allTowersE.push_back(allE); 0672 m_cluster_allTowersEta.push_back(allEta); 0673 m_cluster_allTowersPhi.push_back(allPhi); 0674 0675 0676 if (isMonteCarlo) { 0677 /* std::cout << "\nEntering cluster isMonteCarlo for cluster "; */ 0678 /* recoCluster->identify(); */ 0679 std::set<PHG4Particle*> all_parts = m_clusterEval->all_truth_primary_particles(recoCluster); 0680 // first add all relevant particles to allTrackIDs for later 0681 for (auto& part : all_parts) { 0682 allTrackIDs.insert(part->get_track_id()); 0683 } 0684 if (all_parts.empty()) { 0685 /* std::cout << "In cluster association, all_parts is empty. recoCluster is "; */ 0686 /* recoCluster->identify(); */ 0687 // noise cluster without any real particles creating the shower 0688 /* m_cluster_maxE_trackID.push_back(-9999); */ 0689 /* m_cluster_maxE_PID.push_back(-9999); */ 0690 /* std::vector<float> all_cluster_trackIDs; */ 0691 /* m_cluster_all_trackIDs.push_back(all_cluster_trackIDs); */ 0692 m_cluster_nParticles.push_back(-9999); 0693 m_cluster_primaryParticlePid.push_back(-9999); 0694 std::vector<float> allPids; 0695 m_cluster_allSecondaryPids.push_back(allPids); 0696 } 0697 else { 0698 /* std::cout << "In cluster association, all_parts has " << all_parts.size() << " entries\n"; */ 0699 PHG4Particle* maxE_part = m_clusterEval->max_truth_primary_particle_by_energy(recoCluster); 0700 assert(maxE_part); 0701 std::vector<float> allPids; 0702 int track_id = maxE_part->get_track_id(); 0703 int embedID = m_truthInfo->isEmbeded(track_id); 0704 genEvent = genEventMap -> get(embedID); 0705 m_theEvent = genEvent -> getEvent(); 0706 /* std::cout << "maxE_part is "; */ 0707 /* maxE_part->identify(); */ 0708 if (maxE_part->get_pid() == 22) { 0709 // "primary" is a photon 0710 if (isDirectPhoton(maxE_part)) { 0711 /* std::cout << "maxE_part is a direct photon, adding primary ID\n"; */ 0712 m_cluster_primaryParticlePid.push_back(22); 0713 } 0714 else { 0715 // it's a decay photon, get the parent form pythia 0716 /* std::cout << "maxE_part is a decay photon\n"; */ 0717 HepMC::GenParticle* parent = GetHepMCParent(maxE_part); 0718 /* std::cout << "Parent is "; */ 0719 /* parent->print(); */ 0720 m_cluster_primaryParticlePid.push_back(parent->pdg_id()); 0721 } 0722 } // end photon check 0723 else { 0724 // if it's not a photon, add the PID directly 0725 /* std::cout << "maxE_part is not a photon, adding it as primary\n"; */ 0726 m_cluster_primaryParticlePid.push_back(maxE_part->get_pid()); 0727 } 0728 // now get the secondaries that produced the shower(s) 0729 std::set<PHG4Shower*> showers = m_clusterEval->all_truth_primary_showers(recoCluster); 0730 /* std::cout << "Cluster " << recoCluster << " has " << showers.size() << " showers\n"; */ 0731 m_cluster_nParticles.push_back(showers.size()); 0732 for (auto& shower : showers) { 0733 int parent_id = shower->get_parent_particle_id(); 0734 PHG4Particle* par = m_truthInfo->GetParticle(parent_id); 0735 /* std::cout << "Shower-producing particle is "; */ 0736 /* par->identify(); */ 0737 allPids.push_back(par->get_pid()); 0738 } 0739 m_cluster_allSecondaryPids.push_back(allPids); 0740 0741 } // end case where all_parts isn't empty 0742 } // end MC-specific stuff 0743 } // end loop over clusters 0744 } 0745 0746 // Next grab the truth particles 0747 if (isMonteCarlo) { 0748 PHG4TruthInfoContainer::Range truthRange = m_truthInfo->GetPrimaryParticleRange(); 0749 PHG4TruthInfoContainer::ConstIterator truthIter; 0750 /* std::cout << "\n\nGreg info: event " << n_events_total << ", starting loop over Geant primary particles\n"; */ 0751 for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) 0752 { 0753 /* std::cout << "Starting new primary. n_primaries=" << n_primaries << "\n"; */ 0754 PHG4Particle *truthPar = truthIter->second; 0755 int track_id = truthPar->get_track_id(); 0756 // only interested in particles which produced a cluster 0757 if (allTrackIDs.count(track_id) == 0) continue; 0758 /* std::cout << "\ttruthPar PID=" << truthPar->get_pid() << "; adding 1 primary\n"; */ 0759 n_primaries++; 0760 /* std::cout << "\nIn main loop, n_primaries=" << n_primaries << ", truthPar is "; */ 0761 /* truthPar->identify(); */ 0762 int embedID = m_truthInfo->isEmbeded(track_id); 0763 genEvent = genEventMap -> get(embedID); 0764 m_theEvent = genEvent -> getEvent(); 0765 int pid = truthPar->get_pid(); 0766 if (abs(pid) > 99) { 0767 /* std::cout << "Geant-decayed hadron. PID=" << pid << "\n"; */ 0768 /* truthPar->identify(); */ 0769 /* std::cout << "IsNeutralHadron=" << IsNeutralHadron(abs(pid)) << ", IsChargedHadron=" << IsChargedHadron(abs(pid)) << "\n"; */ 0770 // if it's not a photon, add it as a primary 0771 /* addPrimaryFromGeant(truthPar); */ 0772 if (IsNeutralHadron(abs(pid))) { 0773 n_neutral_hadrons++; 0774 n_neutral_hadrons_geant++; 0775 } 0776 if (IsChargedHadron(abs(pid))) { 0777 n_charged_hadrons++; 0778 n_charged_hadrons_geant++; 0779 } 0780 /* PrintParent(truthPar); */ 0781 if (pid == 111 || pid == 221) { 0782 /* if (false) { */ 0783 std::vector<PHG4Particle*> daughters = GetAllDaughters(truthPar); 0784 if (daughters.size() == 2) { 0785 if (!withinAcceptance(truthPar)) continue; 0786 /* std::cout << "Greg info: found a geant primary pi0/eta with 2 daughters:\n"; */ 0787 /* PrintParent(truthPar); */ 0788 FillParent(truthPar); 0789 PHG4Particle* decay1 = nullptr; 0790 PHG4Particle* decay2 = nullptr; 0791 GetBestDaughters(truthPar, decay1, decay2); 0792 /* std::cout << "After GetBestDaughters: decay1=" << decay1 << "; decay2=" << decay2 << "\n"; */ 0793 /* std::cout << "Filling decay1 (geant decay)\n"; */ 0794 FillDecay("decay1", decay1, truthPar); 0795 FillDecay("decay2", decay2, truthPar); 0796 } 0797 } 0798 } 0799 if (abs(pid) < 20) n_leptons++; 0800 if (pid == 22) { 0801 n_primary_photons++; 0802 // if it is a photon, we have to check whether it's truly primary, 0803 // or a decay photon handled in pythia 0804 if (isDirectPhoton(truthPar)) { 0805 // if it's really a direct photon, add it normally 0806 /* addPrimaryFromGeant(truthPar); */ 0807 if (!withinAcceptance(truthPar)) continue; 0808 FillParent(truthPar); 0809 /* std::cout << "Filling decay1 (direct photon)\n"; */ 0810 FillDecay("decay1", nullptr, truthPar); 0811 FillDecay("decay2", nullptr, truthPar); 0812 n_direct_photons++; 0813 } 0814 else { 0815 // "primary" photon is actually a decay photon from a decay handled by pythia 0816 /* std::cout << "\tIt's a decay photon, subtracting 1 primary\n"; */ 0817 n_primaries--; 0818 // find the parent in pythia and add that as the true primary 0819 HepMC::GenParticle* parent = GetHepMCParent(truthPar); 0820 /* std::cout << "\tnot a direct photon -- parent is\n"; */ 0821 /* parent->print(); */ 0822 // Check if we've already added this primary via another "primary" decay photon 0823 if (primaryBarcodes.count(parent->barcode()) != 0) { 0824 /* std::cout << "Already added HepMC primary " << parent << ":\n"; */ 0825 /* parent->print(); */ 0826 continue; 0827 } 0828 primaryBarcodes.insert(parent->barcode()); 0829 /* std::cout << "\tAdding parent as 1 primary\n"; */ 0830 n_primaries++; 0831 pid = parent->pdg_id(); 0832 if (IsNeutralHadron(abs(pid))) { 0833 n_neutral_hadrons++; 0834 n_neutral_hadrons_pythia++; 0835 n_pythia_decayed_hadrons++; 0836 } 0837 if (IsChargedHadron(abs(pid))) { 0838 n_charged_hadrons++; 0839 n_charged_hadrons_pythia++; 0840 n_pythia_decayed_hadrons++; 0841 } 0842 if (parent->status() == 2 && abs(parent->pdg_id()) >= 100) { 0843 // truthPar is a decay photon and parent is the primary hadron 0844 /* addPrimaryFromPythia((*parent)); */ 0845 if (pid == 111 || pid == 221) { 0846 /* if (!withinAcceptance(parent)) continue; */ 0847 /* std::cout << "Adding particles via FillTruth (pythia decay)\n"; */ 0848 FillTruth(parent); 0849 } 0850 } 0851 else { 0852 // weird case?? 0853 std::cout << "\nGreg info: weird combination of photon or parent status or PID... photon:\n"; 0854 truthPar->identify(); 0855 std::cout << "Parent:\n"; 0856 parent->print(); 0857 std::cout << "\n"; 0858 } 0859 } // end decay photon case 0860 } // end photon check 0861 /* std::cout << "\tEnding primary. n_primaries=" << n_primaries << "\n"; */ 0862 } // end loop over primary particles 0863 } // end grabbing truth info 0864 0865 0866 clusters_Towers -> Fill(); 0867 if (isMonteCarlo) { 0868 truth_particles -> Fill(); 0869 delete caloevalstack; 0870 } 0871 0872 // Below is some example code showing cluster->particle matching 0873 0874 /* // Start with clusters, and find the corresponding truth particle(s) that */ 0875 /* // produced the shower */ 0876 /* RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters(); */ 0877 /* RawClusterContainer::ConstIterator clusterIter; */ 0878 /* for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++) */ 0879 /* { */ 0880 /* RawCluster *recoCluster = clusterIter -> second; */ 0881 /* CLHEP::Hep3Vector vertex = CLHEP::Hep3Vector(mcVtx->get_x(), mcVtx->get_y(), mcVtx->get_z()); */ 0882 /* CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex); */ 0883 /* float clusE = E_vec_cluster.mag(); */ 0884 /* // Ignore clusters below 300MeV */ 0885 /* if (clusE < 0.3) continue; */ 0886 0887 /* // all particles that contributed energy to this cluster */ 0888 /* std::set<PHG4Particle*> all_parts = clustereval->all_truth_primary_particles(recoCluster); */ 0889 /* // need at least one particle; if not, cluster is just noise */ 0890 /* if (all_parts.size() < 1) continue; */ 0891 /* // particle which contributed the most energy to this cluster */ 0892 /* PHG4Particle* max_part = clustereval->max_truth_primary_particle_by_energy(recoCluster); */ 0893 /* /1* if (max_part->get_pid() > 0) continue; *1/ */ 0894 /* std::cout << "Greg info: printing cluster info:\n"; */ 0895 /* std::cout << Form("id=%d, (E, phi, eta) = (%f, %f, %f)\n", recoCluster->get_id(), clusE, E_vec_cluster.phi(), E_vec_cluster.pseudoRapidity()); */ 0896 /* /1* recoCluster->identify(); *1/ */ 0897 /* std::cout << "Corresponding particle info: (all_parts size=" << all_parts.size() << ", max_part=" << max_part << ")\n"; */ 0898 /* max_part->identify(); */ 0899 /* std::set<RawCluster*> max_part_clusters = clustereval->all_clusters_from(max_part); */ 0900 /* std::cout << "max_part contributed to " << max_part_clusters.size() << " clusters; printing those clusters\n"; */ 0901 /* for (auto& cl : max_part_clusters) { */ 0902 /* cl->identify(); */ 0903 /* } */ 0904 /* std::cout << "\n"; */ 0905 /* } */ 0906 0907 0908 // Below shows how to loop over all the different HepMC events in this Fun4All event 0909 0910 /* PHHepMCGenEventMap::Iter genEventIter; */ 0911 /* for (genEventIter = genEventMap->begin(); genEventIter != genEventMap->end(); genEventIter++) { */ 0912 /* int embedID = genEventIter->first; */ 0913 /* // For HIJING embedded samples, the PYTHIA event has a positive embedID */ 0914 /* if (hasHIJING && embedID < 1) continue; */ 0915 /* if (embedID < 1) continue; */ 0916 /* genEvent = genEventIter->second; */ 0917 /* /1* std::cout << "Greg info: embedID=" << embedID << "; printing genEvent\n"; *1/ */ 0918 /* /1* genEvent->identify(); *1/ */ 0919 /* theEvent = genEvent->getEvent(); */ 0920 /* /1* std::cout << "Greg info: printing theEvent\n"; *1/ */ 0921 /* /1* theEvent->print(); *1/ */ 0922 /* /1* std::cout << Form("Greg info: embedID=%d, genEvent->is_simulated()=%d\n", embedID, genEvent->is_simulated()); *1/ */ 0923 /* // now loop over PYTHIA particles */ 0924 /* for(HepMC::GenEvent::particle_const_iterator p = theEvent -> particles_begin(); p != theEvent -> particles_end(); ++p) */ 0925 /* { */ 0926 /* assert(*p); */ 0927 /* // look for photons only */ 0928 /* if ((*p)->pdg_id() == 22) { */ 0929 /* // only want final-state particles, i.e. status = 1 */ 0930 /* if ((*p)->status() != 1) continue; */ 0931 /* // get the photon's parent */ 0932 /* HepMC::GenVertex* prod_vtx = (*p)->production_vertex(); */ 0933 /* if (prod_vtx->particles_in_size() == 1) { */ 0934 /* HepMC::GenVertex::particles_in_const_iterator parent = prod_vtx->particles_in_const_begin(); */ 0935 /* assert((*parent)); */ 0936 /* // parent's status and pdg_id tell us whether it's a direct photon */ 0937 /* // or a decay photon */ 0938 /* if ((*parent)->status() > 2 && abs((*parent)->pdg_id()) < 100) { */ 0939 /* // direct photon */ 0940 /* /1* addPrimaryFromPythia((*p)); *1/ */ 0941 /* n_pythia_direct_photons++; */ 0942 /* } */ 0943 /* else if ((*parent)->status() == 2 && abs((*parent)->pdg_id()) >= 100) { */ 0944 /* // decay photon */ 0945 /* n_pythia_decay_photons++; */ 0946 /* if (withinAcceptance((*p))) n_pythia_decay_photons_in_acceptance++; */ 0947 /* std::pair<int,int> embedBarcode = std::pair<int,int>(embedID, (*parent)->barcode()); */ 0948 /* if (! vector_contains(embedBarcode, pythia_primary_barcodes) ) { */ 0949 /* pythia_primary_barcodes.push_back(embedBarcode); */ 0950 /* n_pythia_decays++; */ 0951 /* if (withinAcceptance((*parent))) n_pythia_decays_in_acceptance++; */ 0952 /* if ( (*parent)->pdg_id() == 111 ) { */ 0953 /* n_pythia_decayed_pi0s++; */ 0954 /* if (withinAcceptance((*parent))) n_pythia_decayed_pi0s_in_acceptance++; */ 0955 /* } */ 0956 /* } */ 0957 0958 /* } */ 0959 /* else { */ 0960 /* // weird case?? */ 0961 /* std::cout << "\nGreg info: weird combination of photon or parent status or PID... photon:\n"; */ 0962 /* (*p)->print(); */ 0963 /* std::cout << "Parent:\n"; */ 0964 /* (*parent)->print(); */ 0965 /* std::cout << "\n"; */ 0966 /* } */ 0967 /* } */ 0968 /* else { */ 0969 /* // weird photon -- they should only have 1 parent */ 0970 /* std::cout << "\nGreg info: in PYTHIA check, found a photon with " << prod_vtx->particles_in_size() << " parent(s). Photon:\n"; */ 0971 /* (*p)->print(); */ 0972 /* for (HepMC::GenVertex::particles_in_const_iterator parent = prod_vtx->particles_in_const_begin(); parent != prod_vtx->particles_in_const_end(); parent++) { */ 0973 /* std::cout << "Parent:\n"; */ 0974 /* (*parent)->print(); */ 0975 /* } */ 0976 /* std::cout << "\n"; */ 0977 /* } */ 0978 /* } // end photon check */ 0979 /* else { */ 0980 /* // check for (final state) non-decayed pi0s */ 0981 /* if ( (*p)->status() == 1 && (*p)->pdg_id() > 100) { */ 0982 /* n_pythia_nondecayed_hadrons++; */ 0983 /* if (withinAcceptance((*p))) n_pythia_nondecayed_hadrons_in_acceptance++; */ 0984 /* if ( (*p)->pdg_id() == 111 ) { */ 0985 /* n_pythia_nondecayed_pi0s++; */ 0986 /* if (withinAcceptance((*p))) n_pythia_nondecayed_pi0s_in_acceptance++; */ 0987 /* } */ 0988 /* } */ 0989 /* } */ 0990 /* } // end PYTHIA loop */ 0991 /* } // end loop over embedID */ 0992 /* std::cout << "PYTHIA: " << n_pythia << " total particles, " << n_pho_pythia << " photons, " << n_pi0_pythia << " pi0s\n"; */ 0993 0994 0995 // Below is the old code I used to find truth photons 0996 // Warning! Likely will not work after uncommenting. Should only be used as an example 0997 0998 // Next check for decays handled by Geant 0999 // Start with primary particles from GEANT 1000 // Look for photons and decide if they are direct or decay photons 1001 /* PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange(); */ 1002 /* PHG4TruthInfoContainer::ConstIterator truthIter; */ 1003 /* /1* std::cout << "\n\nGreg info: starting loop over Geant primary particles\n"; *1/ */ 1004 /* for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) */ 1005 /* { */ 1006 /* PHG4Particle *truthPar = truthIter->second; */ 1007 /* int embedID = truthinfo->isEmbeded(truthPar->get_track_id()); */ 1008 /* if (hasHIJING && embedID < 1) continue; */ 1009 /* if (embedID < 1) continue; */ 1010 /* int geant_barcode = truthPar->get_barcode(); */ 1011 /* int pid = truthPar->get_pid(); */ 1012 /* // look for photons only */ 1013 /* if (pid == 22) { */ 1014 /* /1* std::cout << "\nFound a primary photon (embedID=" << embedID << ", barcode=" << geant_barcode << "); details:\t"; *1/ */ 1015 /* /1* truthPar->identify(); *1/ */ 1016 /* genEvent = genEventMap -> get(embedID); */ 1017 /* theEvent = genEvent -> getEvent(); */ 1018 /* // is it a direct photon? */ 1019 /* if (isDirectPhoton(truthPar, theEvent)) { */ 1020 /* /1* std::cout << "\tbarcode " << geant_barcode << " is a direct photon... adding to list\n"; *1/ */ 1021 /* n_direct_photons++; */ 1022 /* addDirectPhoton(truthPar, truthinfo); */ 1023 /* TLorentzVector truth_momentum; */ 1024 /* truth_momentum.SetPxPyPzE(truthPar->get_px(), truthPar->get_py(), truthPar->get_pz(), truthPar->get_e()); */ 1025 /* float eta = truth_momentum.PseudoRapidity(); */ 1026 /* if (abs(eta) <= 1.1) n_direct_photons_in_acceptance++; */ 1027 /* } */ 1028 /* else { */ 1029 /* // decay photon */ 1030 /* /1* std::cout << "\tbarcode " << geant_barcode << " is a decay photon... check parent info\n"; *1/ */ 1031 /* n_decay_photons++; */ 1032 /* if (withinAcceptance(truthPar)) n_decay_photons_in_acceptance++; */ 1033 /* addDecayPhoton(truthPar, truthinfo, theEvent); */ 1034 /* // find the parent and add it as a primary */ 1035 /* // primary photon means geant doesn't know about the parent */ 1036 /* // so match this photon to the pythia photon */ 1037 /* HepMC::GenParticle* pho = getGenParticle(geant_barcode, theEvent); */ 1038 /* /1* assert(pho); *1/ */ 1039 /* if (!pho) { */ 1040 /* // problem -- we couldn't find the corresponding pythia particle */ 1041 /* std::cout << "\t\tGreg info: skipping Geant primary with barcode " << geant_barcode << " because corresponding pythia particle could not be found; more details:\t"; */ 1042 /* truthPar->identify(); */ 1043 /* continue; */ 1044 /* } */ 1045 /* // now find the parent in pythia */ 1046 /* HepMC::GenVertex* prod_vtx = pho->production_vertex(); */ 1047 /* if (prod_vtx->particles_in_size() == 1) { */ 1048 /* HepMC::GenVertex::particles_in_const_iterator parent = prod_vtx->particles_in_const_begin(); */ 1049 /* assert((*parent)); */ 1050 /* /1* std::cout << "\t\tParent barcode=" << (*parent)->barcode() << ", id=" << (*parent)->pdg_id() << "... "; *1/ */ 1051 /* std::pair<int,int> embedBarcode(embedID, (*parent)->barcode()); */ 1052 /* if (! vector_contains(embedBarcode, primaryBarcodes) ) { */ 1053 /* /1* std::cout << "adding to primaries.\n"; *1/ */ 1054 /* primaryBarcodes.push_back(embedBarcode); */ 1055 /* n_primary++; */ 1056 /* if (withinAcceptance((*parent))) n_primary_in_acceptance++; */ 1057 /* /1* n_pythia_decays++; *1/ */ 1058 /* /1* std::cout << "adding primary from pythia. PID=" << (*parent)->pdg_id() << "\n"; *1/ */ 1059 /* addPrimaryHadronFromPythia((*parent), embedID); */ 1060 /* } */ 1061 /* else { */ 1062 /* /1* std::cout << "already added this primary.\n"; *1/ */ 1063 /* } */ 1064 /* } */ 1065 /* else { */ 1066 /* std::cout << "\t\tGreg info: pythia-decayed photon with " << prod_vtx->particles_in_size() << " parents. Skipping...\n"; */ 1067 /* } */ 1068 /* } */ 1069 /* } // end photon check */ 1070 /* else { */ 1071 /* if (pid > 100) { */ 1072 /* n_geant_primary_hadrons++; */ 1073 /* if (withinAcceptance(truthPar)) n_geant_primary_hadrons_in_acceptance++; */ 1074 /* } */ 1075 /* if (pid == 111) { */ 1076 /* n_geant_primary_pi0s++; */ 1077 /* if (withinAcceptance(truthPar)) n_geant_primary_pi0s_in_acceptance++; */ 1078 /* } */ 1079 /* } */ 1080 /* } // end loop over geant primary particles */ 1081 1082 1083 // Now loop over GEANT secondary particles, taking only the daughters 1084 // of primary particles we already found 1085 /* truthRange = truthinfo->GetSecondaryParticleRange(); */ 1086 /* /1* std::cout << "Greg info: starting loop over Geant secondary particles\n"; *1/ */ 1087 /* for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) */ 1088 /* { */ 1089 /* PHG4Particle *truthPar = truthIter->second; */ 1090 /* int embedID = truthinfo->isEmbeded(truthPar->get_track_id()); */ 1091 /* if (hasHIJING && embedID < 1) continue; */ 1092 /* if (embedID < 1) continue; */ 1093 /* // only looking for decay photons and their corresponding primaries */ 1094 /* if (truthPar->get_pid() == 22) { */ 1095 /* // get the parent */ 1096 /* PHG4Particle* parent = truthinfo->GetParticle(truthPar->get_parent_id()); */ 1097 /* // make sure parent is a primary */ 1098 /* if ( truthinfo->is_primary(parent) ) { */ 1099 /* // make sure the "parent" isn't also a photon, i.e. the */ 1100 /* // same photon before and after an interaction */ 1101 /* if (parent->get_pid() == 22) continue; */ 1102 1103 /* /1* std::cout << "\tFound a secondary photon (barcode=" << truthPar->get_barcode() << ")\n"; *1/ */ 1104 /* /1* std::cout << "\t\tParent is primary (barcode=" << parent->get_barcode() << ", id=" << parent->get_pid() << "), "; *1/ */ 1105 /* n_decay_photons++; */ 1106 /* if (withinAcceptance(truthPar)) n_decay_photons_in_acceptance++; */ 1107 /* n_geant_decay_photons++; */ 1108 /* if (withinAcceptance(truthPar)) n_geant_decay_photons_in_acceptance++; */ 1109 /* addDecayPhoton(truthPar, truthinfo, theEvent); */ 1110 /* int embedID = truthinfo->isEmbeded(parent->get_track_id()); */ 1111 /* std::pair<int,int> embedBarcode(embedID, parent->get_barcode()); */ 1112 /* /1* if (parent->get_pid() == 111) { *1/ */ 1113 /* /1* std::cout << "Found a primary pi0 with embedID " << embedID << ", barcode " << parent->get_barcode() << "\n"; *1/ */ 1114 /* /1* } *1/ */ 1115 /* if (! vector_contains(embedBarcode, primaryBarcodes) ) { */ 1116 /* /1* std::cout << "adding to primaries.\n"; *1/ */ 1117 /* primaryBarcodes.push_back(embedBarcode); */ 1118 /* n_primary++; */ 1119 /* if (withinAcceptance(parent)) n_primary_in_acceptance++; */ 1120 /* n_geant_decays++; */ 1121 /* if (withinAcceptance(parent)) n_geant_decays_in_acceptance++; */ 1122 1123 /* /1* std::cout << "\tadded primary with embedID " << embedID << ", barcode " << parent->get_barcode() << "\n"; *1/ */ 1124 /* /1* std::cout << "Adding primary from geant. PID=" << parent->get_pid() << "\n"; *1/ */ 1125 /* addPrimaryHadronFromGeant(parent, truthinfo); */ 1126 /* /1* std::cout << "\nFound a geant-decayed primary:\n"; *1/ */ 1127 /* /1* parent->identify(); *1/ */ 1128 /* /1* std::cout << "Daughter photon is\n"; *1/ */ 1129 /* /1* truthPar->identify(); *1/ */ 1130 /* /1* std::cout << "\n"; *1/ */ 1131 /* } */ 1132 /* else { */ 1133 /* /1* std::cout << "already added this primary.\n"; *1/ */ 1134 /* } */ 1135 /* } */ 1136 /* /1* else { *1/ */ 1137 /* /1* std::cout << "\t\tParent is *not* primary. Skipping!\n"; *1/ */ 1138 /* /1* } *1/ */ 1139 /* } //end photon check */ 1140 /* } // end loop over geant secondaries */ 1141 1142 1143 return Fun4AllReturnCodes::EVENT_OK; 1144 } 1145 1146 //____________________________________________________________________________.. 1147 int pythiaEMCalAna::ResetEvent(PHCompositeNode *topNode) 1148 { 1149 1150 m_tower_energy.clear(); 1151 m_eta_center.clear(); 1152 m_phi_center.clear(); 1153 m_cluster_ID.clear(); 1154 m_cluster_e.clear(); 1155 m_cluster_eta.clear(); 1156 m_cluster_phi.clear(); 1157 m_cluster_ecore.clear(); 1158 m_cluster_chi2.clear(); 1159 m_cluster_prob.clear(); 1160 m_cluster_nTowers.clear(); 1161 m_cluster_allTowersE.clear(); 1162 m_cluster_allTowersEta.clear(); 1163 m_cluster_allTowersPhi.clear(); 1164 m_cluster_nParticles.clear(); 1165 m_cluster_primaryParticlePid.clear(); 1166 m_cluster_allSecondaryPids.clear(); 1167 /* m_cluster_maxE_trackID.clear(); */ 1168 /* m_cluster_maxE_PID.clear(); */ 1169 /* m_cluster_all_trackIDs.clear(); */ 1170 m_truth_Parent_Barcode.clear(); 1171 m_truth_Parent_Pid.clear(); 1172 m_truth_Parent_M.clear(); 1173 m_truth_Parent_E.clear(); 1174 m_truth_Parent_Eta.clear(); 1175 m_truth_Parent_Phi.clear(); 1176 m_truth_Parent_Pt.clear(); 1177 m_truth_Parent_xF.clear(); 1178 m_truth_Parent_EndVtx_x.clear(); 1179 m_truth_Parent_EndVtx_y.clear(); 1180 m_truth_Parent_EndVtx_z.clear(); 1181 /* m_truth_Parent_EndVtx_t.clear(); */ 1182 m_truth_Parent_EndVtx_r.clear(); 1183 /* m_truth_Parent_TotalNDaughters.clear(); */ 1184 /* m_truth_Parent_AllDaughterPids.clear(); */ 1185 /* m_truth_Parent_AllDaughterEnergyFractions.clear(); */ 1186 m_truth_Parent_TotalNClusters.clear(); 1187 m_truth_Parent_AllClusterIDs.clear(); 1188 m_truth_Parent_AllClusterEnergyFractions.clear(); 1189 m_truth_Decay1_Barcode.clear(); 1190 m_truth_Decay1_Pid.clear(); 1191 m_truth_Decay1_M.clear(); 1192 m_truth_Decay1_E.clear(); 1193 m_truth_Decay1_Eta.clear(); 1194 m_truth_Decay1_Phi.clear(); 1195 m_truth_Decay1_Pt.clear(); 1196 m_truth_Decay1_xF.clear(); 1197 m_truth_Decay1_EndVtx_x.clear(); 1198 m_truth_Decay1_EndVtx_y.clear(); 1199 m_truth_Decay1_EndVtx_z.clear(); 1200 /* m_truth_Decay1_EndVtx_t.clear(); */ 1201 m_truth_Decay1_EndVtx_r.clear(); 1202 m_truth_Decay1_TotalNClusters.clear(); 1203 m_truth_Decay1_BestClusterID.clear(); 1204 m_truth_Decay1_BestClusterEfraction.clear(); 1205 /* m_truth_Decay1_AllClusterIDs.clear(); */ 1206 /* m_truth_Decay1_AllClusterEnergyFractions.clear(); */ 1207 m_truth_Decay2_Barcode.clear(); 1208 m_truth_Decay2_Pid.clear(); 1209 m_truth_Decay2_M.clear(); 1210 m_truth_Decay2_E.clear(); 1211 m_truth_Decay2_Eta.clear(); 1212 m_truth_Decay2_Phi.clear(); 1213 m_truth_Decay2_Pt.clear(); 1214 m_truth_Decay2_xF.clear(); 1215 m_truth_Decay2_EndVtx_x.clear(); 1216 m_truth_Decay2_EndVtx_y.clear(); 1217 m_truth_Decay2_EndVtx_z.clear(); 1218 /* m_truth_Decay2_EndVtx_t.clear(); */ 1219 m_truth_Decay2_EndVtx_r.clear(); 1220 m_truth_Decay2_TotalNClusters.clear(); 1221 m_truth_Decay2_BestClusterID.clear(); 1222 m_truth_Decay2_BestClusterEfraction.clear(); 1223 /* m_truth_Decay2_AllClusterIDs.clear(); */ 1224 /* m_truth_Decay2_AllClusterEnergyFractions.clear(); */ 1225 /* m_truth_Cluster1_Id.clear(); */ 1226 /* m_truth_Cluster1_E.clear(); */ 1227 /* m_truth_Cluster1_Ecore.clear(); */ 1228 /* m_truth_Cluster1_Eta.clear(); */ 1229 /* m_truth_Cluster1_Phi.clear(); */ 1230 /* m_truth_Cluster1_Pt.clear(); */ 1231 /* m_truth_Cluster1_xF.clear(); */ 1232 /* m_truth_Cluster1_Chi2.clear(); */ 1233 /* m_truth_Cluster1_Decay1EnergyFraction.clear(); */ 1234 /* m_truth_Cluster1_Decay2EnergyFraction.clear(); */ 1235 /* m_truth_Cluster2_Id.clear(); */ 1236 /* m_truth_Cluster2_E.clear(); */ 1237 /* m_truth_Cluster2_Ecore.clear(); */ 1238 /* m_truth_Cluster2_Eta.clear(); */ 1239 /* m_truth_Cluster2_Phi.clear(); */ 1240 /* m_truth_Cluster2_Pt.clear(); */ 1241 /* m_truth_Cluster2_xF.clear(); */ 1242 /* m_truth_Cluster2_Chi2.clear(); */ 1243 /* m_truth_Cluster2_Decay1EnergyFraction.clear(); */ 1244 /* m_truth_Cluster2_Decay2EnergyFraction.clear(); */ 1245 allTrackIDs.clear(); 1246 primaryBarcodes.clear(); 1247 1248 return Fun4AllReturnCodes::EVENT_OK; 1249 } 1250 1251 //____________________________________________________________________________.. 1252 int pythiaEMCalAna::EndRun(const int runnumber) 1253 { 1254 std::cout << "pythiaEMCalAna::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl; 1255 return Fun4AllReturnCodes::EVENT_OK; 1256 } 1257 1258 //____________________________________________________________________________.. 1259 int pythiaEMCalAna::End(PHCompositeNode *topNode) 1260 { 1261 std::cout << "pythiaEMCalAna::End(PHCompositeNode *topNode) This is the End..." << std::endl; 1262 std::cout << Form("Total # events: %ld\n", n_events_total); 1263 std::cout << Form("\t# minbias: %ld (%.2f%% total)\n", n_events_minbias, (float)(100.0*n_events_minbias/n_events_total)); 1264 std::cout << Form("\t# w/ vertex: %ld (%.2f%% total, %.2f%% minbias)\n", n_events_with_vertex, (float)(100.0*n_events_with_vertex/n_events_total), (float)(100.0*n_events_with_vertex/n_events_minbias)); 1265 std::cout << Form("\t# w/ good vertex: %ld (%.2f%% total, %.2f%% minbias, %.2f%% w/ vertex)\n", n_events_with_good_vertex, (float)(100.0*n_events_with_good_vertex/n_events_total), (float)(100.0*n_events_with_good_vertex/n_events_minbias), (float)(100.0*n_events_with_good_vertex/n_events_with_vertex)); 1266 std::cout << Form("\t# w/ positive calo E: %ld (%.2f%% total, %.2f%% minbias, %.2f%% w/ vertex, %.2f%% w/ good vertex)\n", n_events_positiveCaloE, (float)(100.0*n_events_positiveCaloE/n_events_total), (float)(100.0*n_events_positiveCaloE/n_events_minbias), (float)(100.0*n_events_positiveCaloE/n_events_with_vertex), (float)(100.0*n_events_positiveCaloE/n_events_with_good_vertex)); 1267 if (isMonteCarlo) { 1268 std::cout << Form("Total primary particles: %ld\n", n_primaries); 1269 std::cout << Form("Total primary hadrons decayed by pythia: %ld\n", n_pythia_decayed_hadrons); 1270 std::cout << Form("Total \"primary\" photons: %ld\n", n_primary_photons); 1271 std::cout << Form("Total direct photons: %ld\n", n_direct_photons); 1272 std::cout << Form("Total neutral hadrons: %ld (%ld geant, %ld pythia)\n", n_neutral_hadrons, n_neutral_hadrons_geant, n_neutral_hadrons_pythia); 1273 std::cout << Form("Total charged hadrons: %ld (%ld geant, %ld pythia)\n", n_charged_hadrons, n_charged_hadrons_geant, n_charged_hadrons_pythia); 1274 std::cout << Form("Total primary leptons: %ld\n", n_leptons); 1275 } 1276 fout -> cd(); 1277 std::cout << "Writing clusters_Towers\n"; 1278 clusters_Towers -> Write(); 1279 if (isMonteCarlo) { 1280 std::cout << "Writing truth_particles\n"; 1281 truth_particles -> Write(); 1282 } 1283 fout -> Close(); 1284 return Fun4AllReturnCodes::EVENT_OK; 1285 } 1286 1287 //____________________________________________________________________________.. 1288 int pythiaEMCalAna::Reset(PHCompositeNode *topNode) 1289 { 1290 std::cout << "pythiaEMCalAna::Reset(PHCompositeNode *topNode) being Reset" << std::endl; 1291 return Fun4AllReturnCodes::EVENT_OK; 1292 } 1293 1294 //____________________________________________________________________________.. 1295 void pythiaEMCalAna::Print(const std::string &what) const 1296 { 1297 std::cout << "pythiaEMCalAna::Print(const std::string &what) const Printing info for " << what << std::endl; 1298 } 1299 1300 /* void pythiaEMCalAna::addPrimaryFromGeant(PHG4Particle* part) { */ 1301 /* TLorentzVector truth_momentum; */ 1302 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 1303 /* if (! withinAcceptance(part)) return; */ 1304 /* int this_id = part->get_track_id(); */ 1305 /* n_primaries++; */ 1306 1307 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id); */ 1308 /* if (!end_vtx) { */ 1309 /* std::cout << "\nGreg info: no end vertex found for Geant-decayed primary:\n"; */ 1310 /* part->identify(); */ 1311 /* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); */ 1312 /* return; */ 1313 /* } */ 1314 1315 /* /1* int embedID = truthInfo->isEmbeded(part->get_track_id()); *1/ */ 1316 /* /1* m_truthBarcode.push_back(std::pair<int,int>(embedID, part->get_barcode())); *1/ */ 1317 /* /1* m_truthParentBarcode.push_back(std::pair<int,int>(embedID, 0)); *1/ */ 1318 /* if (part->get_pid() == 22) { */ 1319 /* /1* n_primary_photons++; *1/ */ 1320 /* // check for direct photons is done in process_event */ 1321 /* n_direct_photons++; */ 1322 /* } */ 1323 /* if (abs(part->get_pid()) < 20) n_leptons++; */ 1324 /* if (abs(part->get_pid()) > 100) { */ 1325 /* if (IsNeutralHadron(abs(part->get_pid()))) { */ 1326 /* n_neutral_hadrons++; */ 1327 /* n_neutral_hadrons_geant++; */ 1328 /* } */ 1329 /* else { */ 1330 /* n_charged_hadrons++; */ 1331 /* n_charged_hadrons_geant++; */ 1332 /* } */ 1333 /* } */ 1334 1335 /* m_truthIsPrimary.push_back(1); */ 1336 /* m_truthTrackID.push_back(this_id); */ 1337 /* m_truthPid.push_back(part->get_pid()); */ 1338 1339 /* m_truthE.push_back(truth_momentum.E()); */ 1340 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 1341 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 1342 /* m_truthPt.push_back(truth_momentum.Pt()); */ 1343 /* m_truthMass.push_back(truth_momentum.M()); */ 1344 1345 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 1346 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 1347 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 1348 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 1349 /* float x = end_vtx->get_x(); */ 1350 /* float y = end_vtx->get_y(); */ 1351 /* float r = sqrt(x*x + y*y); */ 1352 /* m_truthEndVtx_r.push_back(r); */ 1353 1354 /* std::vector<float> allClusterIDs; */ 1355 1356 /* std::set<RawCluster*> clusters = m_clusterEval->all_clusters_from(part); */ 1357 /* for (auto& cl : clusters) { */ 1358 /* allClusterIDs.push_back(cl->get_id()); */ 1359 /* } */ 1360 /* m_truth_all_clusterIDs.push_back(allClusterIDs); */ 1361 /* } */ 1362 1363 1364 /* void pythiaEMCalAna::addPrimaryFromPythia(HepMC::GenParticle* part) { */ 1365 /* /1* std::cout << "Entering pythiaEMCalAna::addPrimaryFromPythia\npart=" << part << "\n"; *1/ */ 1366 /* HepMC::FourVector mom = part->momentum(); */ 1367 /* TLorentzVector truth_momentum; */ 1368 /* truth_momentum.SetPxPyPzE(mom.px(), mom.py(), mom.pz(), mom.e()); */ 1369 /* if (! withinAcceptance(part)) return; */ 1370 1371 /* primaryBarcodes.insert(part->barcode()); */ 1372 /* n_primaries++; */ 1373 /* n_pythia_decayed_hadrons++; */ 1374 /* if (IsNeutralHadron(abs(part->pdg_id()))) { */ 1375 /* n_neutral_hadrons++; */ 1376 /* n_neutral_hadrons_pythia++; */ 1377 /* } */ 1378 /* else { */ 1379 /* n_charged_hadrons++; */ 1380 /* n_charged_hadrons_pythia++; */ 1381 /* } */ 1382 1383 /* m_truthIsPrimary.push_back(1); */ 1384 /* // no track ID available from pythia */ 1385 /* m_truthTrackID.push_back(-999); */ 1386 /* m_truthPid.push_back(part->pdg_id()); */ 1387 1388 /* m_truthE.push_back(truth_momentum.E()); */ 1389 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 1390 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 1391 /* m_truthPt.push_back(truth_momentum.Pt()); */ 1392 /* m_truthMass.push_back(truth_momentum.M()); */ 1393 1394 /* /1* std::cout << "Getting end vertex\n"; *1/ */ 1395 /* HepMC::GenVertex* end_vtx = part->end_vertex(); */ 1396 /* assert(end_vtx); */ 1397 /* /1* std::cout << "Found end vertex in pythia\n"; *1/ */ 1398 /* HepMC::FourVector position = end_vtx->position(); */ 1399 1400 /* m_truthEndVtx_x.push_back(position.x()); */ 1401 /* m_truthEndVtx_y.push_back(position.y()); */ 1402 /* m_truthEndVtx_z.push_back(position.z()); */ 1403 /* m_truthEndVtx_t.push_back(position.t()); */ 1404 /* float x = position.x(); */ 1405 /* float y = position.y(); */ 1406 /* float r = sqrt(x*x + y*y); */ 1407 /* /1* std::cout << Form("pythiaEMCalAna::addPrimaryFromPythia(): Vertex r=%f, perp=%f\n", r, position.perp()); *1/ */ 1408 /* m_truthEndVtx_r.push_back(r); */ 1409 1410 /* // no cluster matching available from pythia */ 1411 /* std::vector<float> allClusterIDs; */ 1412 /* m_truth_all_clusterIDs.push_back(allClusterIDs); */ 1413 /* } */ 1414 1415 bool pythiaEMCalAna::isDirectPhoton(PHG4Particle* part) { 1416 /* std::cout << "\tGreg info: entering isDirectPhoton(). "; */ 1417 /* std::cout << "G4 particle is " << part << ", barcode " << part->get_barcode() << "\n";//; printing info\n"; */ 1418 /* part->identify(); */ 1419 1420 // Only bother with this check if we have PYTHIA information 1421 // Otherwise assume any photon is a "direct" photon, since we can't tell from HIJING alone 1422 if (!hasPYTHIA) return (part->get_pid() == 22); 1423 HepMC::GenParticle* genpart = getGenParticle(part->get_barcode()); 1424 if (!genpart) { 1425 std::cout << "\t\tGreg info: in isDirectPhoton(), could not find pythia particle with barcode " << part->get_barcode() << "; returning true. (This may be an error!)\n"; 1426 return true; 1427 } 1428 assert(genpart); 1429 /* std::cout << "\tFound corresponding pythia particle; printing info\t"; */ 1430 /* genpart->print(); */ 1431 /* else std::cout << "Greg info: found corresponding pythia photon (" << genpart << ")\n"; */ 1432 HepMC::GenVertex* prod_vtx = genpart->production_vertex(); 1433 if (prod_vtx->particles_in_size() == 1) { 1434 HepMC::GenVertex::particles_in_const_iterator parent = prod_vtx->particles_in_const_begin(); 1435 assert((*parent)); 1436 // parent's status and pdg_id tell us whether it's a direct photon 1437 // or a decay photon 1438 if ((*parent)->status() > 2 && abs((*parent)->pdg_id()) < 100) { 1439 // direct photon 1440 /* std::cout << "\nGreg info: found a direct photon. Photon:\n"; */ 1441 /* genpart->print(); */ 1442 bool printHistory = false; 1443 if (printHistory) { 1444 // print history of all ancestors 1445 // generation 0 is genpart, gen -1 is its parent, etc. 1446 int generation = -1; 1447 std::cout << "\tGeneration " << generation << " -- "; 1448 (*parent)->print(); 1449 while (true) { 1450 generation--; 1451 HepMC::GenVertex::particles_in_const_iterator parentparent; 1452 HepMC::GenVertex* parent_prod_vtx = (*parent)->production_vertex(); 1453 if (parent_prod_vtx) { 1454 parentparent = parent_prod_vtx->particles_in_const_begin(); 1455 std::cout << "\tGeneration " << generation << ": "; 1456 (*parentparent)->print(); 1457 parent = parentparent; 1458 } 1459 else break; 1460 } 1461 parent = prod_vtx->particles_in_const_begin(); // reset parent to genpart's actual parent 1462 } 1463 return true; 1464 } 1465 else if ((*parent)->status() == 2 && abs((*parent)->pdg_id()) >= 100) { 1466 // decay photon 1467 /* std::cout << "\nGreg info: found a decay photon. Photon:\n"; */ 1468 /* genpart->print(); */ 1469 /* std::cout << "Parent:\n"; */ 1470 /* (*parent)->print(); */ 1471 /* std::cout << "\n"; */ 1472 return false; 1473 } 1474 else { 1475 // weird case?? 1476 std::cout << "\nGreg info: weird combination of photon or parent status or PID... photon:\n"; 1477 genpart->print(); 1478 std::cout << "Parent:\n"; 1479 (*parent)->print(); 1480 std::cout << "\n"; 1481 return false; 1482 } 1483 } // single parent check 1484 else { 1485 // weird photon -- they should only have 1 parent 1486 std::cout << "Greg info: found a photon with " << prod_vtx->particles_in_size() << " parent(s). Photon:\n"; 1487 genpart->print(); 1488 for (HepMC::GenVertex::particles_in_const_iterator parent = prod_vtx->particles_in_const_begin(); parent != prod_vtx->particles_in_const_end(); parent++) { 1489 std::cout << "Parent:\n"; 1490 (*parent)->print(); 1491 } 1492 /* std::cout << "\n"; */ 1493 return false; 1494 } 1495 } 1496 1497 void pythiaEMCalAna::FillTruthParticle(std::string which, PHG4Particle* par) { 1498 std::vector<float>* Barcode = nullptr; 1499 std::vector<float>* Pid = nullptr; 1500 std::vector<float>* M = nullptr; 1501 std::vector<float>* E = nullptr; 1502 std::vector<float>* Eta = nullptr; 1503 std::vector<float>* Phi = nullptr; 1504 std::vector<float>* Pt = nullptr; 1505 std::vector<float>* xF = nullptr; 1506 std::vector<float>* EndVtx_x = nullptr; 1507 std::vector<float>* EndVtx_y = nullptr; 1508 std::vector<float>* EndVtx_z = nullptr; 1509 /* std::vector<float>* EndVtx_t = nullptr; */ 1510 std::vector<float>* EndVtx_r = nullptr; 1511 1512 if (which == "parent") { 1513 Barcode = &m_truth_Parent_Barcode; 1514 Pid = &m_truth_Parent_Pid; 1515 M = &m_truth_Parent_M; 1516 E = &m_truth_Parent_E; 1517 Eta = &m_truth_Parent_Eta; 1518 Phi = &m_truth_Parent_Phi; 1519 Pt = &m_truth_Parent_Pt; 1520 xF = &m_truth_Parent_xF; 1521 EndVtx_x = &m_truth_Parent_EndVtx_x; 1522 EndVtx_y = &m_truth_Parent_EndVtx_y; 1523 EndVtx_z = &m_truth_Parent_EndVtx_z; 1524 /* EndVtx_t = &m_truth_Parent_EndVtx_t; */ 1525 EndVtx_r = &m_truth_Parent_EndVtx_r; 1526 } 1527 if (which == "decay1") { 1528 /* std::cout << "In FillTruthParticle: setting decay1\n"; */ 1529 Barcode = &m_truth_Decay1_Barcode; 1530 Pid = &m_truth_Decay1_Pid; 1531 M = &m_truth_Decay1_M; 1532 E = &m_truth_Decay1_E; 1533 Eta = &m_truth_Decay1_Eta; 1534 Phi = &m_truth_Decay1_Phi; 1535 Pt = &m_truth_Decay1_Pt; 1536 xF = &m_truth_Decay1_xF; 1537 EndVtx_x = &m_truth_Decay1_EndVtx_x; 1538 EndVtx_y = &m_truth_Decay1_EndVtx_y; 1539 EndVtx_z = &m_truth_Decay1_EndVtx_z; 1540 /* EndVtx_t = &m_truth_Decay1_EndVtx_t; */ 1541 EndVtx_r = &m_truth_Decay1_EndVtx_r; 1542 } 1543 if (which == "decay2") { 1544 Barcode = &m_truth_Decay2_Barcode; 1545 Pid = &m_truth_Decay2_Pid; 1546 M = &m_truth_Decay2_M; 1547 E = &m_truth_Decay2_E; 1548 Eta = &m_truth_Decay2_Eta; 1549 Phi = &m_truth_Decay2_Phi; 1550 Pt = &m_truth_Decay2_Pt; 1551 xF = &m_truth_Decay2_xF; 1552 EndVtx_x = &m_truth_Decay2_EndVtx_x; 1553 EndVtx_y = &m_truth_Decay2_EndVtx_y; 1554 EndVtx_z = &m_truth_Decay2_EndVtx_z; 1555 /* EndVtx_t = &m_truth_Decay2_EndVtx_t; */ 1556 EndVtx_r = &m_truth_Decay2_EndVtx_r; 1557 } 1558 1559 if (!par) { 1560 /* std::cout << "In FillTruthParticle: pushing back for nullptr\n"; */ 1561 Barcode->push_back(-9999); 1562 Pid->push_back(-1); 1563 M->push_back(-9999); 1564 E->push_back(-9999); 1565 Eta->push_back(-9999); 1566 Phi->push_back(-9999); 1567 Pt->push_back(-9999); 1568 xF->push_back(-9999); 1569 EndVtx_x->push_back(-9999); 1570 EndVtx_y->push_back(-9999); 1571 EndVtx_z->push_back(-9999); 1572 /* EndVtx_t->push_back(-9999); */ 1573 EndVtx_r->push_back(-9999); 1574 return; 1575 } 1576 else { 1577 if (which == "parent" && !withinAcceptance(par)) return; 1578 TLorentzVector truth_momentum; 1579 truth_momentum.SetPxPyPzE(par->get_px(), par->get_py(), par->get_pz(), par->get_e()); 1580 int this_id = par->get_track_id(); 1581 1582 /* std::cout << "In FillTruthParticle: pushing back barcode " << par->get_barcode() << ", E=" << par->get_e() << " for " << which << "\n"; */ 1583 Barcode->push_back(par->get_barcode()); 1584 Pid->push_back(par->get_pid()); 1585 M->push_back(truth_momentum.M()); 1586 E->push_back(truth_momentum.E()); 1587 Eta->push_back(truth_momentum.PseudoRapidity()); 1588 Phi->push_back(truth_momentum.Phi()); 1589 Pt->push_back(truth_momentum.Pt()); 1590 float xf = 2.0*truth_momentum.Pz()/200.0; 1591 xF->push_back(xf); 1592 1593 PHG4VtxPoint* end_vtx = getG4EndVtx(this_id); 1594 if (!end_vtx) { 1595 /* std::cout << "\nGreg info: no end vertex found for Geant-decayed particle:\n"; */ 1596 /* par->identify(); */ 1597 /* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); */ 1598 EndVtx_x->push_back(-9999); 1599 EndVtx_y->push_back(-9999); 1600 EndVtx_z->push_back(-9999); 1601 /* EndVtx_t->push_back(-9999); */ 1602 EndVtx_r->push_back(-9999); 1603 /* return; */ 1604 } 1605 else { 1606 EndVtx_x->push_back(end_vtx->get_x()); 1607 EndVtx_y->push_back(end_vtx->get_y()); 1608 EndVtx_z->push_back(end_vtx->get_z()); 1609 /* EndVtx_t->push_back(end_vtx->get_t()); */ 1610 float x = end_vtx->get_x(); 1611 float y = end_vtx->get_y(); 1612 float r = sqrt(x*x + y*y); 1613 EndVtx_r->push_back(r); 1614 } 1615 } 1616 } 1617 1618 void pythiaEMCalAna::FillAllClustersFromParent(PHG4Particle* par) { 1619 std::set<RawCluster*> cl_set = m_clusterEval->all_clusters_from(par); 1620 std::vector<float> allIDs; 1621 std::vector<float> allEfractions; 1622 if (cl_set.empty()) { 1623 std::cout << "Greg info: in FillAllClustersFromParent, cl_set is empty. Parent is\n"; 1624 par->identify(); 1625 m_truth_Parent_TotalNClusters.push_back(0); 1626 m_truth_Parent_AllClusterIDs.push_back(allIDs); 1627 m_truth_Parent_AllClusterEnergyFractions.push_back(allEfractions); 1628 return; 1629 } 1630 m_truth_Parent_TotalNClusters.push_back(cl_set.size()); 1631 for (auto& cl : cl_set) { 1632 allIDs.push_back(cl->get_id()); 1633 allEfractions.push_back(cl->get_ecore()/par->get_e()); 1634 } 1635 m_truth_Parent_AllClusterIDs.push_back(allIDs); 1636 m_truth_Parent_AllClusterEnergyFractions.push_back(allEfractions); 1637 return; 1638 } 1639 1640 void pythiaEMCalAna::FillParent(PHG4Particle* par) { 1641 FillTruthParticle("parent", par); 1642 FillAllClustersFromParent(par); 1643 /* std::vector<PHG4Particle*> daughters = GetAllDaughters(par); */ 1644 /* m_truth_Parent_TotalNDaughters.push_back(daughters.size()); */ 1645 /* std::vector<float> allPids; */ 1646 /* std::vector<float> allEfractions; */ 1647 /* for (auto daughter : daughters) { */ 1648 /* allPids.push_back(daughter->get_pid()); */ 1649 /* allEfractions.push_back(daughter->get_e()/par->get_e()); */ 1650 /* } */ 1651 /* m_truth_Parent_AllDaughterPids.push_back(allPids); */ 1652 /* m_truth_Parent_AllDaughterEnergyFractions.push_back(allEfractions); */ 1653 } 1654 1655 HepMC::GenParticle* pythiaEMCalAna::GetHepMCParent(PHG4Particle* par) { 1656 HepMC::GenParticle* p = getGenParticle(par->get_barcode()); 1657 HepMC::GenVertex* prod_vtx = p->production_vertex(); 1658 HepMC::GenVertex::particles_in_const_iterator parent; 1659 // should only have 1 parent; otherwise print an error 1660 if (prod_vtx->particles_in_size() == 1) { 1661 parent = prod_vtx->particles_in_const_begin(); 1662 assert((*parent)); 1663 return (*parent); 1664 } 1665 else { 1666 // weird photon -- they should only have 1 parent 1667 std::cout << "\nGreg info: in GetHepMCParent, found a photon with " << prod_vtx->particles_in_size() << " parent(s). Photon:\n"; 1668 p->print(); 1669 for (parent = prod_vtx->particles_in_const_begin(); parent != prod_vtx->particles_in_const_end(); parent++) { 1670 std::cout << "Parent:\n"; 1671 (*parent)->print(); 1672 } 1673 std::cout << "\n"; 1674 return nullptr; 1675 } 1676 } 1677 1678 PHG4Particle* pythiaEMCalAna::getG4Particle(int barcode) { 1679 PHG4TruthInfoContainer::Range truthRange = m_truthInfo->GetPrimaryParticleRange(); 1680 PHG4TruthInfoContainer::ConstIterator truthIter; 1681 /* std::cout << "\n\nGreg info: starting loop over Geant primary particles\n"; */ 1682 for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) 1683 { 1684 PHG4Particle *p = truthIter->second; 1685 if ( p->get_barcode() == barcode ) { 1686 // found the right particle 1687 return p; 1688 } 1689 } 1690 // reached end of loop... if we still haven't found the right particle, 1691 // we have a problem 1692 std::cout << "\t\tGreg info: in getG4Particle(), could not find correct generated particle!\n"; 1693 return nullptr; 1694 } 1695 1696 void pythiaEMCalAna::FillTruth(HepMC::GenParticle* par) { 1697 if (!par) { 1698 std::cout << "In FillTruth: par is empty!\n"; 1699 m_truth_Parent_Barcode.push_back(-9999); 1700 m_truth_Parent_Pid.push_back(-1); 1701 m_truth_Parent_M.push_back(-9999); 1702 m_truth_Parent_E.push_back(-9999); 1703 m_truth_Parent_Eta.push_back(-9999); 1704 m_truth_Parent_Phi.push_back(-9999); 1705 m_truth_Parent_Pt.push_back(-9999); 1706 m_truth_Parent_xF.push_back(-9999); 1707 m_truth_Parent_EndVtx_x.push_back(-9999); 1708 m_truth_Parent_EndVtx_y.push_back(-9999); 1709 m_truth_Parent_EndVtx_z.push_back(-9999); 1710 /* m_truth_Parent_EndVtx_t.push_back(-9999); */ 1711 m_truth_Parent_EndVtx_r.push_back(-9999); 1712 m_truth_Parent_TotalNClusters.push_back(-9999); 1713 std::vector<float> allIDs; 1714 std::vector<float> allEfractions; 1715 m_truth_Parent_AllClusterIDs.push_back(allIDs); 1716 m_truth_Parent_AllClusterEnergyFractions.push_back(allEfractions); 1717 return; 1718 } 1719 else { 1720 if (! withinAcceptance(par)) { 1721 /* std::cout << "In FillTruth: HepMC particle is outside acceptance. Details:\n"; */ 1722 /* par->print(); */ 1723 /* HepMC::FourVector mom = par->momentum(); */ 1724 /* std::cout << "(eta, E) = (" << mom.pseudoRapidity() << ", " << mom.e() << ")\n"; */ 1725 return; 1726 } 1727 primaryBarcodes.insert(par->barcode()); 1728 HepMC::FourVector mom = par->momentum(); 1729 TLorentzVector truth_momentum; 1730 truth_momentum.SetPxPyPzE(mom.px(), mom.py(), mom.pz(), mom.e()); 1731 1732 HepMC::GenVertex* end_vtx = par->end_vertex(); 1733 if (!end_vtx) { 1734 std::cout << "\nGreg info: no end vertex found for Pythia-decayed particle:\n"; 1735 par->print(); 1736 std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); 1737 return; 1738 } 1739 HepMC::FourVector position = end_vtx->position(); 1740 1741 m_truth_Parent_Barcode.push_back(par->barcode()); 1742 m_truth_Parent_Pid.push_back(par->pdg_id()); 1743 m_truth_Parent_M.push_back(truth_momentum.M()); 1744 m_truth_Parent_E.push_back(truth_momentum.E()); 1745 m_truth_Parent_Eta.push_back(truth_momentum.PseudoRapidity()); 1746 m_truth_Parent_Phi.push_back(truth_momentum.Phi()); 1747 m_truth_Parent_Pt.push_back(truth_momentum.Pt()); 1748 float xf = 2.0*truth_momentum.Pz()/200.0; 1749 m_truth_Parent_xF.push_back(xf); 1750 1751 m_truth_Parent_EndVtx_x.push_back(position.x()); 1752 m_truth_Parent_EndVtx_y.push_back(position.y()); 1753 m_truth_Parent_EndVtx_z.push_back(position.z()); 1754 /* EndVtx_t.push_back(end_vtx->get_t()); */ 1755 float x = position.x(); 1756 float y = position.y(); 1757 float r = sqrt(x*x + y*y); 1758 m_truth_Parent_EndVtx_r.push_back(r); 1759 1760 int NClusters = 0; 1761 std::vector<float> allIDs; 1762 std::vector<float> allEfractions; 1763 bool first = false; 1764 bool second = false; 1765 1766 HepMC::GenVertex::particles_out_const_iterator out; 1767 /* std::cout << "Greg info: In FillParent(HepMC), end_vtx has " << end_vtx->particles_out_size() << " outgoing particles. Adding their info (parent is)\n"; */ 1768 /* par->print(); */ 1769 for (out = end_vtx->particles_out_const_begin(); out != end_vtx->particles_out_const_end(); out++) { 1770 PHG4Particle* out_g4 = getG4Particle((*out)->barcode()); 1771 std::set<RawCluster*> cl_set = m_clusterEval->all_clusters_from(out_g4); 1772 if (cl_set.empty()) { 1773 /* std::cout << "Greg info: in FillTruth, cl_set is empty. Daughter is\n"; */ 1774 /* out_g4->identify(); */ 1775 } 1776 else { 1777 for (auto& cl : cl_set) { 1778 NClusters++; 1779 allIDs.push_back(cl->get_id()); 1780 allEfractions.push_back(cl->get_ecore()/truth_momentum.E()); 1781 } 1782 } 1783 if (!first) { 1784 first = true; 1785 /* std::cout << "Filling decay1\n"; */ 1786 FillDecay("decay1", out_g4, nullptr); 1787 } 1788 else { 1789 if (!second) { 1790 second = true; 1791 FillDecay("decay2", out_g4, nullptr); 1792 } 1793 } 1794 } // end loop over outgoing particles 1795 m_truth_Parent_TotalNClusters.push_back(NClusters); 1796 m_truth_Parent_AllClusterIDs.push_back(allIDs); 1797 m_truth_Parent_AllClusterEnergyFractions.push_back(allEfractions); 1798 } 1799 } 1800 1801 std::vector<PHG4Particle*> pythiaEMCalAna::GetAllDaughters(PHG4Particle* parent) { 1802 std::vector<PHG4Particle*> vec; 1803 int parent_id = parent->get_track_id(); 1804 /* std::cout << "In GetAllDaughters: parent is\n"; */ 1805 /* parent->identify(); */ 1806 PHG4TruthInfoContainer::Range truthRange = m_truthInfo->GetSecondaryParticleRange(); 1807 PHG4TruthInfoContainer::ConstIterator truthIter; 1808 for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) 1809 { 1810 PHG4Particle *decayPar = truthIter->second; 1811 int this_parent_id = decayPar->get_parent_id(); 1812 if (this_parent_id == parent_id) { 1813 vec.push_back(decayPar); 1814 } 1815 } 1816 return vec; 1817 } 1818 1819 void pythiaEMCalAna::PrintParent(PHG4Particle* par) { 1820 std::vector<PHG4Particle*> daughters = GetAllDaughters(par); 1821 std::set<RawCluster*> cl_set = m_clusterEval->all_clusters_from(par); 1822 std::cout << "In PrintParent: parent is \n"; 1823 par->identify(); 1824 std::cout << Form("Parent has %ld daughter particles and %ld associated clusters.\n", daughters.size(), cl_set.size()); 1825 return; 1826 } 1827 1828 void pythiaEMCalAna::GetBestDaughters(PHG4Particle* parent, PHG4Particle* &decay1, PHG4Particle* &decay2) { 1829 /* std::cout << "\nIn GetBestDaughters: parent is\n"; */ 1830 /* parent->identify(); */ 1831 std::vector<PHG4Particle*> daughters = GetAllDaughters(parent); 1832 /* std::cout << "Parent has " << daughters.size() << " daughters.\n"; */ 1833 if (daughters.size() == 0) { 1834 decay1 = nullptr; 1835 decay2 = nullptr; 1836 return; 1837 } 1838 if (daughters.size() == 1) { 1839 decay1 = daughters.at(0); 1840 /* std::cout << "In GetBestDaughters: only daughter is" << decay1 << ":\n"; */ 1841 /* decay1->identify(); */ 1842 decay2 = nullptr; 1843 return; 1844 } 1845 std::vector<float> Efractions; 1846 /* for (std::vector<PHG4Particle*>::iterator it = daughters.begin(); it != daughters.end(); it++) { */ 1847 for (auto daughter : daughters) { 1848 Efractions.push_back(daughter->get_e()/parent->get_e()); 1849 } 1850 float maxE = 0.0; 1851 float submaxE = 0.0; 1852 std::vector<float>::iterator max; 1853 max = std::max_element(Efractions.begin(), Efractions.end()); 1854 maxE = *max; 1855 unsigned int max_index = std::distance(Efractions.begin(), max); 1856 decay1 = daughters.at(max_index); 1857 if (false) std::cout << "In GetBestDaughters: highest energy daughter is " << decay1 << " with energy fraction " << maxE << ":\n"; 1858 /* decay1->identify(); */ 1859 // set the max to 0 so we can use max_element to find the submax 1860 Efractions.at(max_index) = 0.0; 1861 max = std::max_element(Efractions.begin(), Efractions.end()); 1862 submaxE = *max; 1863 unsigned int submax_index = std::distance(Efractions.begin(), max); 1864 decay2 = daughters.at(submax_index); 1865 if (false) std::cout << "In GetBestDaughters: second highest energy daughter is " << decay2 << " with energy fraction " << submaxE << ":\n"; 1866 /* decay2->identify(); */ 1867 } 1868 1869 /* void pythiaEMCalAna::GetBestDaughters(PHG4Particle* parent, PHG4Particle* &decay1, PHG4Particle* &decay2) { */ 1870 /* std::cout << "In GetBestDaughters: parent is\n"; */ 1871 /* parent->identify(); */ 1872 /* bool first = false; */ 1873 /* bool second = false; */ 1874 /* PHG4TruthInfoContainer::Range truthRange = m_truthInfo->GetSecondaryParticleRange(); */ 1875 /* PHG4TruthInfoContainer::ConstIterator truthIter; */ 1876 /* for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) */ 1877 /* { */ 1878 /* PHG4Particle *decayPar = truthIter->second; */ 1879 /* int this_parent_id = decayPar->get_parent_id(); */ 1880 /* if (this_parent_id == parent_id) { */ 1881 /* if (!first) { */ 1882 /* std::cout << "In GetDaughters: found first. decayPar=" << decayPar << "\n"; */ 1883 /* decay1 = decayPar; */ 1884 /* first = true; */ 1885 /* std::cout << "decay1=" << decay1 << "\n"; */ 1886 /* std::cout << "Found decay1 with primary id " << decay1->get_primary_id() << ": "; */ 1887 /* decay1->identify(); */ 1888 /* } */ 1889 /* else { */ 1890 /* if (!second) { */ 1891 /* std::cout << "In GetDaughters: found second. decayPar=" << decayPar << "\n"; */ 1892 /* decay2 = decayPar; */ 1893 /* second = true; */ 1894 /* std::cout << "decay2=" << decay2 << "\n"; */ 1895 /* std::cout << "Found decay2 with primary id " << decay2->get_primary_id() << ": "; */ 1896 /* decay2->identify(); */ 1897 /* /1* break; *1/ */ 1898 /* } */ 1899 /* else { */ 1900 /* std::cout << "Greg info: found a parent with more than 2 decay particles. Parent is \n"; */ 1901 /* parent->identify(); */ 1902 /* std::cout << "Decay1 is \n"; */ 1903 /* decay1->identify(); */ 1904 /* std::cout << "Decay2 is \n"; */ 1905 /* decay2->identify(); */ 1906 /* std::cout << "decayPar is \n"; */ 1907 /* decayPar->identify(); */ 1908 /* break; */ 1909 /* } */ 1910 /* } */ 1911 /* } */ 1912 /* } */ 1913 /* if (!first || !second) { */ 1914 /* std::cout << "Greg info: got to end of secondary particles range without finding two decays!\n"; */ 1915 /* if (!first) std::cout << "\tFound 0 decays. Parent is\n"; */ 1916 /* else std::cout << "\tFound 1 decay. Parent is\n"; */ 1917 /* parent->identify(); */ 1918 /* } */ 1919 /* } */ 1920 1921 void pythiaEMCalAna::FillDecay(std::string which, PHG4Particle* decay, PHG4Particle* parent) { 1922 /* std::cout << "Entering FillDecay(" << which << ")\n"; */ 1923 FillTruthParticle(which, decay); 1924 std::vector<float>* TotalNClusters = nullptr; 1925 std::vector<float>* BestClusterID = nullptr; 1926 std::vector<float>* BestClusterEfraction = nullptr; 1927 /* std::vector<std::vector<float>> AllClusterIDs; */ 1928 /* std::vector<std::vector<float>> AllClusterEnergyFractions; */ 1929 if (which == "decay1") { 1930 TotalNClusters = &m_truth_Decay1_TotalNClusters; 1931 BestClusterID = &m_truth_Decay1_BestClusterID; 1932 BestClusterEfraction = &m_truth_Decay1_BestClusterEfraction; 1933 /* AllClusterIDs = m_truth_Decay1_AllClusterIDs; */ 1934 /* AllClusterEnergyFractions = m_truth_Decay1_AllClusterEnergyFractions; */ 1935 } 1936 if (which == "decay2") { 1937 TotalNClusters = &m_truth_Decay2_TotalNClusters; 1938 BestClusterID = &m_truth_Decay2_BestClusterID; 1939 BestClusterEfraction = &m_truth_Decay2_BestClusterEfraction; 1940 /* AllClusterIDs = m_truth_Decay2_AllClusterIDs; */ 1941 /* AllClusterEnergyFractions = m_truth_Decay2_AllClusterEnergyFractions; */ 1942 } 1943 if (!decay) { 1944 TotalNClusters->push_back(0); 1945 BestClusterID->push_back(-9999); 1946 BestClusterEfraction->push_back(-9999); 1947 /* std::vector<float> allIDs; */ 1948 /* std::vector<float> allEfractions; */ 1949 /* AllClusterIDs->push_back(allIDs); */ 1950 /* AllClusterEnergyFractions->push_back(allEfractions); */ 1951 return; 1952 } 1953 1954 std::set<RawCluster*> cl_set = m_clusterEval->all_clusters_from(decay); 1955 if (cl_set.empty()) { 1956 /* std::cout << "Greg info: in FillDecay, cl_set is empty. Decay is " << decay << ":\n"; */ 1957 /* decay->identify(); */ 1958 TotalNClusters->push_back(-9999); 1959 } 1960 else { 1961 TotalNClusters->push_back(cl_set.size()); 1962 } 1963 RawCluster* best_cl = m_clusterEval->best_cluster_from(decay); 1964 if (!best_cl) { 1965 /* std::cout << "Greg info: in FillDecay, best_cluster_from failed\n"; */ 1966 /* std::cout << "decay is "; */ 1967 /* decay->identify(); */ 1968 if (parent) best_cl = FindBestCluster(decay, parent); 1969 else { 1970 /* std::cout << "No parent given, and best_cluster_from failed!\n"; */ 1971 BestClusterID->push_back(-9999); 1972 BestClusterEfraction->push_back(-9999); 1973 return; 1974 } 1975 /* std::cout << "Result of FindBestCluster is "; */ 1976 /* CLHEP::Hep3Vector cl3vec = best_cl->get_position(); */ 1977 /* std::cout << Form("id=%d, E=%f, eta=%f, phi=%f\n", best_cl->get_id(), best_cl->get_energy(), cl3vec.pseudoRapidity(), cl3vec.phi()); */ 1978 1979 } 1980 else { 1981 /* std::cout << "Greg info: in FillDecay, best_cluster_from succeeded\n"; */ 1982 } 1983 BestClusterID->push_back(best_cl->get_id()); 1984 BestClusterEfraction->push_back(best_cl->get_ecore()/decay->get_e()); 1985 1986 /* std::vector<float> allIDs; */ 1987 /* std::vector<float> allEfractions; */ 1988 /* if (cl_set.empty()) { */ 1989 /* allIDs->push_back(best_cl->get_id()); */ 1990 /* allEfractions->push_back(best_cl->get_ecore()/decay->get_e()); */ 1991 /* } */ 1992 /* else { */ 1993 /* for (auto& cl : cl_set) { */ 1994 /* allIDs->push_back(cl->get_id()); */ 1995 /* allEfractions->push_back(cl->get_ecore()/decay->get_e()); */ 1996 /* } */ 1997 /* } */ 1998 /* AllClusterIDs->push_back(allIDs); */ 1999 /* AllClusterEnergyFractions->push_back(allEfractions); */ 2000 } 2001 2002 RawCluster* pythiaEMCalAna::FindBestCluster(PHG4Particle* decay, PHG4Particle* parent) { 2003 RawCluster* best_cl = nullptr; 2004 std::set<RawCluster*> cl_set = m_clusterEval->all_clusters_from(parent); 2005 if (cl_set.empty()) { 2006 std::cout << "Greg info: in FindBestCluster, cl_set is empty. Parent is\n"; 2007 parent->identify(); 2008 return best_cl; 2009 } 2010 /* std::cout << "In FindBestCluster: parent has " << cl_set.size() << " associated clusters.\n"; */ 2011 TLorentzVector decay_vec; 2012 decay_vec.SetPxPyPzE(decay->get_px(), decay->get_py(), decay->get_pz(), decay->get_e()); 2013 TLorentzVector cl_vec; 2014 float min_dR = 999.9; 2015 for (auto& cl : cl_set) { 2016 CLHEP::Hep3Vector cl3vec = cl->get_position(); 2017 cl_vec.SetPtEtaPhiM(cl3vec.perp(), cl3vec.pseudoRapidity(), cl3vec.phi(), 0); 2018 float dR = decay_vec.DeltaR(cl_vec); 2019 min_dR = std::min(min_dR, dR); 2020 if (dR == min_dR) best_cl = cl; 2021 /* std::cout << "In FindBestCluster: min_dR=" << min_dR << " from cluster id " << best_cl->get_id() << "\n"; */ 2022 } 2023 /* std::cout << "Done in FindBestCluster. min_dR=" << min_dR << "\n"; */ 2024 return best_cl; 2025 } 2026 2027 HepMC::GenParticle* pythiaEMCalAna::getGenParticle(int barcode) { 2028 /* std::cout << "Greg info: in getGenParticle, looking for barcode " << barcode << "\n"; */ 2029 /* std::cout << "m_theEvent = " << m_theEvent << "\n"; */ 2030 for(HepMC::GenEvent::particle_const_iterator p=m_theEvent->particles_begin(); p!=m_theEvent->particles_end(); ++p) 2031 { 2032 /* std::cout << "In loop, (*p)=" << (*p) << "\n"; */ 2033 assert(*p); 2034 /* std::cout << "\tpythia barcode " << (*p)->barcode() << "\n"; */ 2035 /* if ((*p)->barcode() == 33) { */ 2036 /* std::cout << "barcode 33 details: \t"; */ 2037 /* (*p)->print(); */ 2038 /* } */ 2039 if ( (*p)->barcode() == barcode ) { 2040 // found the right particle 2041 /* std::cout << "\tpythia barcode " << (*p)->barcode() << ", returning (*p)\n"; */ 2042 return (*p); 2043 } 2044 } 2045 // reached end of loop... if we still haven't found the right particle, 2046 // we have a problem 2047 std::cout << "\t\tGreg info: in getGenParticle(), could not find correct generated particle!\n"; 2048 return nullptr; 2049 } 2050 2051 bool pythiaEMCalAna::withinAcceptance(PHG4Particle* part) { 2052 float max_eta = 1.1; 2053 float min_E = 1.0; 2054 TLorentzVector truth_momentum; 2055 truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); 2056 float eta = truth_momentum.PseudoRapidity(); 2057 if (abs(eta) > max_eta) return false; 2058 if (truth_momentum.E() < min_E) return false; 2059 else return true; 2060 } 2061 2062 bool pythiaEMCalAna::withinAcceptance(HepMC::GenParticle* part) { 2063 float max_eta = 1.1; 2064 float min_E = 1.0; 2065 HepMC::FourVector mom = part->momentum(); 2066 TLorentzVector truth_momentum; 2067 truth_momentum.SetPxPyPzE(mom.px(), mom.py(), mom.pz(), mom.e()); 2068 float eta = truth_momentum.PseudoRapidity(); 2069 if (abs(eta) > max_eta) return false; 2070 if (truth_momentum.E() < min_E) return false; 2071 else return true; 2072 } 2073 2074 PHG4VtxPoint* pythiaEMCalAna::getG4EndVtx(int id) { 2075 /* std::cout << "Entering getG4EndVtx\n"; */ 2076 PHG4VtxPoint* end_vtx = nullptr; 2077 PHG4TruthInfoContainer::Range truthRange = m_truthInfo->GetSecondaryParticleRange(); 2078 PHG4TruthInfoContainer::ConstIterator truthIter; 2079 for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) 2080 { 2081 PHG4Particle *decayPar = truthIter->second; 2082 /* std::cout << "In getG4EndVtx: decayPar = " << decayPar << "\n"; */ 2083 int parent_id = decayPar->get_parent_id(); 2084 /* std::cout << "Parent id is " << parent_id << "\n"; */ 2085 if (parent_id == id) { 2086 end_vtx = m_truthInfo->GetVtx(decayPar->get_vtx_id()); 2087 /* std::cout << "In getG4EndVtx: assigning end_vtx\n"; */ 2088 break; 2089 } 2090 } 2091 if (!end_vtx) { 2092 // couldn't find any daughter particles in secondary list 2093 // check primary list instead?? 2094 /* std::cout << "Could not find any daughters in secondary particle list\n"; */ 2095 truthRange = m_truthInfo->GetPrimaryParticleRange(); 2096 for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++) 2097 { 2098 PHG4Particle *decayPar = truthIter->second; 2099 int parent_id = decayPar->get_parent_id(); 2100 if (parent_id == id) { 2101 std::cout << "Greg info: in getG4EndVtx, found daughter in *primary* particle range!\n"; 2102 end_vtx = m_truthInfo->GetVtx(decayPar->get_vtx_id()); 2103 break; 2104 } 2105 } 2106 } 2107 /* std::cout << "In getG4EndVtx: returning end_vtx = " << end_vtx << "\n"; */ 2108 return end_vtx; 2109 } 2110 2111 // Old stuff 2112 2113 /* void pythiaEMCalAna::addDirectPhoton(PHG4Particle* part, PHG4TruthInfoContainer* truthInfo) { */ 2114 /* /1* std::cout << "Found a direct photon!\n"; *1/ */ 2115 /* TLorentzVector truth_momentum; */ 2116 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 2117 /* if (! withinAcceptance(part)) return; */ 2118 /* /1* std::cout << "\tIs within acceptance!\n"; *1/ */ 2119 /* /1* part->identify(); *1/ */ 2120 /* int this_id = part->get_track_id(); */ 2121 2122 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); */ 2123 /* if (!end_vtx) { */ 2124 /* std::cout << "\nGreg info: no end vertex found for direct photon:\n"; */ 2125 /* part->identify(); */ 2126 /* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); */ 2127 /* return; */ 2128 /* } */ 2129 2130 /* int embedID = truthInfo->isEmbeded(part->get_track_id()); */ 2131 /* m_truthBarcode.push_back(std::pair<int,int>(embedID, part->get_barcode())); */ 2132 /* m_truthParentBarcode.push_back(std::pair<int,int>(embedID, 0)); */ 2133 /* m_truthIsPrimary.push_back(1); */ 2134 /* m_truthPid.push_back(part->get_pid()); */ 2135 2136 /* /1* TLorentzVector truth_momentum; *1/ */ 2137 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2138 /* m_truthE.push_back(truth_momentum.E()); */ 2139 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2140 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2141 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2142 /* m_truthMass.push_back(truth_momentum.M()); */ 2143 2144 /* // part->get_vtx_id() gives the *production* vertex, not the end vertex */ 2145 /* /1* PHG4VtxPoint* end_vtx = truthInfo->GetVtx(part->get_vtx_id()); // DOES NOT WORK *1/ */ 2146 /* // get the end vertex from one of the daughter particles */ 2147 /* /1* int this_id = part->get_track_id(); *1/ */ 2148 /* /1* std::cout << "this_id = " << this_id << "; getting end vertex\n"; *1/ */ 2149 /* /1* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); *1/ */ 2150 /* /1* std::cout << "end_vtx = " << end_vtx << "\n"; *1/ */ 2151 /* assert(end_vtx); */ 2152 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 2153 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 2154 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 2155 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 2156 /* float x = end_vtx->get_x(); */ 2157 /* float y = end_vtx->get_y(); */ 2158 /* float r = sqrt(x*x + y*y); */ 2159 /* m_truthEndVtx_r.push_back(r); */ 2160 2161 /* /1* primaryBarcodes.push_back(part->get_barcode()); *1/ */ 2162 /* } */ 2163 2164 /* void pythiaEMCalAna::addDecayPhoton(PHG4Particle* part, PHG4TruthInfoContainer* truthInfo, HepMC::GenEvent* theEvent) { */ 2165 /* TLorentzVector truth_momentum; */ 2166 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 2167 /* if (! withinAcceptance(part)) return; */ 2168 /* int this_id = part->get_track_id(); */ 2169 2170 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); */ 2171 /* if (!end_vtx) { */ 2172 /* /1* std::cout << "\nGreg info: no end vertex found for decay photon:\n"; *1/ */ 2173 /* /1* part->identify(); *1/ */ 2174 /* /1* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); *1/ */ 2175 /* return; */ 2176 /* } */ 2177 2178 /* int parent_barcode; */ 2179 /* // If Geant handled the decay, we can get the parent directly */ 2180 /* PHG4Particle* geant_parent = truthInfo->GetParticle(part->get_parent_id()); */ 2181 /* if (geant_parent) { */ 2182 /* parent_barcode = geant_parent->get_barcode(); */ 2183 /* } */ 2184 /* else { */ 2185 /* // pythia handled the decay, so get the parent from there */ 2186 /* HepMC::GenParticle* pho = getGenParticle(part->get_barcode(), theEvent); */ 2187 /* assert(pho); */ 2188 /* HepMC::GenVertex* prod_vtx = pho->production_vertex(); */ 2189 /* if (prod_vtx->particles_in_size() == 1) { */ 2190 /* HepMC::GenVertex::particles_in_const_iterator pythia_parent = prod_vtx->particles_in_const_begin(); */ 2191 /* assert((*pythia_parent)); */ 2192 /* parent_barcode = (*pythia_parent)->barcode(); */ 2193 /* } */ 2194 /* else { */ 2195 /* std::cout << "Greg info: pythia-decayed photon with " << prod_vtx->particles_in_size() << " parents. Skipping...\n"; */ 2196 /* return; */ 2197 /* } */ 2198 /* } // done getting parent barcode */ 2199 2200 /* /1* std::cout << "Adding decay photon with barcode " << part->get_barcode() << "; geant_parent is " << geant_parent << "\n"; *1/ */ 2201 /* int embedID = truthInfo->isEmbeded(part->get_track_id()); */ 2202 /* m_truthBarcode.push_back(std::pair<int,int>(embedID, part->get_barcode())); */ 2203 /* m_truthParentBarcode.push_back(std::pair<int,int>(embedID, parent_barcode)); */ 2204 /* m_truthIsPrimary.push_back(0); */ 2205 /* m_truthPid.push_back(part->get_pid()); */ 2206 2207 /* /1* TLorentzVector truth_momentum; *1/ */ 2208 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2209 /* m_truthE.push_back(truth_momentum.E()); */ 2210 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2211 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2212 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2213 /* m_truthMass.push_back(truth_momentum.M()); */ 2214 2215 /* // part->get_vtx_id() gives the *production* vertex, not the end vertex */ 2216 /* /1* PHG4VtxPoint* end_vtx = truthInfo->GetVtx(part->get_vtx_id()); // DOES NOT WORK *1/ */ 2217 /* // get the end vertex from one of the daughter particles */ 2218 /* /1* int this_id = part->get_track_id(); *1/ */ 2219 /* /1* std::cout << "this_id = " << this_id << "; getting end vertex\n"; *1/ */ 2220 /* /1* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); *1/ */ 2221 /* /1* std::cout << "end_vtx = " << end_vtx << "\n"; *1/ */ 2222 /* assert(end_vtx); */ 2223 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 2224 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 2225 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 2226 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 2227 /* float x = end_vtx->get_x(); */ 2228 /* float y = end_vtx->get_y(); */ 2229 /* float r = sqrt(x*x + y*y); */ 2230 /* m_truthEndVtx_r.push_back(r); */ 2231 2232 /* /1* secondaryBarcodes.push_back(part->get_barcode()); *1/ */ 2233 /* } */ 2234 2235 /* void pythiaEMCalAna::addPrimaryHadronFromPythia(HepMC::GenParticle* part, int embedID) { */ 2236 /* // case where PYTHIA handles the decay and Geant never knows about the primary */ 2237 /* HepMC::FourVector mom = part->momentum(); */ 2238 /* TLorentzVector truth_momentum; */ 2239 /* truth_momentum.SetPxPyPzE(mom.px(), mom.py(), mom.pz(), mom.e()); */ 2240 /* if (! withinAcceptance(part)) return; */ 2241 /* HepMC::GenVertex* end_vtx = part->end_vertex(); */ 2242 /* if (!end_vtx) { */ 2243 /* std::cout << "\nGreg info: no end vertex found for Pythia-decayed primary:\n"; */ 2244 /* part->print(); */ 2245 /* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); */ 2246 /* return; */ 2247 /* } */ 2248 2249 /* m_truthBarcode.push_back(std::pair<int,int>(embedID, part->barcode())); */ 2250 /* m_truthParentBarcode.push_back(std::pair<int,int>(embedID, 0)); */ 2251 /* m_truthIsPrimary.push_back(1); */ 2252 /* m_truthPid.push_back(part->pdg_id()); */ 2253 /* if (part->pdg_id() == 22) std::cout << "Greg info: primary in addPrimaryHadronFromPythia is a photon!\n"; */ 2254 2255 /* /1* TLorentzVector truth_momentum; *1/ */ 2256 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2257 /* m_truthE.push_back(truth_momentum.E()); */ 2258 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2259 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2260 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2261 /* m_truthMass.push_back(truth_momentum.M()); */ 2262 2263 /* assert(end_vtx); */ 2264 /* HepMC::FourVector position = end_vtx->position(); */ 2265 2266 /* m_truthEndVtx_x.push_back(position.x()); */ 2267 /* m_truthEndVtx_y.push_back(position.y()); */ 2268 /* m_truthEndVtx_z.push_back(position.z()); */ 2269 /* m_truthEndVtx_t.push_back(position.t()); */ 2270 /* float x = position.x(); */ 2271 /* float y = position.y(); */ 2272 /* float r = sqrt(x*x + y*y); */ 2273 /* /1* std::cout << Form("pythiaEMCalAna::addPrimaryFromPythia(): Vertex r=%f, perp=%f\n", r, position.perp()); *1/ */ 2274 /* m_truthEndVtx_r.push_back(r); */ 2275 2276 /* /1* primaryBarcodes.push_back(part->barcode()); *1/ */ 2277 /* } */ 2278 2279 /* void pythiaEMCalAna::addPrimaryHadronFromGeant(PHG4Particle* part, PHG4TruthInfoContainer* truthInfo) { */ 2280 /* // case where Geant handles the decay */ 2281 /* TLorentzVector truth_momentum; */ 2282 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 2283 /* if (! withinAcceptance(part)) return; */ 2284 /* int this_id = part->get_track_id(); */ 2285 2286 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); */ 2287 /* if (!end_vtx) { */ 2288 /* std::cout << "\nGreg info: no end vertex found for Geant-decayed primary:\n"; */ 2289 /* part->identify(); */ 2290 /* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); */ 2291 /* return; */ 2292 /* } */ 2293 2294 /* int embedID = truthInfo->isEmbeded(part->get_track_id()); */ 2295 /* m_truthBarcode.push_back(std::pair<int,int>(embedID, part->get_barcode())); */ 2296 /* m_truthParentBarcode.push_back(std::pair<int,int>(embedID, 0)); */ 2297 /* m_truthIsPrimary.push_back(1); */ 2298 /* m_truthPid.push_back(part->get_pid()); */ 2299 /* if (part->get_pid() == 22) std::cout << "Greg info: primary in addPrimaryHadronFromGeant is a photon!\n"; */ 2300 2301 /* /1* TLorentzVector truth_momentum; *1/ */ 2302 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2303 /* m_truthE.push_back(truth_momentum.E()); */ 2304 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2305 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2306 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2307 /* m_truthMass.push_back(truth_momentum.M()); */ 2308 2309 /* assert(end_vtx); */ 2310 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 2311 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 2312 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 2313 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 2314 /* float x = end_vtx->get_x(); */ 2315 /* float y = end_vtx->get_y(); */ 2316 /* float r = sqrt(x*x + y*y); */ 2317 /* m_truthEndVtx_r.push_back(r); */ 2318 2319 /* /1* primaryBarcodes.push_back(part->get_barcode()); *1/ */ 2320 /* } */ 2321 2322 /* void pythiaEMCalAna::addSecondaryFromPythia(HepMC::GenParticle* part) { */ 2323 /* /1* std::cout << "Entering pythiaEMCalAna::addSecondaryFromPythia\npart=" << part << "\n"; *1/ */ 2324 /* HepMC::FourVector mom = part->momentum(); */ 2325 /* TLorentzVector truth_momentum; */ 2326 /* truth_momentum.SetPxPyPzE(mom.px(), mom.py(), mom.pz(), mom.e()); */ 2327 /* float eta = truth_momentum.PseudoRapidity(); */ 2328 /* if (abs(eta) > 1.1) return; */ 2329 2330 /* HepMC::GenVertex* prod_vtx = (*p)->production_vertex(); */ 2331 /* HepMC::GenVertex::particles_in_const_iterator parent; */ 2332 /* if (prod_vtx->particles_in_size() == 1) { */ 2333 /* parent = prod_vtx->particles_in_const_begin(); */ 2334 /* } */ 2335 /* else { */ 2336 /* std::cout << "Greg info: found a weird decay photon in pythia. Has " << prod_vtx->particles_in_size() << " parents. Skipping this photon!\n"; */ 2337 /* return; */ 2338 /* } */ 2339 2340 /* m_truthBarcode.push_back(part->barcode()); */ 2341 /* m_truthParentBarcode.push_back((*parent)->barcode()); */ 2342 /* m_truthIsPrimary.push_back(0); */ 2343 /* m_truthPid.push_back(part->pdg_id()); */ 2344 2345 /* /1* std::cout << "Getting momentum\n"; *1/ */ 2346 /* /1* HepMC::FourVector mom = part->momentum(); *1/ */ 2347 /* /1* TLorentzVector truth_momentum; *1/ */ 2348 /* /1* truth_momentum.SetPxPyPzE(mom.px(), mom.py(), mom.pz(), mom.e()); *1/ */ 2349 /* m_truthE.push_back(truth_momentum.E()); */ 2350 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2351 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2352 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2353 /* m_truthMass.push_back(truth_momentum.M()); */ 2354 2355 /* /1* std::cout << "Getting end vertex\n"; *1/ */ 2356 /* HepMC::GenVertex* end_vtx = part->end_vertex(); */ 2357 /* assert(end_vtx); */ 2358 /* /1* std::cout << "Found end vertex in pythia\n"; *1/ */ 2359 /* HepMC::FourVector position = end_vtx->position(); */ 2360 2361 /* m_truthEndVtx_x.push_back(position.x()); */ 2362 /* m_truthEndVtx_y.push_back(position.y()); */ 2363 /* m_truthEndVtx_z.push_back(position.z()); */ 2364 /* m_truthEndVtx_t.push_back(position.t()); */ 2365 /* float x = position.x(); */ 2366 /* float y = position.y(); */ 2367 /* float r = sqrt(x*x + y*y); */ 2368 /* /1* std::cout << Form("pythiaEMCalAna::addPrimaryFromPythia(): Vertex r=%f, perp=%f\n", r, position.perp()); *1/ */ 2369 /* m_truthEndVtx_r.push_back(r); */ 2370 2371 /* primaryBarcodes.push_back(part->barcode()); */ 2372 /* } */ 2373 2374 /* void pythiaEMCalAna::addPrimaryFromGeant(PHG4Particle* part, PHG4TruthInfoContainer* truthInfo) { */ 2375 /* /1* std::cout << "Entering pythiaEMCalAna::addPrimaryFromGeant\npart=" << part << "\n"; *1/ */ 2376 /* /1* part->identify(); *1/ */ 2377 /* /1* std::cout << "Barcode is " << part->get_barcode() << "\n"; *1/ */ 2378 /* /1* int vtx_id = part->get_vtx_id(); *1/ */ 2379 /* /1* std::cout << "vertex id is " << vtx_id << "\n"; *1/ */ 2380 /* /1* PHG4VtxPoint* vtxpt = truthInfo->GetVtx(vtx_id); *1/ */ 2381 /* /1* vtxpt->identify(); *1/ */ 2382 /* TLorentzVector truth_momentum; */ 2383 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 2384 /* float eta = truth_momentum.PseudoRapidity(); */ 2385 /* if (abs(eta) > 1.1) return; */ 2386 /* int this_id = part->get_track_id(); */ 2387 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); */ 2388 /* if (!end_vtx) { */ 2389 /* /1* std::cout << "\nGreg info: no end vertex found for primary particle:\n"; *1/ */ 2390 /* /1* part->identify(); *1/ */ 2391 /* /1* std::cout << Form("Mass=%f, pT=%f\n\n", truth_momentum.M(), truth_momentum.Pt()); *1/ */ 2392 /* return; */ 2393 /* } */ 2394 2395 /* m_truthBarcode.push_back(part->get_barcode()); */ 2396 /* m_truthParentBarcode.push_back(0); */ 2397 /* m_truthIsPrimary.push_back(1); */ 2398 /* m_truthPid.push_back(part->get_pid()); */ 2399 2400 /* /1* TLorentzVector truth_momentum; *1/ */ 2401 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2402 /* m_truthE.push_back(truth_momentum.E()); */ 2403 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2404 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2405 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2406 /* m_truthMass.push_back(truth_momentum.M()); */ 2407 2408 /* // part->get_vtx_id() gives the *production* vertex, not the end vertex */ 2409 /* /1* PHG4VtxPoint* end_vtx = truthInfo->GetVtx(part->get_vtx_id()); // DOES NOT WORK *1/ */ 2410 /* // get the end vertex from one of the daughter particles */ 2411 /* /1* int this_id = part->get_track_id(); *1/ */ 2412 /* /1* std::cout << "this_id = " << this_id << "; getting end vertex\n"; *1/ */ 2413 /* /1* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); *1/ */ 2414 /* /1* std::cout << "end_vtx = " << end_vtx << "\n"; *1/ */ 2415 /* assert(end_vtx); */ 2416 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 2417 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 2418 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 2419 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 2420 /* float x = end_vtx->get_x(); */ 2421 /* float y = end_vtx->get_y(); */ 2422 /* float r = sqrt(x*x + y*y); */ 2423 /* m_truthEndVtx_r.push_back(r); */ 2424 2425 /* primaryBarcodes.push_back(part->get_barcode()); */ 2426 /* } */ 2427 2428 /* void pythiaEMCalAna::addSecondaryFromGeant(PHG4Particle* part, PHG4TruthInfoContainer* truthInfo) { */ 2429 /* /1* std::cout << "Entering pythiaEMCalAna::addSecondary(), geant particle is " << part << "\n"; *1/ */ 2430 /* TLorentzVector truth_momentum; */ 2431 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 2432 /* float eta = truth_momentum.PseudoRapidity(); */ 2433 /* if (abs(eta) > 1.1) return; */ 2434 /* /1* std::cout << "Done checking eta\n"; *1/ */ 2435 /* int this_id = part->get_track_id(); */ 2436 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); */ 2437 /* if (!end_vtx) { */ 2438 /* /1* std::cout << "\n\nGreg info: no end vertex found for secondary particle:\n"; *1/ */ 2439 /* /1* part->identify(); *1/ */ 2440 /* /1* std::cout << "Particle eta=" << eta << "\n"; *1/ */ 2441 /* return; */ 2442 /* } */ 2443 /* /1* std::cout << "Done checking end_vtx\n"; *1/ */ 2444 2445 /* m_truthBarcode.push_back(part->get_barcode()); */ 2446 /* PHG4Particle* parent = truthInfo->GetParticle(part->get_parent_id()); */ 2447 /* m_truthParentBarcode.push_back(parent->get_barcode()); */ 2448 /* m_truthIsPrimary.push_back(0); */ 2449 /* m_truthPid.push_back(part->get_pid()); */ 2450 /* /1* std::cout << "Done adding pid to vector\n"; *1/ */ 2451 2452 /* /1* TLorentzVector truth_momentum; *1/ */ 2453 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2454 /* m_truthE.push_back(truth_momentum.E()); */ 2455 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2456 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2457 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2458 /* m_truthMass.push_back(truth_momentum.M()); */ 2459 /* /1* std::cout << "Done adding kinematics to vector\n"; *1/ */ 2460 2461 /* /1* std::cout << "Asserting end_vtx= " << end_vtx << "\n"; *1/ */ 2462 /* assert(end_vtx); */ 2463 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 2464 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 2465 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 2466 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 2467 /* float x = end_vtx->get_x(); */ 2468 /* float y = end_vtx->get_y(); */ 2469 /* float r = sqrt(x*x + y*y); */ 2470 /* m_truthEndVtx_r.push_back(r); */ 2471 /* /1* std::cout << "Done adding vertex to vector\n"; *1/ */ 2472 2473 /* secondaryBarcodes.push_back(part->get_barcode()); */ 2474 /* } */ 2475 2476 /* void pythiaEMCalAna::addSecondaryWithoutParent(PHG4Particle* part, PHG4TruthInfoContainer* truthInfo, HepMC::GenParticle* genParent) { */ 2477 /* /1* std::cout << "Entering pythiaEMCalAna::addSecondaryWithoutParent(), geant particle is " << part << ", pythia particle is " << genParent << "\n"; *1/ */ 2478 /* TLorentzVector truth_momentum; */ 2479 /* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); */ 2480 /* float eta = truth_momentum.PseudoRapidity(); */ 2481 /* if (abs(eta) > 1.1) return; */ 2482 /* /1* std::cout << "Done checking eta\n"; *1/ */ 2483 /* int this_id = part->get_track_id(); */ 2484 /* PHG4VtxPoint* end_vtx = getG4EndVtx(this_id, truthInfo); */ 2485 /* if (!end_vtx) { */ 2486 /* /1* std::cout << "\n\nGreg info: no end vertex found for secondary particle:\n"; *1/ */ 2487 /* /1* part->identify(); *1/ */ 2488 /* /1* std::cout << "Particle eta=" << eta << "\n"; *1/ */ 2489 /* return; */ 2490 /* } */ 2491 /* /1* std::cout << "Done checking end_vtx\n"; *1/ */ 2492 2493 /* m_truthBarcode.push_back(part->get_barcode()); */ 2494 /* m_truthParentBarcode.push_back(genParent->barcode()); */ 2495 /* m_truthIsPrimary.push_back(0); */ 2496 /* m_truthPid.push_back(part->get_pid()); */ 2497 /* /1* std::cout << "Done adding pid to vector\n"; *1/ */ 2498 2499 /* /1* TLorentzVector truth_momentum; *1/ */ 2500 /* /1* truth_momentum.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e()); *1/ */ 2501 /* m_truthE.push_back(truth_momentum.E()); */ 2502 /* m_truthEta.push_back(truth_momentum.PseudoRapidity()); */ 2503 /* m_truthPhi.push_back(truth_momentum.Phi()); */ 2504 /* m_truthPt.push_back(truth_momentum.Pt()); */ 2505 /* m_truthMass.push_back(truth_momentum.M()); */ 2506 /* /1* std::cout << "Done adding kinematics to vector\n"; *1/ */ 2507 2508 /* /1* std::cout << "Asserting end_vtx= " << end_vtx << "\n"; *1/ */ 2509 /* assert(end_vtx); */ 2510 /* m_truthEndVtx_x.push_back(end_vtx->get_x()); */ 2511 /* m_truthEndVtx_y.push_back(end_vtx->get_y()); */ 2512 /* m_truthEndVtx_z.push_back(end_vtx->get_z()); */ 2513 /* m_truthEndVtx_t.push_back(end_vtx->get_t()); */ 2514 /* float x = end_vtx->get_x(); */ 2515 /* float y = end_vtx->get_y(); */ 2516 /* float r = sqrt(x*x + y*y); */ 2517 /* m_truthEndVtx_r.push_back(r); */ 2518 /* /1* std::cout << "Done adding vertex to vector\n"; *1/ */ 2519 2520 /* secondaryBarcodes.push_back(part->get_barcode()); */ 2521 /* } */ 2522 2523 2524 /* bool pythiaEMCalAna::vector_contains(std::pair<int,int> val, std::vector<std::pair<int,int>> vec) { */ 2525 /* unsigned int s = vec.size(); */ 2526 /* if ( s == 0 ) return false; */ 2527 /* unsigned int i; */ 2528 /* for (i=0; i<s; i++) { */ 2529 /* if (vec.at(i) == val) break; */ 2530 /* } */ 2531 /* if ( i == s ) return false; */ 2532 /* else return true; */ 2533 /* } */ 2534
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|