File indexing completed on 2025-08-05 08:12:42
0001
0002 #include <fun4all/Fun4AllBase.h>
0003 #include <iostream>
0004 #include "TruthCaloTree.h"
0005
0006
0007 #include <fun4all/Fun4AllServer.h>
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <phool/getClass.h>
0010 #include <phool/PHCompositeNode.h>
0011 #include <g4main/PHG4TruthInfoContainer.h>
0012 #include <g4main/PHG4Particle.h>
0013 #include <g4main/PHG4Hit.h>
0014 #include <g4main/PHG4VtxPoint.h>
0015
0016
0017 #include <calobase/RawTowerv2.h>
0018 #include <calobase/RawTowerGeom.h>
0019 #include <calobase/RawTowerContainer.h>
0020 #include <calobase/RawCluster.h>
0021 #include <calobase/RawClusterContainer.h>
0022
0023
0024 #include "TLorentzVector.h"
0025 #include "TMath.h"
0026 #include "TROOT.h"
0027 #include "TH2.h"
0028
0029 TruthCaloTree::TruthCaloTree(const std::string &name) : SubsysReco("TRUTH_CALO_TREE")
0030 {
0031 _foutname = name;
0032 _hcal_sim_name = "TOWER_SIM_HCALOUT";
0033 _hcalIN_sim_name = "TOWER_SIM_HCALIN";
0034 _EMcal_sim_name = "TOWER_SIM_CEMC";
0035
0036 }
0037
0038 int TruthCaloTree::reset_tree_vars() {
0039
0040 _b_event = -99;
0041 _b_tower_sim_n = -99;
0042 _b_EMcal_sim_n = -99;
0043 _b_hcalIN_sim_n = -99;
0044 _nBH = -99;
0045
0046 for (int i = 0; i <nTowers; i++){
0047
0048 _b_tower_sim_E[i] = -99;
0049 _b_tower_sim_eta[i] = -99;
0050 _b_tower_sim_phi[i] = -99;
0051 _b_tower_sim_ieta[i] = -99;
0052 _b_tower_sim_iphi[i] = -99;
0053
0054 _b_EMcal_sim_E[i] = -99;
0055 _b_EMcal_sim_eta[i] = -99;
0056 _b_EMcal_sim_phi[i] = -99;
0057 _b_EMcal_sim_iphi[i] = -99;
0058 _b_EMcal_sim_ieta[i] = -99;
0059
0060 _b_hcalIN_sim_E[i] = -99;
0061 _b_hcalIN_sim_eta[i] = -99;
0062 _b_hcalIN_sim_phi[i] = -99;
0063 _b_hcalIN_sim_iphi[i] = -99;
0064 _b_hcalIN_sim_ieta[i] = -99;
0065 }
0066
0067 _b_n_truth = -99;
0068 n_child = -99;
0069 n_vertex = -99;
0070
0071 for (int i = 0; i < nTruth; i++) {
0072 _b_truthenergy[i] = -99;
0073 _b_trutheta[i] = -99;
0074 _b_truthphi[i] = -99;
0075 _b_truthpt[i] = -99;
0076 _b_truthp[i] = -99;
0077 _b_truthpid[i] = -99;
0078 _b_truth_trackid[i] = -99;
0079 vertex_id[i] = -99;
0080 vertex_x[i] = -99;
0081 vertex_y[i] = -99;
0082 vertex_z[i] = -99;
0083 _BH_e[i] = -99;
0084 _BH_px[i] = -99;
0085 _BH_py[i] = -99;
0086 _BH_pz[i] = -99;
0087 _BH_track_id[i] = -99;
0088 }
0089
0090 for (int i = 0; i < nTruth; i++) {
0091 child_vertex_id[i] = -99;
0092 child_parent_id[i] = -99;
0093 child_pid[i] = -99;
0094 child_px[i] = -99;
0095 child_py[i] = -99;
0096 child_pz[i] = -99;
0097 child_energy[i] = -99;
0098 }
0099
0100 return 1;
0101 }
0102
0103 int TruthCaloTree::GetNodes(PHCompositeNode *topNode) {
0104
0105 blackhole = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BH_1");
0106 if (!blackhole) std::cout << "No blackhole" << std::endl;
0107
0108 if (_debug) std::cout<<"GettingNodes..."<<std::endl;
0109 _geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0110 if(!_geomOH) std::cout<<"No TOWERGeOM_HCALOUT"<<std::endl;
0111
0112 _geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0113 if(!_geomIH) std::cout<<"No TOWERGeOM_HCALIN"<<std::endl;
0114
0115 _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0116 if(!_geomEM) std::cout<<"No TOWERGeOM_CEMC"<<std::endl;
0117
0118 _towersSimOH = findNode::getClass<RawTowerContainer>(topNode, _hcal_sim_name);
0119 if (!_towersSimOH) std::cout<<"No TOWER_SIM_HCALOUT Node"<<std::endl;
0120
0121 _towersSimIH = findNode::getClass<RawTowerContainer>(topNode, _hcalIN_sim_name);
0122 if (!_towersSimIH) std::cout<<"No TOWER_SIM_HCALIN Node"<<std::endl;
0123
0124 _towersSimEM = findNode::getClass<RawTowerContainer>(topNode, _EMcal_sim_name);
0125 if (!_towersSimEM) std::cout<<"No TOWER_SIM_CEMC Node"<<std::endl;
0126
0127 truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0128 if (!truthinfo) std::cout << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"<< std::endl;
0129
0130 return 1;
0131 }
0132
0133 int TruthCaloTree::Init(PHCompositeNode *topNode) {
0134 _ievent = 0;
0135 _b_event = -1;
0136
0137 if (_debug) std::cout<<"Initiating..."<<std::endl;
0138
0139 reset_tree_vars();
0140
0141 _file = new TFile(_foutname.c_str(), "RECREATE");
0142
0143 _tree = new TTree("T", "keep on giving tree");
0144
0145 _tree->Branch("n_truth",&_b_n_truth,"n_truth/I");
0146 _tree->Branch("truthenergy",_b_truthenergy,"truthenergy[n_truth]/F");
0147 _tree->Branch("trutheta",_b_trutheta,"trutheta[n_truth]/F");
0148 _tree->Branch("truthphi",_b_truthphi,"truthphi[n_truth]/F");
0149 _tree->Branch("truthpt",_b_truthpt,"truthpt[n_truth]/F");
0150 _tree->Branch("truthp",_b_truthp,"truthp[n_truth]/F");
0151 _tree->Branch("truthpid",_b_truthpid,"truthpid[n_truth]/I");
0152 _tree->Branch("truth_track_id",_b_truth_trackid,"truth_track_id[n_truth]/I");
0153
0154 _tree->Branch("EMcal_sim_n", &_b_EMcal_sim_n, "EMcal_sim_n/I");
0155 _tree->Branch("EMcal_sim_E", _b_EMcal_sim_E, "EMcal_sim_E[EMcal_sim_n]/F");
0156 _tree->Branch("EMcal_sim_eta", _b_EMcal_sim_eta, "EMcal_sim_eta[EMcal_sim_n]/F");
0157 _tree->Branch("EMcal_sim_phi", _b_EMcal_sim_phi, "EMcal_sim_phi[EMcal_sim_n]/F");
0158 _tree->Branch("EMcal_sim_iphi", _b_EMcal_sim_iphi, "EMcal_sim_iphi[EMcal_sim_n]/I");
0159 _tree->Branch("EMcal_sim_ieta", _b_EMcal_sim_ieta, "EMcal_sim_ieta[EMcal_sim_n]/I");
0160
0161 _tree->Branch("hcalIN_sim_n", &_b_hcalIN_sim_n, "hcalIN_sim_n/I");
0162 _tree->Branch("hcalIN_sim_E", _b_hcalIN_sim_E, "hcalIN_sim_E[hcalIN_sim_n]/F");
0163 _tree->Branch("hcalIN_sim_eta", _b_hcalIN_sim_eta, "hcalIN_sim_eta[hcalIN_sim_n]/F");
0164 _tree->Branch("hcalIN_sim_phi", _b_hcalIN_sim_phi, "hcalIN_sim_phi[hcalIN_sim_n]/F");
0165 _tree->Branch("hcalIN_sim_iphi", _b_hcalIN_sim_iphi, "hcalIN_sim_iphi[hcalIN_sim_n]/I");
0166 _tree->Branch("hcalIN_sim_ieta", _b_hcalIN_sim_ieta, "hcalIN_sim_ieta[hcalIN_sim_n]/I");
0167
0168 _tree->Branch("tower_sim_n",&_b_tower_sim_n, "tower_sim_n/I");
0169 _tree->Branch("tower_sim_E",_b_tower_sim_E, "tower_sim_E[tower_sim_n]/F");
0170 _tree->Branch("tower_sim_eta",_b_tower_sim_eta, "tower_sim_eta[tower_sim_n]/F");
0171 _tree->Branch("tower_sim_phi",_b_tower_sim_phi, "tower_sim_phi[tower_sim_n]/F");
0172 _tree->Branch("tower_sim_ieta",_b_tower_sim_ieta, "tower_sim_ieta[tower_sim_n]/I");
0173 _tree->Branch("tower_sim_iphi",_b_tower_sim_iphi, "tower_sim_iphi[tower_sim_n]/I");
0174
0175 _tree->Branch("n_vertex",&n_vertex,"n_vertex/I");
0176 _tree->Branch("vertex_id",vertex_id,"vertex_id[n_vertex]/I");
0177 _tree->Branch("vertex_x",vertex_x,"vertex_x[n_vertex]/F");
0178 _tree->Branch("vertex_y",vertex_y,"vertex_y[n_vertex]/F");
0179 _tree->Branch("vertex_z",vertex_z,"vertex_z[n_vertex]/F");
0180
0181 _tree->Branch("n_child",&n_child,"n_child/I");
0182 _tree->Branch("child_vertex_id",child_vertex_id,"child_vertex_id[n_child]/I");
0183 _tree->Branch("child_parent_id",child_parent_id,"child_parent_id[n_child]/I");
0184 _tree->Branch("child_pid",child_pid,"child_pid[n_child]/I");
0185 _tree->Branch("child_px",child_px,"child_px[n_child]/F");
0186 _tree->Branch("child_py",child_py,"child_py[n_child]/F");
0187 _tree->Branch("child_pz",child_pz,"child_pz[n_child]/F");
0188 _tree->Branch("child_energy",child_energy,"child_energy[n_child]/F");
0189
0190 _tree->Branch("BH_n",&_nBH,"BH_n/I");
0191 _tree->Branch("BH_px",_BH_px,"BH_px[BH_n]/F");
0192 _tree->Branch("BH_py",_BH_py,"BH_py[BH_n]/F");
0193 _tree->Branch("BH_pz",_BH_pz,"BH_pz[BH_n]/F");
0194 _tree->Branch("BH_track_id",_BH_track_id,"BH_track_id[BH_n]/I");
0195 _tree->Branch("BH_e",_BH_e,"BH_e[BH_n]/F");
0196
0197
0198 return 0;
0199 }
0200
0201 int TruthCaloTree::process_event(PHCompositeNode *topNode) {
0202
0203 GetNodes(topNode);
0204
0205 _b_event = _ievent;
0206 if (_ievent %5000==0) std::cout<<"Event: "<<_ievent<<std::endl;
0207
0208
0209 if (_debug) std::cout<<"Processing Event: "<< _b_event<<std::endl;
0210 if (_debug) std::cout << "hello";
0211
0212
0213
0214
0215
0216 _nBH = 0;
0217 PHG4HitContainer::ConstRange bh_hit_range = blackhole->getHits();
0218 for (PHG4HitContainer::ConstIterator hit_iter = bh_hit_range.first;
0219 hit_iter != bh_hit_range.second; hit_iter++)
0220 {
0221 PHG4Hit *this_hit = hit_iter->second;
0222 assert(this_hit);
0223 _BH_e[_nBH] = this_hit->get_edep();
0224 _BH_px[_nBH] = this_hit->get_px(0);
0225 _BH_py[_nBH] = this_hit->get_py(0);
0226 _BH_pz[_nBH] = this_hit->get_pz(0);
0227 PHG4Particle *p = truthinfo->GetParticle(this_hit->get_trkid());
0228 _BH_track_id[_nBH] = p->get_track_id();
0229 _nBH++;
0230 }
0231
0232
0233
0234
0235
0236 _b_tower_sim_n=0;
0237 RawTowerContainer::ConstRange begin_end_sim = _towersSimOH->getTowers();
0238 if (_debug) std::cout<<"Got the iterator"<<std::endl;
0239
0240 for (RawTowerContainer::ConstIterator rtiter = begin_end_sim.first; rtiter != begin_end_sim.second; ++rtiter) {
0241 if ((_b_tower_sim_n%10==0)&&(_debug)) std::cout<<"At sim tower: "<< _b_tower_sim_n<<std::endl;
0242
0243 RawTower *tower = rtiter->second;
0244 RawTowerGeom *tower_geom = _geomOH->get_tower_geometry(tower->get_key());
0245 _b_tower_sim_E [ _b_tower_sim_n ] = tower->get_energy();
0246 _b_tower_sim_eta [ _b_tower_sim_n ] = tower_geom->get_eta();
0247 _b_tower_sim_phi [ _b_tower_sim_n ] = tower_geom->get_phi();
0248 _b_tower_sim_ieta[ _b_tower_sim_n ] = tower_geom->get_bineta();
0249 _b_tower_sim_iphi[ _b_tower_sim_n ] = tower_geom->get_binphi();
0250 _b_tower_sim_n++;
0251 if(_b_tower_sim_n >= nTowers){
0252 std::cout << __FILE__ << " ERROR: _b_tower_sim_n has hit cap of " << nTowers << "!!!" << std::endl;
0253 }
0254 }
0255
0256
0257
0258
0259
0260 _b_hcalIN_sim_n = 0;
0261 RawTowerContainer::ConstRange begin_end_simIN = _towersSimIH->getTowers();
0262 for (RawTowerContainer::ConstIterator rtiter = begin_end_simIN.first; rtiter != begin_end_simIN.second; ++rtiter) {
0263 RawTower *tower = rtiter->second;
0264 RawTowerGeom *tower_geom = _geomIH->get_tower_geometry(tower->get_key());
0265 _b_hcalIN_sim_E [ _b_hcalIN_sim_n ] = tower->get_energy();
0266 _b_hcalIN_sim_eta [ _b_hcalIN_sim_n ] = tower_geom->get_eta();
0267 _b_hcalIN_sim_phi [ _b_hcalIN_sim_n ] = tower_geom->get_phi();
0268 _b_hcalIN_sim_ieta [ _b_hcalIN_sim_n ] = tower_geom->get_bineta();
0269 _b_hcalIN_sim_iphi [ _b_hcalIN_sim_n ] = tower_geom->get_binphi();
0270 _b_hcalIN_sim_n++;
0271
0272 }
0273
0274
0275
0276
0277
0278 _b_EMcal_sim_n = 0;
0279 RawTowerContainer::ConstRange begin_end_simEM = _towersSimEM->getTowers();
0280 for (RawTowerContainer::ConstIterator rtiter = begin_end_simEM.first; rtiter != begin_end_simEM.second; ++rtiter) {
0281 RawTower *tower = rtiter->second;
0282 RawTowerGeom *tower_geom = _geomEM->get_tower_geometry(tower->get_key());
0283 _b_EMcal_sim_E [ _b_EMcal_sim_n ] = tower->get_energy();
0284 _b_EMcal_sim_eta [ _b_EMcal_sim_n ] = tower_geom->get_eta();
0285 _b_EMcal_sim_phi [ _b_EMcal_sim_n ] = tower_geom->get_phi();
0286 _b_EMcal_sim_ieta [ _b_EMcal_sim_n ] = tower_geom->get_bineta();
0287 _b_EMcal_sim_iphi [ _b_EMcal_sim_n ] = tower_geom->get_binphi();
0288 _b_EMcal_sim_n++;
0289
0290 }
0291
0292
0293
0294
0295
0296
0297 n_child = 0;
0298
0299 std::set<int> vertex;
0300
0301 if (!vertex.empty()) vertex.clear();
0302 PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0303 _b_n_truth = 0;
0304 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0305 iter != range.second; ++iter) {
0306 const PHG4Particle *truth = iter->second;
0307 if (!truthinfo->is_primary(truth)) continue;
0308 _b_truthp[_b_n_truth] = sqrt(truth->get_px() * truth->get_px()
0309 + truth->get_py() * truth->get_py()
0310 + truth->get_pz() * truth->get_pz());
0311 _b_truthpt[_b_n_truth] = sqrt(truth->get_px() * truth->get_px()
0312 + truth->get_py() * truth->get_py());
0313 _b_truthenergy[_b_n_truth] = truth->get_e();
0314 _b_truthphi[_b_n_truth] = atan2(truth->get_py(),truth->get_px());
0315 _b_trutheta[_b_n_truth] = atanh(truth->get_pz() / _b_truthenergy[_b_n_truth]);
0316
0317
0318 if (_b_trutheta[_b_n_truth] != _b_trutheta[_b_n_truth])
0319 _b_trutheta[_b_n_truth] = -99;
0320 _b_truthpid[_b_n_truth] = truth->get_pid();
0321 _b_truth_trackid[_b_n_truth] = truth->get_track_id();
0322
0323
0324 _b_n_truth++;
0325 if (_b_n_truth >= 50000) break;
0326 }
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353 PHG4TruthInfoContainer::Range second_range = truthinfo->GetSecondaryParticleRange();
0354
0355 for (PHG4TruthInfoContainer::ConstIterator iter = second_range.first; iter != second_range.second; ++iter) {
0356 const PHG4Particle *child = iter->second;
0357 if (child->get_parent_id() > 0) {
0358 PHG4Particle *parent = truthinfo->GetParticle(child->get_parent_id());
0359 if (parent->get_e() > 2.0) {
0360 vertex.insert(child->get_vtx_id());
0361 child_pid[n_child] = child->get_pid();
0362 child_parent_id[n_child] = child->get_parent_id();
0363 child_vertex_id[n_child] = child->get_vtx_id();
0364 child_px[n_child] = child->get_px();
0365 child_py[n_child] = child->get_py();
0366 child_pz[n_child] = child->get_pz();
0367 child_energy[n_child] = child->get_e();
0368 n_child++;
0369 }
0370 }
0371 }
0372
0373
0374 PHG4TruthInfoContainer::VtxRange vtxrange = truthinfo->GetVtxRange();
0375 n_vertex = 0;
0376 for (PHG4TruthInfoContainer::ConstVtxIterator iter = vtxrange.first; iter != vtxrange.second; ++iter) {
0377 PHG4VtxPoint *vtx = iter->second;
0378 if (vertex.find(vtx->get_id()) != vertex.end()) {
0379 vertex_x[n_vertex] = vtx->get_x();
0380 vertex_y[n_vertex] = vtx->get_y();
0381 vertex_z[n_vertex] = vtx->get_z();
0382 vertex_id[n_vertex] = vtx->get_id();
0383 n_vertex++;
0384 }
0385 }
0386
0387 _file->cd();
0388 _tree->Fill();
0389 _ievent++;
0390
0391 return 0;
0392 }
0393
0394 int TruthCaloTree::End(PHCompositeNode *topNode)
0395 {
0396 if (_debug) std::cout<<"Writing File"<<std::endl;
0397 _file->cd();
0398 _file->Write();
0399 _file->Close();
0400
0401 return 0;
0402 }
