All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
evgen::NDKGen Class Reference

A module to check the results from the Monte Carlo generator. More...

Inheritance diagram for evgen::NDKGen:

Public Member Functions

 NDKGen (fhicl::ParameterSet const &pset)
 
virtual ~NDKGen ()
 

Private Member Functions

void produce (art::Event &evt) override
 
void beginJob () override
 
void beginRun (art::Run &run) override
 
void endJob () override
 
std::string ParticleStatus (int StatusCode)
 
std::string ReactionChannel (int ccnc, int mode)
 
void FillHistograms (simb::MCTruth const &mc)
 

Private Attributes

std::string fNdkFile
 
std::ifstream fEventFile
 
TStopwatch fStopwatch
 
std::string fNDKModuleLabel
 keep track of how long it takes to run the job More...
 
CLHEP::HepRandomEngine & fEngine
 art-managed random-number engine More...
 
TH1F * fGenerated [6]
 Spectra as generated. More...
 
TH1F * fVertexX
 vertex location of generated events in x More...
 
TH1F * fVertexY
 vertex location of generated events in y More...
 
TH1F * fVertexZ
 vertex location of generated events in z More...
 
TH2F * fVertexXY
 vertex location in xy More...
 
TH2F * fVertexXZ
 vertex location in xz More...
 
TH2F * fVertexYZ
 vertex location in yz More...
 
TH1F * fDCosX
 direction cosine in x More...
 
TH1F * fDCosY
 direction cosine in y More...
 
TH1F * fDCosZ
 direction cosine in z More...
 
TH1F * fMuMomentum
 momentum of outgoing muons More...
 
TH1F * fMuDCosX
 direction cosine of outgoing mu in x More...
 
TH1F * fMuDCosY
 direction cosine of outgoing mu in y More...
 
TH1F * fMuDCosZ
 direction cosine of outgoing mu in z More...
 
TH1F * fEMomentum
 momentum of outgoing electrons More...
 
TH1F * fEDCosX
 direction cosine of outgoing e in x More...
 
TH1F * fEDCosY
 direction cosine of outgoing e in y More...
 
TH1F * fEDCosZ
 direction cosine of outgoing e in z More...
 
TH1F * fCCMode
 CC interaction mode. More...
 
TH1F * fNCMode
 CC interaction mode. More...
 
TH1F * fECons
 histogram to determine if energy is conserved in the event More...
 

Detailed Description

A module to check the results from the Monte Carlo generator.

Definition at line 49 of file NDKGen_module.cc.

Constructor & Destructor Documentation

evgen::NDKGen::NDKGen ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 108 of file NDKGen_module.cc.

109  : EDProducer{pset}
110  , fNdkFile{pset.get<std::string>("NdkFile")}
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"))
std::string fNdkFile
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
std::ifstream fEventFile
evgen::NDKGen::~NDKGen ( )
virtual

Definition at line 129 of file NDKGen_module.cc.

130  {
131  fStopwatch.Stop();
132  }
TStopwatch fStopwatch

Member Function Documentation

void evgen::NDKGen::beginJob ( )
overrideprivate

Definition at line 136 of file NDKGen_module.cc.

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  }
TH1F * fMuMomentum
momentum of outgoing muons
process_name opflash particleana ie ie ie z
TH2F * fVertexXY
vertex location in xy
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
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexXZ
vertex location in xz
TH1F * fVertexY
vertex location of generated events in y
TH1F * fCCMode
CC interaction mode.
TH1F * fEDCosY
direction cosine of outgoing e in y
process_name opflash particleana ie ie y
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
TH1F * fEDCosX
direction cosine of outgoing e in x
art::ServiceHandle< art::TFileService > tfs
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosY
direction cosine in y
TH1F * fVertexX
vertex location of generated events in x
TH1F * fDCosX
direction cosine in x
void evgen::NDKGen::beginRun ( art::Run &  run)
overrideprivate

Definition at line 197 of file NDKGen_module.cc.

198  {
199  art::ServiceHandle<geo::Geometry const> geo;
200  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
201  }
void evgen::NDKGen::endJob ( )
overrideprivate

Definition at line 204 of file NDKGen_module.cc.

205  {
206  fEventFile.close();
207  }
std::ifstream fEventFile
void evgen::NDKGen::FillHistograms ( simb::MCTruth const &  mc)
private

look for the outgoing lepton in the particle stack just interested in the first one

