Calorimetery/ReadAssociationsAndContributions/ReadInAssociationsAndContributions.cxx

0001 // ----------------------------------------------------------------------------
0002 // 'ReadInAssociationsAndContributions.cxx'
0003 // Derek Anderson
0004 // 06.10.2024
0005 //
0006 // Quick ROOT macro to read in both the
0007 // associations of and contributions
0008 // to a reconstructed cluster.
0009 //
0010 // Usage:
0011 //   root -b -q ReadInAssociationsAndContributions.cxx++
0012 // ----------------------------------------------------------------------------
0013
0014 #define READINASSOCIATIONSANDCONTRIBUTIONS_CXX
0015
0016 // c++ utilities
0017 #include <cmath>
0018 #include <string>
0019 #include <limits>
0020 #include <optional>
0021 #include <iostream>
0022 // root libraries
0023 #include <TSystem.h>
0024 // podio libraries
0025 #include <podio/Frame.h>
0026 #include <podio/CollectionBase.h>
0027 #include <podio/ROOTFrameReader.h>
0028 // edm4eic types
0029 #include <edm4eic/Cluster.h>
0030 #include <edm4eic/ClusterCollection.h>
0031 #include <edm4eic/MCRecoClusterParticleAssociation.h>
0032 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
0033 #include <edm4eic/CalorimeterHit.h>
0034 // edm4hep types
0035 #include <edm4hep/SimCalorimeterHit.h>
0036 #include <edm4hep/SimCalorimeterHitCollection.h>
0037 #include <edm4hep/CaloHitContribution.h>
0038
0039
0040
0041 // user options ---------------------------------------------------------------
0042
0043 struct Options {
0044   std::string in_file;       // input file
0045   std::string out_file;      // output file
0046   std::string clusters;      // name of collection of clusters
0047   std::string associations;  // name of collection of cluster-particle associations
0048   std::string sim_hits;      // name of collection of mc calorimeter hits
0049 } DefaultOptions = {
0050   "../normal/forQuickAssociationAnswers_normalMode_nEvt500.py8s18x275minq100ang0025div1.run0file0.podio.root",
0051   "test.root",
0052   "LFHCALClusters",
0053   "LFHCALClusterAssociations",
0054   "LFHCALHits"
0055 };
0056
0057
0058
0059 // macro body -----------------------------------------------------------------
0060
0061 void ReadInAssociationsAndContributions(const Options& opt = DefaultOptions) {
0062
0063   // announce start of macro
0064   std::cout << "\n  Beginning cluster-track projection matching macro!" << std::endl;
0065
0066   // open file w/ frame reader
0067   podio::ROOTFrameReader reader = podio::ROOTFrameReader();
0068   reader.openFile( opt.in_file );
0069   std::cout << "    Opened ROOT-based frame reader." << std::endl;
0070
0071   // get no. of frames and annoucne
0072   const uint64_t nFrames = reader.getEntries(podio::Category::Event);
0073   std::cout << "    Starting frame loop: " << reader.getEntries(podio::Category::Event) << " frames to process." << std::endl;
0074
0075   // iterate through frames (i.e. events in this case)
0076   for (uint64_t iFrame = 0; iFrame < nFrames; ++iFrame) {
0077
0078     // announce progress
0079     std::cout << "      Processing frame " << iFrame + 1 << "/" << nFrames << "..." << std::endl;
0080
0081     // grab frame
0082     auto frame = podio::Frame( reader.readNextEntry(podio::Category::Event) );
0083
0084     // grab collections
0085     auto& clusters     = frame.get<edm4eic::ClusterCollection>( opt.clusters );
0086     auto& associations = frame.get<edm4eic::MCRecoClusterParticleAssociationCollection>( opt.associations );
0087     auto& sim_hits     = frame.get<edm4hep::SimCalorimeterHitCollection>( opt.sim_hits );
0088
0089     // look at associations ---------------------------------------------------
0090
0091     // loop over associations
0092     for (size_t iAssoc = 0; edm4eic::MCRecoClusterParticleAssociation assoc : associations) {
0093
0094       // grab associated cluster, particle
0095       auto cluster  = assoc.getRec();
0096       auto particle = assoc.getSim();
0097
0098       // print IDs of both objects
0099       std::cout << "        [Assoc. #" << iAssoc << "]: (cluster, particle) ID = (" << assoc.getRecID() << ", " << assoc.getSimID() << ")" << std::endl;
0100
0101       // now print some additional info from each
0102       std::cout << "          Cluster:  energy = " << cluster.getEnergy() << ", nHits = " << cluster.getNhits() << "\n"
0103                 << "          Particle: energy = " << particle.getEnergy() << ", pdg = " << particle.getPDG() << ", gen status = " << particle.getGeneratorStatus()
0104                 << std::endl;
0105       ++iAssoc;
0106
0107     }  // end association loop
0108
0109     // look at contributions --------------------------------------------------
0110
0111     // loop over clusters
0112     for (size_t iClust = 0; edm4eic::Cluster cluster : clusters) {
0113
0114       // loop through reconstructed cluster hits (cells)
0115       for (size_t iRec = 0; edm4eic::CalorimeterHit rec : cluster.getHits()) {
0116
0117         // grab cell ID
0118         const uint64_t recCellID = rec.getCellID();
0119
0120         // loop over sim hits
0121         for (size_t iSim = 0; edm4hep::SimCalorimeterHit sim : sim_hits) {
0122
0123           // match sim-to-reco hit based on cell ID
0124           const uint64_t simCellID  = sim.getCellID();
0125           const bool     isSameCell = (recCellID == simCellID);
0126           if (!isSameCell) continue;
0127
0128           // now loop over contributions
0129           for (size_t iContrib = 0; edm4hep::CaloHitContribution contrib : sim.getContributions()) {
0130
0131             // grab corresponding particle
0132             auto particle = contrib.getParticle();
0133
0134             // print IDs of all 5 objects in hierarchy
0135             std::cout << "        [Contrib #" << iContrib << "] (cluster, reco hit, sim hit, contrib, particle) ID = ("
0136                       << cluster.getObjectID().index  << ", "
0137                       << rec.getObjectID().index      << ", "
0138                       << sim.getObjectID().index      << ", "
0139                       << contrib.getObjectID().index  << ", "
0140                       << particle.getObjectID().index << ")"
0141                       << std::endl;
0142
0143             // and now print some info from each layer
0144             std::cout << "          Cluster:  energy = " << cluster.getEnergy() << ", nHits = << " << cluster.getNhits() << "\n"
0145                       << "          Rec Hit:  energy = " << rec.getEnergy() << ", time = " << rec.getTime() << "\n"
0146                       << "          Sim Hit:  energy = " << sim.getEnergy() << "\n"
0147                       << "          Contrib:  energy = " << contrib.getEnergy() << ", time = " << contrib.getTime() << "\n"
0148                       << "          Particle: energy = " << particle.getEnergy() << ", pdg = " << particle.getPDG() << ", gen status = " << particle.getGeneratorStatus()
0149                       << std::endl;
0150             ++iContrib;
0151
0152           }  // end contrib loop
0153           ++iSim;
0154
0155         }  // end sim hit loop
0156         ++iRec;
0157
0158       }  // end reco hit (cell) loop
0159       ++iClust;
0160
0161     }  // end cluster loop
0162   }  // end frame loop
0163   std::cout << "    Finished frame loop"  << std::endl;
0164
0165   // announce end and exit
0166   std::cout << "  End of macro!\n" << std::endl;
0167   return;
0168
0169 }
0170
0171 // end ------------------------------------------------------------------------