All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
evgen::GENIEGen Class Reference

A module to check the results from the Monte Carlo generator. More...

Inheritance diagram for evgen::GENIEGen:

Public Member Functions

 GENIEGen (fhicl::ParameterSet const &pset)
 
virtual ~GENIEGen ()
 
void produce (art::Event &evt)
 
void beginJob ()
 
void beginRun (art::Run &run)
 
void beginSubRun (art::SubRun &sr)
 
void endSubRun (art::SubRun &sr)
 

Private Member Functions

std::string ParticleStatus (int StatusCode)
 
std::string ReactionChannel (int ccnc, int mode)
 
void FillHistograms (simb::MCTruth mc)
 

Private Attributes

evgb::GENIEHelper * fGENIEHelp
 GENIEHelper object. More...
 
bool fDefinedVtxHistRange
 
std::vector< double > fVtxPosHistRange
 use defined hist range; it is useful to have for asymmetric ranges like in DP FD. More...
 
int fPassEmptySpills
 whether or not to kill evnets with no interactions More...
 
TStopwatch fStopwatch
 
double fGlobalTimeOffset
 keep track of how long it takes to run the job More...
 
double fRandomTimeOffset
 The start of a simulated "beam gate". More...
 
::sim::BeamType_t fBeamType
 The width of a simulated "beam gate". More...
 
double fPrevTotPOT
 The type of beam. More...
 
double fPrevTotGoodPOT
 Total good POT from subruns previous to current subrun. More...
 
TH1F * fGenerated [6]
 Spectra as generated. More...
 
TH1F * fVertexX
 vertex location of generated events in x More...
 
TH1F * fVertexY
 vertex location of generated events in y More...
 
TH1F * fVertexZ
 vertex location of generated events in z More...
 
TH2F * fVertexXY
 vertex location in xy More...
 
TH2F * fVertexXZ
 vertex location in xz More...
 
TH2F * fVertexYZ
 vertex location in yz More...
 
TH1F * fDCosX
 direction cosine in x More...
 
TH1F * fDCosY
 direction cosine in y More...
 
TH1F * fDCosZ
 direction cosine in z More...
 
TH1F * fMuMomentum
 momentum of outgoing muons More...
 
TH1F * fMuDCosX
 direction cosine of outgoing mu in x More...
 
TH1F * fMuDCosY
 direction cosine of outgoing mu in y More...
 
TH1F * fMuDCosZ
 direction cosine of outgoing mu in z More...
 
TH1F * fEMomentum
 momentum of outgoing electrons More...
 
TH1F * fEDCosX
 direction cosine of outgoing e in x More...
 
TH1F * fEDCosY
 direction cosine of outgoing e in y More...
 
TH1F * fEDCosZ
 direction cosine of outgoing e in z More...
 
TH1F * fCCMode
 CC interaction mode. More...
 
TH1F * fNCMode
 CC interaction mode. More...
 
TH1F * fDeltaE
 difference in neutrino energy from MCTruth::Enu() vs TParticle More...
 
TH1F * fECons
 histogram to determine if energy is conserved in the event More...
 

Detailed Description

A module to check the results from the Monte Carlo generator.

Note on random number generator

GENIE uses a TRandom generator for its purposes. Since art's RandomNumberGenerator service only provides CLHEP::HepRandomEngine, the standard LArSoft/art mechanism for handling the random stream can't be used. GENIEHelper, interface to GENIE provided by nugen, creates a TRandom that GENIE can use. It initializes it with a random seed read from RandomSeed configuration parameter. This and all the other parameters are inherited from the art module (that is, GENIEGen) configuration. LArSoft meddles with this mechanism to provide support for the standard "Seed" parameter and NuRandomService service.

Configuration parameters

As custom, if the random seed is not provided by the configuration, one is fetched from NuRandomService (if available), with the behaviour in lar::util::FetchRandomSeed().

Definition at line 81 of file GENIEGen_module.cc.

Constructor & Destructor Documentation

evgen::GENIEGen::GENIEGen ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 149 of file GENIEGen_module.cc.

