Particle Ancestry Map

larsim/MergeSimSources/MergeSimSources.cxx

@@ -201,6 +201,17 @@ void sim::MergeSimSourcesUtility::MergeAuxDetHits(std::vector<sim::AuxDetHit>& d
   std::transform(begin(src), end(src), back_inserter(dest), offsetAuxDetHitID);
 }
 
+void sim::MergeSimSourcesUtility::MergeParticleAncestryMaps(
+  std::vector<sim::ParticleAncestryMap>& dest,
+  const sim::ParticleAncestryMap& src,
+  std::size_t source_index) const
+{
+  const int offset = fG4TrackIDOffsets.at(source_index);
+
+  dest.reserve(dest.size() + 1);
+  dest.push_back(offsetParticleAncestryMapTrackID(src, offset));
+}
+
 void sim::MergeSimSourcesUtility::UpdateG4TrackIDRange(std::pair<int, int> newrange,
                                                        size_t source_index)
 {
@@ -277,3 +288,21 @@ sim::AuxDetHit sim::MergeSimSourcesUtility::offsetAuxDetHitTrackID(sim::AuxDetHi
     adh.GetExitMomentumZ(),
   };
 } // sim::MergeSimSourcesUtility::offsetAuxDetHitTrackID()
+
+sim::ParticleAncestryMap sim::MergeSimSourcesUtility::offsetParticleAncestryMapTrackID(
+  sim::ParticleAncestryMap const& pam,
+  int offset)
+{
+  std::map<int, std::set<int>> newMap;
+
+  for (auto const& [ancestor, descendants] : pam.GetMap()) {
+    const int newAnc = (ancestor >= 0) ? ancestor + offset : ancestor - offset;
+
+    for (auto const& descendant : descendants) {
+      const int newDesc = (descendant >= 0) ? descendant + offset : descendant - offset;
+      newMap[newAnc].insert(newDesc);
+    }
+  }
+
+  return sim::ParticleAncestryMap(newMap);
+}

larsim/MergeSimSources/MergeSimSources.h

@@ -10,12 +10,13 @@
  * Typically just merges vectors/maps/etc together. But, if anything as a G4 trackID, applies
  * a user-defined offset to those IDs.
  *
-*/
+ */
 
 #include "nusimdata/SimulationBase/MCParticle.h"
 
 #include "lardataobj/Simulation/AuxDetHit.h"
 #include "lardataobj/Simulation/AuxDetSimChannel.h"
+#include "lardataobj/Simulation/ParticleAncestryMap.h"
 #include "lardataobj/Simulation/SimChannel.h"
 #include "lardataobj/Simulation/SimEnergyDeposit.h"
 #include "lardataobj/Simulation/SimPhotons.h"
@@ -57,6 +58,10 @@ namespace sim {
                          const std::vector<sim::AuxDetHit>&,
                          size_t) const;
 
+    void MergeParticleAncestryMaps(std::vector<sim::ParticleAncestryMap>&,
+                                   const sim::ParticleAncestryMap&,
+                                   size_t) const;
+
     const std::vector<std::vector<size_t>>& GetMCParticleListMap() { return fMCParticleListMap; }
 
   private:
@@ -71,6 +76,10 @@ namespace sim {
 
     static sim::AuxDetHit offsetAuxDetHitTrackID(sim::AuxDetHit const&, int);
 
+    static sim::ParticleAncestryMap offsetParticleAncestryMapTrackID(
+      sim::ParticleAncestryMap const&,
+      int);
+
   }; //end MergeSimSourcesUtility class
 
 } //end namespace sim

larsim/MergeSimSources/MergeSimSources_module.cc

@@ -112,6 +112,11 @@ class sim::MergeSimSources : public art::EDProducer {
                                "Other"} // default
     };
 
+    fhicl::Atom<bool> FillParticleAncestryMaps{
+      fhicl::Name{"FillParticleAncestryMaps"},
+      fhicl::Comment{"whether to merge particle ancestry maps"},
+      true};
+
   }; // struct Config
 
   using Parameters = art::EDProducer::Table<Config>;
@@ -134,6 +139,7 @@ class sim::MergeSimSources : public art::EDProducer {
   std::vector<std::string> const fEnergyDepositionInstances;
   bool const fFillAuxDetHits;
   std::vector<std::string> const fAuxDetHitsInstanceLabels;
+  bool const fFillParticleAncestryMaps;
 
   static std::string const ReflectedLabel;
 
@@ -176,6 +182,7 @@ sim::MergeSimSources::MergeSimSources(Parameters const& params)
   , fEnergyDepositionInstances(params().EnergyDepositInstanceLabels())
   , fFillAuxDetHits(params().FillAuxDetHits())
   , fAuxDetHitsInstanceLabels(params().AuxDetHitsInstanceLabels())
+  , fFillParticleAncestryMaps(params().FillParticleAncestryMaps())
 {
 
   if (fInputSourcesLabels.size() != fTrackIDOffsets.size()) {
@@ -223,6 +230,8 @@ sim::MergeSimSources::MergeSimSources(Parameters const& params)
       }
     }
 
+    if (fFillParticleAncestryMaps) { consumes<sim::ParticleAncestryMap>(tag); }
+
   } // for input labels
 
   if (fFillMCParticles) {
@@ -256,6 +265,8 @@ sim::MergeSimSources::MergeSimSources(Parameters const& params)
       produces<std::vector<sim::AuxDetHit>>(auxdethit_inst);
   }
 
+  if (fFillParticleAncestryMaps) { produces<std::vector<sim::ParticleAncestryMap>>(); }
+
   dumpConfiguration();
 }
 
@@ -271,6 +282,7 @@ void sim::MergeSimSources::produce(art::Event& e)
   auto tpassn =
     std::make_unique<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo>>();
   auto adCol = std::make_unique<std::vector<sim::AuxDetSimChannel>>();
+  auto pamCol = std::make_unique<std::vector<sim::ParticleAncestryMap>>();
 
   using edeps_t = std::vector<sim::SimEnergyDeposit>;
   std::vector<edeps_t> edepCols;
@@ -349,6 +361,11 @@ void sim::MergeSimSources::produce(art::Event& e)
         MergeUtility.MergeAuxDetHits(
           auxdethitCol, e.getProduct<aux_det_hits_t>(auxdethit_tag), i_source);
       }
+
+      if (fFillParticleAncestryMaps) {
+        auto const& input_pamCol = e.getProduct<sim::ParticleAncestryMap>(input_label);
+        MergeUtility.MergeParticleAncestryMaps(*pamCol, input_pamCol, i_source);
+      }
     }
   }
 
@@ -383,6 +400,8 @@ void sim::MergeSimSources::produce(art::Event& e)
       e.put(std::make_unique<aux_det_hits_t>(move(auxdethitCol)), auxdethit_inst);
     }
   }
+
+  if (fFillParticleAncestryMaps) { e.put(std::move(pamCol)); }
 }
 
 void sim::MergeSimSources::dumpConfiguration() const
@@ -426,6 +445,8 @@ void sim::MergeSimSources::dumpConfiguration() const
     log << ")";
   }
 
+  if (fFillParticleAncestryMaps) log << "\n - filling ParticleAncestryMaps";
+
 } // sim::MergeSimSources::dumpConfiguration()
 
 DEFINE_ART_MODULE(sim::MergeSimSources)