Implement digi mixing for tracker

CommonMC/src/ProtonBunchTimeMCFromProtonBunchTime_module.cc

@@ -0,0 +1,52 @@
+// Convert proton bunch time marker from data to MC data product
+// To be used for mixing simulated signal into real data from detector
+//
+// Ed Callaghan, 2024
+
+// art
+#include <art/Framework/Core/EDProducer.h>
+#include <art/Framework/Principal/Event.h>
+
+// canvas
+#include <canvas/Utilities/InputTag.h>
+
+// fhiclcpp
+#include <fhiclcpp/types/Atom.h>
+
+// mu2e
+#include <Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh>
+#include <Offline/RecoDataProducts/inc/ProtonBunchTime.hh>
+
+namespace mu2e{
+  class ProtonBunchTimeMCFromProtonBunchTime: public art::EDProducer{
+    public:
+      struct Config{
+        fhicl::Atom<art::InputTag> pbttag{
+          fhicl::Name("PBTTag"),
+          fhicl::Comment("ProtonBunchTime provenance"),
+          "" // default, eventually to specify "from data"
+        };
+      };
+
+      using Parameters = art::EDProducer::Table<Config>;
+      explicit ProtonBunchTimeMCFromProtonBunchTime(Parameters const& config);
+    protected:
+      /**/
+    private:
+      void produce(art::Event& event) override;
+      art::InputTag pbttag_;
+  };
+
+  ProtonBunchTimeMCFromProtonBunchTime::ProtonBunchTimeMCFromProtonBunchTime(Parameters const& config): EDProducer(config), pbttag_(config().pbttag()){
+    this->produces<ProtonBunchTimeMC>();
+  }
+
+  void ProtonBunchTimeMCFromProtonBunchTime::produce(art::Event& event){
+    auto pbt = event.getHandle<ProtonBunchTime>(pbttag_);
+    std::unique_ptr<ProtonBunchTimeMC> pbtmc(new ProtonBunchTimeMC);
+    pbtmc->pbtime_ = pbt->pbtime_;
+    event.put(std::move(pbtmc));
+  }
+}
+
+DEFINE_ART_MODULE(mu2e::ProtonBunchTimeMCFromProtonBunchTime);

EventMixing/inc/Mu2eProductMixer.hh

@@ -41,7 +41,9 @@
 #include "Offline/MCDataProducts/inc/CosmicLivetime.hh"
 #include "Offline/MCDataProducts/inc/SimTimeOffset.hh"
 #include "Offline/MCDataProducts/inc/PhysicalVolumeInfoMultiCollection.hh"
-
+#include "Offline/RecoDataProducts/inc/StrawDigi.hh"
+#include "Offline/MCDataProducts/inc/StrawDigiMC.hh"
+#include "Offline/DataProducts/inc/EventWindowMarker.hh"
 
 
 //================================================================