150  : EDProducer{pset}
151  , fGENIEHelp(0)
152  , fDefinedVtxHistRange (pset.get< bool >("DefinedVtxHistRange"))
153  , fVtxPosHistRange (pset.get< std::vector<double> >("VtxPosHistRange"))
154  , fPassEmptySpills (pset.get< bool >("PassEmptySpills"))
155  , fGlobalTimeOffset(pset.get< double >("GlobalTimeOffset",0))
156  , fRandomTimeOffset(pset.get< double >("RandomTimeOffset",1600.)) // BNB default value
157  , fBeamType(::sim::kBNB)
158  {
159  fStopwatch.Start();
160 
161  produces< std::vector<simb::MCTruth> >();
162  produces< std::vector<simb::MCFlux> >();
163  produces< std::vector<simb::GTruth> >();
165  produces< sumdata::POTSummary, art::InSubRun >();
166  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
167  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
168  produces< std::vector<sim::BeamGateInfo> >();
169 
170  std::string beam_type_name = pset.get<std::string>("BeamName");
171 
172  if(beam_type_name == "numi")
173 
175 
176  else if(beam_type_name == "booster")
177 
179 
180  else
181 
182  fBeamType = ::sim::kUnknown;
183 
184  art::ServiceHandle<geo::Geometry const> geo;
185 
186  signed int temp_seed; // the seed read by GENIEHelper is a signed integer...
187  fhicl::ParameterSet GENIEconfig(pset);
188  if (!GENIEconfig.get_if_present("RandomSeed", temp_seed)) { // TODO use has_key() when it becomes available
189  // no RandomSeed specified; check for the LArSoft-style "Seed" instead:
190  // obtain the random seed from a service,
191  // unless overridden in configuration with key "Seed"
192  unsigned int seed;
193  if (!GENIEconfig.get_if_present("Seed", seed))
194  seed = art::ServiceHandle<rndm::NuRandomService>()->getSeed();
195 
196  // The seed is not passed to RandomNumberGenerator,
197  // since GENIE uses a TRandom generator that is owned by the GENIEHelper.
198  // Instead, we explicitly configure the random seed for GENIEHelper:
199  GENIEconfig.put("RandomSeed", seed);
200  } // if no RandomSeed present
201 
202  fGENIEHelp = new evgb::GENIEHelper(GENIEconfig,
203  geo->ROOTGeoManager(),
204  geo->ROOTFile(),
205  geo->TotalMass(pset.get< std::string>("TopVolume").c_str()));
206 
207  }
produces< sumdata::RunData, art::InRun >()
::sim::BeamType_t fBeamType
The width of a simulated &quot;beam gate&quot;.
NuMI.
Definition: BeamTypes.h:12
int fPassEmptySpills
whether or not to kill evnets with no interactions
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
double fGlobalTimeOffset
keep track of how long it takes to run the job
unsigned int seed
BNB.
Definition: BeamTypes.h:11
TStopwatch fStopwatch
std::vector< double > fVtxPosHistRange
use defined hist range; it is useful to have for asymmetric ranges like in DP FD. ...
double fRandomTimeOffset
The start of a simulated &quot;beam gate&quot;.
evgen::GENIEGen::~GENIEGen ( )
virtual

Definition at line 210 of file GENIEGen_module.cc.

211  {
212  if(fGENIEHelp) delete fGENIEHelp;
213  fStopwatch.Stop();
214  mf::LogInfo("GENIEProductionTime") << "real time to produce file: " << fStopwatch.RealTime();
215  }
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
TStopwatch fStopwatch

Member Function Documentation

void evgen::GENIEGen::beginJob ( )

Definition at line 218 of file GENIEGen_module.cc.

