Calorimetery/MatchClustersToTrackProjections/MatchProjectionsAndClusters.cxx

0001 // ----------------------------------------------------------------------------
0002 // 'MatchProjectionsAndClusters.cxx'
0003 // Derek Anderson
0004 // 05.20.2024
0005 //
0006 // Quick ROOT macro to match calorimeter clusters
0007 // to track projections.
0008 //
0009 // Usage:
0010 //   root -b -q MatchProjectionsAndClusters.cxx++
0011 // ----------------------------------------------------------------------------
0012
0013 #define MATCHPROJECTIONSANDCLUSTERS_CXX
0014
0015 // c++ utilities
0016 #include <cmath>
0017 #include <string>
0018 #include <limits>
0019 #include <optional>
0020 #include <iostream>
0021 // root libraries
0022 #include <TSystem.h>
0023 // podio libraries
0024 #include <podio/Frame.h>
0025 #include <podio/CollectionBase.h>
0026 #include <podio/ROOTFrameReader.h>
0027 // edm4eic types
0028 #include <edm4eic/Cluster.h>
0029 #include <edm4eic/ClusterCollection.h>
0030 #include <edm4eic/TrackPoint.h>
0031 #include <edm4eic/TrackSegment.h>
0032 #include <edm4eic/TrackSegmentCollection.h>
0033 // edm4hep types
0034 #include <edm4hep/Vector3f.h>
0035 #include <edm4hep/utils/vector_utils.h>
0036
0037
0038
0039 // user options ---------------------------------------------------------------
0040
0041 struct Options {
0042   std::string in_file;      // input file
0043   std::string out_file;     // output file
0044   std::string clusters;     // name of collection of clusters to match to
0045   std::string projections;  // name of collection of track projections
0046   uint64_t    system;       // system id of calo to match to
0047 } DefaultOptions = {
0048   "normal/forQuickAssociationAnswers_normalMode_nEvt500.py8s18x275minq100ang0025div1.run0file0.podio.root",
0049   "test.root",
0050   "HcalEndcapNClusters",
0051   "CalorimeterTrackProjections",
0052   113
0053 };
0054
0055
0056
0057 // macro body -----------------------------------------------------------------
0058
0059 void MatchProjectionsAndClusters(const Options& opt = DefaultOptions) {
0060
0061   // announce start of macro
0062   std::cout << "\n  Beginning cluster-track projection matching macro!" << std::endl;
0063
0064   // open file w/ frame reader
0065   podio::ROOTFrameReader reader = podio::ROOTFrameReader();
0066   reader.openFile( opt.in_file );
0067   std::cout << "    Opened ROOT-based frame reader." << std::endl;
0068
0069   // get no. of frames and annoucne
0070   const uint64_t nFrames = reader.getEntries(podio::Category::Event);
0071   std::cout << "    Starting frame loop: " << reader.getEntries(podio::Category::Event) << " frames to process." << std::endl;
0072
0073   // iterate through frames (i.e. events in this case)
0074   uint64_t nClustTotal   = 0;
0075   uint64_t nClustMatched = 0;
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& segments = frame.get<edm4eic::TrackSegmentCollection>( opt.projections );
0087
0088     // loop over clusters
0089     for (size_t iClust = 0; edm4eic::Cluster cluster : clusters) {
0090
0091       // grab eta/phi of cluster
0092       const double etaClust  = edm4hep::utils::eta( cluster.getPosition() );
0093       const double phiClust  = edm4hep::utils::angleAzimuthal( cluster.getPosition() );
0094
0095       // match based on eta/phi dstiance
0096       double distMatch = std::numeric_limits<double>::max();
0097
0098       // loop over projections to find matching one
0099       std::optional<edm4eic::TrackPoint> match;
0100       for (edm4eic::TrackSegment segment : segments) {
0101         for (edm4eic::TrackPoint projection : segment.getPoints()) {
0102
0103           // ignore if not pointing to calo or at face of calo
0104           const bool isInSystem = (projection.system == opt.system);
0105           const bool isAtFace   = (projection.surface == 1);
0106           if (!isInSystem || !isAtFace) continue;
0107
0108           // grab eta/phi of projection
0109           const double etaProject = edm4hep::utils::eta( projection.position );
0110           const double phiProject = edm4hep::utils::angleAzimuthal( projection.position );
0111
0112           // get distance to projection
0113           const double distProject = std::hypot(
0114             etaClust - etaProject,
0115             phiClust - phiProject
0116           );
0117
0118           // if smallest distance found, update variables accordingly
0119           if (distProject < distMatch) {
0120             distMatch = distProject;
0121             match     = projection;
0122           }
0123
0124         }  // end point loop
0125       }  // end segment loop
0126
0127       // do analysis if match was found...
0128       if ( match.has_value() ) {
0129         edm4eic::TrackPoint matchedProject = match.value();
0130         std::cout << "        Found match! match = " << matchedProject << std::endl;
0131         ++nClustMatched;
0132       }
0133       ++nClustTotal;
0134
0135     }  // end cluster loop
0136   }  // end frame loop
0137   std::cout << "    Finished frame loop: " << nClustMatched << "/" << nClustTotal << " clusters matched." << std::endl;
0138
0139   // announce end and exit
0140   std::cout << "  End of macro!\n" << std::endl;
0141   return;
0142
0143 }
0144
0145 // end ------------------------------------------------------------------------