34 LOG(
"GHepUtils",
pWARN) <<
"Null event!";
40 Interaction * interaction =
event->Summary();
42 const ProcessInfo & proc = interaction->ProcInfo();
43 const InitialState & init = interaction->InitState();
44 const XclsTag & xcls = interaction->ExclTag();
45 const Kinematics & kine = interaction->Kine();
46 const Target & tgt = init.Tgt();
48 bool is_cc = proc.IsWeakCC();
49 bool is_nc = proc.IsWeakNC();
50 bool is_charm = xcls.IsCharmEvent();
51 bool is_qel = proc.IsQuasiElastic();
52 bool is_dis = proc.IsDeepInelastic();
53 bool is_res = proc.IsResonant();
54 bool is_coh_pr = proc.IsCoherentProduction();
55 bool is_ve = proc.IsNuElectronElastic();
56 bool is_mec = proc.IsMEC();
57 bool is_imd = proc.IsInverseMuDecay();
58 bool is_ask = proc.IsSingleKaon();
59 bool is_diff = proc.IsDiffractive();
60 bool is_p = tgt.HitNucIsSet() ? tgt.HitNucPdg()==
kPdgProton :
false;
61 bool is_n = tgt.HitNucIsSet() ? tgt.HitNucPdg()==
kPdgNeutron :
false;
64 bool W_gt_2 = kine.KVSet(
kKVW) ? (kine.W() > 2.0) :
false;
68 if (is_qel && !is_charm && is_cc && is_nu ) evtype = 1;
69 else if (is_qel && !is_charm && is_nc && is_nu && is_p ) evtype = 51;
70 else if (is_qel && !is_charm && is_nc && is_nu && is_n ) evtype = 52;
71 else if (is_qel && !is_charm && is_cc && is_nubar ) evtype = -1;
72 else if (is_qel && !is_charm && is_nc && is_nubar && is_p) evtype = -51;
73 else if (is_qel && !is_charm && is_nc && is_nubar && is_n) evtype = -52;
76 else if (is_mec && !is_charm && is_cc && is_nu ) evtype = 2;
77 else if (is_mec && !is_charm && is_cc && is_nubar ) evtype = -2;
81 else if (is_qel && is_charm && is_cc && is_nu ) evtype = 26;
82 else if (is_qel && is_charm && is_cc && is_nubar ) evtype = -26;
86 else if ( is_imd ) evtype = 9;
87 else if ( is_ve ) evtype = 59;
91 else if (is_coh_pr && is_cc && is_nu ) evtype = 16;
92 else if (is_coh_pr && is_cc && is_nubar) evtype = -16;
93 else if (is_coh_pr && is_nc && is_nu ) evtype = 36;
94 else if (is_coh_pr && is_nc && is_nubar) evtype = -36;
99 else if (is_dis && W_gt_2 && is_cc && is_nu ) evtype = 26;
100 else if (is_dis && W_gt_2 && is_nc && is_nu ) evtype = 46;
101 else if (is_dis && W_gt_2 && is_cc && is_nubar) evtype = -26;
102 else if (is_dis && W_gt_2 && is_nc && is_nubar) evtype = -46;
106 else if ( is_res || (is_dis && !W_gt_2) || is_ask ) {
108 LOG(
"GHepUtils",
pNOTICE) <<
"Current event is RES or DIS with W<2";
113 int nn=0, np=0, npi0=0, npip=0, npim=0, nKp=0, nKm=0, nK0=0, neta=0, nlambda=0, ngamma=0;
114 bool nuclear_target = init.Tgt().IsNucleus();
116 TIter event_iter(event);
117 GHepParticle * p = 0;
119 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
122 int ghep_pdgc = p->Pdg();
123 int ghep_fm = p->FirstMother();
124 int ghep_fmpdgc = (ghep_fm==-1) ? 0 : event->Particle(ghep_fm)->Pdg();
132 bool parent_included = (ghep_fmpdgc==
kPdgPi0 || ghep_fmpdgc==
kPdgEta);
136 (!nuclear_target && decayed) ||
139 if(!count_it)
continue;
143 if(ghep_pdgc ==
kPdgPiP) npip++;
144 if(ghep_pdgc ==
kPdgPiM) npim++;
145 if(ghep_pdgc ==
kPdgPi0) npi0++;
146 if(ghep_pdgc ==
kPdgEta) neta++;
147 if(ghep_pdgc ==
kPdgKP) nKp++;
148 if(ghep_pdgc ==
kPdgKM) nKm++;
149 if(ghep_pdgc ==
kPdgK0) nK0++;
156 <<
"Num of primary particles: \n p = " << np <<
", n = " << nn
157 <<
", pi+ = " << npip <<
", pi- = " << npim <<
", pi0 = " << npi0
158 <<
", eta = " << neta
159 <<
", K+ = " << nKp <<
", K- = " << nKm <<
", K0 = " << nK0
160 <<
", Lambda's = " << nlambda
161 <<
", gamma's = " << ngamma;
164 int npi = npi0 + npip + npim;
165 int nK = nK0 + nKp + nKm;
166 int neKL = neta + nK + nlambda;
168 bool is_radiative_dec = (nnuc==1) && (npi==0) && (ngamma==1);
174 if (is_res && is_nu && is_cc && is_n && is_radiative_dec) evtype = 17;
175 else if (is_res && is_nu && is_nc && is_n && is_radiative_dec) evtype = 38;
176 else if (is_res && is_nu && is_nc && is_p && is_radiative_dec) evtype = 39;
178 else if (is_res && is_nubar && is_cc && is_p && is_radiative_dec) evtype = -17;
179 else if (is_res && is_nubar && is_nc && is_n && is_radiative_dec) evtype = -38;
180 else if (is_res && is_nubar && is_nc && is_p && is_radiative_dec) evtype = -39;
187 else if (is_nu && is_cc && is_p && np==1 && nn==0 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 11;
188 else if (is_nu && is_cc && is_n && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 12;
189 else if (is_nu && is_cc && is_n && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 13;
192 else if (is_nu && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 31;
193 else if (is_nu && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 32;
194 else if (is_nu && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = 33;
195 else if (is_nu && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 34;
198 else if (is_nubar && is_cc && is_n && np==0 && nn==1 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -11;
199 else if (is_nubar && is_cc && is_p && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -12;
200 else if (is_nubar && is_cc && is_p && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -13;
203 else if (is_nubar && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -31;
204 else if (is_nubar && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -32;
205 else if (is_nubar && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -33;
206 else if (is_nubar && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = -34;
212 else if (is_res && is_nu && is_cc && is_n && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 22;
213 else if (is_res && is_nu && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 42;
214 else if (is_res && is_nu && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 43;
216 else if (is_res && is_nubar && is_cc && is_p && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -22;
217 else if (is_res && is_nubar && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -42;
218 else if (is_res && is_nubar && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -43;
224 else if (is_res && is_nu && is_cc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 23;
225 else if (is_res && is_nu && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 44;
226 else if (is_res && is_nu && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 45;
228 else if (is_res && is_nubar && is_cc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -23;
229 else if (is_res && is_nubar && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -44;
230 else if (is_res && is_nubar && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -45;
236 else if (is_ask && is_nu && is_cc && is_n && nn==1 && np==0 && nKp==1 && neKL==1) evtype = 18;
237 else if (is_ask && is_nu && is_cc && is_n && nn==0 && np==1 && nK0==1 && neKL==1) evtype = 19;
238 else if (is_ask && is_nu && is_cc && is_p && nn==0 && np==1 && nKp==1 && neKL==1) evtype = 20;
250 else if (is_nu && is_cc && npi>1) evtype = 21;
251 else if (is_nu && is_nc && npi>1) evtype = 41;
252 else if (is_nubar && is_cc && npi>1) evtype = -21;
253 else if (is_nubar && is_nc && npi>1) evtype = -41;
262 <<
"Rare RES/low-W DIS final state: Bundled-in with multi-pi events";
264 if (is_nu && is_cc) evtype = 21;
265 else if (is_nu && is_nc) evtype = 41;
266 else if (is_nubar && is_cc) evtype = -21;
267 else if (is_nubar && is_nc) evtype = -41;
272 else if ( is_diff && is_cc ) {
273 if ( is_nu ) evtype = 15;
274 else if ( is_nubar ) evtype = -15;
276 else if ( is_diff && is_nc ) {
277 if ( is_nu ) evtype = 35;
278 else if ( is_nubar ) evtype = -35;
bool IsNeutrino(int pdgc)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
bool IsAntiNeutrino(int pdgc)
enum genie::EGHepStatus GHepStatus_t