CaloDebugTools/HCalDebugTreeMaker/FlatHCalDebugTreeMakerProcessor.cc

0001 // ----------------------------------------------------------------------------
0002 // 'FlatDebugHCalTreeMaker.cc'
0003 // Derek Anderson
0004 // 03.09.2024
0005 //
0006 // A JANA plugin to construct a flat tree of HCal
0007 // cells and clusters, generated particles, and
0008 // all MC particles.
0009 // ----------------------------------------------------------------------------
0010
0011 #define FLATHCALDEBUGTREEMAKERPROCESSOR_CC
0012
0013 // plugin definition
0014 #include "FlatHCalDebugTreeMakerProcessor.h"
0015
0016 // The following just makes this a JANA plugin
0017 extern "C" {
0018     void InitPlugin(JApplication *app) {
0019         InitJANAPlugin(app);
0020         app->Add(new FlatHCalDebugTreeMakerProcessor);
0021     }
0022 }
0023
0024
0025
0026 // public methods -------------------------------------------------------------
0027
0028 void FlatHCalDebugTreeMakerProcessor::InitWithGlobalRootLock(){
0029
0030   //  grab output file
0031   auto rootfile_svc = GetApplication() -> GetService<RootFile_service>();
0032   auto rootfile     = rootfile_svc     -> GetHistFile();
0033
0034   // create directory
0035   m_pluginDir = rootfile -> mkdir(m_config.sPlugin.data());
0036   m_pluginDir -> cd();
0037
0038   // reset tree variables
0039   ResetVariables();
0040
0041   // initialize event index
0042   m_evtIndex = 0;
0043
0044   // initialize output trees
0045   InitializeDecoder();
0046   InitializeTree();
0047   return;
0048
0049 }  // end 'InitWithGlobalLock()'
0050
0051
0052
0053 void FlatHCalDebugTreeMakerProcessor::ProcessSequential(const std::shared_ptr<const JEvent>& event) {
0054
0055   // reset tree variables
0056   ResetVariables();
0057
0058   // loop over generated particles
0059   size_t nGenPar = 0;
0060   for (auto gen : m_genPars()) {
0061
0062     // only accept truth particles
0063     const int  typeGen = gen -> getType();
0064     const bool isTruth = (typeGen == 1);
0065     if (!isTruth) continue;
0066
0067     FillGenVariables( gen );
0068     ++nGenPar;
0069
0070   }  // end generated particle loop
0071
0072   // loop over mc particles
0073   size_t nMCPar = 0;
0074   for (auto mc : m_mcPars()) {
0075
0076     FillMCVariables( mc );
0077     ++nMCPar;
0078
0079   }  // end mc particle loop
0080
0081   // loop over bhcal hits
0082   size_t nCell = 0;
0083   for (auto cell : m_recHits()) {
0084
0085     FillCellVariables( cell );
0086     ++nCell;
0087
0088   }  // end bhcal hit loop
0089
0090   // loop over clusters
0091   int64_t nClust = 0;
0092   for (auto cluster : m_clusters()) {
0093
0094     FillClusterVariables( cluster, nClust );
0095     ++nClust;
0096
0097   }  // end bhcal cluster loop
0098
0099   // set event variables
0100   m_nGen   = nGenPar;
0101   m_nMC    = nMCPar;
0102   m_nCells = nCell;
0103   m_nClust = nClust;
0104
0105   // fill output tree
0106   m_outTree -> Fill();
0107
0108   // increment event index and exit routine
0109   ++m_evtIndex;
0110   return;
0111
0112 }  // end 'ProcessSequential(std::shared_ptr<JEvent>&)'
0113
0114
0115
0116 void FlatHCalDebugTreeMakerProcessor::FinishWithGlobalRootLock() {
0117
0118   // clean up variables
0119   ResetVariables();
0120   return;
0121
0122 }  // end 'FinishWithGlobalRootLock()'
0123
0124
0125
0126 // internal methods -----------------------------------------------------------
0127
0128 void FlatHCalDebugTreeMakerProcessor::InitializeDecoder() {
0129
0130   // grab detector
0131   auto detector = GetApplication() -> GetService<DD4hep_service>() -> detector();
0132
0133   // make sure readout is available
0134   dd4hep::IDDescriptor descriptor;
0135   try {
0136     descriptor  = detector -> readout(m_config.sSimHits.data()).idSpec();
0137   } catch (const std::runtime_error &err) {
0138     throw std::runtime_error("PANIC: readout class is not in output!");
0139   }
0140
0141   // grab decoder and test
0142   std::vector<short> indices(m_config.sFields.size());
0143   try {
0144     m_decoder = descriptor.decoder();
0145     for (size_t iField = 0; iField < m_config.sFields.size(); iField++) {
0146       indices[iField] = m_decoder -> index(m_config.sFields[iField]);
0147     }
0148   } catch (const std::runtime_error &err) {
0149     throw std::runtime_error("PANIC: something went wrong grabbing the decoder!");
0150   }
0151   return;
0152
0153 }  // end 'InitializeDecoder()'
0154
0155
0156
0157 void FlatHCalDebugTreeMakerProcessor::InitializeTree() {
0158
0159   // switch to output directory
0160   m_pluginDir -> cd();
0161
0162   // instantiate output tree
0163   m_outTree = new TTree("FlatDebugTree", "Flat Tree with Particles, Cells, and Clusters");
0164   m_outTree -> SetDirectory(m_pluginDir);
0165
0166   // set branches
0167   m_outTree -> Branch("EvtIndex",     &m_evtIndex, "EvtIndex/I");
0168   m_outTree -> Branch("EvtNGenPar",   &m_nGen,     "EvtNGenPar/I");
0169   m_outTree -> Branch("EvtNMCPar",    &m_nMC,      "EvtNMCPar/I");
0170   m_outTree -> Branch("EvtNCell",     &m_nCells,   "EvtNCell/I");
0171   m_outTree -> Branch("EvtNClust",    &m_nClust,   "EvtNClust/I");
0172   m_outTree -> Branch("GenType",      &m_genType);
0173   m_outTree -> Branch("GenPDG",       &m_genPDG);
0174   m_outTree -> Branch("GenE",         &m_genEne);
0175   m_outTree -> Branch("GenPhi",       &m_genPhi);
0176   m_outTree -> Branch("GenEta",       &m_genEta);
0177   m_outTree -> Branch("GenMass",      &m_genMass);
0178   m_outTree -> Branch("MCGenStat",    &m_mcGenStat);
0179   m_outTree -> Branch("MCSimStat",    &m_mcSimStat);
0180   m_outTree -> Branch("MCPDG",        &m_mcPDG);
0181   m_outTree -> Branch("MCE",          &m_mcEne);
0182   m_outTree -> Branch("MCPhi",        &m_mcPhi);
0183   m_outTree -> Branch("MCEta",        &m_mcEta);
0184   m_outTree -> Branch("MCMass",       &m_mcMass);
0185   m_outTree -> Branch("MCStartVx",    &m_mcStartVX);
0186   m_outTree -> Branch("MCStartVy",    &m_mcStartVY);
0187   m_outTree -> Branch("MCStartVz",    &m_mcStartVZ);
0188   m_outTree -> Branch("MCStopVx",     &m_mcStopVX);
0189   m_outTree -> Branch("MCStopVy",     &m_mcStopVY);
0190   m_outTree -> Branch("MCStopVz",     &m_mcStopVZ);
0191   m_outTree -> Branch("MCTime",       &m_mcTime);
0192   m_outTree -> Branch("CellEne",      &m_cellEne);
0193   m_outTree -> Branch("CellPosX",     &m_cellRX);
0194   m_outTree -> Branch("CellPosY",     &m_cellRY);
0195   m_outTree -> Branch("CellPosZ",     &m_cellRZ);
0196   m_outTree -> Branch("CellTime",     &m_cellTime);
0197   m_outTree -> Branch("CellID",       &m_cellID);
0198   m_outTree -> Branch("CellIndexA",   &m_cellIndexA);
0199   m_outTree -> Branch("CellIndexB",   &m_cellIndexB);
0200   m_outTree -> Branch("CellIndexC",   &m_cellIndexC);
0201   m_outTree -> Branch("CellIndexE",   &m_cellIndexE);
0202   m_outTree -> Branch("CellIndexF",   &m_cellIndexF);
0203   m_outTree -> Branch("CellIndexG",   &m_cellIndexG);
0204   m_outTree -> Branch("ClustIndex",   &m_clustIndex);
0205   m_outTree -> Branch("ClustNCells",  &m_clustNCells);
0206   m_outTree -> Branch("ClustEne",     &m_clustEne);
0207   m_outTree -> Branch("ClustEta",     &m_clustEta);
0208   m_outTree -> Branch("ClustPhi",     &m_clustPhi);
0209   m_outTree -> Branch("ClustPosX",    &m_clustRX);
0210   m_outTree -> Branch("ClustPosY",    &m_clustRY);
0211   m_outTree -> Branch("ClustPosZ",    &m_clustRZ);
0212   m_outTree -> Branch("ClustTime",    &m_clustTime);
0213   return;
0214
0215 } // end 'InitializeTrees()'
0216
0217
0218
0219 void FlatHCalDebugTreeMakerProcessor::ResetVariables() {
0220
0221   // reset event variables
0222   m_nGen   = 0;
0223   m_nMC    = 0;
0224   m_nCells = 0;
0225   m_nClust = 0;
0226
0227   // reset gen particle variables
0228   m_genType.clear();
0229   m_genPDG.clear();
0230   m_genEne.clear();
0231   m_genPhi.clear();
0232   m_genEta.clear();
0233   m_genMass.clear();
0234
0235   // reset event variables
0236   m_mcGenStat.clear();
0237   m_mcSimStat.clear();
0238   m_mcPDG.clear();
0239   m_mcEne.clear();
0240   m_mcPhi.clear();
0241   m_mcEta.clear();
0242   m_mcMass.clear();
0243   m_mcStartVX.clear();
0244   m_mcStartVY.clear();
0245   m_mcStartVZ.clear();
0246   m_mcStopVX.clear();
0247   m_mcStopVY.clear();
0248   m_mcStopVZ.clear();
0249   m_mcTime.clear();
0250
0251   // reset cell members
0252   m_cellEne.clear();
0253   m_cellRX.clear();
0254   m_cellRY.clear();
0255   m_cellRZ.clear();
0256   m_cellTime.clear();
0257   m_cellIndexA.clear();
0258   m_cellIndexB.clear();
0259   m_cellIndexC.clear();
0260   m_cellIndexE.clear();
0261   m_cellIndexF.clear();
0262   m_cellIndexG.clear();
0263
0264   // reset cluster variables
0265   m_clustNCells.clear();
0266   m_clustEne.clear();
0267   m_clustEta.clear();
0268   m_clustPhi.clear();
0269   m_clustRX.clear();
0270   m_clustRY.clear();
0271   m_clustRZ.clear();
0272   m_clustTime.clear();
0273   return;
0274
0275 }  // end 'ResetVariables()'
0276
0277
0278
0279 void FlatHCalDebugTreeMakerProcessor::FillMCVariables(const edm4hep::MCParticle* mc) {
0280
0281   // grab particle kinematics
0282   ROOT::Math::PxPyPzM4D<float> pMC(
0283     mc -> getMomentum().x,
0284     mc -> getMomentum().y,
0285     mc -> getMomentum().z,
0286     mc -> getMass()
0287   );
0288
0289   // set output variables
0290   m_mcGenStat.push_back( mc -> getGeneratorStatus() );
0291   m_mcSimStat.push_back( mc -> getSimulatorStatus() );
0292   m_mcPDG.push_back( mc -> getPDG() );
0293   m_mcEne.push_back( pMC.E() );
0294   m_mcPhi.push_back( pMC.Phi() );
0295   m_mcEta.push_back( pMC.Eta() );
0296   m_mcMass.push_back( mc -> getMass() );
0297   m_mcStartVX.push_back( mc -> getVertex().x );
0298   m_mcStartVY.push_back( mc -> getVertex().y );
0299   m_mcStartVZ.push_back( mc -> getVertex().z );
0300   m_mcStopVX.push_back( mc -> getEndpoint().x );
0301   m_mcStopVY.push_back( mc -> getEndpoint().y );
0302   m_mcStopVZ.push_back( mc -> getEndpoint().z );
0303   m_mcTime.push_back( mc -> getTime() );
0304   return;
0305
0306 }  // end 'FillMCVariables(edm4hep::MCParticle*)'
0307
0308
0309
0310 void FlatHCalDebugTreeMakerProcessor::FillGenVariables(const edm4eic::ReconstructedParticle* gen) {
0311
0312   // get particle kinematics
0313   ROOT::Math::PxPyPzE4D<float> pGen(
0314     gen -> getMomentum().x,
0315     gen -> getMomentum().y,
0316     gen -> getMomentum().z,
0317     gen -> getEnergy()
0318   );
0319
0320   // set output variables
0321   m_genType.push_back( gen -> getType() );
0322   m_genPDG.push_back( gen -> getPDG() );
0323   m_genEne.push_back( pGen.E() );
0324   m_genEta.push_back( pGen.Eta() );
0325   m_genPhi.push_back( pGen.Phi() );
0326   m_genMass.push_back( gen -> getMass() );
0327   return;
0328
0329 }  // end 'FillGenVariables(const edm4eic::ReconstructedParticle*)'
0330
0331
0332
0333 void FlatHCalDebugTreeMakerProcessor::FillCellVariables(const edm4eic::CalorimeterHit* cell) {
0334
0335   // create vector to hold indices
0336   std::vector<short> indices(m_const.nFieldMax, -1);
0337
0338   // get hit indices
0339   const int64_t cellID = cell -> getCellID();
0340   GetCellIndices(cellID, indices);
0341
0342   // calculate time
0343   const double time = GetTime(cell -> getTime());
0344
0345   // set output variables
0346   m_cellEne.push_back( cell -> getEnergy() );
0347   m_cellRX.push_back( cell -> getPosition().x );
0348   m_cellRY.push_back( cell -> getPosition().y );
0349   m_cellRZ.push_back( cell -> getPosition().z );
0350   m_cellTime.push_back( time );
0351   m_cellID.push_back( cellID );
0352   m_cellIndexA.push_back( indices[0] );
0353   m_cellIndexB.push_back( indices[1] );
0354   m_cellIndexC.push_back( indices[2] );
0355   m_cellIndexE.push_back( indices[2] );
0356   m_cellIndexF.push_back( indices[3] );
0357   m_cellIndexG.push_back( indices[5] );
0358   return;
0359
0360 }  // end 'FillCellVariables(edm4eic::CalorimeterHit*)'
0361
0362
0363
0364 void FlatHCalDebugTreeMakerProcessor::FillClusterVariables(const edm4eic::Cluster* clust, const int64_t iClust) {
0365
0366   // calculate eta
0367   const double theta = clust -> getIntrinsicTheta();
0368   const double eta   = GetEta(theta);
0369
0370   m_clustIndex.push_back( iClust );
0371   m_clustNCells.push_back( clust -> getNhits() );
0372   m_clustEne.push_back( clust -> getEnergy() );
0373   m_clustEta.push_back( eta );
0374   m_clustPhi.push_back( clust -> getIntrinsicPhi() );
0375   m_clustRX.push_back( clust -> getPosition().x );
0376   m_clustRY.push_back( clust -> getPosition().y );
0377   m_clustRZ.push_back( clust -> getPosition().z );
0378   m_clustTime.push_back( clust -> getTime() );
0379   return;
0380
0381 }  // end 'FillClusterVariables(edm4eic::Cluster*, int64_t)'
0382
0383
0384
0385 void FlatHCalDebugTreeMakerProcessor::GetCellIndices(const int64_t cellID, std::vector<short>& indices) {
0386
0387   for (size_t iField = 0; iField < m_config.sFields.size(); iField++) {
0388     if (iField < m_const.nFieldMax) {
0389       indices.at(iField) = m_decoder -> get(cellID, m_config.sFields[iField].data());
0390     }
0391   }
0392   return;
0393
0394 }  // end 'GetCellIndices(int64_t, vector<short>&)'
0395
0396
0397
0398 double FlatHCalDebugTreeMakerProcessor::GetTime(const double tIn) {
0399
0400   // set max time
0401   const double maxTime = std::numeric_limits<double>::max();
0402
0403   // if above, return max
0404   double tOut = tIn;
0405   if (tOut > maxTime) {
0406     tOut = maxTime;
0407   }
0408   return tOut;
0409
0410 }  // end 'GetTime(double)'
0411
0412
0413
0414 double FlatHCalDebugTreeMakerProcessor::GetEta(const double theta) {
0415
0416   const double expo = std::tan(theta / 2.);
0417   const double eta  = -1. * std::log(expo);
0418   return eta;
0419
0420 }  // end 'GetEta(double)'
0421
0422 // end ------------------------------------------------------------------------