128 #include "libxml/parser.h" 129 #include "libxml/xmlmemory.h" 136 #include <TObjString.h> 162 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 170 using std::ostringstream;
174 using std::setprecision;
177 using std::setiosflags;
180 using namespace genie;
196 int HAProbeFSI (
int,
int,
int,
double [],
int [],
int,
int,
int);
298 const double e_h = 1.3;
304 int brFSPrimLept = 0;
310 bool brFromSea =
false;
312 bool brIsQel =
false;
313 bool brIsRes =
false;
314 bool brIsDis =
false;
315 bool brIsCoh =
false;
316 bool brIsMec =
false;
317 bool brIsDfr =
false;
318 bool brIsImd =
false;
319 bool brIsSingleK =
false;
320 bool brIsImdAnh =
false;
321 bool brIsNuEL =
false;
325 bool brIsCharmPro =
false;
326 bool brIsAMNuGamma =
false;
328 int brCodeNuance = 0;
333 double brKineQ2s = 0;
407 TTree * s_tree =
new TTree(
"gst",
"GENIE Summary Event Tree");
411 s_tree->Branch(
"iev", &brIev,
"iev/I" );
412 s_tree->Branch(
"neu", &brNeutrino,
"neu/I" );
413 s_tree->Branch(
"fspl", &brFSPrimLept,
"fspl/I" );
414 s_tree->Branch(
"tgt", &brTarget,
"tgt/I" );
415 s_tree->Branch(
"Z", &brTargetZ,
"Z/I" );
416 s_tree->Branch(
"A", &brTargetA,
"A/I" );
417 s_tree->Branch(
"hitnuc", &brHitNuc,
"hitnuc/I" );
418 s_tree->Branch(
"hitqrk", &brHitQrk,
"hitqrk/I" );
419 s_tree->Branch(
"resid", &brResId,
"resid/I" );
420 s_tree->Branch(
"sea", &brFromSea,
"sea/O" );
421 s_tree->Branch(
"qel", &brIsQel,
"qel/O" );
422 s_tree->Branch(
"mec", &brIsMec,
"mec/O" );
423 s_tree->Branch(
"res", &brIsRes,
"res/O" );
424 s_tree->Branch(
"dis", &brIsDis,
"dis/O" );
425 s_tree->Branch(
"coh", &brIsCoh,
"coh/O" );
426 s_tree->Branch(
"dfr", &brIsDfr,
"dfr/O" );
427 s_tree->Branch(
"imd", &brIsImd,
"imd/O" );
428 s_tree->Branch(
"imdanh", &brIsImdAnh,
"imdanh/O" );
429 s_tree->Branch(
"singlek", &brIsSingleK,
"singlek/O" );
430 s_tree->Branch(
"nuel", &brIsNuEL,
"nuel/O" );
431 s_tree->Branch(
"em", &brIsEM,
"em/O" );
432 s_tree->Branch(
"cc", &brIsCC,
"cc/O" );
433 s_tree->Branch(
"nc", &brIsNC,
"nc/O" );
434 s_tree->Branch(
"charm", &brIsCharmPro,
"charm/O" );
435 s_tree->Branch(
"amnugamma", &brIsAMNuGamma,
"amnugamma/O" );
436 s_tree->Branch(
"neut_code", &brCodeNeut,
"neut_code/I" );
437 s_tree->Branch(
"nuance_code", &brCodeNuance,
"nuance_code/I" );
438 s_tree->Branch(
"wght", &brWeight,
"wght/D" );
439 s_tree->Branch(
"xs", &brKineXs,
"xs/D" );
440 s_tree->Branch(
"ys", &brKineYs,
"ys/D" );
441 s_tree->Branch(
"ts", &brKineTs,
"ts/D" );
442 s_tree->Branch(
"Q2s", &brKineQ2s,
"Q2s/D" );
443 s_tree->Branch(
"Ws", &brKineWs,
"Ws/D" );
444 s_tree->Branch(
"x", &brKineX,
"x/D" );
445 s_tree->Branch(
"y", &brKineY,
"y/D" );
446 s_tree->Branch(
"t", &brKineT,
"t/D" );
447 s_tree->Branch(
"Q2", &brKineQ2,
"Q2/D" );
448 s_tree->Branch(
"W", &brKineW,
"W/D" );
449 s_tree->Branch(
"EvRF", &brEvRF,
"EvRF/D" );
450 s_tree->Branch(
"Ev", &brEv,
"Ev/D" );
451 s_tree->Branch(
"pxv", &brPxv,
"pxv/D" );
452 s_tree->Branch(
"pyv", &brPyv,
"pyv/D" );
453 s_tree->Branch(
"pzv", &brPzv,
"pzv/D" );
454 s_tree->Branch(
"En", &brEn,
"En/D" );
455 s_tree->Branch(
"pxn", &brPxn,
"pxn/D" );
456 s_tree->Branch(
"pyn", &brPyn,
"pyn/D" );
457 s_tree->Branch(
"pzn", &brPzn,
"pzn/D" );
458 s_tree->Branch(
"El", &brEl,
"El/D" );
459 s_tree->Branch(
"pxl", &brPxl,
"pxl/D" );
460 s_tree->Branch(
"pyl", &brPyl,
"pyl/D" );
461 s_tree->Branch(
"pzl", &brPzl,
"pzl/D" );
462 s_tree->Branch(
"pl", &brPl,
"pl/D" );
463 s_tree->Branch(
"cthl", &brCosthl,
"cthl/D" );
464 s_tree->Branch(
"nfp", &brNfP,
"nfp/I" );
465 s_tree->Branch(
"nfn", &brNfN,
"nfn/I" );
466 s_tree->Branch(
"nfpip", &brNfPip,
"nfpip/I" );
467 s_tree->Branch(
"nfpim", &brNfPim,
"nfpim/I" );
468 s_tree->Branch(
"nfpi0", &brNfPi0,
"nfpi0/I" );
469 s_tree->Branch(
"nfkp", &brNfKp,
"nfkp/I" );
470 s_tree->Branch(
"nfkm", &brNfKm,
"nfkm/I" );
471 s_tree->Branch(
"nfk0", &brNfK0,
"nfk0/I" );
472 s_tree->Branch(
"nfem", &brNfEM,
"nfem/I" );
473 s_tree->Branch(
"nfother", &brNfOther,
"nfother/I" );
474 s_tree->Branch(
"nip", &brNiP,
"nip/I" );
475 s_tree->Branch(
"nin", &brNiN,
"nin/I" );
476 s_tree->Branch(
"nipip", &brNiPip,
"nipip/I" );
477 s_tree->Branch(
"nipim", &brNiPim,
"nipim/I" );
478 s_tree->Branch(
"nipi0", &brNiPi0,
"nipi0/I" );
479 s_tree->Branch(
"nikp", &brNiKp,
"nikp/I" );
480 s_tree->Branch(
"nikm", &brNiKm,
"nikm/I" );
481 s_tree->Branch(
"nik0", &brNiK0,
"nik0/I" );
482 s_tree->Branch(
"niem", &brNiEM,
"niem/I" );
483 s_tree->Branch(
"niother", &brNiOther,
"niother/I" );
484 s_tree->Branch(
"ni", &brNi,
"ni/I" );
485 s_tree->Branch(
"pdgi", brPdgi,
"pdgi[ni]/I" );
486 s_tree->Branch(
"resc", brResc,
"resc[ni]/I" );
487 s_tree->Branch(
"Ei", brEi,
"Ei[ni]/D" );
488 s_tree->Branch(
"pxi", brPxi,
"pxi[ni]/D" );
489 s_tree->Branch(
"pyi", brPyi,
"pyi[ni]/D" );
490 s_tree->Branch(
"pzi", brPzi,
"pzi[ni]/D" );
491 s_tree->Branch(
"nf", &brNf,
"nf/I" );
492 s_tree->Branch(
"pdgf", brPdgf,
"pdgf[nf]/I" );
493 s_tree->Branch(
"Ef", brEf,
"Ef[nf]/D" );
494 s_tree->Branch(
"pxf", brPxf,
"pxf[nf]/D" );
495 s_tree->Branch(
"pyf", brPyf,
"pyf[nf]/D" );
496 s_tree->Branch(
"pzf", brPzf,
"pzf[nf]/D" );
497 s_tree->Branch(
"pf", brPf,
"pf[nf]/D" );
498 s_tree->Branch(
"cthf", brCosthf,
"cthf[nf]/D" );
499 s_tree->Branch(
"vtxx", &brVtxX,
"vtxx/D" );
500 s_tree->Branch(
"vtxy", &brVtxY,
"vtxy/D" );
501 s_tree->Branch(
"vtxz", &brVtxZ,
"vtxz/D" );
502 s_tree->Branch(
"vtxt", &brVtxT,
"vtxt/D" );
503 s_tree->Branch(
"sumKEf", &brSumKEf,
"sumKEf/D" );
504 s_tree->Branch(
"calresp0", &brCalResp0,
"calresp0/D" );
510 er_tree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
513 LOG(
"gntpc",
pERROR) <<
"Null input GHEP event tree";
516 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
520 er_tree->SetBranchAddress(
"gmcrec", &mcrec);
522 LOG(
"gntpc",
pERROR) <<
"Null MC record";
527 Long64_t nmax = (
gOptN<0) ?
528 er_tree->GetEntries() : TMath::Min( er_tree->GetEntries(),
gOptN );
530 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
534 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
536 TLorentzVector pdummy(0,0,0,0);
539 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
540 er_tree->GetEntry(
iev);
549 bool is_unphysical =
event.IsUnphysical();
551 LOG(
"gntpc",
pINFO) <<
"Skipping unphysical event";
602 TLorentzVector * vtx =
event.Vertex();
614 bool is_em = proc_info.
IsEM();
615 bool is_weakcc = proc_info.
IsWeakCC();
616 bool is_weaknc = proc_info.
IsWeakNC();
617 bool is_mec = proc_info.
IsMEC();
620 if (!hitnucl && neutrino) {
621 assert(is_coh || is_imd || is_imdanh || is_nuel | is_amnugamma);
625 int qrk = (is_dis) ? tgt.
HitQrkPdg() : 0;
626 bool seaq = (is_dis) ? tgt.
HitSeaQrk() :
false;
640 double weight =
event.Weight();
646 bool get_selected =
true;
647 double xs = kine.
x (get_selected);
648 double ys = kine.
y (get_selected);
649 double ts = (is_coh || is_dfr) ? kine.
t (get_selected) : -1;
650 double Q2s = kine.
Q2(get_selected);
651 double Ws = kine.
W (get_selected);
654 <<
"[Select] Q2 = " << Q2s <<
", W = " << Ws
655 <<
", x = " << xs <<
", y = " << ys <<
", t = " << ts;
661 const TLorentzVector & k1 = (neutrino) ? *(neutrino->
P4()) : pdummy;
662 const TLorentzVector & k2 = (fsl) ? *(fsl->
P4()) : pdummy;
663 const TLorentzVector &
p1 = (hitnucl) ? *(hitnucl->
P4()) : pdummy;
666 TLorentzVector q = k1-k2;
667 double Q2 = -1 * q.M2();
669 double v = (hitnucl) ? q.Energy() : -1;
673 x = (hitnucl) ? 0.5*Q2/(M*v) : -1;
674 y = (hitnucl) ? v/k1.Energy() : -1;
676 W2 = (hitnucl) ? M*M + 2*M*v - Q2 : -1;
684 W2 = M*M + 2*M*v -
Q2;
689 double t = (is_coh || is_dfr) ? kine.
t (get_selected) : -1;
692 TLorentzVector k1_rf = k1;
694 k1_rf.Boost(-1.*p1.BoostVector());
706 <<
"[Calc] Q2 = " << Q2 <<
", W = " << W
707 <<
", x = " << x <<
", y = " << y <<
", t = " <<
t;
712 bool study_hadsyst = (is_qel || is_res || is_dis || is_coh || is_dfr || is_mec || is_singlek);
715 TObjArrayIter piter(&event);
724 LOG(
"gntpc",
pDEBUG) <<
"Extracting final state hadronic system";
726 vector<int> final_had_syst;
727 while( (p = (
GHepParticle *) piter.Next()) && study_hadsyst)
732 if(event.Particle(ip)->FirstMother()==0)
continue;
741 if(event.Particle(igmom)->Pdg() !=
kPdgPi0) { final_had_syst.push_back(ip); }
744 final_had_syst.push_back(ip);
750 int fd_pdgc =
event.Particle(ifd)->Pdg();
753 final_had_syst.push_back(ip);
758 if(
count(final_had_syst.begin(), final_had_syst.end(), -1) > 0) {
771 LOG(
"gntpc",
pDEBUG) <<
"Extracting primary hadronic system";
774 TObjArrayIter piter_prim(&event);
776 vector<int> prim_had_syst;
782 for( vector<int>::const_iterator hiter = final_had_syst.begin();
783 hiter != final_had_syst.end(); ++hiter) {
785 prim_had_syst.push_back(*hiter);
811 int ist_comp = p->
Status();
818 prim_had_syst.push_back(ip);
825 int ist_comp = p->
Status();
831 prim_had_syst.push_back(ip);
838 int ist_comp = p->
Status();
845 prim_had_syst.push_back(ip);
852 int ist_comp = p->
Status();
859 prim_had_syst.push_back(ip);
877 if(
count(prim_had_syst.begin(), prim_had_syst.end(), -1) > 0) {
886 brNeutrino = (neutrino) ? neutrino->
Pdg() : 0;
887 brFSPrimLept = (fsl) ? fsl->
Pdg() : 0;
888 brTarget = target->
Pdg();
891 brHitNuc = (hitnucl) ? hitnucl->
Pdg() : 0;
901 brIsSingleK = is_singlek;
907 brIsCharmPro = charm;
908 brIsAMNuGamma= is_amnugamma;
920 brEvRF = k1_rf.Energy();
925 brEn = (hitnucl) ? p1.Energy() : 0;
926 brPxn = (hitnucl) ? p1.Px() : 0;
927 brPyn = (hitnucl) ? p1.Py() : 0;
928 brPzn = (hitnucl) ? p1.Pz() : 0;
934 brCosthl = TMath::Cos( k2.Vect().Angle(k1.Vect()) );
947 brNi = prim_had_syst.size();
948 for(
int j=0;
j<brNi;
j++) {
949 p =
event.Particle(prim_had_syst[
j]);
951 brPdgi[
j] = p->
Pdg();
970 <<
"Counting in primary hadronic system: idx = " << prim_had_syst[
j]
971 <<
" -> " << p->
Name();
976 <<
", N(n):" << brNiN
977 <<
", N(pi+):" << brNiPip
978 <<
", N(pi-):" << brNiPim
979 <<
", N(pi0):" << brNiPi0
980 <<
", N(K+,K-,K0):" << brNiKp+brNiKm+brNiK0
981 <<
", N(gamma,e-,e+):" << brNiEM
982 <<
", N(etc):" << brNiOther <<
"\n";
996 brSumKEf = (fsl) ? fsl->
KinE() : 0;
999 brNf = final_had_syst.size();
1000 for(
int j=0;
j<brNf;
j++) {
1001 p =
event.Particle(final_had_syst[
j]);
1004 int hpdg = p->
Pdg();
1006 double hKE = p->
KinE();
1007 double hpx = p->
Px();
1008 double hpy = p->
Py();
1009 double hpz = p->
Pz();
1010 double hp =
TMath::Sqrt(hpx*hpx + hpy*hpy + hpz*hpz);
1011 double hm = p->
Mass();
1012 double hcth = TMath::Cos( p->
P4()->Vect().Angle(k1.Vect()) );
1024 if ( hpdg ==
kPdgProton ) { brNfP++; brCalResp0 += hKE; }
1025 else if ( hpdg ==
kPdgAntiProton ) { brNfP++; brCalResp0 += (hE + 2*hm);}
1026 else if ( hpdg ==
kPdgNeutron ) { brNfN++; brCalResp0 += hKE; }
1027 else if ( hpdg ==
kPdgAntiNeutron ) { brNfN++; brCalResp0 += (hE + 2*hm);}
1028 else if ( hpdg ==
kPdgPiP ) { brNfPip++; brCalResp0 += hKE; }
1029 else if ( hpdg ==
kPdgPiM ) { brNfPim++; brCalResp0 += hKE; }
1030 else if ( hpdg ==
kPdgPi0 ) { brNfPi0++; brCalResp0 += (e_h *
hE); }
1031 else if ( hpdg ==
kPdgKP ) { brNfKp++; brCalResp0 += hKE; }
1032 else if ( hpdg ==
kPdgKM ) { brNfKm++; brCalResp0 += hKE; }
1033 else if ( hpdg ==
kPdgK0 ) { brNfK0++; brCalResp0 += hKE; }
1034 else if ( hpdg ==
kPdgAntiK0 ) { brNfK0++; brCalResp0 += hKE; }
1035 else if ( hpdg ==
kPdgGamma ) { brNfEM++; brCalResp0 += (e_h *
hE); }
1036 else if ( hpdg ==
kPdgElectron ) { brNfEM++; brCalResp0 += (e_h *
hE); }
1037 else if ( hpdg ==
kPdgPositron ) { brNfEM++; brCalResp0 += (e_h *
hE); }
1038 else { brNfOther++; brCalResp0 += hKE; }
1041 <<
"Counting in f/s system from hadronic vtx: idx = " << final_had_syst[
j]
1042 <<
" -> " << p->
Name();
1047 <<
", N(n):" << brNfN
1048 <<
", N(pi+):" << brNfPip
1049 <<
", N(pi-):" << brNfPim
1050 <<
", N(pi0):" << brNfPi0
1051 <<
", N(K+,K-,K0):" << brNfKp+brNfKm+brNfK0
1052 <<
", N(gamma,e-,e+):" << brNfEM
1053 <<
", N(etc):" << brNfOther <<
"\n";
1069 TFolder * genv = (TFolder*)
fin.Get(
"genv");
1070 TFolder * gconfig = (TFolder*)
fin.Get(
"gconfig");
1072 genv ->
Write(
"genv");
1073 gconfig ->
Write(
"gconfig");
1090 tree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
1093 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
1097 tree->SetBranchAddress(
"gmcrec", &mcrec);
1103 output <<
"<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>";
1105 output <<
"<!-- generated by GENIE gntpc utility -->";
1106 output << endl <<
endl;
1107 output <<
"<genie_event_list version=\"1.00\">" <<
endl;
1110 Long64_t nmax = (
gOptN<0) ?
1111 tree->GetEntries() : TMath::Min(tree->GetEntries(),
gOptN);
1113 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
1116 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
1119 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
1120 tree->GetEntry(
iev);
1131 output << endl <<
endl;
1132 output <<
" <!-- GENIE GHEP event -->" <<
endl;
1133 output <<
" <ghep np=\"" <<
event.GetEntries()
1134 <<
"\" unphysical=\"" 1135 << (
event.IsUnphysical() ?
"true" :
"false") <<
"\">" << endl;
1136 output << setiosflags(ios::scientific);
1140 output <<
" <!-- event weight -->";
1141 output <<
" <wgt> " <<
event.Weight() <<
" </wgt>";
1144 output <<
" <!-- cross sections -->";
1145 output <<
" <xsec_evnt> " <<
event.XSec() <<
" </xsec_evnt>";
1146 output <<
" <xsec_kine> " <<
event.DiffXSec() <<
" </xsec_kine>";
1149 output <<
" <!-- event vertex -->";
1150 output <<
" <vx> " <<
event.Vertex()->X() <<
" </vx>";
1151 output <<
" <vy> " <<
event.Vertex()->Y() <<
" </vy>";
1152 output <<
" <vz> " <<
event.Vertex()->Z() <<
" </vz>";
1153 output <<
" <vt> " <<
event.Vertex()->T() <<
" </vt>";
1157 output <<
" <!-- particle list -->" <<
endl;
1160 TIter event_iter(&event);
1161 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
1167 output <<
" <p idx=\"" << i <<
"\" type=\"" << type <<
"\">" <<
endl;
1169 output <<
" <pdg> " << p->
Pdg() <<
" </pdg>";
1170 output <<
" <ist> " << p->
Status() <<
" </ist>";
1173 output <<
" <mother> " 1174 <<
" <fst> " << setfill(
' ') << setw(3) << p->
FirstMother() <<
" </fst> " 1175 <<
" <lst> " << setfill(
' ') << setw(3) << p->
LastMother() <<
" </lst> " 1179 output <<
" <daughter> " 1180 <<
" <fst> " << setfill(
' ') << setw(3) << p->
FirstDaughter() <<
" </fst> " 1181 <<
" <lst> " << setfill(
' ') << setw(3) << p->
LastDaughter() <<
" </lst> " 1185 output <<
" <px> " << setfill(
' ') << setw(20) << p->
Px() <<
" </px>";
1186 output <<
" <py> " << setfill(
' ') << setw(20) << p->
Py() <<
" </py>";
1187 output <<
" <pz> " << setfill(
' ') << setw(20) << p->
Pz() <<
" </pz>";
1188 output <<
" <E> " << setfill(
' ') << setw(20) << p->
E() <<
" </E> ";
1191 output <<
" <x> " << setfill(
' ') << setw(20) << p->
Vx() <<
" </x> ";
1192 output <<
" <y> " << setfill(
' ') << setw(20) << p->
Vy() <<
" </y> ";
1193 output <<
" <z> " << setfill(
' ') << setw(20) << p->
Vz() <<
" </z> ";
1194 output <<
" <t> " << setfill(
' ') << setw(20) << p->
Vt() <<
" </t> ";
1206 output <<
" <rescatter> " << p->
RescatterCode() <<
" </rescatter>";
1210 output <<
" </p>" <<
endl;
1213 output <<
" </ghep>" <<
endl;
1219 output << endl <<
endl;
1220 output <<
"<genie_event_list version=\"1.00\">";
1225 LOG(
"gntpc",
pINFO) <<
"\nDone converting GENIE's GHEP ntuple";
1236 tree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
1239 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
1243 tree->SetBranchAddress(
"gmcrec", &mcrec);
1246 Long64_t nmax = (
gOptN<0) ?
1247 tree->GetEntries() : TMath::Min(tree->GetEntries(),
gOptN);
1249 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
1252 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
1260 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
1261 tree->GetEntry(
iev);
1271 stripped_event -> AttachSummary (nullint);
1272 stripped_event -> SetWeight (event.Weight());
1273 stripped_event -> SetVertex (*event.Vertex());
1282 p->Pdg(), ist, -1,-1,-1,-1, *p->P4(), *p->X4());
1297 LOG(
"gntpc",
pINFO) <<
"\nDone converting GENIE's GHEP ntuple";
1308 tree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
1311 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
1319 tree->SetBranchAddress(
"gmcrec", &mcrec);
1321 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 1323 tree->SetBranchAddress(
"flux", &flux_info);
1326 <<
"\n Flux drivers are not enabled." 1327 <<
"\n No flux pass-through information will be written-out in the rootracker file" 1328 <<
"\n If this isn't what you are supposed to be doing then build GENIE by adding " 1329 <<
"--with-flux-drivers in the configuration step.";
1336 Long64_t nmax = (
gOptN<0) ?
1337 tree->GetEntries() : TMath::Min(tree->GetEntries(),
gOptN);
1339 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
1342 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
1345 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
1346 tree->GetEntry(
iev);
1355 TIter event_iter(&event);
1372 LOG(
"gntpc",
pNOTICE) <<
"NEUT-like event type = " << evtype;
1379 LOG(
"gntpc",
pNOTICE) <<
"NUANCE-like event type = " << evtype;
1392 <<
event.Vertex()->X() <<
" " 1393 <<
event.Vertex()->Y() <<
" " 1394 <<
event.Vertex()->Z() <<
" " 1395 <<
event.Vertex()->T() <<
endl;
1407 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
1411 int ghep_pdgc = p->
Pdg();
1424 bool is_pi0_dec =
false;
1429 for(
int jd = ghep_fd; jd <= ghep_ld; jd++) {
1431 pi0dv.push_back(event.Particle(jd)->Pdg());
1434 sort(pi0dv.begin(), pi0dv.end());
1441 int ghep_fmpdgc = (ghep_fm==-1) ? 0 : event.Particle(ghep_fm)->Pdg();
1457 tracks.push_back(iparticle);
1464 for(
int iloop=0; iloop<=1; iloop++)
1466 for(vector<int>::const_iterator
ip = tracks.begin();
ip != tracks.end(); ++
ip)
1469 p =
event.Particle(iparticle);
1471 int ghep_pdgc = p->
Pdg();
1477 if(iloop==0 && fs)
continue;
1478 if(iloop==1 && !fs)
continue;
1493 default: ist = -999;
break;
1499 int pdgc = ghep_pdgc;
1511 double R = rnd->
RndGen().Rndm();
1517 const TLorentzVector *
p4 = p->
P4();
1524 double dcosx = (P>0) ? Px/P : -999;
1525 double dcosy = (P>0) ? Py/P : -999;
1526 double dcosz = (P>0) ? Pz/P : -999;
1542 <<
"Adding $track corrsponding to GHEP particle at position: " << iparticle
1543 <<
" (tracker status code: " << ist <<
")";
1545 output <<
"$ track " << pdgc <<
" " << E <<
" " 1546 << dcosx <<
" " << dcosy <<
" " << dcosz <<
" " 1619 <<
event.Weight() <<
" " 1620 <<
event.Probability()
1622 output <<
"$ info " <<
event.Vertex()->X() <<
" " 1623 <<
event.Vertex()->Y() <<
" " 1624 <<
event.Vertex()->Z() <<
" " 1625 <<
event.Vertex()->T()
1634 quark_id = 10 * quark_pdg + sorv;
1644 output <<
"$ info " <<
event.GetEntries() <<
endl;
1645 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) )
1653 << p->
X4()->X() <<
" " << p->
X4()->Y() <<
" " << p->
X4()->Z() <<
" " << p->
X4()->T() <<
" " 1654 << p->
P4()->Px() <<
" " << p->
P4()->Py() <<
" " << p->
P4()->Pz() <<
" " << p->
P4()->E() <<
" ";
1665 int rescat_code = -1;
1666 bool have_rescat_code =
false;
1668 if(gFileMinorVrs >= 5) {
1669 if(gFileRevisVrs >= 1) {
1670 have_rescat_code =
true;
1674 if(have_rescat_code) {
1750 LOG(
"gntpc",
pINFO) <<
"\nDone converting GENIE's GHEP ntuple";
1761 TBits* brEvtFlags = 0;
1762 TObjString* brEvtCode = 0;
1771 int brStdHepPdg [
kNPmax];
1772 int brStdHepStatus[
kNPmax];
1773 int brStdHepRescat[
kNPmax];
1774 double brStdHepX4 [
kNPmax][4];
1775 double brStdHepP4 [
kNPmax][4];
1776 double brStdHepPolz [
kNPmax][3];
1785 TObjString* brNuFileName = 0;
1790 int brNuParentDecMode;
1791 double brNuParentDecP4 [4];
1792 double brNuParentDecX4 [4];
1793 double brNuParentProP4 [4];
1794 double brNuParentProX4 [4];
1795 int brNuParentProNVtx;
1846 int brNumiFluxEvtno;
1847 double brNumiFluxNdxdz;
1848 double brNumiFluxNdydz;
1849 double brNumiFluxNpz;
1850 double brNumiFluxNenergy;
1851 double brNumiFluxNdxdznea;
1852 double brNumiFluxNdydznea;
1853 double brNumiFluxNenergyn;
1854 double brNumiFluxNwtnear;
1855 double brNumiFluxNdxdzfar;
1856 double brNumiFluxNdydzfar;
1857 double brNumiFluxNenergyf;
1858 double brNumiFluxNwtfar;
1859 int brNumiFluxNorig;
1860 int brNumiFluxNdecay;
1875 int brNumiFluxNtype;
1876 double brNumiFluxVx;
1877 double brNumiFluxVy;
1878 double brNumiFluxVz;
1879 double brNumiFluxPdpx;
1880 double brNumiFluxPdpy;
1881 double brNumiFluxPdpz;
1882 double brNumiFluxPpdxdz;
1883 double brNumiFluxPpdydz;
1884 double brNumiFluxPppz;
1885 double brNumiFluxPpenergy;
1886 int brNumiFluxPpmedium;
1887 int brNumiFluxPtype;
1888 double brNumiFluxPpvx;
1889 double brNumiFluxPpvy;
1890 double brNumiFluxPpvz;
1891 double brNumiFluxMuparpx;
1892 double brNumiFluxMuparpy;
1893 double brNumiFluxMuparpz;
1894 double brNumiFluxMupare;
1895 double brNumiFluxNecm;
1896 double brNumiFluxNimpwt;
1897 double brNumiFluxXpoint;
1898 double brNumiFluxYpoint;
1899 double brNumiFluxZpoint;
1900 double brNumiFluxTvx;
1901 double brNumiFluxTvy;
1902 double brNumiFluxTvz;
1903 double brNumiFluxTpx;
1904 double brNumiFluxTpy;
1905 double brNumiFluxTpz;
1906 double brNumiFluxTptype;
1907 double brNumiFluxTgen;
1911 double brNumiFluxTgptype;
1912 double brNumiFluxTgppx;
1914 double brNumiFluxTgppy;
1916 double brNumiFluxTgppz;
1918 double brNumiFluxTprivx;
1919 double brNumiFluxTprivy;
1920 double brNumiFluxTprivz;
1921 double brNumiFluxBeamx;
1922 double brNumiFluxBeamy;
1923 double brNumiFluxBeamz;
1924 double brNumiFluxBeampx;
1925 double brNumiFluxBeampy;
1926 double brNumiFluxBeampz;
1932 TTree * rootracker_tree =
new TTree(
"gRooTracker",
"GENIE event tree rootracker format");
1942 rootracker_tree->Branch(
"EvtFlags",
"TBits", &brEvtFlags, 32000, 1);
1943 rootracker_tree->Branch(
"EvtCode",
"TObjString", &brEvtCode, 32000, 1);
1944 rootracker_tree->Branch(
"EvtNum", &brEvtNum,
"EvtNum/I");
1945 rootracker_tree->Branch(
"EvtXSec", &brEvtXSec,
"EvtXSec/D");
1946 rootracker_tree->Branch(
"EvtDXSec", &brEvtDXSec,
"EvtDXSec/D");
1947 rootracker_tree->Branch(
"EvtWght", &brEvtWght,
"EvtWght/D");
1948 rootracker_tree->Branch(
"EvtProb", &brEvtProb,
"EvtProb/D");
1949 rootracker_tree->Branch(
"EvtVtx", brEvtVtx,
"EvtVtx[4]/D");
1950 rootracker_tree->Branch(
"StdHepN", &brStdHepN,
"StdHepN/I");
1951 rootracker_tree->Branch(
"StdHepPdg", brStdHepPdg,
"StdHepPdg[StdHepN]/I");
1952 rootracker_tree->Branch(
"StdHepStatus", brStdHepStatus,
"StdHepStatus[StdHepN]/I");
1953 rootracker_tree->Branch(
"StdHepRescat", brStdHepRescat,
"StdHepRescat[StdHepN]/I");
1954 rootracker_tree->Branch(
"StdHepX4", brStdHepX4,
"StdHepX4[StdHepN][4]/D");
1955 rootracker_tree->Branch(
"StdHepP4", brStdHepP4,
"StdHepP4[StdHepN][4]/D");
1956 rootracker_tree->Branch(
"StdHepPolz", brStdHepPolz,
"StdHepPolz[StdHepN][3]/D");
1957 rootracker_tree->Branch(
"StdHepFd", brStdHepFd,
"StdHepFd[StdHepN]/I");
1958 rootracker_tree->Branch(
"StdHepLd", brStdHepLd,
"StdHepLd[StdHepN]/I");
1959 rootracker_tree->Branch(
"StdHepFm", brStdHepFm,
"StdHepFm[StdHepN]/I");
1960 rootracker_tree->Branch(
"StdHepLm", brStdHepLm,
"StdHepLm[StdHepN]/I");
1963 rootracker_tree->Branch(
"EvtNum", &brEvtNum,
"EvtNum/I");
1964 rootracker_tree->Branch(
"EvtWght", &brEvtWght,
"EvtWght/D");
1965 rootracker_tree->Branch(
"EvtVtx", brEvtVtx,
"EvtVtx[4]/D");
1966 rootracker_tree->Branch(
"StdHepN", &brStdHepN,
"StdHepN/I");
1967 rootracker_tree->Branch(
"StdHepPdg", brStdHepPdg,
"StdHepPdg[StdHepN]/I");
1968 rootracker_tree->Branch(
"StdHepX4", brStdHepX4,
"StdHepX4[StdHepN][4]/D");
1969 rootracker_tree->Branch(
"StdHepP4", brStdHepP4,
"StdHepP4[StdHepN][4]/D");
1976 rootracker_tree->Branch(
"G2NeutEvtCode", &brNeutCode,
"G2NeutEvtCode/I");
1978 rootracker_tree->Branch(
"NuFileName",
"TObjString", &brNuFileName, 32000, 1);
1979 rootracker_tree->Branch(
"NuParentPdg", &brNuParentPdg,
"NuParentPdg/I");
1980 rootracker_tree->Branch(
"NuParentDecMode", &brNuParentDecMode,
"NuParentDecMode/I");
1981 rootracker_tree->Branch(
"NuParentDecP4", brNuParentDecP4,
"NuParentDecP4[4]/D");
1982 rootracker_tree->Branch(
"NuParentDecX4", brNuParentDecX4,
"NuParentDecX4[4]/D");
1983 rootracker_tree->Branch(
"NuParentProP4", brNuParentProP4,
"NuParentProP4[4]/D");
1984 rootracker_tree->Branch(
"NuParentProX4", brNuParentProX4,
"NuParentProX4[4]/D");
1985 rootracker_tree->Branch(
"NuParentProNVtx", &brNuParentProNVtx,
"NuParentProNVtx/I");
1987 rootracker_tree->Branch(
"NuFluxEntry", &brNuFluxEntry,
"NuFluxEntry/L");
1988 rootracker_tree->Branch(
"NuIdfd", &brNuIdfd,
"NuIdfd/I");
1989 rootracker_tree->Branch(
"NuCospibm", &brNuCospibm,
"NuCospibm/F");
1990 rootracker_tree->Branch(
"NuCospi0bm", &brNuCospi0bm,
"NuCospi0bm/F");
1991 rootracker_tree->Branch(
"NuGipart", &brNuGipart,
"NuGipart/I");
1992 rootracker_tree->Branch(
"NuGpos0", brNuGpos0,
"NuGpos0[3]/F");
1993 rootracker_tree->Branch(
"NuGvec0", brNuGvec0,
"NuGvec0[3]/F");
1994 rootracker_tree->Branch(
"NuGamom0", &brNuGamom0,
"NuGamom0/F");
1996 rootracker_tree->Branch(
"NuXnu", brNuXnu,
"NuXnu[2]/F");
1997 rootracker_tree->Branch(
"NuRnu", &brNuRnu,
"NuRnu/F");
1998 rootracker_tree->Branch(
"NuNg", &brNuNg,
"NuNg/I");
1999 rootracker_tree->Branch(
"NuGpid", brNuGpid,
"NuGpid[NuNg]/I");
2000 rootracker_tree->Branch(
"NuGmec", brNuGmec,
"NuGmec[NuNg]/I");
2001 rootracker_tree->Branch(
"NuGv", brNuGv,
"NuGv[NuNg][3]/F");
2002 rootracker_tree->Branch(
"NuGp", brNuGp,
"NuGp[NuNg][3]/F");
2003 rootracker_tree->Branch(
"NuGcosbm", brNuGcosbm,
"NuGcosbm[NuNg]/F");
2004 rootracker_tree->Branch(
"NuGmat", brNuGmat,
"NuGmat[NuNg]/I");
2005 rootracker_tree->Branch(
"NuGdistc", brNuGdistc,
"NuGdistc[NuNg]/F");
2006 rootracker_tree->Branch(
"NuGdistal", brNuGdistal,
"NuGdistal[NuNg]/F");
2007 rootracker_tree->Branch(
"NuGdistti", brNuGdistti,
"NuGdistti[NuNg]/F");
2008 rootracker_tree->Branch(
"NuGdistfe", brNuGdistfe,
"NuGdistfe[NuNg]/F");
2009 rootracker_tree->Branch(
"NuNorm", &brNuNorm,
"NuNorm/F");
2010 rootracker_tree->Branch(
"NuEnusk", &brNuEnusk,
"NuEnusk/F");
2011 rootracker_tree->Branch(
"NuNormsk", &brNuNormsk,
"NuNormsk/F");
2012 rootracker_tree->Branch(
"NuAnorm", &brNuAnorm,
"NuAnorm/F");
2013 rootracker_tree->Branch(
"NuVersion", &brNuVersion,
"NuVersion/F");
2014 rootracker_tree->Branch(
"NuNtrig", &brNuNtrig,
"NuNtrig/I");
2015 rootracker_tree->Branch(
"NuTuneid", &brNuTuneid,
"NuTuneid/I");
2016 rootracker_tree->Branch(
"NuPint", &brNuPint,
"NuPint/I");
2017 rootracker_tree->Branch(
"NuBpos", brNuBpos,
"NuBpos[2]/F");
2018 rootracker_tree->Branch(
"NuBtilt", brNuBtilt,
"NuBtilt[2]/F");
2019 rootracker_tree->Branch(
"NuBrms", brNuBrms,
"NuBrms[2]/F");
2020 rootracker_tree->Branch(
"NuEmit", brNuEmit,
"NuEmit[2]/F");
2021 rootracker_tree->Branch(
"NuAlpha", brNuAlpha,
"NuAlpha[2]/F");
2022 rootracker_tree->Branch(
"NuHcur", brNuHcur,
"NuHcur[3]/F");
2023 rootracker_tree->Branch(
"NuRand", &brNuRand,
"NuRand/I");
2031 rootracker_tree->Branch(
"NumiFluxRun", &brNumiFluxRun,
"NumiFluxRun/I");
2032 rootracker_tree->Branch(
"NumiFluxEvtno", &brNumiFluxEvtno,
"NumiFluxEvtno/I");
2033 rootracker_tree->Branch(
"NumiFluxNdxdz", &brNumiFluxNdxdz,
"NumiFluxNdxdz/D");
2034 rootracker_tree->Branch(
"NumiFluxNdydz", &brNumiFluxNdydz,
"NumiFluxNdydz/D");
2035 rootracker_tree->Branch(
"NumiFluxNpz", &brNumiFluxNpz,
"NumiFluxNpz/D");
2036 rootracker_tree->Branch(
"NumiFluxNenergy", &brNumiFluxNenergy,
"NumiFluxNenergy/D");
2037 rootracker_tree->Branch(
"NumiFluxNdxdznea", &brNumiFluxNdxdznea,
"NumiFluxNdxdznea/D");
2038 rootracker_tree->Branch(
"NumiFluxNdydznea", &brNumiFluxNdydznea,
"NumiFluxNdydznea/D");
2039 rootracker_tree->Branch(
"NumiFluxNenergyn", &brNumiFluxNenergyn,
"NumiFluxNenergyn/D");
2040 rootracker_tree->Branch(
"NumiFluxNwtnear", &brNumiFluxNwtnear,
"NumiFluxNwtnear/D");
2041 rootracker_tree->Branch(
"NumiFluxNdxdzfar", &brNumiFluxNdxdzfar,
"NumiFluxNdxdzfar/D");
2042 rootracker_tree->Branch(
"NumiFluxNdydzfar", &brNumiFluxNdydzfar,
"NumiFluxNdydzfar/D");
2043 rootracker_tree->Branch(
"NumiFluxNenergyf", &brNumiFluxNenergyf,
"NumiFluxNenergyf/D");
2044 rootracker_tree->Branch(
"NumiFluxNwtfar", &brNumiFluxNwtfar,
"NumiFluxNwtfar/D");
2045 rootracker_tree->Branch(
"NumiFluxNorig", &brNumiFluxNorig,
"NumiFluxNorig/I");
2046 rootracker_tree->Branch(
"NumiFluxNdecay", &brNumiFluxNdecay,
"NumiFluxNdecay/I");
2047 rootracker_tree->Branch(
"NumiFluxNtype", &brNumiFluxNtype,
"NumiFluxNtype/I");
2048 rootracker_tree->Branch(
"NumiFluxVx", &brNumiFluxVx,
"NumiFluxVx/D");
2049 rootracker_tree->Branch(
"NumiFluxVy", &brNumiFluxVy,
"NumiFluxVy/D");
2050 rootracker_tree->Branch(
"NumiFluxVz", &brNumiFluxVz,
"NumiFluxVz/D");
2051 rootracker_tree->Branch(
"NumiFluxPdpx", &brNumiFluxPdpx,
"NumiFluxPdpx/D");
2052 rootracker_tree->Branch(
"NumiFluxPdpy", &brNumiFluxPdpy,
"NumiFluxPdpy/D");
2053 rootracker_tree->Branch(
"NumiFluxPdpz", &brNumiFluxPdpz,
"NumiFluxPdpz/D");
2054 rootracker_tree->Branch(
"NumiFluxPpdxdz", &brNumiFluxPpdxdz,
"NumiFluxPpdxdz/D");
2055 rootracker_tree->Branch(
"NumiFluxPpdydz", &brNumiFluxPpdydz,
"NumiFluxPpdydz/D");
2056 rootracker_tree->Branch(
"NumiFluxPppz", &brNumiFluxPppz,
"NumiFluxPppz/D");
2057 rootracker_tree->Branch(
"NumiFluxPpenergy", &brNumiFluxPpenergy,
"NumiFluxPpenergy/D");
2058 rootracker_tree->Branch(
"NumiFluxPpmedium", &brNumiFluxPpmedium,
"NumiFluxPpmedium/I");
2059 rootracker_tree->Branch(
"NumiFluxPtype", &brNumiFluxPtype,
"NumiFluxPtype/I");
2060 rootracker_tree->Branch(
"NumiFluxPpvx", &brNumiFluxPpvx,
"NumiFluxPpvx/D");
2061 rootracker_tree->Branch(
"NumiFluxPpvy", &brNumiFluxPpvy,
"NumiFluxPpvy/D");
2062 rootracker_tree->Branch(
"NumiFluxPpvz", &brNumiFluxPpvz,
"NumiFluxPpvz/D");
2063 rootracker_tree->Branch(
"NumiFluxMuparpx", &brNumiFluxMuparpx,
"NumiFluxMuparpx/D");
2064 rootracker_tree->Branch(
"NumiFluxMuparpy", &brNumiFluxMuparpy,
"NumiFluxMuparpy/D");
2065 rootracker_tree->Branch(
"NumiFluxMuparpz", &brNumiFluxMuparpz,
"NumiFluxMuparpz/D");
2066 rootracker_tree->Branch(
"NumiFluxMupare", &brNumiFluxMupare,
"NumiFluxMupare/D");
2067 rootracker_tree->Branch(
"NumiFluxNecm", &brNumiFluxNecm,
"NumiFluxNecm/D");
2068 rootracker_tree->Branch(
"NumiFluxNimpwt", &brNumiFluxNimpwt,
"NumiFluxNimpwt/D");
2069 rootracker_tree->Branch(
"NumiFluxXpoint", &brNumiFluxXpoint,
"NumiFluxXpoint/D");
2070 rootracker_tree->Branch(
"NumiFluxYpoint", &brNumiFluxYpoint,
"NumiFluxYpoint/D");
2071 rootracker_tree->Branch(
"NumiFluxZpoint", &brNumiFluxZpoint,
"NumiFluxZpoint/D");
2072 rootracker_tree->Branch(
"NumiFluxTvx", &brNumiFluxTvx,
"NumiFluxTvx/D");
2073 rootracker_tree->Branch(
"NumiFluxTvy", &brNumiFluxTvy,
"NumiFluxTvy/D");
2074 rootracker_tree->Branch(
"NumiFluxTvz", &brNumiFluxTvz,
"NumiFluxTvz/D");
2075 rootracker_tree->Branch(
"NumiFluxTpx", &brNumiFluxTpx,
"NumiFluxTpx/D");
2076 rootracker_tree->Branch(
"NumiFluxTpy", &brNumiFluxTpy,
"NumiFluxTpy/D");
2077 rootracker_tree->Branch(
"NumiFluxTpz", &brNumiFluxTpz,
"NumiFluxTpz/D");
2078 rootracker_tree->Branch(
"NumiFluxTptype", &brNumiFluxTptype,
"NumiFluxTptype/I");
2079 rootracker_tree->Branch(
"NumiFluxTgen", &brNumiFluxTgen,
"NumiFluxTgen/I");
2080 rootracker_tree->Branch(
"NumiFluxTgptype", &brNumiFluxTgptype,
"NumiFluxTgptype/I");
2081 rootracker_tree->Branch(
"NumiFluxTgppx", &brNumiFluxTgppx,
"NumiFluxTgppx/D");
2082 rootracker_tree->Branch(
"NumiFluxTgppy", &brNumiFluxTgppy,
"NumiFluxTgppy/D");
2083 rootracker_tree->Branch(
"NumiFluxTgppz", &brNumiFluxTgppz,
"NumiFluxTgppz/D");
2084 rootracker_tree->Branch(
"NumiFluxTprivx", &brNumiFluxTprivx,
"NumiFluxTprivx/D");
2085 rootracker_tree->Branch(
"NumiFluxTprivy", &brNumiFluxTprivy,
"NumiFluxTprivy/D");
2086 rootracker_tree->Branch(
"NumiFluxTprivz", &brNumiFluxTprivz,
"NumiFluxTprivz/D");
2087 rootracker_tree->Branch(
"NumiFluxBeamx", &brNumiFluxBeamx,
"NumiFluxBeamx/D");
2088 rootracker_tree->Branch(
"NumiFluxBeamy", &brNumiFluxBeamy,
"NumiFluxBeamy/D");
2089 rootracker_tree->Branch(
"NumiFluxBeamz", &brNumiFluxBeamz,
"NumiFluxBeamz/D");
2090 rootracker_tree->Branch(
"NumiFluxBeampx", &brNumiFluxBeampx,
"NumiFluxBeampx/D");
2091 rootracker_tree->Branch(
"NumiFluxBeampy", &brNumiFluxBeampy,
"NumiFluxBeampy/D");
2092 rootracker_tree->Branch(
"NumiFluxBeampz", &brNumiFluxBeampz,
"NumiFluxBeampz/D");
2099 gtree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
2102 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
2106 gtree->SetBranchAddress(
"gmcrec", &mcrec);
2117 LOG(
"gntpc",
pINFO) <<
"Found T2KMetaData!";
2122 <<
"Could not find T2KMetaData attached to the event tree!";
2126 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 2129 gtree->SetBranchAddress(
"flux", &jnubeam_flux_info);
2133 gtree->SetBranchAddress(
"flux", &gnumi_flux_info);
2137 <<
"\n Flux drivers are not enabled." 2138 <<
"\n No flux pass-through information will be written-out in the rootracker file" 2139 <<
"\n If this isn't what you are supposed to be doing then build GENIE by adding " 2140 <<
"--with-flux-drivers in the configuration step.";
2144 Long64_t nmax = (
gOptN<0) ?
2145 gtree->GetEntries() : TMath::Min(gtree->GetEntries(),
gOptN);
2147 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
2150 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
2153 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
2154 gtree->GetEntry(
iev);
2162 LOG(
"gntpc",
pINFO) << *interaction;
2163 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 2165 if(jnubeam_flux_info) {
2166 LOG(
"gntpc",
pINFO) << *jnubeam_flux_info;
2168 LOG(
"gntpc",
pINFO) <<
"No JNUBEAM flux info associated with this event";
2176 if(brEvtFlags)
delete brEvtFlags;
2178 if(brEvtCode)
delete brEvtCode;
2185 for(
int k=0; k<4; k++) {
2190 brStdHepPdg [
i] = 0;
2191 brStdHepStatus[
i] = -1;
2192 brStdHepRescat[
i] = -1;
2193 for(
int k=0; k<4; k++) {
2194 brStdHepX4 [
i][k] = 0;
2195 brStdHepP4 [
i][k] = 0;
2197 for(
int k=0; k<3; k++) {
2198 brStdHepPolz [
i][k] = 0;
2206 brNuParentDecMode = 0;
2207 for(
int k=0; k<4; k++) {
2208 brNuParentDecP4 [k] = 0;
2209 brNuParentDecX4 [k] = 0;
2210 brNuParentProP4 [k] = 0;
2211 brNuParentProX4 [k] = 0;
2213 brNuParentProNVtx = 0;
2217 brNuCospibm = -999999.;
2218 brNuCospi0bm = -999999.;
2220 brNuGamom0 = -999999.;
2221 for(
int k=0; k< 3; k++){
2222 brNuGvec0[k] = -999999.;
2223 brNuGpos0[k] = -999999.;
2226 for(
int k=0; k<2; k++) {
2227 brNuXnu[k] = brNuBpos[k] = brNuBtilt[k] = brNuBrms[k] = brNuEmit[k] = brNuAlpha[k] = -999999.;
2229 for(
int k=0; k<3; k++) brNuHcur[k] = -999999.;
2231 for(
int k=0; k<3; k++){
2232 brNuGv[np][k] = -999999.;
2233 brNuGp[np][k] = -999999.;
2235 brNuGpid[np] = -999999;
2236 brNuGmec[np] = -999999;
2237 brNuGmat[np] = -999999;
2238 brNuGcosbm[np] = -999999.;
2239 brNuGdistc[np] = -999999.;
2240 brNuGdistal[np] = -999999.;
2241 brNuGdistti[np] = -999999.;
2242 brNuGdistfe[np] = -999999.;
2246 brNuNorm = -999999.;
2247 brNuEnusk = -999999.;
2248 brNuNormsk = -999999.;
2249 brNuAnorm = -999999.;
2250 brNuVersion= -999999.;
2251 brNuNtrig = -999999;
2252 brNuTuneid = -999999;
2255 if(brNuFileName)
delete brNuFileName;
2262 brEvtFlags =
new TBits(*event.EventFlags());
2263 brEvtCode =
new TObjString(event.Summary()->AsString().c_str());
2266 brEvtDXSec = (1
E+38/
units::cm2) * event.DiffXSec();
2267 brEvtWght =
event.Weight();
2268 brEvtProb =
event.Probability();
2269 brEvtVtx[0] =
event.Vertex()->X();
2270 brEvtVtx[1] =
event.Vertex()->Y();
2271 brEvtVtx[2] =
event.Vertex()->Z();
2272 brEvtVtx[3] =
event.Vertex()->T();
2276 TIter event_iter(&event);
2277 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2283 brStdHepPdg [iparticle] = p->
Pdg();
2284 brStdHepStatus[iparticle] = (
int) p->
Status();
2286 brStdHepX4 [iparticle][0] = p->
X4()->X();
2287 brStdHepX4 [iparticle][1] = p->
X4()->Y();
2288 brStdHepX4 [iparticle][2] = p->
X4()->Z();
2289 brStdHepX4 [iparticle][3] = p->
X4()->T();
2290 brStdHepP4 [iparticle][0] = p->
P4()->Px();
2291 brStdHepP4 [iparticle][1] = p->
P4()->Py();
2292 brStdHepP4 [iparticle][2] = p->
P4()->Pz();
2293 brStdHepP4 [iparticle][3] = p->
P4()->E();
2305 brStdHepN = iparticle;
2315 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 2320 if(jnubeam_flux_info) {
2322 brNuParentDecMode = jnubeam_flux_info->
mode;
2324 brNuParentDecP4 [0] = jnubeam_flux_info->
ppi * jnubeam_flux_info->
npi[0];
2325 brNuParentDecP4 [1] = jnubeam_flux_info->
ppi * jnubeam_flux_info->
npi[1];
2326 brNuParentDecP4 [2] = jnubeam_flux_info->
ppi * jnubeam_flux_info->
npi[2];
2328 TMath::Power(pdglib->
Find(brNuParentPdg)->Mass(), 2.)
2329 + TMath::Power(jnubeam_flux_info->
ppi, 2.)
2331 brNuParentDecX4 [0] = jnubeam_flux_info->
xpi[0];
2332 brNuParentDecX4 [1] = jnubeam_flux_info->
xpi[1];
2333 brNuParentDecX4 [2] = jnubeam_flux_info->
xpi[2];
2334 brNuParentDecX4 [3] = 0;
2336 brNuParentProP4 [0] = jnubeam_flux_info->
ppi0 * jnubeam_flux_info->
npi0[0];
2337 brNuParentProP4 [1] = jnubeam_flux_info->
ppi0 * jnubeam_flux_info->
npi0[1];
2338 brNuParentProP4 [2] = jnubeam_flux_info->
ppi0 * jnubeam_flux_info->
npi0[2];
2340 TMath::Power(pdglib->
Find(brNuParentPdg)->Mass(), 2.)
2341 + TMath::Power(jnubeam_flux_info->
ppi0, 2.)
2343 brNuParentProX4 [0] = jnubeam_flux_info->
xpi0[0];
2344 brNuParentProX4 [1] = jnubeam_flux_info->
xpi0[1];
2345 brNuParentProX4 [2] = jnubeam_flux_info->
xpi0[2];
2346 brNuParentProX4 [3] = 0;
2348 brNuParentProNVtx = jnubeam_flux_info->
nvtx0;
2351 brNuFluxEntry = jnubeam_flux_info->
fluxentry;
2352 brNuIdfd = jnubeam_flux_info->
idfd;
2353 brNuCospibm = jnubeam_flux_info->
cospibm;
2354 brNuCospi0bm = jnubeam_flux_info->
cospi0bm;
2355 brNuGipart = jnubeam_flux_info->
gipart;
2356 brNuGamom0 = jnubeam_flux_info->
gamom0;
2357 for(
int k=0; k<3; k++){
2358 brNuGpos0[k] = (double) jnubeam_flux_info->
gpos0[k];
2359 brNuGvec0[k] = (
double) jnubeam_flux_info->
gvec0[k];
2362 brNuXnu[0] = (double) jnubeam_flux_info->
xnu;
2363 brNuXnu[1] = (
double) jnubeam_flux_info->
ynu;
2364 brNuRnu = (double) jnubeam_flux_info->
rnu;
2365 for(
int k=0; k<2; k++){
2366 brNuBpos[k] = (double) jnubeam_flux_info->
bpos[k];
2367 brNuBtilt[k] = (
double) jnubeam_flux_info->
btilt[k];
2368 brNuBrms[k] = (double) jnubeam_flux_info->
brms[k];
2369 brNuEmit[k] = (
double) jnubeam_flux_info->
emit[k];
2370 brNuAlpha[k] = (double) jnubeam_flux_info->
alpha[k];
2372 for(
int k=0; k<3; k++) brNuHcur[k] = jnubeam_flux_info->
hcur[k];
2373 for(
int np = 0; np < flux::fNgmax; np++){
2374 brNuGv[np][0] = jnubeam_flux_info->
gvx[np];
2375 brNuGv[np][1] = jnubeam_flux_info->
gvy[np];
2376 brNuGv[np][2] = jnubeam_flux_info->
gvz[np];
2377 brNuGp[np][0] = jnubeam_flux_info->
gpx[np];
2378 brNuGp[np][1] = jnubeam_flux_info->
gpy[np];
2379 brNuGp[np][2] = jnubeam_flux_info->
gpz[np];
2380 brNuGpid[np] = jnubeam_flux_info->
gpid[np];
2381 brNuGmec[np] = jnubeam_flux_info->
gmec[np];
2382 brNuGcosbm[np] = jnubeam_flux_info->
gcosbm[np];
2383 brNuGmat[np] = jnubeam_flux_info->
gmat[np];
2384 brNuGdistc[np] = jnubeam_flux_info->
gdistc[np];
2385 brNuGdistal[np] = jnubeam_flux_info->
gdistal[np];
2386 brNuGdistti[np] = jnubeam_flux_info->
gdistti[np];
2387 brNuGdistfe[np] = jnubeam_flux_info->
gdistfe[np];
2389 brNuNg = jnubeam_flux_info->
ng;
2390 brNuNorm = jnubeam_flux_info->
norm;
2391 brNuEnusk = jnubeam_flux_info->
Enusk;
2392 brNuNormsk = jnubeam_flux_info->
normsk;
2393 brNuAnorm = jnubeam_flux_info->
anorm;
2394 brNuVersion= jnubeam_flux_info->
version;
2395 brNuNtrig = jnubeam_flux_info->
ntrig;
2396 brNuTuneid = jnubeam_flux_info->
tuneid;
2397 brNuPint = jnubeam_flux_info->
pint;
2398 brNuRand = jnubeam_flux_info->
rand;
2399 brNuFileName =
new TObjString(jnubeam_flux_info->
fluxfilename.c_str());
2408 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__ 2410 if(gnumi_flux_info) {
2411 brNumiFluxRun = gnumi_flux_info->
run;
2412 brNumiFluxEvtno = gnumi_flux_info->
evtno;
2413 brNumiFluxNdxdz = gnumi_flux_info->
ndxdz;
2414 brNumiFluxNdydz = gnumi_flux_info->
ndydz;
2415 brNumiFluxNpz = gnumi_flux_info->
npz;
2416 brNumiFluxNenergy = gnumi_flux_info->
nenergy;
2417 brNumiFluxNdxdznea = gnumi_flux_info->
ndxdznea;
2418 brNumiFluxNdydznea = gnumi_flux_info->
ndydznea;
2419 brNumiFluxNenergyn = gnumi_flux_info->
nenergyn;
2420 brNumiFluxNwtnear = gnumi_flux_info->
nwtnear;
2421 brNumiFluxNdxdzfar = gnumi_flux_info->
ndxdzfar;
2422 brNumiFluxNdydzfar = gnumi_flux_info->
ndydzfar;
2423 brNumiFluxNenergyf = gnumi_flux_info->
nenergyf;
2424 brNumiFluxNwtfar = gnumi_flux_info->
nwtfar;
2425 brNumiFluxNorig = gnumi_flux_info->
norig;
2426 brNumiFluxNdecay = gnumi_flux_info->
ndecay;
2427 brNumiFluxNtype = gnumi_flux_info->
ntype;
2428 brNumiFluxVx = gnumi_flux_info->
vx;
2429 brNumiFluxVy = gnumi_flux_info->
vy;
2430 brNumiFluxVz = gnumi_flux_info->
vz;
2431 brNumiFluxPdpx = gnumi_flux_info->
pdpx;
2432 brNumiFluxPdpy = gnumi_flux_info->
pdpy;
2433 brNumiFluxPdpz = gnumi_flux_info->
pdpz;
2434 brNumiFluxPpdxdz = gnumi_flux_info->
ppdxdz;
2435 brNumiFluxPpdydz = gnumi_flux_info->
ppdydz;
2436 brNumiFluxPppz = gnumi_flux_info->
pppz;
2437 brNumiFluxPpenergy = gnumi_flux_info->
ppenergy;
2438 brNumiFluxPpmedium = gnumi_flux_info->
ppmedium;
2439 brNumiFluxPtype = gnumi_flux_info->
ptype;
2440 brNumiFluxPpvx = gnumi_flux_info->
ppvx;
2441 brNumiFluxPpvy = gnumi_flux_info->
ppvy;
2442 brNumiFluxPpvz = gnumi_flux_info->
ppvz;
2443 brNumiFluxMuparpx = gnumi_flux_info->
muparpx;
2444 brNumiFluxMuparpy = gnumi_flux_info->
muparpy;
2445 brNumiFluxMuparpz = gnumi_flux_info->
muparpz;
2446 brNumiFluxMupare = gnumi_flux_info->
mupare;
2447 brNumiFluxNecm = gnumi_flux_info->
necm;
2448 brNumiFluxNimpwt = gnumi_flux_info->
nimpwt;
2449 brNumiFluxXpoint = gnumi_flux_info->
xpoint;
2450 brNumiFluxYpoint = gnumi_flux_info->
ypoint;
2451 brNumiFluxZpoint = gnumi_flux_info->
zpoint;
2452 brNumiFluxTvx = gnumi_flux_info->
tvx;
2453 brNumiFluxTvy = gnumi_flux_info->
tvy;
2454 brNumiFluxTvz = gnumi_flux_info->
tvz;
2455 brNumiFluxTpx = gnumi_flux_info->
tpx;
2456 brNumiFluxTpy = gnumi_flux_info->
tpy;
2457 brNumiFluxTpz = gnumi_flux_info->
tpz;
2458 brNumiFluxTptype = gnumi_flux_info->
tptype;
2459 brNumiFluxTgen = gnumi_flux_info->
tgen;
2460 brNumiFluxTgptype = gnumi_flux_info->
tgptype;
2461 brNumiFluxTgppx = gnumi_flux_info->
tgppx;
2462 brNumiFluxTgppy = gnumi_flux_info->
tgppy;
2463 brNumiFluxTgppz = gnumi_flux_info->
tgppz;
2464 brNumiFluxTprivx = gnumi_flux_info->
tprivx;
2465 brNumiFluxTprivy = gnumi_flux_info->
tprivy;
2466 brNumiFluxTprivz = gnumi_flux_info->
tprivz;
2467 brNumiFluxBeamx = gnumi_flux_info->
beamx;
2468 brNumiFluxBeamy = gnumi_flux_info->
beamy;
2469 brNumiFluxBeamz = gnumi_flux_info->
beamz;
2470 brNumiFluxBeampx = gnumi_flux_info->
beampx;
2471 brNumiFluxBeampy = gnumi_flux_info->
beampy;
2472 brNumiFluxBeampz = gnumi_flux_info->
beampz;
2478 rootracker_tree->Fill();
2484 double pot = gtree->GetWeight();
2485 rootracker_tree->SetWeight(pot);
2489 TFolder * genv = (TFolder*)
fin.Get(
"genv");
2490 TFolder * gconfig = (TFolder*)
fin.Get(
"gconfig");
2492 genv ->
Write(
"genv");
2493 gconfig ->
Write(
"gconfig");
2501 LOG(
"gntpc",
pINFO) <<
"\nDone converting GENIE's GHEP ntuple";
2524 tree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
2527 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
2531 tree->SetBranchAddress(
"gmcrec", &mcrec);
2538 TFile
fout(
"ghad.root",
"recreate");
2539 TTree * ghad =
new TTree(
"ghad",
"");
2540 ghad->Branch(
"i", &brIev,
"i/I" );
2541 ghad->Branch(
"W", &brW,
"W/D" );
2542 ghad->Branch(
"n", &brN,
"n/I" );
2543 ghad->Branch(
"pdg", brPdg,
"pdg[n]/I" );
2544 ghad->Branch(
"E", brE,
"E[n]/D" );
2545 ghad->Branch(
"px", brPx,
"px[n]/D" );
2546 ghad->Branch(
"py", brPy,
"py[n]/D" );
2547 ghad->Branch(
"pz", brPz,
"pz[n]/D" );
2551 Long64_t nmax = (
gOptN<0) ?
2552 tree->GetEntries() : TMath::Min(tree->GetEntries(),
gOptN);
2554 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
2557 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
2560 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
2561 tree->GetEntry(
iev);
2570 for(
int k=0; k<
kNPmax; k++) {
2582 const Interaction * interaction =
event.Summary();
2590 bool pass = is_cc && (is_dis || is_res);
2596 int ccnc = is_cc ? 1 : 0;
2600 if (init_state.
IsNuP ()) im = 1;
2601 else if (init_state.
IsNuN ()) im = 2;
2602 else if (init_state.
IsNuBarP ()) im = 3;
2603 else if (init_state.
IsNuBarN ()) im = 4;
2615 int nupdg = neutrino->
Pdg();
2616 int fslpdg = fsl->
Pdg();
2617 int A = target->
A();
2618 int Z = target->
Z();
2620 const TLorentzVector & k1 = *(neutrino->
P4());
2621 const TLorentzVector & k2 = *(fsl->
P4());
2627 GHepParticle * hadsyst =
event.FinalStateHadronicSystem();
2629 ph = *(hadsyst->
P4());
2633 ph = *(hadres->
P4());
2637 bool get_selected =
true;
2638 double x = kine.
x (get_selected);
2639 double y = kine.
y (get_selected);
2640 double W = kine.
W (get_selected);
2644 TIter event_iter(&event);
2647 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2653 if (pdg ==
kPdgIndep ) { hadmod=13; ihadmom=
i; }
2658 << nupdg <<
"\t" << ccnc <<
"\t" << im <<
"\t" 2659 << A <<
"\t" << Z <<
endl;
2660 output << inttyp <<
"\t" << x <<
"\t" << y <<
"\t" << W <<
"\t" 2663 << k1.Px() <<
"\t" << k1.Py() <<
"\t" << k1.Pz() <<
"\t" 2664 << k1.Energy() <<
"\t" << k1.M() <<
endl;
2666 << k2.Px() <<
"\t" << k2.Py() <<
"\t" << k2.Pz() <<
"\t" 2667 << k2.Energy() <<
"\t" << k2.M() <<
endl;
2669 << ph.Px() <<
"\t" << ph.Py() <<
"\t" << ph.Pz() <<
"\t" 2670 << ph.Energy() <<
"\t" << ph.M() <<
endl;
2676 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2678 if(i<ihadmom)
continue;
2688 int mom_pdg = mom->
Pdg();
2690 if(!skip) { hadv.push_back(i); }
2705 vector<int>::const_iterator hiter = hadv.begin();
2706 for( ; hiter != hadv.end(); ++hiter) {
2709 int pdg = particle->
Pdg();
2710 double px = particle->
P4()->Px();
2711 double py = particle->
P4()->Py();
2712 double pz = particle->
P4()->Pz();
2713 double E = particle->
P4()->Energy();
2714 double m = particle->
P4()->M();
2716 << px <<
"\t" << py <<
"\t" << pz <<
"\t" 2717 << E <<
"\t" << m <<
endl;
2741 ghad->Write(
"ghad");
2746 LOG(
"gntpc",
pINFO) <<
"\nDone converting GENIE's GHEP ntuple";
2789 TTree * tEvtTree =
new TTree(
"ginuke",
"GENIE INuke Summary Tree");
2794 tEvtTree->Branch(
"iev", &brIEv,
"iev/I" );
2795 tEvtTree->Branch(
"probe", &brProbe,
"probe/I" );
2796 tEvtTree->Branch(
"tgt" , &brTarget,
"tgt/I" );
2797 tEvtTree->Branch(
"ke", &brKE,
"ke/D" );
2798 tEvtTree->Branch(
"e", &brE,
"e/D" );
2799 tEvtTree->Branch(
"p", &brP,
"p/D" );
2800 tEvtTree->Branch(
"A", &brTgtA,
"A/I" );
2801 tEvtTree->Branch(
"Z", &brTgtZ,
"Z/I" );
2802 tEvtTree->Branch(
"vtxx", &brVtxX,
"vtxx/D" );
2803 tEvtTree->Branch(
"vtxy", &brVtxY,
"vtxy/D" );
2804 tEvtTree->Branch(
"vtxz", &brVtxZ,
"vtxz/D" );
2805 tEvtTree->Branch(
"probe_fsi", &brProbeFSI,
"probe_fsi/I" );
2806 tEvtTree->Branch(
"dist", &brDist,
"dist/D" );
2807 tEvtTree->Branch(
"nh", &brNh,
"nh/I" );
2808 tEvtTree->Branch(
"pdgh", brPdgh,
"pdgh[nh]/I" );
2809 tEvtTree->Branch(
"Eh", brEh,
"Eh[nh]/D" );
2810 tEvtTree->Branch(
"ph", brPh,
"ph[nh]/D" );
2811 tEvtTree->Branch(
"pxh", brPxh,
"pxh[nh]/D" );
2812 tEvtTree->Branch(
"pyh", brPyh,
"pyh[nh]/D" );
2813 tEvtTree->Branch(
"pzh", brPzh,
"pzh[nh]/D" );
2814 tEvtTree->Branch(
"cth", brCosth,
"cth[nh]/D" );
2815 tEvtTree->Branch(
"mh", brMh,
"mh[nh]/D" );
2816 tEvtTree->Branch(
"np", &brNp,
"np/I" );
2817 tEvtTree->Branch(
"nn", &brNn,
"nn/I" );
2818 tEvtTree->Branch(
"npip", &brNpip,
"npip/I" );
2819 tEvtTree->Branch(
"npim", &brNpim,
"npim/I" );
2820 tEvtTree->Branch(
"npi0", &brNpi0,
"npi0/I" );
2824 TTree * er_tree = 0;
2826 er_tree = dynamic_cast <TTree *> (
fin.Get(
"gtree") );
2829 LOG(
"gntpc",
pERROR) <<
"Null input tree";
2832 LOG(
"gntpc",
pINFO) <<
"Input tree header: " << *thdr;
2836 er_tree->SetBranchAddress(
"gmcrec", &mcrec);
2838 LOG(
"gntpc",
pERROR) <<
"Null MC record";
2843 Long64_t nmax = (
gOptN<0) ?
2844 er_tree->GetEntries() : TMath::Min( er_tree->GetEntries(),
gOptN );
2846 LOG(
"gntpc",
pERROR) <<
"Number of events = 0";
2849 LOG(
"gntpc",
pNOTICE) <<
"*** Analyzing: " << nmax <<
" events";
2851 for(Long64_t
iev = 0;
iev < nmax;
iev++) {
2853 er_tree->GetEntry(
iev);
2881 brProbe = probe -> Pdg();
2882 brTarget = target -> Pdg();
2883 brKE = probe -> KinE();
2885 brP = probe -> P4()->Vect().Mag();
2888 brVtxX = probe -> Vx();
2889 brVtxY = probe -> Vy();
2890 brVtxZ = probe -> Vz();
2891 brProbeFSI = probe -> RescatterCode();
2893 assert(rescattered_hadron);
2898 double x = rescattered_hadron->
Vx();
2899 double y = rescattered_hadron->
Vy();
2900 double z = rescattered_hadron->
Vz();
2901 double d2 = TMath::Power(brVtxX-x,2) +
2902 TMath::Power(brVtxY-y,2) +
2903 TMath::Power(brVtxZ-z,2);
2915 TIter event_iter(&event);
2916 while ( (p = dynamic_cast<GHepParticle *>(event_iter.Next())) ) {
2920 brPdgh[
i] = p->
Pdg();
2922 brPxh [
i] = p->
Px();
2923 brPyh [
i] = p->
Py();
2924 brPzh [
i] = p->
Pz();
2927 +brPzh[i]*brPzh[i]);
2928 brCosth[
i] = brPzh[
i]/brPh[
i];
2929 brMh [
i] = p->
Mass();
2942 int tempProbeFSI = brProbeFSI;
2943 brProbeFSI =
HAProbeFSI(tempProbeFSI, brProbe, brNh, brEh, brPdgh, brNpip, brNpim, brNpi0);
2959 LOG(
"gntpc",
pINFO) <<
"\nDone converting GENIE's GHEP ntuple";
2975 LOG(
"gntpc",
pINFO) <<
"Reading input filename";
2979 <<
"Unspecified input filename - Exiting";
2989 <<
"The input ROOT file [" 2997 LOG(
"gntpc",
pINFO) <<
"Reading output file format";
3014 LOG(
"gntpc",
pFATAL) <<
"Unknown output file format (" << fmt <<
")";
3020 LOG(
"gntpc",
pFATAL) <<
"Unspecified output file format";
3027 LOG(
"gntpc",
pINFO) <<
"Reading output filename";
3031 <<
"Unspecified output filename - Using default";
3037 LOG(
"gntpc",
pINFO) <<
"Reading number of events to analyze";
3041 <<
"Unspecified number of events to analyze - Use all";
3047 LOG(
"gntpc",
pINFO) <<
"Reading format version number";
3053 <<
"Unspecified version number - Use latest";
3064 LOG(
"gntpc",
pINFO) <<
"Reading random number seed";
3067 LOG(
"gntpc",
pINFO) <<
"Unspecified random number seed - Using default";
3075 LOG(
"gntpc",
pNOTICE) <<
"Number of events to be converted = " <<
gOptN;
3099 unsigned int L = inpname.length();
3103 if(inpname.substr(L-4, L).find(
"root") != string::npos) {
3104 inpname.erase(L-4, L);
3108 size_t pos = inpname.find(
"ghep.");
3109 if(pos != string::npos) {
3110 inpname.erase(pos, pos+4);
3114 name << inpname <<
ext;
3116 return gSystem->BaseName(name.str().c_str());
3139 string thisfile = basedir +
string(
"/src/Apps/gNtpConv.cxx");
3140 string cmd =
"less " + thisfile;
3142 gSystem->Exec(cmd.c_str());
3146 int HAProbeFSI(
int probe_fsi,
int probe_pdg,
int numh,
double E_had[],
int pdg_had[],
int numpip,
int numpim,
int numpi0)
3151 for(
int i=0;
i<numh;
i++)
3152 { energy += E_had[
i]; }
3156 if (probe_fsi==3 && numh==1)
3158 else if (energy==E_had[0] && numh==1)
3160 else if (
pdg::IsPion(probe_pdg) && numpip+numpi0+numpim==0)
3162 else if ( (
pdg::IsNucleon(probe_pdg) && numpip+numpi0+numpim==0 && numh>2 )
3163 || (probe_pdg==
kPdgGamma && energy!=E_had[0] && numpip+numpi0+numpim==0))
3165 else if ( numpip+numpi0+numpim > (
pdg::IsPion(probe_pdg) ? 1 : 0) )
3169 for(
int i = 0;
i < numh;
i++)
3171 if ( (
pdg::IsPion(probe_pdg) && ( probe_pdg==pdg_had[
i] ))
3183 for (
int iter = 0; iter < numh; iter++)
3185 if (probe_pdg==211 && pdg_had[iter]==-211) { index=8; undef=
false; }
3186 else if (probe_pdg==-211 && pdg_had[iter]==211) { index=8; undef=
false; }
3189 if (undef) { index=0; }
bool IsResonant(void) const
void RandGen(long int seed)
static void SetPrintLevel(int print_level)
int main(int argc, char **argv)
int NeutReactionCode(const GHepRecord *evrec)
double W(bool selected=false) const
int RescatterCode(void) const
int GenieMajorVrsNum(string tag)
bool IsParticle(int pdgc)
not ion or pseudo-particle
bool HitSeaQrk(void) const
bool IsWeakCC(void) const
bool IsNuBarN(void) const
is anti-neutrino + neutron?
NtpMCRecHeader hdr
record header
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
double E(void) const
Get energy.
void AddDarkMatter(double mass, double med_ratio)
static const double kNucleonMass
double Q2(const Interaction *const i)
static RandomGen * Instance()
Access instance.
void CustomizeFilename(string filename)
const TLorentzVector * P4(void) const
int FirstDaughter(void) const
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.
double PolzPolarAngle(void) const
int HitQrkPdg(void) const
bool IsInverseMuDecay(void) const
void ReadFromCommandLine(int argc, char **argv)
bool IsQuasiElastic(void) const
enum genie::EGHepStatus GHepStatus_t
bool IsNuN(void) const
is neutrino + neutron?
int NuanceReactionCode(const GHepRecord *evrec)
bool gOptCopyJobMeta
copy MC job metadata (gconfig, genv TFolders)
int IonPdgCodeToA(int pdgc)
Generated/set kinematical variables for an event.
double x(bool selected=false) const
bool IsDiffractive(void) const
double Mass(void) const
Mass that corresponds to the PDG code.
bool IsIMDAnnihilation(void) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
string gOptInpFileName
input file name
double Pz(void) const
Get Pz.
string AsString(void) const
GHepStatus_t Status(void) const
double Energy(void) const
Get energy.
Contains minimal information for tagging exclusive processes.
void ConvertToGINuke(void)
double Px(void) const
Get Px.
double y(bool selected=false) const
int LastMother(void) const
bool IsCharmEvent(void) const
double Vt(void) const
Get production time.
bool CheckRootFilename(string filename)
bool IsSingleKaon(void) const
int FirstMother(void) const
enum EGNtpcFmt GNtpcFmt_t
#define P(a, b, c, d, e, x)
string Name(void) const
Name that corresponds to the PDG code.
string DefaultOutputFile(void)
int GenieMinorVrsNum(string tag)
Summary information for an interaction.
int GeantToPdg(int geant_code)
int LastDaughter(void) const
bool IsWeakNC(void) const
bool IsNuP(void) const
is neutrino + proton?
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void ConvertToGRooTracker(void)
bool IsNuElectronElastic(void) const
static constexpr double L
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
bool IsAMNuGamma(void) const
Long64_t gOptN
number of events to process
const Kinematics & Kine(void) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
void Save(void)
get the even tree
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Resonance_t Resonance(void) const
int HAProbeFSI(int, int, int, double[], int[], int, int, int)
GNtpcFmt_t gOptOutFileFormat
output file format id
bool PolzIsSet(void) const
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
bool IsDeepInelastic(void) const
void Initialize(void)
add event
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
static RunOpt * Instance(void)
int GenieRevisVrsNum(string tag)
int gOptVersion
output file format version
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Singleton class to load & serve a TDatabasePDG.
bool HitQrkIsSet(void) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
double Vz(void) const
Get production z.
const TLorentzVector * X4(void) const
const int kPdgAntiNeutron
void ConvertToGHepMock(void)
bool IsPseudoParticle(int pdgc)
const XclsTag & ExclTag(void) const
void ConvertToGTracker(void)
bool IsNuBarP(void) const
is anti-neutrino + proton?
TRandom3 & RndGen(void) const
rnd number generator for generic usage
TParticlePDG * Find(int pdgc)
assert(nhit_max >=nhit_nbins)
virtual void AddParticle(const GHepParticle &p)
const InitialState & InitState(void) const
int LatestFormatVersionNumber(void)
const ProcessInfo & ProcInfo(void) const
double t(bool selected=false) const
int IonPdgCodeToZ(int pdgc)
long int gOptRanSeed
random number seed
void MesgThresholds(string inpfile)
double Vy(void) const
Get production y.
Command line argument parser.
double Q2(bool selected=false) const
const Target & Tgt(void) const
double PolzAzimuthAngle(void) const
void Clear(Option_t *opt="")
const int kPdgHadronicSyst
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
bool OptionExists(char opt)
was option set?
void GetCommandLineArgs(int argc, char **argv)
double Vx(void) const
Get production x.
Initial State information.
bool IsCoherent(void) const
double Py(void) const
Get Py.
string gOptOutFileName
output file name