Merge branch 'main' into main

Author: ehrlich-uva

fcl/prolog.fcl

@@ -286,7 +286,7 @@ EventNtupleMaker : {
   FillTimeClusterInfo : false
   CaloClustersTag : "CaloClusterMaker"
   CaloHitsTag : "CaloHitMaker"
-  CaloRecoDigisTag : ""
+  CaloRecoDigisTag : "CaloRecoDigiMaker"
   CaloDigisTag : "SelectReco"
   FillCaloClusters : true
   FillCaloHits : true
@@ -315,6 +315,7 @@ EventNtupleMaker : {
   PrimaryParticleTag : "compressRecoMCs"
   KalSeedMCAssns : "SelectReco"
   CaloClusterMCTag : "compressRecoMCs"
+  CaloHitMCTag : "compressRecoMCs"
   InfoMCStructHelper : {
     SimParticleCollectionTag : "compressRecoMCs"
     MinGoodMomFraction : 0.9

helper/ntuplehelper.py

@@ -246,12 +246,11 @@ def list_all_branches(self, export_to_md=False):
         else:
             print("## Calorimeter Branches\n")
             print("These branches are vectors of calorimeter clusters/hits/recodigis/digis that happened during the event.")
-            print("The branch is empty if there are no calo cluster during the event.\n")
             print("While each branch can be read independently, i.e. all the hits of the event, each element contains indexes to the other branches for parentage link.")
             print("The cluster element contains the vector \'hits_\' containing the indexes of the hits branch belonging to this cluster.")
             print("Similarly, each hit contains the indexes of its two recodigis (left and right channels) and the index of its parent cluster.")
             print("Example: cluster 3 has hits_ = {12, 13, 14, 15}. Each of those hits will have \'clusterIdx_\' = 3.")
-            print("Calo hits have the crystal (x,y,z) position saved. Frame origin is center of tracker.")
+            print("Calo digis and hits have the crystal (x,y,z) position saved. Frame origin is center of tracker.")
             print("| branch | structure | explanation | leaf information |")
             print("|--------|-----------|-------------|------------------|")
         for branch in self.calo_branches:

inc/CaloClusterInfoMC.hh

@@ -7,6 +7,7 @@
 namespace mu2e
 {
   struct CaloClusterInfoMC {
+    int nhits; // number of hitmc associated with this cluster
     int nsim; // # of sim particles associated with this cluster
     float etot; // total true energy from all particles in this cluster
     float tavg; // average time over all particles

inc/CaloDigiInfo.hh

@@ -5,22 +5,32 @@
 #define CaloDigiInfo_HH
 
 #include <vector>
+#include "Offline/RecoDataProducts/inc/CaloDigi.hh"
 
 namespace mu2e
 {
   struct CaloDigiInfo {
 
     int                 SiPMID_;           // Offline SiPM ID number
-    int                 diskID_;           // Offline disk ID number
     int                 t0_;               // Time of first waveform sample
     std::vector<int>    waveform_;         // Waveform
     int                 peakpos_;          // Waveform index of peak
+
+    // not in dataproduct
+    int                 diskID_;           // Offline disk ID number
     int                 caloRecoDigiIdx_;  // RecoDigi index
-    float               posX_;             // Crystal x position in detector frame
-    float               posY_;             // Crystal y position in detector frame
+    int                 peakval_;          // Waveform value at peak
+    XYZVectorF          crystalPos_;       // Crystal position
 
-    CaloDigiInfo() : SiPMID_(-1), t0_(0), waveform_(), peakpos_(0), caloRecoDigiIdx_(-1), posX_(0.0), posY_(0.0) {}
+    CaloDigiInfo() : SiPMID_(-1), t0_(0), waveform_(), peakpos_(0), diskID_(-1), caloRecoDigiIdx_(-1), peakval_(0), crystalPos_() {}
     void reset() { *this = CaloDigiInfo(); }
+    
+    bool operator== (const CaloDigi& other) const {
+      return SiPMID_          == other.SiPMID()  &&
+             t0_              == other.t0()      &&
+             peakpos_         == other.peakpos() &&
+             waveform_.size() == other.waveform().size(); //no need to check the entire waveform
+    }
   };
 }
 #endif

inc/CaloHitInfo.hh

@@ -6,6 +6,7 @@
 
 #include <vector>
 #include "EventNtuple/inc/RootVectors.hh"
+#include "Offline/RecoDataProducts/inc/CaloHit.hh"
 
 namespace mu2e
 {
@@ -18,12 +19,23 @@ namespace mu2e
     float             eDep_;                // Hit energy
     float             eDepErr_;             // Hit energy error
     std::vector<int>  recoDigis_;           // vector of branch indices of reco digis
-    int               clusterIdx_;          // Cluster index
 
+    // not in dataproduct
+    int               clusterIdx_;          // Cluster index
     XYZVectorF        crystalPos_;          // Crystal position
 
     CaloHitInfo() : crystalId_(0), nSiPMs_(0), time_(0.0), timeErr_(0.0), eDep_(0.0), eDepErr_(0.0), recoDigis_(), clusterIdx_(-1), crystalPos_() {}
     void reset() { *this = CaloHitInfo(); }
+
+    bool operator== (const CaloHit& other) const {
+      const double eps = 1e-9;
+      return crystalId_       == other.crystalID()        &&
+             nSiPMs_          == other.nSiPMs()           &&
+             std::abs(time_    - other.time()) < eps      &&
+             std::abs(timeErr_ - other.timeErr()) < eps   &&
+             std::abs(eDep_    - other.energyDep()) < eps &&
+             std::abs(eDepErr_ - other.energyDepErr()) < eps;
+    }
   };
 }
 #endif

inc/CaloHitInfoMC.hh

@@ -4,10 +4,13 @@
 #ifndef CaloHitInfoMC_HH
 #define CaloHitInfoMC_HH
 #include "Offline/MCDataProducts/inc/MCRelationship.hh"
+#include "Offline/MCDataProducts/inc/CaloHitMC.hh"
+
 namespace mu2e
 {
   struct CaloHitInfoMC {
 
+    int crystalID_; // Crystal ID of the hit
     int nsim; // # of sim particles associated with this hit
     float eDep; // total (corrected) deposited energy in the hit
     float eDepG4; // total G4 deposited energy in the hit
@@ -20,9 +23,25 @@ namespace mu2e
     std::vector<MCRelationship> simRels; // relationship to the particle that deposited the most energy in the calo Hit
     int clusterIdx_; // Cluster index
 
-    CaloHitInfoMC() : nsim(0), eDep(0.0), eDepG4(0.0), eprimary(0.0), tprimary(0.0), clusterIdx_(-1){}
+    CaloHitInfoMC() : crystalID_(0), nsim(0), eDep(0.0), eDepG4(0.0), eprimary(0.0), tprimary(0.0), clusterIdx_(-1){}
 
     void reset() { *this = CaloHitInfoMC(); }
+
+    bool operator== (const CaloHitMC& other) const {
+      const double eps = 1e-9;
+      auto const& edeps = other.energyDeposits();
+      bool same1 = crystalID_     == other.crystalID()             &&
+                   nsim           == int(edeps.size())             &&
+                   std::abs(eDep   - other.totalEnergyDep()) < eps &&
+                   std::abs(eDepG4 - other.totalEnergyDepG4()) < eps;
+      if (nsim > 0){ //Only compare first edep if exists
+        return same1                                                &&
+               std::abs(eprimary - edeps.front().energyDep()) < eps &&
+               std::abs(tprimary - edeps.front().time()) < eps;
+      } else {
+        return same1;
+      }
+    }
   };
 }
 #endif

inc/CaloRecoDigiInfo.hh

@@ -7,6 +7,7 @@
 #include "CLHEP/Vector/ThreeVector.h"
 #include "CLHEP/Units/PhysicalConstants.h"
 #include <vector>
+#include "Offline/RecoDataProducts/inc/CaloRecoDigi.hh"
 
 namespace mu2e
 {
@@ -20,10 +21,24 @@ namespace mu2e
     unsigned            ndf_;              // Fit NDF
     bool                pileUp_;           // Flag indicating if the hit is affected by pile-up
     int                 caloDigiIdx_;      // Index to the associated CaloDigi
+
+    // not in dataproduct
     int                 caloHitIdx_;       // Hit index
+    int                 SiPMID_;           // Offline SiPM ID number 
 
-    CaloRecoDigiInfo() : eDep_(0.0), eDepErr_(0.0), time_(0.0), timeErr_(0.0), chi2_(0.0), ndf_(0), pileUp_(false), caloDigiIdx_(-1), caloHitIdx_(-1) {}
+    CaloRecoDigiInfo() : eDep_(0.0), eDepErr_(0.0), time_(0.0), timeErr_(0.0), chi2_(0.0), ndf_(0), pileUp_(false), caloDigiIdx_(-1), caloHitIdx_(-1), SiPMID_(0) {}
     void reset() { *this = CaloRecoDigiInfo(); }
+    
+    bool operator== (const CaloRecoDigi& other) const {
+      const double eps = 1e-9;
+      return std::abs(eDep_    - other.energyDep())    < eps &&
+             std::abs(eDepErr_ - other.energyDepErr()) < eps &&
+             std::abs(time_    - other.time())    < eps &&
+             std::abs(timeErr_ - other.timeErr()) < eps &&
+             std::abs(chi2_    - other.chi2())    < eps &&
+             ndf_             == other.ndf()            &&
+             pileUp_          == other.pileUp();
+    }
   };
 }
 #endif

inc/InfoMCStructHelper.hh

@@ -75,7 +75,7 @@ namespace mu2e {
       void fillVDInfo(KalSeed const& kseed, const KalSeedMC& kseedmc, std::vector<std::vector<MCStepInfo>>& all_vdinfos);
       void fillHitInfoMCs(const KalSeed& kseed, const KalSeedMC& kseedmc, std::vector<std::vector<TrkStrawHitInfoMC>>& all_tshinfomcs);
       void fillCaloClusterInfoMC(CaloClusterMC const& ccmc, std::vector<CaloClusterInfoMC>& ccimc);
-      void fillCaloHitInfoMC(CaloHitMC const& chmc, std::vector<CaloHitInfoMC>& chimc, int clusterIdx);
+      void fillCaloHitInfoMC(CaloHitMC const& chmc, std::vector<CaloHitInfoMC>& chimc, int clusterIdx = -1);
       void fillCaloDigiMCInfo(CaloShowerSim const& shower, std::vector<CaloDigiMCInfo>& calodigimc);
       void fillCaloDigiSimInfos(CaloShowerSim const& shower, std::vector<SimInfo>& cdsis);
       void fillCaloSimInfos(CaloClusterMC const& ccmc, std::vector<SimInfo>& csis);