20 #include "TDatabasePDG.h" 23 #include "TStopwatch.h" 25 #include "CLHEP/Random/RandFlat.h" 33 #include "art_root_io/TFileService.h" 105 ,
fNdkFile{pset.get<std::string>(
"NdkFile")}
113 produces<std::vector<simb::MCTruth>>();
114 produces<sumdata::RunData, art::InRun>();
127 fGenerated[0] = tfs->make<TH1F>(
"fGenerated_necc",
"", 100, 0.0, 20.0);
128 fGenerated[1] = tfs->make<TH1F>(
"fGenerated_nebcc",
"", 100, 0.0, 20.0);
129 fGenerated[2] = tfs->make<TH1F>(
"fGenerated_nmcc",
"", 100, 0.0, 20.0);
130 fGenerated[3] = tfs->make<TH1F>(
"fGenerated_nmbcc",
"", 100, 0.0, 20.0);
131 fGenerated[4] = tfs->make<TH1F>(
"fGenerated_nnc",
"", 100, 0.0, 20.0);
132 fGenerated[5] = tfs->make<TH1F>(
"fGenerated_nbnc",
"", 100, 0.0, 20.0);
134 fDCosX = tfs->make<TH1F>(
"fDCosX",
";dx/ds", 200, -1., 1.);
135 fDCosY = tfs->make<TH1F>(
"fDCosY",
";dy/ds", 200, -1., 1.);
136 fDCosZ = tfs->make<TH1F>(
"fDCosZ",
";dz/ds", 200, -1., 1.);
138 fMuMomentum = tfs->make<TH1F>(
"fMuMomentum",
";p_{#mu} (GeV/c)", 500, 0., 50.);
139 fMuDCosX = tfs->make<TH1F>(
"fMuDCosX",
";dx/ds;", 200, -1., 1.);
140 fMuDCosY = tfs->make<TH1F>(
"fMuDCosY",
";dy/ds;", 200, -1., 1.);
141 fMuDCosZ = tfs->make<TH1F>(
"fMuDCosZ",
";dz/ds;", 200, -1., 1.);
143 fEMomentum = tfs->make<TH1F>(
"fEMomentum",
";p_{e} (GeV/c)", 500, 0., 50.);
144 fEDCosX = tfs->make<TH1F>(
"fEDCosX",
";dx/ds;", 200, -1., 1.);
145 fEDCosY = tfs->make<TH1F>(
"fEDCosY",
";dy/ds;", 200, -1., 1.);
146 fEDCosZ = tfs->make<TH1F>(
"fEDCosZ",
";dz/ds;", 200, -1., 1.);
148 fCCMode = tfs->make<TH1F>(
"fCCMode",
";CC Interaction Mode;", 5, 0., 5.);
149 fCCMode->GetXaxis()->SetBinLabel(1,
"QE");
150 fCCMode->GetXaxis()->SetBinLabel(2,
"Res");
151 fCCMode->GetXaxis()->SetBinLabel(3,
"DIS");
152 fCCMode->GetXaxis()->SetBinLabel(4,
"Coh");
153 fCCMode->GetXaxis()->SetBinLabel(5,
"kInverseMuDecay");
154 fCCMode->GetXaxis()->CenterLabels();
156 fNCMode = tfs->make<TH1F>(
"fNCMode",
";NC Interaction Mode;", 5, 0., 5.);
157 fNCMode->GetXaxis()->SetBinLabel(1,
"QE");
158 fNCMode->GetXaxis()->SetBinLabel(2,
"Res");
159 fNCMode->GetXaxis()->SetBinLabel(3,
"DIS");
160 fNCMode->GetXaxis()->SetBinLabel(4,
"Coh");
161 fNCMode->GetXaxis()->SetBinLabel(5,
"kNuElectronElastic");
162 fNCMode->GetXaxis()->CenterLabels();
164 fECons = tfs->make<TH1F>(
"fECons",
";#Delta E(#nu,lepton);", 500, -5., 5.);
170 int xdiv = TMath::Nint(2 * x / 5.);
171 int ydiv = TMath::Nint(2 * y / 5.);
172 int zdiv = TMath::Nint(2 * z / 5.);
174 fVertexX = tfs->make<TH1F>(
"fVertexX",
";x (cm)", xdiv, -
x,
x);
175 fVertexY = tfs->make<TH1F>(
"fVertexY",
";y (cm)", ydiv, -
y,
y);
176 fVertexZ = tfs->make<TH1F>(
"fVertexZ",
";z (cm)", zdiv, -0.2 *
z,
z);
178 fVertexXY = tfs->make<TH2F>(
"fVertexXY",
";x (cm);y (cm)", xdiv, -
x,
x, ydiv, -
y,
y);
179 fVertexXZ = tfs->make<TH2F>(
"fVertexXZ",
";z (cm);x (cm)", zdiv, -0.2 *
z,
z, xdiv, -
x,
x);
180 fVertexYZ = tfs->make<TH2F>(
"fVertexYZ",
";z (cm);y (cm)", zdiv, -0.2 *
z,
z, ydiv, -
y,
y);
200 std::cout << std::endl;
201 std::cout <<
"------------------------------------------------------------------------------" 203 std::cout <<
"event : " << evt.
id().
event() << std::endl;
205 std::string name, k, dollar;
211 int Idx, Ist, PDG, Mother1, Mother2, Daughter1, Daughter2;
212 double Px, Py, Pz,
E, m;
213 std::string p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14;
216 std::string primary(
"primary");
217 int FirstMother = -1;
223 TLorentzVector Neutrino;
224 TLorentzVector Lepton;
225 TLorentzVector Target;
227 TLorentzVector Hadron4mom;
231 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
237 double const fvCut{5.0};
247 auto const world = tpc.GetCenter();
248 if (minx > world.X() - tpc.HalfWidth()) minx = world.X() - tpc.HalfWidth();
249 if (maxx < world.X() + tpc.HalfWidth()) maxx = world.X() + tpc.HalfWidth();
250 if (miny > world.Y() - tpc.HalfHeight()) miny = world.Y() - tpc.HalfHeight();
251 if (maxy < world.Y() + tpc.HalfHeight()) maxy = world.Y() + tpc.HalfHeight();
252 if (minz > world.Z() - tpc.Length() / 2.) minz = world.Z() - tpc.Length() / 2.;
253 if (maxz < world.Z() + tpc.Length() / 2.) maxz = world.Z() + tpc.Length() / 2.;
257 double X0 = 0.0 + flat.fire(minx + fvCut, maxx - fvCut);
258 double Y0 = 0.0 + flat.fire(miny + fvCut, maxy - fvCut);
259 double Z0 = 0.0 + flat.fire(minz + fvCut, maxz - fvCut);
261 std::cout <<
"NDKGen_module: X, Y, Z of vtx: " << X0 <<
", " << Y0 <<
", " << Z0 << std::endl;
265 if (!
fEventFile.good()) std::cout <<
"NdkFile: Problem reading Ndk file" << std::endl;
269 if (k.find(
"** Event:") != std::string::npos) {
270 std::istringstream
in;
274 in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >>
276 std::cout <<
"Genie Evt " << GenieEvt <<
" art evt " << evt.
id().
event() <<
"\n";
279 if (GenieEvt + 1 != static_cast<int>(evt.
id().
event()))
283 if (!k.compare(0, 25,
"GENIE Interaction Summary"))
285 if (k.compare(0, 1,
"|") || k.compare(1, 2,
" "))
287 if (k.find(
"Fin-Init") != std::string::npos)
continue;
288 if (k.find(
"Ar") != std::string::npos)
continue;
289 if (k.find(
"Cl") != std::string::npos)
continue;
290 if (k.find(
"HadrBlob") != std::string::npos)
continue;
291 if (k.find(
"NucBindE") != std::string::npos)
continue;
292 if (k.find(
"FLAGS") != std::string::npos)
break;
293 if (k.find(
"Vertex") != std::string::npos)
break;
295 if (!k.compare(26, 1,
"1"))
298 std::istringstream
in;
302 in >> p1 >> Idx >> p2 >> Name >> p3 >> Ist >> p4 >> PDG >> p5 >> Mother1 >> p6 >>
303 Mother2 >> p7 >> Daughter1 >> p8 >> Daughter2 >> p9 >> Px >> p10 >> Py >> p11 >> Pz >>
304 p12 >> E >> p13 >> m >> p14;
305 if (Ist != 1)
continue;
307 std::cout <<
"PDG = " << PDG << std::endl;
309 const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
310 const TParticlePDG* definition = databasePDG->GetParticle(PDG);
311 Mass = definition->Mass();
312 if (E - Mass < 0.001)
continue;
318 P = std::sqrt(pow(E, 2.) - pow(Mass, 2.));
319 std::cout <<
"Momentum = " << P << std::endl;
321 TLorentzVector pos(X0, Y0, Z0, 0);
325 TLorentzVector mom(Px, Py, Pz, E);
336 std::cout <<
"NDKGen.cxx: Putting " << truth.
NParticles() <<
" tracks on stack." << std::endl;
337 truthcol->push_back(truth);
338 evt.
put(std::move(truthcol));
344 switch (StatusCode) {
345 case 0:
return "kIStInitialState";
346 case 1:
return "kIStFinalState";
347 case 11:
return "kIStNucleonTarget";
348 default:
return "Status Unknown";
355 std::string ReactionChannelName =
" ";
358 ReactionChannelName =
"kCC";
360 ReactionChannelName =
"kNC";
362 std::cout <<
"Current mode unknown!! " << std::endl;
365 ReactionChannelName +=
"_kQE";
367 ReactionChannelName +=
"_kRes";
369 ReactionChannelName +=
"_kDIS";
371 ReactionChannelName +=
"_kCoh";
373 ReactionChannelName +=
"_kNuElectronElastic";
375 ReactionChannelName +=
"_kInverseMuDecay";
377 std::cout <<
"interaction mode unknown!! " << std::endl;
379 return ReactionChannelName;
407 if (
id == -1) abort();
434 std::cout <<
"-----------> Particles in the Stack = " << mc.
NParticles() << std::endl;
435 std::cout << std::setiosflags(
std::ios::left) << std::setw(20) <<
"PARTICLE" 436 << std::setiosflags(
std::ios::left) << std::setw(32) <<
"STATUS" << std::setw(18)
437 <<
"E (GeV)" << std::setw(18) <<
"m (GeV/c2)" << std::setw(18) <<
"Ek (GeV)" 441 const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
447 if (
part.PdgCode() == 18040)
449 else if (
part.PdgCode() != -99999) {
450 name = databasePDG->GetParticle(
part.PdgCode())->GetName();
453 int code =
part.StatusCode();
455 double mass =
part.Mass();
457 double Ek = (energy - mass);
459 std::cout << std::setiosflags(
std::ios::left) << std::setw(20) << name
460 << std::setiosflags(
std::ios::left) << std::setw(32) << status << std::setw(18)
461 << energy << std::setw(18) << mass << std::setw(18) << Ek << std::endl;
470 if (
abs(
part.PdgCode()) == 11) {
478 else if (
abs(
part.PdgCode()) == 13) {
double E(const int i=0) const
NDKGen(fhicl::ParameterSet const &pset)
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
base_engine_t & createEngine(seed_t seed)
TH1F * fMuMomentum
momentum of outgoing muons
TH2F * fVertexXY
vertex location in xy
const simb::MCNeutrino & GetNeutrino() const
TH1F * fGenerated[6]
Spectra as generated.
double Py(const int i=0) const
TH1F * fMuDCosY
direction cosine of outgoing mu in y
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
TH1F * fECons
histogram to determine if energy is conserved in the event
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
void SetOrigin(simb::Origin_t origin)
TH1F * fEDCosZ
direction cosine of outgoing e in z
const simb::MCParticle & Nu() const
TH1F * fEMomentum
momentum of outgoing electrons
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.
TH2F * fVertexXZ
vertex location in xz
Geometry information for a single TPC.
double Px(const int i=0) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
TH1F * fVertexY
vertex location of generated events in y
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
TH1F * fCCMode
CC interaction mode.
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
TH1F * fEDCosY
direction cosine of outgoing e in y
void produce(art::Event &evt) override
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
A module to check the results from the Monte Carlo generator.
#define DEFINE_ART_MODULE(klass)
double P(const int i=0) const
TH1F * fMuDCosX
direction cosine of outgoing mu in x
TH2F * fVertexYZ
vertex location in yz
TH1F * fVertexZ
vertex location of generated events in z
TH1F * fDCosZ
direction cosine in z
std::string ReactionChannel(int ccnc, int mode)
std::string fNDKModuleLabel
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
const simb::MCParticle & GetParticle(int i) const
TH1F * fEDCosX
direction cosine of outgoing e in x
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)
void Add(simb::MCParticle const &part)
void FillHistograms(simb::MCTruth const &mc)
double Pz(const int i=0) const
double Vz(const int i=0) const
EventNumber_t event() const
std::string ParticleStatus(int StatusCode)
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.
void beginRun(art::Run &run) override
Event generator information.
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosY
direction cosine in y
Event generator information.
Namespace collecting geometry-related classes utilities.
Event Generation using GENIE, cosmics or single particles.
TH1F * fVertexX
vertex location of generated events in x
TH1F * fDCosX
direction cosine in x
double Vy(const int i=0) const
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
The data type to uniquely identify a cryostat.