CaloDebugTools/CheckClusterSplitting/CheckClusterSplittingProcessor.h

0001 // ----------------------------------------------------------------------------
0002 // 'CheckClusterSplittingProcessor.h'
0003 // Derek Anderson
0004 // 04.11.2024
0005 //
0006 // EICrecon plugin to collect some metrics on cluster splitting
0007 // across several calorimeters.
0008 // ----------------------------------------------------------------------------
0009
0010 #ifndef CHECKCLUSTERSPLITTINGPROCESSOR_H
0011 #define CHECKCLUSTERSPLITTINGPROCESSOR_H
0012
0013 // c++ utilities
0014 #include <map>
0015 #include <cmath>
0016 #include <limits>
0017 #include <string>
0018 #include <vector>
0019 #include <cstdlib>
0020 #include <utility>
0021 #include <algorithm>
0022 // root libraries
0023 #include <TH1.h>
0024 #include <TH2.h>
0025 #include <TFile.h>
0026 #include <TDirectory.h>
0027 // edm4hep types
0028 #include <edm4hep/Vector3f.h>
0029 #include <edm4hep/utils/vector_utils.h>
0030 #include <edm4hep/MCParticle.h>
0031 #include <edm4hep/MCParticleCollection.h>
0032 // edm4eic types
0033 #include <edm4eic/CalorimeterHit.h>
0034 #include <edm4eic/Cluster.h>
0035 #include <edm4eic/ClusterCollection.h>
0036 #include <edm4eic/TrackPoint.h>
0037 #include <edm4eic/TrackSegment.h>
0038 #include <edm4eic/TrackSegmentCollection.h>
0039 #include <edm4eic/ReconstructedParticle.h>
0040 #include <edm4eic/ReconstructedParticleCollection.h>
0041 // dd4hep utilities
0042 #include "DD4hep/Detector.h"
0043 #include "DD4hep/DetElement.h"
0044 // jana libraries
0045 #include <JANA/JEventProcessorSequentialRoot.h>
0046 // eicrecon services
0047 #include "services/rootfile/RootFile_service.h"
0048 #include "services/geometry/dd4hep/DD4hep_service.h"
0049
0050 // make common namespaces implicit
0051 using namespace std;
0052
0053
0054
0055 // convenience types ----------------------------------------------------------
0056
0057 // binning types
0058 typedef tuple<int64_t, float, float> Binning;
0059 typedef pair<string, Binning> AxisDef;
0060
0061 // histogram types
0062 typedef tuple<string, AxisDef> HistDef1D;
0063 typedef tuple<string, AxisDef, AxisDef> HistDef2D;
0064
0065 // definition vectors
0066 typedef vector<AxisDef> VecAxisDef;
0067 typedef vector<HistDef1D> VecHistDef1D;
0068 typedef vector<HistDef2D> VecHistDef2D;
0069
0070 // other types
0071 typedef pair<float, float> FPair;
0072 typedef vector<edm4eic::TrackPoint> TrkVec;
0073 typedef vector<vector<vector<TH1D*>>> TH1Vec;
0074 typedef vector<vector<vector<TH2D*>>> TH2Vec;
0075
0076
0077
0078 // histogram utilities --------------------------------------------------------
0079
0080 namespace Hist {
0081
0082   // indices/accessors --------------------------------------------------------
0083
0084   // calorimeters
0085   enum Calo {NHCal, NECal, BHCal, BECal, PHCal, PECal};
0086
0087   // histogram axes
0088   enum EVar {EvtNTr, EvtNCl, EvtPE, EvtLE, EvtTE, EvtLT, EvtLP, EvtTP, EvtNA90, EvtNB90, EvtNR};
0089   enum TVar {TrkRX, TrkRY, TrkRZ, TrkH, TrkF, TrkP, TrkDR, TrkDRS, TrkEP, TrkDep};
0090   enum CVar {CalNH, CalNP, CalRX, CalRY, CalRZ, CalH, CalF, CalE, CalEps, CalRho, CalEP, CalSig, CalN90, CalP90, CalDR90};
0091
0092
0093   // event histogram accessors
0094   enum EType {EvtAll};
0095   enum E1D   {EvtNTrk, EvtNClust, EvtParEne, EvtLeadEne, EvtTotEne, EvtLeadDivTot, EvtLeadDivPar, EvtTotDivPar, EvtNClAb90, EvtNClBe90, EvtNClRec};
0096
0097   // projection histogram accessors
0098   enum TType {TrkAll, TrkMatch};
0099   enum T1D   {TrkPosX, TrkPosY, TrkPosZ, TrkEta, TrkPhi, TrkMom, TrkDeltaR, TrkDeltaRScale, TrkEOverP, TrkDeposit};
0100   enum T2D   {TrkPosYvsX, TrkMomVsEta, TrkEtaVsPhi, TrkEPvsDR, TrkEPvsDRSig};
0101
0102   // cluster histogram accessors
0103   enum CType {CalAll, CalAbove90, CalBelow90};
0104   enum C1D   {CalNHit, CalNProj, CalPosX, CalPosY, CalPosZ, CalEta, CalPhi, CalEne, CalFraction, CalPurity, CalEOverP, CalSigma, CalNAdd90, CalPAdd90, CalDRAdd90};
0105   enum C2D   {CalPosYvsX, CalEneVsEta, CalEneVsPhi, CalEtaVsPhi, CalPurVsEne, CalN90VsEP, CalP90VsEP, CalDR90VsEP};
0106
0107   // histogram content definitions ------------------------------------------
0108
0109   // event histogram content
0110   struct EvtContent {
0111
0112     // data
0113     int64_t nTrk   = 0;
0114     int64_t nClust = 0;
0115     int64_t nCl90  = 0;
0116     int64_t nClB90 = 0;
0117     int64_t nClRec = 0;
0118     double  ePar   = 0.;
0119     double  eLead  = 0.;
0120     double  eTot   = 0.;
0121     double  lOverT = 0.;
0122     double  lOverP = 0.;
0123     double  tOverP = 0.;
0124
0125   };  // end EvtContent
0126
0127   // projection histogram content
0128   struct TrkContent {
0129
0130     // data
0131     double rx      = 0.;
0132     double ry      = 0.;
0133     double rz      = 0.;
0134     double eta     = 0.;
0135     double phi     = 0.;
0136     double p       = 0.;
0137     double dr      = 0.;
0138     double drSig   = 0.;
0139     double eOverP  = 0.;
0140     double deposit = 0.;
0141
0142     // get info from projection
0143     void GetInfo(edm4eic::TrackPoint& point) {
0144       rx  = point.position.x;
0145       ry  = point.position.y;
0146       rz  = point.position.z;
0147       eta = edm4hep::utils::eta( point.position );
0148       phi = atan2( point.position.y, point.position.x );
0149       p   = edm4hep::utils::magnitude( point.momentum );
0150       return;
0151     }  // end 'GetInfo(edm4eic::TrackPoint&)'
0152   };  // end TrkContent
0153
0154   // cluster histogram content
0155   struct CalContent {
0156
0157     // data
0158     int64_t nHit     = 0;
0159     int64_t nProj    = 0;
0160     int64_t nAdd90   = 0;
0161     double  rx       = 0.;
0162     double  ry       = 0.;
0163     double  rz       = 0.;
0164     double  eta      = 0.;
0165     double  phi      = 0.;
0166     double  ene      = 0.;
0167     double  frac     = 0.;
0168     double  purity   = 0.;
0169     double  eOverP   = 0.;
0170     double  sigma    = 0.;
0171     double  perAdd90 = 0.;
0172     double  drAdd90  = 0.;
0173
0174     // get info from cluster
0175     void GetInfo(edm4eic::Cluster& cluster) {
0176       nHit = cluster.getHits().size();
0177       rx   = cluster.getPosition().x;
0178       ry   = cluster.getPosition().y;
0179       rz   = cluster.getPosition().z;
0180       eta  = edm4hep::utils::eta( cluster.getPosition() );
0181       phi  = atan2( cluster.getPosition().y, cluster.getPosition().x );
0182       ene  = cluster.getEnergy();
0183       return;
0184     }  // end 'GetInfo(edm4eic::Cluster&)'
0185   };  // end CalContent
0186
0187 }  // end Hist namespace
0188
0189
0190
0191 // CheckClusterSplittingProcessor definition ----------------------------------
0192
0193 class CheckClusterSplittingProcessor : public JEventProcessorSequentialRoot {
0194
0195   // struct to hold user options
0196   struct Config {
0197     float  drStep;
0198     float  drMax;
0199     string parCollect;
0200     string projCollect;
0201     vector<pair<string, string>> vecCalos;
0202     vector<double> vecAvgEP;
0203     vector<double> vecSigEP;
0204   } m_config = {
0205     .drStep      = 0.05,
0206     .drMax       = 3.0,
0207     .parCollect  = "GeneratedParticles",
0208     .projCollect = "CalorimeterTrackProjections",
0209     .vecCalos    = {
0210       make_pair("HcalEndcapNClusters", "HcalEndcapN"),
0211       make_pair("EcalEndcapNClusters", "EcalEndcapN"),
0212       make_pair("HcalBarrelClusters",  "HcalBarrel"),
0213       make_pair("EcalBarrelClusters",  "EcalBarrelImaging"),
0214       make_pair("LFHCALClusters",      "LFHCAL"),
0215       make_pair("EcalEndcapPClusters", "EcalEndcapP")
0216     },
0217     .vecAvgEP = {
0218       1.0,
0219       1.0,
0220       1.0,
0221       1.0,
0222       1.0,
0223       1.0
0224     },
0225     .vecSigEP = {
0226       1.0,
0227       1.0,
0228       1.0,
0229       1.0,
0230       1.0,
0231       1.0
0232     }
0233   };  // end Config
0234
0235   public:
0236
0237     // ctor
0238     CheckClusterSplittingProcessor() { SetTypeName(NAME_OF_THIS); }
0239
0240     // public methods
0241     void InitWithGlobalRootLock() override;
0242     void ProcessSequential(const shared_ptr<const JEvent>& event) override;
0243     void FinishWithGlobalRootLock() override;
0244
0245   private:
0246
0247     // private methods
0248     void  BuildEvtHistograms();
0249     void  BuildTrkHistograms();
0250     void  BuildCalHistograms();
0251     void  FillEvtHistograms(const int calo, const int type, const Hist::EvtContent& content);
0252     void  FillTrkHistograms(const int calo, const int type, const Hist::TrkContent& content);
0253     void  FillCalHistograms(const int calo, const int type, const Hist::CalContent& content);
0254     void  SaveHistograms();
0255     void  TryToRecover(const double eStart, const double ePar, bool& didRecover, float& drRecover, int32_t& nRecover);
0256     void  GetClustersAndDr(const edm4eic::ClusterCollection& others, const Hist::CalContent& reference, const int idRef);
0257     void  GetProjections(const edm4eic::TrackSegmentCollection& projections, const int calo);
0258     int   GetCaloID(const int iCalo);
0259     FPair GetClusterWidths(const edm4eic::Cluster& cluster);
0260     FPair GetEtaRange(const edm4eic::Cluster& cluster);
0261     FPair GetPhiRange(const edm4eic::Cluster& cluster);
0262
0263     // output histograms
0264     TH1Vec m_vecEvtH1D;
0265     TH1Vec m_vecTrkH1D;
0266     TH1Vec m_vecCalH1D;
0267     TH2Vec m_vecTrkH2D;
0268     TH2Vec m_vecCalH2D;
0269
0270     // for matching tracks to clusters
0271     TrkVec m_vecTrkProj;
0272
0273     // for adding clusters during split recovery
0274     map<int, bool>   m_mapClustToChecked;
0275     map<int, double> m_mapClustToEne;
0276     map<int, double> m_mapClustToDr;
0277
0278 };  // end CheckClusterSplittingProcessor
0279
0280 #endif
0281
0282 // end ------------------------------------------------------------------------