218  {
219  fGENIEHelp->Initialize();
220 
221  fPrevTotPOT = 0.;
222  fPrevTotGoodPOT = 0.;
223 
224  // Get access to the TFile service.
225  art::ServiceHandle<art::TFileService const> tfs;
226 
227  fGenerated[0] = tfs->make<TH1F>("fGenerated_necc","", 100, 0.0, 20.0);
228  fGenerated[1] = tfs->make<TH1F>("fGenerated_nebcc","", 100, 0.0, 20.0);
229  fGenerated[2] = tfs->make<TH1F>("fGenerated_nmcc","", 100, 0.0, 20.0);
230  fGenerated[3] = tfs->make<TH1F>("fGenerated_nmbcc","", 100, 0.0, 20.0);
231  fGenerated[4] = tfs->make<TH1F>("fGenerated_nnc","", 100, 0.0, 20.0);
232  fGenerated[5] = tfs->make<TH1F>("fGenerated_nbnc","", 100, 0.0, 20.0);
233 
234  fDCosX = tfs->make<TH1F>("fDCosX", ";dx/ds", 200, -1., 1.);
235  fDCosY = tfs->make<TH1F>("fDCosY", ";dy/ds", 200, -1., 1.);
236  fDCosZ = tfs->make<TH1F>("fDCosZ", ";dz/ds", 200, -1., 1.);
237 
238  fMuMomentum = tfs->make<TH1F>("fMuMomentum", ";p_{#mu} (GeV/c)", 500, 0., 50.);
239  fMuDCosX = tfs->make<TH1F>("fMuDCosX", ";dx/ds;", 200, -1., 1.);
240  fMuDCosY = tfs->make<TH1F>("fMuDCosY", ";dy/ds;", 200, -1., 1.);
241  fMuDCosZ = tfs->make<TH1F>("fMuDCosZ", ";dz/ds;", 200, -1., 1.);
242 
243  fEMomentum = tfs->make<TH1F>("fEMomentum", ";p_{e} (GeV/c)", 500, 0., 50.);
244  fEDCosX = tfs->make<TH1F>("fEDCosX", ";dx/ds;", 200, -1., 1.);
245  fEDCosY = tfs->make<TH1F>("fEDCosY", ";dy/ds;", 200, -1., 1.);
246  fEDCosZ = tfs->make<TH1F>("fEDCosZ", ";dz/ds;", 200, -1., 1.);
247 
248  fCCMode = tfs->make<TH1F>("fCCMode", ";CC Interaction Mode;", 4, 0., 4.);
249  fCCMode->GetXaxis()->SetBinLabel(1, "QE");
250  fCCMode->GetXaxis()->SetBinLabel(2, "Res");
251  fCCMode->GetXaxis()->SetBinLabel(3, "DIS");
252  fCCMode->GetXaxis()->SetBinLabel(4, "Coh");
253  fCCMode->GetXaxis()->CenterLabels();
254 
255  fNCMode = tfs->make<TH1F>("fNCMode", ";NC Interaction Mode;", 4, 0., 4.);
256  fNCMode->GetXaxis()->SetBinLabel(1, "QE");
257  fNCMode->GetXaxis()->SetBinLabel(2, "Res");
258  fNCMode->GetXaxis()->SetBinLabel(3, "DIS");
259  fNCMode->GetXaxis()->SetBinLabel(4, "Coh");
260  fNCMode->GetXaxis()->CenterLabels();
261 
262  fDeltaE = tfs->make<TH1F>("fDeltaE", ";#Delta E_{#nu} (GeV);", 200, -1., 1.);
263  fECons = tfs->make<TH1F>("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.);
264 
265  if (fDefinedVtxHistRange == false)
266  {
267  art::ServiceHandle<geo::Geometry const> geo;
268  double x = 2.1*geo->DetHalfWidth();
269  double y = 2.1*geo->DetHalfHeight();
270  double z = 2.*geo->DetLength();
271  int xdiv = TMath::Nint(2*x/5.);
272  int ydiv = TMath::Nint(2*y/5.);
273  int zdiv = TMath::Nint(2*z/5.);
274 
275  fVertexX = tfs->make<TH1F>("fVertexX", ";x (cm)", xdiv, -0.1*x, x);
276  fVertexY = tfs->make<TH1F>("fVertexY", ";y (cm)", ydiv, -y, y);
277  fVertexZ = tfs->make<TH1F>("fVertexZ", ";z (cm)", zdiv, -0.1*z, z);
278 
279  fVertexXY = tfs->make<TH2F>("fVertexXY", ";x (cm);y (cm)", xdiv, -0.1*x, x, ydiv, -y, y);
280  fVertexXZ = tfs->make<TH2F>("fVertexXZ", ";z (cm);x (cm)", zdiv, -0.2*z, z, xdiv, -0.1*x, x);
281  fVertexYZ = tfs->make<TH2F>("fVertexYZ", ";z (cm);y (cm)", zdiv, -0.2*z, z, ydiv, -y, y);
282  }
283  else
284  {
285  int xdiv = TMath::Nint((fVtxPosHistRange[1]-fVtxPosHistRange[0])/5.);
286  int ydiv = TMath::Nint((fVtxPosHistRange[3]-fVtxPosHistRange[2])/5.);
287  int zdiv = TMath::Nint((fVtxPosHistRange[5]-fVtxPosHistRange[4])/5.);
288 
289  fVertexX = tfs->make<TH1F>("fVertexX", ";x (cm)", xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1]);
290  fVertexY = tfs->make<TH1F>("fVertexY", ";y (cm)", ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]);
291  fVertexZ = tfs->make<TH1F>("fVertexZ", ";z (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5]);
292 
293  fVertexXY = tfs->make<TH2F>("fVertexXY", ";x (cm);y (cm)", xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1], ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]);
294  fVertexXZ = tfs->make<TH2F>("fVertexXZ", ";z (cm);x (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5], xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1]);
295  fVertexYZ = tfs->make<TH2F>("fVertexYZ", ";z (cm);y (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5], ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]);
296  }
297 
298  }
TH1F * fVertexX
vertex location of generated events in x
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
TH2F * fVertexXZ
vertex location in xz
process_name opflash particleana ie ie ie z
process_name opflash particleana ie x
TH1F * fDeltaE
difference in neutrino energy from MCTruth::Enu() vs TParticle
TH2F * fVertexXY
vertex location in xy
TH1F * fMuDCosY
direction cosine of outgoing mu in y
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosZ
direction cosine in z
TH1F * fEDCosY
direction cosine of outgoing e in y
process_name opflash particleana ie ie y
TH1F * fCCMode
CC interaction mode.
double fPrevTotGoodPOT
Total good POT from subruns previous to current subrun.
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
TH1F * fEDCosX
direction cosine of outgoing e in x
TH1F * fMuDCosX
direction cosine of outgoing mu in x
TH1F * fGenerated[6]
Spectra as generated.
double fPrevTotPOT
The type of beam.
TH1F * fDCosY
direction cosine in y
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexYZ
vertex location in yz
TH1F * fECons
histogram to determine if energy is conserved in the event
TH1F * fEDCosZ
direction cosine of outgoing e in z
TH1F * fDCosX
direction cosine in x
art::ServiceHandle< art::TFileService > tfs
TH1F * fMuMomentum
momentum of outgoing muons
TH1F * fVertexY
vertex location of generated events in y
std::vector< double > fVtxPosHistRange
use defined hist range; it is useful to have for asymmetric ranges like in DP FD. ...
TH1F * fVertexZ
vertex location of generated events in z
void evgen::GENIEGen::beginRun ( art::Run &  run)

