Merge pull request #114 from brownd1978/simhits

Add content to MCSim branch vector
Mu2e · Nov 22, 2023 · 40522cf · 40522cf
2 parents a4905fc + 38c6bfa
commit 40522cf
Show file tree

Hide file tree

Showing 6 changed files with 89 additions and 48 deletions.
diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl
@@ -97,13 +97,15 @@ dioLLWeight : {
 }
 genCountLogger: { module_type: GenEventCountReader }
 
-AllOpt : { fillMC : true
+AllOpt : {
+  fillMC : true
   trkqual : "TrkQual"
   fillTrkQual : true
   trkpid : "TrkPID"
   fillTrkPID : true
   fillHits : true
   genealogyDepth : -1
+  matchDepth : -1
 }
 
 DeM : { input : "KK"

diff --git a/fcl/prolog_trigger.fcl b/fcl/prolog_trigger.fcl
@@ -6,13 +6,6 @@ BEGIN_PROLOG
 TrkAnaTrigger : {
   TrkAnaTT : {
     @table::TrkAnaTreeMaker
-      candidate : {
-        @table::DeM
-          options : {
-            @table::AllOpt
-            fillHits : true
-          }
-      }
     supplements : []
     ExtraMCStepCollectionTags : []
     diagLevel : 2
@@ -37,5 +30,33 @@ TrkAnaTrigger : {
     }
   }
 }
+TrkAnaTrigger : {
+  @table::TrkAnaTrigger
+  TrkAnaDeMTT : {
+    @table::TrkAnaTrigger.TrkAnaTT
+    candidate : {
+      branch: "trk"
+      suffix : "DeM"
+      options : {
+        @table::AllOpt
+        fillHits : true
+      }
+    }
+  }
+
+  TrkAnaDePTT : {
+    @table::TrkAnaTrigger.TrkAnaTT
+    candidate : {
+      branch: "trk"
+      suffix : "DeP"
+      options : {
+        @table::AllOpt
+        fillHits : true
+      }
+    }
+  }
+
+
+}
 
 END_PROLOG
diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh
@@ -60,7 +60,7 @@ namespace mu2e {
       void fillTrkInfoMC(const KalSeed& kseed, const KalSeedMC& kseedmc, TrkInfoMC& trkinfomc);
       void fillTrkInfoMCDigis(const KalSeed& kseed, const KalSeedMC& kseedmc, TrkInfoMC& trkinfomc);
       void fillHitInfoMC(const KalSeedMC& kseedmc, TrkStrawHitInfoMC& tshinfomc, const TrkStrawHitMC& tshmc);
-      void fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector<SimInfo>& siminfos, int n_generations);
+      void fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector<SimInfo>& siminfos, int n_generations, int n_match);
       void fillVDInfo(KalSeed const& kseed, const KalSeedMC& kseedmc, std::vector<MCStepInfo>& vdinfos);
       void fillHitInfoMCs(const KalSeedMC& kseedmc, std::vector<TrkStrawHitInfoMC>& tshinfomcs);
       void fillCaloClusterInfoMC(CaloClusterMC const& ccmc, CaloClusterInfoMC& ccimc);

diff --git a/inc/SimInfo.hh b/inc/SimInfo.hh
@@ -13,13 +13,18 @@ namespace mu2e
 // general info about the SimParticle which was simulated
   struct SimInfo {
     bool valid = false; // true/false as to whether the data is valid
+    Int_t nhits = -1; // number of matched hits in the track generated by this particle
+    Int_t nactive = -1; // number of active matched hits in the track generated by this particle
+    Int_t rank = -1; // rank of this particle, by # of associated hits
     Int_t pdg = -1; // true PDG code
+    Int_t proc = -1; // process code
     Int_t gen = -1; // generator code
-    Float_t time = -1.0;  // time of this step
-    XYZVectorF mom = XYZVectorF();  // particle momentum at the start of this step
-    XYZVectorF pos = XYZVectorF();  // particle position at the start of this step
+    Float_t time = -1.0;  // Origin time of the SimParticle
+    Int_t index; // index into the SimInfo vector
+    XYZVectorF mom = XYZVectorF();  // origin momentumof the SimParticle
+    XYZVectorF pos = XYZVectorF();  // origin position of the SimParticle
     MCRelationship prirel = MCRelationship::none; // relationship to the event primary particles
-    MCRelationship trkrel = MCRelationship::none; // relationship to the particle that created the track
+    MCRelationship trkrel = MCRelationship::none; // relationship to the particle that created hits in the track
     void reset() { *this = SimInfo(); }
   };
 }

diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc
@@ -123,39 +123,51 @@ namespace mu2e {
     tshinfomc.doca = -1*dperp;
   }
 
