MC-Calibrations/src/JESMCTreeGen.cc

0001 #include "JESMCTreeGen.h"
0002
0003 // Vertex.
0004 #include <globalvertex/GlobalVertex.h>
0005 #include <globalvertex/GlobalVertexMap.h>
0006
0007 // Tower.
0008 #include <calobase/TowerInfo.h>
0009 #include <calobase/TowerInfov1.h>
0010 #include <calobase/TowerInfov2.h>
0011 #include <calobase/TowerInfov3.h>
0012 #include <calobase/TowerInfoContainer.h>
0013 #include <calobase/TowerInfoContainerv1.h>
0014 #include <calobase/TowerInfoContainerv2.h>
0015 #include <calobase/TowerInfoContainerv3.h>
0016 #include <calobase/TowerInfoContainerv4.h>
0017 #include <calobase/TowerInfoDefs.h>
0018
0019 // Cluster.
0020 #include <calobase/RawCluster.h>
0021 #include <calobase/RawClusterv1.h>
0022 #include <calobase/RawClusterContainer.h>
0023 #include <calobase/RawClusterUtility.h>
0024 #include <calobase/RawTowerDefs.h>
0025 #include <CLHEP/Vector/ThreeVector.h>
0026
0027 // Jet.
0028 #include <jetbase/Jetv1.h>
0029 #include <jetbase/Jetv2.h>
0030 #include <jetbase/JetContainer.h>
0031
0032 // Truth particle.
0033 #include <g4main/PHG4Particle.h>
0034 #include <g4main/PHG4TruthInfoContainer.h>
0035 #include <g4main/PHG4VtxPoint.h>
0036 #include <TDatabasePDG.h>
0037
0038 // Fun4All.
0039 #include <fun4all/Fun4AllReturnCodes.h>
0040 #include <fun4all/Fun4AllServer.h>
0041
0042 // Nodes.
0043 #include <phool/PHCompositeNode.h>
0044 #include <phool/getClass.h>
0045 #include <phool/phool.h>
0046 #include <phool/PHIODataNode.h>
0047 #include <phool/PHNode.h>
0048 #include <phool/PHNodeIterator.h>
0049 #include <phool/PHObject.h>
0050
0051 // General.
0052 #include <cstdint>
0053 #include <iostream>
0054 #include <utility>
0055
0056 class PHCompositeNode;
0057 class TowerInfoContainer;
0058 class JESMCTreeGen;
0059
0060 JESMCTreeGen::JESMCTreeGen(const std::string &name, const std::string &outfilename)
0061   :SubsysReco(name)
0062 {
0063   verbosity = 0;
0064   foutname = outfilename;
0065   std::cout << "JESMCTreeGen::JESMCTreeGen(const std::string &name) Calling ctor" << std::endl;
0066 }
0067
0068 int JESMCTreeGen::Init(PHCompositeNode * /*topNode*/) {
0069   if (verbosity > 0) std::cout << "Processing initialization: CaloEmulatorTreeMaker::Init(PHCompositeNode *topNode)" << std::endl;
0070   file = new TFile( foutname.c_str(), "RECREATE");
0071   tree = new TTree("ttree","TTree for JES calibration");
0072
0073   // Vertex information.
0074   if (doZVertexReadout) Initialize_z_vertex();
0075   // Tower information.
0076   if (doTowerInfoReadout) {
0077     Initialize_calo_tower();
0078     if (doTowerInfoReadout_RetowerEMCal) Initialize_calo_tower_RetowerEMCal();
0079   }
0080   // Cluster information.
0081   if (doClusterReadout) {
0082     Initialize_cluster();
0083     if (doClusterIsoReadout) Initialize_clusterIso();
0084     if (doClusterCNNReadout) Initialize_clusterCNN();
0085   }
0086   // Jet information.
0087   if (doRecoJetReadout) {
0088     for (int ir = 0; ir < static_cast<int>(doRecoJet_radius.size()); ir++) {
0089       Initialize_jet("recojet", doRecoJet_radius[ir]);
0090     }
0091     if (doSeedJetRawReadout) Initialize_jet("seedjetraw", 0.2);
0092     if (doSeedJetSubReadout) Initialize_jet("seedjetsub", 0.2);
0093   }
0094   if (doTruthJetReadout) {
0095     for (int ir = 0; ir < static_cast<int>(doTruthJet_radius.size()); ir++) {
0096       Initialize_jet("truthjet", doTruthJet_radius[ir]); // Only 0.2, 0.3, 0.4, 0.5 are produced.
0097     }
0098   }
0099   // Truth particle information.
0100   if (doTruthParticleReadout) Initialize_truthparticle();
0101
0102   ievent = 0;
0103   return Fun4AllReturnCodes::EVENT_OK;
0104 }
0105
0106 int JESMCTreeGen::process_event(PHCompositeNode *topNode) {
0107   if (verbosity >= 0) {
0108     if (ievent%100 == 0) std::cout << "Processing event " << ievent << std::endl;
0109   }
0110
0111   // Vertex information.
0112   if (doZVertexReadout) Fill_z_vertex(topNode);
0113   // Tower information.
0114   if (doTowerInfoReadout) {
0115     Fill_calo_tower(topNode, "CEMC");
0116     Fill_calo_tower(topNode, "HCALIN");
0117     Fill_calo_tower(topNode, "HCALOUT");
0118     if (doTowerInfoReadout_RetowerEMCal) Fill_emcal_retower(topNode);
0119   }
0120   // Cluster information.
0121   if (doClusterReadout) Fill_cluster(topNode);
0122   // Jet information.
0123   if (doRecoJetReadout) {
0124     for (int ir = 0; ir < static_cast<int>(doRecoJet_radius.size()); ir++) {
0125       Fill_recojet(topNode, "recojet", "subtracted", doRecoJet_radius[ir]);
0126     }
0127     if (doSeedJetRawReadout) Fill_recojet(topNode, "TowerInfo_HIRecoSeedsRaw", "seedjetraw", 0.2);
0128     if (doSeedJetSubReadout) Fill_recojet(topNode, "TowerInfo_HIRecoSeedsSub", "seedjetsub", 0.2);
0129   }
0130   if (doTruthJetReadout) {
0131     for (int ir = 0; ir < static_cast<int>(doTruthJet_radius.size()); ir++) {
0132       Fill_truthjet(topNode, doTruthJet_radius[ir]); // Only 0.2 - 0.8 are produced.
0133     }
0134   }
0135   // Truth particle information.
0136   if (doTruthParticleReadout) Fill_truthparticle(topNode);
0137
0138   tree->Fill();
0139   ievent++;
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143 int JESMCTreeGen::ResetEvent(PHCompositeNode * /*topNode*/) {
0144   if (Verbosity() > 1) std::cout << "Resetting the tree branches" << std::endl;
0145
0146   if (doClusterReadout) Reset_cluster();
0147   if (doRecoJetReadout || doTruthJetReadout) Reset_jet();
0148   if (doTruthParticleReadout) Reset_truthparticle();
0149
0150   return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 int JESMCTreeGen::End(PHCompositeNode * /*topNode*/) {
0154   if (verbosity > 0) std::cout << "Processing end: CaloEmulatorTreeMaker::End(PHCompositeNode *topNode)" << std::endl;
0155   file->cd();
0156   tree->Write();
0157   file->Close();
0158   delete file;
0159   if (verbosity > 0) std::cout << "CaloEmulatorTreeMaker complete." << std::endl;
0160   return Fun4AllReturnCodes::EVENT_OK;
0161 }
0162
0163 ////////// ********** Initialize functions ********** //////////
0164 void JESMCTreeGen::Initialize_z_vertex() {
0165   tree->Branch("z_vertex", &z_vertex, "z_vertex/F");
0166 }
0167
0168 void JESMCTreeGen::Initialize_calo_tower() {
0169   tree->Branch("emcal_tower_e", emcal_tower_e, "emcal_tower_e[96][256]/F"); // EMCal tower is original tower here.
0170   tree->Branch("ihcal_tower_e", ihcal_tower_e, "ihcal_tower_e[24][64]/F");
0171   tree->Branch("ohcal_tower_e", ohcal_tower_e, "ohcal_tower_e[24][64]/F");
0172   tree->Branch("emcal_tower_time", emcal_tower_time, "emcal_tower_time[96][256]/F"); // unit is sample.
0173   tree->Branch("ihcal_tower_time", ihcal_tower_time, "ihcal_tower_time[24][64]/F");
0174   tree->Branch("ohcal_tower_time", ohcal_tower_time, "ohcal_tower_time[24][64]/F");
0175   tree->Branch("emcal_tower_status", emcal_tower_status, "emcal_tower_status[96][256]/I"); // 0 is good, 1 is bad.
0176   tree->Branch("ihcal_tower_status", ihcal_tower_status, "ihcal_tower_status[24][64]/I");
0177   tree->Branch("ohcal_tower_status", ohcal_tower_status, "ohcal_tower_status[24][64]/I");
0178 }
0179
0180 void JESMCTreeGen::Initialize_calo_tower_RetowerEMCal() {
0181   tree->Branch("emcal_retower_e", emcal_retower_e, "emcal_retower_e[24][64]/F");
0182 }
0183
0184 void JESMCTreeGen::Initialize_cluster() {
0185   tree->Branch("cluster_e", &cluster_e);
0186   tree->Branch("cluster_pt", &cluster_pt);
0187   tree->Branch("cluster_eta", &cluster_eta);
0188   tree->Branch("cluster_phi", &cluster_phi);
0189   tree->Branch("cluster_e_core", &cluster_e_core);
0190   tree->Branch("cluster_pt_core", &cluster_pt_core);
0191   tree->Branch("cluster_eta_core", &cluster_eta_core);
0192   tree->Branch("cluster_phi_core", &cluster_phi_core);
0193   tree->Branch("cluster_probability", &cluster_probability);
0194 }
0195
0196 void JESMCTreeGen::Initialize_clusterIso() {
0197   tree->Branch("cluster_iso", &cluster_iso);
0198 }
0199
0200 void JESMCTreeGen::Initialize_clusterCNN() {
0201   tree->Branch("cluster_e_cnn", &cluster_e_cnn);
0202   tree->Branch("cluster_probability_cnn", &cluster_probability_cnn);
0203 }
0204
0205 void JESMCTreeGen::Initialize_jet(std::string jet_prefix , float radius) {
0206   // Prefix for the jet.
0207   std::string vector_prefix = jet_prefix + "0" + std::to_string((int)(radius * 10)) + "_";
0208   // Jet general information.
0209   std::string jet_e_vectorname = vector_prefix + "e"; jet_e_map[jet_e_vectorname] = std::vector<float>(); tree->Branch(jet_e_vectorname.c_str(), &jet_e_map[jet_e_vectorname]);
0210   std::string jet_et_vectorname = vector_prefix + "et"; jet_et_map[jet_et_vectorname] = std::vector<float>(); tree->Branch(jet_et_vectorname.c_str(), &jet_et_map[jet_et_vectorname]);
0211   std::string jet_pt_vectorname = vector_prefix + "pt"; jet_pt_map[jet_pt_vectorname] = std::vector<float>(); tree->Branch(jet_pt_vectorname.c_str(), &jet_pt_map[jet_pt_vectorname]);
0212   std::string jet_eta_vectorname = vector_prefix + "eta"; jet_eta_map[jet_eta_vectorname] = std::vector<float>(); tree->Branch(jet_eta_vectorname.c_str(), &jet_eta_map[jet_eta_vectorname]);
0213   std::string jet_phi_vectorname = vector_prefix + "phi"; jet_phi_map[jet_phi_vectorname] = std::vector<float>(); tree->Branch(jet_phi_vectorname.c_str(), &jet_phi_map[jet_phi_vectorname]);
0214   // Jet tower information or truth jet particle information.
0215   if (jet_prefix == "truthjet") {
0216     std::string jet_truthparticle_e_vectorname = vector_prefix + "truthparticle_e"; jet_truthparticle_e_map[jet_truthparticle_e_vectorname] = std::vector<float>(); tree->Branch(jet_truthparticle_e_vectorname.c_str(), &jet_truthparticle_e_map[jet_truthparticle_e_vectorname]);
0217     std::string jet_truthparticle_pt_vectorname = vector_prefix + "truthparticle_pt"; jet_truthparticle_pt_map[jet_truthparticle_pt_vectorname] = std::vector<float>(); tree->Branch(jet_truthparticle_pt_vectorname.c_str(), &jet_truthparticle_pt_map[jet_truthparticle_pt_vectorname]);
0218     std::string jet_truthparticle_eta_vectorname = vector_prefix + "truthparticle_eta"; jet_truthparticle_eta_map[jet_truthparticle_eta_vectorname] = std::vector<float>(); tree->Branch(jet_truthparticle_eta_vectorname.c_str(), &jet_truthparticle_eta_map[jet_truthparticle_eta_vectorname]);
0219     std::string jet_truthparticle_phi_vectorname = vector_prefix + "truthparticle_phi"; jet_truthparticle_phi_map[jet_truthparticle_phi_vectorname] = std::vector<float>(); tree->Branch(jet_truthparticle_phi_vectorname.c_str(), &jet_truthparticle_phi_map[jet_truthparticle_phi_vectorname]);
0220     std::string jet_truthparticle_pid_vectorname = vector_prefix + "truthparticle_pid"; jet_truthparticle_pid_map[jet_truthparticle_pid_vectorname] = std::vector<int>(); tree->Branch(jet_truthparticle_pid_vectorname.c_str(), &jet_truthparticle_pid_map[jet_truthparticle_pid_vectorname]);
0221   } else {
0222     // EMCal tower information. (Retower)
0223     std::string jet_emcal_ncomponent_vectorname = vector_prefix + "emcal_ncomponent"; jet_tower_ncomponent_map[jet_emcal_ncomponent_vectorname] = std::vector<int>(); tree->Branch(jet_emcal_ncomponent_vectorname.c_str(), &jet_tower_ncomponent_map[jet_emcal_ncomponent_vectorname]);
0224     std::string jet_emcal_tower_e_vectorname = vector_prefix + "emcal_tower_e"; jet_tower_e_map[jet_emcal_tower_e_vectorname] = std::vector<float>(); tree->Branch(jet_emcal_tower_e_vectorname.c_str(), &jet_tower_e_map[jet_emcal_tower_e_vectorname]);
0225     std::string jet_emcal_tower_ieta_vectorname = vector_prefix + "emcal_tower_ieta"; jet_tower_ieta_map[jet_emcal_tower_ieta_vectorname] = std::vector<int>(); tree->Branch(jet_emcal_tower_ieta_vectorname.c_str(), &jet_tower_ieta_map[jet_emcal_tower_ieta_vectorname]);
0226     std::string jet_emcal_tower_iphi_vectorname = vector_prefix + "emcal_tower_iphi"; jet_tower_iphi_map[jet_emcal_tower_iphi_vectorname] = std::vector<int>(); tree->Branch(jet_emcal_tower_iphi_vectorname.c_str(), &jet_tower_iphi_map[jet_emcal_tower_iphi_vectorname]);
0227     // iHCaL tower information.
0228     std::string jet_ihcal_ncomponent_vectorname = vector_prefix + "ihcal_ncomponent"; jet_tower_ncomponent_map[jet_ihcal_ncomponent_vectorname] = std::vector<int>(); tree->Branch(jet_ihcal_ncomponent_vectorname.c_str(), &jet_tower_ncomponent_map[jet_ihcal_ncomponent_vectorname]);
0229     std::string jet_ihcal_tower_e_vectorname = vector_prefix + "ihcal_tower_e"; jet_tower_e_map[jet_ihcal_tower_e_vectorname] = std::vector<float>(); tree->Branch(jet_ihcal_tower_e_vectorname.c_str(), &jet_tower_e_map[jet_ihcal_tower_e_vectorname]);
0230     std::string jet_ihcal_tower_ieta_vectorname = vector_prefix + "ihcal_tower_ieta"; jet_tower_ieta_map[jet_ihcal_tower_ieta_vectorname] = std::vector<int>(); tree->Branch(jet_ihcal_tower_ieta_vectorname.c_str(), &jet_tower_ieta_map[jet_ihcal_tower_ieta_vectorname]);
0231     std::string jet_ihcal_tower_iphi_vectorname = vector_prefix + "ihcal_tower_iphi"; jet_tower_iphi_map[jet_ihcal_tower_iphi_vectorname] = std::vector<int>(); tree->Branch(jet_ihcal_tower_iphi_vectorname.c_str(), &jet_tower_iphi_map[jet_ihcal_tower_iphi_vectorname]);
0232     // oHCaL tower information.
0233     std::string jet_ohcal_ncomponent_vectorname = vector_prefix + "ohcal_ncomponent"; jet_tower_ncomponent_map[jet_ohcal_ncomponent_vectorname] = std::vector<int>(); tree->Branch(jet_ohcal_ncomponent_vectorname.c_str(), &jet_tower_ncomponent_map[jet_ohcal_ncomponent_vectorname]);
0234     std::string jet_ohcal_tower_e_vectorname = vector_prefix + "ohcal_tower_e"; jet_tower_e_map[jet_ohcal_tower_e_vectorname] = std::vector<float>(); tree->Branch(jet_ohcal_tower_e_vectorname.c_str(), &jet_tower_e_map[jet_ohcal_tower_e_vectorname]);
0235     std::string jet_ohcal_tower_ieta_vectorname = vector_prefix + "ohcal_tower_ieta"; jet_tower_ieta_map[jet_ohcal_tower_ieta_vectorname] = std::vector<int>(); tree->Branch(jet_ohcal_tower_ieta_vectorname.c_str(), &jet_tower_ieta_map[jet_ohcal_tower_ieta_vectorname]);
0236     std::string jet_ohcal_tower_iphi_vectorname = vector_prefix + "ohcal_tower_iphi"; jet_tower_iphi_map[jet_ohcal_tower_iphi_vectorname] = std::vector<int>(); tree->Branch(jet_ohcal_tower_iphi_vectorname.c_str(), &jet_tower_iphi_map[jet_ohcal_tower_iphi_vectorname]);
0237   }
0238 }
0239
0240 void JESMCTreeGen::Initialize_truthparticle() {
0241   tree->Branch("truthparticle_e", &truthparticle_e);
0242   tree->Branch("truthparticle_pt", &truthparticle_pt);
0243   tree->Branch("truthparticle_eta", &truthparticle_eta);
0244   tree->Branch("truthparticle_phi", &truthparticle_phi);
0245   tree->Branch("truthparticle_pid", &truthparticle_pid);
0246 }
0247
0248 ////////// ********** Fill functions ********** //////////
0249 void JESMCTreeGen::Fill_z_vertex(PHCompositeNode *topNode) {
0250   z_vertex = -999;
0251   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0252   if (vertexmap) {
0253     if (!vertexmap->empty()) {
0254       GlobalVertex *vtx = vertexmap->begin()->second;
0255       if (vtx) {
0256         z_vertex = vtx->get_z();
0257       }
0258     } else {
0259       std::cout << "GlobalVertexMap is empty" << std::endl;
0260     }
0261   } else {
0262     std::cout << "GlobalVertexMap is missing" << std::endl;
0263   }
0264 }
0265
0266 void JESMCTreeGen::Fill_calo_tower(PHCompositeNode *topNode, std::string calorimeter) {
0267   std::string tower_info_container_name = "TOWERINFO_CALIB_" + calorimeter;
0268   TowerInfoContainer *_towers_calo = findNode::getClass<TowerInfoContainer>(topNode, tower_info_container_name);
0269   if (_towers_calo) {
0270     for (int channel = 0; channel < (int)_towers_calo->size(); ++channel) {
0271         TowerInfo *_tower = _towers_calo->get_tower_at_channel(channel);
0272         unsigned int towerkey = _towers_calo->encode_key(channel);
0273         int etabin = _towers_calo->getTowerEtaBin(towerkey);
0274         int phibin = _towers_calo->getTowerPhiBin(towerkey);
0275       float towE = _tower->get_energy();
0276       float towT = _tower->get_time_float();
0277       int towS = 1;
0278       if (!_tower->get_isHot() && !_tower->get_isNoCalib() && !_tower->get_isNotInstr() && !_tower->get_isBadChi2()) {
0279         towS = 0;
0280       }
0281         if (calorimeter == "CEMC") {
0282         emcal_tower_e[etabin][phibin] = towE;
0283         emcal_tower_time[etabin][phibin] = towT;
0284         emcal_tower_status[etabin][phibin] = towS;
0285       } else if (calorimeter == "HCALIN") {
0286         ihcal_tower_e[etabin][phibin] = towE;
0287         ihcal_tower_time[etabin][phibin] = towT;
0288         ihcal_tower_status[etabin][phibin] = towS;
0289       } else if (calorimeter == "HCALOUT") {
0290         ohcal_tower_e[etabin][phibin] = towE;
0291         ohcal_tower_time[etabin][phibin] = towT;
0292         ohcal_tower_status[etabin][phibin] = towS;
0293       }
0294     }
0295   } else {
0296     std::cout << "TowerInfoContainer for " << calorimeter << " is missing" << std::endl;
0297   }
0298 }
0299
0300 void JESMCTreeGen::Fill_emcal_retower(PHCompositeNode *topNode) {
0301   TowerInfoContainer *_towers_calo = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0302   if (_towers_calo) {
0303     for (int channel = 0; channel < (int)_towers_calo->size(); ++channel) {
0304         TowerInfo *_tower = _towers_calo->get_tower_at_channel(channel);
0305         unsigned int towerkey = _towers_calo->encode_key(channel);
0306         int etabin = _towers_calo->getTowerEtaBin(towerkey);
0307         int phibin = _towers_calo->getTowerPhiBin(towerkey);
0308       float towE = _tower->get_energy();
0309       emcal_retower_e[etabin][phibin] = towE;
0310     }
0311   } else {
0312     std::cout << "TowerInfoContainer EMCal retower node is missing" << std::endl;
0313   }
0314 }
0315
0316 void JESMCTreeGen::Fill_cluster(PHCompositeNode *topNode) {
0317   std::string _clustername;
0318   _clustername = "CLUSTERINFO_CEMC";
0319   RawClusterContainer *_clusters = findNode::getClass<RawClusterContainer>(topNode, _clustername.c_str());
0320   if (_clusters) {
0321     RawClusterContainer::ConstRange clusterEnd = _clusters->getClusters();
0322       for (RawClusterContainer::ConstIterator clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++) {
0323         RawCluster *_cluster = clusterIter->second;
0324       CLHEP::Hep3Vector vertex(0, 0, 0); // Vertex is assumed to be at (0, 0, 0). This is a placeholder for vertex correction.
0325       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*_cluster, vertex);
0326       CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*_cluster, vertex);
0327       float clusE = E_vec_cluster_Full.mag();
0328       if (clusE < 1) continue; // Skip clusters with energy less than 1 GeV.
0329       float clusEcore = E_vec_cluster.mag();
0330       float clus_eta = E_vec_cluster_Full.pseudoRapidity();
0331       float clus_phi = E_vec_cluster_Full.phi();
0332       float clus_pt = E_vec_cluster_Full.perp();
0333       float clus_coreeta = E_vec_cluster.pseudoRapidity();
0334       float clus_corephi = E_vec_cluster.phi();
0335       float clus_corept = E_vec_cluster.perp();
0336       float clus_prob = _cluster->get_prob();
0337       cluster_e.push_back(clusE);
0338         cluster_pt.push_back(clus_pt);
0339         cluster_eta.push_back(clus_eta);
0340       cluster_phi.push_back(clus_phi);
0341       cluster_e_core.push_back(clusEcore);
0342       cluster_pt_core.push_back(clus_corept);
0343       cluster_eta_core.push_back(clus_coreeta);
0344       cluster_phi_core.push_back(clus_corephi);
0345       cluster_probability.push_back(clus_prob);
0346       if (doClusterIsoReadout) {
0347         float clus_iso = _cluster->get_et_iso(4, false, true);
0348         cluster_iso.push_back(clus_iso);
0349       }
0350       }
0351   } else {
0352     std::cout << "RawClusterContainer is missing" << std::endl;
0353   }
0354   if (doClusterCNNReadout) {
0355     std::string _clustername_cnn;
0356     _clustername_cnn = "CLUSTERINFO_CEMC_CNN";
0357     RawClusterContainer *_clusters_cnn = findNode::getClass<RawClusterContainer>(topNode, _clustername_cnn.c_str());
0358     if (_clusters_cnn) {
0359       RawClusterContainer::ConstRange clusterEnd_cnn = _clusters_cnn->getClusters();
0360       for (RawClusterContainer::ConstIterator clusterIter = clusterEnd_cnn.first; clusterIter != clusterEnd_cnn.second; clusterIter++) {
0361         RawCluster *_cluster = clusterIter->second;
0362         float clus_prob_cnn = _cluster->get_prob();
0363         float clusE = _cluster->get_energy();
0364         cluster_probability_cnn.push_back(clus_prob_cnn);
0365         cluster_e_cnn.push_back(clusE);
0366       }
0367     } else {
0368       std::cout << "RawClusterContainer for CNN cluster is missing" << std::endl;
0369     }
0370   }
0371 }
0372
0373 void JESMCTreeGen::Fill_recojet(PHCompositeNode *topNode, std::string jet_prefix, std::string node_prefix, float radius) {
0374   // Set vector name prefix for the jet map.
0375   std::string vector_prefix = jet_prefix + "0" + std::to_string((int)(radius * 10)) + "_";
0376   // Get the jet container.
0377   std::string nodename = "AntiKt_" + node_prefix + "_r0" + std::to_string((int)(radius * 10));
0378   JetContainer* _jets = findNode::getClass<JetContainer>(topNode, nodename);
0379   if (!_jets) std::cout << "JetContainer for " << nodename << " is missing" << std::endl;
0380   // Get the towerinfo container.
0381   TowerInfoContainer *_emcal_towers = nullptr;
0382   TowerInfoContainer *_ihcal_towers = nullptr;
0383   TowerInfoContainer *_ohcal_towers = nullptr;
0384   if (node_prefix == "TowerInfo_HIRecoSeedsRaw") {
0385     _emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0386     _ihcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0387     _ohcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0388   } else {
0389     _emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0390     _ihcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0391     _ohcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0392   }
0393   if (!_emcal_towers || !_ihcal_towers || !_ohcal_towers) std::cout << "TowerInfoContainer is missing" << std::endl;
0394   // Fill the jet information.
0395   for (auto _jet : *_jets) {
0396     if (_jet->get_pt() < 1) continue; // Skip jets with pt less than 1 GeV. (same cut in jet reconstruction)
0397     jet_e_map[vector_prefix+"e"].push_back(_jet->get_e());
0398     jet_et_map[vector_prefix+"et"].push_back(_jet->get_et());
0399     jet_pt_map[vector_prefix+"pt"].push_back(_jet->get_pt());
0400     jet_eta_map[vector_prefix+"eta"].push_back(_jet->get_eta());
0401     jet_phi_map[vector_prefix+"phi"].push_back(_jet->get_phi());
0402     int ncomponent_emcal = 0;
0403     int ncomponent_ihcal = 0;
0404     int ncomponent_ohcal = 0;
0405     for (auto comp: _jet->get_comp_vec()) {
0406         unsigned int channel = comp.second;
0407       if (comp.first == 14 || comp.first == 29 || comp.first == 25 || comp.first == 28) {
0408         TowerInfo *_tower = _emcal_towers->get_tower_at_channel(channel);
0409           unsigned int calokey = _emcal_towers->encode_key(channel);
0410           int jet_tower_ieta = _emcal_towers->getTowerEtaBin(calokey);
0411           int jet_tower_iphi = _emcal_towers->getTowerPhiBin(calokey);
0412           float jet_tower_e = _tower->get_energy();
0413         jet_tower_e_map[vector_prefix+"emcal_tower_e"].push_back(jet_tower_e);
0414         jet_tower_ieta_map[vector_prefix+"emcal_tower_ieta"].push_back(jet_tower_ieta);
0415         jet_tower_iphi_map[vector_prefix+"emcal_tower_iphi"].push_back(jet_tower_iphi);
0416         ncomponent_emcal++;
0417         } else if (comp.first == 15 ||  comp.first == 30 || comp.first == 26) {
0418         TowerInfo *_tower = _ihcal_towers->get_tower_at_channel(channel);
0419           unsigned int calokey = _ihcal_towers->encode_key(channel);
0420           int jet_tower_ieta = _ihcal_towers->getTowerEtaBin(calokey);
0421           int jet_tower_iphi = _ihcal_towers->getTowerPhiBin(calokey);
0422           float jet_tower_e = _tower->get_energy();
0423         jet_tower_e_map[vector_prefix+"ihcal_tower_e"].push_back(jet_tower_e);
0424         jet_tower_ieta_map[vector_prefix+"ihcal_tower_ieta"].push_back(jet_tower_ieta);
0425         jet_tower_iphi_map[vector_prefix+"ihcal_tower_iphi"].push_back(jet_tower_iphi);
0426         ncomponent_ihcal++;
0427         } else if (comp.first == 16 || comp.first == 31 || comp.first == 27) {
0428         TowerInfo *_tower = _ohcal_towers->get_tower_at_channel(channel);
0429           unsigned int calokey = _ohcal_towers->encode_key(channel);
0430           int jet_tower_ieta = _ohcal_towers->getTowerEtaBin(calokey);
0431           int jet_tower_iphi = _ohcal_towers->getTowerPhiBin(calokey);
0432           float jet_tower_e = _tower->get_energy();
0433         jet_tower_e_map[vector_prefix+"ohcal_tower_e"].push_back(jet_tower_e);
0434         jet_tower_ieta_map[vector_prefix+"ohcal_tower_ieta"].push_back(jet_tower_ieta);
0435         jet_tower_iphi_map[vector_prefix+"ohcal_tower_iphi"].push_back(jet_tower_iphi);
0436         ncomponent_ohcal++;
0437       }
0438       }
0439     jet_tower_ncomponent_map[vector_prefix+"emcal_ncomponent"].push_back(ncomponent_emcal);
0440     jet_tower_ncomponent_map[vector_prefix+"ihcal_ncomponent"].push_back(ncomponent_ihcal);
0441     jet_tower_ncomponent_map[vector_prefix+"ohcal_ncomponent"].push_back(ncomponent_ohcal);
0442   }
0443 }
0444
0445 void JESMCTreeGen::Fill_truthjet(PHCompositeNode *topNode, float radius) {
0446   // Set vector name prefix for the jet map.
0447   std::string vector_prefix = "truthjet0" + std::to_string((int)(radius * 10)) + "_";
0448   // Get the jet container.
0449   std::string nodename = "AntiKt_Truth_r0" + std::to_string((int)(radius * 10));
0450   JetContainer* _truth_jets = findNode::getClass<JetContainer>(topNode, nodename);
0451   if (!_truth_jets) std::cout << "JetContainer for " << nodename << " is missing" << std::endl;
0452   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0453   if (_truth_jets && truthinfo) {
0454     for (auto _truthjet : *_truth_jets) {
0455       if (_truthjet->get_eta() < -2 || _truthjet->get_eta() > 2 || _truthjet->get_pt() < 5 || _truthjet->get_pt() > 100) continue; // Cut for truth jets.
0456       jet_e_map[vector_prefix+"e"].push_back(_truthjet->get_e());
0457       jet_et_map[vector_prefix+"et"].push_back(_truthjet->get_et());
0458       jet_pt_map[vector_prefix+"pt"].push_back(_truthjet->get_pt());
0459       jet_eta_map[vector_prefix+"eta"].push_back(_truthjet->get_eta());
0460       jet_phi_map[vector_prefix+"phi"].push_back(_truthjet->get_phi());
0461       for (auto comp: _truthjet->get_comp_vec()) {
0462           PHG4Particle* _truthparticle = truthinfo->GetParticle(comp.second);
0463         float particle_e = _truthparticle->get_e();
0464         float particle_pt = sqrt(_truthparticle->get_px() * _truthparticle->get_px() + _truthparticle->get_py() * _truthparticle->get_py());
0465         float particle_eta = atanh(_truthparticle->get_pz() / sqrt(_truthparticle->get_px()*_truthparticle->get_px()+_truthparticle->get_py()*_truthparticle->get_py()+_truthparticle->get_pz()*_truthparticle->get_pz()));
0466         float particle_phi = atan2(_truthparticle->get_py(), _truthparticle->get_px());
0467         int particle_pid = _truthparticle->get_pid();
0468         jet_truthparticle_e_map[vector_prefix+"truthparticle_e"].push_back(particle_e);
0469         jet_truthparticle_pt_map[vector_prefix+"truthparticle_pt"].push_back(particle_pt);
0470         jet_truthparticle_eta_map[vector_prefix+"truthparticle_eta"].push_back(particle_eta);
0471         jet_truthparticle_phi_map[vector_prefix+"truthparticle_phi"].push_back(particle_phi);
0472         jet_truthparticle_pid_map[vector_prefix+"truthparticle_pid"].push_back(particle_pid);
0473       }
0474       }
0475   } else {
0476     std::cout << "JetContainer for " << nodename << " is missing, or PHG4TruthInfoContainer is missing" << std::endl;
0477   }
0478 }
0479
0480 void JESMCTreeGen::Fill_truthparticle(PHCompositeNode *topNode) {
0481   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0482   if (truthinfo) {
0483     PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0484     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) {
0485       PHG4Particle *truth = iter->second;
0486       if (!truthinfo->is_primary(truth)) continue;
0487       float particle_eta = atanh(truth->get_pz() / sqrt(truth->get_px()*truth->get_px()+truth->get_py()*truth->get_py()+truth->get_pz()*truth->get_pz()));
0488       if (fabs(particle_eta) > 1.6) continue;
0489       float particle_e = truth->get_e();
0490       float particle_pt = sqrt(truth->get_px() * truth->get_px() + truth->get_py() * truth->get_py());
0491       float particle_phi = atan2(truth->get_py(), truth->get_px());
0492       int particle_pid = truth->get_pid();
0493       truthparticle_e.push_back(particle_e);
0494       truthparticle_pt.push_back(particle_pt);
0495       truthparticle_eta.push_back(particle_eta);
0496       truthparticle_phi.push_back(particle_phi);
0497       truthparticle_pid.push_back(particle_pid);
0498     }
0499   } else {
0500     std::cout << "PHG4TruthInfoContainer is missing" << std::endl;
0501   }
0502 }
0503
0504 ////////// ********** Reset functions ********** //////////
0505 void JESMCTreeGen::Reset_cluster() {
0506   cluster_e.clear();
0507   cluster_pt.clear();
0508   cluster_eta.clear();
0509   cluster_phi.clear();
0510   cluster_e_core.clear();
0511   cluster_pt_core.clear();
0512   cluster_eta_core.clear();
0513   cluster_phi_core.clear();
0514   cluster_probability.clear();
0515   if (doClusterIsoReadout) {
0516     cluster_iso.clear();
0517   }
0518   if (doClusterCNNReadout) {
0519     cluster_e_cnn.clear();
0520     cluster_probability_cnn.clear();
0521   }
0522 }
0523
0524 void JESMCTreeGen::Reset_jet() {
0525   for (auto& entry : jet_e_map) entry.second.clear();
0526   for (auto& entry : jet_et_map) entry.second.clear();
0527   for (auto& entry : jet_pt_map) entry.second.clear();
0528   for (auto& entry : jet_eta_map) entry.second.clear();
0529   for (auto& entry : jet_phi_map) entry.second.clear();
0530   if (doRecoJetReadout) {
0531     for (auto& entry : jet_tower_ncomponent_map) entry.second.clear();
0532     for (auto& entry : jet_tower_e_map) entry.second.clear();
0533     for (auto& entry : jet_tower_ieta_map) entry.second.clear();
0534     for (auto& entry : jet_tower_iphi_map) entry.second.clear();
0535   }
0536   if (doTruthJetReadout) {
0537     for (auto& entry : jet_truthparticle_e_map) entry.second.clear();
0538     for (auto& entry : jet_truthparticle_pt_map) entry.second.clear();
0539     for (auto& entry : jet_truthparticle_eta_map) entry.second.clear();
0540     for (auto& entry : jet_truthparticle_phi_map) entry.second.clear();
0541     for (auto& entry : jet_truthparticle_pid_map) entry.second.clear();
0542   }
0543 }
0544
0545 void JESMCTreeGen::Reset_truthparticle() {
0546   truthparticle_e.clear();
0547   truthparticle_pt.clear();
0548   truthparticle_eta.clear();
0549   truthparticle_phi.clear();
0550   truthparticle_pid.clear();
0551 }