Project

General

Profile

find_bad_nuance.C

script to process events and identify problematic ones and print them out - Robert Hatcher, 11/22/2017 04:18 PM

 
1
//
2
// Prints the events in the input GHEP event file
3
//   *for events that return a nuance code of 0"
4
//
5
// Usage:
6
// shell% genie -b -q find_bad_nuance.C\(\"icarus.12345.ghep.root\"\)
7
//
8
// simple dumper by: Costas Andreopoulos, Dec 13, 2006
9
// adapted by Robert Hatcher 2017-11-22
10
//
11
#include <iostream>
12
using namespace std;
13
#include "TFile.h"
14
#include "TTree.h"
15
#include "GENIE/Ntuple/NtpMCEventRecord.h"
16
#include "GENIE/EVGCore/EventRecord.h"
17
#include "GENIE/Interaction/Interaction.h"
18
#include "GENIE/Interaction/ProcessInfo.h"
19
#include "GENIE/Interaction/InitialState.h"
20
#include "GENIE/GHEP/GHepParticle.h"
21
#include "GENIE/PDG/PDGCodes.h"
22
#include "GENIE/GHEP/GHepUtils.h"   // NuanceReactionCode() function
23

    
24
// forward declare our own version
25
namespace genie {
26
  namespace utils {
27
    namespace ghep {
28
      int NuanceReactionCodeModified(const genie::GHepRecord * event);
29
    }
30
  }
31
}
32

    
33
void find_bad_nuance(const char * filename)
34
{
35

    
36
 // open the ROOT file and get the event TTree
37
 TFile inpfile(filename,"READ");
38
 TTree * tree = dynamic_cast <TTree *> (inpfile.Get("gtree"));
39

    
40
 // set the branch address
41
 genie::NtpMCEventRecord * mcrec = 0;
42
 tree->SetBranchAddress("gmcrec", &mcrec);
43

    
44
 int nchecked=0;
45
 int nmec=0;
46
 int nres=0;
47
 int nother=0;
48

    
49
 // loop over event tree
50
 for(int i = 0; i< tree->GetEntries(); i++) {
51
    tree->GetEntry(i);
52

    
53
    genie::NtpMCRecHeader rec_header = mcrec->hdr;
54
    genie::EventRecord &  event      = *(mcrec->event);
55

    
56
    // Nuance function needs  const GHepRecord*
57
    //   genie::EventRecord derives from GHepRecord
58
    // we'll make additional variables so we can cut-and-paste
59
    // other code in to this routine
60
    genie::GHepRecord* recptr = &event;
61

    
62
    int nuance = genie::utils::ghep::NuanceReactionCode(recptr);
63
    ++nchecked;
64

    
65
    if ( nuance == 0 ) {
66

    
67
      genie::Interaction * interaction = recptr->Summary();
68
      const genie::ProcessInfo &  proc = interaction->ProcInfo();
69
      const genie::InitialState & init = interaction->InitState();
70
      if ( proc.IsMEC() ) {
71
        ++nmec;
72
        cout << rec_header << " <==== is MEC" << endl;
73
      } else {
74
        // the more complicated unhandled case
75

    
76
        // print-out
77
        cout << rec_header;
78
        cout << event;
79

    
80
        int nuance2 = genie::utils::ghep::NuanceReactionCodeModified(recptr);
81
        if ( proc.IsResonant() ) {
82
          ++nres;
83
        } else {
84
          ++nother;
85
        }
86
      }
87
    }
88

    
89
    mcrec->Clear();
90
 }
91

    
92
 inpfile.Close();
93

    
94
 cout << "Final results: " << nchecked << " checked, "
95
      << nmec << " MEC, " << nres << " Resonant, "
96
      << nother << " other" << endl;
97
}
98

    
99

    
100
//____________________________________________________________________________
101
int genie::utils::ghep::NuanceReactionCodeModified(const genie::GHepRecord * event)
102
{
103
// Josh Spitz, Costas Andreopoulos
104
//
105
  if(!event) {
106
    //LOG("GHepUtils", pWARN) << "Null event!";
107
    return 0;
108
  }
109

    
110
  int evtype = 0;
111

    
112
  genie::Interaction * interaction = event->Summary();
113

    
114
  const genie::ProcessInfo &  proc = interaction->ProcInfo();
115
  const genie::InitialState & init = interaction->InitState();
116
  if      (proc.IsQuasiElastic()   && proc.IsWeakCC()) evtype =  1;
117
  else if (proc.IsQuasiElastic()   && proc.IsWeakNC()) evtype =  2;
118
  else if (proc.IsDeepInelastic()  && proc.IsWeakCC()) evtype = 91;
119
  else if (proc.IsDeepInelastic()  && proc.IsWeakNC()) evtype = 92;
120
  else if (proc.IsCoherent()       && proc.IsWeakNC()) evtype = 96;
121
  else if (proc.IsCoherent()       && proc.IsWeakCC()) evtype = 97;
122
  else if (proc.IsNuElectronElastic())                 evtype = 98;
123
  else if (proc.IsInverseMuDecay())                    evtype = 99;
124
  else if (proc.IsResonant()) {
125
     int nn=0, np=0, npi0=0, npip=0, npim=0;
126
     bool nuclear_target = init.Tgt().IsNucleus();
127
     genie::GHepStatus_t matched_ist = (nuclear_target) ?
128
                kIStHadronInTheNucleus : kIStStableFinalState;
129

    
130
     TIter event_iter(event);
131
     genie::GHepParticle * p = 0;
132

    
133
     while ( (p = dynamic_cast<genie::GHepParticle *>(event_iter.Next())) )
134
     {
135
         genie::GHepStatus_t ghep_ist = (genie::GHepStatus_t) p->Status();
136
         if(ghep_ist != matched_ist) continue;
137

    
138
         int ghep_pdgc = p->Pdg();
139
         if(ghep_pdgc == genie::kPdgProton ) np++;
140
         if(ghep_pdgc == genie::kPdgNeutron) nn++;
141
         if(ghep_pdgc == genie::kPdgPi0)     npi0++;
142
         if(ghep_pdgc == genie::kPdgPiP)     npip++;
143
         if(ghep_pdgc == genie::kPdgPiM)     npim++;
144
     }
145
     if(proc.IsWeakCC() && init.IsNuP()) {
146
         // v p -> l- p pi+
147
         if(np==1 && nn==0 && npip==1 && npi0==0 && npim==0) evtype = 3;
148
     }
149
     if(proc.IsWeakCC() && init.IsNuN()) {
150
         // v n -> l- p pi0
151
         if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 4;
152
         // v n -> l- n pi+
153
         if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 5;
154
     }
155
     if(proc.IsWeakNC() && init.IsNuP()) {
156
         // v p -> v p pi0
157
         if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 6;
158
         // v p -> v n pi+
159
         if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 7;
160
     }
161
     if(proc.IsWeakNC() && init.IsNuN()) {
162
         // v n -> v n pi0
163
         if(np==0 && nn==1 && npip==0 && npi0==1 && npim==0) evtype = 8;
164
         // v n -> v p pi-
165
         if(np==1 && nn==0 && npip==0 && npi0==0 && npim==1) evtype = 9;
166
     }
167
     if(proc.IsWeakCC() && init.IsNuBarN()) {
168
         // \bar{v} n -> l+ n pi-
169
         if(np==1 && nn==0 && npip==1 && npi0==0 && npim==0) evtype = 10;
170
     }
171
     if(proc.IsWeakCC() && init.IsNuBarP()) {
172
         // \bar{v} p -> l+ n pi0
173
         if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 11;
174
         // \bar{v} p -> l+ p pi-
175
         if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 12;
176
     }
177
     if(proc.IsWeakNC() && init.IsNuBarP()) {
178
         // \bar{v} p -> \bar{v} p pi0
179
         if(np==1 && nn==0 && npip==0 && npi0==1 && npim==0) evtype = 13;
180
         // \bar{v} p -> \bar{v} n pi+
181
         if(np==0 && nn==1 && npip==1 && npi0==0 && npim==0) evtype = 14;
182
      }
183
      if(proc.IsWeakNC() && init.IsNuBarN()) {
184
         // \bar{v} n -> \bar{v} n pi0
185
         if(np==0 && nn==1 && npip==0 && npi0==1 && npim==0) evtype = 15;
186
         // \bar{v} n -> \bar{v} p pi-
187
         if(np==1 && nn==0 && npip==0 && npi0==0 && npim==1) evtype = 16;
188
      }
189
      cout << "FOR CONVENIENCE: IsResonant " << proc.IsResonant() << endl
190
           << " np " << np << " nn " << nn
191
           << " npip " << npip << " npim " << npim
192
           << " npi0 " << npi0 << endl
193
           << " IsWeakNC() " << proc.IsWeakNC()
194
           << " IsWeakCC() " << proc.IsWeakCC() << endl
195
           << " IsNuP() " << init.IsNuP()
196
           << " IsNuN() " << init.IsNuN() << endl
197
           << " IsNuBarP() " << init.IsNuBarP()
198
           << " IsNuBarN() " << init.IsNuBarN() << endl;
199
  } else {
200
    cout << "FOR CONVENIENCE: IsResonant " << proc.IsResonant()
201
         << " ... fell through " << endl;
202
  }
203

    
204
  return evtype;
205
}
206
//____________________________________________________________________________