Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:13:41 

0001 #include "IsolatedTrackAnalysis.h"
0002 
0003 // Centrality includes
0004 #include <centrality/CentralityInfo.h>
0005 
0006 // Tracking includes
0007 #include <trackbase_historic/SvtxTrack.h>
0008 #include <trackbase_historic/SvtxTrackMap.h>
0009 #include <trackbase_historic/SvtxTrackState.h>
0010 #include <trackbase_historic/SvtxVertexMap.h>
0011 
0012 #include <trackbase/TrkrCluster.h>
0013 #include <trackbase/TrkrClusterv3.h>
0014 #include <trackbase/TrkrClusterv4.h>
0015 #include <trackbase/TrkrHit.h>
0016 #include <trackbase/TrkrHitSetContainer.h>
0017 #include <trackbase/TrkrDefs.h>
0018 #include <trackbase/ActsGeometry.h>
0019 #include <trackbase/ClusterErrorPara.h>
0020 
0021 #include <trackbase/TpcDefs.h>
0022 
0023 #include <trackbase/TrkrClusterContainer.h>
0024 #include <trackbase/TrkrHitSet.h>
0025 #include <trackbase/TrkrClusterHitAssoc.h>
0026 #include <trackbase/TrkrClusterIterationMapv1.h>
0027 
0028 // Cluster includes
0029 #include <calobase/RawCluster.h>
0030 #include <calobase/RawClusterContainer.h>
0031 #include <calobase/RawClusterUtility.h>
0032 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0033 
0034 // Tower includes
0035 #include <calobase/RawTower.h>
0036 #include <calobase/RawTowerContainer.h>
0037 #include <calobase/RawTowerGeom.h>
0038 #include <calobase/RawTowerGeomContainer.h>
0039 #include <calobase/RawTowerDefs.h>
0040 
0041 // G4 truth includes
0042 #include <g4eval/SvtxEvalStack.h>
0043 #include <g4main/PHG4Particle.h>
0044 #include <g4main/PHG4TruthInfoContainer.h>
0045 #include <g4main/PHG4Hit.h>
0046 #include <g4main/PHG4HitContainer.h>
0047 #include <g4main/PHG4VtxPoint.h>
0048 
0049 // HepMC truth includes
0050 #pragma GCC diagnostic push
0051 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0052 #include <HepMC/GenEvent.h>
0053 #pragma GCC diagnostic pop
0054 #include <HepMC/GenVertex.h>
0055 #include <phhepmc/PHHepMCGenEvent.h>
0056 #include <phhepmc/PHHepMCGenEventMap.h>
0057 
0058 // Fun4All includes
0059 #include <fun4all/Fun4AllReturnCodes.h>
0060 #include <phool/PHCompositeNode.h>
0061 #include <phool/getClass.h>
0062 
0063 // Acts includes
0064 #include <Acts/Definitions/Algebra.hpp>
0065 
0066 // ROOT includes
0067 #include <TFile.h>
0068 #include <TTree.h>
0069 
0070 // System includes
0071 #include <iostream>
0072 #include <cmath>
0073 #include <string>
0074 #include <vector>
0075 #include <set>
0076 
0077 //____________________________________________________________________________..
0078 IsolatedTrackAnalysis::IsolatedTrackAnalysis(const std::string &name, const std::string &fileName):
0079   SubsysReco(name),
0080   m_outputFileName(fileName),
0081   m_minEMClusterEnergy(0.1),
0082   m_minIHClusterEnergy(0.012),
0083   m_minOHClusterEnergy(0.025),
0084   m_minCemcTowerEnergy(0.04),
0085   m_minHcalTowerEnergy(0.006),
0086   m_minSimTowerEnergy(0.0),
0087   m_analyzeTracks(true),
0088   m_analyzeClusters(true),
0089   m_analyzeTowers(true),
0090   m_analyzeSimTowers(true),
0091   m_analyzeHepMCTruth(true),
0092   m_analyzeG4Truth(true),
0093   m_analyzeCentrality(true),
0094   m_analyzeAddTruth(true)
0095 {
0096   initializeTrees();
0097 }
0098 
0099 //____________________________________________________________________________..
0100 IsolatedTrackAnalysis::~IsolatedTrackAnalysis()
0101 {
0102   delete m_tracktree;
0103   delete m_clustertree;
0104   delete m_towertree;
0105   delete m_simtowertree;
0106   delete m_hepmctree;
0107   delete m_g4tree;
0108   delete m_centraltree;
0109   delete m_addtruthtree;
0110 }
0111 
0112 //____________________________________________________________________________..
0113 int IsolatedTrackAnalysis::Init(PHCompositeNode *topNode)
0114 {
0115   m_outputFile = new TFile(m_outputFileName.c_str(), "RECREATE");
0116 
0117   counter = 0;
0118 
0119   return Fun4AllReturnCodes::EVENT_OK;
0120 }
0121 
0122 //____________________________________________________________________________..
0123 int IsolatedTrackAnalysis::InitRun(PHCompositeNode *topNode)
0124 {
0125   if(m_analyzeTracks){
0126     RawTowerGeomContainer_Cylinderv1* cemcGeomContainer = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0127     RawTowerGeomContainer_Cylinderv1* ihcalGeomContainer = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0128     RawTowerGeomContainer_Cylinderv1* ohcalGeomContainer = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0129 
0130     if(!cemcGeomContainer){
0131       std::cout << "ERROR: TOWERGEOM_CEMC not found" << std::endl;
0132     }
0133     if(!ihcalGeomContainer){
0134       std::cout << "ERROR: TOWERGEOM_HCALIN not found" << std::endl;
0135     }
0136     if(!ohcalGeomContainer){
0137       std::cout << "ERROR: TOWERGEOM_HCALOUT not found" << std::endl;
0138     }
0139 
0140     m_cemcRadius = cemcGeomContainer->get_radius();
0141     m_ihcalRadius = ihcalGeomContainer->get_radius();
0142     m_ohcalRadius = ohcalGeomContainer->get_radius();
0143   }
0144 
0145   return Fun4AllReturnCodes::EVENT_OK;
0146 }
0147 
0148 //____________________________________________________________________________..
0149 int IsolatedTrackAnalysis::process_event(PHCompositeNode *topNode)
0150 {
0151   if(m_analyzeTracks){ getTracks(topNode); }
0152   if(m_analyzeClusters){ getClusters(topNode); }
0153   if(m_analyzeTowers){ getTowers(topNode); }
0154   if(m_analyzeSimTowers){ getSimTowers(topNode); }
0155   if(m_analyzeHepMCTruth){ getHepMCTruth(topNode); }
0156   if(m_analyzeG4Truth){ getG4Truth(topNode); }
0157   if(m_analyzeCentrality){ getCentrality(topNode); }
0158   if(m_analyzeAddTruth){ getAddTruth(topNode); }
0159 
0160   counter++;
0161 
0162   //if(counter % 100 == 0)
0163     std::cout << counter << " events processed" << std::endl;
0164 
0165   return Fun4AllReturnCodes::EVENT_OK;
0166 }
0167 
0168 //____________________________________________________________________________..
0169 int IsolatedTrackAnalysis::ResetEvent(PHCompositeNode *topNode)
0170 {
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173 
0174 //____________________________________________________________________________..
0175 int IsolatedTrackAnalysis::End(PHCompositeNode *topNode)
0176 {
0177   m_outputFile->cd();
0178   m_tracktree->Write();
0179   m_clustertree->Write();
0180   m_towertree->Write();
0181   m_simtowertree->Write();
0182   m_hepmctree->Write();
0183   m_g4tree->Write();
0184   m_centraltree->Write();
0185   m_addtruthtree->Write();
0186   m_outputFile->Close();
0187   return Fun4AllReturnCodes::EVENT_OK;
0188 }
0189 
0190 
0191 //////////////////////////////////////////////////////////////////////////////////////////////////
0192 
0193 ////////////////////////////////////
0194 // Function extracting centrality //
0195 ////////////////////////////////////
0196 
0197 void IsolatedTrackAnalysis::getCentrality(PHCompositeNode *topNode) 
0198 {
0199   central = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0200   centrality = central->get_centile(CentralityInfo::PROP::bimp);
0201   m_centraltree->Fill();
0202 }
0203 
0204 ////////////////////////////////
0205 // Function extracting tracks //
0206 ////////////////////////////////
0207 
0208 void IsolatedTrackAnalysis::getTracks(PHCompositeNode *topNode)
0209 {
0210   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0211 
0212   // Check whether SvtxTrackMap is present in the node tree or not
0213   if (!trackmap)
0214   {
0215     std::cout << PHWHERE
0216          << "SvtxTrackMap node is missing, can't collect tracks"
0217          << std::endl;
0218     return;
0219   }
0220 
0221   SvtxVertexMap *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0222 
0223   // Check whether SvtxTrackMap is present in the node tree or not
0224   if (!vertexmap)
0225   {
0226     std::cout << PHWHERE
0227          << "SvtxVertexMap node is missing, can't collect tracks"
0228          << std::endl;
0229     return;
0230   }
0231 
0232   // EvalStack for truth track matching
0233   if(!m_svtxEvalStack){ m_svtxEvalStack = new SvtxEvalStack(topNode); }
0234   
0235   m_svtxEvalStack->next_event(topNode);
0236 
0237   // Get the track evaluator
0238   SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
0239 
0240   // Get the range for primary tracks
0241   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0242 
0243   m_trkmult = 0;
0244   // Looping over tracks in an event
0245   for(auto entry : *trackmap)
0246   {
0247     SvtxTrack *track = entry.second;
0248 
0249     m_tr_p[m_trkmult] = track->get_p();
0250     m_tr_pt[m_trkmult] = track->get_pt();
0251     m_tr_eta[m_trkmult] = track->get_eta();
0252     m_tr_phi[m_trkmult] = track->get_phi();
0253   
0254     m_tr_charge[m_trkmult] = track->get_charge();
0255     m_tr_chisq[m_trkmult] = track->get_chisq();
0256     m_tr_ndf[m_trkmult] = track->get_ndf();
0257     m_tr_silicon_hits[m_trkmult] = track->get_silicon_seed()->size_cluster_keys();
0258 
0259     float dca_xy, dca_z, dca_xy_error, dca_z_error;
0260     calculateDCA(track, vertexmap, dca_xy, dca_z, dca_xy_error, dca_z_error);
0261 
0262     m_tr_dca_xy[m_trkmult] = dca_xy;
0263     m_tr_dca_xy_error[m_trkmult] = dca_xy_error;
0264     m_tr_dca_z[m_trkmult] = dca_z;
0265     m_tr_dca_z_error[m_trkmult] = dca_z_error;
0266 
0267     m_tr_x[m_trkmult] = track->get_x();
0268     m_tr_y[m_trkmult] = track->get_y();
0269     m_tr_z[m_trkmult] = track->get_z();
0270 
0271     m_tr_vertex_id[m_trkmult] = track->get_vertex_id();
0272 
0273     // Project tracks to CEMC
0274     if(track->find_state(m_cemcRadius) != track->end_states()){
0275       m_tr_cemc_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_cemcRadius)->second);
0276       m_tr_cemc_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_cemcRadius)->second);
0277     }
0278     else{
0279       m_tr_cemc_eta[m_trkmult] = -999;
0280       m_tr_cemc_phi[m_trkmult] = -999;
0281     }
0282     
0283     // Project tracks to inner HCAL
0284     if(track->find_state(m_ihcalRadius) != track->end_states()){
0285       m_tr_ihcal_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_ihcalRadius)->second);
0286       m_tr_ihcal_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_ihcalRadius)->second);
0287     }
0288     else{
0289       m_tr_ihcal_eta[m_trkmult] = -999;
0290       m_tr_ihcal_phi[m_trkmult] = -999;
0291     }
0292 
0293     // Project tracks to outer HCAL
0294     if(track->find_state(m_ohcalRadius) != track->end_states()){
0295       m_tr_ohcal_eta[m_trkmult] = calculateProjectionEta( track->find_state(m_ohcalRadius)->second);
0296       m_tr_ohcal_phi[m_trkmult] = calculateProjectionPhi( track->find_state(m_ohcalRadius)->second);
0297     }
0298     else{
0299       m_tr_ohcal_eta[m_trkmult] = -999;
0300       m_tr_ohcal_phi[m_trkmult] = -999;
0301     }
0302     
0303     // Matching SvtxTracks to G4 truth
0304     //std::cout << "before truth matching" << std::endl;
0305     PHG4Particle *truthtrack = trackeval->max_truth_particle_by_nclusters(track);
0306     //std::cout << "after truth matching" << std::endl;
0307 
0308     if(truthtrack){
0309       //std::cout << "found truth track" << std::endl;
0310       m_tr_truth_pid[m_trkmult] =truthtrack->get_pid();
0311       m_tr_truth_is_primary[m_trkmult] =truthinfo->is_primary(truthtrack);
0312   
0313       m_tr_truth_e[m_trkmult] =truthtrack->get_e();
0314 
0315       float px = truthtrack->get_px();
0316       float py = truthtrack->get_py();
0317       float pz = truthtrack->get_pz();
0318 
0319       m_tr_truth_pt[m_trkmult] = sqrt(px*px+py*py);
0320       m_tr_truth_phi[m_trkmult] = atan2(py, px);
0321       m_tr_truth_eta[m_trkmult] = atanh(pz/sqrt(px*px+py*py+pz*pz));
0322       m_tr_truth_track_id[m_trkmult] = truthtrack->get_track_id();
0323     }
0324     else{
0325       //std::cout << "truth track null" << std::endl;
0326       m_tr_truth_pid[m_trkmult] = -999;
0327       m_tr_truth_is_primary[m_trkmult] = -999;
0328       m_tr_truth_e[m_trkmult] = -999;
0329       m_tr_truth_pt[m_trkmult] = -999;
0330       m_tr_truth_eta[m_trkmult] = -999;
0331       m_tr_truth_phi[m_trkmult] = -999;
0332       m_tr_truth_track_id[m_trkmult] = -999;
0333     }
0334 
0335     m_trkmult++;
0336   }
0337 
0338   m_vtxmult = 0;
0339   for(auto entry : *vertexmap)
0340   {
0341     SvtxVertex *vertex = entry.second;
0342     
0343     m_vertex_id[m_vtxmult] = vertex->get_id();
0344     m_vx[m_vtxmult] = vertex->get_x();
0345     m_vy[m_vtxmult] = vertex->get_y();
0346     m_vz[m_vtxmult] = vertex->get_z();
0347     m_vtxmult++;
0348   }
0349 
0350   m_tracktree->Fill();
0351 }
0352 
0353 //////////////////////////////////
0354 // Function extracting clusters //
0355 //////////////////////////////////
0356 void IsolatedTrackAnalysis::getClusters(PHCompositeNode *topNode){
0357   
0358   RawClusterContainer *cemcContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC"); 
0359   
0360   // Check whether CLUSTER_CEMC is present in the node tree or not
0361   if (!cemcContainer)
0362   {
0363     std::cout << PHWHERE
0364          << "CLUSTER_CEMC node is missing, can't collect EMCal clusters"
0365          << std::endl;
0366     return;
0367   }
0368 
0369   RawClusterContainer::Map cemcMap = cemcContainer->getClustersMap();
0370   
0371   m_clsmult_cemc = 0;
0372   for(auto entry : cemcMap){
0373     RawCluster* cluster = entry.second;
0374 
0375     if(cluster->get_energy() > m_minEMClusterEnergy){
0376       // calculating eta
0377       // we need the nominal eta, with respect to the origin, not the vertex!
0378       CLHEP::Hep3Vector origin(0, 0, 0);
0379       CLHEP::Hep3Vector cluster_vector = RawClusterUtility::GetECoreVec(*cluster, origin);
0380        
0381       m_cl_cemc_e[m_clsmult_cemc] = cluster->get_energy();
0382       m_cl_cemc_eta[m_clsmult_cemc] = cluster_vector.pseudoRapidity();
0383       m_cl_cemc_phi[m_clsmult_cemc] = cluster->get_phi();
0384       m_cl_cemc_r[m_clsmult_cemc] = cluster->get_r();
0385       m_cl_cemc_z[m_clsmult_cemc] = cluster->get_z();
0386 
0387       m_clsmult_cemc++;  
0388     }
0389   }
0390 
0391   RawClusterContainer *ihcalContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN"); 
0392   
0393    // Check whether CLUSTER_HCALIN is present in the node tree or not
0394   if (!ihcalContainer)
0395   {
0396     std::cout << PHWHERE
0397          << "CLUSTER_HCALIN node is missing, can't collect EMCal clusters"
0398          << std::endl;
0399     return;
0400   }
0401 
0402   RawClusterContainer::Map ihcalMap = ihcalContainer->getClustersMap();
0403   
0404   m_clsmult_ihcal = 0;
0405   for(auto entry : ihcalMap){
0406     RawCluster* cluster = entry.second;
0407 
0408     if(cluster->get_energy() > m_minIHClusterEnergy){
0409       // calculating eta
0410       // we need the nominal eta, with respect to the origin, not the vertex!
0411       CLHEP::Hep3Vector origin(0, 0, 0);
0412       CLHEP::Hep3Vector cluster_vector = RawClusterUtility::GetECoreVec(*cluster, origin);
0413       
0414       m_cl_ihcal_e[m_clsmult_ihcal] = cluster->get_energy();
0415       m_cl_ihcal_eta[m_clsmult_ihcal] = cluster_vector.pseudoRapidity();
0416       m_cl_ihcal_phi[m_clsmult_ihcal] = cluster->get_phi();
0417       m_cl_ihcal_r[m_clsmult_ihcal] = cluster->get_r();
0418       m_cl_ihcal_z[m_clsmult_ihcal] = cluster->get_z();
0419         
0420       m_clsmult_ihcal++;
0421     }
0422   } 
0423   
0424   RawClusterContainer *ohcalContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT");   
0425 
0426    // Check whether CLUSTER_HCALOUT is present in the node tree or not
0427   if (!ohcalContainer)
0428   {
0429     std::cout << PHWHERE
0430          << "CLUSTER_HCALOUT node is missing, can't collect EMCal clusters"
0431          << std::endl;
0432     return;
0433   }
0434 
0435   RawClusterContainer::Map ohcalMap = ohcalContainer->getClustersMap();
0436   
0437   m_clsmult_ohcal = 0;
0438   for(auto entry : ohcalMap){
0439     RawCluster* cluster = entry.second;
0440 
0441     if(cluster->get_energy() > m_minOHClusterEnergy){
0442       // calculating eta
0443       // we need the nominal eta, with respect to the origin, not the vertex!
0444       CLHEP::Hep3Vector origin(0, 0, 0);
0445       CLHEP::Hep3Vector cluster_vector = RawClusterUtility::GetECoreVec(*cluster, origin);
0446  
0447       m_cl_ohcal_e[m_clsmult_ohcal] = cluster->get_energy();
0448       m_cl_ohcal_eta[m_clsmult_ohcal] = cluster_vector.pseudoRapidity();
0449       m_cl_ohcal_phi[m_clsmult_ohcal] = cluster->get_phi();
0450       m_cl_ohcal_r[m_clsmult_ohcal] = cluster->get_r();
0451       m_cl_ohcal_z[m_clsmult_ohcal] = cluster->get_z();
0452     
0453       m_clsmult_ohcal++;  
0454     }   
0455   }
0456   
0457   m_clustertree->Fill();
0458 
0459 }
0460 
0461 ////////////////////////////////
0462 // Function extracting towers //
0463 ////////////////////////////////
0464 
0465 void IsolatedTrackAnalysis::getTowers(PHCompositeNode *topNode) {
0466 
0467   m_twrmult_cemc = 0;
0468 
0469   RawTowerContainer *cemcTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0470   RawTowerGeomContainer *cemcGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0471 
0472   if(cemcTowerContainer && cemcGeom)
0473   {
0474     RawTowerContainer::ConstRange range = cemcTowerContainer->getTowers();
0475     for(RawTowerContainer::ConstIterator it = range.first; it != range.second; it++){
0476       RawTower *tower = it->second;
0477 
0478         RawTowerGeom *tower_geom = cemcGeom->get_tower_geometry(tower->get_key()); 
0479         if(tower->get_energy() > m_minCemcTowerEnergy)
0480         {
0481           m_twr_cemc_e[m_twrmult_cemc] = tower->get_energy();
0482         m_twr_cemc_eta[m_twrmult_cemc] = tower_geom->get_eta();
0483           m_twr_cemc_phi[m_twrmult_cemc] = tower_geom->get_phi();
0484         m_twr_cemc_ieta[m_twrmult_cemc] = tower_geom->get_bineta();
0485           m_twr_cemc_iphi[m_twrmult_cemc] = tower_geom->get_binphi();
0486           m_twrmult_cemc++;
0487         }
0488     }
0489   }
0490 
0491 
0492   m_twrmult_ihcal = 0;
0493 
0494   RawTowerContainer *ihcalTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0495   RawTowerGeomContainer *ihcalGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0496 
0497   if(ihcalTowerContainer && ihcalGeom)
0498   {
0499     RawTowerContainer::ConstRange range = ihcalTowerContainer->getTowers();
0500     for(RawTowerContainer::ConstIterator it = range.first; it != range.second; it++){
0501       RawTower *tower = it->second;
0502 
0503         RawTowerGeom *tower_geom = ihcalGeom->get_tower_geometry(tower->get_key()); 
0504         if(tower->get_energy() > m_minHcalTowerEnergy)
0505         {
0506           m_twr_ihcal_e[m_twrmult_ihcal] = tower->get_energy();
0507         m_twr_ihcal_eta[m_twrmult_ihcal] = tower_geom->get_eta();
0508           m_twr_ihcal_phi[m_twrmult_ihcal] = tower_geom->get_phi();
0509         m_twr_ihcal_ieta[m_twrmult_ihcal] = tower_geom->get_bineta();
0510           m_twr_ihcal_iphi[m_twrmult_ihcal] = tower_geom->get_binphi();
0511           m_twrmult_ihcal++;
0512         }
0513     }
0514   }
0515 
0516 
0517   m_twrmult_ohcal = 0;
0518 
0519   RawTowerContainer *ohcalTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0520   RawTowerGeomContainer *ohcalGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0521 
0522   if(ohcalTowerContainer && ohcalGeom)
0523   {
0524     RawTowerContainer::ConstRange range = ohcalTowerContainer->getTowers();
0525     for(RawTowerContainer::ConstIterator it = range.first; it != range.second; it++){
0526       RawTower *tower = it->second;
0527 
0528         RawTowerGeom *tower_geom = ohcalGeom->get_tower_geometry(tower->get_key()); 
0529         if(tower->get_energy() > m_minHcalTowerEnergy)
0530         {
0531           m_twr_ohcal_e[m_twrmult_ohcal] = tower->get_energy();
0532         m_twr_ohcal_eta[m_twrmult_ohcal] = tower_geom->get_eta();
0533           m_twr_ohcal_phi[m_twrmult_ohcal] = tower_geom->get_phi();
0534         m_twr_ohcal_ieta[m_twrmult_ohcal] = tower_geom->get_bineta();
0535           m_twr_ohcal_iphi[m_twrmult_ohcal] = tower_geom->get_binphi();
0536           m_twrmult_ohcal++;
0537         }
0538     }
0539   }
0540 
0541   m_towertree->Fill();
0542  
0543 }
0544 
0545 ////////////////////////////////
0546 // Function extracting towers //
0547 ////////////////////////////////
0548 
0549 void IsolatedTrackAnalysis::getSimTowers(PHCompositeNode *topNode) {
0550 
0551   m_simtwrmult_cemc = 0;
0552 
0553   RawTowerContainer *cemcTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_SIM_CEMC");
0554   RawTowerGeomContainer *cemcGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0555 
0556   if(cemcTowerContainer && cemcGeom)
0557   {
0558     RawTowerContainer::ConstRange range = cemcTowerContainer->getTowers();
0559     for(RawTowerContainer::ConstIterator it = range.first; it != range.second; it++){
0560       RawTower *tower = it->second;
0561 
0562         RawTowerGeom *tower_geom = cemcGeom->get_tower_geometry(tower->get_key()); 
0563         if(tower->get_energy() > m_minSimTowerEnergy)
0564         {
0565           m_simtwr_cemc_e[m_simtwrmult_cemc] = tower->get_energy();
0566         m_simtwr_cemc_eta[m_simtwrmult_cemc] = tower_geom->get_eta();
0567           m_simtwr_cemc_phi[m_simtwrmult_cemc] = tower_geom->get_phi();
0568         m_simtwr_cemc_ieta[m_simtwrmult_cemc] = tower_geom->get_bineta();
0569           m_simtwr_cemc_iphi[m_simtwrmult_cemc] = tower_geom->get_binphi();
0570           m_simtwrmult_cemc++;
0571         }
0572     }
0573   }
0574 
0575 
0576   m_simtwrmult_ihcal = 0;
0577 
0578   RawTowerContainer *ihcalTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_SIM_HCALIN");
0579   RawTowerGeomContainer *ihcalGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0580 
0581   if(ihcalTowerContainer && ihcalGeom)
0582   {
0583     RawTowerContainer::ConstRange range = ihcalTowerContainer->getTowers();
0584     for(RawTowerContainer::ConstIterator it = range.first; it != range.second; it++){
0585       RawTower *tower = it->second;
0586 
0587         RawTowerGeom *tower_geom = ihcalGeom->get_tower_geometry(tower->get_key()); 
0588         if(tower->get_energy() > m_minSimTowerEnergy)
0589         {
0590           m_simtwr_ihcal_e[m_simtwrmult_ihcal] = tower->get_energy();
0591         m_simtwr_ihcal_eta[m_simtwrmult_ihcal] = tower_geom->get_eta();
0592           m_simtwr_ihcal_phi[m_simtwrmult_ihcal] = tower_geom->get_phi();
0593         m_simtwr_ihcal_ieta[m_simtwrmult_ihcal] = tower_geom->get_bineta();
0594           m_simtwr_ihcal_iphi[m_simtwrmult_ihcal] = tower_geom->get_binphi();
0595           m_simtwrmult_ihcal++;
0596         }
0597     }
0598   }
0599 
0600   m_simtwrmult_ohcal = 0;
0601 
0602   RawTowerContainer *ohcalTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_SIM_HCALOUT");
0603   RawTowerGeomContainer *ohcalGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0604 
0605   if(ohcalTowerContainer && ohcalGeom)
0606   {
0607     RawTowerContainer::ConstRange range = ohcalTowerContainer->getTowers();
0608     for(RawTowerContainer::ConstIterator it = range.first; it != range.second; it++){
0609       RawTower *tower = it->second;
0610 
0611         RawTowerGeom *tower_geom = ohcalGeom->get_tower_geometry(tower->get_key()); 
0612         if(tower->get_energy() > m_minSimTowerEnergy)
0613         {
0614           m_simtwr_ohcal_e[m_simtwrmult_ohcal] = tower->get_energy();
0615         m_simtwr_ohcal_eta[m_simtwrmult_ohcal] = tower_geom->get_eta();
0616           m_simtwr_ohcal_phi[m_simtwrmult_ohcal] = tower_geom->get_phi();
0617         m_simtwr_ohcal_ieta[m_simtwrmult_ohcal] = tower_geom->get_bineta();
0618           m_simtwr_ohcal_iphi[m_simtwrmult_ohcal] = tower_geom->get_binphi();
0619           m_simtwrmult_ohcal++;
0620         }
0621     }
0622   }
0623 
0624   m_simtowertree->Fill(); 
0625 }
0626 
0627 
0628 
0629 ////////////////////////////////////////////
0630 // Functions extracting truth information //
0631 ////////////////////////////////////////////
0632 
0633 void IsolatedTrackAnalysis::getHepMCTruth(PHCompositeNode *topNode) {
0634 
0635   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0636   if (!hepmceventmap) std::cout << PHWHERE << "HEPMC event map node is missing, can't collected HEPMC truth particles"<< std::endl;
0637 
0638   for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
0639          eventIter != hepmceventmap->end(); ++eventIter) {
0640      
0641     /// Get the event
0642     PHHepMCGenEvent *hepmcevent = eventIter->second;
0643       
0644     // To fill TTree, require that the event be the primary event (embedding_id > 0)
0645     if (hepmcevent && hepmcevent->get_embedding_id() == 0)
0646     {
0647       /// Get the event characteristics, inherited from HepMC classes
0648       HepMC::GenEvent *truthevent = hepmcevent->getEvent();
0649       if (!truthevent) {
0650         std::cout << PHWHERE
0651              << "no evt pointer under phhepmvgeneventmap found "
0652              << std::endl;
0653       }
0654 
0655       /// Loop over all the truth particles and get their information
0656       m_hepmc = 0;
0657       for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
0658            iter != truthevent->particles_end(); ++iter) {
0659 
0660         m_hepmc_e[m_hepmc] = (*iter)->momentum().e();
0661         m_hepmc_pid[m_hepmc] = (*iter)->pdg_id();
0662         m_hepmc_eta[m_hepmc] = (*iter)->momentum().pseudoRapidity();
0663         m_hepmc_phi[m_hepmc] = (*iter)->momentum().phi();
0664         m_hepmc_pt[m_hepmc] = sqrt((*iter)->momentum().px() * (*iter)->momentum().px() 
0665                             + (*iter)->momentum().py() * (*iter)->momentum().py());
0666         /// Fill the truth tree
0667         m_hepmc++;
0668         if (m_hepmc >= 20000) break;
0669       }
0670       break;
0671     } 
0672   }
0673 
0674   m_hepmctree->Fill();
0675 
0676 }
0677 
0678 void IsolatedTrackAnalysis::getG4Truth(PHCompositeNode *topNode) {
0679 
0680   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0681   if (!truthinfo) std::cout << PHWHERE << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"<< std::endl;
0682 
0683   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0684   m_g4 = 0;
0685   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0686        iter != range.second; ++iter) {
0687 
0688     // Get truth particle
0689     const PHG4Particle *truth = iter->second;
0690     if (!truthinfo->is_primary(truth)) continue;
0691 
0692     /// Get this particles momentum, etc.
0693     m_g4_pt[m_g4] = sqrt(truth->get_px() * truth->get_px()
0694                             + truth->get_py() * truth->get_py());
0695     m_g4_e[m_g4] = truth->get_e();
0696     m_g4_phi[m_g4] = atan2(truth->get_py(), truth->get_px());
0697     m_g4_eta[m_g4] = atanh(truth->get_pz() / truth->get_e());
0698     if (m_g4_eta[m_g4] != m_g4_eta[m_g4]) m_g4_eta[m_g4] = -999; // check for nans
0699     m_g4_pid[m_g4] = truth->get_pid();
0700     m_g4_track_id[m_g4] = truth->get_track_id();
0701     m_g4_parent_id[m_g4] = truth->get_parent_id();
0702     m_g4++;
0703     if (m_g4 >= 20000) break;
0704   }
0705 
0706   PHG4TruthInfoContainer::Range second_range = truthinfo->GetSecondaryParticleRange();
0707 
0708   for (PHG4TruthInfoContainer::ConstIterator siter = second_range.first;
0709         siter != second_range.second; ++siter) {
0710     if (m_g4 >= 19999) break;
0711     // Get photons from pi0 decays 
0712     const PHG4Particle *truth = siter->second;
0713     if (truth->get_pid() == 22) {
0714       PHG4Particle *parent = truthinfo->GetParticle(truth->get_parent_id());
0715       if (parent->get_pid() == 111 && parent->get_track_id() > 0) {
0716         m_g4_pt[m_g4] = sqrt(truth->get_px() * truth->get_px()
0717                           + truth->get_py() * truth->get_py());
0718         m_g4_e[m_g4] = truth->get_e();
0719         m_g4_phi[m_g4] = atan2(truth->get_py(), truth->get_px());
0720         m_g4_eta[m_g4] = atanh(truth->get_pz() / truth->get_e());
0721         if (m_g4_eta[m_g4] != m_g4_eta[m_g4]) m_g4_eta[m_g4] = -999; // check for nans
0722         m_g4_pid[m_g4] = truth->get_pid();
0723         m_g4_track_id[m_g4] = truth->get_track_id();
0724         m_g4_parent_id[m_g4] = truth->get_parent_id();
0725         m_g4++;
0726       }
0727     }
0728   }
0729 
0730   m_g4tree->Fill();
0731 
0732 }
0733 
0734 ///////////////////////////////////////////////////////
0735 // Functions extracting additional truth information //
0736 ///////////////////////////////////////////////////////
0737 
0738 void IsolatedTrackAnalysis::getAddTruth(PHCompositeNode *topNode) {
0739 
0740   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0741   if (!truthinfo) std::cout << PHWHERE << "PHG4TruthInfoContainer node is missing, can't collect additional truth particles"<< std::endl;
0742   PHG4TruthInfoContainer::Range second_range = truthinfo->GetSecondaryParticleRange();
0743 
0744   n_child = 0;
0745   std::set<int> vertex;
0746   if (!vertex.empty()) vertex.clear();
0747   for (PHG4TruthInfoContainer::ConstIterator iter = second_range.first; iter != second_range.second; ++iter) {
0748     const PHG4Particle *child = iter->second;
0749     if (child->get_parent_id() > 0) {
0750       vertex.insert(child->get_vtx_id());
0751       child_pid[n_child] = child->get_pid();
0752       child_parent_id[n_child] = child->get_parent_id();
0753       child_vertex_id[n_child] = child->get_vtx_id();
0754       child_px[n_child] = child->get_px();
0755       child_py[n_child] = child->get_py();
0756       child_pz[n_child] = child->get_pz();
0757       child_energy[n_child] = child->get_e();
0758       n_child++;
0759       if (n_child >= 100000) break;
0760     }
0761   }
0762   
0763   PHG4TruthInfoContainer::VtxRange vtxrange = truthinfo->GetVtxRange();
0764   n_vertex = 0;
0765   for (PHG4TruthInfoContainer::ConstVtxIterator iter = vtxrange.first; iter != vtxrange.second; ++iter) {
0766     PHG4VtxPoint *vtx = iter->second;
0767     if (vertex.find(vtx->get_id()) != vertex.end()) {
0768       vertex_x[n_vertex] = vtx->get_x();
0769       vertex_y[n_vertex] = vtx->get_y();
0770       vertex_z[n_vertex] = vtx->get_z();
0771       vertex_id[n_vertex] = vtx->get_id();
0772       n_vertex++;
0773       if (n_vertex >= 100000) break;
0774     }
0775   }
0776 
0777   PHG4HitContainer* blackhole = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BH_1");
0778   if (!blackhole) std::cout << "No blackhole" << std::endl;
0779 
0780   _nBH = 0;
0781   PHG4HitContainer::ConstRange bh_hit_range = blackhole->getHits();
0782   for (PHG4HitContainer::ConstIterator hit_iter = bh_hit_range.first;
0783          hit_iter != bh_hit_range.second; hit_iter++) {
0784     PHG4Hit *this_hit = hit_iter->second;
0785     assert(this_hit);
0786     _BH_e[_nBH] = this_hit->get_edep();
0787     _BH_px[_nBH] = this_hit->get_px(0);
0788     _BH_py[_nBH] = this_hit->get_py(0);
0789     _BH_pz[_nBH] = this_hit->get_pz(0);
0790     PHG4Particle *p = truthinfo->GetParticle(this_hit->get_trkid());
0791     _BH_track_id[_nBH] = p->get_track_id();
0792     _nBH++;
0793   }
0794 
0795   m_addtruthtree->Fill();
0796 }
0797 
0798 ////////////////////////////
0799 // Track helper functions //
0800 ////////////////////////////
0801 
0802 float IsolatedTrackAnalysis::calculateProjectionEta(SvtxTrackState* projectedState){
0803   float x = projectedState->get_x();// - initialState->get_x();
0804   float y = projectedState->get_y();// - initialState->get_y();
0805   float z = projectedState->get_z();// - initialState->get_z();
0806 
0807   float theta = acos(z / sqrt(x*x + y*y + z*z) );
0808 
0809   return -log( tan(theta/2.0) );
0810 }
0811 
0812 float IsolatedTrackAnalysis::calculateProjectionPhi(SvtxTrackState* projectedState){
0813   float x = projectedState->get_x();// - initialState->get_x();
0814   float y = projectedState->get_y();// - initialState->get_y();
0815  
0816   return atan2(y,x);
0817 }
0818 
0819 // From SvtxEval.cc
0820 void IsolatedTrackAnalysis::calculateDCA(SvtxTrack* track, SvtxVertexMap* vertexmap,
0821           float& dca3dxy, float& dca3dz,
0822           float& dca3dxysigma, float& dca3dzsigma)
0823 {
0824   Acts::Vector3 pos(track->get_x(),
0825         track->get_y(),
0826         track->get_z());
0827   Acts::Vector3 mom(track->get_px(),
0828         track->get_py(),
0829         track->get_pz());
0830 
0831   auto vtxid = track->get_vertex_id();
0832   auto svtxVertex = vertexmap->get(vtxid);
0833   if(!svtxVertex)
0834     { return; }
0835   Acts::Vector3 vertex(svtxVertex->get_x(),
0836            svtxVertex->get_y(),
0837            svtxVertex->get_z());
0838 
0839   pos -= vertex;
0840 
0841   Acts::ActsSymMatrix<3> posCov;
0842   for(int i = 0; i < 3; ++i)
0843     {
0844       for(int j = 0; j < 3; ++j)
0845   {
0846     posCov(i, j) = track->get_error(i, j);
0847   } 
0848     }
0849   
0850   Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
0851   float phi = atan2(r(1), r(0));
0852   
0853   Acts::RotationMatrix3 rot;
0854   Acts::RotationMatrix3 rot_T;
0855   rot(0,0) = cos(phi);
0856   rot(0,1) = -sin(phi);
0857   rot(0,2) = 0;
0858   rot(1,0) = sin(phi);
0859   rot(1,1) = cos(phi);
0860   rot(1,2) = 0;
0861   rot(2,0) = 0;
0862   rot(2,1) = 0;
0863   rot(2,2) = 1;
0864   
0865   rot_T = rot.transpose();
0866 
0867   Acts::Vector3 pos_R = rot * pos;
0868   Acts::ActsSymMatrix<3> rotCov = rot * posCov * rot_T;
0869 
0870   dca3dxy = pos_R(0);
0871   dca3dz = pos_R(2);
0872   dca3dxysigma = sqrt(rotCov(0,0));
0873   dca3dzsigma = sqrt(rotCov(2,2));
0874   
0875 }
0876 
0877 
0878 ///////////////////////////////////
0879 // Function to set tree branches //
0880 ///////////////////////////////////
0881 
0882 void IsolatedTrackAnalysis::initializeTrees()
0883 {
0884   //////////////////
0885   // Central tree //
0886   //////////////////
0887   m_centraltree = new TTree("centraltree", "Tree of centralities");
0888   m_centraltree->Branch("centrality", &centrality, "centrality/F");
0889   
0890   ////////////////
0891   // Track tree //
0892   ////////////////
0893   m_tracktree = new TTree("tracktree", "Tree of svtx tracks");
0894   
0895   m_tracktree->Branch("m_trkmult", &m_trkmult, "m_trkmult/I");
0896 
0897   // Basic properties
0898   m_tracktree->Branch("m_tr_p",  m_tr_p, "m_tr_p[m_trkmult]/F");
0899   m_tracktree->Branch("m_tr_pt",  m_tr_pt, "m_tr_pt[m_trkmult]/F");
0900   m_tracktree->Branch("m_tr_eta", m_tr_eta, "m_tr_eta[m_trkmult]/F");
0901   m_tracktree->Branch("m_tr_phi", m_tr_phi, "m_tr_phi[m_trkmult]/F");
0902   m_tracktree->Branch("m_tr_charge", m_tr_charge, "m_tr_charge[m_trkmult]/I");
0903   m_tracktree->Branch("m_tr_chisq",  m_tr_chisq, "m_tr_chisq[m_trkmult]/F");
0904   m_tracktree->Branch("m_tr_ndf",    m_tr_ndf, "m_tr_ndf[m_trkmult]/I");
0905   m_tracktree->Branch("m_tr_silicon_hits",    m_tr_silicon_hits, "m_tr_silicon_hits[m_trkmult]/I");
0906  
0907   // Distance of closest approach
0908   m_tracktree->Branch("m_tr_dca_xy", m_tr_dca_xy, "m_tr_dca_xy[m_trkmult]/F");
0909   m_tracktree->Branch("m_tr_dca_xy_error", m_tr_dca_xy_error, "m_tr_dc_xy_error[m_trkmult]/F");
0910   m_tracktree->Branch("m_tr_dca_z",  m_tr_dca_z, "m_tr_dca_z[m_trkmult]/F");
0911   m_tracktree->Branch("m_tr_dca_z_error", m_tr_dca_z_error, "m_tr_dca_z_error[m_trkmult]/F");
0912 
0913   // Initial point of track
0914   m_tracktree->Branch("m_tr_x", m_tr_x, "m_tr_x[m_trkmult]/F");
0915   m_tracktree->Branch("m_tr_y", m_tr_y, "m_tr_y[m_trkmult]/F");
0916   m_tracktree->Branch("m_tr_z", m_tr_z, "m_tr_z[m_trkmult]/F");
0917 
0918   // Vertex id of track
0919   m_tracktree->Branch("m_tr_vertex_id", m_tr_vertex_id, "m_tr_vertex_id[m_trkmult]/I");
0920   
0921   // Vertex ids and positions, also stored on track tree
0922   m_tracktree->Branch("m_vtxmult", &m_vtxmult, "m_vtxmult/I");
0923   m_tracktree->Branch("m_vertex_id", m_vertex_id, "m_vertex_id[m_vtxmult]/I");
0924   m_tracktree->Branch("m_vx", m_vx, "m_vx[m_vtxmult]/F");
0925   m_tracktree->Branch("m_vy", m_vy, "m_vy[m_vtxmult]/F");
0926   m_tracktree->Branch("m_vz", m_vz, "m_vz[m_vtxmult]/F");
0927 
0928   // Projection of track to calorimeters
0929   // CEMC
0930   m_tracktree->Branch("m_tr_cemc_eta", m_tr_cemc_eta, "m_tr_cemc_eta[m_trkmult]/F");
0931   m_tracktree->Branch("m_tr_cemc_phi", m_tr_cemc_phi, "m_tr_cemc_phi[m_trkmult]/F");
0932   // Inner HCAL
0933   m_tracktree->Branch("m_tr_ihcal_eta", m_tr_ihcal_eta, "m_tr_ihcal_eta[m_trkmult]/F");
0934   m_tracktree->Branch("m_tr_ihcal_phi", m_tr_ihcal_phi, "m_tr_ihcal_phi[m_trkmult]/F");
0935   /// Outer HCAL
0936   m_tracktree->Branch("m_tr_ohcal_eta", m_tr_ohcal_eta, "m_tr_ohcal_eta[m_trkmult]/F");
0937   m_tracktree->Branch("m_tr_ohcal_phi", m_tr_ohcal_phi, "m_tr_ohcal_phi[m_trkmult]/F");
0938 
0939   // Matched truth track
0940   m_tracktree->Branch("m_tr_truth_pid", m_tr_truth_pid, "m_tr_truth_pid[m_trkmult]/I");
0941   m_tracktree->Branch("m_tr_truth_is_primary",m_tr_truth_is_primary,"m_tr_truth_is_primary[m_trkmult]/I");
0942   m_tracktree->Branch("m_tr_truth_e", m_tr_truth_e, "m_tr_truth_e[m_trkmult]/F");
0943   m_tracktree->Branch("m_tr_truth_pt", m_tr_truth_pt, "m_tr_truth_pt[m_trkmult]/F");
0944   m_tracktree->Branch("m_tr_truth_eta", m_tr_truth_eta, "m_tr_truth_eta[m_trkmult]/F");
0945   m_tracktree->Branch("m_tr_truth_phi", m_tr_truth_phi, "m_tr_truth_phi[m_trkmult]/F");
0946   m_tracktree->Branch("m_tr_truth_track_id", m_tr_truth_track_id, "m_tr_truth_track_id[m_trkmult]/I");
0947 
0948   //////////////////
0949   // Cluster tree //
0950   //////////////////
0951   m_clustertree = new TTree("clustertree", "Tree of raw clusters");
0952  
0953   // CEMC clusters
0954   m_clustertree->Branch("m_clsmult_cemc", &m_clsmult_cemc, "m_clsmult_cemc/I");
0955   m_clustertree->Branch("m_cl_cemc_e", m_cl_cemc_e, "m_cl_cemc_e[m_clsmult_cemc]/F");
0956   m_clustertree->Branch("m_cl_cemc_eta", m_cl_cemc_eta, "m_cl_cemc_eta[m_clsmult_cemc]/F");
0957   m_clustertree->Branch("m_cl_cemc_phi", m_cl_cemc_phi, "m_cl_cemc_phi[m_clsmult_cemc]/F");
0958   m_clustertree->Branch("m_cl_cemc_r", m_cl_cemc_r, "m_cl_cemc_r[m_clsmult_cemc]/F");
0959   m_clustertree->Branch("m_cl_cemc_z", m_cl_cemc_z, "m_cl_cemc_z[m_clsmult_cemc]/F");
0960 
0961   // Inner HCAL clusters
0962   m_clustertree->Branch("m_clsmult_ihcal", &m_clsmult_ihcal, "m_clsmult_ihcal/I");
0963   m_clustertree->Branch("m_cl_ihcal_e", m_cl_ihcal_e, "m_cl_ihcal_e[m_clsmult_ihcal]/F");
0964   m_clustertree->Branch("m_cl_ihcal_eta", m_cl_ihcal_eta, "m_cl_ihcal_eta[m_clsmult_ihcal]/F");
0965   m_clustertree->Branch("m_cl_ihcal_phi", m_cl_ihcal_phi, "m_cl_ihcal_phi[m_clsmult_ihcal]/F");
0966   m_clustertree->Branch("m_cl_ihcal_r", m_cl_ihcal_r, "m_cl_ihcal_r[m_clsmult_ihcal]/F");
0967   m_clustertree->Branch("m_cl_ihcal_z", m_cl_ihcal_z, "m_cl_ihcal_z[m_clsmult_ihcal]/F");
0968 
0969   // Outer HCAL clusters
0970   m_clustertree->Branch("m_clsmult_ohcal", &m_clsmult_ohcal, "m_clsmult_ohcal/I");
0971   m_clustertree->Branch("m_cl_ohcal_e", m_cl_ohcal_e, "m_cl_ohcal_e[m_clsmult_ohcal]/F");
0972   m_clustertree->Branch("m_cl_ohcal_eta", m_cl_ohcal_eta, "m_cl_ohcal_eta[m_clsmult_ohcal]/F");
0973   m_clustertree->Branch("m_cl_ohcal_phi", m_cl_ohcal_phi, "m_cl_ohcal_phi[m_clsmult_ohcal]/F");
0974   m_clustertree->Branch("m_cl_ohcal_r", m_cl_ohcal_r, "m_cl_ohcal_r[m_clsmult_ohcal]/F");
0975   m_clustertree->Branch("m_cl_ohcal_z", m_cl_ohcal_z, "m_cl_ohcal_z[m_clsmult_ohcal]/F");
0976 
0977   ////////////////
0978   // Tower tree //
0979   ////////////////
0980   m_towertree = new TTree("towertree", "Tree of raw towers");
0981 
0982   // CEMC towers
0983   m_towertree->Branch("m_twrmult_cemc", &m_twrmult_cemc, "m_twrmult_cemc/I");
0984   m_towertree->Branch("m_twr_cemc_e", m_twr_cemc_e, "m_twr_cemc_e[m_twrmult_cemc]/F");
0985   m_towertree->Branch("m_twr_cemc_eta", m_twr_cemc_eta, "m_twr_cemc_eta[m_twrmult_cemc]/F");
0986   m_towertree->Branch("m_twr_cemc_phi", m_twr_cemc_phi, "m_twr_cemc_phi[m_twrmult_cemc]/F");
0987   m_towertree->Branch("m_twr_cemc_ieta", m_twr_cemc_ieta, "m_twr_cemc_ieta[m_twrmult_cemc]/I");
0988   m_towertree->Branch("m_twr_cemc_iphi", m_twr_cemc_iphi, "m_twr_cemc_iphi[m_twrmult_cemc]/I");
0989 
0990   // Inner HCAL towers
0991   m_towertree->Branch("m_twrmult_ihcal", &m_twrmult_ihcal, "m_twrmult_ihcal/I");
0992   m_towertree->Branch("m_twr_ihcal_e", m_twr_ihcal_e, "m_twr_ihcal_e[m_twrmult_ihcal]/F");
0993   m_towertree->Branch("m_twr_ihcal_eta", m_twr_ihcal_eta, "m_twr_ihcal_eta[m_twrmult_ihcal]/F");
0994   m_towertree->Branch("m_twr_ihcal_phi", m_twr_ihcal_phi, "m_twr_ihcal_phi[m_twrmult_ihcal]/F");
0995   m_towertree->Branch("m_twr_ihcal_ieta", m_twr_ihcal_ieta, "m_twr_ihcal_ieta[m_twrmult_ihcal]/I");
0996   m_towertree->Branch("m_twr_ihcal_iphi", m_twr_ihcal_iphi, "m_twr_ihcal_iphi[m_twrmult_ihcal]/I");
0997 
0998   // Outer HCAL towers
0999   m_towertree->Branch("m_twrmult_ohcal", &m_twrmult_ohcal, "m_twrmult_ohcal/I");
1000   m_towertree->Branch("m_twr_ohcal_e", m_twr_ohcal_e, "m_twr_ohcal_e[m_twrmult_ohcal]/F");
1001   m_towertree->Branch("m_twr_ohcal_eta", m_twr_ohcal_eta, "m_twr_ohcal_eta[m_twrmult_ohcal]/F");
1002   m_towertree->Branch("m_twr_ohcal_phi", m_twr_ohcal_phi, "m_twr_ohcal_phi[m_twrmult_ohcal]/F");
1003   m_towertree->Branch("m_twr_ohcal_ieta", m_twr_ohcal_ieta, "m_twr_ohcal_ieta[m_twrmult_ohcal]/I");
1004   m_towertree->Branch("m_twr_ohcal_iphi", m_twr_ohcal_iphi, "m_twr_ohcal_iphi[m_twrmult_ohcal]/I");
1005 
1006   ////////////////////
1007   // Sim tower tree //
1008   ////////////////////
1009   m_simtowertree = new TTree("simtowertree", "Tree of raw simtowers");
1010 
1011   // CEMC simtowers
1012   m_simtowertree->Branch("m_simtwrmult_cemc", &m_simtwrmult_cemc, "m_simtwrmult_cemc/I");
1013   m_simtowertree->Branch("m_simtwr_cemc_e", m_simtwr_cemc_e, "m_simtwr_cemc_e[m_simtwrmult_cemc]/F");
1014   m_simtowertree->Branch("m_simtwr_cemc_eta", m_simtwr_cemc_eta, "m_simtwr_cemc_eta[m_simtwrmult_cemc]/F");
1015   m_simtowertree->Branch("m_simtwr_cemc_phi", m_simtwr_cemc_phi, "m_simtwr_cemc_phi[m_simtwrmult_cemc]/F");
1016   m_simtowertree->Branch("m_simtwr_cemc_ieta", m_simtwr_cemc_ieta, "m_simtwr_cemc_ieta[m_simtwrmult_cemc]/I");
1017   m_simtowertree->Branch("m_simtwr_cemc_iphi", m_simtwr_cemc_iphi, "m_simtwr_cemc_iphi[m_simtwrmult_cemc]/I");
1018 
1019   // Inner HCAL simtowers
1020   m_simtowertree->Branch("m_simtwrmult_ihcal", &m_simtwrmult_ihcal, "m_simtwrmult_ihcal/I");
1021   m_simtowertree->Branch("m_simtwr_ihcal_e", m_simtwr_ihcal_e, "m_simtwr_ihcal_e[m_simtwrmult_ihcal]/F");
1022   m_simtowertree->Branch("m_simtwr_ihcal_eta", m_simtwr_ihcal_eta, "m_simtwr_ihcal_eta[m_simtwrmult_ihcal]/F");
1023   m_simtowertree->Branch("m_simtwr_ihcal_phi", m_simtwr_ihcal_phi, "m_simtwr_ihcal_phi[m_simtwrmult_ihcal]/F");
1024   m_simtowertree->Branch("m_simtwr_ihcal_ieta", m_simtwr_ihcal_ieta, "m_simtwr_ihcal_ieta[m_simtwrmult_ihcal]/I");
1025   m_simtowertree->Branch("m_simtwr_ihcal_iphi", m_simtwr_ihcal_iphi, "m_simtwr_ihcal_iphi[m_simtwrmult_ihcal]/I");
1026 
1027   // Outer HCAL simtowers
1028   m_simtowertree->Branch("m_simtwrmult_ohcal", &m_simtwrmult_ohcal, "m_simtwrmult_ohcal/I");
1029   m_simtowertree->Branch("m_simtwr_ohcal_e", m_simtwr_ohcal_e, "m_simtwr_ohcal_e[m_simtwrmult_ohcal]/F");
1030   m_simtowertree->Branch("m_simtwr_ohcal_eta", m_simtwr_ohcal_eta, "m_simtwr_ohcal_eta[m_simtwrmult_ohcal]/F");
1031   m_simtowertree->Branch("m_simtwr_ohcal_phi", m_simtwr_ohcal_phi, "m_simtwr_ohcal_phi[m_simtwrmult_ohcal]/F");
1032   m_simtowertree->Branch("m_simtwr_ohcal_ieta", m_simtwr_ohcal_ieta, "m_simtwr_ohcal_ieta[m_simtwrmult_ohcal]/I");
1033   m_simtowertree->Branch("m_simtwr_ohcal_iphi", m_simtwr_ohcal_iphi, "m_simtwr_ohcal_iphi[m_simtwrmult_ohcal]/I");
1034 
1035   /////////////////////////
1036   // HepMC particle tree //
1037   /////////////////////////
1038 
1039   m_hepmctree = new TTree("hepmctree","Tree of truth hepmc info and particles");
1040   m_hepmctree->Branch("m_hepmc",&m_hepmc, "m_hepmc/I");
1041   m_hepmctree->Branch("m_hepmc_pid", m_hepmc_pid, "m_hepmc_pid[m_hepmc]/I");
1042   m_hepmctree->Branch("m_hepmc_e", m_hepmc_e, "m_hepmc_e[m_hepmc]/F");
1043   m_hepmctree->Branch("m_hepmc_pt", m_hepmc_pt, "m_hepmc_pt[m_hepmc]/F");
1044   m_hepmctree->Branch("m_hepmc_eta", m_hepmc_eta, "m_hepmc_eta[m_hepmc]/F");
1045   m_hepmctree->Branch("m_hepmc_phi", m_hepmc_phi, "m_hepmc_phi[m_hepmc]/F");
1046 
1047   //////////////////////
1048   // G4 particle tree //
1049   //////////////////////
1050 
1051   m_g4tree = new TTree("g4tree","Tree of truth G4 particles");
1052   m_g4tree->Branch("m_g4",&m_g4, "m_g4/I");
1053   m_g4tree->Branch("m_g4_pid", m_g4_pid, "m_g4_pid[m_g4]/I");
1054   m_g4tree->Branch("m_g4_e", m_g4_e, "m_g4_e[m_g4]/F");
1055   m_g4tree->Branch("m_g4_pt", m_g4_pt, "m_g4_pt[m_g4]/F");
1056   m_g4tree->Branch("m_g4_eta", m_g4_eta, "m_g4_eta[m_g4]/F");
1057   m_g4tree->Branch("m_g4_phi", m_g4_phi, "m_g4_phi[m_g4]/F");
1058   m_g4tree->Branch("m_g4_track_id", m_g4_track_id, "m_g4_track_id[m_g4]/I");
1059   m_g4tree->Branch("m_g4_parent_id", m_g4_parent_id, "m_g4_parent_id[m_g4]/I");
1060 
1061   ///////////////////////////
1062   // Additional truth tree //
1063   ///////////////////////////
1064 
1065   m_addtruthtree = new TTree("addtruthtree", "Tree of additional truth information");
1066   
1067   // child information
1068   m_addtruthtree->Branch("n_child",&n_child,"n_child/I");
1069   m_addtruthtree->Branch("child_pid",child_pid,"child_pid[n_child]/I");
1070   m_addtruthtree->Branch("child_parent_id",child_parent_id,"child_parent_id[n_child]/I");
1071   m_addtruthtree->Branch("child_vertex_id",child_vertex_id,"child_vertex_id[n_child]/I");
1072   m_addtruthtree->Branch("child_px",child_px,"child_px[n_child]/F");
1073   m_addtruthtree->Branch("child_py",child_py,"child_py[n_child]/F");
1074   m_addtruthtree->Branch("child_pz",child_pz,"child_pz[n_child]/F");
1075   m_addtruthtree->Branch("child_energy",child_energy,"child_energy[n_child]/F");
1076 
1077   // vertex information
1078   m_addtruthtree->Branch("n_vertex",&n_vertex,"n_vertex/I");
1079   m_addtruthtree->Branch("vertex_id",vertex_id,"vertex_id[n_vertex]/I");
1080   m_addtruthtree->Branch("vertex_x",vertex_x,"vertex_x[n_vertex]/F");
1081   m_addtruthtree->Branch("vertex_y",vertex_y,"vertex_y[n_vertex]/F");
1082   m_addtruthtree->Branch("vertex_z",vertex_z,"vertex_z[n_vertex]/F");
1083 
1084   // blackhole information
1085   m_addtruthtree->Branch("_nBH",&_nBH,"_nBH/I");
1086   m_addtruthtree->Branch("BH_track_id",_BH_track_id,"BH_track_id[_nBH]/I");
1087   m_addtruthtree->Branch("BH_px",_BH_px,"BH_px[_nBH]/F");
1088   m_addtruthtree->Branch("BH_py",_BH_py,"BH_py[_nBH]/F");
1089   m_addtruthtree->Branch("BH_pz",_BH_pz,"BH_pz[_nBH]/F");
1090   m_addtruthtree->Branch("BH_e",_BH_e,"BH_e[_nBH]/F");
1091 
1092   
1093 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!