Merge pull request #1341 from sophiemiddleton/pion_filters

pions

EventGenerator/fcl/prolog.fcl

@@ -171,14 +171,21 @@ GenFilter : {
         isNull : true
 }
 
+PionFilter : {
+        module_type : PionFilter
+        tmin : 350
+        tmax : 1e8
+        isNull : false
+}
+
 genFilters : # in case we want to add more here
 {
    filters : {
         GenFilter: @local::GenFilter
+        PionFilter: @local::PionFilter
    }
 }
 
-
 EventGenerator : { @table::EventGenerator
     producers : {
         CeEndpointGun               : { @table::CeEndpointGun               }

EventGenerator/src/PionFilter_module.cc

@@ -0,0 +1,84 @@
+// Purpose: Event filter for RPC simulations
+// author: S Middleton  2024
+#include "art/Framework/Core/EDFilter.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "art_root_io/TFileService.h"
+#include "art/Framework/Principal/Handle.h"
+
+// Mu2e includes.
+#include "Offline/MCDataProducts/inc/StageParticle.hh"
+#include "Offline/SeedService/inc/SeedService.hh"
+#include "Offline/GeometryService/inc/GeomHandle.hh"
+#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh"
+#include "Offline/GlobalConstantsService/inc/PhysicsParams.hh"
+#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh"
+#include "Offline/TrackerGeom/inc/Tracker.hh"
+#include <iostream>
+#include <string>
+
+using namespace std;
+namespace mu2e {
+
+  class PionFilter : public art::EDFilter {
+    public:
+      struct Config {
+        using Name=fhicl::Name;
+        using Comment=fhicl::Comment;
+        fhicl::Atom<int> diagLevel{Name("diagLevel"),0};
+        fhicl::Atom<double> tmin{Name("tmin"),0};
+        fhicl::Atom<double> tmax{Name("tmax"),1e6};
+        fhicl::Atom<bool> isNull{Name("isNull"),true};
+      };
+      explicit PionFilter(const art::EDFilter::Table<Config>& config);
+      virtual bool filter(art::Event& event) override;
+
+    private:
+      const SimParticleCollection* SimCol_;
+      int diagLevel_;
+      double tmin_;
+      double tmax_;
+      bool isNull_;
+      double totalweight = 0;
+      double selectedweight = 0;
+  };
+
+  PionFilter::PionFilter(const art::EDFilter::Table<Config>& config) :
+     EDFilter{config}
+    , diagLevel_{config().diagLevel()}
+    , tmin_{config().tmin()}
+    , tmax_{config().tmax()}
+    , isNull_{config().isNull()}
+  {
+  }
+
+  bool PionFilter::filter(art::Event& evt) {
+      if(isNull_) return true;
+      bool passed = false;
+      std::vector<art::Handle<SimParticleCollection>> vah = evt.getMany<SimParticleCollection>();
+      for (auto const& ah : vah) { //always one collection
+        for(const auto& aParticle : *ah){
+          art::Ptr<SimParticle> pp(ah, aParticle.first.asUint());
+          float _endglobaltime = pp->endGlobalTime();
+          if( pp->stoppingCode() == ProcessCode::mu2eKillerVolume and std::abs(pp->pdgId()) == PDGCode::pi_plus){
+            const PhysicsParams& gc = *GlobalConstantsHandle<PhysicsParams>();
+            totalweight += exp(-1*pp->endProperTime() / gc.getParticleLifetime(pp->pdgId()));
+          }
+          if( pp->stoppingCode() == ProcessCode::mu2eKillerVolume and (std::abs(pp->pdgId()) == PDGCode::pi_plus and _endglobaltime > tmin_ and _endglobaltime < tmax_ )){
+            passed = true;
+            const PhysicsParams& gc = *GlobalConstantsHandle<PhysicsParams>();
+            selectedweight += exp(-1*pp->endProperTime() / gc.getParticleLifetime(pp->pdgId()));
+          }
+        }
+      }
+    if(diagLevel_ > 0 ){
+      std::cout<<"Total weight for all stops "<<totalweight<<std::endl;
+      std::cout<<"Selected weight for chosen stops "<<selectedweight<<std::endl;
+    }
+    return passed;
+  }
+}
+
+using mu2e::PionFilter;
+DEFINE_ART_MODULE(PionFilter)

EventGenerator/src/RPCGun_module.cc

@@ -95,6 +95,7 @@ namespace mu2e {
     TH1F* _hPosiPx {nullptr};
     TH1F* _hPosiPy {nullptr};
     TH1F* _hPosiPz {nullptr};
+    TH1F* _hTime {nullptr};
     TH1F* _hMee;
     TH2F* _hMeeVsE;
     TH1F* _hMeeOverE;                   // M(ee)/E(gamma)
@@ -140,6 +141,7 @@ namespace mu2e {
           _hPosiPx  = tfdir.make<TH1F>("hPosiPx" , "Produced positron momentum Px", 140,  -140. , 140.);
           _hPosiPy  = tfdir.make<TH1F>("hPosiPy" , "Produced positron momentum Py", 140,  -140. , 140.);
           _hPosiPz  = tfdir.make<TH1F>("hPosiPz" , "Produced positron momentum Pz", 140,  -140. , 140.);
+          _hTime    = tfdir.make<TH1F>("hTime" , "pion time",100,0,1700);
           _hMee      = tfdir.make<TH1F>("hMee"     , "M(e+e-) "           , 200,0.,200.);
           _hMeeVsE   = tfdir.make<TH2F>("hMeeVsE"  , "M(e+e-) vs E"       , 200,0.,200.,200,0,200);
           _hMeeOverE = tfdir.make<TH1F>("hMeeOverE", "M(e+e-)/E "         , 200, 0.,1);
@@ -152,14 +154,12 @@ namespace mu2e {
       const PhysicsParams& gc = *GlobalConstantsHandle<PhysicsParams>();
       double weight = 0.;
       double tau = part->endProperTime() / gc.getParticleLifetime(part->pdgId());
-
       while(part->parent().isNonnull()) { //while particle has a parent which is not null
         part = part->parent();
-        if ( abs(part->pdgId()) == PDGCode::pi_plus ) { //if pion
+        if ( std::abs(part->pdgId()) == PDGCode::pi_plus ) { //if pion
             tau += part->endProperTime() / gc.getParticleLifetime(part->pdgId());
           }
       }
-
     weight = exp(-tau);
     return weight;
   }
@@ -169,20 +169,17 @@ namespace mu2e {
     auto output{std::make_unique<StageParticleCollection>()};
     const auto simh = event.getValidHandle<SimParticleCollection>(simsToken_);
     const auto pis = stoppedPiMinusList(simh);
-
     if(pis.empty()) {
       throw   cet::exception("BADINPUT")
         <<"RPCGun::produce(): no suitable stopped pion in the input SimParticleCollection\n";
     }
-
     unsigned int randIn = randomFlat_.fireInt(pis.size());
     double time_weight = 1.0;
     if(pionDecayOff_) time_weight = MakeEventWeight(pis[randIn]);
     std::unique_ptr<EventWeight> pw(new EventWeight(time_weight));
     event.put(std::move(pw));
     addParticles(output.get(), pis[randIn]);
     event.put(std::move(output));
-
   }
 
   void RPCGun::addParticles(StageParticleCollection* output,
@@ -230,7 +227,7 @@ namespace mu2e {
           _hPosiPx ->Fill(momp.vect().x());
           _hPosiPy ->Fill(momp.vect().y());
           _hPosiPz ->Fill(momp.vect().z());
-
+          _hTime->Fill(pistop->endGlobalTime());
           double mee = (mome+momp).m();
           _hMee->Fill(mee);
           _hMeeVsE->Fill(energy,mee);