Merge pull request #376 from RobMina/codex/crv-pulse-hit-index

CRV pulse hit index to replace pulses vector

Author: AndrewEdmonds11

doc/branches.md

@@ -112,10 +112,13 @@ Calorimeter hits store crystal position in (x,y,z) with frame origin at tracker
 The calorimeter MC branches provide MC-truth information at different levels of the calorimeter reconstruction hierarchy.
 
 **Alignment with Reconstruction Branches:**
-MC information branches are aligned with their corresponding reco branches:
+Most calorimeter MC branches are index-aligned with their reco counterparts (same size and element index). One important exception:
+
+> **Warning — `calohitsmc` is not aligned with `calohits`.** Reco hits come from `CaloHitMaker`; MC hits come from `compressRecoMCs`, which repacks `CaloHitMC` objects when copying cluster MC truth. The two branches can differ in size and order. Offline pairs reco and MC hits via `CaloHitMCTruthAssn` at reconstruction time, but that association is not kept in production art files read by EventNtuple. Each `calohitsmc` element therefore carries `caloHitIdx_` (index into `calohits`, or `-1`). RooUtil resolves MC hits through `caloHitIdx_`; do **not** assume `calohitsmc[i]` corresponds to `calohits[i]`.
+
 - `caloclustersmc` ← aligned with `caloclusters` (same size & indices)
-- `calohitsmc` ← aligned with `calohits` (same size & indices)
-- `calodigismc` ← aligned with `calodigis` (same size & indices; raw digi MC truth) (also: `calorecodigis` too)
+- `calohitsmc` ← **not** index-aligned with `calohits` (see warning above); use `caloHitIdx_` to find the reco hit
+- `calodigismc` ← aligned with `calodigis` when both branches are filled (same size & indices; raw digi MC truth)
 - `calodigisim` ← unique SimParticles from raw digis
 - `calomcsim` ← unique SimParticles from MC clusters
 