-  void InfoMCStructHelper::fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector<SimInfo>& siminfos, int n_generations) {
-    auto trkprimaryptr = kseedmc.simParticle().simParticle(_spcH);
-    auto trkprimary = trkprimaryptr->originParticle();
-
-    auto current_sim_particle_ptr = trkprimaryptr;
-    auto current_sim_particle = trkprimary;
-    if (n_generations == -1) { // means do all generations
+  void InfoMCStructHelper::fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector<SimInfo>& siminfos, int n_generations, int n_match) {
+    // interpret -1 as no llimit
+    if (n_generations == -1) {
       n_generations = std::numeric_limits<int>::max();
     }
-
-    for (int i_generation = 0; i_generation < n_generations; ++i_generation) {
-      SimInfo sim_info;
-      fillSimInfo(current_sim_particle, sim_info);
-      sim_info.trkrel = MCRelationship(current_sim_particle_ptr, trkprimaryptr);
-
-      auto bestprimarysp = primary.primarySimParticles().front();
-      MCRelationship bestrel;
-      for(auto const& spp : primary.primarySimParticles()){
-        MCRelationship mcrel(current_sim_particle_ptr, spp);
-        if(mcrel > bestrel){
-          bestrel = mcrel;
-          bestprimarysp = spp;
+    if (n_match == -1 ) {
+      n_match = std::numeric_limits<int>::max();
+    }
+    for(int imatch = 0 ; imatch < std::min(n_match,static_cast<int>(kseedmc.simParticles().size())); ++imatch) {
+      auto trkprimaryptr = kseedmc.simParticle(imatch).simParticle(_spcH);
+      auto trkprimary = trkprimaryptr->originParticle();
+      auto current_sim_particle_ptr = trkprimaryptr;
+      auto current_sim_particle = trkprimary;
+
+      for (int i_generation = 0; i_generation < n_generations; ++i_generation) {
+        SimInfo sim_info;
+        fillSimInfo(current_sim_particle, sim_info);
+        sim_info.trkrel = MCRelationship(current_sim_particle_ptr, trkprimaryptr);
+        sim_info.rank = imatch;
+
+        auto bestprimarysp = primary.primarySimParticles().front();
+        MCRelationship bestrel;
+        for(auto const& spp : primary.primarySimParticles()){
+          MCRelationship mcrel(current_sim_particle_ptr, spp);
+          if(mcrel > bestrel){
+            bestrel = mcrel;
+            bestprimarysp = spp;
+          }
+        }
+        sim_info.prirel = bestrel;
+        // only count hits for direct contributors
+        if(i_generation == 0){
+          sim_info.nhits = kseedmc.simParticle(imatch)._nhits;
+          sim_info.nactive = kseedmc.simParticle(imatch)._nactive;
+        }
+        // record the index this object will have
+        sim_info.index = siminfos.size();
+        siminfos.push_back(sim_info);
+        if (current_sim_particle.parent().isNonnull()) {
+          current_sim_particle_ptr = current_sim_particle.parent();
+          current_sim_particle = current_sim_particle_ptr->originParticle();
+        }
+        else {
+          break; // this particle doesn't have a parent
         }
-      }
-      sim_info.prirel = bestrel;
-
-      siminfos.push_back(sim_info);
-      if (current_sim_particle.parent().isNonnull()) {
-        current_sim_particle_ptr = current_sim_particle.parent();
-        current_sim_particle = current_sim_particle_ptr->originParticle();
-      }
-      else {
-        break; // this particle doesn't have a parent
       }
     }
 
@@ -173,12 +185,13 @@ namespace mu2e {
         }
       }
       if (!already_added) {
+        auto trkprimaryptr = kseedmc.simParticle().simParticle(_spcH);
         sim_info.trkrel = MCRelationship(spp, trkprimaryptr);
         sim_info.prirel = MCRelationship(spp, spp);
+        sim_info.index = siminfos.size();
         siminfos.push_back(sim_info);
       }
     }
-
   }
 
 
@@ -192,12 +205,11 @@ namespace mu2e {
   }
 
   void InfoMCStructHelper::fillSimInfo(const SimParticle& sp, SimInfo& siminfo) {
-
     GeomHandle<DetectorSystem> det;
-
     siminfo.valid = true;
+    if(sp.genParticle().isNonnull())siminfo.gen = sp.genParticle()->generatorId().id();
+    siminfo.proc = sp.creationCode();
     siminfo.pdg = sp.pdgId();
-    siminfo.gen = sp.creationCode();
     siminfo.time = sp.startGlobalTime();
     siminfo.mom = XYZVectorF(sp.startMomentum());
     siminfo.pos = XYZVectorF(det->toDetector(sp.startPosition()));

diff --git a/src/TrkAnaTreeMaker_module.cc b/src/TrkAnaTreeMaker_module.cc
@@ -113,6 +113,7 @@ namespace mu2e {
         fhicl::Atom<bool> filltrkpid{Name("fillTrkPID"), Comment("Switch to turn on filling of the full TrkPIDInfo for this set of tracks"), false};
         fhicl::Atom<bool> required{Name("required"), Comment("True/false if you require this type of track in the event"), false};
         fhicl::Atom<int> genealogyDepth{Name("genealogyDepth"), Comment("The depth of the genealogy information you want to keep"), 1};
+        fhicl::Atom<int> matchDepth{Name("matchDepth"), Comment("The depth into the MC true particle matching you want to keep"), 1};
       };
 
       struct BranchConfig {
@@ -854,7 +855,7 @@ namespace mu2e {
           _infoMCStructHelper.fillTrkInfoMC(kseed, kseedmc, _allMCTIs.at(i_branch));
           auto& mcvdis = _allMCVDInfos.at(i_branch);
           _infoMCStructHelper.fillVDInfo(kseed, kseedmc, mcvdis);
-          _infoMCStructHelper.fillAllSimInfos(kseedmc, primary, _allMCSimTIs.at(i_branch), branchConfig.options().genealogyDepth());
+          _infoMCStructHelper.fillAllSimInfos(kseedmc, primary, _allMCSimTIs.at(i_branch), branchConfig.options().genealogyDepth(), branchConfig.options().matchDepth());
 
           if(_conf.diag() > 1 || (_conf.fillhits() && branchConfig.options().fillhits())){
             _infoMCStructHelper.fillHitInfoMCs(kseedmc, _allTSHIMCs.at(i_branch));