From c88cc148b62dece100b6ee3a8701039305c55e0a Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 21 Jun 2021 16:06:30 -0500
Subject: [PATCH 1/5] Fix IPA material for stop filtering

---
 Analyses/src/StoppedParticlesDumper_module.cc | 37 +++++++++++++------
 .../geom/protonAbsorber_cylindrical_v04.txt   |  2 +-
 Mu2eG4/src/ConstructMaterials.cc              |  2 +-
 3 files changed, 28 insertions(+), 13 deletions(-)

diff --git a/Analyses/src/StoppedParticlesDumper_module.cc b/Analyses/src/StoppedParticlesDumper_module.cc
index 236c7686b6..5f5f9f8a2c 100644
--- a/Analyses/src/StoppedParticlesDumper_module.cc
+++ b/Analyses/src/StoppedParticlesDumper_module.cc
@@ -50,21 +50,26 @@ namespace mu2e {
       float z;
       float t;
       float tau; // proper time, for stopped pion weights
+      int pdg;
+      int origin;
+      int term;
 
-      StopInfo() : x(), y(), z(), t(), tau() {}
+      StopInfo() : x(), y(), z(), t(), tau(), pdg(0), origin(0), term(0) {}
 
       StopInfo(const art::Ptr<SimParticle>& p, const VspMC& spMCcolls, float tt)
         : x(p->endPosition().x())
         , y(p->endPosition().y())
         , z(p->endPosition().z())
         , t(p->endGlobalTime())
-        , tau(tt)
-      {
-        if(!p->endDefined()) {
-          throw cet::exception("BADINPUTS")
-            <<"StoppedParticlesDumper: input SimParticle does not have end defined!\n";
-        }
-      }
+	, tau(tt)
+	, pdg(p->pdgId())
+	, origin(p->creationCode())
+	, term(p->stoppingCode()) {
+	  if(!p->endDefined()) {
+	    throw cet::exception("BADINPUTS")
+	      <<"StoppedParticlesDumper: input SimParticle does not have end defined!\n";
+	  }
+	}
     };
 
   }// namespace
@@ -101,6 +106,12 @@ namespace mu2e {
           false
           };
 
