CaloDebugTools/HCalDebugTreeMaker/RelationalHCalDebugTreeMakerProcessor.cc

0001 // ----------------------------------------------------------------------------
0002 // 'RelationalDebugHCalTreeMaker.cc'
0003 // Derek Anderson
0004 // 03.09.2024
0005 //
0006 // A JANA plugin to construct a tree of HCal clusters
0007 // and related cells and particles.
0008 // ----------------------------------------------------------------------------
0009
0010 #define RELATIONALHCALDEBUGTREEMAKERPROCESSOR_CC
0011
0012 // plugin definition
0013 #include "RelationalHCalDebugTreeMakerProcessor.h"
0014
0015 // The following just makes this a JANA plugin
0016 extern "C" {
0017     void InitPlugin(JApplication *app) {
0018         InitJANAPlugin(app);
0019         app->Add(new RelationalHCalDebugTreeMakerProcessor);
0020     }
0021 }
0022
0023
0024
0025 // public methods -------------------------------------------------------------
0026
0027 void RelationalHCalDebugTreeMakerProcessor::InitWithGlobalRootLock(){
0028
0029   //  grab output file
0030   auto rootfile_svc = GetApplication() -> GetService<RootFile_service>();
0031   auto rootfile     = rootfile_svc     -> GetHistFile();
0032
0033   // create directory
0034   m_pluginDir = rootfile -> mkdir(m_config.sPlugin.data());
0035   m_pluginDir -> cd();
0036   // reset tree variables
0037   ResetVariables();
0038
0039   // initialize event index
0040   m_evtIndex = 0;
0041
0042   // initialize output trees
0043   InitializeDecoder();
0044   InitializeTree();
0045   return;
0046
0047 }  // end 'InitWithGlobalLock()'
0048
0049
0050
0051 void RelationalHCalDebugTreeMakerProcessor::ProcessSequential(const std::shared_ptr<const JEvent>& event) {
0052
0053   // reset tree variables
0054   ResetVariables();
0055
0056   // loop over clusters
0057   int64_t nClust = 0;
0058   for (auto cluster : m_clusters()) {
0059
0060     // loop over associated mc particles
0061     int64_t nAssocToClust = 0;
0062     for (auto assoc : m_assocs()) {
0063
0064       // consider only current cluster
0065       const bool isSameClust = (assoc -> getRecID() == nClust);
0066       if (!isSameClust) continue;
0067
0068       // grab associated particle
0069       auto mcPar = assoc -> getSim();
0070
0071       FillAssocVariables( mcPar, nClust );
0072       ++nAssocToClust;
0073
0074     }  // end association loop
0075
0076     // grab consituent hits
0077     const auto clustHits = cluster -> getHits();
0078
0079     // associate each cell, contributing particle with corresponding cluster
0080     int64_t nCellsInClust    = 0;
0081     int64_t nContribToClust = 0;
0082     for (auto clustHit : clustHits) {
0083
0084       // get cell ID
0085       const uint64_t cellID = clustHit.getCellID();
0086
0087       // get contributing particles
0088       for (auto simHit : m_simHits()) {
0089
0090         // identify sim hit based on cell ID
0091         const bool isSameCell = (cellID == simHit -> getCellID());
0092         if (!isSameCell) continue;
0093
0094         // loop over hit contributions
0095         const auto contribs = simHit -> getContributions();
0096         for (auto contrib : contribs) {
0097
0098           // grab contributing particle
0099           auto mcPar = contrib.getParticle();
0100
0101           FillContribVariables( mcPar, nClust, cellID );
0102           ++nContribToClust;
0103
0104         }  // end contribution loop
0105       }  // end sim hit loop
0106
0107       FillCellVariables( clustHit, nClust );
0108       ++nCellsInClust;
0109
0110     }  // end cell loop
0111
0112
0113     FillClusterVariables( cluster, nClust, nCellsInClust, nAssocToClust, nContribToClust );
0114     ++nClust;
0115
0116   }  // end cluster loop
0117
0118   // set event variables
0119   m_nClust = nClust;
0120
0121   // fill output tree
0122   m_outTree -> Fill();
0123
0124   // increment event index and exit routine
0125   ++m_evtIndex;
0126   return;
0127
0128 }  // end 'ProcessSequential(std::shared_ptr<JEvent>&)'
0129
0130
0131
0132 void RelationalHCalDebugTreeMakerProcessor::FinishWithGlobalRootLock() {
0133
0134   // clean up variables
0135   ResetVariables();
0136   return;
0137
0138 }  // end 'FinishWithGlobalRootLock()'
0139
0140
0141
0142 // internal methods -----------------------------------------------------------
0143
0144 void RelationalHCalDebugTreeMakerProcessor::InitializeDecoder() {
0145
0146   // grab detector
0147   auto detector = GetApplication() -> GetService<DD4hep_service>() -> detector();
0148
0149   // make sure readout is available
0150   dd4hep::IDDescriptor descriptor;
0151   try {
0152     descriptor  = detector -> readout(m_config.sSimHits.data()).idSpec();
0153   } catch (const std::runtime_error &err) {
0154     throw std::runtime_error("PANIC: readout class is not in output!");
0155   }
0156
0157   // grab decoder and test
0158   std::vector<short> indices(m_config.sFields.size());
0159   try {
0160     m_decoder = descriptor.decoder();
0161     for (size_t iField = 0; iField < m_config.sFields.size(); iField++) {
0162       indices[iField] = m_decoder -> index(m_config.sFields[iField]);
0163     }
0164   } catch (const std::runtime_error &err) {
0165     throw std::runtime_error("PANIC: something went wrong grabbing the decoder!");
0166   }
0167   return;
0168
0169 }  // end 'InitializeDecoder()'
0170
0171
0172
0173 void RelationalHCalDebugTreeMakerProcessor::InitializeTree() {
0174
0175   // switch to output directory
0176   m_pluginDir -> cd();
0177
0178   // instantiate output tree
0179   m_outTree = new TTree("RelationalDebugTree", "Tree with Cluster-Cell/Particle Relations ");
0180   m_outTree -> SetDirectory(m_pluginDir);
0181
0182   // set branches
0183   m_outTree -> Branch("EvtIndex",          &m_evtIndex, "EvtIndex/I");
0184   m_outTree -> Branch("EvtNClust",         &m_nClust,   "EvtNClust/I");
0185   m_outTree -> Branch("ClustIndex",        &m_clustIndex);
0186   m_outTree -> Branch("ClustNCell",        &m_clustNCells);
0187   m_outTree -> Branch("ClustNAssocPar",    &m_clustNAssoc);
0188   m_outTree -> Branch("ClustNContribPar",  &m_clustNContrib);
0189   m_outTree -> Branch("ClustEne",          &m_clustEne);
0190   m_outTree -> Branch("ClustEta",          &m_clustEta);
0191   m_outTree -> Branch("ClustPhi",          &m_clustPhi);
0192   m_outTree -> Branch("ClustPosX",         &m_clustRX);
0193   m_outTree -> Branch("ClustPosY",         &m_clustRY);
0194   m_outTree -> Branch("ClustPosZ",         &m_clustRZ);
0195   m_outTree -> Branch("ClustTime",         &m_clustTime);
0196   m_outTree -> Branch("CelllID",           &m_cellID);
0197   m_outTree -> Branch("CellClustIndex",    &m_cellClustIndex);
0198   m_outTree -> Branch("CellEne",           &m_cellEne);
0199   m_outTree -> Branch("CellPosX",          &m_cellRX);
0200   m_outTree -> Branch("CellPosY",          &m_cellRY);
0201   m_outTree -> Branch("CellPosZ",          &m_cellRZ);
0202   m_outTree -> Branch("CellTime",          &m_cellTime);
0203   m_outTree -> Branch("CellIndexA",        &m_cellIndexA);
0204   m_outTree -> Branch("CellIndexB",        &m_cellIndexB);
0205   m_outTree -> Branch("CellIndexC",        &m_cellIndexC);
0206   m_outTree -> Branch("CellIndexE",        &m_cellIndexE);
0207   m_outTree -> Branch("CellIndexF",        &m_cellIndexF);
0208   m_outTree -> Branch("CellIndexG",        &m_cellIndexG);
0209   m_outTree -> Branch("AssocClustIndex",   &m_assocClustIndex);
0210   m_outTree -> Branch("AssocGenStat",      &m_assocGenStat);
0211   m_outTree -> Branch("AssocSimStat",      &m_assocSimStat);
0212   m_outTree -> Branch("AssocPDG",          &m_assocPDG);
0213   m_outTree -> Branch("AssocEne",          &m_assocEne);
0214   m_outTree -> Branch("AssocPhi",          &m_assocPhi);
0215   m_outTree -> Branch("AssocEta",          &m_assocEta);
0216   m_outTree -> Branch("AssocMass",         &m_assocMass);
0217   m_outTree -> Branch("AssocStartVX",      &m_assocStartVX);
0218   m_outTree -> Branch("AssocStartVY",      &m_assocStartVY);
0219   m_outTree -> Branch("AssocStartVZ",      &m_assocStartVZ);
0220   m_outTree -> Branch("AssocStopVX",       &m_assocStopVX);
0221   m_outTree -> Branch("AssocStopVY",       &m_assocStopVY);
0222   m_outTree -> Branch("AssocStopVZ",       &m_assocStopVZ);
0223   m_outTree -> Branch("AssocTime",         &m_assocTime);
0224   m_outTree -> Branch("ContribClustIndex", &m_contribClustIndex);
0225   m_outTree -> Branch("ContribCellID",     &m_contribCellID);
0226   m_outTree -> Branch("ContribGenStat",    &m_contribGenStat);
0227   m_outTree -> Branch("ContribSimStat",    &m_contribSimStat);
0228   m_outTree -> Branch("ContribPDG",        &m_contribPDG);
0229   m_outTree -> Branch("ContribE",          &m_contribEne);
0230   m_outTree -> Branch("ContribPhi",        &m_contribPhi);
0231   m_outTree -> Branch("ContribEta",        &m_contribEta);
0232   m_outTree -> Branch("ContribMass",       &m_contribMass);
0233   m_outTree -> Branch("ContribStartVX",    &m_contribStartVX);
0234   m_outTree -> Branch("ContribStartVY",    &m_contribStartVY);
0235   m_outTree -> Branch("ContribStartVZ",    &m_contribStartVZ);
0236   m_outTree -> Branch("ContribStopVX",     &m_contribStopVX);
0237   m_outTree -> Branch("ContribStopVY",     &m_contribStopVY);
0238   m_outTree -> Branch("ContribStopVZ",     &m_contribStopVZ);
0239   m_outTree -> Branch("ContribTime",       &m_contribTime);
0240   return;
0241
0242 } // end 'InitializeTrees()'
0243
0244
0245
0246 void RelationalHCalDebugTreeMakerProcessor::ResetVariables() {
0247
0248   // reset event variables
0249   m_nClust = 0;
0250
0251   // reset cluster variables
0252   m_clustNCells.clear();
0253   m_clustNContrib.clear();
0254   m_clustNAssoc.clear();
0255   m_clustEne.clear();
0256   m_clustEta.clear();
0257   m_clustPhi.clear();
0258   m_clustRX.clear();
0259   m_clustRY.clear();
0260   m_clustRZ.clear();
0261   m_clustTime.clear();
0262
0263   // reset cells in cluster variables
0264   m_cellID.clear();
0265   m_cellClustIndex.clear();
0266   m_cellEne.clear();
0267   m_cellRX.clear();
0268   m_cellRY.clear();
0269   m_cellRZ.clear();
0270   m_cellTime.clear();
0271   m_cellIndexA.clear();
0272   m_cellIndexB.clear();
0273   m_cellIndexC.clear();
0274   m_cellIndexE.clear();
0275   m_cellIndexF.clear();
0276   m_cellIndexG.clear();
0277
0278   // reset mc particle associated to cluster variables
0279   m_assocClustIndex.clear();
0280   m_assocGenStat.clear();
0281   m_assocSimStat.clear();
0282   m_assocPDG.clear();
0283   m_assocEne.clear();
0284   m_assocPhi.clear();
0285   m_assocEta.clear();
0286   m_assocMass.clear();
0287   m_assocStartVX.clear();
0288   m_assocStartVY.clear();
0289   m_assocStartVZ.clear();
0290   m_assocStopVX.clear();
0291   m_assocStopVY.clear();
0292   m_assocStopVZ.clear();
0293   m_assocTime.clear();
0294
0295   // reset mc particle contributing to cluster variables
0296   m_contribClustIndex.clear();
0297   m_contribCellID.clear();
0298   m_contribGenStat.clear();
0299   m_contribSimStat.clear();
0300   m_contribPDG.clear();
0301   m_contribEne.clear();
0302   m_contribPhi.clear();
0303   m_contribEta.clear();
0304   m_contribMass.clear();
0305   m_contribStartVX.clear();
0306   m_contribStartVY.clear();
0307   m_contribStartVZ.clear();
0308   m_contribStopVX.clear();
0309   m_contribStopVY.clear();
0310   m_contribStopVZ.clear();
0311   m_contribTime.clear();
0312   return;
0313
0314 }  // end 'ResetVariables()'
0315
0316
0317
0318 void RelationalHCalDebugTreeMakerProcessor::FillClusterVariables(
0319   const edm4eic::Cluster* clust,
0320   const int64_t iClust,
0321   const int64_t nCells,
0322   const int64_t nAssoc,
0323   const int64_t nContrib
0324 ) {
0325
0326   // calculate eta
0327   const double theta = clust -> getIntrinsicTheta();
0328   const double eta   = GetEta(theta);
0329
0330   m_clustIndex.push_back( iClust );
0331   m_clustNCells.push_back( nCells );
0332   m_clustNAssoc.push_back( nAssoc );
0333   m_clustNContrib.push_back( nContrib );
0334   m_clustEne.push_back( clust -> getEnergy() );
0335   m_clustEta.push_back( eta );
0336   m_clustPhi.push_back( clust -> getIntrinsicPhi() );
0337   m_clustRX.push_back( clust -> getPosition().x );
0338   m_clustRY.push_back( clust -> getPosition().y );
0339   m_clustRZ.push_back( clust -> getPosition().z );
0340   m_clustTime.push_back( clust -> getTime() );
0341   return;
0342
0343 }  // end 'FillClusterVariables(edm4eic::Cluster*, int64_t)'
0344
0345
0346
0347 void RelationalHCalDebugTreeMakerProcessor::FillCellVariables(const edm4eic::CalorimeterHit cell, const int64_t iClust) {
0348
0349   // create vector to hold indices
0350   std::vector<short> indices(m_const.nFieldMax, -1);
0351
0352   // get hit indices
0353   const int64_t cellID = cell.getCellID();
0354   GetCellIndices(cellID, indices);
0355
0356   // calculate time
0357   const double time = GetTime(cell.getTime());
0358
0359   // set output variables
0360   m_cellID.push_back( cellID );
0361   m_cellClustIndex.push_back( iClust );
0362   m_cellEne.push_back( cell.getEnergy() );
0363   m_cellRX.push_back( cell.getPosition().x );
0364   m_cellRY.push_back( cell.getPosition().y );
0365   m_cellRZ.push_back( cell.getPosition().z );
0366   m_cellTime.push_back( time );
0367   m_cellIndexA.push_back( indices[0] );
0368   m_cellIndexB.push_back( indices[1] );
0369   m_cellIndexC.push_back( indices[2] );
0370   m_cellIndexE.push_back( indices[3] );
0371   m_cellIndexF.push_back( indices[4] );
0372   m_cellIndexG.push_back( indices[5] );
0373   return;
0374
0375 }  // end 'FillCellVariables(edm4eic::CalorimeterHit, int64_t)'
0376
0377
0378
0379 void RelationalHCalDebugTreeMakerProcessor::FillAssocVariables(const edm4hep::MCParticle assoc, const int64_t iClust) {
0380
0381   // grab particle kinematics
0382   ROOT::Math::PxPyPzM4D<float> pMC(
0383     assoc.getMomentum().x,
0384     assoc.getMomentum().y,
0385     assoc.getMomentum().z,
0386     assoc.getMass()
0387   );
0388
0389   // set output variables
0390   m_assocClustIndex.push_back( iClust );
0391   m_assocGenStat.push_back( assoc.getGeneratorStatus() );
0392   m_assocSimStat.push_back( assoc.getSimulatorStatus() );
0393   m_assocPDG.push_back( assoc.getPDG() );
0394   m_assocEne.push_back( pMC.E() );
0395   m_assocPhi.push_back( pMC.Phi() );
0396   m_assocEta.push_back( pMC.Eta() );
0397   m_assocMass.push_back( assoc.getMass() );
0398   m_assocStartVX.push_back( assoc.getVertex().x );
0399   m_assocStartVY.push_back( assoc.getVertex().y );
0400   m_assocStartVZ.push_back( assoc.getVertex().z );
0401   m_assocStopVX.push_back( assoc.getEndpoint().x );
0402   m_assocStopVY.push_back( assoc.getEndpoint().y );
0403   m_assocStopVZ.push_back( assoc.getEndpoint().z );
0404   m_assocTime.push_back( assoc.getTime() );
0405   return;
0406
0407 }  // end 'FillAssocVariables(edm4hep::MCParticle, int64_t)'
0408
0409
0410
0411 void RelationalHCalDebugTreeMakerProcessor::FillContribVariables(
0412   const edm4hep::MCParticle mc,
0413   const int64_t iClust,
0414   const int64_t cellID
0415 ) {
0416
0417   // grab particle kinematics
0418   ROOT::Math::PxPyPzM4D<float> pMC(
0419     mc.getMomentum().x,
0420     mc.getMomentum().y,
0421     mc.getMomentum().z,
0422     mc.getMass()
0423   );
0424
0425   // set output variables
0426   m_contribClustIndex.push_back( iClust );
0427   m_contribCellID.push_back( cellID );
0428   m_contribGenStat.push_back( mc.getGeneratorStatus() );
0429   m_contribSimStat.push_back( mc.getSimulatorStatus() );
0430   m_contribPDG.push_back( mc.getPDG() );
0431   m_contribEne.push_back( pMC.E() );
0432   m_contribPhi.push_back( pMC.Phi() );
0433   m_contribEta.push_back( pMC.Eta() );
0434   m_contribMass.push_back( mc.getMass() );
0435   m_contribStartVX.push_back( mc.getVertex().x );
0436   m_contribStartVY.push_back( mc.getVertex().y );
0437   m_contribStartVZ.push_back( mc.getVertex().z );
0438   m_contribStopVX.push_back( mc.getEndpoint().x );
0439   m_contribStopVY.push_back( mc.getEndpoint().y );
0440   m_contribStopVZ.push_back( mc.getEndpoint().z );
0441   m_contribTime.push_back( mc.getTime() );
0442   return;
0443
0444 }  // end 'FillContribVariables(edm4hep::MCParticle, int64_t, int64_t)'
0445
0446
0447
0448 void RelationalHCalDebugTreeMakerProcessor::GetCellIndices(const int64_t cellID, std::vector<short>& indices) {
0449
0450   for (size_t iField = 0; iField < m_config.sFields.size(); iField++) {
0451     if (iField < m_const.nFieldMax) {
0452       indices.at(iField) = m_decoder -> get(cellID, m_config.sFields[iField].data());
0453     }
0454   }
0455   return;
0456
0457 }  // end 'GetCellIndices(int64_t, vector<short>)'
0458
0459
0460
0461 double RelationalHCalDebugTreeMakerProcessor::GetTime(const double tIn) {
0462
0463   // set max time
0464   const double maxTime = std::numeric_limits<double>::max();
0465
0466   // if above, return max
0467   double tOut = tIn;
0468   if (tOut > maxTime) {
0469     tOut = maxTime;
0470   }
0471   return tOut;
0472
0473 }  // end 'GetTime(double)'
0474
0475
0476
0477 double RelationalHCalDebugTreeMakerProcessor::GetEta(const double theta) {
0478
0479   const double expo = std::tan(theta / 2.);
0480   const double eta  = -1. * std::log(expo);
0481   return eta;
0482
0483 }  // end 'GetEta(double)'
0484
0485 // end ------------------------------------------------------------------------
0486