35 #include <TObjString.h> 36 #include <TLorentzVector.h> 40 using std::ostringstream;
46 bool using_new_version =
false;
47 bool event_printout =
false;
92 const int kStdHepIdxPx = 0;
93 const int kStdHepIdxPy = 1;
94 const int kStdHepIdxPz = 2;
95 const int kStdHepIdxE = 3;
101 TFile
file(filename,
"READ");
102 TTree *
tree = (TTree *) file.Get(
"gRooTracker");
110 TObjString* EvtCode = 0;
121 double StdHepX4 [
kNPmax][4];
122 double StdHepP4 [
kNPmax][4];
123 double StdHepPolz [
kNPmax][3];
131 double NuParentDecP4 [4];
132 double NuParentDecX4 [4];
133 double NuParentProP4 [4];
134 double NuParentProX4 [4];
141 TBranch * brEvtFlags = tree ->
GetBranch (
"EvtFlags");
142 TBranch * brEvtCode = tree ->
GetBranch (
"EvtCode");
143 TBranch * brEvtNum = tree ->
GetBranch (
"EvtNum");
144 TBranch * brEvtXSec = tree ->
GetBranch (
"EvtXSec");
145 TBranch * brEvtDXSec = tree ->
GetBranch (
"EvtDXSec");
146 TBranch * brEvtWght = tree ->
GetBranch (
"EvtWght");
147 TBranch * brEvtProb = tree ->
GetBranch (
"EvtProb");
148 TBranch * brEvtVtx = tree ->
GetBranch (
"EvtVtx");
149 TBranch * brStdHepN = tree ->
GetBranch (
"StdHepN");
150 TBranch * brStdHepPdg = tree ->
GetBranch (
"StdHepPdg");
151 TBranch * brStdHepStatus = tree ->
GetBranch (
"StdHepStatus");
152 TBranch * brStdHepRescat = (using_new_version) ? tree ->
GetBranch (
"StdHepRescat") : 0;
153 TBranch * brStdHepX4 = tree ->
GetBranch (
"StdHepX4");
154 TBranch * brStdHepP4 = tree ->
GetBranch (
"StdHepP4");
155 TBranch * brStdHepPolz = tree ->
GetBranch (
"StdHepPolz");
156 TBranch * brStdHepFd = tree ->
GetBranch (
"StdHepFd");
157 TBranch * brStdHepLd = tree ->
GetBranch (
"StdHepLd");
158 TBranch * brStdHepFm = tree ->
GetBranch (
"StdHepFm");
159 TBranch * brStdHepLm = tree ->
GetBranch (
"StdHepLm");
160 TBranch * brG2NeutEvtCode = (using_new_version) ? tree ->
GetBranch (
"G2NeutEvtCode") : 0;
161 TBranch * brNuParentPdg = tree ->
GetBranch (
"NuParentPdg");
162 TBranch * brNuParentDecMode = tree ->
GetBranch (
"NuParentDecMode");
163 TBranch * brNuParentDecP4 = tree ->
GetBranch (
"NuParentDecP4");
164 TBranch * brNuParentDecX4 = tree ->
GetBranch (
"NuParentDecX4");
165 TBranch * brNuParentProP4 = tree ->
GetBranch (
"NuParentProP4");
166 TBranch * brNuParentProX4 = tree ->
GetBranch (
"NuParentProX4");
167 TBranch * brNuParentProNVtx = tree ->
GetBranch (
"NuParentProNVtx");
169 brEvtFlags -> SetAddress ( &EvtFlags );
170 brEvtCode -> SetAddress ( &EvtCode );
171 brEvtNum -> SetAddress ( &EvtNum );
172 brEvtXSec -> SetAddress ( &EvtXSec );
173 brEvtDXSec -> SetAddress ( &EvtDXSec );
174 brEvtWght -> SetAddress ( &EvtWght );
175 brEvtProb -> SetAddress ( &EvtProb );
176 brEvtVtx -> SetAddress ( EvtVtx );
177 brStdHepN -> SetAddress ( &StdHepN );
178 brStdHepPdg -> SetAddress ( StdHepPdg );
179 brStdHepStatus -> SetAddress ( StdHepStatus );
180 if(using_new_version) {
181 brStdHepRescat -> SetAddress ( StdHepRescat );
183 brStdHepX4 -> SetAddress ( StdHepX4 );
184 brStdHepP4 -> SetAddress ( StdHepP4 );
185 brStdHepPolz -> SetAddress ( StdHepPolz );
186 brStdHepFd -> SetAddress ( StdHepFd );
187 brStdHepLd -> SetAddress ( StdHepLd );
188 brStdHepFm -> SetAddress ( StdHepFm );
189 brStdHepLm -> SetAddress ( StdHepLm );
190 if(using_new_version) {
191 brG2NeutEvtCode -> SetAddress ( &G2NeutEvtCode );
193 brNuParentPdg -> SetAddress ( &NuParentPdg );
194 brNuParentDecMode -> SetAddress ( &NuParentDecMode );
195 brNuParentDecP4 -> SetAddress ( NuParentDecP4 );
196 brNuParentDecX4 -> SetAddress ( NuParentDecX4 );
197 brNuParentProP4 -> SetAddress ( NuParentProP4 );
198 brNuParentProX4 -> SetAddress ( NuParentProX4 );
199 brNuParentProNVtx -> SetAddress ( &NuParentProNVtx );
206 outfilename << filename <<
".reaction_codes";
208 ofstream
outfile(outfilename.str().c_str());
210 outfile <<
"# NEUT reaction code for GENIE file: " << filename <<
endl;
219 int n = tree->GetEntries();
220 printf(
"Number of entries: %d", n);
222 for(
int i=0;
i < tree->GetEntries();
i++) {
228 printf(
"\n\n ** Current entry: %d \n",
i);
235 printf(
"\n -----------------------------------------------------------------------------------------------------------------");
236 printf(
"\n Event code : %s", EvtCode->String().Data());
237 printf(
"\n Event x-section : %10.5f * 1E-38* cm^2", EvtXSec);
238 printf(
"\n Event kinematics x-section : %10.5f * 1E-38 * cm^2/{K^n}", EvtDXSec);
239 printf(
"\n Event weight : %10.8f", EvtWght);
240 printf(
"\n Event vertex : x = %8.2f mm, y = %8.2f mm, z = %8.2f mm", EvtVtx[0], EvtVtx[1], EvtVtx[2]);
241 printf(
"\n * Particle list:");
242 printf(
"\n --------------------------------------------------------------------------------------------------------------------------");
243 printf(
"\n | Idx | Ist | PDG | Rescat | Mother | Daughter | Px | Py | Pz | E | x | y | z |");
244 printf(
"\n | | | | | | |(GeV/c) |(GeV/c) |(GeV/c) | (GeV) | (fm) | (fm) | (fm) |");
245 printf(
"\n --------------------------------------------------------------------------------------------------------------------------");
247 for(
int ip=0;
ip<StdHepN;
ip++) {
248 printf(
"\n | %3d | %3d | %10d | %6d | %3d | %3d | %3d | %3d | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f | %6.3f |",
249 ip, StdHepStatus[
ip], StdHepPdg[ip], StdHepRescat[ip],
250 StdHepFm[ip], StdHepLm[ip], StdHepFd[ip], StdHepLd[ip],
251 StdHepP4[ip][0], StdHepP4[ip][1], StdHepP4[ip][2], StdHepP4[ip][3],
252 StdHepX4[ip][0], StdHepX4[ip][1], StdHepX4[ip][2]);
254 printf(
"\n --------------------------------------------------------------------------------------------------------------------------");
255 printf(
"\n * Flux Info:");
256 printf(
"\n Parent hadron pdg code : %d", NuParentPdg);
257 printf(
"\n Parent hadron decay mode : %d", NuParentDecMode);
258 printf(
"\n Parent hadron 4p at decay : Px = %6.3f GeV/c, Py = %6.3f GeV/c, Pz = %6.3f GeV/c, E = %6.3f GeV",
259 NuParentDecP4[0], NuParentDecP4[1], NuParentDecP4[2], NuParentDecP4[3]);
260 printf(
"\n Parent hadron 4p at prod. : Px = %6.3f GeV/c, Py = %6.3f GeV/c, Pz = %6.3f GeV/c, E = %6.3f GeV",
261 NuParentProP4[0], NuParentProP4[1], NuParentProP4[2], NuParentProP4[3]);
262 printf(
"\n Parent hadron 4x at decay : x = %6.3f m, y = %6.3f m, z = %6.3f m, t = %6.3f s",
263 NuParentDecX4[0], NuParentDecX4[1], NuParentDecX4[2], NuParentDecX4[3]);
264 printf(
"\n Parent hadron 4x at prod. : x = %6.3f m, y = %6.3f m, z = %6.3f m, t = %6.3f s",
265 NuParentProX4[0], NuParentProX4[1], NuParentProX4[2], NuParentProX4[3]);
266 printf(
"\n -------------------------------------------------------------------------------------------------------------------------- \n");
276 string genie_code = EvtCode->GetString().Data();
277 bool is_cc = (genie_code.find(
"Weak[CC]") != string::npos);
278 bool is_nc = (genie_code.find(
"Weak[NC]") != string::npos);
279 bool is_charm = (genie_code.find(
"charm") != string::npos);
280 bool is_qel = (genie_code.find(
"QES") != string::npos);
281 bool is_dis = (genie_code.find(
"DIS") != string::npos);
282 bool is_res = (genie_code.find(
"RES") != string::npos);
283 bool is_cohpi = (genie_code.find(
"COH") != string::npos);
284 bool is_ve = (genie_code.find(
"NuEEL") != string::npos);
285 bool is_imd = (genie_code.find(
"IMD") != string::npos);
289 int fsl_pos = StdHepFd[inu_pos];
294 int nu_code = StdHepPdg[inu_pos];
295 bool is_nu = (nu_code == kPdgNuE || nu_code == kPdgNuMu || nu_code ==
kPdgNuTau );
296 bool is_nubar = (nu_code == kPdgAntiNuE || nu_code == kPdgAntiNuMu || nu_code ==
kPdgAntiNuTau);
298 int target_code = StdHepPdg[tgt_pos];
299 bool nuclear_target = (target_code > 1000000000);
301 bool hitnuc_set =
false;
304 for(
int ip = 1;
ip <= 2;
ip++)
306 int ghep_pdgc = StdHepPdg[
ip];
307 if(ghep_pdgc == kPdgProton) {
313 if(ghep_pdgc == kPdgNeutron) {
327 TLorentzVector p4v ( StdHepP4 [inu_pos] [kStdHepIdxPx],
328 StdHepP4 [inu_pos] [kStdHepIdxPy],
329 StdHepP4 [inu_pos] [kStdHepIdxPz],
330 StdHepP4 [inu_pos] [kStdHepIdxE] );
331 TLorentzVector p4l ( StdHepP4 [fsl_pos] [kStdHepIdxPx],
332 StdHepP4 [fsl_pos] [kStdHepIdxPy],
333 StdHepP4 [fsl_pos] [kStdHepIdxPz],
334 StdHepP4 [fsl_pos] [kStdHepIdxE] );
335 TLorentzVector p4n ( StdHepP4 [hitnuc_pos] [kStdHepIdxPx],
336 StdHepP4 [hitnuc_pos] [kStdHepIdxPy],
337 StdHepP4 [hitnuc_pos] [kStdHepIdxPz],
338 StdHepP4 [hitnuc_pos] [kStdHepIdxE] );
340 TLorentzVector q = p4v - p4l;
341 TLorentzVector
w = p4n + q;
349 if (is_qel && !is_charm && is_cc && is_nu ) evtype = 1;
350 else if (is_qel && !is_charm && is_nc && is_nu && is_p ) evtype = 51;
351 else if (is_qel && !is_charm && is_nc && is_nu && is_n ) evtype = 52;
352 else if (is_qel && !is_charm && is_cc && is_nubar ) evtype = -1;
353 else if (is_qel && !is_charm && is_nc && is_nubar && is_p) evtype = -51;
354 else if (is_qel && !is_charm && is_nc && is_nubar && is_n) evtype = -52;
358 else if (is_qel && is_charm && is_cc && is_nu ) evtype = 25;
359 else if (is_qel && is_charm && is_cc && is_nubar ) evtype = -25;
363 else if ( is_imd ) evtype = 9;
364 else if ( is_ve ) evtype = 59;
368 else if (is_cohpi && is_cc && is_nu ) evtype = 16;
369 else if (is_cohpi && is_cc && is_nubar) evtype = -16;
370 else if (is_cohpi && is_nc && is_nu ) evtype = 36;
371 else if (is_cohpi && is_nc && is_nubar) evtype = -36;
376 else if (is_dis && W_gt_2 && is_cc && is_nu ) evtype = 26;
377 else if (is_dis && W_gt_2 && is_nc && is_nu ) evtype = 46;
378 else if (is_dis && W_gt_2 && is_cc && is_nubar) evtype = -26;
379 else if (is_dis && W_gt_2 && is_nc && is_nubar) evtype = -46;
383 else if ( is_res || (is_dis && !W_gt_2) ) {
390 int nn=0, np=0, npi0=0, npip=0, npim=0, nKp=0, nKm=0, nK0=0, neta=0, nlambda=0, ngamma=0;
392 for(
int ip = 0;
ip < StdHepN;
ip++)
394 int ghep_ist = StdHepStatus[
ip];
395 int ghep_pdgc = StdHepPdg[
ip];
396 int ghep_fm = StdHepFm[
ip];
397 int ghep_fmpdgc = (ghep_fm==-1) ? 0 : StdHepPdg[ghep_fm];
404 bool decayed = (ghep_ist==kIStDecayedState && (ghep_pdgc==kPdgPi0 || ghep_pdgc==
kPdgEta));
405 bool parent_included = (ghep_fmpdgc==kPdgPi0 || ghep_fmpdgc==
kPdgEta);
409 (!nuclear_target && decayed) ||
410 (!nuclear_target && ghep_ist==kIStStableFinalState && !parent_included);
412 if(!count_it)
continue;
414 if(ghep_pdgc == kPdgProton ) np++;
415 if(ghep_pdgc == kPdgNeutron) nn++;
416 if(ghep_pdgc == kPdgPiP) npip++;
417 if(ghep_pdgc == kPdgPiM) npim++;
418 if(ghep_pdgc == kPdgPi0) npi0++;
419 if(ghep_pdgc == kPdgEta) neta++;
420 if(ghep_pdgc == kPdgKP) nKp++;
421 if(ghep_pdgc == kPdgKM) nKm++;
422 if(ghep_pdgc == kPdgK0) nK0++;
423 if(ghep_pdgc == kPdgAntiK0) nK0++;
424 if(ghep_pdgc == kPdgLambda) nlambda++;
425 if(ghep_pdgc == kPdgAntiLambda) nlambda++;
426 if(ghep_pdgc == kPdgGamma) ngamma++;
429 cout <<
"Num of primary particles: \n p = " << np <<
", n = " << nn
430 <<
", pi+ = " << npip <<
", pi- = " << npim <<
", pi0 = " << npi0
431 <<
", eta = " << neta
432 <<
", K+ = " << nKp <<
", K- = " << nKm <<
", K0 = " << nK0
433 <<
", Lambda's = " << nlambda
434 <<
", gamma's = " << ngamma
438 int npi = npi0 + npip + npim;
439 int nK = nK0 + nKp + nKm;
440 int neKL = neta + nK + nlambda;
442 bool is_radiative_dec = (nnuc==1) && (npi==0) && (ngamma==1);
448 if (is_res && is_nu && is_cc && is_n && is_radiative_dec) evtype = 17;
449 else if (is_res && is_nu && is_nc && is_n && is_radiative_dec) evtype = 38;
450 else if (is_res && is_nu && is_nc && is_p && is_radiative_dec) evtype = 39;
452 else if (is_res && is_nubar && is_cc && is_p && is_radiative_dec) evtype = -17;
453 else if (is_res && is_nubar && is_nc && is_n && is_radiative_dec) evtype = -38;
454 else if (is_res && is_nubar && is_nc && is_p && is_radiative_dec) evtype = -39;
461 else if (is_nu && is_cc && is_p && np==1 && nn==0 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 11;
462 else if (is_nu && is_cc && is_n && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 12;
463 else if (is_nu && is_cc && is_n && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 13;
466 else if (is_nu && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 31;
467 else if (is_nu && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = 32;
468 else if (is_nu && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = 33;
469 else if (is_nu && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = 34;
472 else if (is_nubar && is_cc && is_n && np==0 && nn==1 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -11;
473 else if (is_nubar && is_cc && is_p && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -12;
474 else if (is_nubar && is_cc && is_p && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -13;
477 else if (is_nubar && is_nc && is_n && np==0 && nn==1 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -31;
478 else if (is_nubar && is_nc && is_p && np==1 && nn==0 && npip==0 && npim==0 && npi0==1 && neKL==0) evtype = -32;
479 else if (is_nubar && is_nc && is_n && np==1 && nn==0 && npip==0 && npim==1 && npi0==0 && neKL==0) evtype = -33;
480 else if (is_nubar && is_nc && is_p && np==0 && nn==1 && npip==1 && npim==0 && npi0==0 && neKL==0) evtype = -34;
486 else if (is_res && is_nu && is_cc && is_n && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 22;
487 else if (is_res && is_nu && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 42;
488 else if (is_res && is_nu && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = 43;
490 else if (is_res && is_nubar && is_cc && is_p && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -22;
491 else if (is_res && is_nubar && is_nc && is_n && np==0 && nn==1 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -42;
492 else if (is_res && is_nubar && is_nc && is_p && np==1 && nn==0 && npi==0 && nK==0 && nlambda==0 && neta==1) evtype = -43;
498 else if (is_res && is_nu && is_cc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 23;
499 else if (is_res && is_nu && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 44;
500 else if (is_res && is_nu && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = 45;
502 else if (is_res && is_nubar && is_cc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -23;
503 else if (is_res && is_nubar && is_nc && is_n && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -44;
504 else if (is_res && is_nubar && is_nc && is_p && nnuc==0 && npi==0 && nK==1 && nlambda==1 && neta==0) evtype = -45;
510 else if (is_nu && is_cc && npi>1) evtype = 21;
511 else if (is_nu && is_nc && npi>1) evtype = 41;
512 else if (is_nubar && is_cc && npi>1) evtype = -21;
513 else if (is_nubar && is_nc && npi>1) evtype = -41;
522 cout <<
"Rare RES/low-W DIS final state: Bundled-in with multi-pi events" <<
endl;
524 if (is_nu && is_cc) evtype = 21;
525 else if (is_nu && is_nc) evtype = 41;
526 else if (is_nubar && is_cc) evtype = -21;
527 else if (is_nubar && is_nc) evtype = -41;
531 cout <<
" *** GENIE event = " <<
i <<
" --> NEUT reaction code = " << evtype <<
endl;
532 if(using_new_version) {
535 cout <<
"NEUT reaction code stored at the rootracker file = " << G2NeutEvtCode <<
endl;
536 assert(evtype == G2NeutEvtCode);
ntuple GetBranch("organID") -> SetAddress(&xx)
string outfilename
knobs that need extra care
printf("%d Experimental points found\n", nlines)
void neut_code_from_rootracker(const char *filename)
assert(nhit_max >=nhit_nbins)