GENIEGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief GENIE neutrino event generator
3 /// \author brebel@fnal.gov, rhatcher@fnal.gov
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <cassert>
8 #include <cstdlib>
9 #include <string>
10 #include <regex>
11 #include <sstream>
12 #include <vector>
13 #include <map>
14 #include <unistd.h>
15 
16 // ROOT includes
17 #include "TStopwatch.h"
18 
19 // Framework includes
24 #include "fhiclcpp/ParameterSet.h"
28 
29 
30 // NOvA includes
31 #include "nugen/EventGeneratorBase/GENIE/GENIEHelper.h"
32 #include "Geometry/Geometry.h"
36 #include "SummaryData/POTSum.h"
37 #include "SummaryData/SpillData.h"
38 #include "SummaryData/RunData.h"
39 #include "SummaryData/SubRunData.h"
40 #include "Utilities/AssociationUtil.h"
41 
42 #include "dk2nu/tree/dk2nu.h"
43 #include "dk2nu/tree/NuChoice.h"
44 #include "dk2nu/genie/GDk2NuFlux.h"
45 #include "dk2nu/tree/dkmeta.h"
46 
47 ///Monte Carlo event generation
48 namespace evgen {
49 
50  /// A module to check the results from the Monte Carlo generator
51  class GENIEGen : public art::EDProducer {
52 
53  public:
54 
55  explicit GENIEGen(fhicl::ParameterSet const &pset);
56  virtual ~GENIEGen();
57 
58  void produce(art::Event& evt);
59  void beginJob();
60  void beginRun(art::Run &run);
61  void endSubRun(art::SubRun &sr);
62 
63  private:
64 
65  evgb::GENIEHelper *fGENIEHelp; ///< GENIEHelper object
66  int fPassEmptySpills; ///< whether or not to allow for events with no interactions to be put into the event
67  TStopwatch fStopwatch; ///< keep track of how long it takes to run the job
68  int fSpillCounter; ///< how many spills so far?
69  double fPOTPerSpill; ///< how many POT to simulate per spill
70  double fEventsPerSpill; ///< how many events to simulate per spill
71  double fTotalExposure; ///< total exposure in POT
72  double fTotalPOTLimit; ///< total POT limit
73  int fCycle; ///< cycle number in the MC generation
74  };
75 }
76 
77 namespace evgen {
78 
79  //___________________________________________________________________________
81  : fGENIEHelp (0)
82  , fPassEmptySpills (pset.get< bool >("PassEmptySpills"))
83  , fSpillCounter (0)
84  , fPOTPerSpill (pset.get< double >("POTPerSpill", 5.0e13))
85  , fEventsPerSpill (pset.get< double >("EventsPerSpill", 0))
86  , fTotalExposure (0)
87  , fTotalPOTLimit (pset.get< double >("TotalPOTLimit"))
88  , fCycle (pset.get< int >("Cycle", 0))
89  {
90  fStopwatch.Start();
91 
92  produces< std::vector<simb::MCTruth> >();
93  produces< std::vector<simb::MCFlux> >();
94  produces< std::vector<simb::GTruth> >();
95  produces< sumdata::SpillData >();
96  produces< sumdata::POTSum, art::InSubRun >();
97  produces< sumdata::SubRunData, art::InSubRun >();
98  produces< sumdata::RunData, art::InRun >();
99  // Associate every truth with the flux it came from
100  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
101  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
102  // Dk2Nu additions
103  produces< std::vector<bsim::Dk2Nu> >();
104  produces< std::vector<bsim::NuChoice> >();
105 #ifdef PUT_DK2NU_ASSN
106  produces< art::Assns<simb::MCTruth, bsim::Dk2Nu> >();
107  produces< art::Assns<simb::MCTruth, bsim::NuChoice> >();
108 #endif
109 
111  fGENIEHelp = new evgb::GENIEHelper(pset,
112  geo->ROOTGeoManager(),
113  geo->ROOTFile(),
114  geo->TotalMass(pset.get< std::string>("TopVolume").c_str()));
115  }
116 
117  //___________________________________________________________________________
119  {
120  fStopwatch.Stop();
121  mf::LogInfo("GENIEGen") << "real time to produce file: "
122  << fStopwatch.RealTime();
123  delete fGENIEHelp; // clean up, and let dtor do its thing
124  }
125 
126  //___________________________________________________________________________
128  {
129  }
130 
131  //___________________________________________________________________________
133  {
134  // grab the geometry object to see what geometry we are using
136 
137  std::unique_ptr<sumdata::RunData>
138  runcol(new sumdata::RunData(geo->DetId(),
139  geo->FileBaseName(),
140  geo->ExtractGDML()));
141 
142  run.put(std::move(runcol));
143 
144  // initialize the GENIEHelper here rather than in beginJob to
145  // avoid problems with the Geometry reloading at a run boundary.
146  // If we ever make more than one run in a single job we will have
147  // to re-evaluate
148  fGENIEHelp->Initialize();
149  fTotalExposure = 0.0;
150 
151  return;
152  }
153 
154  //___________________________________________________________________________
156  {
157  std::unique_ptr< sumdata::POTSum > p(new sumdata::POTSum);
158 
159  // p->totpot = fGENIEHelp->TotalExposure();
160  // p->totgoodpot = fGENIEHelp->TotalExposure();
161  p->totpot = fTotalExposure;
162  p->totgoodpot = fTotalExposure;
163  p->totspills = fSpillCounter;
164  p->goodspills = fSpillCounter;
165  p->Print(std::cout);
166 
167  sr.put(std::move(p));
168 
169  // store the cycle information
171  std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
172  sr.put(std::move(sd));
173  }
174 
175  //___________________________________________________________________________
177  {
178  // A temporary value is needed to store the spill exposure that GENIEHelper uses. GENIEGen
179  // needs to remember the number after GENIEHelper has reset it to zero for the purposes of
180  // updating fTotalExposure.
181  double SpillExpTemp = 0.0;
182 
183  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
184  std::unique_ptr< std::vector<simb::MCFlux> > fluxcol (new std::vector<simb::MCFlux >);
185  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol (new std::vector<simb::GTruth >);
186  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > tgtassn(new art::Assns<simb::MCTruth, simb::GTruth>);
187  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > assns(new art::Assns<simb::MCTruth, simb::MCFlux>);
188 
189  std::unique_ptr< std::vector<bsim::Dk2Nu> >
190  dk2nucol(new std::vector<bsim::Dk2Nu>);
191  std::unique_ptr< std::vector<bsim::NuChoice> >
192  nuchoicecol(new std::vector<bsim::NuChoice>);
193 
194  std::unique_ptr< art::Assns<simb::MCTruth, bsim::Dk2Nu> >
196  std::unique_ptr< art::Assns<simb::MCTruth, bsim::NuChoice> >
198 
199  bool passedPOTLimit = false;
200 
201  // keep throwing spills until we get one with something in it
202  while ( truthcol->size() < 1 ) {
203 
204  // POT limit. Maybe we passed the limit in a previous event, so now we
205  // have to keep skipping until the end of the file. Or maybe we're in the
206  // process of breaking out of the inner loop, and need to break again.
207  if(fTotalPOTLimit > 0 &&
208  fTotalExposure + fGENIEHelp->SpillExposure() >= fTotalPOTLimit){
209  passedPOTLimit = true;
210  break;
211  }
212 
213  while ( ! fGENIEHelp->Stop() ) {
214 
215  simb::MCTruth truth;
217  simb::GTruth gTruth;
218 
219  // GENIEHelper returns a false in the sample method if
220  // either no neutrino was generated, or the interaction
221  // occurred beyond the detector's z extent - ie something we
222  // would never see anyway.
223  if ( fGENIEHelp->Sample(truth, flux, gTruth ) ) {
224 
225  // POT limit
226  if(fTotalPOTLimit > 0 &&
227  fTotalExposure + fGENIEHelp->SpillExposure() > fTotalPOTLimit){
228  // Top the total exposure up to the right value
229  SpillExpTemp = fTotalPOTLimit - fTotalExposure;
230  break;
231  }
232 
233  // When running in "POT per Spill" mode, this will prevent the last neutrino (that was produced after
234  // the POT limit had been passed) from being put into the event.
235  //
236  // This is only a temporary fix to correct the problem of GENIEHelper always making one more neutrino
237  // than it should.
238  if((fGENIEHelp->SpillExposure() <= fPOTPerSpill && fEventsPerSpill == 0) || fEventsPerSpill > 0) {
239  truthcol ->push_back(truth);
240  gtruthcol->push_back(gTruth);
241  fluxcol ->push_back(flux);
242 
243  util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns,
244  fluxcol->size()-1, fluxcol->size());
245 
246  util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *tgtassn,
247  gtruthcol->size()-1, gtruthcol->size());
248 
249  genie::GFluxI* fdriver = fGENIEHelp->GetFluxDriver(true);
250  genie::flux::GDk2NuFlux* dk2nuDriver =
251  dynamic_cast<genie::flux::GDk2NuFlux*>(fdriver);
252  if ( dk2nuDriver ) {
253  const bsim::Dk2Nu& dk2nuObj = dk2nuDriver->GetDk2Nu();
254  dk2nucol ->push_back(dk2nuObj);
255  const bsim::NuChoice& nuchoiceObj = dk2nuDriver->GetNuChoice();
256  nuchoicecol->push_back(nuchoiceObj);
257 
258 #ifdef PUT_DK2NU_ASSN
259  evgb::util::CreateAssn(*this, evt, *truthcol, *dk2nucol, *dk2nuassn,
260  dk2nucol->size()-1, dk2nucol->size());
261  evgb::util::CreateAssn(*this, evt, *truthcol, *nuchoicecol, *nuchoiceassn,
262  nuchoicecol->size()-1, nuchoicecol->size());
263 #endif
264  }
265 
266  // update exposure
267  SpillExpTemp = fGENIEHelp->SpillExposure();
268  }
269 
270  } // end if genie was able to make an event
271 
272  } // end event generation loop
273 
274  // check to see if we are to pass empty spills
275  if ( truthcol->size() < 1 && fPassEmptySpills ) {
276  LOG_DEBUG("GENIEGen")
277  << "no events made for this spill "
278  << " continue on until we get a spill with events";
279  break;
280  }
281 
282  } // end loop while no interactions are made
283 
284  // put the collections in the event
285  evt.put(std::move(truthcol));
286  evt.put(std::move(fluxcol));
287  evt.put(std::move(gtruthcol));
288  evt.put(std::move(assns));
289  evt.put(std::move(tgtassn));
290 
291  // in the constructor we said these were produced ...
292  // ... so we have to put them in the record, even if just empty
293  evt.put(std::move(dk2nucol));
294  evt.put(std::move(nuchoicecol));
295 #ifdef PUT_DK2NU_ASSN
296  evt.put(std::move(dk2nuassn));
297  evt.put(std::move(nuchoiceassn));
298 #endif
299 
300  ++fSpillCounter;
301 
302  // If fEventsPerSpill > 0, then we should update by what we actually used.
303  // Otherwise, we are in "POT per Spill" mode so we should update by the cutoff
304  // value, fPOTPerSpill.
305 
306  double ThisSpillPoT = fPOTPerSpill;
307  if(fEventsPerSpill > 0 || passedPOTLimit) ThisSpillPoT=SpillExpTemp;
308  fTotalExposure += ThisSpillPoT;
309 
310  std::unique_ptr<sumdata::SpillData> sd(new sumdata::SpillData);
311  sd->spillpot = ThisSpillPoT;
312  sd->goodbeam = 1;
313 
314  // Fill in some typical values to hopefully keep good beam cuts happy
315  // without having to special-case for MC.
316  sd->deltaspilltimensec = 0;
317  sd->hornI = -199.5;
318  sd->posx = 0.9;
319  sd->posy = 1.3;
320  sd->widthx = 1.2;
321  sd->widthy = 1.2;
322 
323  genie::GFluxI* fdriver = fGENIEHelp->GetFluxDriver(true);
324  genie::flux::GDk2NuFlux* dk2nuDriver = dynamic_cast<genie::flux::GDk2NuFlux*>(fdriver);
325  // if we're using dk2nu flux, set SpillData from the flux metadata
326  if ( dk2nuDriver ) {
327  // get metadata for last generated dk2nu ray
328  // all rays used for a single simulation run have the same horn current and beam position
329  const bsim::DkMeta &dkmeta = dk2nuDriver->GetDkMeta();
330  //strip out alphabetic characters and convert to a float
331  //only works if the format is something like 200i or -200i
332  std::string stripped_horncfg = std::regex_replace(dkmeta.horncfg, std::regex(R"([A-Za-z])"), "");
333 
334  sd->hornI = std::stof(stripped_horncfg);
335 
336  std::cout<<"sd->hornI: "<<sd->hornI<<std::endl;
337 
338  if (sd->hornI < 0){
339  sd->isRHC = true;
340  }else if (sd->hornI == 0){
341  sd->is0HC = true;
342  }else {
343  sd->isRHC = false;
344  sd->is0HC = false;
345  }
346 
347  std::cout<<"sd->isRHC : "<<sd->isRHC<<std::endl;
348  std::cout<<"sd->is0HC : "<<sd->is0HC<<std::endl;
349 
350  sd->posx = dkmeta.beam0x;
351  sd->posy = dkmeta.beam0y;
352  sd->widthx = dkmeta.beamhwidth*10; //convert to mm
353  sd->widthy = dkmeta.beamvwidth*10;
354  }
355 
356  //based on data
357 /*
358  sd->bposx[0]=-0.8;
359  sd->bposx[1]=-0.2;
360  sd->bposx[2]=-0.04;
361  sd->bposx[3]=0.05;
362  sd->bposx[4]=-0.008;
363  sd->bposx[5]=-0.007;
364 
365  for(int i=0;i<6;i++){
366  sd->bposy[i]=0.5;
367  //assuming no slip-stacking
368  sd->intx[i]=ThisSpillPoT/6;
369  sd->inty[i]=ThisSpillPoT/6;
370  }
371  */
372  evt.put(std::move(sd));
373  }
374 
375 }// namespace
376 
377 namespace evgen { DEFINE_ART_MODULE(GENIEGen) }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
double fEventsPerSpill
how many events to simulate per spill
int fSpillCounter
how many spills so far?
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void beginRun(art::Run &run)
Double_t beamvwidth
vertical width of beam
Definition: dkmeta.h:110
const char * p
Definition: xmltok.h:285
Double_t beamhwidth
horizontal width of beam
Definition: dkmeta.h:109
void endSubRun(art::SubRun &sr)
const bsim::NuChoice & GetNuChoice(void)
Definition: GDk2NuFlux.h:103
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
Loaders::FluxType flux
object containing MC flux information
Definition: Run.h:31
Double_t beam0y
y of beam center at start
Definition: dkmeta.h:107
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double sd(Eigen::VectorXd x)
double fTotalExposure
total exposure in POT
void produce(art::Event &evt)
GENIEGen(fhicl::ParameterSet const &pset)
int fPassEmptySpills
whether or not to allow for events with no interactions to be put into the event
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
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
int evt
caf::StandardRecord * sr
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
std::string horncfg
horn config e.g. "FHC/185A/LE/h1xoff=1mm"
Definition: dkmeta.h:99
const bsim::DkMeta & GetDkMeta(void)
Definition: GDk2NuFlux.h:105
bsim::DkMeta * dkmeta
int fCycle
cycle number in the MC generation
Double_t beam0x
x of beam center at start
Definition: dkmeta.h:106
double fPOTPerSpill
how many POT to simulate per spill
ProductID put(std::unique_ptr< PROD > &&)
double fTotalPOTLimit
total POT limit
An implementation of the GENIE GFluxI interface ("flux driver") encapsulating reading/processing the ...
Definition: GDk2NuFlux.h:69
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
Event generator information.
Definition: MCTruth.h:32
TStopwatch fStopwatch
keep track of how long it takes to run the job
Helper for AttenCurve.
Definition: Path.h:10
const bsim::Dk2Nu & GetDk2Nu(void)
Definition: GDk2NuFlux.h:104
Module to generate only pions from cosmic rays.
size_t size() const
Encapsulate the geometry of one entire detector (near, far, ndos)
A module to check the results from the Monte Carlo generator.
std::string FileBaseName() const
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
enum BeamMode string