@@ -124,7 +127,7 @@ Each MC element contains `simParticleIds` to index into the corresponding `*mcsi
 | branch | structure | explanation | key fields | leaf information |
 |--------|-----------|-------------|-----------|------------------|
 | caloclustersmc |  Vector branch |   MC-truth information for calorimeter clusters (aligned with caloclusters)| `nsim` (# of sim particles), `etot` (total true energy), `tavg` (avg time), `eprimary` (primary particle energy), `tprimary` (primary time), `simParticleIds`, `simRels` (MCRelationship), `hits_`, `prel` (primary to event primary relationship) | [see CaloClusterInfoMC.hh](../inc/CaloClusterInfoMC.hh)
-| calohitsmc |  Vector branch |   MC-truth information for calorimeter hits (aligned with calohits)| `nsim`, `eDep`, `eDepG4`, `eprimary`, `tprimary`, `eDeps`, `tDeps`, `momentumIns`, `simParticleIds`, `simRels`, `caloClusterIdx_`, etc. | [see CaloHitInfoMC.hh](../inc/CaloHitInfoMC.hh)
+| calohitsmc |  Vector branch |   MC-truth for calorimeter hits (**not** index-aligned with `calohits`; use `caloHitIdx_`)| `nsim`, `eDep`, `eDepG4`, `eprimary`, `tprimary`, `eDeps`, `tDeps`, `momentumIns`, `simParticleIds`, `simRels`, `clusterIdx_`, `caloHitIdx_`, etc. | [see CaloHitInfoMC.hh](../inc/CaloHitInfoMC.hh)
 | calodigismc |  Vector branch |   MC-truth information for raw calorimeter digis (aligned with calodigis)| `nsim`, `eDep`, `eDepG4`, `eprimary`, `tprimary`, `eDeps`, `tDeps`, `momentumIns`, `simParticleIds`, `simRels`, `caloHitIdx_`, `crystalID_`, `diskID_`, `energyCorr_`, `timeCorr_`, `posX_`, `posY_` | [see CaloDigiMCInfo.hh](../inc/CaloDigiMCInfo.hh)
 | calodigisim |  Vector branch |   unique SimParticles in all raw digis (genealogy information)| SimParticle genealogy info | [see SimInfo.hh](../inc/SimInfo.hh)
 | calomcsim |  Vector branch |   unique SimParticles in all MC clusters (genealogy information)| SimParticle genealogy info | [see SimInfo.hh](../inc/SimInfo.hh)
@@ -140,6 +143,8 @@ The branch is empty if there are no CRV hit during the event.
 | crvcoincs |  Vector branch |   information about a cluster of CRV coincidence triplets| [see CrvHitInfoReco.hh](../inc/CrvHitInfoReco.hh)
 | crvcoincsmc |  Vector branch |   information about the MC track which most likely caused the CRV coincidence triplets| [see CrvHitInfoMC.hh](../inc/CrvHitInfoMC.hh)
 | crvcoincsmcplane |  Vector branch |   information about the point where the MC trajectory crosses the xz plane of CRV-T| [see CrvPlaneInfoMC.hh](../inc/CrvPlaneInfoMC.hh)
+| crvpulses |  Vector branch |   information about CRV reco pulses, including ROC/FEB/FEB channel IDs. `crvHitIndex` gives the index in `crvcoincs`, or -1 if the pulse was not clustered into a CRV hit| [see CrvPulseInfoReco.hh](../inc/CrvPulseInfoReco.hh)
+| crvpulsesmc |  Vector branch |   MC-truth information corresponding entry-by-entry to `crvpulses`| [see CrvHitInfoMC.hh](../inc/CrvHitInfoMC.hh)
 ## Trigger Branches
 
 Trigger branches store trigger decision information for each event and track hypothesis. Each trigger branch corresponds to a different trigger path or decision criterion.

fcl/prolog.fcl

@@ -361,6 +361,7 @@ EventNtupleMaker : {
     recoPulsesTag : "SelectReco"
     stepsTag : "compressRecoMCs"
     fillPulses : false
+    keepUnclusteredPulses : false
     fillDigis : false
     digisTag : "SelectReco"
     planeY : @local::crvPlaneY.CRV_T

inc/CaloHitInfoMC.hh

@@ -22,8 +22,9 @@ namespace mu2e
     std::vector<int> simParticleIds; // list of simparticle ids
     std::vector<MCRelationship> simRels; // relationship to the particle that deposited the most energy in the calo Hit
     int clusterIdx_; // Cluster index
+    int caloHitIdx_; // index into calohits branch, -1 if unset; calohitsmc is NOT index-aligned with calohits
 
-    CaloHitInfoMC() : crystalID_(0), 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), caloHitIdx_(-1){}
 
     void reset() { *this = CaloHitInfoMC(); }
 

inc/CrvHitInfoReco.hh

@@ -3,7 +3,6 @@
 
 // CrvHitInfoReco: information about a cluster of CRV coincidence triplets
 
-#include "EventNtuple/inc/CrvPulseInfoReco.hh"
 #include "EventNtuple/inc/RootVectors.hh"
 #include "Offline/DataProducts/inc/CRVId.hh"
 namespace CLHEP {class Hep3Vector; }
@@ -30,10 +29,9 @@ namespace mu2e
     int               nHits = -1;      //number of coincidence hits in this cluster
     int               nLayers = -1;      //number of coincidence layers in this cluster
     float             angle = -999;   //coincidence direction
-    std::vector<CrvPulseInfoReco> crvPulses;
 
     CrvHitInfoReco(){}
-    CrvHitInfoReco(int sectorType, bool hitPosAndTimeCalculated, CLHEP::Hep3Vector hpos, float timeWindowStart, float timeWindowEnd, float timeAvg, float PEs, std::array<float, CRVId::nLayers> PEsPerLayer, std::array<float, CRVId::nLayers * CRVId::nSidesPerBar> sidePEsPerLayer, std::array<size_t, CRVId::nSidesPerBar> sidePulses, std::array<float, CRVId::nSidesPerBar> sidePEs, std::array<float, CRVId::nSidesPerBar> sideTimes, int nCoincidenceHits, int nCoincidenceLayers, float coincidenceAngle, std::vector<CrvPulseInfoReco> crvPulses);
+    CrvHitInfoReco(int sectorType, bool hitPosAndTimeCalculated, CLHEP::Hep3Vector hpos, float timeWindowStart, float timeWindowEnd, float timeAvg, float PEs, std::array<float, CRVId::nLayers> PEsPerLayer, std::array<float, CRVId::nLayers * CRVId::nSidesPerBar> sidePEsPerLayer, std::array<size_t, CRVId::nSidesPerBar> sidePulses, std::array<float, CRVId::nSidesPerBar> sidePEs, std::array<float, CRVId::nSidesPerBar> sideTimes, int nCoincidenceHits, int nCoincidenceLayers, float coincidenceAngle);
   };
   typedef std::vector<CrvHitInfoReco> CrvHitInfoRecoCollection;  //this is the reco vector which will be stored in the main TTree
 }

inc/CrvInfoHelper.hh

@@ -42,10 +42,17 @@ namespace mu2e
           CrvPlaneInfoMCCollection &MCInfoPlane, double crvPlaneY,
           art::Handle<PrimaryParticle> const& primary);
 
+      void FillCrvPulseHitIndices(
+          art::Handle<CrvCoincidenceClusterCollection> const& crvCoincidences,
+          art::Handle<CrvRecoPulseCollection> const& crvRecoPulses,
+          std::vector<int> &pulseHitIndices);
+
       void FillCrvPulseInfoCollections(
           art::Handle<CrvRecoPulseCollection> const& crvRecoPulses,
           art::Handle<CrvDigiMCCollection> const& crvDigiMCs,
           art::Handle<EventWindowMarker> const& ewmh,
+          const std::vector<int> &pulseHitIndices,
+          bool keepUnclusteredPulses,
           CrvPulseInfoRecoCollection &recoInfo, CrvHitInfoMCCollection &MCInfo);
 
       void FillCrvDigiInfoCollections(

inc/CrvPulseInfoReco.hh

@@ -15,15 +15,19 @@ namespace mu2e
     int    barId = -1;           //CRV counter ID
     int    sectorId = -1;        //CRV sector ID
     int    SiPMId = -1;          //SiPMId number
+    int    ROC = -1;             //ROC number
+    int    FEB = -1;             //FEB number
+    int    FEBchannel = -1;      //FEB channel number
     int    PEs = -1;             //PEs using pulse integral
     int    PEsPulseHeight = -1;  //PEs using pulse height
     float  pulseHeight = -1;     //Pulse height
     float  pulseBeta = -1;       //Pulse beta
     float  pulseFitChi2 = -1;    //Pulse Fit chi2
     float  time = -1;            //Time
+    int    crvHitIndex = -1;     //index of the CRV hit containing this pulse, or -1 if unclustered
 
     CrvPulseInfoReco(){}
-    CrvPulseInfoReco(CLHEP::Hep3Vector ppos, int barId, int sectorId, int SiPMId, int PEs, int PEsPulseHeight, float pulseHeight, float pulseBeta, float pulseFitChi2, float time);
+    CrvPulseInfoReco(CLHEP::Hep3Vector ppos, int barId, int sectorId, int SiPMId, int ROC, int FEB, int FEBchannel, int PEs, int PEsPulseHeight, float pulseHeight, float pulseBeta, float pulseFitChi2, float time, int crvHitIndex);
   };
 
   typedef std::vector<CrvPulseInfoReco> CrvPulseInfoRecoCollection;  //this is the reco vector which will be stored in the main TTree

rooutil/examples/PrintEvents.C

@@ -212,7 +212,7 @@ void PrintEvents(std::string filename) {
     // crvpulses branch
     if (event.crvpulses != nullptr) {
       for (const auto& crvpulse : *(event.crvpulses)) {
-        std::cout << "crvpulse: " << crvpulse.pos.x() << "," << crvpulse.barId << "," << crvpulse.time << std::endl;
+        std::cout << "crvpulse: " << crvpulse.pos.x() << "," << crvpulse.barId << "," << crvpulse.ROC << "," << crvpulse.FEB << "," << crvpulse.FEBchannel << "," << crvpulse.time << "," << crvpulse.crvHitIndex << std::endl;
       }
     }
 

rooutil/inc/Event.hh

@@ -36,6 +36,7 @@
 #include "EventNtuple/inc/TrkStrawHitCalibInfo.hh"
 
 #include "EventNtuple/inc/MVAResultInfo.hh"
+#include <stdexcept>
 
 #include "EventNtuple/inc/TrigInfo.hh"
 #include "EventNtuple/inc/MCStepInfo.hh"
@@ -155,6 +156,8 @@ namespace rooutil {
       for (int i_crv_coinc = 0; i_crv_coinc < nCrvCoincs(); ++i_crv_coinc) {
         CrvCoinc crv_coinc(&(crvcoincs->at(i_crv_coinc)));
         if (crvcoincsmc != nullptr) {
+          if (i_crv_coinc >= static_cast<int>(crvcoincsmc->size()))
+            throw std::runtime_error("CRV coincidence cluster count mismatch between reco and MC collections");
           crv_coinc.mc = &(crvcoincsmc->at(i_crv_coinc));
         }
         crv_coincs.emplace_back(crv_coinc);
@@ -185,8 +188,15 @@ namespace rooutil {
               CaloHit calohit;
               calohit.reco = &(calohits->at(calohit_Idx)); // passing the addresses of the underlying structs
 
-              if (calohitsmc != nullptr) { // 1-to-1 matching of calohits and calohitsmc
-                calohit.mc = &(calohitsmc->at(calohit_Idx));
+              // WARNING: calohitsmc is NOT index-aligned with calohits (different art products).
+              // Do not use calohitsmc->at(calohit_Idx); look up via caloHitIdx_ instead.
+              if (calohitsmc != nullptr) {
+                for (int i_mc = 0; i_mc < static_cast<int>(calohitsmc->size()); ++i_mc) {
+                  if (calohitsmc->at(i_mc).caloHitIdx_ == static_cast<int>(calohit_Idx)) {
+                    calohit.mc = &(calohitsmc->at(i_mc));
+                    break;
+                  }
+                }
               }
 
               calo_cluster.hits.emplace_back(calohit);

src/CrvHitInfoReco.cc

@@ -1,7 +1,7 @@
 #include "EventNtuple/inc/CrvHitInfoReco.hh"
 #include "CLHEP/Vector/ThreeVector.h"
 namespace mu2e {
-  CrvHitInfoReco::CrvHitInfoReco(int sectorType, bool hitPosAndTimeCalculated, CLHEP::Hep3Vector hpos, float timeWindowStart, float timeWindowEnd, float timeAvg, float PEs, std::array<float, CRVId::nLayers> PEsPerLayer, std::array<float, CRVId::nLayers * CRVId::nSidesPerBar> sidePEsPerLayer, std::array<size_t, CRVId::nSidesPerBar> sidePulses, std::array<float, CRVId::nSidesPerBar> sidePEs, std::array<float, CRVId::nSidesPerBar> sideTimes, int nCoincidenceHits, int nCoincidenceLayers, float coincidenceAngle, CrvPulseInfoRecoCollection crvPulses) :
+  CrvHitInfoReco::CrvHitInfoReco(int sectorType, bool hitPosAndTimeCalculated, CLHEP::Hep3Vector hpos, float timeWindowStart, float timeWindowEnd, float timeAvg, float PEs, std::array<float, CRVId::nLayers> PEsPerLayer, std::array<float, CRVId::nLayers * CRVId::nSidesPerBar> sidePEsPerLayer, std::array<size_t, CRVId::nSidesPerBar> sidePulses, std::array<float, CRVId::nSidesPerBar> sidePEs, std::array<float, CRVId::nSidesPerBar> sideTimes, int nCoincidenceHits, int nCoincidenceLayers, float coincidenceAngle) :
     sectorType(sectorType),
     hitPosAndTimeCalculated(hitPosAndTimeCalculated),
     pos(hpos),
@@ -16,7 +16,6 @@ namespace mu2e {
     sideTimes(sideTimes),
     nHits(nCoincidenceHits),
     nLayers(nCoincidenceLayers),
-    angle(coincidenceAngle),
-    crvPulses(crvPulses)
+    angle(coincidenceAngle)
   {}
 }