Definition at line 301 of file GENIEGen_module.cc.

302  {
303  art::ServiceHandle<geo::Geometry const> geo;
304  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
305  }
void evgen::GENIEGen::beginSubRun ( art::SubRun &  sr)

Definition at line 308 of file GENIEGen_module.cc.

309  {
310 
311  fPrevTotPOT = fGENIEHelp->TotalExposure();
312  fPrevTotGoodPOT = fGENIEHelp->TotalExposure();
313 
314  return;
315  }
double fPrevTotGoodPOT
Total good POT from subruns previous to current subrun.
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
double fPrevTotPOT
The type of beam.
void evgen::GENIEGen::endSubRun ( art::SubRun &  sr)

Definition at line 318 of file GENIEGen_module.cc.

319  {
320 
321  auto p = std::make_unique<sumdata::POTSummary>();
322 
323  p->totpot = fGENIEHelp->TotalExposure() - fPrevTotPOT;
324  p->totgoodpot = fGENIEHelp->TotalExposure() - fPrevTotGoodPOT;
325 
326  sr.put(std::move(p));
327 
328  return;
329  }
pdgs p
Definition: selectors.fcl:22
double fPrevTotGoodPOT
Total good POT from subruns previous to current subrun.
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
double fPrevTotPOT
The type of beam.
void evgen::GENIEGen::FillHistograms ( simb::MCTruth  mc)
private

< fill the vertex histograms from the neutrino - that is always particle 0 in the list

