All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NDKGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 //
4 // NDK neutrino event generator
5 //
6 // echurch@fnal.gov
7 //
8 ////////////////////////////////////////////////////////////////////////
9 #include <cstdlib>
10 #include <fstream>
11 #include <iomanip>
12 #include <iostream>
13 #include <memory>
14 #include <sstream>
15 #include <stdio.h>
16 #include <string>
17 #include <vector>
18 
19 // ROOT includes
20 #include "TH1.h"
21 #include "TH2.h"
22 #include "TDatabasePDG.h"
23 #include "TStopwatch.h"
24 
25 #include "CLHEP/Random/RandFlat.h"
26 
27 // Framework includes
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"
35 
36 // art extensions
37 #include "nurandom/RandomUtils/NuRandomService.h"
38 #include "nusimdata/SimulationBase/MCTruth.h"
39 #include "nusimdata/SimulationBase/MCParticle.h"
40 #include "nusimdata/SimulationBase/MCNeutrino.h"
41 
42 // LArSoft includes
45 
46 
47 namespace evgen {
48  /// A module to check the results from the Monte Carlo generator
49  class NDKGen : public art::EDProducer {
50  public:
51  explicit NDKGen(fhicl::ParameterSet const &pset);
52  virtual ~NDKGen();
53 
54  private:
55 
56  void produce(art::Event& evt) override;
57  void beginJob() override;
58  void beginRun(art::Run& run) override;
59  void endJob() override;
60 
61  std::string ParticleStatus(int StatusCode);
62  std::string ReactionChannel(int ccnc,int mode);
63 
64  void FillHistograms(simb::MCTruth const& mc);
65 
66  std::string fNdkFile;
67  std::ifstream fEventFile;
68  TStopwatch fStopwatch; ///keep track of how long it takes to run the job
69 
70  std::string fNDKModuleLabel;
71  CLHEP::HepRandomEngine& fEngine; ///< art-managed random-number engine
72 
73  TH1F* fGenerated[6]; ///< Spectra as generated
74 
75  TH1F* fVertexX; ///< vertex location of generated events in x
76  TH1F* fVertexY; ///< vertex location of generated events in y
77  TH1F* fVertexZ; ///< vertex location of generated events in z
78 
79  TH2F* fVertexXY; ///< vertex location in xy
80  TH2F* fVertexXZ; ///< vertex location in xz
81  TH2F* fVertexYZ; ///< vertex location in yz
82 
83  TH1F* fDCosX; ///< direction cosine in x
84  TH1F* fDCosY; ///< direction cosine in y
85  TH1F* fDCosZ; ///< direction cosine in z
86 
87  TH1F* fMuMomentum; ///< momentum of outgoing muons
88  TH1F* fMuDCosX; ///< direction cosine of outgoing mu in x
89  TH1F* fMuDCosY; ///< direction cosine of outgoing mu in y
90  TH1F* fMuDCosZ; ///< direction cosine of outgoing mu in z
91 
92  TH1F* fEMomentum; ///< momentum of outgoing electrons
93  TH1F* fEDCosX; ///< direction cosine of outgoing e in x
94  TH1F* fEDCosY; ///< direction cosine of outgoing e in y
95  TH1F* fEDCosZ; ///< direction cosine of outgoing e in z
96 
97  TH1F* fCCMode; ///< CC interaction mode
98  TH1F* fNCMode; ///< CC interaction mode
99 
100  TH1F* fECons; ///< histogram to determine if energy is conserved in the event
101 
102  };
103 } // namespace
104 
105 namespace evgen{
106 
107  //____________________________________________________________________________
108  NDKGen::NDKGen(fhicl::ParameterSet const& pset)
109  : EDProducer{pset}
110  , fNdkFile{pset.get<std::string>("NdkFile")}
111  , fEventFile{fNdkFile}
112  // create a default random engine; obtain the random seed from NuRandomService,
113  // unless overridden in configuration with key "Seed"
114  , fEngine(art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*this, pset, "Seed"))
115  {
116  fStopwatch.Start();
117 
118  produces< std::vector<simb::MCTruth> >();
120 
121  if(!fEventFile.good()) {
122  throw cet::exception("NDKGen")
123  << "Could not open file: " << fNdkFile << '\n';
124  }
125 
126  }
127 
128  //____________________________________________________________________________
130  {
131  fStopwatch.Stop();
132  }
133 
134 //___________________________________________________________________________
135 
137  {
138  art::ServiceHandle<art::TFileService const> tfs;
139 
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);
146 
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.);
150 
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.);
155 
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.);
160 
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();
168 
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();
176 
177  fECons = tfs->make<TH1F>("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.);
178 
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.);
186 
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);
190 
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);
194  }
195 
196  //____________________________________________________________________________
197  void NDKGen::beginRun(art::Run& run)
198  {
199  art::ServiceHandle<geo::Geometry const> geo;
200  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
201  }
202 
203  //____________________________________________________________________________
205  {
206  fEventFile.close();
207  }
208 
209  //____________________________________________________________________________
210  void NDKGen::produce(art::Event& evt)
211  {
212 
213  std::cout << std::endl;
214  std::cout<<"------------------------------------------------------------------------------"<<std::endl;
215  //std::cout << "run : " << evt.Header().Run() << std::endl;
216  //std::cout << "subrun : " << evt.Header().Subrun() << std::endl;
217  //std::cout << "event : " << evt.Header().Event() << std::endl;
218  std::cout << "event : " << evt.id().event() << std::endl;
219 
220  // TODO: fill more quantities out, as below.
221  /*
222  double X; // vertex position from Ndk
223  double Y; // vertex position from Ndk
224  double Z; // vertex position from Ndk
225  double PDGCODE = -9999.;
226  double CHANNEL = -9999.;
227  int channel = -9999;
228  double energy = 0.; // in MeV from Ndk
229  double cosx = 0.;
230  double cosy = 0.;
231  double cosz = 0.;
232 
233  int partnumber = 0;
234  */
235 
236  std::string name, k, dollar;
237 
238 
239  // event dump format on file output by the two commands ....
240  // gevgen_ndcy -g 1000180400 -m 8 -n 400 -o ndk
241  // gevdump -f ndk.1000.ghep.root > ndk.out
242  std::string Name;
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;
246 
247 
248  int trackid = -1; // set track id to -i as these are all primary particles and have id <= 0
249  std::string primary("primary");
250  int FirstMother = -1;
251  double Mass = -9999;
252  int Status = -9999;
253 
254  double P; // momentum of MCParticle IN GeV/c
255 
256  // TODO: Could perhaps imagine using these in NDk.
257  /*
258  int targetnucleusPdg = -9999;
259  int hitquarkPdg = -9999;
260  double Q2 = -9999;
261  */
262  TLorentzVector Neutrino;
263  TLorentzVector Lepton;
264  TLorentzVector Target;
265  TLorentzVector q;
266  TLorentzVector Hadron4mom;
267 
268  // TODO: Could perhaps imagine using these in NDk.
269  /*
270  int Tpdg = 0; // for target
271  double Tmass = 0;
272  int Tstatus = 11;
273  double Tcosx, Tcosy, Tcosz, Tenergy;
274  */
275 
276  TLorentzVector Tpos;
277 
278 
279  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
280  simb::MCTruth truth;
281 
282  art::ServiceHandle<geo::Geometry const> geo;
283  CLHEP::RandFlat flat(fEngine);
284 
285  double const fvCut{5.0}; // force vtx to be this far from any wall.
286 
287  // Find boundary of active volume
288  double minx = 1e9;
289  double maxx = -1e9;
290  double miny = 1e9;
291  double maxy = -1e9;
292  double minz = 1e9;
293  double maxz = -1e9;
294  for (size_t i = 0; i<geo->NTPC(); ++i)
295  {
296  double local[3] = {0.,0.,0.};
297  double world[3] = {0.,0.,0.};
298  const geo::TPCGeo &tpc = geo->TPC(i);
299  tpc.LocalToWorld(local,world);
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.;
312  }
313 
314  // Assign vertice position
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 );
318 
319  std::cout << "NDKGen_module: X, Y, Z of vtx: " << X0 << ", "<< Y0 << ", "<< Z0 << std::endl;
320 
321  int GenieEvt = -999;
322 
323  if(!fEventFile.good())
324  std::cout << "NdkFile: Problem reading Ndk file" << std::endl;
325 
326  while(getline(fEventFile,k)){
327 
328  if (k.find("** Event:")!= std::string::npos) {
329  std::istringstream in;
330  in.clear();
331  in.str(k);
332  std::string dummy;
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";
335  }
336 
337  if (GenieEvt+1 != static_cast<int>(evt.id().event()))
338  continue;
339  else {
340 
341  if (!k.compare(0,25,"GENIE Interaction Summary")) // testing for new event.
342  break;
343  if (k.compare(0,1,"|") || k.compare(1,2," ")) continue; // uninteresting line if it doesn't start with "|" and if second and third characters aren't spaces.
344  if (k.find("Fin-Init") != std::string::npos) continue; // Meh.
345  if (k.find("Ar") != std::string::npos) continue; // Meh.
346  if (k.find("Cl") != std::string::npos) continue; // ignore chlorine nucleus in nnbar events
347  if (k.find("HadrBlob") != std::string::npos) continue; // Meh.
348  if (k.find("NucBindE") != std::string::npos) continue; // Meh. atmo
349  if (k.find("FLAGS") != std::string::npos) break; // Event end. gevgen_ndcy
350  if (k.find("Vertex") != std::string::npos) break; // Event end. atmo
351 
352  // if (!k.compare(26,1,"3") || !k.compare(26,1,"1")) ; // New event or stable particles.
353  if (!k.compare(26,1,"1")) // New event or stable particles.
354  {
355 
356  std::istringstream in;
357  in.clear();
358  in.str(k);
359 
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;
361  //std::cout<<std::setprecision(9)<<dollar<<" "<<name<<" "<<PDGCODE<<" "<<energy<<" "<<cosx<<" "<<cosy<<" "<<cosz<<" "<<partnumber<<std::endl;
362  if (Ist!=1) continue;
363 
364  std::cout << "PDG = " << PDG << std::endl;
365 
366  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
367  const TParticlePDG* definition = databasePDG->GetParticle(PDG);
368  Mass = definition->Mass(); // GeV
369  if (E-Mass < 0.001) continue; // KE is too low.
370 
371  // if(partnumber == -1)
372  Status = 1;
373 
374  simb::MCParticle mcpart(trackid,
375  PDG,
376  primary,
377  FirstMother,
378  Mass,
379  Status
380  );
381 
382  P = std::sqrt(pow(E,2.) - pow(Mass,2.)); // GeV/c
383  std::cout << "Momentum = " << P << std::endl;
384 
385  TLorentzVector pos(X0, Y0, Z0, 0);
386 
387  Tpos = pos; // for target
388 
389  TLorentzVector mom(Px,Py,Pz, E);
390 
391  mcpart.AddTrajectoryPoint(pos,mom);
392  truth.Add(mcpart);
393 
394 
395  }// loop over particles in an event
396  truth.SetOrigin(simb::kUnknown);
397 
398  //if (!k.compare(1,1,"FLAGS")) // end of event
399  // break;
400 
401  }
402  } // end while loop
403 
404  /////////////////////////////////
405  std::cout << "NDKGen.cxx: Putting " << truth.NParticles() << " tracks on stack." << std::endl;
406  truthcol->push_back(truth);
407  //FillHistograms(truth);
408  evt.put(std::move(truthcol));
409 
410  return;
411  }
412 
413 // //......................................................................
414  std::string NDKGen::ParticleStatus(int StatusCode)
415  {
416  int code = StatusCode;
417  std::string ParticleStatusName;
418 
419  switch(code)
420  {
421  case 0:
422  ParticleStatusName = "kIStInitialState";
423  break;
424  case 1:
425  ParticleStatusName = "kIStFinalState";
426  break;
427  case 11:
428  ParticleStatusName = "kIStNucleonTarget";
429  break;
430  default:
431  ParticleStatusName = "Status Unknown";
432  }
433  return ParticleStatusName;
434  }
435 
436 
437 // //......................................................................
438  std::string NDKGen::ReactionChannel(int ccnc,int mode)
439  {
440  std::string ReactionChannelName=" ";
441 
442  if(ccnc==0)
443  ReactionChannelName = "kCC";
444  else if(ccnc==1)
445  ReactionChannelName = "kNC";
446  else std::cout<<"Current mode unknown!! "<<std::endl;
447 
448  if(mode==0)
449  ReactionChannelName += "_kQE";
450  else if(mode==1)
451  ReactionChannelName += "_kRes";
452  else if(mode==2)
453  ReactionChannelName += "_kDIS";
454  else if(mode==3)
455  ReactionChannelName += "_kCoh";
456  else if(mode==4)
457  ReactionChannelName += "_kNuElectronElastic";
458  else if(mode==5)
459  ReactionChannelName += "_kInverseMuDecay";
460  else std::cout<<"interaction mode unknown!! "<<std::endl;
461 
462  return ReactionChannelName;
463  }
464 
465 // //......................................................................
466  void NDKGen::FillHistograms(simb::MCTruth const& mc)
467  {
468  // Decide which histograms to put the spectrum in
469  int id = -1;
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;
476  else return;
477  }
478  else {
479  fNCMode->Fill(mc.GetNeutrino().Mode());
480  if (mc.GetNeutrino().Nu().PdgCode() > 0) id = 4;
481  else id = 5;
482  }
483  if (id==-1) abort();
484 
485  // Fill the specta histograms
486  fGenerated[id]->Fill(mc.GetNeutrino().Nu().E() );
487 
488  //< fill the vertex histograms from the neutrino - that is always
489  //< particle 0 in the list
490  simb::MCNeutrino mcnu = mc.GetNeutrino();
491  const simb::MCParticle nu = mcnu.Nu();
492 
493  fVertexX->Fill(nu.Vx());
494  fVertexY->Fill(nu.Vy());
495  fVertexZ->Fill(nu.Vz());
496 
497  fVertexXY->Fill(nu.Vx(), nu.Vy());
498  fVertexXZ->Fill(nu.Vz(), nu.Vx());
499  fVertexYZ->Fill(nu.Vz(), nu.Vy());
500 
501  double mom = nu.P();
502  if(std::abs(mom) > 0.){
503  fDCosX->Fill(nu.Px()/mom);
504  fDCosY->Fill(nu.Py()/mom);
505  fDCosZ->Fill(nu.Pz()/mom);
506  }
507 
508 
509 // MF_LOG_DEBUG("GENIEInteractionInformation")
510 // << std::endl
511 // << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode())
512 // << std::endl
513 // << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl
514 // << std::setiosflags(std::ios::left)
515 // << std::setw(20) << "PARTICLE"
516 // << std::setiosflags(std::ios::left)
517 // << std::setw(32) << "STATUS"
518 // << std::setw(18) << "E (GeV)"
519 // << std::setw(18) << "m (GeV/c2)"
520 // << std::setw(18) << "Ek (GeV)"
521 // << std::endl << std::endl;
522 
523 // const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
524 
525 // // Loop over the particle stack for this event
526 // for(int i = 0; i < mc.NParticles(); ++i){
527 // simb::MCParticle part(mc.GetParticle(i));
528 // std::string name = databasePDG->GetParticle(part.PdgCode())->GetName();
529 // int code = part.StatusCode();
530 // std::string status = ParticleStatus(code);
531 // double mass = part.Mass();
532 // double energy = part.E();
533 // double Ek = (energy-mass); // Kinetic Energy (GeV)
534 // if(status=="kIStFinalStB4Interactions")
535 // MF_LOG_DEBUG("GENIEFinalState")
536 // << std::setiosflags(std::ios::left) << std::setw(20) << name
537 // << std::setiosflags(std::ios::left) << std::setw(32) <<status
538 // << std::setw(18)<< energy
539 // << std::setw(18)<< mass
540 // << std::setw(18)<< Ek <<std::endl;
541 // else
542 // MF_LOG_DEBUG("GENIEFinalState")
543 // << std::setiosflags(std::ios::left) << std::setw(20) << name
544 // << std::setiosflags(std::ios::left) << std::setw(32) << status
545 // << std::setw(18) << energy
546 // << std::setw(18) << mass <<std::endl;
547 
548  std::cout << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode()) << std::endl;
549  std::cout << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl;
550  std::cout << std::setiosflags(std::ios::left)
551  << std::setw(20) << "PARTICLE"
552  << std::setiosflags(std::ios::left)
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;
558 
559  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
560 
561  // Loop over the particle stack for this event
562  for(int i = 0; i < mc.NParticles(); ++i){
563  simb::MCParticle part(mc.GetParticle(i));
564  std::string name;
565  if (part.PdgCode() == 18040)
566  name = "Ar40 18040";
567  else if (part.PdgCode() != -99999 )
568  {
569  name = databasePDG->GetParticle(part.PdgCode())->GetName();
570  }
571 
572  int code = part.StatusCode();
573  std::string status = ParticleStatus(code);
574  double mass = part.Mass();
575  double energy = part.E();
576  double Ek = (energy-mass); // Kinetic Energy (GeV)
577 
578  std::cout << std::setiosflags(std::ios::left) << std::setw(20) << name
579  << std::setiosflags(std::ios::left) << std::setw(32) <<status
580  << std::setw(18)<< energy
581  << std::setw(18)<< mass
582  << std::setw(18)<< Ek <<std::endl;
583  }
584 
585  if(mc.GetNeutrino().CCNC() == simb::kCC){
586 
587  ///look for the outgoing lepton in the particle stack
588  ///just interested in the first one
589  for(int i = 0; i < mc.NParticles(); ++i){
590  simb::MCParticle part(mc.GetParticle(i));
591  if(abs(part.PdgCode()) == 11){
592  fEMomentum->Fill(part.P());
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());
597  break;
598  }
599  else if(abs(part.PdgCode()) == 13){
600  fMuMomentum->Fill(part.P());
601  fMuDCosX->Fill(part.Px()/part.P());
602  fMuDCosY->Fill(part.Py()/part.P());
603  fMuDCosZ->Fill(part.Pz()/part.P());
604  fECons->Fill(nu.E() - part.E());
605  break;
606  }
607  }// end loop over particles
608  }//end if CC interaction
609 
610  return;
611  }
612 
613  DEFINE_ART_MODULE(NDKGen)
614 
615 } // namespace
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
Definition: gen_protons.fcl:45
std::string fNdkFile
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexXZ
vertex location in xz
Geometry information for a single TPC.
Definition: TPCGeo.h:38
TH1F * fVertexY
vertex location of generated events in y
produces< sumdata::RunData, art::InRun >()
process_name opdaq physics producers generator physics producers generator Y0
Definition: gen_protons.fcl:45
process_name E
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
TStopwatch fStopwatch
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
then local
T abs(T value)
A module to check the results from the Monte Carlo generator.
Charged-current interactions.
Definition: IPrediction.h:38
virtual ~NDKGen()
standard_singlep gaussian distribution X0
Definition: multigen.fcl:8
process_name opflash particleana ie ie y
const char mode
Definition: noise_ana.cxx:20
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
walls no left
Definition: selectors.fcl:105
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
void beginJob() override
std::ifstream fEventFile
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 endJob() override
float mass
Definition: dedx.py:47
void FillHistograms(simb::MCTruth const &mc)
then echo fcl name
std::string ParticleStatus(int StatusCode)
art::ServiceHandle< art::TFileService > tfs
void beginRun(art::Run &run) override
TCEvent evt
Definition: DataStructs.cxx:8
pdgs k
Definition: selectors.fcl:22
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosY
direction cosine in y
BEGIN_PROLOG SN nu
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.
Definition: TPCGeo.h:563
TH1F * fDCosX
direction cosine in x
physics associatedGroupsWithLeft p1
art framework interface to geometry description
BEGIN_PROLOG could also be cout