LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
NucleonDecay_module.cc
Go to the documentation of this file.
1 // Class: NucleonDecay
3 // Module Type: producer
4 // GENIE nucleon decay generator
5 //
6 // Converted from gNucleonDecayEvGen.cxx by
7 // tjyang@fnal.gov
8 //
9 // 2016 PDG numbering scheme in pp.8-10 of
10 // http://www-pdg.lbl.gov/2016/listings/rpp2016-list-p.pdf (tau1 through tau60)
11 //
13 
19 #include "fhiclcpp/ParameterSet.h"
21 
22 // GENIE includes
23 #include "Framework/Algorithm/AlgFactory.h"
24 #include "Framework/EventGen/EventRecord.h"
25 #include "Framework/EventGen/EventRecordVisitorI.h"
26 #include "Framework/GHEP/GHepParticle.h"
27 #include "Framework/ParticleData/PDGLibrary.h"
28 #include "Framework/Utils/AppInit.h"
29 #include "Physics/NucleonDecay/NucleonDecayMode.h"
30 #include "Physics/NucleonDecay/NucleonDecayUtils.h"
31 
32 // larsoft includes
35 
37 
41 
42 // c++ includes
43 #include <memory>
44 #include <string>
45 
46 #include "CLHEP/Random/RandFlat.h"
47 
48 namespace evgen {
49  class NucleonDecay;
50 }
51 
53 public:
54  explicit NucleonDecay(fhicl::ParameterSet const& p);
55  // The destructor generated by the compiler is fine for classes
56  // without bare pointers or other resource use.
57 
58  // Plugins should not be copied or assigned.
59  NucleonDecay(NucleonDecay const&) = delete;
60  NucleonDecay(NucleonDecay&&) = delete;
61  NucleonDecay& operator=(NucleonDecay const&) = delete;
62  NucleonDecay& operator=(NucleonDecay&&) = delete;
63 
64  // Required functions.
65  void produce(art::Event& e) override;
66 
67  // Selected optional functions.
68  void beginRun(art::Run& run) override;
69 
70 private:
71  // Declare member data here.
72  const genie::EventRecordVisitorI* mcgen;
73  genie::NucleonDecayMode_t gOptDecayMode = genie::kNDNull; // nucleon decay mode
74  int dpdg = 0;
75  CLHEP::RandFlat flatDist;
76 };
77 
79  : art::EDProducer{p}
80  // create a default random engine; obtain the random seed from NuRandomService,
81  // unless overridden in configuration with key "Seed"
83  -> registerAndSeedEngine(createEngine(0), p, "Seed")}
84 {
85  string sname = "genie::EventGenerator";
86  string sconfig = "NucleonDecay";
87  // GENIE v3 needs a tune (even if irrelevant)
88  evgb::SetEventGeneratorListAndTune("Default", "Default");
89 
90  genie::PDGLibrary::Instance(); //Ensure Messenger is started first in GENIE.
91 
92  genie::AlgFactory* algf = genie::AlgFactory::Instance();
93  mcgen = dynamic_cast<const genie::EventRecordVisitorI*>(algf->GetAlgorithm(sname, sconfig));
94  if (!mcgen) {
95  throw cet::exception("NucleonDecay") << "Couldn't instantiate the nucleon decay generator";
96  }
97  int fDecayMode = p.get<int>("DecayMode");
98  gOptDecayMode = (genie::NucleonDecayMode_t)fDecayMode;
99 
100  if (p.get<int>("DecayedNucleon") > 0) { dpdg = p.get<int>("DecayedNucleon"); }
101  else {
102  dpdg = genie::utils::nucleon_decay::DecayedNucleonPdgCode(gOptDecayMode);
103  }
104 
105  produces<std::vector<simb::MCTruth>>();
106  produces<sumdata::RunData, art::InRun>();
107 
108  unsigned int seed = art::ServiceHandle<rndm::NuRandomService>()->getSeed();
109  genie::utils::app_init::RandGen(seed);
110 }
111 
113 {
114  // Implementation of required member function here.
115  genie::EventRecord* event = new genie::EventRecord;
116  int target = 1000180400; //Only use argon target
117  int decay = (int)gOptDecayMode;
118  genie::Interaction* interaction = genie::Interaction::NDecay(target, decay, dpdg);
119  event->AttachSummary(interaction);
120 
121  // Simulate decay
122  mcgen->ProcessEventRecord(event);
123 
124  // genie::Interaction *inter = event->Summary();
125  // const genie::InitialState &initState = inter->InitState();
126  // std::cout<<"initState = "<<initState.AsString()<<std::endl;
127  // const genie::ProcessInfo &procInfo = inter->ProcInfo();
128  // std::cout<<"procInfo = "<<procInfo.AsString()<<std::endl;
129  MF_LOG_DEBUG("NucleonDecay") << "Generated event: " << *event;
130 
131  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
132  simb::MCTruth truth;
133 
135 
136  // Find boundary of active volume
137  double minx = 1e9;
138  double maxx = -1e9;
139  double miny = 1e9;
140  double maxy = -1e9;
141  double minz = 1e9;
142  double maxz = -1e9;
143  for (auto const& tpc : geo->Iterate<geo::TPCGeo>(geo::CryostatID{0})) {
144  if (minx > tpc.MinX()) minx = tpc.MinX();
145  if (maxx < tpc.MaxX()) maxx = tpc.MaxX();
146  if (miny > tpc.MinY()) miny = tpc.MinY();
147  if (maxy < tpc.MaxY()) maxy = tpc.MaxY();
148  if (minz > tpc.MinZ()) minz = tpc.MinZ();
149  if (maxz < tpc.MaxZ()) maxz = tpc.MaxZ();
150  }
151 
152  // Assign vertice position
153  double X0 = flatDist.fire(minx, maxx);
154  double Y0 = flatDist.fire(miny, maxy);
155  double Z0 = flatDist.fire(minz, maxz);
156 
157  TIter partitr(event);
158  genie::GHepParticle* part = 0;
159  // GHepParticles return units of GeV/c for p. the V_i are all in fermis
160  // and are relative to the center of the struck nucleus.
161  // add the vertex X/Y/Z to the V_i for status codes 0 and 1
162  int trackid = 0;
163  std::string primary("primary");
164 
165  while ((part = dynamic_cast<genie::GHepParticle*>(partitr.Next()))) {
166 
167  simb::MCParticle tpart(
168  trackid, part->Pdg(), primary, part->FirstMother(), part->Mass(), part->Status());
169 
170  TLorentzVector pos(X0, Y0, Z0, 0);
171  TLorentzVector mom(part->Px(), part->Py(), part->Pz(), part->E());
172  tpart.AddTrajectoryPoint(pos, mom);
173  if (part->PolzIsSet()) {
174  TVector3 polz;
175  part->GetPolarization(polz);
176  tpart.SetPolarization(polz);
177  }
178  tpart.SetRescatter(part->RescatterCode());
179  truth.Add(tpart);
180 
181  ++trackid;
182  } // end loop to convert GHepParticles to MCParticles
183  truth.SetOrigin(simb::kUnknown);
184  truthcol->push_back(truth);
185  //FillHistograms(truth);
186  e.put(std::move(truthcol));
187 
188  delete event;
189 }
190 
192 {
194  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
195 }
196 
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
base_engine_t & createEngine(seed_t seed)
genie::NucleonDecayMode_t gOptDecayMode
constexpr auto fullRun()
NucleonDecay(fhicl::ParameterSet const &p)
CLHEP::RandFlat flatDist
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
const genie::EventRecordVisitorI * mcgen
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Geometry information for a single TPC.
Definition: TPCGeo.h:36
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
Functions for transforming GENIE objects into ART objects (and back)
Particle class.
void beginRun(art::Run &run) override
Definition: Run.h:37
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
TString part[npart]
Definition: Style.C:32
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
long seed
Definition: chem4.cc:67
void produce(art::Event &e) override
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
NucleonDecay & operator=(NucleonDecay const &)=delete
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
void SetEventGeneratorListAndTune(const std::string &evtlistname="", const std::string &tunename="${GENIE_XSEC_TUNE}")
Definition: GENIE2ART.cxx:128
#define MF_LOG_DEBUG(id)
cout<< "-> Edep in the target
Definition: analysis.C:53
Definition: MVAAlg.h:12
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:203
Event generator information.
Definition: MCTruth.h:32
Float_t e
Definition: plot.C:35
Namespace collecting geometry-related classes utilities.
Event Generation using GENIE, cosmics or single particles.
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Event finding and building.
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192