@@ -102,6 +104,10 @@ namespace mu2e {
       fhicl::Table<CollectionMixerConfig> crvStepMixer { fhicl::Name("crvStepMixer") };
       fhicl::Table<CollectionMixerConfig> extMonSimHitMixer { fhicl::Name("extMonSimHitMixer") };
       fhicl::Table<CollectionMixerConfig> eventIDMixer { fhicl::Name("eventIDMixer") };
+      fhicl::Table<CollectionMixerConfig> strawDigiMixer { fhicl::Name("strawDigiMixer") };
+      fhicl::Table<CollectionMixerConfig> strawDigiADCWaveformMixer { fhicl::Name("strawDigiADCWaveformMixer") };
+      fhicl::Table<CollectionMixerConfig> strawDigiMCMixer { fhicl::Name("strawDigiMCMixer") };
+      fhicl::Table<CollectionMixerConfig> eventWindowMarkerMixer { fhicl::Name("eventWindowMarkerMixer") };
       fhicl::OptionalTable<CosmicLivetimeMixerConfig> cosmicLivetimeMixer { fhicl::Name("cosmicLivetimeMixer") };
       fhicl::OptionalTable<VolumeInfoMixerConfig> volumeInfoMixer { fhicl::Name("volumeInfoMixer") };
       fhicl::OptionalAtom<art::InputTag> simTimeOffset { fhicl::Name("simTimeOffset"), fhicl::Comment("Simulation time offset to apply (optional)") };
@@ -148,6 +154,23 @@ namespace mu2e {
                           ExtMonFNALSimHitCollection& out,
                           art::PtrRemapper const& remap);
 
+    bool mixStrawDigis(std::vector<StrawDigiCollection const*> const& in,
+                       StrawDigiCollection& out,
+                       art::PtrRemapper const& remap);
+
+    bool mixStrawDigiADCWaveforms(std::vector<StrawDigiADCWaveformCollection const*> const& in,
+                       StrawDigiADCWaveformCollection& out,
+                       art::PtrRemapper const& remap);
+
+    bool mixStrawDigiMCs(std::vector<StrawDigiMCCollection const*> const& in,
+                       StrawDigiMCCollection& out,
+                       art::PtrRemapper const& remap);
+
+    bool mixEventWindowMarkers(std::vector<EventWindowMarker const*> const& in,
+                       EventWindowMarker& out,
+                       art::PtrRemapper const& remap);
+
+
     bool mixEventIDs(std::vector<art::EventIDSequence const*> const &in,
                      art::EventIDSequence& out,
                      art::PtrRemapper const& remap);

EventMixing/src/MixDigis_module.cc

@@ -0,0 +1,113 @@
+// Mix digi-level data products, as opposed to simulation-level data products,
+// into art events
+//
+// Ed Callaghan, 2023
+
+// stl
+
+// art
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/SubRun.h"
+#include "art/Framework/IO/ProductMix/MixHelper.h"
+#include "art/Framework/Modules/MixFilter.h"
+#include "art_root_io/RootIOPolicy.h"
+#include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/Comment.h"
+#include "fhiclcpp/types/Name.h"
+#include "fhiclcpp/types/Sequence.h"
+#include "fhiclcpp/types/Table.h"
+#include "fhiclcpp/types/TupleAs.h"
+
+// cetlib_except
+#include "cetlib_except/exception.h"
+
+// mu2e
+#include "Offline/EventMixing/inc/Mu2eProductMixer.hh"
+
+namespace mu2e{
+  // detail class for art::MixFilter<>
+  class MixDigisDetail{
+    public:
+      struct Mu2eConfig{
+        // configuration is simple --- which data products to mix in
+        fhicl::Table<Mu2eProductMixer::Config> products{
+          fhicl::Name("products"),
+          fhicl::Comment("Products to be mixed, as form of a mixingMap, i.e. tuples of InputTags to output instance names.")
+
+        };
+      };
+
+      struct Config{
+        fhicl::Table<Mu2eConfig> mu2e{ fhicl::Name("mu2e") };
+      };
+
+      using Parameters = art::MixFilterTable<Config>;
+      explicit MixDigisDetail(const Parameters& parameters,
+                                  art::MixHelper& helper);
+      size_t nSecondaries();
+
+      void processEventIDs(const art::EventIDSequence& seq);
+      void beginSubRun(const art::SubRun& subrun);
+      void endSubRun(art::SubRun& subrun);
+      void startEvent(const art::Event& event);
+      void finalizeEvent(art::Event& event);
+
+    protected:
+      /**/
+
+    private:
+      Mu2eProductMixer _mixer;
+      art::EventIDSequence evids;
+  };
+
+  // implementation
+  MixDigisDetail::MixDigisDetail(const Parameters& parameters,
+                                         art::MixHelper& helper)
+      : _mixer{parameters().mu2e().products(), helper}{
+    helper.produces<art::EventIDSequence>();
+  }
+
+  // always overlay simulated event onto one event window
+  size_t MixDigisDetail::nSecondaries(){
+    size_t rv = 1;
+    return rv;
+  }
+
+  // always keep track of secondary EventIDs
+  void MixDigisDetail::processEventIDs(const art::EventIDSequence& seq){
+    // intent is to overlay signals onto a nominally-complete event,
+    // so we should only mix in one event
+    if (seq.size() != 1){
+      throw cet::exception("MIX") << "mu2e::DixDigisDetail: mixing more than one digi frame" << std::endl;
+    }
+    _mixer.processEventIDs(seq);
+    evids = seq;
+  }
+
+  // forward beginSubRun
+  void MixDigisDetail::beginSubRun(const art::SubRun& subrun){
+    _mixer.beginSubRun(subrun);
+  }
+
+  // forward endSubRun
+  void MixDigisDetail::endSubRun(art::SubRun& subrun){
+    _mixer.endSubRun(subrun);
+  }
+
+  // forward startEvent
+  void MixDigisDetail::startEvent(const art::Event& event){
+    _mixer.startEvent(event);
+  }
+
+  void MixDigisDetail::finalizeEvent(art::Event& event){
+    // "manually" mix in EventIDs, as these do not preexist as data products
+    auto seq = std::make_unique<art::EventIDSequence>();
+    seq->swap(evids);
+    event.put(std::move(seq));
+  }
+
+  // define the module
+  typedef art::MixFilter<MixDigisDetail,art::RootIOPolicy> MixDigis;
+}
+
+DEFINE_ART_MODULE(mu2e::MixDigis);

EventMixing/src/Mu2eProductMixer.cc

@@ -6,6 +6,7 @@
 #include <algorithm>
 #include <iterator>
 #include <iostream>
+#include <string>
 
 #include "cetlib_except/exception.h"
 
@@ -96,6 +97,26 @@ namespace mu2e {
         (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixEventIDs, *this);
     }
 
+    for(const auto& e: conf.strawDigiMixer().mixingMap()) {
+      helper.declareMixOp
+        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixStrawDigis, *this);
+    }
+
+    for(const auto& e: conf.strawDigiADCWaveformMixer().mixingMap()) {
+      helper.declareMixOp
+        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixStrawDigiADCWaveforms, *this);
+    }
+
+    for(const auto& e: conf.strawDigiMCMixer().mixingMap()) {
+      helper.declareMixOp
+        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixStrawDigiMCs, *this);
+    }
+
+    for(const auto& e: conf.eventWindowMarkerMixer().mixingMap()) {
+      helper.declareMixOp
+        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixEventWindowMarkers, *this);
+    }
+
     //----------------------------------------------------------------
     // VolumeInfo handling
 
