GENIENeutronOscGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief GENIE neutron oscillation event generator.
3 /// Multiple functions are taken from GENIE source code
4 /// src/Apps/gNNBarOscEvGen.cxx, and credits go to the
5 /// author there.
6 /// Please ignore the warnings of 'hides overloaded virtual
7 /// function [-Woverloaded-virtual]' when compiling, as it is
8 /// harmless has been fixed in more recent releases of GENIE.
9 /// It will disappear when NOvA progresses to use more recent
10 /// versions of GENIE in the future.
11 /// \author junting@utexas.edu
12 /// \date 2018/01/26
13 ////////////////////////////////////////////////////////////////////////
14 
15 #include <cassert>
16 #include <cstdlib>
17 #include <string>
18 #include <sstream>
19 #include <vector>
20 #include <map>
21 #include <unistd.h>
22 
23 // Both TStopwatch.h and EventRecord.h have these warnings. They're not our
24 // fault, so hide them for now.
25 #pragma GCC diagnostic push
26 #pragma GCC diagnostic ignored "-Woverloaded-virtual"
27 
28 // ROOT includes
29 #include "TStopwatch.h"
30 
31 // GENIE includes
32 #include "GENIE/Framework/EventGen/EventRecord.h"
33 #pragma GCC diagnostic pop
34 
35 #include "GENIE/Framework/EventGen/EventRecordVisitorI.h"
36 #include "GENIE/Framework/Algorithm/AlgFactory.h"
37 #include "GENIE/Framework/Numerical/RandomGen.h"
38 #include "GENIE/Framework/GHEP/GHepParticle.h"
39 #include "GENIE/Physics/NNBarOscillation/NNBarOscUtils.h"
40 
41 // Framework includes
46 #include "fhiclcpp/ParameterSet.h"
50 
51 // NOvA includes
52 #include "Geometry/Geometry.h"
56 #include "SummaryData/POTSum.h"
57 #include "SummaryData/SpillData.h"
58 #include "SummaryData/RunData.h"
59 #include "SummaryData/SubRunData.h"
60 #include "Utilities/AssociationUtil.h"
61 
62 
63 namespace evgen {
64 
66 
67  public:
68 
69  explicit GENIENeutronOscGen(fhicl::ParameterSet const &pset);
70  virtual ~GENIENeutronOscGen();
71 
72  void produce(art::Event& evt);
73  void beginJob();
74  void beginRun(art::Run &run);
75  void endSubRun(art::SubRun &sr);
77  int selectAnnihilationMode(int pdgCode);
78  void fillMCTruth(const genie::EventRecord* record, simb::MCTruth &truth);
80 
81  private:
82 
83  TStopwatch stopwatch_; ///< keep track of how long it takes to run the job
84  int cycle_; ///< cycle number in the MC generation
85  int target_; ///< pdg code of the target
86  genie::NNBarOscMode_t gOptDecayMode_; ///< neutron oscillation mode
88  double vertexX_; // m
89  double vertexY_; // m
90  double vertexZ_; // m
91  double vertexT_; // s
94  double fdWidth_;
95  double fdHeight_;
96  double fdLength_;
97  };
98 }
99 
100 namespace evgen {
101 
102  //___________________________________________________________________________
104  : cycle_(pset.get<int>("Cycle", 0))
105  , target_(pset.get<int>("Target", 1000060120))
106  , randomVertexPosition_(pset.get<bool>("RandomVertexPosition", false))
107  , vertexX_(pset.get<double>("VertexX", 0.))
108  , vertexY_(pset.get<double>("VertexY", 0.))
109  , vertexZ_(pset.get<double>("VertexZ", 30.))
110  , vertexT_(pset.get<double>("VertexT", 225.e-6))
111  {
112  stopwatch_.Start();
113 
114  gOptDecayMode_ = (genie::NNBarOscMode_t) pset.get<int>("AnnihilationMode", 0);
116  if (!validMode) {
117  mf::LogError("GENIENeutronOscGen") << "You need to specify a valid annihilation mode (0-16)";
118  exit(0);
119  }
120 
122  fdWidth_ = 2. * geo->DetHalfWidth() / 100.; // m
123  fdHeight_ = 2. * geo->DetHalfHeight() / 100.; // m
124  fdLength_ = geo->DetLength() / 100.; // m
125 
127 
129  rnd_->SetSeed(0);
130 
131  produces<std::vector<simb::MCTruth>>();
132  produces<sumdata::SubRunData, art::InSubRun>();
133  produces<sumdata::RunData, art::InRun>();
134  }
135 
136  //___________________________________________________________________________
138  {
139  // delete mcgen_;
140 
141  stopwatch_.Stop();
142  mf::LogInfo("GENIENeutronOscGen") << "real time to produce file: "
143  << stopwatch_.RealTime();
144  }
145 
146  //___________________________________________________________________________
148  {
149  }
150 
151  //___________________________________________________________________________
153  {
155 
156  std::unique_ptr<sumdata::RunData>
157  runcol(new sumdata::RunData(geo->DetId(),
158  geo->FileBaseName(),
159  geo->ExtractGDML()));
160 
161  run.put(std::move(runcol));
162 
163  return;
164  }
165 
166  //___________________________________________________________________________
168  {
170  std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(cycle_));
171  sr.put(std::move(sd));
172  }
173 
174  //___________________________________________________________________________
176  {
177  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
178  simb::MCTruth truth;
179 
181  genie::Interaction* interaction = genie::Interaction::NOsc(target_, decay);
183  event->AttachSummary(interaction);
184  mcgen_->ProcessEventRecord(event);
185  // mf::LogInfo("GENIENeutronOscGen") << "Generated event: " << *event;
186 
187  event->SetVertex(vertexX_, vertexY_, vertexZ_, vertexT_); // in si unit, m and s
188  if (randomVertexPosition_) {
190  }
191 
192  fillMCTruth(event, truth);
193  truthcol->push_back(truth);
194 
195  evt.put(std::move(truthcol));
196  }
197 
198  //___________________________________________________________________________
200  {
201  double randomNumber = rnd_->RndNum().Rndm();
202  double randomX = (randomNumber - 0.5) * fdWidth_; // m
203 
204  randomNumber = rnd_->RndNum().Rndm();
205  double randomY = (randomNumber - 0.5) * fdHeight_; // m
206 
207  randomNumber = rnd_->RndNum().Rndm();
208  double randomZ = randomNumber * fdLength_; // m
209 
210  record->SetVertex(randomX, randomY, randomZ, record->Vertex()->T());
211  }
212 
213  //___________________________________________________________________________
215  {
216  string sname = "genie::EventGenerator";
217  string sconfig = "NeutronOsc";
219  mcgen_ = dynamic_cast<const genie::EventRecordVisitorI *> (algf->GetAlgorithm(sname,sconfig));
220 
221  if (!mcgen_) {
222  mf::LogError("GENIENeutronOscGen") << "Couldn't instantiate the neutron oscillation generator";
223  // genie::gAbortingInErr = true;
224  exit(1);
225  }
226  }
227 
228  //___________________________________________________________________________
230  {
232  int mode;
233 
234  std::string pdgString = std::to_string(static_cast<long long>(pdgCode));
235  if (pdgString.size() != 10) {
236  mf::LogError("GENIENeutronOscGen") << "Expecting PDG code to be a 10-digit integer; instead, it's the following: " << pdgString;
237  // genie::gAbortingInErr = true;
238  exit(1);
239  }
240 
241  int nNucleons = std::stoi(pdgString.substr(6,3)) - 1;
242  int nProtons = std::stoi(pdgString.substr(3,3));
243 
244  double protonFrac = ((double) nProtons) / ((double) nNucleons);
245  double neutronFrac = 1 - protonFrac;
246 
247  const int nModes = 16;
248  double br [nModes] = { 0.010, 0.080, 0.100, 0.220,
249  0.360, 0.160, 0.070, 0.020,
250  0.015, 0.065, 0.110, 0.280,
251  0.070, 0.240, 0.100, 0.100 };
252 
253  for (int i = 0; i < nModes; i++) {
254  if (i < 7)
255  br[i] *= protonFrac;
256  else
257  br[i] *= neutronFrac;
258  }
259 
260  double randomNumber = rnd_->RndNum().Rndm();
261  double threshold = 0.;
262  for (int i = 0; i < nModes; i++) {
263  threshold += br[i];
264  if (randomNumber < threshold) {
265  mode = i + 1;
266  return mode;
267  }
268  }
269 
270  mf::LogError("GENIENeutronOscGen") << "Random selection of final state failed!";
271 
272  // genie::gAbortingInErr = true;
273  exit(1);
274  }
275 
276  else {
277  int mode = (int) gOptDecayMode_;
278  return mode;
279  }
280  }
281 
282  //___________________________________________________________________________
284  {
285  TIter partitr(record);
287  int trackid = 0;
288  std::string primary("primary");
289 
290  while( (part = dynamic_cast<genie::GHepParticle*>(partitr.Next())) ) {
291  if (part->Status() != genie::kIStStableFinalState ) {
292  continue;
293  }
294 
295  simb::MCParticle tpart(trackid,
296  part->Pdg(),
297  primary,
298  part->FirstMother(),
299  part->Mass(),
300  part->Status());
301 
302  TLorentzVector* vertex = record->Vertex();
303  double vtx[4] = {
304  100. * (part->Vx() * 1.e-15 + vertex->X()), // cm
305  100. * (part->Vy() * 1.e-15 + vertex->Y()), // cm
306  100. * (part->Vz() * 1.e-15 + vertex->Z()), // cm
307  part->Vt() + vertex->T() * 1.e9 // ns
308  };
309 
310  TLorentzVector position(vtx[0], vtx[1], vtx[2], vtx[3]);
311  TLorentzVector momentum(part->Px(), part->Py(), part->Pz(), part->E());
312  tpart.AddTrajectoryPoint(position, momentum);
313  truth.Add(tpart);
314  trackid++;
315  }
316  }
317 
318 }
319 
TStopwatch stopwatch_
keep track of how long it takes to run the job
GENIENeutronOscGen(fhicl::ParameterSet const &pset)
Definition: event.h:34
void setRandomEventVertexPosition(genie::EventRecord *record)
virtual void ProcessEventRecord(GHepRecord *event_rec) const =0
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
int cycle_
cycle number in the MC generation
double DetLength() const
double Mass(void) const
Mass that corresponds to the PDG code.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
virtual void AttachSummary(Interaction *interaction)
Definition: GHepRecord.cxx:143
object containing MC flux information
Definition: Run.h:31
virtual TLorentzVector * Vertex(void) const
Definition: GHepRecord.h:141
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
double Vt(void) const
Get production time.
Definition: GHepParticle.h:98
int Pdg(void) const
Definition: GHepParticle.h:64
int FirstMother(void) const
Definition: GHepParticle.h:67
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double sd(Eigen::VectorXd x)
Summary information for an interaction.
Definition: Interaction.h:56
TString part[npart]
Definition: Style.C:32
enum genie::ENNBarOscMode NNBarOscMode_t
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
double DetHalfHeight() const
Definition: run.py:1
genie::NNBarOscMode_t gOptDecayMode_
neutron oscillation mode
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators & other numerical methods
Definition: RandomGen.h:78
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double DetHalfWidth() const
Definition: decay.py:1
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
double Vz(void) const
Get production z.
Definition: GHepParticle.h:97
void fillMCTruth(const genie::EventRecord *record, simb::MCTruth &truth)
ProductID put(std::unique_ptr< PROD > &&)
virtual void SetVertex(double x, double y, double z, double t)
Definition: GHepRecord.cxx:863
exit(0)
bool IsValidMode(NNBarOscMode_t ndm)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
double Vy(void) const
Get production y.
Definition: GHepParticle.h:96
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
Event generator information.
Definition: MCTruth.h:32
const genie::EventRecordVisitorI * mcgen_
Float_t e
Definition: plot.C:35
Helper for AttenCurve.
Definition: Path.h:10
Module to generate only pions from cosmic rays.
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
Encapsulate the geometry of one entire detector (near, far, ndos)
static constexpr Double_t sr
Definition: Munits.h:164
static Interaction * NOsc(int tgt, int annihilation_mode=-1)
int target_
pdg code of the target
std::string FileBaseName() const
double Vx(void) const
Get production x.
Definition: GHepParticle.h:95
void SetSeed(long int seed)
Definition: RandomGen.cxx:90
double Py(void) const
Get Py.
Definition: GHepParticle.h:90