22 #include "TDatabasePDG.h"
23 #include "TStopwatch.h"
25 #include "CLHEP/Random/RandFlat.h"
28 #include "art/Framework/Core/EDProducer.h"
29 #include "art/Framework/Core/ModuleMacros.h"
30 #include "art/Framework/Principal/Event.h"
31 #include "art/Framework/Principal/Run.h"
32 #include "art/Framework/Services/Registry/ServiceHandle.h"
33 #include "art_root_io/TFileService.h"
34 #include "fhiclcpp/ParameterSet.h"
37 #include "nurandom/RandomUtils/NuRandomService.h"
38 #include "nusimdata/SimulationBase/MCTruth.h"
39 #include "nusimdata/SimulationBase/MCParticle.h"
40 #include "nusimdata/SimulationBase/MCNeutrino.h"
49 class NDKGen :
public art::EDProducer {
51 explicit NDKGen(fhicl::ParameterSet
const &pset);
58 void beginRun(art::Run& run)
override;
110 , fNdkFile{pset.get<std::string>(
"NdkFile")}
111 , fEventFile{fNdkFile}
114 ,
fEngine(art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*
this, pset,
"Seed"))
118 produces< std::vector<simb::MCTruth> >();
121 if(!fEventFile.good()) {
122 throw cet::exception(
"NDKGen")
123 <<
"Could not open file: " << fNdkFile <<
'\n';
138 art::ServiceHandle<art::TFileService const>
tfs;
140 fGenerated[0] = tfs->make<TH1F>(
"fGenerated_necc",
"", 100, 0.0, 20.0);
141 fGenerated[1] = tfs->make<TH1F>(
"fGenerated_nebcc",
"", 100, 0.0, 20.0);
142 fGenerated[2] = tfs->make<TH1F>(
"fGenerated_nmcc",
"", 100, 0.0, 20.0);
143 fGenerated[3] = tfs->make<TH1F>(
"fGenerated_nmbcc",
"", 100, 0.0, 20.0);
144 fGenerated[4] = tfs->make<TH1F>(
"fGenerated_nnc",
"", 100, 0.0, 20.0);
145 fGenerated[5] = tfs->make<TH1F>(
"fGenerated_nbnc",
"", 100, 0.0, 20.0);
147 fDCosX = tfs->make<TH1F>(
"fDCosX",
";dx/ds", 200, -1., 1.);
148 fDCosY = tfs->make<TH1F>(
"fDCosY",
";dy/ds", 200, -1., 1.);
149 fDCosZ = tfs->make<TH1F>(
"fDCosZ",
";dz/ds", 200, -1., 1.);
151 fMuMomentum = tfs->make<TH1F>(
"fMuMomentum",
";p_{#mu} (GeV/c)", 500, 0., 50.);
152 fMuDCosX = tfs->make<TH1F>(
"fMuDCosX",
";dx/ds;", 200, -1., 1.);
153 fMuDCosY = tfs->make<TH1F>(
"fMuDCosY",
";dy/ds;", 200, -1., 1.);
154 fMuDCosZ = tfs->make<TH1F>(
"fMuDCosZ",
";dz/ds;", 200, -1., 1.);
156 fEMomentum = tfs->make<TH1F>(
"fEMomentum",
";p_{e} (GeV/c)", 500, 0., 50.);
157 fEDCosX = tfs->make<TH1F>(
"fEDCosX",
";dx/ds;", 200, -1., 1.);
158 fEDCosY = tfs->make<TH1F>(
"fEDCosY",
";dy/ds;", 200, -1., 1.);
159 fEDCosZ = tfs->make<TH1F>(
"fEDCosZ",
";dz/ds;", 200, -1., 1.);
161 fCCMode = tfs->make<TH1F>(
"fCCMode",
";CC Interaction Mode;", 5, 0., 5.);
162 fCCMode->GetXaxis()->SetBinLabel(1,
"QE");
163 fCCMode->GetXaxis()->SetBinLabel(2,
"Res");
164 fCCMode->GetXaxis()->SetBinLabel(3,
"DIS");
165 fCCMode->GetXaxis()->SetBinLabel(4,
"Coh");
166 fCCMode->GetXaxis()->SetBinLabel(5,
"kInverseMuDecay");
167 fCCMode->GetXaxis()->CenterLabels();
169 fNCMode = tfs->make<TH1F>(
"fNCMode",
";NC Interaction Mode;", 5, 0., 5.);
170 fNCMode->GetXaxis()->SetBinLabel(1,
"QE");
171 fNCMode->GetXaxis()->SetBinLabel(2,
"Res");
172 fNCMode->GetXaxis()->SetBinLabel(3,
"DIS");
173 fNCMode->GetXaxis()->SetBinLabel(4,
"Coh");
174 fNCMode->GetXaxis()->SetBinLabel(5,
"kNuElectronElastic");
175 fNCMode->GetXaxis()->CenterLabels();
177 fECons = tfs->make<TH1F>(
"fECons",
";#Delta E(#nu,lepton);", 500, -5., 5.);
179 art::ServiceHandle<geo::Geometry const> geo;
180 double x = 2.1*geo->DetHalfWidth();
181 double y = 2.1*geo->DetHalfHeight();
182 double z = 2.*geo->DetLength();
183 int xdiv = TMath::Nint(2*x/5.);
184 int ydiv = TMath::Nint(2*y/5.);
185 int zdiv = TMath::Nint(2*z/5.);
187 fVertexX = tfs->make<TH1F>(
"fVertexX",
";x (cm)", xdiv, -
x,
x);
188 fVertexY = tfs->make<TH1F>(
"fVertexY",
";y (cm)", ydiv, -
y,
y);
189 fVertexZ = tfs->make<TH1F>(
"fVertexZ",
";z (cm)", zdiv, -0.2*
z,
z);
191 fVertexXY = tfs->make<TH2F>(
"fVertexXY",
";x (cm);y (cm)", xdiv, -
x,
x, ydiv, -
y,
y);
192 fVertexXZ = tfs->make<TH2F>(
"fVertexXZ",
";z (cm);x (cm)", zdiv, -0.2*
z,
z, xdiv, -
x,
x);
193 fVertexYZ = tfs->make<TH2F>(
"fVertexYZ",
";z (cm);y (cm)", zdiv, -0.2*
z,
z, ydiv, -
y,
y);
199 art::ServiceHandle<geo::Geometry const> geo;
200 run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
214 std::cout<<
"------------------------------------------------------------------------------"<<std::endl;
218 std::cout <<
"event : " << evt.id().event() << std::endl;
236 std::string
name,
k, dollar;
243 int Idx, Ist,
PDG, Mother1, Mother2, Daughter1 ,Daughter2;
244 double Px, Py, Pz,
E,
m ;
245 std::string
p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14;
249 std::string primary(
"primary");
250 int FirstMother = -1;
262 TLorentzVector Neutrino;
263 TLorentzVector Lepton;
264 TLorentzVector Target;
266 TLorentzVector Hadron4mom;
279 std::unique_ptr< std::vector<simb::MCTruth> > truthcol(
new std::vector<simb::MCTruth>);
282 art::ServiceHandle<geo::Geometry const> geo;
285 double const fvCut{5.0};
294 for (
size_t i = 0; i<geo->NTPC(); ++i)
296 double local[3] = {0.,0.,0.};
297 double world[3] = {0.,0.,0.};
300 if (minx>world[0]-geo->DetHalfWidth(i))
301 minx = world[0]-geo->DetHalfWidth(i);
302 if (maxx<world[0]+geo->DetHalfWidth(i))
303 maxx = world[0]+geo->DetHalfWidth(i);
304 if (miny>world[1]-geo->DetHalfHeight(i))
305 miny = world[1]-geo->DetHalfHeight(i);
306 if (maxy<world[1]+geo->DetHalfHeight(i))
307 maxy = world[1]+geo->DetHalfHeight(i);
308 if (minz>world[2]-geo->DetLength(i)/2.)
309 minz = world[2]-geo->DetLength(i)/2.;
310 if (maxz<world[2]+geo->DetLength(i)/2.)
311 maxz = world[2]+geo->DetLength(i)/2.;
315 double X0 = 0.0 + flat.fire( minx+fvCut , maxx-fvCut );
316 double Y0 = 0.0 + flat.fire( miny+fvCut , maxy-fvCut );
317 double Z0 = 0.0 + flat.fire( minz+fvCut , maxz-fvCut );
319 std::cout <<
"NDKGen_module: X, Y, Z of vtx: " << X0 <<
", "<< Y0 <<
", "<< Z0 << std::endl;
324 std::cout <<
"NdkFile: Problem reading Ndk file" << std::endl;
328 if (k.find(
"** Event:")!= std::string::npos) {
329 std::istringstream
in;
333 in>> dummy>> dummy>> dummy >> dummy>> dummy>> dummy>> dummy >> dummy>> dummy>> dummy >> GenieEvt;
334 std::cout<<
"Genie Evt "<< GenieEvt <<
" art evt "<<evt.id().event()<<
"\n";
337 if (GenieEvt+1 != static_cast<int>(evt.id().event()))
341 if (!k.compare(0,25,
"GENIE Interaction Summary"))
343 if (k.compare(0,1,
"|") || k.compare(1,2,
" "))
continue;
344 if (k.find(
"Fin-Init") != std::string::npos)
continue;
345 if (k.find(
"Ar") != std::string::npos)
continue;
346 if (k.find(
"Cl") != std::string::npos)
continue;
347 if (k.find(
"HadrBlob") != std::string::npos)
continue;
348 if (k.find(
"NucBindE") != std::string::npos)
continue;
349 if (k.find(
"FLAGS") != std::string::npos)
break;
350 if (k.find(
"Vertex") != std::string::npos)
break;
353 if (!k.compare(26,1,
"1"))
356 std::istringstream
in;
360 in>>p1>> Idx >>p2>> Name >>p3>> Ist >>p4>> PDG >>p5>>Mother1 >> p6 >> Mother2 >>p7>> Daughter1 >>p8>> Daughter2 >>p9>>Px>>p10>>Py>>p11>>Pz>>p12>>E>>p13>> m>>p14;
362 if (Ist!=1)
continue;
364 std::cout <<
"PDG = " << PDG << std::endl;
366 const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
367 const TParticlePDG* definition = databasePDG->GetParticle(PDG);
368 Mass = definition->Mass();
369 if (E-Mass < 0.001)
continue;
374 simb::MCParticle mcpart(trackid,
382 P = std::sqrt(pow(E,2.) - pow(Mass,2.));
383 std::cout <<
"Momentum = " << P << std::endl;
385 TLorentzVector pos(X0, Y0, Z0, 0);
389 TLorentzVector mom(Px,Py,Pz, E);
391 mcpart.AddTrajectoryPoint(pos,mom);
396 truth.SetOrigin(simb::kUnknown);
405 std::cout <<
"NDKGen.cxx: Putting " << truth.NParticles() <<
" tracks on stack." << std::endl;
406 truthcol->push_back(truth);
408 evt.put(std::move(truthcol));
416 int code = StatusCode;
422 ParticleStatusName =
"kIStInitialState";
425 ParticleStatusName =
"kIStFinalState";
428 ParticleStatusName =
"kIStNucleonTarget";
431 ParticleStatusName =
"Status Unknown";
440 std::string ReactionChannelName=
" ";
443 ReactionChannelName =
"kCC";
445 ReactionChannelName =
"kNC";
446 else std::cout<<
"Current mode unknown!! "<<std::endl;
449 ReactionChannelName +=
"_kQE";
451 ReactionChannelName +=
"_kRes";
453 ReactionChannelName +=
"_kDIS";
455 ReactionChannelName +=
"_kCoh";
457 ReactionChannelName +=
"_kNuElectronElastic";
459 ReactionChannelName +=
"_kInverseMuDecay";
460 else std::cout<<
"interaction mode unknown!! "<<std::endl;
462 return ReactionChannelName;
470 if (mc.GetNeutrino().CCNC()==
simb::kCC) {
471 fCCMode->Fill(mc.GetNeutrino().Mode());
472 if (mc.GetNeutrino().Nu().PdgCode() == 12)
id = 0;
473 else if (mc.GetNeutrino().Nu().PdgCode() == -12)
id = 1;
474 else if (mc.GetNeutrino().Nu().PdgCode() == 14)
id = 2;
475 else if (mc.GetNeutrino().Nu().PdgCode() == -14)
id = 3;
479 fNCMode->Fill(mc.GetNeutrino().Mode());
480 if (mc.GetNeutrino().Nu().PdgCode() > 0)
id = 4;
486 fGenerated[id]->Fill(mc.GetNeutrino().Nu().E() );
490 simb::MCNeutrino mcnu = mc.GetNeutrino();
491 const simb::MCParticle
nu = mcnu.Nu();
503 fDCosX->Fill(nu.Px()/mom);
504 fDCosY->Fill(nu.Py()/mom);
505 fDCosZ->Fill(nu.Pz()/mom);
549 std::cout <<
"-----------> Particles in the Stack = " << mc.NParticles() << std::endl;
551 << std::setw(20) <<
"PARTICLE"
553 << std::setw(32) <<
"STATUS"
554 << std::setw(18) <<
"E (GeV)"
555 << std::setw(18) <<
"m (GeV/c2)"
556 << std::setw(18) <<
"Ek (GeV)"
557 << std::endl << std::endl;
559 const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
562 for(
int i = 0; i < mc.NParticles(); ++i){
563 simb::MCParticle part(mc.GetParticle(i));
565 if (part.PdgCode() == 18040)
567 else if (part.PdgCode() != -99999 )
569 name = databasePDG->GetParticle(part.PdgCode())->GetName();
572 int code = part.StatusCode();
574 double mass = part.Mass();
575 double energy = part.E();
576 double Ek = (energy-
mass);
580 << std::setw(18)<< energy
581 << std::setw(18)<< mass
582 << std::setw(18)<< Ek <<std::endl;
585 if(mc.GetNeutrino().CCNC() ==
simb::kCC){
589 for(
int i = 0; i < mc.NParticles(); ++i){
590 simb::MCParticle part(mc.GetParticle(i));
591 if(
abs(part.PdgCode()) == 11){
593 fEDCosX->Fill(part.Px()/part.P());
594 fEDCosY->Fill(part.Py()/part.P());
595 fEDCosZ->Fill(part.Pz()/part.P());
596 fECons->Fill(nu.E() - part.E());
599 else if(
abs(part.PdgCode()) == 13){
604 fECons->Fill(nu.E() - part.E());
NDKGen(fhicl::ParameterSet const &pset)
TH1F * fMuMomentum
momentum of outgoing muons
process_name opflash particleana ie ie ie z
TH2F * fVertexXY
vertex location in xy
process_name stream1 can override from command line with o or output services DetectorPropertiesService services DetectorPropertiesService services DetectorPropertiesService services DetectorPropertiesService physics analyzers pmtresponse NeutronTrackingCut services LArG4Parameters gaussian physics producers generator PDG
TH1F * fGenerated[6]
Spectra as generated.
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
process_name opflash particleana ie x
TH1F * fEDCosZ
direction cosine of outgoing e in z
process_name opdaq physics producers generator physics producers generator physics producers generator Z0
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexXZ
vertex location in xz
Geometry information for a single TPC.
TH1F * fVertexY
vertex location of generated events in y
produces< sumdata::RunData, art::InRun >()
process_name opdaq physics producers generator physics producers generator Y0
TH1F * fCCMode
CC interaction mode.
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
TH1F * fEDCosY
direction cosine of outgoing e in y
void produce(art::Event &evt) override
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
A module to check the results from the Monte Carlo generator.
Charged-current interactions.
standard_singlep gaussian distribution X0
process_name opflash particleana ie ie y
TH1F * fMuDCosX
direction cosine of outgoing mu in x
fEngine(art::ServiceHandle< rndm::NuRandomService >() ->createEngine(*this, pset,"Seed"))
TH2F * fVertexYZ
vertex location in yz
TH1F * fVertexZ
vertex location of generated events in z
TH1F * fDCosZ
direction cosine in z
BEGIN_PROLOG vertical distance to the surface Name
std::string ReactionChannel(int ccnc, int mode)
std::string fNDKModuleLabel
keep track of how long it takes to run the job
TH1F * fEDCosX
direction cosine of outgoing e in x
if &&[-z"$BASH_VERSION"] then echo Attempting to switch to bash bash shellSwitch exit fi &&["$1"= 'shellSwitch'] shift declare a IncludeDirectives for Dir in
void FillHistograms(simb::MCTruth const &mc)
std::string ParticleStatus(int StatusCode)
art::ServiceHandle< art::TFileService > tfs
void beginRun(art::Run &run) override
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosY
direction cosine in y
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
TH1F * fVertexX
vertex location of generated events in x
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
TH1F * fDCosX
direction cosine in x
physics associatedGroupsWithLeft p1
art framework interface to geometry description
BEGIN_PROLOG could also be cout