look for the outgoing lepton in the particle stack just interested in the first one

Definition at line 482 of file GENIEGen_module.cc.

483  {
484  // Decide which histograms to put the spectrum in
485  int id = -1;
486  if (mc.GetNeutrino().CCNC()==simb::kCC) {
487  fCCMode->Fill(mc.GetNeutrino().Mode());
488  if (mc.GetNeutrino().Nu().PdgCode() == 12) id = 0;
489  else if (mc.GetNeutrino().Nu().PdgCode() == -12) id = 1;
490  else if (mc.GetNeutrino().Nu().PdgCode() == 14) id = 2;
491  else if (mc.GetNeutrino().Nu().PdgCode() == -14) id = 3;
492  else return;
493  }
494  else {
495  fNCMode->Fill(mc.GetNeutrino().Mode());
496  if (mc.GetNeutrino().Nu().PdgCode() > 0) id = 4;
497  else id = 5;
498  }
499  if (id==-1) abort();
500 
501  // Fill the specta histograms
502  fGenerated[id]->Fill(mc.GetNeutrino().Nu().E() );
503 
504  ///< fill the vertex histograms from the neutrino - that is always
505  ///< particle 0 in the list
506  simb::MCNeutrino mcnu = mc.GetNeutrino();
507  const simb::MCParticle nu = mcnu.Nu();
508 
509  fVertexX->Fill(nu.Vx());
510  fVertexY->Fill(nu.Vy());
511  fVertexZ->Fill(nu.Vz());
512 
513  fVertexXY->Fill(nu.Vx(), nu.Vy());
514  fVertexXZ->Fill(nu.Vz(), nu.Vx());
515  fVertexYZ->Fill(nu.Vz(), nu.Vy());
516 
517  double mom = nu.P();
518  if(std::abs(mom) > 0.){
519  fDCosX->Fill(nu.Px()/mom);
520  fDCosY->Fill(nu.Py()/mom);
521  fDCosZ->Fill(nu.Pz()/mom);
522  }
523 
524 
525  MF_LOG_DEBUG("GENIEInteractionInformation")
526  << std::endl
527  << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode())
528  << std::endl
529  << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl
530  << std::setiosflags(std::ios::left)
531  << std::setw(20) << "PARTICLE"
532  << std::setiosflags(std::ios::left)
533  << std::setw(32) << "STATUS"
534  << std::setw(18) << "E (GeV)"
535  << std::setw(18) << "m (GeV/c2)"
536  << std::setw(18) << "Ek (GeV)"
537  << std::endl << std::endl;
538 
539  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
540 
541  // Loop over the particle stack for this event
542  for(int i = 0; i < mc.NParticles(); ++i){
543  simb::MCParticle part(mc.GetParticle(i));
544  std::string name = databasePDG->GetParticle(part.PdgCode())->GetName();
545  int code = part.StatusCode();
546  std::string status = ParticleStatus(code);
547  double mass = part.Mass();
548  double energy = part.E();
549  double Ek = (energy-mass); // Kinetic Energy (GeV)
550  if(status=="kIStStableFinalState"||status=="kIStHadronInTheNucleus")
551  MF_LOG_DEBUG("GENIEFinalState")
552  << std::setiosflags(std::ios::left) << std::setw(20) << name
553  << std::setiosflags(std::ios::left) << std::setw(32) <<status
554  << std::setw(18)<< energy
555  << std::setw(18)<< mass
556  << std::setw(18)<< Ek <<std::endl;
557  else
558  MF_LOG_DEBUG("GENIEFinalState")
559  << std::setiosflags(std::ios::left) << std::setw(20) << name
560  << std::setiosflags(std::ios::left) << std::setw(32) << status
561  << std::setw(18) << energy
562  << std::setw(18) << mass <<std::endl;
563  }
564 
565 
566  if(mc.GetNeutrino().CCNC() == simb::kCC){
567 
568  ///look for the outgoing lepton in the particle stack
569  ///just interested in the first one
570  for(int i = 0; i < mc.NParticles(); ++i){
571  simb::MCParticle part(mc.GetParticle(i));
572  if(abs(part.PdgCode()) == 11){
573  fEMomentum->Fill(part.P());
574  fEDCosX->Fill(part.Px()/part.P());
575  fEDCosY->Fill(part.Py()/part.P());
576  fEDCosZ->Fill(part.Pz()/part.P());
577  fECons->Fill(nu.E() - part.E());
578  break;
579  }
580  else if(abs(part.PdgCode()) == 13){
581  fMuMomentum->Fill(part.P());
582  fMuDCosX->Fill(part.Px()/part.P());
583  fMuDCosY->Fill(part.Py()/part.P());
584  fMuDCosZ->Fill(part.Pz()/part.P());
585  fECons->Fill(nu.E() - part.E());
586  break;
587  }
588  }// end loop over particles
589  }//end if CC interaction
590 
591  return;
592  }
TH1F * fVertexX
vertex location of generated events in x
std::string ParticleStatus(int StatusCode)
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
TH2F * fVertexXZ
vertex location in xz
TH2F * fVertexXY
vertex location in xy
std::string ReactionChannel(int ccnc, int mode)
TH1F * fMuDCosY
direction cosine of outgoing mu in y
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosZ
direction cosine in z
T abs(T value)
Charged-current interactions.
Definition: IPrediction.h:38
TH1F * fEDCosY
direction cosine of outgoing e in y
TH1F * fCCMode
CC interaction mode.
TH1F * fEDCosX
direction cosine of outgoing e in x
TH1F * fMuDCosX
direction cosine of outgoing mu in x
TH1F * fGenerated[6]
Spectra as generated.
walls no left
Definition: selectors.fcl:105
TH1F * fDCosY
direction cosine in y
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexYZ
vertex location in yz
TH1F * fECons
histogram to determine if energy is conserved in the event
float mass
Definition: dedx.py:47
TH1F * fEDCosZ
direction cosine of outgoing e in z
TH1F * fDCosX
direction cosine in x
then echo fcl name
TH1F * fMuMomentum
momentum of outgoing muons
BEGIN_PROLOG SN nu
TH1F * fVertexY
vertex location of generated events in y
TH1F * fVertexZ
vertex location of generated events in z
std::string evgen::GENIEGen::ParticleStatus ( int  StatusCode)
private

