LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
NUANCEGen_module.cc
Go to the documentation of this file.
1 //
3 //
4 // NUANCE neutrino event generator
5 //
6 // brebel@fnal.gov
7 // saima@ksu.edu
8 //
10 #ifndef EVGEN_NUANCEGEN_H
11 #define EVGEN_NUANCEGEN_H
12 
13 #include <cstdlib>
14 #include <string>
15 #include <iostream>
16 #include <iomanip>
17 #include <sstream>
18 #include <vector>
19 #include <map>
20 #include <memory>
21 #include <unistd.h>
22 #include <stdio.h>
23 #include <fstream>
24 #include <memory> // std::unique_ptr<>
25 
26 // ROOT includes
27 #include "TH1.h"
28 #include "TH2.h"
29 #include "TDatabasePDG.h"
30 #include "TSystem.h"
31 #include "TStopwatch.h"
32 
33 // Framework includes
37 #include "fhiclcpp/ParameterSet.h"
44 
45 // LArSoft includes
50 
51 // class TH1F;
52 // class TH2F;
53 
54 // namespace simb { class MCTruth; }
55 
56 namespace evgen {
58  class NUANCEGen : public art::EDProducer {
59  public:
60  explicit NUANCEGen(fhicl::ParameterSet const &pset);
61  virtual ~NUANCEGen();
62 
63  void produce(art::Event& evt);
64  void beginJob();
65  void beginRun(art::Run& run);
66  void reconfigure(fhicl::ParameterSet const& p);
67  void endJob();
68 
69  private:
70 
71  std::string ParticleStatus(int StatusCode);
72  std::string ReactionChannel(int ccnc,int mode);
73 
75 
76  std::string fNuanceFile;
78  std::unique_ptr<std::ifstream> fEventFile;
79  TStopwatch fStopwatch;
80 
81 
82  std::string fNUANCEModuleLabel;
83 
84  TH1F* fGenerated[6];
85 
86  TH1F* fVertexX;
87  TH1F* fVertexY;
88  TH1F* fVertexZ;
89 
90  TH2F* fVertexXY;
91  TH2F* fVertexXZ;
92  TH2F* fVertexYZ;
93 
94  TH1F* fDCosX;
95  TH1F* fDCosY;
96  TH1F* fDCosZ;
97 
98  TH1F* fMuMomentum;
99  TH1F* fMuDCosX;
100  TH1F* fMuDCosY;
101  TH1F* fMuDCosZ;
102 
103  TH1F* fEMomentum;
104  TH1F* fEDCosX;
105  TH1F* fEDCosY;
106  TH1F* fEDCosZ;
107 
108  TH1F* fCCMode;
109  TH1F* fNCMode;
110 
111  // for c2: fDeltaE is no longer used
112  //TH1F* fDeltaE; ///< difference in neutrino energy from MCTruth::Enu() vs TParticle
113  TH1F* fECons;
114 
115  };
116 }
117 
118 namespace evgen{
119 
120  //____________________________________________________________________________
122  {
123  this->reconfigure(pset);
124  fStopwatch.Start();
125 
126  produces< std::vector<simb::MCTruth> >();
127  produces< sumdata::RunData, art::InRun >();
128 
129  mf::LogInfo("NUANCEGen") << "Reading events from '" << fNuanceFile << "'";
130  fEventFile = std::make_unique<std::ifstream>(fNuanceFile.c_str());
131  if (!(fEventFile->good())) {
133  << "Failed to open input file '" << fNuanceFile << "'.\n";
134  }
135  }
136 
137  //____________________________________________________________________________
139  {
140  fStopwatch.Stop();
141  }
142 
143  //____________________________________________________________________________
145  {
146  fNuanceFile =(p.get< std::string >("NuanceFile"));
147  fBeamVerticalAngle =(p.get< double >("BeamVerticalAngle"));
148  return;
149  }
150 //___________________________________________________________________________
151 
153 
154  // Get access to the TFile service.
156 
157  fGenerated[0] = tfs->make<TH1F>("fGenerated_necc","", 100, 0.0, 20.0);
158  fGenerated[1] = tfs->make<TH1F>("fGenerated_nebcc","", 100, 0.0, 20.0);
159  fGenerated[2] = tfs->make<TH1F>("fGenerated_nmcc","", 100, 0.0, 20.0);
160  fGenerated[3] = tfs->make<TH1F>("fGenerated_nmbcc","", 100, 0.0, 20.0);
161  fGenerated[4] = tfs->make<TH1F>("fGenerated_nnc","", 100, 0.0, 20.0);
162  fGenerated[5] = tfs->make<TH1F>("fGenerated_nbnc","", 100, 0.0, 20.0);
163 
164  fDCosX = tfs->make<TH1F>("fDCosX", ";dx/ds", 200, -1., 1.);
165  fDCosY = tfs->make<TH1F>("fDCosY", ";dy/ds", 200, -1., 1.);
166  fDCosZ = tfs->make<TH1F>("fDCosZ", ";dz/ds", 200, -1., 1.);
167 
168  fMuMomentum = tfs->make<TH1F>("fMuMomentum", ";p_{#mu} (GeV/c)", 500, 0., 50.);
169  fMuDCosX = tfs->make<TH1F>("fMuDCosX", ";dx/ds;", 200, -1., 1.);
170  fMuDCosY = tfs->make<TH1F>("fMuDCosY", ";dy/ds;", 200, -1., 1.);
171  fMuDCosZ = tfs->make<TH1F>("fMuDCosZ", ";dz/ds;", 200, -1., 1.);
172 
173  fEMomentum = tfs->make<TH1F>("fEMomentum", ";p_{e} (GeV/c)", 500, 0., 50.);
174  fEDCosX = tfs->make<TH1F>("fEDCosX", ";dx/ds;", 200, -1., 1.);
175  fEDCosY = tfs->make<TH1F>("fEDCosY", ";dy/ds;", 200, -1., 1.);
176  fEDCosZ = tfs->make<TH1F>("fEDCosZ", ";dz/ds;", 200, -1., 1.);
177 
178  fCCMode = tfs->make<TH1F>("fCCMode", ";CC Interaction Mode;", 5, 0., 5.);
179  fCCMode->GetXaxis()->SetBinLabel(1, "QE");
180  fCCMode->GetXaxis()->SetBinLabel(2, "Res");
181  fCCMode->GetXaxis()->SetBinLabel(3, "DIS");
182  fCCMode->GetXaxis()->SetBinLabel(4, "Coh");
183  fCCMode->GetXaxis()->SetBinLabel(5, "kInverseMuDecay");
184  fCCMode->GetXaxis()->CenterLabels();
185 
186  fNCMode = tfs->make<TH1F>("fNCMode", ";NC Interaction Mode;", 5, 0., 5.);
187  fNCMode->GetXaxis()->SetBinLabel(1, "QE");
188  fNCMode->GetXaxis()->SetBinLabel(2, "Res");
189  fNCMode->GetXaxis()->SetBinLabel(3, "DIS");
190  fNCMode->GetXaxis()->SetBinLabel(4, "Coh");
191  fNCMode->GetXaxis()->SetBinLabel(5, "kNuElectronElastic");
192  fNCMode->GetXaxis()->CenterLabels();
193 
194  //fDeltaE = tfs->make<TH1F>("fDeltaE", ";#Delta E_{#nu} (GeV);", 200, -1., 1.);
195  fECons = tfs->make<TH1F>("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.);
196 
198  double x = 2.1*geo->DetHalfWidth();
199  double y = 2.1*geo->DetHalfHeight();
200  double z = 2.*geo->DetLength();
201  int xdiv = TMath::Nint(2*x/5.);
202  int ydiv = TMath::Nint(2*y/5.);
203  int zdiv = TMath::Nint(2*z/5.);
204 
205  fVertexX = tfs->make<TH1F>("fVertexX", ";x (cm)", xdiv, -x, x);
206  fVertexY = tfs->make<TH1F>("fVertexY", ";y (cm)", ydiv, -y, y);
207  fVertexZ = tfs->make<TH1F>("fVertexZ", ";z (cm)", zdiv, -0.2*z, z);
208 
209  fVertexXY = tfs->make<TH2F>("fVertexXY", ";x (cm);y (cm)", xdiv, -x, x, ydiv, -y, y);
210  fVertexXZ = tfs->make<TH2F>("fVertexXZ", ";z (cm);x (cm)", zdiv, -0.2*z, z, xdiv, -x, x);
211  fVertexYZ = tfs->make<TH2F>("fVertexYZ", ";z (cm);y (cm)", zdiv, -0.2*z, z, ydiv, -y, y);
212 
213  }
214 
215  //____________________________________________________________________________
217  {
218 
219  // grab the geometry object to see what geometry we are using
221  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetectorName()));
222 
223  run.put(std::move(runcol));
224 
225  return;
226  }
227 
228  //____________________________________________________________________________
230  {
231  fEventFile->close();
232  }
233  //____________________________________________________________________________
235  {
236 
237  std::cout << std::endl;
238  std::cout<<"------------------------------------------------------------------------------"<<std::endl;
239 // std::cout << "run : " << evt.Header().Run() << std::endl;
240 // std::cout << "subrun : " << evt.Header().Subrun() << std::endl;
241 // std::cout << "event : " << evt.Header().Event() << std::endl;
242  std::cout << "event : " << evt.id().event() << std::endl;
243 
244  double X = 0.; // vertex position from Nuance
245  double Y = 0.; // vertex position from Nuance
246  double Z = 0.; // vertex position from Nuance
247  double PDGCODE = -9999.;
248  double CHANNEL = -9999.;
249  int channel = -9999;
250  double energy = 0.; // in MeV from Nuance
251  double cosx = 0.;
252  double cosy = 0.;
253  double cosz = 0.;
254  std::string name, k, dollar;
255  int partnumber = 0;
256 
257  int trackid = -1; // set track id to -i as these are all primary particles and have id <= 0
258  std::string primary("primary");
259  int FirstMother = -1;
260  double Mass = -9999;
261  int Status = -9999;
262 
263 
264  int ccnc = -9999;
265  int mode = -9999;
266  int targetnucleusPdg = -9999;
267  int hitquarkPdg = -9999;
268 
269  double P; // momentum of MCParticle in GeV/c
270  TLorentzVector Neutrino;
271  TLorentzVector Lepton;
272  TLorentzVector Target;
273  TLorentzVector q;
274  TLorentzVector Hadron4mom;
275 
276  double InvariantMass = -9999;
277  double x = -9999;
278  double y = -9999;
279  double Q2 = -9999;
280 
281  int Tpdg = 0; // for target
282  double Tmass = 0;
283  int Tstatus = 11;
284  double Tcosx = 0., Tcosy = 0., Tcosz = 0., Tenergy = 0.;
285  TLorentzVector Tpos;
286  double M = 0;
287 
288  // rotated direction cosines
289  double CosX = 0;
290  double CosY = 0;
291  double CosZ = 0;
292 
293 
294  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
295  simb::MCTruth truth;
296 
297  if(!fEventFile->good())
298  std::cout << "NuanceFile: Problem reading nuance file" << std::endl;
299 
300  while(getline(*fEventFile,k)){
301  std::istringstream in;
302  in.clear();
303  in.str(k);
304 
305  in>>dollar>>name>>PDGCODE>>energy>>cosx>>cosy>>cosz>>partnumber;
306  //std::cout<<std::setprecision(9)<<dollar<<" "<<name<<" "<<PDGCODE<<" "<<energy<<" "<<cosx<<" "<<cosy<<" "<<cosz<<" "<<partnumber<<std::endl;
307 
308  //get the nuance channel number
309  if(name == "nuance"){
310  CHANNEL = PDGCODE;
311  channel = int (CHANNEL);
312  channel = abs(channel)+ simb::kNuanceOffset;
313  //std::cout << "channel = " << channel << std::endl;
314 
315  //set the interaction type; CC or NC
316  if ( abs(channel) == simb::kCCQE
317  || ( abs(channel) >= simb::kResCCNuProtonPiPlus && abs(channel) <= simb::kResCCNuNeutronPiPlus )
318  || ( abs(channel) >= simb::kResCCNuBarNeutronPiMinus && abs(channel) <= simb::kResCCNuBarProtonPiMinus )
319  || ( abs(channel) >= simb::kResCCNuDeltaPlusPiPlus && abs(channel) <= simb::kResCCNuDelta2PlusPiMinus )
320  || ( abs(channel) >= simb::kResCCNuBarDelta0PiMinus && abs(channel) <= simb::kResCCNuBarDeltaMinusPiPlus )
321  || ( abs(channel) >= simb::kResCCNuProtonRhoPlus && abs(channel) <= simb::kResCCNuNeutronRhoPlus )
322  || ( abs(channel) >= simb::kResCCNuBarNeutronRhoMinus && abs(channel) <= simb::kResCCNuBarNeutronRho0 )
323  || ( abs(channel) >= simb::kResCCNuSigmaPlusKaonPlus && abs(channel) <= simb::kResCCNuSigmaPlusKaon0 )
324  || ( abs(channel) >= simb::kResCCNuBarSigmaMinusKaon0 && abs(channel) <= simb::kResCCNuBarSigma0Kaon0 )
325  || abs(channel) == simb::kResCCNuProtonEta
326  || abs(channel) == simb::kResCCNuBarNeutronEta
327  || abs(channel) == simb::kResCCNuKaonPlusLambda0
328  || abs(channel) == simb::kResCCNuBarKaon0Lambda0
329  || ( abs(channel) >= simb::kResCCNuProtonPiPlusPiMinus && abs(channel) <= simb::kResCCNuProtonPi0Pi0 )
330  || ( abs(channel) >= simb::kResCCNuBarNeutronPiPlusPiMinus && abs(channel) <= simb::kResCCNuBarNeutronPi0Pi0 )
331  || abs(channel) == simb::kCCDIS || abs(channel) == simb::kCCQEHyperon
332  || abs(channel) == simb::kCCCOH || abs(channel) == simb::kInverseMuDecay )
333  ccnc = simb::kCC;
334 
335  else if ( abs(channel) != simb::kUnUsed1 && abs(channel) != simb::kUnUsed2 )
336  ccnc = simb::kNC;
337 
338  //set the interaction mode; QE, Res, DIS, Coh, kNuElectronElastic, kInverseMuDecay
339  if ( abs(channel) == simb::kCCQE || abs(channel) == simb::kNCQE || abs(channel) == simb::kCCQEHyperon )
340  mode = simb::kQE;
341  else if ( abs(channel) >= simb::kResCCNuProtonPiPlus && abs(channel) <= simb::kResCCNuBarProtonPi0Pi0 )
342  mode = simb::kRes;
343  else if ( abs(channel) == simb::kCCDIS || abs(channel) == simb::kNCDIS )
344  mode = simb::kDIS;
345  else if ( abs(channel) == simb::kNCCOH || abs(channel) == simb::kCCCOH )
346  mode = simb::kCoh;
347  else if ( abs(channel) == simb::kNuElectronElastic )
348  mode = 4;
349  else if ( abs(channel) == simb::kInverseMuDecay )
350  mode = 5;
351 
352  } // if(name == "nuance")
353 
354 
355  //get the nuance vertex position:
356  if(name == "vertex"){
357  X = PDGCODE;
358  Y = energy;
359  Z = cosx;
360  //std::cout << "vertex from nuance = " << X << " " << Y << " " << Z << std::endl;
361  }
362 
363  double PI = 3.14159265;
364 
365  CosX = cosx;
366  CosY = (TMath::Cos(fBeamVerticalAngle*PI/180)*cosy)+(TMath::Sin(fBeamVerticalAngle*PI/180)*cosz);
367  CosZ = (-TMath::Sin(fBeamVerticalAngle*PI/180)*cosy)+(TMath::Cos(fBeamVerticalAngle*PI/180)*cosz);
368 
369  //get the target info
370  if(name == "track" && (PDGCODE == 2212 || PDGCODE == 2112 || PDGCODE == 18040 || PDGCODE == 11) && partnumber == -1){
371  Tpdg = int (PDGCODE);
372 
373  if ( PDGCODE == 18040 ){
374  Tmass = 37.5593438; // GeV
375  }
376 
377  else {
378  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
379  const TParticlePDG* definition = databasePDG->GetParticle(Tpdg);
380  Tmass = definition->Mass();
381  }
382 
383  Tenergy = energy;
384  Tcosx = CosX;
385  Tcosy = CosY;
386  Tcosz = CosZ;
387  }
388 
389  //get mcparticles other than target
390  if(name == "track" && (
391  partnumber == 0 ||
392  (partnumber == -1 && (PDGCODE != 2212 && PDGCODE != 2112 && PDGCODE != 18040 && PDGCODE != 11))
393  )
394  ){
395 
396  int pdgcode = int (PDGCODE);
397  //std::cout << "pdgcode = " << pdgcode << std::endl;
398 
399  if ( PDGCODE == 18040 )
400  Mass = 37.5593438; // GeV
401 
402  else {
403  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
404  const TParticlePDG* definition = databasePDG->GetParticle(pdgcode);
405  Mass = definition->Mass(); // GeV
406  }
407 
408  if(partnumber == -1)
409  Status = 0;
410 
411  if(partnumber == 0)
412  Status = 1;
413 
414  simb::MCParticle mcpart(trackid,
415  pdgcode,
416  primary,
417  FirstMother,
418  Mass,
419  Status
420  );
421 
422  P = std::sqrt(pow(energy/1000,2.) - pow(Mass,2.)); // GeV/c
423  //std::cout << "Momentum = " << P << std::endl;
424 
426 
427  double X0 = X + geo->DetHalfWidth();
428  double Y0 = Y;
429  double Z0 = Z + 0.5*geo->DetLength();
430 
431  TLorentzVector pos(X0, Y0, Z0, 0);
432  Tpos = pos; // for target
433 
434  TLorentzVector mom(CosX*P, CosY*P, CosZ*P, energy/1000);
435 
436  mcpart.AddTrajectoryPoint(pos,mom);
437  truth.Add(mcpart);
438 
439 
440  if(name == "track" && (abs(pdgcode) == 14 || abs(pdgcode) == 12) && partnumber == -1)
441  Neutrino.SetPxPyPzE(CosX*P, CosY*P, CosZ*P, energy/1000);
442 
443  if(name == "track" && (abs(pdgcode) == 13 || abs(pdgcode) == 11 || abs(pdgcode) == 14 || abs(pdgcode) == 12) && partnumber == 0)
444  Lepton.SetPxPyPzE(CosX*P, CosY*P, CosZ*P, energy/1000);
445 
446  }// loop over particles in an event. (excluding target)
447 
448  if(name == "end"){
449  break;
450  }
451 
452  } // end while loop
453 
454 
456 
457  //adding the target to mcparticle
458  simb::MCParticle mcpart(trackid,
459  Tpdg,
460  primary,
461  FirstMother,
462  Tmass,
463  Tstatus
464  );
465 
466  TLorentzVector Tmom;
467  Tmom.SetPxPyPzE(Tcosx, Tcosy, Tcosz, Tenergy/1000); // target momentum coordinates make a unit vector, may be modified later, does not effect the event simulation
468 
469  mcpart.AddTrajectoryPoint(Tpos,Tmom);
470  truth.Add(mcpart);
471 
472 // std::cout<<"Neutrino 4-Momentum = "
473 // <<Neutrino.Px()<<" "
474 // << Neutrino.Py()<<" "
475 // <<Neutrino.Pz()<<" "
476 // <<Neutrino.E()<<std::endl;
477 
478 // std::cout << "Target 4-Momentum = "
479 // << Target.Px() << " "
480 // << Target.Py() << " "
481 // << Target.Pz()<< " "
482 // << Target.E() << std::endl;
483 
484 // std::cout << "4-momentum of Lepton = "
485 // << Lepton.Px() << " "
486 // << Lepton.Py() << " "
487 // << Lepton.Pz() << " "
488 // << Lepton.E() << std::endl;
489 
490  M = Tmass;
491  q = Neutrino - Lepton;
492  Q2 = -(q*q);
493 
494  double v = q.E();
495  x = Q2/ (2*M*v);
496  y = v/ Neutrino.E();
497  double W2 = M*M + 2*M*v - Q2;
498  InvariantMass = TMath::Sqrt(TMath::Max(0.,W2));
499 
500  if(mode == simb::kCoh){
501  x = -1;
502  y = -1;
503  InvariantMass = -1;
504  }
505 
507  truth.SetNeutrino(ccnc, mode, channel,
508  targetnucleusPdg,
509  Tpdg,
510  hitquarkPdg,
511  //InvariantMass, x, y, Q2
512  InvariantMass, x, y, Q2
513  );
514 
515  std::cout << truth.GetNeutrino() << std::endl;
516 
517  truthcol->push_back(truth);
518  FillHistograms(truth);
519  evt.put(std::move(truthcol));
520 
521  return;
522  }
523 
524 // //......................................................................
525  std::string NUANCEGen::ParticleStatus(int StatusCode)
526  {
527  int code = StatusCode;
528  std::string ParticleStatusName;
529 
530  switch(code)
531  {
532  case 0:
533  ParticleStatusName = "kIStInitialState";
534  break;
535  case 1:
536  ParticleStatusName = "kIStFinalState";
537  break;
538  case 11:
539  ParticleStatusName = "kIStNucleonTarget";
540  break;
541  default:
542  ParticleStatusName = "Status Unknown";
543  }
544  return ParticleStatusName;
545  }
546 
547 
548 // //......................................................................
549  std::string NUANCEGen::ReactionChannel(int ccnc,int mode)
550  {
551  std::string ReactionChannelName=" ";
552 
553  if(ccnc==0)
554  ReactionChannelName = "kCC";
555  else if(ccnc==1)
556  ReactionChannelName = "kNC";
557  else std::cout<<"Current mode unknown!! "<<std::endl;
558 
559  if(mode==0)
560  ReactionChannelName += "_kQE";
561  else if(mode==1)
562  ReactionChannelName += "_kRes";
563  else if(mode==2)
564  ReactionChannelName += "_kDIS";
565  else if(mode==3)
566  ReactionChannelName += "_kCoh";
567  else if(mode==4)
568  ReactionChannelName += "_kNuElectronElastic";
569  else if(mode==5)
570  ReactionChannelName += "_kInverseMuDecay";
571  else std::cout<<"interaction mode unknown!! "<<std::endl;
572 
573  return ReactionChannelName;
574  }
575 
576 // //......................................................................
578  {
579  // Decide which histograms to put the spectrum in
580  int id = -1;
581  if (mc.GetNeutrino().CCNC()==simb::kCC) {
582  fCCMode->Fill(mc.GetNeutrino().Mode());
583  if (mc.GetNeutrino().Nu().PdgCode() == 12) id = 0;
584  else if (mc.GetNeutrino().Nu().PdgCode() == -12) id = 1;
585  else if (mc.GetNeutrino().Nu().PdgCode() == 14) id = 2;
586  else if (mc.GetNeutrino().Nu().PdgCode() == -14) id = 3;
587  else return;
588  }
589  else {
590  fNCMode->Fill(mc.GetNeutrino().Mode());
591  if (mc.GetNeutrino().Nu().PdgCode() > 0) id = 4;
592  else id = 5;
593  }
594  if (id==-1) abort();
595 
596  // Fill the specta histograms
597  fGenerated[id]->Fill(mc.GetNeutrino().Nu().E() );
598 
599  //< fill the vertex histograms from the neutrino - that is always
600  //< particle 0 in the list
601  simb::MCNeutrino mcnu = mc.GetNeutrino();
602  const simb::MCParticle nu = mcnu.Nu();
603 
604  fVertexX->Fill(nu.Vx());
605  fVertexY->Fill(nu.Vy());
606  fVertexZ->Fill(nu.Vz());
607 
608  fVertexXY->Fill(nu.Vx(), nu.Vy());
609  fVertexXZ->Fill(nu.Vz(), nu.Vx());
610  fVertexYZ->Fill(nu.Vz(), nu.Vy());
611 
612  double mom = nu.P();
613  if(std::abs(mom) > 0.){
614  fDCosX->Fill(nu.Px()/mom);
615  fDCosY->Fill(nu.Py()/mom);
616  fDCosZ->Fill(nu.Pz()/mom);
617  }
618 
619 
620 // LOG_DEBUG("GENIEInteractionInformation")
621 // << std::endl
622 // << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode())
623 // << std::endl
624 // << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl
625 // << std::setiosflags(std::ios::left)
626 // << std::setw(20) << "PARTICLE"
627 // << std::setiosflags(std::ios::left)
628 // << std::setw(32) << "STATUS"
629 // << std::setw(18) << "E (GeV)"
630 // << std::setw(18) << "m (GeV/c2)"
631 // << std::setw(18) << "Ek (GeV)"
632 // << std::endl << std::endl;
633 
634 // const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
635 
636 // // Loop over the particle stack for this event
637 // for(int i = 0; i < mc.NParticles(); ++i){
638 // simb::MCParticle part(mc.GetParticle(i));
639 // std::string name = databasePDG->GetParticle(part.PdgCode())->GetName();
640 // int code = part.StatusCode();
641 // std::string status = ParticleStatus(code);
642 // double mass = part.Mass();
643 // double energy = part.E();
644 // double Ek = (energy-mass); // Kinetic Energy (GeV)
645 // if(status=="kIStFinalStB4Interactions")
646 // LOG_DEBUG("GENIEFinalState")
647 // << std::setiosflags(std::ios::left) << std::setw(20) << name
648 // << std::setiosflags(std::ios::left) << std::setw(32) <<status
649 // << std::setw(18)<< energy
650 // << std::setw(18)<< mass
651 // << std::setw(18)<< Ek <<std::endl;
652 // else
653 // LOG_DEBUG("GENIEFinalState")
654 // << std::setiosflags(std::ios::left) << std::setw(20) << name
655 // << std::setiosflags(std::ios::left) << std::setw(32) << status
656 // << std::setw(18) << energy
657 // << std::setw(18) << mass <<std::endl;
658 
659  std::cout << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode()) << std::endl;
660  std::cout << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl;
661  std::cout << std::setiosflags(std::ios::left)
662  << std::setw(20) << "PARTICLE"
663  << std::setiosflags(std::ios::left)
664  << std::setw(32) << "STATUS"
665  << std::setw(18) << "E (GeV)"
666  << std::setw(18) << "m (GeV/c2)"
667  << std::setw(18) << "Ek (GeV)"
668  << std::endl << std::endl;
669 
670  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
671 
672  // Loop over the particle stack for this event
673  for(int i = 0; i < mc.NParticles(); ++i){
675  std::string name;
676  if (part.PdgCode() == 18040)
677  name = "Ar40 18040";
678  else if (part.PdgCode() != -99999 )
679  {
680  name = databasePDG->GetParticle(part.PdgCode())->GetName();
681  }
682 
683  int code = part.StatusCode();
684  std::string status = ParticleStatus(code);
685  double mass = part.Mass();
686  double energy = part.E();
687  double Ek = (energy-mass); // Kinetic Energy (GeV)
688 
689  std::cout << std::setiosflags(std::ios::left) << std::setw(20) << name
690  << std::setiosflags(std::ios::left) << std::setw(32) <<status
691  << std::setw(18)<< energy
692  << std::setw(18)<< mass
693  << std::setw(18)<< Ek <<std::endl;
694  }
695 
696 
697  if(mc.GetNeutrino().CCNC() == simb::kCC){
698 
701  for(int i = 0; i < mc.NParticles(); ++i){
703  if(abs(part.PdgCode()) == 11){
704  fEMomentum->Fill(part.P());
705  fEDCosX->Fill(part.Px()/part.P());
706  fEDCosY->Fill(part.Py()/part.P());
707  fEDCosZ->Fill(part.Pz()/part.P());
708  fECons->Fill(nu.E() - part.E());
709  break;
710  }
711  else if(abs(part.PdgCode()) == 13){
712  fMuMomentum->Fill(part.P());
713  fMuDCosX->Fill(part.Px()/part.P());
714  fMuDCosY->Fill(part.Py()/part.P());
715  fMuDCosZ->Fill(part.Pz()/part.P());
716  fECons->Fill(nu.E() - part.E());
717  break;
718  }
719  }// end loop over particles
720  }//end if CC interaction
721 
722  return;
723  }
724 
725 }
726 
727 namespace evgen{
728 
730 
731 }
732 
733 #endif
Float_t x
Definition: compare.C:6
double E(const int i=0) const
Definition: MCParticle.h:237
TH1F * fVertexX
vertex location of generated events in x
neutral current quasi-elastic
Definition: MCNeutrino.h:101
TH2F * fVertexYZ
vertex location in yz
int PdgCode() const
Definition: MCParticle.h:216
int CCNC() const
Definition: MCNeutrino.h:152
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:74
std::string fNuanceFile
double Py(const int i=0) const
Definition: MCParticle.h:235
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
charged current deep inelastic scatter
Definition: MCNeutrino.h:138
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:261
NUANCEGen(fhicl::ParameterSet const &pset)
resonant charged current, nubar p -> l+ p pi-
Definition: MCNeutrino.h:111
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:78
neutrino electron elastic scatter
Definition: MCNeutrino.h:144
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
TH2F * fVertexXY
vertex location in xy
Float_t y
Definition: compare.C:6
TH1F * fDCosY
direction cosine in y
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:150
TH1F * fEMomentum
momentum of outgoing electrons
Double_t z
Definition: plot.C:279
Float_t Y
Definition: plot.C:39
double Px(const int i=0) const
Definition: MCParticle.h:234
TH1F * fECons
histogram to determine if energy is conserved in the event
offset to account for adding in Nuance codes to this enum
Definition: MCNeutrino.h:99
TH1F * fGenerated[6]
Spectra as generated.
TH1F * fNCMode
CC interaction mode.
TH2F * fVertexXZ
vertex location in xz
int NParticles() const
Definition: MCTruth.h:72
void SetNeutrino(int CCNC, int mode, int interactionType, int target, int nucleon, int quark, double w, double x, double y, double qsqr)
Definition: MCTruth.cxx:30
charged current quasi-elastic
Definition: MCNeutrino.h:100
std::string ParticleStatus(int StatusCode)
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:148
void reconfigure(fhicl::ParameterSet const &p)
void Add(simb::MCParticle &part)
Definition: MCTruth.h:77
TH1F * fDCosZ
direction cosine in z
constexpr double PI
Definition: Run.h:30
Float_t Z
Definition: plot.C:39
TH1F * fCCMode
CC interaction mode.
TH1F * fDCosX
direction cosine in x
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
TH1F * fEDCosX
direction cosine of outgoing e in x
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
resonant charged current, nu n -> l- n pi+
Definition: MCNeutrino.h:104
void FillHistograms(simb::MCTruth mc)
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
TH1F * fMuDCosY
direction cosine of outgoing mu in y
std::string ReactionChannel(int ccnc, int mode)
double P(const int i=0) const
Definition: MCParticle.h:238
T get(std::string const &key) const
Definition: ParameterSet.h:231
TString part[npart]
Definition: Style.C:32
double energy
Definition: plottest35.C:25
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
std::unique_ptr< std::ifstream > fEventFile
TH1F * fEDCosZ
direction cosine of outgoing e in z
void produce(art::Event &evt)
TH1F * fMuMomentum
momentum of outgoing muons
charged current deep inelastic scatter
Definition: MCNeutrino.h:137
resonant charged current, nu p -> l- p pi+
Definition: MCNeutrino.h:102
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:73
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
double Vx(const int i=0) const
Definition: MCParticle.h:225
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
T * make(ARGS...args) const
ifstream in
Definition: comparison.C:7
charged current coherent pion
Definition: MCNeutrino.h:143
inverse muon decay
Definition: MCNeutrino.h:145
A module to check the results from the Monte Carlo generator.
double Pz(const int i=0) const
Definition: MCParticle.h:236
resonant charged current, nubar n -> l+ n pi-
Definition: MCNeutrino.h:109
TH1F * fMuDCosX
direction cosine of outgoing mu in x
double Vz(const int i=0) const
Definition: MCParticle.h:227
EventNumber_t event() const
Definition: EventID.h:117
TH1F * fEDCosY
direction cosine of outgoing e in y
Event generator information.
Definition: MCTruth.h:30
Event generator information.
Definition: MCNeutrino.h:18
Float_t X
Definition: plot.C:39
Namespace collecting geometry-related classes utilities.
int Mode() const
Definition: MCNeutrino.h:153
Event Generation using GENIE, cosmics or single particles.
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
void beginRun(art::Run &run)
TH1F * fVertexY
vertex location of generated events in y
EventID id() const
Definition: Event.h:56
double Vy(const int i=0) const
Definition: MCParticle.h:226
std::string fNUANCEModuleLabel
keep track of how long it takes to run the job
art framework interface to geometry description
TH1F * fVertexZ
vertex location of generated events in z
Beam neutrinos.
Definition: MCTruth.h:21