@@ -397,6 +418,43 @@ namespace mu2e {
     return true;
   }
 
+  bool Mu2eProductMixer::mixStrawDigis(std::vector<StrawDigiCollection const*> const& in,
+                     StrawDigiCollection& out,
+                     art::PtrRemapper const& remap)
+  {
+    art::flattenCollections(in, out);
+    return true;
+  }
+
+  bool Mu2eProductMixer::mixStrawDigiADCWaveforms(std::vector<StrawDigiADCWaveformCollection const*> const& in,
+                     StrawDigiADCWaveformCollection& out,
+                     art::PtrRemapper const& remap)
+  {
+    art::flattenCollections(in, out);
+    return true;
+  }
+
+  bool Mu2eProductMixer::mixStrawDigiMCs(std::vector<StrawDigiMCCollection const*> const& in,
+                     StrawDigiMCCollection& out,
+                     art::PtrRemapper const& remap)
+  {
+    art::flattenCollections(in, out);
+    return true;
+  }
+
+  bool Mu2eProductMixer::mixEventWindowMarkers(std::vector<EventWindowMarker const*> const& in,
+                     EventWindowMarker& out,
+                     art::PtrRemapper const& remap){
+    // assert that only one EventWindowMarker be present
+    if (in.size() != 1){
+      std::string msg = "attempting to mix more than 1 EventWindowMarker: ";
+      msg += std::to_string(in.size()) + " present";
+      throw cet::exception("BADINPUT") << msg << std::endl;
+    }
+    out = *in[0];
+    return true;
+  }
+
   //----------------------------------------------------------------
   bool Mu2eProductMixer::mixEventIDs(std::vector<art::EventIDSequence const*> const &in,
                                      art::EventIDSequence& out,

MCDataProducts/inc/StrawDigiMC.hh

@@ -1,7 +1,7 @@
 #ifndef MCDataProducts_StrawDigiMC_hh
 #define MCDataProducts_StrawDigiMC_hh
 //
-//  Summary of MC information used to create a StrawDigi.  Everything is referenced by the thresold digitization end
+//  Summary of MC information used to create a StrawDigi.  Everything is referenced by the threshold digitization end
 //
 // Original author David Brown, LBNL
 //
@@ -20,28 +20,43 @@
 #include <iostream>
 #include <vector>
 #include <array>
+#include <string>
 
 // Mu2e includes
+#include <Offline/GeneralUtilities/inc/EnumToStringSparse.hh>
 
 namespace mu2e {
+    // enum equipped with std::string descriptions
+    class StrawDigiProvenanceDetail{
+      public:
+        enum enum_type {unknown=0, Simulation, Mixed, External};
+        static std::string const& typeName();
+        static std::map<enum_type, std::string> const& names();
+    };
+    using StrawDigiProvenance = EnumToStringSparse<StrawDigiProvenanceDetail>;
 
   class StrawDigiMC{
 
   public:
+
     typedef art::Ptr<StrawGasStep> SGSP;
     typedef std::array<SGSP,StrawEnd::nends> SGSPA;
     typedef std::array<XYZVectorF,StrawEnd::nends> PA;
     typedef std::array<float,StrawEnd::nends> FA;
 
     StrawDigiMC();
     // construct from hitlets
-    StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs);
+    StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs, StrawDigiProvenance::enum_type=StrawDigiProvenance::Simulation);
 
     // use compuater copy construcors
     StrawDigiMC(const StrawDigiMC& rhs, SGSPA sgsp ); // update the Ptrs
+    StrawDigiMC(const StrawDigiMC& rhs, StrawDigiProvenance::enum_type provenance ); // update validity
+
     // Accessors
     StrawId strawId() const { return _strawid; }
 
+    StrawDigiProvenance provenance() const { return _provenance; }
+
     SGSP const&  strawGasStep(StrawEnd strawend) const { return _sgspa[strawend]; }
     SGSPA const&  strawGasSteps() const { return _sgspa; }
     StrawEnd earlyEnd() const { return (_wtime[StrawEnd::cal] < _wtime[StrawEnd::hv]) ?  StrawEnd::cal : StrawEnd::hv; }
@@ -68,6 +83,7 @@ namespace mu2e {
     FA _ctime; // times of the trigger clusters
     FA _wtime; // times at the wire ends of the signals which fired the TDC.
     SGSPA _sgspa; // StrawGasStep that triggered each end
+    StrawDigiProvenance _provenance; // level of association with any true MC events
   };
 
   inline std::ostream& operator<<( std::ostream& ost,

MCDataProducts/src/StrawDigiMC.cc

@@ -15,19 +15,41 @@ using namespace std;
 
 namespace mu2e {
 
+  std::string const& StrawDigiProvenanceDetail::typeName(){
+    static const std::string rv = "Provenance";
+    return rv;
+  }
+
+  static const std::map<StrawDigiProvenanceDetail::enum_type, std::string> nam{
+    std::make_pair(StrawDigiProvenance::unknown,    "unknown"),
+    std::make_pair(StrawDigiProvenance::Simulation, "Simulation"),
+    std::make_pair(StrawDigiProvenance::Mixed,      "Mixed"),
+    std::make_pair(StrawDigiProvenance::External,   "External")
+  };
+
+  std::map<StrawDigiProvenance::enum_type, std::string> const& StrawDigiProvenanceDetail::names(){
+    return nam;
+  }
+
   // Default constructor is required for persistable classes
   StrawDigiMC::StrawDigiMC()
     : _strawid(StrawId::_invalid)
+    , _provenance(StrawDigiProvenance::Simulation)
   {}
 
-  StrawDigiMC::StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs):
+  StrawDigiMC::StrawDigiMC(StrawId sid, PA cpos, FA ctime, FA wetime, SGSPA sgs, StrawDigiProvenance::enum_type provenance):
     _strawid(sid), _cpos(cpos), _ctime(ctime), _wtime(wetime), _sgspa(sgs)
+    , _provenance(provenance)
   {}
 
   StrawDigiMC::StrawDigiMC(const StrawDigiMC& rhs, SGSPA sgspa ) : StrawDigiMC(rhs)  {
     _sgspa = sgspa;
   }
 
+  StrawDigiMC::StrawDigiMC(const StrawDigiMC& rhs, StrawDigiProvenance::enum_type provenance): StrawDigiMC(rhs){
+    _provenance = provenance;
+  }
+
   bool StrawDigiMC::isCrossTalk(StrawEnd strawend) const {
     bool retval(false);
     if(!_sgspa[strawend].isNull()){