Merge pull request #1300 from brownd1978/surfacestep

Introduce SurfaceStep

CommonMC/fcl/TestSurfaceSteps.fcl

@@ -0,0 +1,29 @@
+#include "Offline/fcl/standardServices.fcl"
+#include "Offline/CommonMC/fcl/prolog.fcl"
+services.DbService.purpose: MDC2020_best
+services.DbService.version: v1_3
+
+process_name: TSS
+source : { module_type : RootInput }
+services : @local::Services.Reco
+physics :
+{
+  producers : {
+    MakeSS : @local::CommonMC.MakeSS
+  }
+  TrigPath : [MakeSS]
+  analyzers : {
+    printSSteps : {
+      module_type : PrintModule
+      surfaceStepPrinter : {
+        verbose : 1
+      }
+    }
+    SSDiag : @local::CommonMC.SSDiag
+  }
+  EndPath : [printSSteps, SSDiag]
+  end_paths : [EndPath]
+  trigger_paths : [TrigPath]
+}
+services.TimeTracker.printSummary: true
+services.TFileService.fileName: "nts.owner.SurfaceStepDiag.version.sequence.root"

CommonMC/fcl/prolog.fcl

@@ -6,26 +6,38 @@ nonProtonGenIds : ["unknown", "particleGun", "cosmicToy", "cosmicDYB", "cosmic",
 muLifeGenIds : [ "StoppedParticleReactionGun", "dioTail", "ExternalRMC", "InternalRMC",
   "MuCapProtonGenTool", "MuCapDeuteronGenTool", "DIOGenTool", "MuCapNeutronGenTool", "MuCapPhotonGenTool", "MuCapGammaRayGenTool" ]
 CommonMC: {
-  producers: {
-    #
-    # time offset producers
-    protonTimeOffset : {
-        module_type: ProtonTimeOffset
-        randPDFparameters : {
-            pulseType : default
-            limitingHalfWidth : 125.0
-            potPulse  : "Offline/Mu2eUtilities/data/potTimingDistribution_20160511.txt"
-            acDipole  : "Offline/Mu2eUtilities/data/acDipoleTransmissionFunction_20160511.txt"
-        }
-    }
-    cosmicTimeOffset : { module_type: CosmicTimeOffset }
-    # Event window marker
+  #
+  # producers
+  #
+  # Event window marker
+  DigiProducers : {
     EWMProducer : {
       module_type : EventWindowMarkerProducer
       RecoFromMCTruth : true
       RecoFromMCTruthErr : 0.5
     }
   }
+
+  protonTimeOffset : {
+    module_type: ProtonTimeOffset
+    randPDFparameters : {
+      pulseType : default
+      limitingHalfWidth : 125.0
+      potPulse  : "Offline/Mu2eUtilities/data/potTimingDistribution_20160511.txt"
+      acDipole  : "Offline/Mu2eUtilities/data/acDipoleTransmissionFunction_20160511.txt"
+    }
+  }
+
+  MakeSS : {
+    module_type : MakeSurfaceSteps
+    debugLevel : 0
+    VDStepPointMCs : "compressDigiMCs:virtualdetector"
+    IPAStepPointMCs : "compressDigiMCs:protonabsorber"
+    TargetStepPointMCs : "compressDigiMCs:stoppingtarget"
+    MaxDistGap : 0.1 # mm
+    MaxTimeGap : 0.1 # ns
+  }
+
   FindMCPrimary : {
     module_type : FindMCPrimary
     debugLevel : 0
@@ -43,6 +55,18 @@ CommonMC: {
     debugLevel : 0
   }
   DigiSim : [ EWMProducer  ]
+  #
+  # Analyzers
+  #
+  SSDiag : {
+    module_type : SurfaceStepDiag
+    SurfaceStepCollection : "MakeSS"
+    VDStepPointMCs : "compressDigiMCs:virtualdetector"
+    IPAStepPointMCs : "compressDigiMCs:protonabsorber"
+    TargetStepPointMCs : "compressDigiMCs:stoppingtarget"
+  }
+
+
 }
 #------------------------------------------------------------------------------
 

CommonMC/src/MakeSurfaceSteps_module.cc

@@ -0,0 +1,158 @@
+//
+//  This module creates SurfaceStep objects from G4 StepPointMCs
+//
+//  Original author: David Brown (LBNL) 7/2024
+//
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Principal/Handle.h"
+#include "cetlib_except/exception.h"
+#include "fhiclcpp/types/Atom.h"
+#include "canvas/Utilities/InputTag.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "Offline/GlobalConstantsService/inc/ParticleData.hh"
+#include "Offline/DataProducts/inc/SurfaceId.hh"
+#include "Offline/MCDataProducts/inc/SurfaceStep.hh"
+#include "Offline/MCDataProducts/inc/StepPointMC.hh"
+#include "Offline/GeometryService/inc/DetectorSystem.hh"
+#include "Offline/GeometryService/inc/GeomHandle.hh"
+
+namespace mu2e {
+
+  class MakeSurfaceSteps : public art::EDProducer {
+    public:
+      using Name=fhicl::Name;
+      using Comment=fhicl::Comment;
+      struct Config {
+        fhicl::Atom<int> debug{ Name("debugLevel"), Comment("Debug Level"), 0};
+        fhicl::Atom<art::InputTag> vdstepmcs { Name("VDStepPointMCs"), Comment("Virtual Detector StepPointMC collection")};
+        fhicl::Atom<art::InputTag> ipastepmcs { Name("IPAStepPointMCs"), Comment("IPAStepPointMC collection")};
+        fhicl::Atom<art::InputTag> ststepmcs { Name("TargetStepPointMCs"), Comment("Stopping target StepPointMC collection")};
+        fhicl::Atom<double> maxdgap{ Name("MaxDistGap"), Comment("Maximum dstance gap between aggregated StepPointMCs")};
+        fhicl::Atom<double> maxtgap{ Name("MaxTimeGap"), Comment("Maximum time gap between aggregated StepPointMCs")};
+      };
+      using Parameters = art::EDProducer::Table<Config>;
+      explicit MakeSurfaceSteps(const Parameters& conf);
+    private:
+      void produce(art::Event& e) override;
+      int debug_;
+      double maxdgap_;
+      double maxtgap_;
+      art::ProductToken<StepPointMCCollection> vdstepmcs_, ipastepmcs_, ststepmcs_;
+      std::map<VirtualDetectorId,SurfaceId> vdmap_; // map of VDIds to surfaceIds
+  };
+
+  MakeSurfaceSteps::MakeSurfaceSteps(const Parameters& config )  :
+    art::EDProducer{config},
+    debug_(config().debug()),
+    maxdgap_(config().maxdgap()),
+    maxtgap_(config().maxtgap()),
+    vdstepmcs_{ consumes<StepPointMCCollection>(config().vdstepmcs())},
+    ipastepmcs_{ consumes<StepPointMCCollection>(config().ipastepmcs())},
+    ststepmcs_{ consumes<StepPointMCCollection>(config().ststepmcs())}
+    {
+      produces <SurfaceStepCollection>();
+      // build the VDId -> SId map by hand. This should come from a service TODO
+      vdmap_[VirtualDetectorId(VirtualDetectorId::TT_FrontHollow)] = SurfaceId("TT_Front");
+      vdmap_[VirtualDetectorId(VirtualDetectorId::TT_Mid)] = SurfaceId("TT_Mid");
+      vdmap_[VirtualDetectorId(VirtualDetectorId::TT_MidInner)] = SurfaceId("TT_Mid");
+      vdmap_[VirtualDetectorId(VirtualDetectorId::TT_Back)] = SurfaceId("TT_Back");
+      vdmap_[VirtualDetectorId(VirtualDetectorId::TT_OutSurf)] = SurfaceId("TT_Outer");
+      vdmap_[VirtualDetectorId(VirtualDetectorId::TT_InSurf)] = SurfaceId("TT_Inner");
+    }
+
+  void MakeSurfaceSteps::produce(art::Event& event) {
+    GeomHandle<DetectorSystem> det;
+    // create output
+    std::unique_ptr<SurfaceStepCollection> ssc(new SurfaceStepCollection);
+    // start with virtual detectors; these are copied directly
+    auto const& vdspmccol_h = event.getValidHandle<StepPointMCCollection>(vdstepmcs_);
+    auto const& vdspmccol = *vdspmccol_h;
+    for(auto const& vdspmc : vdspmccol) {
+      // only some VDs are kept
+      auto isid = vdmap_.find(vdspmc.virtualDetectorId());
+      if(isid != vdmap_.end())ssc->emplace_back(isid->second,vdspmc,det); // no aggregation of VD hits
+    }
+    auto nvdsteps = ssc->size();
+    if(debug_ > 0)std::cout << "Added " << nvdsteps << " VD steps " << std::endl;
+    // now IPA
+    // colate adjacent IPA steppointMCs from the same particle. The following assumes the StepPointMCs from the same SimParticle are adjacent and (roughly) time-ordered. This is observationally true currently, but if G4 ever evolves so that its not, this code will need reworking
+    auto const& ipaspmccol_h = event.getValidHandle<StepPointMCCollection>(ipastepmcs_);
+    auto const& ipaspmccol = *ipaspmccol_h;
+    if(debug_ > 0)std::cout << "IPA StepPointMC collection has " << ipaspmccol.size() << " Entries" << std::endl;
+    // match steps by their sim particle. There is only 1 volume for IPA
+    SurfaceStep ipastep;
+    for(auto const& spmc : ipaspmccol) {
+      bool added(false);
+      // decide if this step is contiguous to existing steps already aggregated
+      // first, test SimParticle (surface is defined by the collection)
+      if(ipastep.simParticle() == spmc.simParticle()){
+        if(debug_ > 2) std::cout << "Same SimParticle" << std::endl;
+        // then time;
+        auto tgap = fabs(ipastep.time()-spmc.time());
+        auto dgap = (ipastep.endPosition() - XYZVectorF(det->toDetector(spmc.position()))).R();
+        if(spmc.time() < ipastep.time())throw cet::exception("Simulation") << " StepPointMC times out-of-order" << std::endl;
+        if(debug_ > 2) std::cout << "Time gap " << tgap << " Distance gap " << dgap << std::endl;
+        if(tgap < maxtgap_ && dgap < maxdgap_){
+          // accumulate this step into the existing surface step
+          if(debug_ > 1)std::cout <<"Added step" << std::endl;
+          ipastep.addStep(spmc,det);
+          added = true;
+        }
+      }
+      if(!added){
+        // if the existing step is valid, save it
+        if(ipastep.surfaceId() != SurfaceIdDetail::unknown && ipastep.simParticle().isNonnull())ssc->push_back(ipastep);
+        // start a new SurfaceStep for this step
+        ipastep = SurfaceStep(SurfaceId(SurfaceIdDetail::IPA),spmc,det);
+        if(debug_ > 1)std::cout <<"New step" << std::endl;
+      }
+    }
+    // save last step
+    if(ipastep.surfaceId() != SurfaceIdDetail::unknown && ipastep.simParticle().isNonnull())ssc->push_back(ipastep);
+    auto nipasteps = ssc->size() - nvdsteps;
+    if(debug_ > 0)std::cout << "Added " << nipasteps << " IPA Steps "<< std::endl;
+    // same for stopping target. Here, the foil number (= volumeid) matters
+    auto const& stspmccol_h = event.getValidHandle<StepPointMCCollection>(ststepmcs_);
+    auto const& stspmccol = *stspmccol_h;
+    if(debug_ > 0)std::cout << "Target StepPointMC collection has " << stspmccol.size() << " Entries" << std::endl;
+    // match steps by their sim particle. There is only 1 volume for IPA
+    SurfaceStep ststep;
+    for(auto const& spmc : stspmccol) {
+      bool added(false);
+      // decide if this step is contiguous to existing steps already aggregated
+      // first, test SimParticle (surface is defined by the collection)
+      if(ststep.simParticle() == spmc.simParticle() && ststep.surfaceId().index()== (int)spmc.volumeId()){
+        if(debug_ > 2) std::cout << "Same SimParticle" << std::endl;
+        // then time;
+        auto tgap = fabs(ststep.time()-spmc.time());
+        auto dgap = (ststep.endPosition() - XYZVectorF(det->toDetector(spmc.position()))).R();
+        if(spmc.time() < ststep.time())throw cet::exception("Simulation") << " StepPointMC times out-of-order" << std::endl;
+        if(debug_ > 2) std::cout << "Time gap " << tgap << " Distance gap " << dgap << std::endl;
+        if(tgap < maxtgap_ && dgap < maxdgap_){
+          // accumulate this step into the existing surface step
+          if(debug_ > 1)std::cout <<"Added step" << std::endl;
+          ststep.addStep(spmc,det);
+          added = true;
+        }
+      }
+      if(!added){
+        // if the existing step is valid, save it
+        if(ststep.surfaceId() != SurfaceIdDetail::unknown && ststep.simParticle().isNonnull())ssc->push_back(ststep);
+        // start a new SurfaceStep for this step
+        // hack around the problem that the ST wires are mixed in with the ST foils
+        SurfaceId stid(SurfaceIdDetail::ST_Foils,spmc.volumeId());
+        if(det->toDetector(spmc.position()).rho() > 75.001)stid = SurfaceId(SurfaceIdDetail::ST_Wires,spmc.volumeId());
+        ststep = SurfaceStep(stid,spmc,det);
+        if(debug_ > 1)std::cout <<"New step" << std::endl;
+      }
+    }
+    if(ststep.surfaceId() != SurfaceIdDetail::unknown && ststep.simParticle().isNonnull())ssc->push_back(ststep);
+    auto nststeps = ssc->size() - nipasteps - nvdsteps;
+    if(debug_ > 0)std::cout << "Added " << nststeps << " stopping target Steps "<< std::endl;
+    // finish
+    event.put(move(ssc));
+  }
+
+}
+DEFINE_ART_MODULE(mu2e::MakeSurfaceSteps)

CommonMC/src/SurfaceStepDiag_module.cc

@@ -0,0 +1,129 @@
+//
+// SurfaceStep diags
+//
+// Original author D. Brown
+//
+// framework
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Core/EDAnalyzer.h"
+#include "Offline/DataProducts/inc/GenVector.hh"
+#include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/OptionalAtom.h"
+#include "art_root_io/TFileService.h"
+#include "TTree.h"
+#include "TH1F.h"
+// data
+#include "Offline/MCDataProducts/inc/SurfaceStep.hh"
+#include "Offline/MCDataProducts/inc/StepPointMC.hh"
+
+namespace mu2e
+{
+  class SurfaceStepDiag : public art::EDAnalyzer {
+    public:
+      using SPM = std::map<art::Ptr<SimParticle>,unsigned>;
+
+      struct Config {
+        using Name = fhicl::Name;
+        using Comment = fhicl::Comment;
+        fhicl::Atom<art::InputTag> SurfaceStepCollection{   Name("SurfaceStepCollection"),   Comment("SurfaceStep collection name") };
+        fhicl::Atom<art::InputTag> vdstepmcs { Name("VDStepPointMCs"), Comment("Virtual Detector StepPointMC collection")};
+        fhicl::Atom<art::InputTag> ipastepmcs { Name("IPAStepPointMCs"), Comment("IPAStepPointMC collection")};
+        fhicl::Atom<art::InputTag> ststepmcs { Name("TargetStepPointMCs"), Comment("Stopping target StepPointMC collection")};
+      };
+
+      explicit SurfaceStepDiag(const art::EDAnalyzer::Table<Config>& config);
+      virtual ~SurfaceStepDiag();
+      virtual void beginJob();
+      virtual void analyze(const art::Event& e);
+    private:
+      art::ProductToken<SurfaceStepCollection> surfsteps_;
+      art::ProductToken<StepPointMCCollection> vdstepmcs_, ipastepmcs_, ststepmcs_;
+      TTree *ssdiag_;
+      int evt_, subrun_, run_, sid_, sindex_;
+      float edep_, path_, time_;
+      XYZVectorF mom_, start_, mid_, end_;
+      TH1F *nvdspmc_, *nipaspmc_, *nstspmc_;
+      TH1F *nvd_, *nipa_, *nstfoil_,* nstwire_;
+  };
+
+  SurfaceStepDiag::SurfaceStepDiag(const art::EDAnalyzer::Table<Config>& config) :
+    art::EDAnalyzer(config),
+    surfsteps_{ consumes<SurfaceStepCollection>(config().SurfaceStepCollection() ) },
+    vdstepmcs_{ consumes<StepPointMCCollection>(config().vdstepmcs())},
+    ipastepmcs_{ consumes<StepPointMCCollection>(config().ipastepmcs())},
+    ststepmcs_{ consumes<StepPointMCCollection>(config().ststepmcs())}
+ {}
+
+  SurfaceStepDiag::~SurfaceStepDiag(){}
+
+  void SurfaceStepDiag::beginJob() {
+    // create diagnostics if requested
+    art::ServiceHandle<art::TFileService> tfs;
+    // detailed diagnostics
+    ssdiag_=tfs->make<TTree>("ssdiag","SurfaceStep diagnostics");
+    ssdiag_->Branch("evt",&evt_,"evt/I");
+    ssdiag_->Branch("subrun",&subrun_,"subrun/I");
+    ssdiag_->Branch("run",&run_,"run/I");
+    ssdiag_->Branch("sid",&sid_,"sid/I");
+    ssdiag_->Branch("sindex",&sindex_,"sindex/I");
+    ssdiag_->Branch("edep",&edep_,"edep/F");
+    ssdiag_->Branch("path",&path_,"path/F");
+    ssdiag_->Branch("time",&time_,"time/F");
+    ssdiag_->Branch("mom.",&mom_);
+    ssdiag_->Branch("start.",&start_);
+    ssdiag_->Branch("mid.",&mid_);
+    ssdiag_->Branch("end.",&end_);
+    nvdspmc_ = tfs->make<TH1F>("nvdspmc","N VD StepPointMCs",100,-0.5,99.5);
+    nipaspmc_ = tfs->make<TH1F>("nipaspmc","N IPA StepPointMCs",100,-0.5,99.5);
+    nstspmc_ = tfs->make<TH1F>("nstspmc","N ST StepPointMCs",100,-0.5,99.5);
+    nvd_ = tfs->make<TH1F>("nvd","N VD SurfaceSteps",100,-0.5,99.5);
+    nipa_ = tfs->make<TH1F>("nipa","N IPA SurfaceSteps",100,-0.5,99.5);
+    nstfoil_ = tfs->make<TH1F>("nstfoil","N ST Foil SurfaceSteps",100,-0.5,99.5);
+    nstwire_ = tfs->make<TH1F>("nstwire","N ST Wire SurfaceSteps",100,-0.5,99.5);
+  }
+
+  void SurfaceStepDiag::analyze(const art::Event& event ) {
+// first histogram StepPointMC precursor sizes
+    auto const& vdspmccol_h = event.getValidHandle<StepPointMCCollection>(vdstepmcs_);
+    auto const& vdspmccol = *vdspmccol_h;
+    nvdspmc_->Fill(vdspmccol.size());
+    auto const& ipaspmccol_h = event.getValidHandle<StepPointMCCollection>(ipastepmcs_);
+    auto const& ipaspmccol = *ipaspmccol_h;
+    nipaspmc_->Fill(ipaspmccol.size());
+    auto const& stspmccol_h = event.getValidHandle<StepPointMCCollection>(ststepmcs_);
+    auto const& stspmccol = *stspmccol_h;
+    nstspmc_->Fill(stspmccol.size());
+
+
+    auto sscolH = event.getValidHandle(surfsteps_);
+    auto const& sscol = *sscolH.product();
+    evt_ = event.id().event();
+    subrun_ = event.id().subRun();
+    run_ = event.id().run();
+    unsigned nvd(0), nipa(0), nstfoil(0), nstwire(0);
+    for(auto const& ss : sscol) {
+      sid_ = ss.surfaceId().id();
+      sindex_ = ss.surfaceId().index();
+      edep_ = ss.energyDeposit();
+      path_ = ss.pathLength();
+      time_ = ss.time();
+      mom_ = ss.momentum();
+      start_ = ss.startPosition();
+      mid_ = ss.midPosition();
+      end_ = ss.endPosition();
+      ssdiag_->Fill();
+      if(ss.surfaceId().id() < SurfaceIdDetail::IPA)nvd++;
+      if(ss.surfaceId().id() == SurfaceIdDetail::IPA)nipa++;
+      if(ss.surfaceId().id() == SurfaceIdDetail::ST_Foils)nstfoil++;
+      if(ss.surfaceId().id() == SurfaceIdDetail::ST_Wires)nstwire++;
+    }
+    nvd_->Fill(nvd);
+    nipa_->Fill(nipa);
+    nstfoil_->Fill(nstfoil);
+    nstwire_->Fill(nstwire);
+  }
+}
+// Part of the magic that makes this class a module.
+using mu2e::SurfaceStepDiag;
+DEFINE_ART_MODULE(SurfaceStepDiag)