Project

General

Profile

UtlEventRec.cxx

Function mapping GENIE vertices to VALOR Interaction codes - Stephen Dennis, 04/25/2017 04:21 PM

 
1
//________________________________________________________________________________________
2
/*
3

4
Copyright (c) 2010-2016, The VALOR Fitting Group, All Rights Reserved.
5

6
*/
7
//________________________________________________________________________________________
8

    
9
#include <vector>
10
#include <iostream>
11

    
12
#include "TObjArray.h"
13
#include "TString.h"
14
#include "TMath.h"
15

    
16
#include "Conventions/GBuild.h"
17
#include "Conventions/Units.h"
18
#include "EVGCore/EventRecord.h"
19
#include "GHEP/GHepParticle.h"
20
#include "PDG/PDGCodes.h"
21
#include "PDG/PDGUtils.h"
22
#include "PDG/PDGLibrary.h"
23

    
24
// valor/dune
25
#include "valor/Conventions/PDG.h"
26
#include "valor/Messenger/Messenger.h"
27
#include "valor/Utilities/UtlSystem.h"
28

    
29
#include "dune/Conventions/BeamConf.h"
30
#include "dune/Conventions/ReacMode.h"
31
#include "dune/Conventions/Sample.h"
32

    
33
#include "dune/GENIEUtilities/UtlEventRec.h"
34

    
35
using namespace valor;
36

    
37
bool valor::util::GetHadronicSystems(const genie::EventRecord * event, vector<int> & prim_had_syst, vector<int> & final_had_syst) 
38
{
39
  // Function adapted from GENIE app gNtpConv (GST).
40
  prim_had_syst .clear();
41
  final_had_syst.clear();
42
  
43
  // Process id
44
  const genie::ProcessInfo &  proc_info = event->Summary()->ProcInfo();
45
  bool is_qel     = proc_info.IsQuasiElastic();
46
  bool is_res     = proc_info.IsResonant();
47
  bool is_dis     = proc_info.IsDeepInelastic();
48
  bool is_coh     = proc_info.IsCoherent();
49
  bool is_dfr     = proc_info.IsDiffractive();
50
  bool is_singlek = proc_info.IsSingleKaon();    
51
  bool is_mec     = proc_info.IsMEC();
52

    
53
  // Extract more info on the hadronic system
54
  // Only for QEL/RES/DIS/COH/MEC events
55
  //
56
  bool study_hadsyst = (is_qel || is_res || is_dis || is_coh || is_dfr || is_mec || is_singlek);
57
  if (!study_hadsyst) {return false;}
58

    
59
  //
60
  TObjArrayIter piter(event);
61
  genie::GHepParticle * p = 0;
62
  int ip=-1;
63

    
64
  //
65
  // Extract the final state system originating from the hadronic vertex 
66
  // (after the intranuclear rescattering step)
67
  //
68

    
69
  LOG("UtlEventRec", pDEBUG) << "Extracting final state hadronic system";
70

    
71
  while( (p = (genie::GHepParticle *) piter.Next()) && study_hadsyst)
72
  {
73
    ip++;
74
    // don't count final state lepton as part hadronic system 
75
    //if(!is_coh && event->Particle(ip)->FirstMother()==0) continue;
76
    if(event->Particle(ip)->FirstMother()==0) continue;
77
    if(genie::pdg::IsPseudoParticle(p->Pdg())) continue;
78
    int pdgc = p->Pdg();
79
    int ist  = p->Status();
80
    if(ist==genie::kIStStableFinalState) {
81
      if (pdgc == genie::kPdgGamma || pdgc == genie::kPdgElectron || pdgc == genie::kPdgPositron)  {
82
        int igmom = p->FirstMother();
83
        if(igmom!=-1) {
84
          // only count e+'s e-'s or gammas not from decay of pi0
85
          if(event->Particle(igmom)->Pdg() != genie::kPdgPi0) { final_had_syst.push_back(ip); }
86
        }
87
      } else {
88
        final_had_syst.push_back(ip);
89
      }
90
    }
91
    // now add pi0's that were decayed as short lived particles
92
    else if(pdgc == genie::kPdgPi0){
93
      int ifd = p->FirstDaughter();
94
      int fd_pdgc = event->Particle(ifd)->Pdg();
95
      // just require that first daughter is one of gamma, e+ or e-  
96
      if(fd_pdgc == genie::kPdgGamma || fd_pdgc == genie::kPdgElectron || fd_pdgc == genie::kPdgPositron){
97
        final_had_syst.push_back(ip);
98
      }
99
    }
100
  }//particle-loop
101

    
102
  //~ if( count(final_had_syst.begin(), final_had_syst.end(), -1) > 0) {
103
    //~ mcrec->Clear();
104
    //~ continue;
105
  //~ }
106

    
107
  //
108
  // Extract info on the primary hadronic system (before any intranuclear rescattering)
109
  // looking for particles with status_code == kIStHadronInTheNucleus 
110
  // An exception is the coherent production and scattering off free nucleon targets 
111
  // (no intranuclear rescattering) in which case primary hadronic system is set to be 
112
  // 'identical' with the final  state hadronic system
113
  //
114

    
115
  LOG("UtlEventRec", pDEBUG) << "Extracting primary hadronic system";
116

    
117
  ip = -1;
118
  TObjArrayIter piter_prim(event);
119
  
120
  if(study_hadsyst) {
121
    // if coherent or free nucleon target set primary states equal to final states
122
    genie::GHepParticle * target = event->Particle(1);
123
    if(!genie::pdg::IsIon(target->Pdg()) || (is_coh)) {
124
      vector<int>::const_iterator hiter = final_had_syst.begin();
125
      for( ; hiter != final_had_syst.end(); ++hiter) {
126
        prim_had_syst.push_back(*hiter);
127
      }
128
    } 
129
    // otherwise loop over all particles and store indices of those which are hadrons
130
    // created within the nucleus
131
    else {
132
      while( (p = (genie::GHepParticle *) piter_prim.Next()) ){
133
        ip++;      
134
        int ist_comp  = p->Status();
135
        if(ist_comp==genie::kIStHadronInTheNucleus) {
136
          prim_had_syst.push_back(ip); 
137
        }
138
      }//particle-loop   
139
      //
140
      // also include gammas from nuclear de-excitations (appearing in the daughter list of the 
141
      // hit nucleus, earlier than the primary hadronic system extracted above)
142
      for(int i = target->FirstDaughter(); i <= target->LastDaughter(); i++) {
143
        if(i<0) continue;
144
        if(event->Particle(i)->Status()==genie::kIStStableFinalState) { prim_had_syst.push_back(i); }
145
      }      
146
    }//freenuc?
147
  }//study_hadsystem?
148
  
149
  return study_hadsyst;
150
}
151
//_________________________________________________________________________________
152
bool valor::util::EventContainsQuarkFlavour(const genie::EventRecord * event, int quark_pdg)
153
{
154
  //
155
  TObjArrayIter piter(event);
156
  genie::GHepParticle * p = 0;
157
  
158
  bool result = false;
159
  
160
  while( (p = (genie::GHepParticle *) piter.Next())) {
161
    int pdg = std::abs(p->Pdg());
162
    pdg = std::abs(pdg);
163
    
164
    // Check if PDG code contains the quark PDG
165
    if (genie::pdg::IsHadron(pdg)) {
166
      int magnitude = 1000;
167
      while(magnitude) {
168
        int digit = (pdg / magnitude) % 10;
169
        if (digit==std::abs(quark_pdg)) {result = true;}
170
        magnitude/=10;
171
      }
172
    }
173
  } //piter
174
  return result;
175
}
176
//_________________________________________________________________________________
177
conv::ReacMode_t valor::util::DetermineReactionMode(const genie::EventRecord * event, bool use_final_state)
178
{
179
  int  nupdg       = event->Probe()->Pdg();
180
  bool is_nue      = (nupdg == genie::kPdgNuE);
181
  bool is_numu     = (nupdg == genie::kPdgNuMu);
182
  bool is_nuebar   = (nupdg == genie::kPdgAntiNuE);
183
  bool is_numubar  = (nupdg == genie::kPdgAntiNuMu);
184

    
185
  bool is_CC       = event->Summary()->ProcInfo().IsWeakCC();
186
  bool is_NC       = event->Summary()->ProcInfo().IsWeakNC();
187
  bool is_charm   = util::EventContainsQuarkFlavour(event, genie::kPdgCQuark);
188
  bool is_strange = util::EventContainsQuarkFlavour(event, genie::kPdgSQuark);
189
  bool is_QE       = event->Summary()->ProcInfo().IsQuasiElastic() && !is_strange && !is_charm;
190
  bool is_MEC      = event->Summary()->ProcInfo().IsMEC();
191
  bool is_Coherent = event->Summary()->ProcInfo().IsCoherent();
192
  bool is_NuEEl    = event->Summary()->ProcInfo().IsNuElectronElastic();
193
  bool is_IMD      = event->Summary()->ProcInfo().IsInverseMuDecay();
194
  bool is_IMDanh   = event->Summary()->ProcInfo().IsIMDAnnihilation();
195

    
196
  /*  if (is_NuEEl) {
197
    LOG("UtlEventRec",pWARN)<<"Ignoring NuEEl event.";
198
    return conv::kmUndefined;
199
    }*/
200
  
201
  vector<int> primary_hadr_system;
202
  vector<int> final_hadr_system;
203
  
204
  GetHadronicSystems(event,primary_hadr_system,final_hadr_system);
205
  
206
  // Set up variables for counting particles
207
  int nP        = 0;
208
  int nN        = 0;
209
  int nPip      = 0;
210
  int nPim      = 0;
211
  int nPi0      = 0;
212
  int nKp       = 0;
213
  int nKm       = 0;
214
  int nK0       = 0;
215
  int nEM       = 0;
216
  int nOtherHad = 0;
217
  int nNucleus  = 0;
218
  
219
  std::vector<const genie::GHepParticle*> unknown_particles;
220
  
221
  const vector<int> & hadr_system = (use_final_state ? final_hadr_system : primary_hadr_system);
222
  // We exclude states with the right number of pions but other mesons.
223
  
224

    
225
  
226
  for (unsigned int ip = 0 ; ip < hadr_system.size(); ip++) {
227
    const genie::GHepParticle * p = event->Particle(hadr_system[ip]);
228
    int pdg = p->Pdg();
229
    if      (pdg == genie::kPdgProton  || pdg == genie::kPdgAntiProton)                           { nP++;       }
230
    else if (pdg == genie::kPdgNeutron || pdg == genie::kPdgAntiNeutron)                          { nN++;       }
231
    else if (pdg == genie::kPdgPiP)                                                               { nPip++;     }
232
    else if (pdg == genie::kPdgPiM)                                                               { nPim++;     }
233
    else if (pdg == genie::kPdgPi0)                                                               { nPi0++;     }
234
    else if (pdg == genie::kPdgKP)                                                                { nKp++;      }
235
    else if (pdg == genie::kPdgKM)                                                                { nKm++;      }
236
    else if (pdg == genie::kPdgK0    || pdg == genie::kPdgAntiK0)                                 { nK0++;      }
237
    else if (pdg == genie::kPdgGamma || pdg == genie::kPdgElectron || pdg == genie::kPdgPositron) { nEM++;      }
238
    else if (genie::pdg::IsHadron(pdg))                                                           { nOtherHad++;}
239
    else if (genie::pdg::IsIon   (pdg))                                                           { nNucleus++; }
240
    else { unknown_particles.push_back(p); }
241
    
242
  }
243
  
244
  if (unknown_particles.size()) {
245
    std::ostringstream s;
246
    for (unsigned int ip = 0 ; ip < unknown_particles.size() ; ip++) {
247
      const genie::GHepParticle * p = unknown_particles[ip];
248
      s<<*p<<std::endl;
249
    }
250
    LOG("UtlEventRec",pWARN)<<unknown_particles.size()<<" unidentified particles in event:\n"<<s.str();
251
  }
252
  
253
  int nPiC = nPip+nPim;
254
  int nK   = nKp+nKm+nK0;
255
  
256
  bool is_1piC     = (!is_Coherent) && (nPiC==1) && (nPi0==0) && (nK==0) && (nOtherHad==0);
257
  bool is_1pi0     = (!is_Coherent) && (nPiC==0) && (nPi0==1) && (nK==0) && (nOtherHad==0); 
258
  bool is_2piC     = (!is_Coherent) && (nPiC==2) && (nPi0==0) && (nK==0) && (nOtherHad==0); 
259
  bool is_2pi0     = (!is_Coherent) && (nPiC==0) && (nPi0==2) && (nK==0) && (nOtherHad==0); 
260
  bool is_1piC1pi0 = (!is_Coherent) && (nPiC==1) && (nPi0==1) && (nK==0) && (nOtherHad==0); 
261
  
262
  bool is_ccnuel   = (is_NuEEl && is_CC);
263
  bool is_ncnuel   = (is_NuEEl && is_NC);
264
  
265
  vector<conv::ReacMode_t> modes;
266
  
267
  // First do exclusive modes
268
  if(is_numu) {
269
    if(is_CC) {
270
      if (is_QE)       { modes.push_back(conv::kmNuMuCCQE    ); } // numu    CC QE
271
      if (is_MEC)      { modes.push_back(conv::kmNuMuCCMEC   ); } // numu    CC MEC
272
      if (is_1piC)     { modes.push_back(conv::kmNuMuCC1piC  ); } // numu    CC 1 pi+/-
273
      if (is_1pi0)     { modes.push_back(conv::kmNuMuCC1pi0  ); } // numu    CC 1 pi0
274
      if (is_2piC)     { modes.push_back(conv::kmNuMuCC2piC  ); } // numu    CC 2 pi+/-
275
      if (is_2pi0)     { modes.push_back(conv::kmNuMuCC2pi0  ); } // numu    CC 2 pi0
276
      if (is_1piC1pi0) { modes.push_back(conv::kmNuMuCCpiCpi0); } // numu    CC 1 pi0 + 1 pi+/-
277
      if (is_Coherent) { modes.push_back(conv::kmNuMuCCcoh   ); } // numu    CC coherent
278
      if (is_IMD)      { modes.push_back(conv::kmNuMuCCIMD   ); } // numu    CC IMD
279
    }
280
    if(is_NC) {
281
      if (is_1piC)     { modes.push_back(conv::kmNuMuNC1piC  ); } // numu    NC 1 pi+/-
282
      if (is_1pi0)     { modes.push_back(conv::kmNuMuNC1pi0  ); } // numu    NC 1 pi0
283
      if (is_Coherent) { modes.push_back(conv::kmNuMuNCcoh   ); } // numu    NC coherent
284
    }
285
    if (is_ccnuel) { modes.push_back(conv::kmNuMuCCELScat   ); } // numu + electron CC
286
    if (is_ncnuel) { modes.push_back(conv::kmNuMuNCELScat   ); } // numu + electron NC
287
  }
288

    
289
  if(is_numubar) {
290
    if(is_CC) {
291
      if (is_QE)       { modes.push_back(conv::kmNuMuBarCCQE    ); } // numubar CC QE
292
      if (is_MEC)      { modes.push_back(conv::kmNuMuBarCCMEC   ); } // numubar CC MEC
293
      if (is_1piC)     { modes.push_back(conv::kmNuMuBarCC1piC  ); } // numubar CC 1 pi+/-
294
      if (is_1pi0)     { modes.push_back(conv::kmNuMuBarCC1pi0  ); } // numubar CC 1 pi0
295
      if (is_2piC)     { modes.push_back(conv::kmNuMuBarCC2piC  ); } // numubar CC 2 pi+/-
296
      if (is_2pi0)     { modes.push_back(conv::kmNuMuBarCC2pi0  ); } // numubar CC 2 pi0
297
      if (is_1piC1pi0) { modes.push_back(conv::kmNuMuBarCCpiCpi0); } // numubar CC 1 pi0 + 1 pi+/-
298
      if (is_Coherent) { modes.push_back(conv::kmNuMuBarCCcoh   ); } // numubar CC coherent
299
    }
300
    if(is_NC) {
301
      if (is_1piC)     { modes.push_back(conv::kmNuMuBarNC1piC ); } // numubar NC 1 pi+/-
302
      if (is_1pi0)     { modes.push_back(conv::kmNuMuBarNC1pi0 ); } // numubar NC 1 pi0
303
      if (is_Coherent) { modes.push_back(conv::kmNuMuBarNCcoh  ); } // numubar NC coherent
304
    }
305
    if (is_ccnuel) { modes.push_back(conv::kmNuMuBarCCELScat   ); } // numubar + electron CC                                                 
306
    if (is_ncnuel) { modes.push_back(conv::kmNuMuBarNCELScat   ); } // numubar + electron NC    
307
  }
308

    
309
  //
310
  if(is_nue) {
311
    if(is_CC) {
312
      if (is_QE)       { modes.push_back(conv::kmNuECCQE       );}     // nue     CC QE
313
      if (is_MEC)      { modes.push_back(conv::kmNuECCMEC      );}     // nue     CC MEC
314
      if (is_1piC)     { modes.push_back(conv::kmNuECC1piC     );}     // nue     CC 1 pi+/-
315
      if (is_1pi0)     { modes.push_back(conv::kmNuECC1pi0     );}     // nue     CC 1 pi0
316
      if (is_2piC)     { modes.push_back(conv::kmNuECC2piC     );}     // nue     CC 2 pi+/-
317
      if (is_2pi0)     { modes.push_back(conv::kmNuECC2pi0     );}     // nue     CC 2 pi0
318
      if (is_1piC1pi0) { modes.push_back(conv::kmNuECCpiCpi0   );}     // nue     CC 1 pi0 + 1 pi+/-
319
      if (is_Coherent) { modes.push_back(conv::kmNuECCcoh      );}     // nue     CC coherent
320
    }
321
    if(is_NC) {
322
      if (is_1piC)     { modes.push_back(conv::kmNuENC1piC     );}     // nue     NC 1 pi+/-
323
      if (is_1pi0)     { modes.push_back(conv::kmNuENC1pi0     );}     // nue     NC 1 pi0
324
      if (is_Coherent) { modes.push_back(conv::kmNuENCcoh      );}     // nue     NC coherent
325
    }
326
    if(is_NuEEl) { modes.push_back(conv::kmNuEELScat     );}     // nue + electron scattering (NC + CC + interference)
327
  }
328
  
329
  //
330
  if(is_nuebar) {
331
    if(is_CC) {
332
      if (is_QE)       { modes.push_back(conv::kmNuEBarCCQE    );}     // nuebar  CC QE
333
      if (is_MEC)      { modes.push_back(conv::kmNuEBarCCMEC   );}     // nuebar  CC MEC
334
      if (is_1piC)     { modes.push_back(conv::kmNuEBarCC1piC  );}     // nuebar  CC 1 pi+/-
335
      if (is_1pi0)     { modes.push_back(conv::kmNuEBarCC1pi0  );}     // nuebar  CC 1 pi0
336
      if (is_2piC)     { modes.push_back(conv::kmNuEBarCC2piC  );}     // nuebar  CC 2 pi+/-
337
      if (is_2pi0)     { modes.push_back(conv::kmNuEBarCC2pi0  );}     // nuebar  CC 2 pi0
338
      if (is_1piC1pi0) { modes.push_back(conv::kmNuEBarCCpiCpi0);}     // nuebar  CC 1 pi0 + 1 pi+/-
339
      if (is_Coherent) { modes.push_back(conv::kmNuEBarCCcoh   );}     // nuebar  CC coherent
340
      if (is_IMDanh)   { modes.push_back(conv::kmNuEBarCCIMD   );}     // nuebar  CC IMD annihilation 
341
    }
342
    if(is_NC) {
343
      if (is_1piC)     { modes.push_back(conv::kmNuEBarNC1piC );}     // nuebar  NC 1 pi+/-
344
      if (is_1pi0)     { modes.push_back(conv::kmNuEBarNC1pi0 );}     // nuebar  NC 1 pi0
345
      if (is_Coherent) { modes.push_back(conv::kmNuEBarNCcoh  );}     // nuebar  NC coherent
346
    }
347
    if(is_NuEEl) { modes.push_back(conv::kmNuEBarELScat     );}     // nuebar + electron scattering (NC + CC + interference)
348
  }
349
  
350
  // If we've got one exclusive mode, we did well.
351
  if (modes.size() == 1) {
352
    return modes[0];
353
  }
354
  // Otherwise look for inclusive
355
  else if (modes.size() == 0) {
356
    if (is_numu) {
357
      if (is_CC) {return conv::kmNuMuCCoth;}
358
      if (is_NC) {return conv::kmNuMuNCoth;}
359
    }
360
    if (is_numubar) {
361
      if (is_CC) {return conv::kmNuMuBarCCoth;}
362
      if (is_NC) {return conv::kmNuMuBarNCoth;}
363
    }
364
    if (is_nue) {
365
      if (is_CC) {return conv::kmNuECCoth;}
366
      if (is_NC) {return conv::kmNuENCoth;}
367
    }
368
    if (is_nuebar) {
369
      if (is_CC) {return conv::kmNuEBarCCoth;}
370
      if (is_NC) {return conv::kmNuEBarNCoth;}
371
    }
372
  }
373
  else {
374
    std::ostringstream mode_names;
375
    for (unsigned int ii = 0 ; ii < modes.size(); ii++) {
376
      mode_names << "\n "<<util::AsString(modes[ii]);
377
    }
378
    LOG("UtlEventRec",pERROR)<<"Event filed into multiple exclusive modes: "
379
      <<mode_names.str()
380
      <<"\nFull record: "<<*event;
381
      util::Die(-1);
382
  }
383

    
384
  return conv::kmUndefined;
385
}
386
//_________________________________________________________________________________