Definition at line 411 of file GENIEGen_module.cc.

412  {
413  int code = StatusCode;
414  std::string ParticleStatusName;
415 
416  switch(code)
417  {
418  case -1:
419  ParticleStatusName = "kIStUndefined";
420  break;
421  case 0:
422  ParticleStatusName = "kIStInitialState";
423  break;
424  case 1:
425  ParticleStatusName = "kIStStableFinalState";
426  break;
427  case 2:
428  ParticleStatusName = "kIStIntermediateState";
429  break;
430  case 3:
431  ParticleStatusName = "kIStDecayedState";
432  break;
433  case 11:
434  ParticleStatusName = "kIStNucleonTarget";
435  break;
436  case 12:
437  ParticleStatusName = "kIStDISPreFragmHadronicState";
438  break;
439  case 13:
440  ParticleStatusName = "kIStPreDecayResonantState";
441  break;
442  case 14:
443  ParticleStatusName = "kIStHadronInTheNucleus";
444  break;
445  case 15:
446  ParticleStatusName = "kIStFinalStateNuclearRemnant";
447  break;
448  case 16:
449  ParticleStatusName = "kIStNucleonClusterTarget";
450  break;
451  default:
452  ParticleStatusName = "Status Unknown";
453  }
454  return ParticleStatusName;
455  }
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
void evgen::GENIEGen::produce ( art::Event &  evt)

Definition at line 332 of file GENIEGen_module.cc.