+      fhicl::Atom<bool> writeCodes {
+        Name("writeParticleCodes"),
+          Comment("Write out particle codes"),
+          false
+          };
+
       fhicl::Sequence<int> decayOffPDGCodes {
         Name("decayOffPDGCodes"),
           Comment("A list of PDG IDs of particles that had their decay process turned off during\n"
@@ -125,7 +136,7 @@ namespace mu2e {
   private:
     bool dumpSimParticleLeaves_;
     art::InputTag input_;
-    bool writeProperTime_;
+    bool writeProperTime_, writeCodes_;
     std::vector<art::InputTag> hitColls_;
 
     std::vector<int> decayOffCodes_;
@@ -143,12 +154,12 @@ namespace mu2e {
     dumpSimParticleLeaves_(conf().dumpSimParticleLeaves()),
     input_(conf().inputCollection()),
     writeProperTime_(conf().writeProperTime()),
+    writeCodes_(conf().writeCodes()),
     nt_()
   {
     if(writeProperTime_) {
       hitColls_ =  conf().hitCollections();
       decayOffCodes_ = conf().decayOffPDGCodes();
-
       // must sort to use binary_search in SimParticleGetTau
       std::sort(decayOffCodes_.begin(), decayOffCodes_.end());
     }
@@ -158,9 +169,13 @@ namespace mu2e {
   void StoppedParticlesDumper::beginJob() {
     art::ServiceHandle<art::TFileService> tfs;
     std::string branchDesc("x/F:y/F:z/F:time/F");
-    if(writeProperTime_) {
+    if(writeProperTime_ || writeCodes_) {
       branchDesc += ":tauNormalized/F";
     }
+    if(writeCodes_){
+      branchDesc += ":pdg/I:origin/I:term/I";
+    }
+      
     nt_ = tfs->make<TTree>( "stops", "Stopped particles ntuple");
     nt_->Branch("stops", &data_, branchDesc.c_str());
   }
diff --git a/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt b/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt
index 5ad48c638a..4db85733ae 100644
--- a/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt
+++ b/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt
@@ -18,7 +18,7 @@ int protonabsorber.verbosityLevel = 0;
 
 //IPA updates
 
-string protonabsorber.materialName = "Polyethylene090";
+string protonabsorber.materialName = "IPAPolyethylene";
 double protonabsorber.OutRadius1   =  300.5; //design is 300 mm inner radius, so add thickness to outer here
 double protonabsorber.OutRadius0   =  300.5;
 double protonabsorber.thickness    =    0.5;
diff --git a/Mu2eG4/src/ConstructMaterials.cc b/Mu2eG4/src/ConstructMaterials.cc
index f0d579c65d..513a55b37f 100644
--- a/Mu2eG4/src/ConstructMaterials.cc
+++ b/Mu2eG4/src/ConstructMaterials.cc
@@ -258,7 +258,7 @@ namespace mu2e {
 
     mat = uniqueMaterialOrThrow( "IPAPolyethylene");
     {
-      G4Material* IPAPolyethylene = new G4Material( mat.name, 0.96*CLHEP::g/CLHEP::cm3, 2);
+      G4Material* IPAPolyethylene = new G4Material( mat.name, 0.90*CLHEP::g/CLHEP::cm3, 2);
       IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_H"), 0.14);
       IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_C"), 0.86);
     }

From a7b93dcdf680b7ea0af866a660654de968269707 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Mon, 21 Jun 2021 18:32:52 -0500
Subject: [PATCH 2/5] Add CosmicLivetime and clean out unused parts.  Allow
 wildcarding the process codes when selecting particles

---
 EventMixing/inc/Mu2eProductMixer.hh      | 13 +++------
 EventMixing/src/Mu2eProductMixer.cc      | 34 +++++-------------------
 Filters/src/ParticleCodeFilter_module.cc | 10 ++++---
 Mu2eG4/src/ConstructMaterials.cc         |  6 ++---
 4 files changed, 20 insertions(+), 43 deletions(-)

diff --git a/EventMixing/inc/Mu2eProductMixer.hh b/EventMixing/inc/Mu2eProductMixer.hh
index 850c4cb9c9..1213d2d36b 100644
--- a/EventMixing/inc/Mu2eProductMixer.hh
+++ b/EventMixing/inc/Mu2eProductMixer.hh
@@ -37,7 +37,7 @@
 #include "MCDataProducts/inc/StrawGasStep.hh"
 #include "MCDataProducts/inc/CrvStep.hh"
 #include "MCDataProducts/inc/ExtMonFNALSimHitCollection.hh"
-#include "MCDataProducts/inc/ProtonBunchIntensity.hh"
+#include "MCDataProducts/inc/CosmicLivetime.hh"
 #include "MCDataProducts/inc/SimParticleTimeMap.hh"
 #include "MCDataProducts/inc/SimTimeOffset.hh"
 #include "MCDataProducts/inc/PhysicalVolumeInfoMultiCollection.hh"
@@ -93,8 +93,7 @@ namespace mu2e {
       fhicl::Table<CollectionMixerConfig> strawGasStepMixer { fhicl::Name("strawGasStepMixer") };
       fhicl::Table<CollectionMixerConfig> crvStepMixer { fhicl::Name("crvStepMixer") };
       fhicl::Table<CollectionMixerConfig> extMonSimHitMixer { fhicl::Name("extMonSimHitMixer") };
-      fhicl::Table<CollectionMixerConfig> protonBunchIntensityMixer { fhicl::Name("protonBunchIntensityMixer") };
-      fhicl::Table<CollectionMixerConfig> protonTimeMapMixer { fhicl::Name("protonTimeMapMixer") };
+      fhicl::Table<CollectionMixerConfig> cosmicLivetimeMixer { fhicl::Name("cosmicLivetimeMixer") };
       fhicl::Table<CollectionMixerConfig> eventIDMixer { fhicl::Name("eventIDMixer") };
       fhicl::OptionalTable<VolumeInfoMixerConfig> volumeInfoMixer { fhicl::Name("volumeInfoMixer") };
       fhicl::Atom<art::InputTag> simTimeOffset { fhicl::Name("simTimeOffset"), fhicl::Comment("Simulation time offset to apply (optional)"), art::InputTag() };
@@ -140,14 +139,10 @@ namespace mu2e {
                           ExtMonFNALSimHitCollection& out,
                           art::PtrRemapper const& remap);
 
-    bool mixProtonBunchIntensity(std::vector<mu2e::ProtonBunchIntensity const*> const &in,
-                                 mu2e::ProtonBunchIntensity& out,
+    bool mixCosmicLivetime(std::vector<mu2e::CosmicLivetime const*> const &in,
+                                 mu2e::CosmicLivetime& out,
                                  art::PtrRemapper const& remap);
 
-    bool mixProtonTimeMap(std::vector<mu2e::SimParticleTimeMap const*> const &in,
-                                mu2e::SimParticleTimeMap& out,
-                                art::PtrRemapper const& remap);
-
     bool mixEventIDs(std::vector<art::EventIDSequence const*> const &in,
                      art::EventIDSequence& out,
                      art::PtrRemapper const& remap);
diff --git a/EventMixing/src/Mu2eProductMixer.cc b/EventMixing/src/Mu2eProductMixer.cc
index 21b5c318d8..f626b4bdb1 100644
--- a/EventMixing/src/Mu2eProductMixer.cc
+++ b/EventMixing/src/Mu2eProductMixer.cc
@@ -88,14 +88,9 @@ namespace mu2e {
         (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixExtMonSimHits, *this);
     }
 
-    for(const auto& e: conf.protonBunchIntensityMixer().mixingMap()) {
+    for(const auto& e: conf.cosmicLivetimeMixer().mixingMap()) {
       helper.declareMixOp
-        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixProtonBunchIntensity, *this);
-    }
-
-    for(const auto& e: conf.protonTimeMapMixer().mixingMap()) {
-      helper.declareMixOp
-        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixProtonTimeMap, *this);
+        (e.inTag, e.resolvedInstanceName(), &Mu2eProductMixer::mixCosmicLivetime, *this);
     }
 
     for(const auto& e: conf.eventIDMixer().mixingMap()) {
@@ -348,29 +343,14 @@ namespace mu2e {
   }
 
   //----------------------------------------------------------------
-  bool Mu2eProductMixer::mixProtonBunchIntensity(std::vector<ProtonBunchIntensity const*> const& in,
-                                                 ProtonBunchIntensity& out,
+  bool Mu2eProductMixer::mixCosmicLivetime(std::vector<CosmicLivetime const*> const& in,
+                                                 CosmicLivetime& out,
                                                  art::PtrRemapper const& remap)
   {
+    if(in.size() > 1)
+      throw cet::exception("BADINPUT")<<"Mu2eProductMixer/evt: can't mix CosmicLiveTime" << std::endl; 
     for(const auto& x: in) {
-      out.add(*x);
-    }
-
-    return true;
-  }
-
-  bool Mu2eProductMixer::mixProtonTimeMap(std::vector<SimParticleTimeMap const*> const& in,
-                                          SimParticleTimeMap& out,
-                                          art::PtrRemapper const& remap)
-  {
-    for(size_t incount = 0; incount < in.size(); ++incount) {
-      auto const& timemap = *in[incount];
-      //std::cout << "Mixing time map " << incount << " size " << timemap.size() << std::endl;
-      for(auto & imap : timemap) {
-        auto newptr = remap(imap.first, simOffsets_[incount]);
-        out[newptr] = imap.second;
-        // do I need to go down the chain?  I think not
-      }
+      out = *x;
     }
 
     return true;
diff --git a/Filters/src/ParticleCodeFilter_module.cc b/Filters/src/ParticleCodeFilter_module.cc
index 0b105cada9..35cb784f16 100644
--- a/Filters/src/ParticleCodeFilter_module.cc
+++ b/Filters/src/ParticleCodeFilter_module.cc
@@ -32,9 +32,11 @@ namespace mu2e {
 
       struct ParticleCodeSelector {
 	PDGCode::type pdgCode_;
-	ProcessCode creationCode_, terminationCode_;
+	int creationCode_, terminationCode_;
 	bool select(SimParticle const& part) const {
-	  return part.pdgId() == pdgCode_ && part.creationCode() == creationCode_ && part.stoppingCode() == terminationCode_;
+	  return part.pdgId() == pdgCode_
+	    && ( creationCode_ < 0 || part.creationCode() == creationCode_) 
+	    && ( terminationCode_ < 0 || part.stoppingCode() == terminationCode_);
 	}
       };
 
@@ -54,8 +56,8 @@ namespace mu2e {
     for(auto const& pconfig : conf().codeConfig()) {
       ParticleCodeSelector pselector;
       pselector.pdgCode_ = static_cast<PDGCode::type>(std::get<0>(pconfig));
-      pselector.creationCode_ = ProcessCode(static_cast<ProcessCode::enum_type>(std::get<1>(pconfig)));
-      pselector.terminationCode_ = ProcessCode(static_cast<ProcessCode::enum_type>(std::get<2>(pconfig)));
+      pselector.creationCode_ = std::get<1>(pconfig);
+      pselector.terminationCode_ = std::get<2>(pconfig);
       if(printLevel_ > 0) std::cout << "Creating selector of PDGcode " << pselector.pdgCode_ 
       << " creation code " << pselector.creationCode_
       << " termination code " << pselector.terminationCode_ << std::endl;
diff --git a/Mu2eG4/src/ConstructMaterials.cc b/Mu2eG4/src/ConstructMaterials.cc
index 513a55b37f..02bb58360e 100644
--- a/Mu2eG4/src/ConstructMaterials.cc
+++ b/Mu2eG4/src/ConstructMaterials.cc
@@ -258,9 +258,9 @@ namespace mu2e {
 
     mat = uniqueMaterialOrThrow( "IPAPolyethylene");
     {
-      G4Material* IPAPolyethylene = new G4Material( mat.name, 0.90*CLHEP::g/CLHEP::cm3, 2);
-      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_H"), 0.14);
-      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_C"), 0.86);
+      G4Material* IPAPolyethylene = new G4Material( mat.name, 0.96*CLHEP::g/CLHEP::cm3, 2);
+      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_H"), 0.12);
+      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_C"), 0.88);
     }
 
     mat = uniqueMaterialOrThrow( "Polyethylene096");

From b97b7dc4813f6b786486b00aa7241949bf566a2f Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Tue, 22 Jun 2021 11:15:45 -0500
Subject: [PATCH 3/5] use typed members in filter config.  Provide references
 for IPA values, update to latest measurement

---
 Filters/src/ParticleCodeFilter_module.cc       | 12 ++++++------
 Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt |  2 +-
 Mu2eG4/src/ConstructMaterials.cc               |  6 +++---
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/Filters/src/ParticleCodeFilter_module.cc b/Filters/src/ParticleCodeFilter_module.cc
index 35cb784f16..bacf83414a 100644
--- a/Filters/src/ParticleCodeFilter_module.cc
+++ b/Filters/src/ParticleCodeFilter_module.cc
@@ -27,16 +27,16 @@ namespace mu2e {
 	using Comment=fhicl::Comment;
 	fhicl::Atom<int> printLevel { Name("PrintLevel"), Comment ("Printout Level"), 0 };
 	fhicl::Atom<art::InputTag> simParticles { Name("SimParticles"), Comment("SimParticle collection") };
-	ParticleCodeConfig codeConfig { fhicl::Name("ParticleCodes") };
+	ParticleCodeConfig codeConfig { fhicl::Name("ParticleCodes"), Comment("particle code to select: PDG, and optionally creation or termination code") };
       };
 
       struct ParticleCodeSelector {
 	PDGCode::type pdgCode_;
-	int creationCode_, terminationCode_;
+	ProcessCode::enum_type creationCode_, terminationCode_;
 	bool select(SimParticle const& part) const {
 	  return part.pdgId() == pdgCode_
-	    && ( creationCode_ < 0 || part.creationCode() == creationCode_) 
-	    && ( terminationCode_ < 0 || part.stoppingCode() == terminationCode_);
+	    && ( creationCode_ == ProcessCode::unknown || part.creationCode() == creationCode_) 
+	    && ( terminationCode_ == ProcessCode::unknown || part.stoppingCode() == terminationCode_);
 	}
       };
 
@@ -56,8 +56,8 @@ namespace mu2e {
     for(auto const& pconfig : conf().codeConfig()) {
       ParticleCodeSelector pselector;
       pselector.pdgCode_ = static_cast<PDGCode::type>(std::get<0>(pconfig));
-      pselector.creationCode_ = std::get<1>(pconfig);
-      pselector.terminationCode_ = std::get<2>(pconfig);
+      pselector.creationCode_ = static_cast<ProcessCode::enum_type>(std::get<1>(pconfig));
+      pselector.terminationCode_ = static_cast<ProcessCode::enum_type>(std::get<2>(pconfig));
       if(printLevel_ > 0) std::cout << "Creating selector of PDGcode " << pselector.pdgCode_ 
       << " creation code " << pselector.creationCode_
       << " termination code " << pselector.terminationCode_ << std::endl;
diff --git a/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt b/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt
index 4db85733ae..f4b91724c7 100644
--- a/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt
+++ b/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt
@@ -21,7 +21,7 @@ int protonabsorber.verbosityLevel = 0;
 string protonabsorber.materialName = "IPAPolyethylene";
 double protonabsorber.OutRadius1   =  300.5; //design is 300 mm inner radius, so add thickness to outer here
 double protonabsorber.OutRadius0   =  300.5;
-double protonabsorber.thickness    =    0.5;
+double protonabsorber.thickness    =    0.511; // measured by S. Krave 6/22/2021
 int    protonabsorber.verbosityLevel = 0;
 bool   protonabsorber.visible      = true;
 bool   protonabsorber.solid        = false;
diff --git a/Mu2eG4/src/ConstructMaterials.cc b/Mu2eG4/src/ConstructMaterials.cc
index 02bb58360e..7d8e64cbdb 100644
--- a/Mu2eG4/src/ConstructMaterials.cc
+++ b/Mu2eG4/src/ConstructMaterials.cc
@@ -258,9 +258,9 @@ namespace mu2e {
 
     mat = uniqueMaterialOrThrow( "IPAPolyethylene");
     {
-      G4Material* IPAPolyethylene = new G4Material( mat.name, 0.96*CLHEP::g/CLHEP::cm3, 2);
-      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_H"), 0.12);
-      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_C"), 0.88);
+      G4Material* IPAPolyethylene = new G4Material( mat.name, 0.954*CLHEP::g/CLHEP::cm3, 2);
+      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_H"), 0.11);
+      IPAPolyethylene->AddMaterial( findMaterialOrThrow("G4_C"), 0.89); // Carbon doped Polytehylene, additional carbon 2-5% from MDS (DeWal DW 402B),  density measured by S. Krave 6/22/2021
     }
 
     mat = uniqueMaterialOrThrow( "Polyethylene096");

From 228d547603153bc8cb64db01b2f97001dd5074f4 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Tue, 22 Jun 2021 14:25:50 -0500
Subject: [PATCH 4/5] adjust comment

---
 Filters/src/ParticleCodeFilter_module.cc | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Filters/src/ParticleCodeFilter_module.cc b/Filters/src/ParticleCodeFilter_module.cc
index bacf83414a..6a06a4ef5b 100644
--- a/Filters/src/ParticleCodeFilter_module.cc
+++ b/Filters/src/ParticleCodeFilter_module.cc
@@ -27,7 +27,8 @@ namespace mu2e {
 	using Comment=fhicl::Comment;
 	fhicl::Atom<int> printLevel { Name("PrintLevel"), Comment ("Printout Level"), 0 };
 	fhicl::Atom<art::InputTag> simParticles { Name("SimParticles"), Comment("SimParticle collection") };
-	ParticleCodeConfig codeConfig { fhicl::Name("ParticleCodes"), Comment("particle code to select: PDG, and optionally creation or termination code") };
+	ParticleCodeConfig codeConfig { Name("ParticleCodes"), Comment("particle code to select: PDG, creation, and termination code \n"
+	"(select 'unknown' to disable testing creation and/or termination codes)") };
       };
 
       struct ParticleCodeSelector {

From e7882297d0d23534dbcbe428245a991ef49fad74 Mon Sep 17 00:00:00 2001
From: David Nathan Brown <dave_brown@lbl.gov>
Date: Tue, 22 Jun 2021 16:28:20 -0500
Subject: [PATCH 5/5] Change to string for processcode

---
 Filters/src/ParticleCodeFilter_module.cc | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/Filters/src/ParticleCodeFilter_module.cc b/Filters/src/ParticleCodeFilter_module.cc
index 6a06a4ef5b..7ece532610 100644
--- a/Filters/src/ParticleCodeFilter_module.cc
+++ b/Filters/src/ParticleCodeFilter_module.cc
@@ -16,11 +16,12 @@
 #include "fhiclcpp/types/Sequence.h"
 #include "DataProducts/inc/PDGCode.hh"
 #include "MCDataProducts/inc/ProcessCode.hh"
+#include <string>
 
 namespace mu2e {
   class ParticleCodeFilter : public art::EDFilter {
     public: 
-      using ParticleCodeConfig = fhicl::Sequence<fhicl::Tuple<int,int,int>>;
+      using ParticleCodeConfig = fhicl::Sequence<fhicl::Tuple<int,std::string,std::string>>;
 
       struct ModuleConfig {
 	using Name=fhicl::Name;
@@ -57,8 +58,8 @@ namespace mu2e {
     for(auto const& pconfig : conf().codeConfig()) {
       ParticleCodeSelector pselector;
       pselector.pdgCode_ = static_cast<PDGCode::type>(std::get<0>(pconfig));
-      pselector.creationCode_ = static_cast<ProcessCode::enum_type>(std::get<1>(pconfig));
-      pselector.terminationCode_ = static_cast<ProcessCode::enum_type>(std::get<2>(pconfig));
+      pselector.creationCode_ = ProcessCode::findByName(std::get<1>(pconfig));
+      pselector.terminationCode_ = ProcessCode::findByName(std::get<2>(pconfig));
       if(printLevel_ > 0) std::cout << "Creating selector of PDGcode " << pselector.pdgCode_ 
       << " creation code " << pselector.creationCode_
       << " termination code " << pselector.terminationCode_ << std::endl;
