17 #include "TDatabasePDG.h" 20 #include "TStopwatch.h" 29 #include "art_root_io/TFileService.h" 50 #include "dk2nu/genie/GDk2NuFlux.h" 51 #include "dk2nu/tree/NuChoice.h" 52 #include "dk2nu/tree/dk2nu.h" 54 #include "GENIE/Framework/EventGen/EventRecord.h" 170 produces<std::vector<simb::MCTruth>>();
171 produces<std::vector<simb::MCFlux>>();
172 produces<std::vector<simb::GTruth>>();
173 produces<sumdata::RunData, art::InRun>();
174 produces<sumdata::POTSummary, art::InSubRun>();
175 produces<art::Assns<simb::MCTruth, simb::MCFlux>>();
176 produces<art::Assns<simb::MCTruth, simb::GTruth>>();
177 produces<std::vector<sim::BeamGateInfo>>();
180 if (pset.get<std::string>(
"FluxType").find(
"dk2nu") != std::string::npos) {
181 produces<std::vector<bsim::Dk2Nu>>();
182 produces<std::vector<bsim::NuChoice>>();
183 produces<art::Assns<simb::MCTruth, bsim::Dk2Nu>>();
184 produces<art::Assns<simb::MCTruth, bsim::NuChoice>>();
187 std::string beam_type_name = pset.get<std::string>(
"BeamName");
189 if (beam_type_name ==
"numi")
193 else if (beam_type_name ==
"booster")
203 signed int temp_seed;
217 GENIEconfig.
put(
"RandomSeed", seed);
220 double detectorMass = 0;
223 if (pset.get<std::string>(
"FluxType").find(
"histogram") == 0 &&
224 pset.get<
double>(
"EventsPerSpill") == 0.0 && pset.get<
double>(
"POTPerSpill") > 0.0) {
227 detectorMass = geo->
TotalMass(pset.get<std::string>(
"TopVolume").c_str());
230 <<
"real time to calculate TotalMass: " << timer.RealTime() <<
" sec";
256 fGenerated[0] = tfs->make<TH1F>(
"fGenerated_necc",
"", 100, 0.0, 20.0);
257 fGenerated[1] = tfs->make<TH1F>(
"fGenerated_nebcc",
"", 100, 0.0, 20.0);
258 fGenerated[2] = tfs->make<TH1F>(
"fGenerated_nmcc",
"", 100, 0.0, 20.0);
259 fGenerated[3] = tfs->make<TH1F>(
"fGenerated_nmbcc",
"", 100, 0.0, 20.0);
260 fGenerated[4] = tfs->make<TH1F>(
"fGenerated_nnc",
"", 100, 0.0, 20.0);
261 fGenerated[5] = tfs->make<TH1F>(
"fGenerated_nbnc",
"", 100, 0.0, 20.0);
263 fDCosX = tfs->make<TH1F>(
"fDCosX",
";dx/ds", 200, -1., 1.);
264 fDCosY = tfs->make<TH1F>(
"fDCosY",
";dy/ds", 200, -1., 1.);
265 fDCosZ = tfs->make<TH1F>(
"fDCosZ",
";dz/ds", 200, -1., 1.);
267 fMuMomentum = tfs->make<TH1F>(
"fMuMomentum",
";p_{#mu} (GeV/c)", 500, 0., 50.);
268 fMuDCosX = tfs->make<TH1F>(
"fMuDCosX",
";dx/ds;", 200, -1., 1.);
269 fMuDCosY = tfs->make<TH1F>(
"fMuDCosY",
";dy/ds;", 200, -1., 1.);
270 fMuDCosZ = tfs->make<TH1F>(
"fMuDCosZ",
";dz/ds;", 200, -1., 1.);
272 fEMomentum = tfs->make<TH1F>(
"fEMomentum",
";p_{e} (GeV/c)", 500, 0., 50.);
273 fEDCosX = tfs->make<TH1F>(
"fEDCosX",
";dx/ds;", 200, -1., 1.);
274 fEDCosY = tfs->make<TH1F>(
"fEDCosY",
";dy/ds;", 200, -1., 1.);
275 fEDCosZ = tfs->make<TH1F>(
"fEDCosZ",
";dz/ds;", 200, -1., 1.);
277 fCCMode = tfs->make<TH1F>(
"fCCMode",
";CC Interaction Mode;", 4, 0., 4.);
278 fCCMode->GetXaxis()->SetBinLabel(1,
"QE");
279 fCCMode->GetXaxis()->SetBinLabel(2,
"Res");
280 fCCMode->GetXaxis()->SetBinLabel(3,
"DIS");
281 fCCMode->GetXaxis()->SetBinLabel(4,
"Coh");
282 fCCMode->GetXaxis()->CenterLabels();
284 fNCMode = tfs->make<TH1F>(
"fNCMode",
";NC Interaction Mode;", 4, 0., 4.);
285 fNCMode->GetXaxis()->SetBinLabel(1,
"QE");
286 fNCMode->GetXaxis()->SetBinLabel(2,
"Res");
287 fNCMode->GetXaxis()->SetBinLabel(3,
"DIS");
288 fNCMode->GetXaxis()->SetBinLabel(4,
"Coh");
289 fNCMode->GetXaxis()->CenterLabels();
291 fDeltaE = tfs->make<TH1F>(
"fDeltaE",
";#Delta E_{#nu} (GeV);", 200, -1., 1.);
292 fECons = tfs->make<TH1F>(
"fECons",
";#Delta E(#nu,lepton);", 500, -5., 5.);
299 int xdiv = TMath::Nint(2 * x / 5.);
300 int ydiv = TMath::Nint(2 * y / 5.);
301 int zdiv = TMath::Nint(2 * z / 5.);
303 fVertexX = tfs->make<TH1F>(
"fVertexX",
";x (cm)", xdiv, -0.1 *
x,
x);
304 fVertexY = tfs->make<TH1F>(
"fVertexY",
";y (cm)", ydiv, -
y,
y);
305 fVertexZ = tfs->make<TH1F>(
"fVertexZ",
";z (cm)", zdiv, -0.1 *
z,
z);
307 fVertexXY = tfs->make<TH2F>(
"fVertexXY",
";x (cm);y (cm)", xdiv, -0.1 *
x,
x, ydiv, -
y,
y);
309 tfs->make<TH2F>(
"fVertexXZ",
";z (cm);x (cm)", zdiv, -0.2 *
z,
z, xdiv, -0.1 *
x,
x);
310 fVertexYZ = tfs->make<TH2F>(
"fVertexYZ",
";z (cm);y (cm)", zdiv, -0.2 *
z,
z, ydiv, -
y,
y);
318 tfs->make<TH1F>(
"fVertexX",
";x (cm)", xdiv,
fVtxPosHistRange[0], fVtxPosHistRange[1]);
320 tfs->make<TH1F>(
"fVertexY",
";y (cm)", ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]);
322 tfs->make<TH1F>(
"fVertexZ",
";z (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5]);
331 fVtxPosHistRange[3]);
339 fVtxPosHistRange[1]);
347 fVtxPosHistRange[3]);
372 auto p = std::make_unique<sumdata::POTSummary>();
385 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
386 std::unique_ptr<std::vector<simb::MCFlux>> fluxcol(
new std::vector<simb::MCFlux>);
387 std::unique_ptr<std::vector<simb::GTruth>> gtruthcol(
new std::vector<simb::GTruth>);
388 std::unique_ptr<art::Assns<simb::MCTruth, simb::MCFlux>> tfassn(
390 std::unique_ptr<art::Assns<simb::MCTruth, simb::GTruth>> tgtassn(
392 std::unique_ptr<std::vector<sim::BeamGateInfo>> gateCollection(
393 new std::vector<sim::BeamGateInfo>);
395 std::unique_ptr<std::vector<bsim::Dk2Nu>> dk2nucol(
new std::vector<bsim::Dk2Nu>);
396 std::unique_ptr<std::vector<bsim::NuChoice>> nuchoicecol(
new std::vector<bsim::NuChoice>);
397 std::unique_ptr<art::Assns<simb::MCTruth, bsim::Dk2Nu>> dk2nuassn(
399 std::unique_ptr<art::Assns<simb::MCTruth, bsim::NuChoice>> nuchoiceassn(
402 genie::flux::GDk2NuFlux* dk2nuDriver =
404 while (truthcol->size() < 1) {
417 truthcol->push_back(truth);
418 fluxcol->push_back(flux);
419 gtruthcol->push_back(gTruth);
427 const bsim::Dk2Nu& dk2nuObj = dk2nuDriver->GetDk2Nu();
428 dk2nucol->push_back(dk2nuObj);
429 const bsim::NuChoice& nuchoiceObj = dk2nuDriver->GetNuChoice();
430 nuchoicecol->push_back(nuchoiceObj);
432 evt, *truthcol, *dk2nucol, *dk2nuassn, dk2nucol->size() - 1, dk2nucol->size());
437 nuchoicecol->size() - 1,
438 nuchoicecol->size());
447 log <<
"Found an interaction that is not represented by the interaction type code in " 452 log <<
"\nGENIE truth record:" 463 MF_LOG_DEBUG(
"GENIEGen") <<
"no events made for this spill but " 464 <<
"passing it on and ending the event anyway";
477 evt.
put(std::move(truthcol));
478 evt.
put(std::move(fluxcol));
479 evt.
put(std::move(gtruthcol));
480 evt.
put(std::move(tfassn));
481 evt.
put(std::move(tgtassn));
482 evt.
put(std::move(gateCollection));
486 evt.
put(std::move(dk2nucol));
487 evt.
put(std::move(nuchoicecol));
488 evt.
put(std::move(dk2nuassn));
489 evt.
put(std::move(nuchoiceassn));
498 int code = StatusCode;
502 case -1: ParticleStatusName =
"kIStUndefined";
break;
503 case 0: ParticleStatusName =
"kIStInitialState";
break;
504 case 1: ParticleStatusName =
"kIStStableFinalState";
break;
505 case 2: ParticleStatusName =
"kIStIntermediateState";
break;
506 case 3: ParticleStatusName =
"kIStDecayedState";
break;
507 case 11: ParticleStatusName =
"kIStNucleonTarget";
break;
508 case 12: ParticleStatusName =
"kIStDISPreFragmHadronicState";
break;
509 case 13: ParticleStatusName =
"kIStPreDecayResonantState";
break;
510 case 14: ParticleStatusName =
"kIStHadronInTheNucleus";
break;
511 case 15: ParticleStatusName =
"kIStFinalStateNuclearRemnant";
break;
512 case 16: ParticleStatusName =
"kIStNucleonClusterTarget";
break;
513 default: ParticleStatusName =
"Status Unknown";
521 std::string ReactionChannelName =
" ";
524 ReactionChannelName =
"kCC";
526 ReactionChannelName =
"kNC";
528 std::cout <<
"Current mode unknown!! " << std::endl;
531 ReactionChannelName +=
"_kQE";
533 ReactionChannelName +=
"_kRes";
535 ReactionChannelName +=
"_kDIS";
537 ReactionChannelName +=
"_kCoh";
539 std::cout <<
"interaction mode unknown!! " << std::endl;
541 return ReactionChannelName;
569 if (
id == -1) abort();
598 <<
"-----------> Particles in the Stack = " << mc.
NParticles() << std::endl
599 << std::setiosflags(
std::ios::left) << std::setw(20) <<
"PARTICLE" 600 << std::setiosflags(
std::ios::left) << std::setw(32) <<
"STATUS" << std::setw(18) <<
"E (GeV)" 601 << std::setw(18) <<
"m (GeV/c2)" << std::setw(18) <<
"Ek (GeV)" << std::endl
604 const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
609 std::string name = databasePDG->GetParticle(
part.PdgCode())->GetName();
610 int code =
part.StatusCode();
612 double mass =
part.Mass();
614 double Ek = (energy - mass);
615 if (status ==
"kIStStableFinalState" || status ==
"kIStHadronInTheNucleus")
618 << std::setiosflags(
std::ios::left) << std::setw(32) << status << std::setw(18) << energy
619 << std::setw(18) << mass << std::setw(18) << Ek << std::endl;
623 << std::setiosflags(
std::ios::left) << std::setw(32) << status << std::setw(18) << energy
624 << std::setw(18) << mass << std::endl;
633 if (
abs(
part.PdgCode()) == 11) {
641 else if (
abs(
part.PdgCode()) == 13) {
TH1F * fVertexX
vertex location of generated events in x
std::string ParticleStatus(int StatusCode)
double E(const int i=0) const
genie::GFluxI * GetFluxDriver(bool base=true)
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
TH2F * fVertexXZ
vertex location in xz
const simb::MCNeutrino & GetNeutrino() const
double Py(const int i=0) const
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.
TH1F * fDeltaE
difference in neutrino energy from MCTruth::Enu() vs TParticle
TH2F * fVertexXY
vertex location in xy
void beginSubRun(art::SubRun &sr)
std::string ReactionChannel(int ccnc, int mode)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void beginRun(art::Run &run)
TH1F * fMuDCosY
direction cosine of outgoing mu in y
const simb::MCParticle & Nu() const
bool Sample(simb::MCTruth &truth, simb::MCFlux &flux, simb::GTruth >ruth)
EDProducer(fhicl::ParameterSet const &pset)
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
double Px(const int i=0) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::string const & ROOTFile() const
Returns the full directory path to the geometry file source.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
offset to account for adding in Nuance codes to this enum
void endSubRun(art::SubRun &sr)
TGeoManager * ROOTGeoManager() const
Access to the ROOT geometry description manager.
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
TH1F * fNCMode
CC interaction mode.
object containing MC flux information
::sim::BeamType_t fBeamType
The width of a simulated "beam gate".
double TotalMass() const
Returns the total mass [kg] of the specified volume (default: world).
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
TH1F * fDCosZ
direction cosine in z
int InteractionType() const
void produce(art::Event &evt)
TH1F * fEDCosY
direction cosine of outgoing e in y
Utility functions to print MC truth information.
#define DEFINE_ART_MODULE(klass)
GENIEGen(fhicl::ParameterSet const &pset)
int fPassEmptySpills
whether or not to kill evnets with no interactions
TH1F * fCCMode
CC interaction mode.
double fPrevTotGoodPOT
Total good POT from subruns previous to current subrun.
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
double P(const int i=0) const
double fGlobalTimeOffset
keep track of how long it takes to run the job
TH1F * fEDCosX
direction cosine of outgoing e in x
bool fDefinedVtxHistRange
TH1F * fMuDCosX
direction cosine of outgoing mu in x
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.
TH1F * fGenerated[6]
Spectra as generated.
Wrapper for generating neutrino interactions with GENIE.
std::vector< double > fVtxPosHistRange
use defined hist range; it is useful to have for asymmetric ranges like in DP FD. ...
double fPrevTotPOT
The type of beam.
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
TH1F * fDCosY
direction cosine in y
const simb::MCParticle & GetParticle(int i) const
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
constexpr auto subRunFragment()
double Vx(const int i=0) const
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexYZ
vertex location in yz
Utility object to perform functions of association.
double TotalExposure() const
TH1F * fECons
histogram to determine if energy is conserved in the event
std::optional< T > get_if_present(std::string const &key) const
double Pz(const int i=0) const
double Vz(const int i=0) const
TH1F * fEDCosZ
direction cosine of outgoing e in z
TH1F * fDCosX
direction cosine in x
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Event generator information.
Event generator information.
Namespace collecting geometry-related classes utilities.
TH1F * fMuMomentum
momentum of outgoing muons
Event Generation using GENIE, cosmics or single particles.
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
TH1F * fVertexY
vertex location of generated events in y
BeamType_t
Defines category of beams to be stored in sim::BeamGateInfo.
void FillHistograms(simb::MCTruth mc)
TH1F * fVertexZ
vertex location of generated events in z
double Vy(const int i=0) const
void put(std::string const &key)
art framework interface to geometry description
A module to check the results from the Monte Carlo generator.
double fRandomTimeOffset
The start of a simulated "beam gate".