Definition at line 466 of file NDKGen_module.cc.

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  }
TH1F * fMuMomentum
momentum of outgoing muons
TH2F * fVertexXY
vertex location in xy
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
TH1F * fEDCosZ
direction cosine of outgoing e in z
TH1F * fEMomentum
momentum of outgoing electrons
TH2F * fVertexXZ
vertex location in xz
TH1F * fVertexY
vertex location of generated events in y
TH1F * fCCMode
CC interaction mode.
TH1F * fEDCosY
direction cosine of outgoing e in y
T abs(T value)
Charged-current interactions.
Definition: IPrediction.h:38
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
walls no left
Definition: selectors.fcl:105
std::string ReactionChannel(int ccnc, int mode)
TH1F * fEDCosX
direction cosine of outgoing e in x
float mass
Definition: dedx.py:47
then echo fcl name
std::string ParticleStatus(int StatusCode)
TH1F * fNCMode
CC interaction mode.
TH1F * fDCosY
direction cosine in y
BEGIN_PROLOG SN nu
TH1F * fVertexX
vertex location of generated events in x
TH1F * fDCosX
direction cosine in x
BEGIN_PROLOG could also be cout
std::string evgen::NDKGen::ParticleStatus ( int  StatusCode)
private

Definition at line 414 of file NDKGen_module.cc.

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  }
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
void evgen::NDKGen::produce ( art::Event &  evt)
overrideprivate

Definition at line 210 of file NDKGen_module.cc.

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  }
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
unsigned int event
Definition: DataStructs.h:634
process_name opdaq physics producers generator physics producers generator physics producers generator Z0
Definition: gen_protons.fcl:45
Geometry information for a single TPC.
Definition: TPCGeo.h:38
process_name opdaq physics producers generator physics producers generator Y0
Definition: gen_protons.fcl:45
process_name E
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
CLHEP::HepRandomEngine & fEngine
art-managed random-number engine
then local
standard_singlep gaussian distribution X0
Definition: multigen.fcl:8
BEGIN_PROLOG vertical distance to the surface Name
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
then echo fcl name
TCEvent evt
Definition: DataStructs.cxx:8
pdgs k
Definition: selectors.fcl:22
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
physics associatedGroupsWithLeft p1
BEGIN_PROLOG could also be cout
std::string evgen::NDKGen::ReactionChannel ( int  ccnc,
int  mode 
)
private

Definition at line 438 of file NDKGen_module.cc.

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  }
const char mode
Definition: noise_ana.cxx:20
BEGIN_PROLOG could also be cout

Member Data Documentation

TH1F* evgen::NDKGen::fCCMode
private

CC interaction mode.

Definition at line 97 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fDCosX
private

direction cosine in x

Definition at line 83 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fDCosY
private

direction cosine in y

Definition at line 84 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fDCosZ
private

direction cosine in z

Definition at line 85 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fECons
private

histogram to determine if energy is conserved in the event

Definition at line 100 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fEDCosX
private

direction cosine of outgoing e in x

Definition at line 93 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fEDCosY
private

direction cosine of outgoing e in y

Definition at line 94 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fEDCosZ
private

direction cosine of outgoing e in z

Definition at line 95 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fEMomentum
private

momentum of outgoing electrons

Definition at line 92 of file NDKGen_module.cc.

CLHEP::HepRandomEngine& evgen::NDKGen::fEngine
private

art-managed random-number engine

Definition at line 71 of file NDKGen_module.cc.

std::ifstream evgen::NDKGen::fEventFile
private

Definition at line 67 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fGenerated[6]
private

Spectra as generated.

Definition at line 73 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fMuDCosX
private

direction cosine of outgoing mu in x

Definition at line 88 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fMuDCosY
private

direction cosine of outgoing mu in y

Definition at line 89 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fMuDCosZ
private

direction cosine of outgoing mu in z

Definition at line 90 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fMuMomentum
private

momentum of outgoing muons

Definition at line 87 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fNCMode
private

CC interaction mode.

Definition at line 98 of file NDKGen_module.cc.

std::string evgen::NDKGen::fNdkFile
private

Definition at line 66 of file NDKGen_module.cc.

std::string evgen::NDKGen::fNDKModuleLabel
private

keep track of how long it takes to run the job

Definition at line 70 of file NDKGen_module.cc.

TStopwatch evgen::NDKGen::fStopwatch
private

Definition at line 68 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fVertexX
private

vertex location of generated events in x

Definition at line 75 of file NDKGen_module.cc.

TH2F* evgen::NDKGen::fVertexXY
private

vertex location in xy

Definition at line 79 of file NDKGen_module.cc.

TH2F* evgen::NDKGen::fVertexXZ
private

vertex location in xz

Definition at line 80 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fVertexY
private

vertex location of generated events in y

Definition at line 76 of file NDKGen_module.cc.

TH2F* evgen::NDKGen::fVertexYZ
private

vertex location in yz

Definition at line 81 of file NDKGen_module.cc.

TH1F* evgen::NDKGen::fVertexZ
private

vertex location of generated events in z

Definition at line 77 of file NDKGen_module.cc.


The documentation for this class was generated from the following file: