GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gEvScan.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gevscan
5 
6 \brief A utility that reads-in a GHEP event tree and performs basic sanity
7  checks / test whether the generated events obey basic conservation laws
8 
9 \syntax gevscan
10  -f ghep_event_file
11  [-o output_error_log_file]
12  [-n nev1[,nev2]]
13  [--add-event-printout-in-error-log]
14  [--max-num-of-errors-shown n]
15  [--event-record-print-level level]
16  [--check-energy-momentum-conservation]
17  [--check-charge-conservation]
18  [--check-for-pseudoparticles-in-final-state]
19  [--check-for-off-mass-shell-particles-in-final-state]
20  [--check-for-num-of-final-state-nucleons-inconsistent-with-target]
21  [--check-vertex-distribution]
22  [--check-decayer-consistency]
23  [--all]
24 
25 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
26  University of Liverpool
27 
28 \created August 13, 2008
29 
30 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
31  For the full text of the license visit http://copyright.genie-mc.org
32 
33 */
34 //____________________________________________________________________________
35 
36 //#define __debug__
37 
38 #include <string>
39 #include <vector>
40 #include <iomanip>
41 #include <sstream>
42 #include <fstream>
43 
44 #include <TSystem.h>
45 #include <TFile.h>
46 #include <TTree.h>
47 #include <TH1D.h>
48 #include <TLorentzVector.h>
49 
63 #include "Framework/Utils/RunOpt.h"
64 
65 using std::ostringstream;
66 using std::ofstream;
67 using std::string;
68 using std::vector;
69 using std::setw;
70 using std::setprecision;
71 using std::setfill;
72 using std::ios;
73 using std::endl;
74 
75 using namespace genie;
76 using namespace genie::constants;
77 
78 void GetCommandLineArgs (int argc, char ** argv);
79 void PrintSyntax (void);
80 bool CheckRootFilename (string filename);
81 
82 // checks
84 void CheckChargeConservation (void);
88 void CheckVertexDistribution (void);
89 void CheckDecayerConsistency (void);
90 
91 // options
92 string gOptInpFilename = "";
93 string gOptOutFilename = "";
94 Long64_t gOptNEvtL = -1;
95 Long64_t gOptNEvtH = -1;
96 int gOptMaxNumErrs = -1;
105 
106 Long64_t gFirstEventNum = -1;
107 Long64_t gLastEventNum = -1;
108 
109 TTree * gEventTree = 0;
111 ofstream gErrLog;
112 
113 //____________________________________________________________________________
114 int main(int argc, char ** argv)
115 {
116  GetCommandLineArgs (argc, argv);
117 
118  // Set GHEP print level
119  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
120 
121  TFile file(gOptInpFilename.c_str(),"READ");
122 
123  NtpMCTreeHeader * thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
124  LOG("gevscan", pINFO) << "Input tree header: " << *thdr;
125  NtpMCFormat_t format = thdr->format;
126  if(format != kNFGHEP) {
127  LOG("gevscan", pERROR)
128  << "*** Unsupported event-tree format : "
129  << NtpMCFormat::AsString(format);
130  file.Close();
131  return 3;
132  }
133 
134  gEventTree = dynamic_cast <TTree *> (file.Get("gtree"));
135  gEventTree->SetBranchAddress("gmcrec", &gMCRec);
136 
137  Long64_t nev = gEventTree->GetEntries();
138  if(gOptNEvtL == -1 && gOptNEvtH == -1) {
139  // read all events
140  gFirstEventNum = 0;
141  gLastEventNum = nev-1;
142  }
143  else {
144  // read a range of events
145  gFirstEventNum = TMath::Max((Long64_t)0, gOptNEvtL);
146  gLastEventNum = TMath::Min(nev-1, gOptNEvtH);
147  if(gLastEventNum - gFirstEventNum < 0) {
148  LOG("gevdump", pFATAL) << "Invalid event range";
149  PrintSyntax();
150  gAbortingInErr = true;
151  exit(1);
152  }
153  }
154 
155 
156  if(gOptOutFilename.size() == 0) {
157  ostringstream logfile;
158  logfile << gOptInpFilename << ".errlog";
159  gOptOutFilename = logfile.str();
160  }
161  if(gOptOutFilename != "none") {
162  gErrLog.open(gOptOutFilename.c_str());
163  gErrLog << "# ..................................................................................." << endl;
164  gErrLog << "# Error log for event file " << gOptInpFilename << endl;
165  gErrLog << "# ..................................................................................." << endl;
166  gErrLog << "# " << endl;
167  }
168 
171  }
174  }
177  }
180  }
183  }
186  }
189  }
190 
191 
192  if(gOptOutFilename != "none") {
193  gErrLog.close();
194  }
195 
196  return 0;
197 }
198 //____________________________________________________________________________
200 {
201  LOG("gevscan", pNOTICE) << "Checking energy/momentum conservation...";
202 
203  if(gErrLog.is_open()) {
204  gErrLog << "# Events failing the energy-momentum conservation test:" << endl;
205  gErrLog << "# " << endl;
206  }
207 
208  int nerr = 0;
209 
210  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
211  {
212  if(gOptMaxNumErrs != -1 && nerr >= gOptMaxNumErrs) break;
213 
214  gEventTree->GetEntry(i);
215 
216  NtpMCRecHeader rec_header = gMCRec->hdr;
217  EventRecord & event = *(gMCRec->event);
218 
219  LOG("gevscan", pINFO) << "Checking event.... " << i;
220 
221  double E_init = 0, E_fin = 0; // E
222  double px_init = 0, px_fin = 0; // px
223  double py_init = 0, py_fin = 0; // py
224  double pz_init = 0, pz_fin = 0; // pz
225 
226  GHepParticle * p = 0;
227  TIter event_iter(&event);
228  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
229 
230  GHepStatus_t ist = p->Status();
231 
232  if(ist == kIStInitialState)
233  {
234  E_init += p->E();
235  px_init += p->Px();
236  py_init += p->Py();
237  pz_init += p->Pz();
238  }
239  if(ist == kIStStableFinalState ||
241  {
242  E_fin += p->E();
243  px_fin += p->Px();
244  py_fin += p->Py();
245  pz_fin += p->Pz();
246  }
247  }//p
248 
249  double epsilon = 1E-3;
250 
251  bool E_conserved = TMath::Abs(E_init - E_fin) < epsilon;
252  bool px_conserved = TMath::Abs(px_init - px_fin) < epsilon;
253  bool py_conserved = TMath::Abs(py_init - py_fin) < epsilon;
254  bool pz_conserved = TMath::Abs(pz_init - pz_fin) < epsilon;
255 
256  bool ok = E_conserved &&
257  px_conserved &&
258  py_conserved &&
259  pz_conserved;
260 
261  if(!ok) {
262  LOG("gevscan", pERROR)
263  << " ** Energy-momentum non-conservation in event: " << i
264  << "\n"
265  << event;
266  if(gErrLog.is_open()) {
267  gErrLog << i;
269  gErrLog << event;
270  }
271  }
272  nerr++;
273  }
274 
275  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
276 
277  }//i
278 
279  if(gErrLog.is_open()) {
280  if(nerr == 0) {
281  gErrLog << "none" << endl;
282  }
283  }
284 
285  LOG("gevscan", pNOTICE)
286  << "Found " << nerr
287  << " events failing the energy/momentum conservation test";
288 }
289 //____________________________________________________________________________
291 {
292  LOG("gevscan", pNOTICE) << "Checking charge conservation...";
293 
294  if(gErrLog.is_open()) {
295  gErrLog << "# Events failing the charge conservation test:" << endl;
296  gErrLog << "# " << endl;
297  }
298 
299  int nerr = 0;
300 
301  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
302  {
303  if(gOptMaxNumErrs != -1 && nerr >= gOptMaxNumErrs) break;
304 
305  gEventTree->GetEntry(i);
306 
307  NtpMCRecHeader rec_header = gMCRec->hdr;
308  EventRecord & event = *(gMCRec->event);
309 
310  LOG("gevscan", pINFO) << "Checking event.... " << i;
311 
312  // Can't run the test for neutrinos scattered off nuclear targets
313  // because of intranuclear rescattering effects and the presence, in the event
314  // record, of a charged nuclear remnant pseudo-particle whose charge is not stored.
315  // To check charge conservation in the primary interaction, use a sample generated
316  // for a free nucleon targets.
317  GHepParticle * nucltgt = event.TargetNucleus();
318  if (nucltgt) {
319  LOG("gevscan", pINFO)
320  << "Event in nuclear target - Skipping test...";
321  }
322  else {
323  double Q_init = 0;
324  double Q_fin = 0;
325 
326  GHepParticle * p = 0;
327  TIter event_iter(&event);
328  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
329 
330  GHepStatus_t ist = p->Status();
331 
332  if(ist == kIStInitialState)
333  {
334  Q_init += p->Charge();
335  }
336  if(ist == kIStStableFinalState)
337  {
338  Q_fin += p->Charge();
339  }
340  }//p
341 
342  double epsilon = 1E-3;
343  bool ok = TMath::Abs(Q_init - Q_fin) < epsilon;
344  if(!ok) {
345  LOG("gevscan", pERROR)
346  << " ** Charge non-conservation in event: " << i
347  << "\n"
348  << event;
349  if(gErrLog.is_open()) {
350  gErrLog << i << endl;
352  gErrLog << event;
353  }
354  }
355  nerr++;
356  }
357 
358  }
359  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
360  }//i
361 
362  if(gErrLog.is_open()) {
363  if(nerr == 0) {
364  gErrLog << "none" << endl;
365  }
366  }
367 
368  LOG("gevscan", pNOTICE)
369  << "Found " << nerr
370  << " events failing the charge conservation test";
371 }
372 //____________________________________________________________________________
374 {
375  LOG("gevscan", pNOTICE)
376  << "Checking for pseudo-particles appearing in final state...";
377 
378  if(gErrLog.is_open()) {
379  gErrLog << "# Events with pseudo-particles in final state:" << endl;
380  gErrLog << "# " << endl;
381  }
382 
383  int nerr = 0;
384 
385  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
386  {
387  if(gOptMaxNumErrs != -1 && nerr >= gOptMaxNumErrs) break;
388 
389  gEventTree->GetEntry(i);
390 
391  NtpMCRecHeader rec_header = gMCRec->hdr;
392  EventRecord & event = *(gMCRec->event);
393 
394  LOG("gevscan", pINFO) << "Checking event.... " << i;
395 
396  GHepParticle * p = 0;
397  TIter event_iter(&event);
398  bool ok = true;
399  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
400 
401  GHepStatus_t ist = p->Status();
402  if(ist != kIStStableFinalState) continue;
403  int pdgc = p->Pdg();
404  if(pdg::IsPseudoParticle(pdgc))
405  {
406  ok = false;
407  break;
408  }
409  }//p
410 
411  if(!ok) {
412  LOG("gevscan", pERROR)
413  << " ** Pseudo-particle final state particle in event: " << i
414  << "\n"
415  << event;
416  if(gErrLog.is_open()) {
417  gErrLog << i << endl;
419  gErrLog << event;
420  }
421  }
422  nerr++;
423  }
424 
425  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
426 
427  }//i
428 
429  if(gErrLog.is_open()) {
430  if(nerr == 0) {
431  gErrLog << "none" << endl;
432  }
433  }
434 
435  LOG("gevscan", pNOTICE)
436  << "Found " << nerr
437  << " events with pseudo-particles in final state";
438 }
439 //____________________________________________________________________________
441 {
442  LOG("gevscan", pNOTICE)
443  << "Checking for off-mass-shell particles appearing in the final state...";
444 
445  if(gErrLog.is_open()) {
446  gErrLog << "# Events with off-mass-shell particles in final state:" << endl;
447  gErrLog << "# " << endl;
448  }
449 
450  int nerr = 0;
451 
452  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
453  {
454  if(gOptMaxNumErrs != -1 && nerr >= gOptMaxNumErrs) break;
455 
456  gEventTree->GetEntry(i);
457 
458  NtpMCRecHeader rec_header = gMCRec->hdr;
459  EventRecord & event = *(gMCRec->event);
460 
461  LOG("gevscan", pINFO) << "Checking event.... " << i;
462 
463  GHepParticle * p = 0;
464  TIter event_iter(&event);
465  bool ok = true;
466  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
467 
468  GHepStatus_t ist = p->Status();
469  if(ist != kIStStableFinalState) continue;
470  if(p->IsOffMassShell())
471  {
472  ok = false;
473  break;
474  }
475  }//p
476 
477  if(!ok) {
478  LOG("gevscan", pERROR)
479  << " ** Off-mass-shell final state particle in event: " << i
480  << "\n"
481  << event;
482  if(gErrLog.is_open()) {
483  gErrLog << i << endl;
485  gErrLog << event;
486  }
487  }
488  nerr++;
489  }
490 
491  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
492 
493  }//i
494 
495  if(gErrLog.is_open()) {
496  if(nerr == 0) {
497  gErrLog << "none" << endl;
498  }
499  }
500 
501  LOG("gevscan", pNOTICE)
502  << "Found " << nerr
503  << " events with off-mass-shell particles in final state";
504 }
505 //____________________________________________________________________________
507 {
508  LOG("gevscan", pNOTICE)
509  << "Checking for number of final state nucleons inconsistent with target...";
510 
511  if(gErrLog.is_open()) {
512  gErrLog << "# Events with number of final state nucleons inconsistent with target:" << endl;
513  gErrLog << "# " << endl;
514  }
515 
516  int nerr = 0;
517 
518  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
519  {
520  if(gOptMaxNumErrs != -1 && nerr >= gOptMaxNumErrs) break;
521 
522  gEventTree->GetEntry(i);
523 
524  NtpMCRecHeader rec_header = gMCRec->hdr;
525  EventRecord & event = *(gMCRec->event);
526 
527  LOG("gevscan", pINFO) << "Checking event.... " << i;
528 
529  // get target nucleus
530  GHepParticle * nucltgt = event.TargetNucleus();
531  if (!nucltgt) {
532  LOG("gevscan", pINFO)
533  << "Event not in nuclear target - Skipping test...";
534  }
535  else {
536  GHepParticle * p = 0;
537 
538  int Z = 0;
539  int N = 0;
540 
541  // get number of spectator nucleons
542  int fd = nucltgt->FirstDaughter();
543  int ld = nucltgt->LastDaughter();
544  for(int d = fd; d <= ld; d++) {
545  p = event.Particle(d);
546  if(!p) continue;
547  int pdgc = p->Pdg();
548  if(pdg::IsIon(pdgc)) {
549  Z = p->Z();
550  N = p->A() - p->Z();
551  }
552  }
553  // add nucleons from the primary interaction
554  TIter event_iter(&event);
555  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
556  GHepStatus_t ist = p->Status();
557  if(ist != kIStHadronInTheNucleus) continue;
558  int pdgc = p->Pdg();
559  if(pdg::IsProton (pdgc)) { Z++; }
560  if(pdg::IsNeutron(pdgc)) { N++; }
561  }//p
562 
563  LOG("gevscan", pINFO)
564  << "Before intranuclear hadron transport: Z = " << Z << ", N = " << N;
565 
566  // count final state nucleons
567  int Zf = 0;
568  int Nf = 0;
569  event_iter.Reset();
570  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
571  GHepStatus_t ist = p->Status();
572  if(ist != kIStStableFinalState) continue;
573  int pdgc = p->Pdg();
574  if(pdg::IsProton (pdgc)) { Zf++; }
575  if(pdg::IsNeutron(pdgc)) { Nf++; }
576  }
577  LOG("gevscan", pINFO)
578  << "In the final state: Z = " << Zf << ", N = " << Nf;
579 
580  bool ok = (Zf <= Z && Nf <= N);
581  if(!ok) {
582  LOG("gevscan", pERROR)
583  << " ** Number of final state nucleons inconsistent with target in event: " << i
584  << "\n"
585  << event;
586  if(gErrLog.is_open()) {
587  gErrLog << i << endl;
589  gErrLog << event;
590  }
591  }
592  nerr++;
593  }
594  } //nucltgt
595 
596  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
597 
598  }//i
599 
600 
601  if(gErrLog.is_open()) {
602  if(nerr == 0) {
603  gErrLog << "none" << endl;
604  }
605  }
606 
607  LOG("gevscan", pNOTICE)
608  << "Found " << nerr
609  << " events with a number of final state nucleons inconsistent with target";
610 }
611 //____________________________________________________________________________
613 {
614  LOG("gevscan", pNOTICE)
615  << "Checking intra-nuclear vertex distribution...";
616 
617  if(gErrLog.is_open()) {
618  gErrLog << "# Intranuclear vertex distribution check:" << endl;
619  gErrLog << "# " << endl;
620  }
621 
622  TH1D * r_distr_mc = new TH1D("r_distr_mc","", 150,0,30); //fm
623  TH1D * r_distr_expected = new TH1D("r_distr_expected","",150,0,30); //fm
624 
625  int Z = -1;
626  int A = -1;
627 
628  // get vertex position distribution
629  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
630  {
631  gEventTree->GetEntry(i);
632 
633  NtpMCRecHeader rec_header = gMCRec->hdr;
634  EventRecord & event = *(gMCRec->event);
635 
636  LOG("gevscan", pINFO) << "Checking event.... " << i;
637 
638  // get target nucleus
639  GHepParticle * nucltgt = event.TargetNucleus();
640  if (!nucltgt) {
641  LOG("gevscan", pINFO)
642  << "Event not in nuclear target - Skipping...";
643  }
644  else {
645  if(Z == -1 && A == -1) {
646  Z = nucltgt->Z();
647  A = nucltgt->A();
648  }
649 
650  // this test is run on a MC sample for a given target
651  if(Z != nucltgt->Z() || A != nucltgt->A()) {
652  LOG("gevscan", pINFO)
653  << "Event not in nuclear target seen first - Skipping...";
654  }
655  else {
656  GHepParticle * probe = event.Particle(0);
657  double r = probe->X4()->Vect().Mag();
658 
659  r_distr_mc->Fill(r);
660  }
661  } //nucltgt
662 
663  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
664 
665  }//i
666 
667  if(A > 1) {
668  // get expected vertex position distribution
669  for(int ir = 1; ir <= r_distr_expected->GetNbinsX(); ir++) {
670  double r = r_distr_expected->GetBinCenter(ir);
671  double rho = utils::nuclear::Density(r,A);
672  double nexp = 4*kPi*r*r*rho;
673  r_distr_expected->SetBinContent(ir,nexp);
674  }
675 
676  // normalize
677  double N = r_distr_mc->GetEntries();
678  r_distr_expected -> Scale (N / r_distr_expected -> Integral());
679 
680  // check consistency
681  double pvalue = r_distr_mc->Chi2Test(r_distr_expected,"WWP");
682  LOG("gevscan", pNOTICE) << "p-value {\\chi^2 test} = " << pvalue;
683 
684  if(gErrLog.is_open()) {
685  if(pvalue < 0.99) {
686  gErrLog << "Problem! p-value = " << pvalue << endl;
687  } else {
688  gErrLog << "OK! p-value = " << pvalue << endl;
689  }
690  }
691 
692 #ifdef __debug__
693  TFile f("./check_vtx.root","recreate");
694  r_distr_mc -> Write();
695  r_distr_expected -> Write();
696  f.Close();
697 #endif
698 
699  }//A
700  else {
701 
702  if(gErrLog.is_open()) {
703  gErrLog << "Can not run test with current sample" << endl;
704  }
705 
706  }
707 
708 }
709 //____________________________________________________________________________
711 {
712 // Check that particles seen in the final state in some events do not appear to
713 // have decayed in other events.
714 // This might happen if, for example, particle decay flags which are applied to
715 // GENIE events do not get applied to intermediate particles appearing in the
716 // PYTHIA hadronization. It might also happen if the decayed particle status is
717 // used incorrectly in some modules (eg intranuke).
718 //
719  LOG("gevscan", pNOTICE)
720  << "Checking decayer consistency...";
721 
722  if(gErrLog.is_open()) {
723  gErrLog << "# Decayer consistency check:" << endl;
724  gErrLog << "# " << endl;
725  }
726 
727  bool allowdup = false;
728  PDGCodeList final_state_particles(allowdup);
729  PDGCodeList decayed_particles(allowdup);
730 
731  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
732  {
733  gEventTree->GetEntry(i);
734 
735  NtpMCRecHeader rec_header = gMCRec->hdr;
736  EventRecord & event = *(gMCRec->event);
737 
738  LOG("gevscan", pINFO) << "Checking event.... " << i;
739 
740  GHepParticle * p = 0;
741  TIter event_iter(&event);
742  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
743  GHepStatus_t ist = p->Status();
744  int pdgc = p->Pdg();
745  if(ist == kIStStableFinalState) { final_state_particles.push_back(pdgc); }
746  if(ist == kIStDecayedState ) { decayed_particles.push_back(pdgc); }
747  }//p
748  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
749  }//i
750 
751  // find particles which appear in both lists
752  PDGCodeList particles_in_both_lists(allowdup);
753 
754  PDGCodeList::const_iterator iter;
755  for(iter = final_state_particles.begin();
756  iter != final_state_particles.end(); ++iter)
757  {
758  int pdgc = *iter;
759  if(decayed_particles.ExistsInPDGCodeList(pdgc))
760  {
761  particles_in_both_lists.push_back(pdgc);
762  }
763  }
764 
765  bool ok = true;
766  ostringstream mesg;
767  if(particles_in_both_lists.size() == 0) {
768  mesg << "OK.\n" << "No particle seen both in the final state and to have decayed.";
769  } else {
770  ok = false;
771  mesg << "Problem!\n" << particles_in_both_lists.size() << " particles seen both final state and to have decayed.";
772  }
773 
774  LOG("gevscan", pNOTICE)
775  << mesg.str();
776  LOG("gevscan", pNOTICE)
777  << "Particles seen in final state: " << final_state_particles;
778  LOG("gevscan", pNOTICE)
779  << "Particles seen to have decayed: " << decayed_particles;
780  LOG("gevscan", pNOTICE)
781  << "Particles seen in both lists: " << particles_in_both_lists;
782 
783  if(gErrLog.is_open()) {
784  gErrLog << mesg.str() << endl;
785  gErrLog << "\nParticles seen in final state:" << final_state_particles << endl;
786  gErrLog << "\nParticles seen to have decayed:" << decayed_particles << endl;
787  gErrLog << "\nParticles seen in both lists:" << particles_in_both_lists << endl;
788  }
789 
790  // find example events
791  if(!ok) {
792  if(gErrLog.is_open()) {
793  gErrLog << "\nExample events: " << endl;
794  }
795  for(iter = particles_in_both_lists.begin();
796  iter != particles_in_both_lists.end(); ++iter)
797  {
798  int pdgc_bothlists = *iter;
799  int iev_decay = -1;
800  int iev_fs = -1;
801  bool have_example = false;
802  for(Long64_t i = gFirstEventNum; i <= gLastEventNum; i++)
803  {
804  have_example = (iev_decay != -1 && iev_fs != -1);
805  if(have_example) break;
806 
807  gEventTree->GetEntry(i);
808  EventRecord & event = *(gMCRec->event);
809  GHepParticle * p = 0;
810  TIter event_iter(&event);
811  while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
812  int pdgc = p->Pdg();
813  if(pdgc != pdgc_bothlists) continue;
814  GHepStatus_t ist = p->Status();
815  if(ist == kIStStableFinalState && iev_fs == -1) { iev_fs = i; }
816  if(ist == kIStDecayedState && iev_decay == -1) { iev_decay = i; }
817  }//p
818  gMCRec->Clear(); // clear out explicitly to prevent memory leak w/Root6
819  }//i
820  if(gErrLog.is_open()) {
821  gErrLog << ">> " << PDGLibrary::Instance()->Find(pdgc_bothlists)->GetName()
822  << ": Decayed in event " << iev_decay
823  << ". Seen in final state in event " << iev_fs << "." << endl;
825  gEventTree->GetEntry(iev_decay);
826  EventRecord & event_dec = *(gMCRec->event);
827  gErrLog << "Event " << iev_decay << ":";
828  gErrLog << event_dec;
829  gEventTree->GetEntry(iev_fs);
830  EventRecord & event_fs = *(gMCRec->event);
831  gErrLog << "Event: " << iev_fs << ":";
832  gErrLog << event_fs;
833  }
834  }
835  }//pdgc
836  }//!ok
837 
838 }
839 //____________________________________________________________________________
840 void GetCommandLineArgs(int argc, char ** argv)
841 {
842  LOG("gevscan", pNOTICE) << "*** Parsing command line arguments";
843 
845 
846  CmdLnArgParser parser(argc,argv);
847 
848  // get input GENIE event sample
849  if( parser.OptionExists('f') ) {
850  LOG("gevscan", pINFO) << "Reading event sample filename";
851  gOptInpFilename = parser.ArgAsString('f');
852  } else {
853  LOG("gevscan", pFATAL)
854  << "Unspecified input filename - Exiting";
855  PrintSyntax();
856  exit(1);
857  }
858 
859  // get output error log
860  if( parser.OptionExists('o') ) {
861  LOG("gevscan", pINFO) << "Reading err log file name";
862  gOptOutFilename = parser.ArgAsString('o');
863  }
864 
865  // number of events
866  if ( parser.OptionExists('n') ) {
867  LOG("gevdump", pINFO) << "Reading number of events to analyze";
868  string nev = parser.ArgAsString('n');
869  if (nev.find(",") != string::npos) {
870  vector<long> vecn = parser.ArgAsLongTokens('n',",");
871  if(vecn.size()!=2) {
872  LOG("gevdump", pFATAL) << "Invalid syntax";
873  PrintSyntax();
874  gAbortingInErr = true;
875  exit(1);
876  }
877  // read a range of events
878  gOptNEvtL = vecn[0];
879  gOptNEvtH = vecn[1];
880  } else {
881  // read single event
882  gOptNEvtL = parser.ArgAsLong('n');
884  }
885  } else {
886  LOG("gevdump", pINFO)
887  << "Unspecified number of events to analyze - Use all";
888  gOptNEvtL = -1;
889  gOptNEvtH = -1;
890  }
891 
893  parser.OptionExists("add-event-printout-in-error-log");
894 
895  if(parser.OptionExists("max-num-of-errors-shown")) {
896  gOptMaxNumErrs = parser.ArgAsInt("max-num-of-errors-shown");
897  gOptMaxNumErrs = TMath::Max(1,gOptMaxNumErrs);
898  }
899 
900  bool all = parser.OptionExists("all");
901 
902  // checks
904  parser.OptionExists("check-energy-momentum-conservation");
906  parser.OptionExists("check-charge-conservation");
908  parser.OptionExists("check-for-num-of-final-state-nucleons-inconsistent-with-target");
910  parser.OptionExists("check-for-pseudoparticles-in-final-state");
912  parser.OptionExists("check-for-off-mass-shell-particles-in-final-state");
914  parser.OptionExists("check-vertex-distribution");
916  parser.OptionExists("check-decayer-consistency");
917 }
918 //____________________________________________________________________________
919 void PrintSyntax(void)
920 {
921  LOG("gevscan", pNOTICE)
922  << "\n\n" << "Syntax:" << "\n"
923  << " gevscan -f sample.root [-n n1[,n2]] [-o errlog] [check names]\n";
924 }
925 //_________________________________________________________________________________
926 bool CheckRootFilename(string filename)
927 {
928  if(filename.size() == 0) return false;
929 
930  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));
931  if (!is_accessible) {
932  LOG("gevscan", pERROR)
933  << "The input ROOT file [" << filename << "] is not accessible";
934  return false;
935  }
936  return true;
937 }
938 //_________________________________________________________________________________
939 
void CheckDecayerConsistency(void)
Definition: gEvScan.cxx:710
static void SetPrintLevel(int print_level)
Definition: GHepRecord.cxx:956
void CheckEnergyMomentumConservation(void)
Definition: gEvScan.cxx:199
int Z(void) const
bool gOptCheckVertexDistribution
Definition: gEvScan.cxx:103
NtpMCRecHeader hdr
record header
Definition: NtpMCRecordI.h:38
#define pERROR
Definition: Messenger.h:59
bool gOptCheckForPseudoParticlesInFinState
Definition: gEvScan.cxx:100
double E(void) const
Get energy.
Definition: GHepParticle.h:91
double Density(double r, int A, double ring=0.)
int FirstDaughter(void) const
Definition: GHepParticle.h:68
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object. Ntuples of this type are intended for feeding GENIE events into other applications (for example the GEANT4 based MC generation framework of an experiment) if no direct interface exists.
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:99
NtpMCEventRecord * gMCRec
Definition: gEvScan.cxx:110
#define pFATAL
Definition: Messenger.h:56
bool ExistsInPDGCodeList(int pdg_code) const
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
A list of PDG codes.
Definition: PDGCodeList.h:32
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
const double epsilon
int Pdg(void) const
Definition: GHepParticle.h:63
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
int gOptMaxNumErrs
Definition: gEvScan.cxx:96
Long64_t gFirstEventNum
Definition: gEvScan.cxx:106
int LastDaughter(void) const
Definition: GHepParticle.h:69
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
void CheckForOffMassShellParticlesInFinState(void)
Definition: gEvScan.cxx:440
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
Long64_t gOptNEvtL
Definition: gEvDump.cxx:74
MINOS-style Ntuple Class to hold an output MC Tree Header.
Long64_t gLastEventNum
Definition: gEvScan.cxx:107
void CheckForPseudoParticlesInFinState(void)
Definition: gEvScan.cxx:373
double Charge(void) const
Chrg that corresponds to the PDG code.
#define pINFO
Definition: Messenger.h:62
bool CheckRootFilename(string filename)
Definition: gEvComp.cxx:2033
static const char * AsString(NtpMCFormat_t fmt)
Definition: NtpMCFormat.h:36
void CheckForNumFinStateNucleonsInconsistentWithTarget(void)
Definition: gEvScan.cxx:506
void CheckVertexDistribution(void)
Definition: gEvScan.cxx:612
bool gOptCheckForNumFinStateNucleonsInconsistentWithTarget
Definition: gEvScan.cxx:102
bool gOptCheckChargeConservation
Definition: gEvScan.cxx:99
bool IsOffMassShell(void) const
void CheckChargeConservation(void)
Definition: gEvScan.cxx:290
bool gOptCheckForOffMassShellParticlesInFinState
Definition: gEvScan.cxx:101
string gOptInpFilename
Definition: gEvDump.cxx:76
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
bool gOptCheckDecayerConsistency
Definition: gEvScan.cxx:104
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
static RunOpt * Instance(void)
Definition: RunOpt.cxx:54
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:42
ofstream gErrLog
Definition: gEvScan.cxx:111
bool gOptCheckEnergyMomentumConservation
Definition: gEvScan.cxx:98
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
bool IsPseudoParticle(int pdgc)
Definition: PDGUtils.cxx:27
bool gOptAddEventPrintoutInErrLog
Definition: gEvScan.cxx:97
Long64_t gOptNEvtH
Definition: gEvDump.cxx:75
hadnt Write("hadnt")
MINOS-style Ntuple Class to hold an MC Event Record Header.
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
Command line argument parser.
enum genie::ENtpMCFormat NtpMCFormat_t
int A(void) const
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
void Clear(Option_t *opt="")
bool gAbortingInErr
Definition: Messenger.cxx:34
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TTree * gEventTree
Definition: gEvScan.cxx:109
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
void PrintSyntax(void)
EventRecord * event
event
string gOptOutFilename
Definition: gEvScan.cxx:93
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
enum genie::EGHepStatus GHepStatus_t
double Py(void) const
Get Py.
Definition: GHepParticle.h:89