333  {
334  std::unique_ptr< std::vector<simb::MCTruth> > truthcol (new std::vector<simb::MCTruth>);
335  std::unique_ptr< std::vector<simb::MCFlux> > fluxcol (new std::vector<simb::MCFlux >);
336  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol (new std::vector<simb::GTruth >);
337  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > tfassn(new art::Assns<simb::MCTruth, simb::MCFlux>);
338  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > tgtassn(new art::Assns<simb::MCTruth, simb::GTruth>);
339  std::unique_ptr< std::vector<sim::BeamGateInfo> > gateCollection(new std::vector<sim::BeamGateInfo>);
340 
341  while(truthcol->size() < 1){
342  while(!fGENIEHelp->Stop()){
343 
344  simb::MCTruth truth;
345  simb::MCFlux flux;
346  simb::GTruth gTruth;
347 
348  // GENIEHelper returns a false in the sample method if
349  // either no neutrino was generated, or the interaction
350  // occurred beyond the detector's z extent - ie something we
351  // would never see anyway.
352  if(fGENIEHelp->Sample(truth, flux, gTruth)){
353 
354  truthcol ->push_back(truth);
355  fluxcol ->push_back(flux);
356  gtruthcol->push_back(gTruth);
357  auto const truthPtr = art::PtrMaker<simb::MCTruth>{evt}(truthcol->size() - 1);
358  tfassn->addSingle(truthPtr, art::PtrMaker<simb::MCFlux>{evt}(fluxcol->size() - 1));
359  tgtassn->addSingle(truthPtr, art::PtrMaker<simb::GTruth>{evt}(gtruthcol->size() - 1));
360 
361  FillHistograms(truth);
362 
363  // check that the process code is not unsupported by GENIE
364  // (see issue #18025 for reference);
365  // if it is, print all the information we can about this truth record
366  if (truth.NeutrinoSet() && (truth.GetNeutrino().InteractionType() == simb::kNuanceOffset)) {
367  mf::LogWarning log("GENIEmissingProcessMapping");
368  log << "Found an interaction that is not represented by the interaction type code in GENIE:"
369  "\nMCTruth record:"
370  "\n"
371  ;
372  sim::dump::DumpMCTruth(log, truth, 2U); // 2 trajectory points per line
373  log <<
374  "\nGENIE truth record:"
375  "\n"
376  ;
377  sim::dump::DumpGTruth(log, gTruth);
378  } // if
379 
380  }// end if genie was able to make an event
381 
382  }// end event generation loop
383 
384  // check to see if we are to pass empty spills
385  if(truthcol->size() < 1 && fPassEmptySpills){
386  MF_LOG_DEBUG("GENIEGen") << "no events made for this spill but "
387  << "passing it on and ending the event anyway";
388  break;
389  }
390 
391  }// end loop while no interactions are made
392 
393  // Create a simulated "beam gate" for these neutrino events.
394  // We're creating a vector of these because, in a
395  // distant-but-possible future, we may be generating more than one
396  // beam gate within a simulated time window.
397  gateCollection->push_back(sim::BeamGateInfo( fGlobalTimeOffset, fRandomTimeOffset, fBeamType ));
398 
399  // put the collections in the event
400  evt.put(std::move(truthcol));
401  evt.put(std::move(fluxcol));
402  evt.put(std::move(gtruthcol));
403  evt.put(std::move(tfassn));
404  evt.put(std::move(tgtassn));
405  evt.put(std::move(gateCollection));
406 
407  return;
408  }
void DumpGTruth(Stream &&out, simb::GTruth const &truth, std::string indent, std::string firstIndent)
Dumps the content of the GENIE truth in the output stream.
Definition: MCDumpers.h:379
offset to account for adding in Nuance codes to this enum
Definition: SREnums.h:85
::sim::BeamType_t fBeamType
The width of a simulated &quot;beam gate&quot;.
int fPassEmptySpills
whether or not to kill evnets with no interactions
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
double fGlobalTimeOffset
keep track of how long it takes to run the job
void DumpMCTruth(Stream &&out, simb::MCTruth const &truth, unsigned int pointsPerLine, std::string indent, std::string firstIndent)
Dumps the content of the specified MC truth in the output stream.
Definition: MCDumpers.h:346
TCEvent evt
Definition: DataStructs.cxx:8
void FillHistograms(simb::MCTruth mc)
double fRandomTimeOffset
The start of a simulated &quot;beam gate&quot;.
std::string evgen::GENIEGen::ReactionChannel ( int  ccnc,
int  mode 
)
private

Definition at line 458 of file GENIEGen_module.cc.

459  {
460  std::string ReactionChannelName=" ";
461 
462  if(ccnc==0)
463  ReactionChannelName = "kCC";
464  else if(ccnc==1)
465  ReactionChannelName = "kNC";
466  else std::cout<<"Current mode unknown!! "<<std::endl;
467 
468  if(mode==0)
469  ReactionChannelName += "_kQE";
470  else if(mode==1)
471  ReactionChannelName += "_kRes";
472  else if(mode==2)
473  ReactionChannelName += "_kDIS";
474  else if(mode==3)
475  ReactionChannelName += "_kCoh";
476  else std::cout<<"interaction mode unknown!! "<<std::endl;
477 
478  return ReactionChannelName;
479  }
const char mode
Definition: noise_ana.cxx:20
BEGIN_PROLOG could also be cout

Member Data Documentation

::sim::BeamType_t evgen::GENIEGen::fBeamType
private

The width of a simulated "beam gate".

Definition at line 108 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fCCMode
private

CC interaction mode.

Definition at line 137 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fDCosX
private

direction cosine in x

Definition at line 123 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fDCosY
private

direction cosine in y

Definition at line 124 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fDCosZ
private

direction cosine in z

Definition at line 125 of file GENIEGen_module.cc.

bool evgen::GENIEGen::fDefinedVtxHistRange
private

Definition at line 100 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fDeltaE
private

difference in neutrino energy from MCTruth::Enu() vs TParticle

Definition at line 140 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fECons
private

histogram to determine if energy is conserved in the event

Definition at line 141 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fEDCosX
private

direction cosine of outgoing e in x

Definition at line 133 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fEDCosY
private

direction cosine of outgoing e in y

Definition at line 134 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fEDCosZ
private

direction cosine of outgoing e in z

Definition at line 135 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fEMomentum
private

momentum of outgoing electrons

Definition at line 132 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fGenerated[6]
private

Spectra as generated.

Definition at line 113 of file GENIEGen_module.cc.

evgb::GENIEHelper* evgen::GENIEGen::fGENIEHelp
private

GENIEHelper object.

Definition at line 99 of file GENIEGen_module.cc.

double evgen::GENIEGen::fGlobalTimeOffset
private

keep track of how long it takes to run the job

Definition at line 106 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fMuDCosX
private

direction cosine of outgoing mu in x

Definition at line 128 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fMuDCosY
private

direction cosine of outgoing mu in y

Definition at line 129 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fMuDCosZ
private

direction cosine of outgoing mu in z

Definition at line 130 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fMuMomentum
private

momentum of outgoing muons

Definition at line 127 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fNCMode
private

CC interaction mode.

Definition at line 138 of file GENIEGen_module.cc.

int evgen::GENIEGen::fPassEmptySpills
private

whether or not to kill evnets with no interactions

Definition at line 103 of file GENIEGen_module.cc.

double evgen::GENIEGen::fPrevTotGoodPOT
private

Total good POT from subruns previous to current subrun.

Definition at line 111 of file GENIEGen_module.cc.

double evgen::GENIEGen::fPrevTotPOT
private

The type of beam.

Total POT from subruns previous to current subrun

Definition at line 110 of file GENIEGen_module.cc.

double evgen::GENIEGen::fRandomTimeOffset
private

The start of a simulated "beam gate".

Definition at line 107 of file GENIEGen_module.cc.

TStopwatch evgen::GENIEGen::fStopwatch
private

Definition at line 104 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fVertexX
private

vertex location of generated events in x

Definition at line 115 of file GENIEGen_module.cc.

TH2F* evgen::GENIEGen::fVertexXY
private

vertex location in xy

Definition at line 119 of file GENIEGen_module.cc.

TH2F* evgen::GENIEGen::fVertexXZ
private

vertex location in xz

Definition at line 120 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fVertexY
private

vertex location of generated events in y

Definition at line 116 of file GENIEGen_module.cc.

TH2F* evgen::GENIEGen::fVertexYZ
private

vertex location in yz

Definition at line 121 of file GENIEGen_module.cc.

TH1F* evgen::GENIEGen::fVertexZ
private

vertex location of generated events in z

Definition at line 117 of file GENIEGen_module.cc.

std::vector< double > evgen::GENIEGen::fVtxPosHistRange
private

use defined hist range; it is useful to have for asymmetric ranges like in DP FD.

Definition at line 101 of file GENIEGen_module.cc.


The documentation for this class was generated from the following file: