4 #include "TTreeReader.h" 5 #include "TTreeReaderValue.h" 12 #include "CAFAna/Core/Var.h" 31 TFile*
inFile = TFile::Open(
file.c_str(),
"READ");
34 TTree*
recTree = (TTree*) inFile->Get(
"recTree");
36 recTree->SetBranchAddress(
"rec", &sr);
46 TFile*
outFile =
new TFile(fout.c_str(),
"RECREATE");
50 unsigned int runNum_FT;
51 unsigned int subrunNum_FT;
53 unsigned int eventNum_FT;
54 unsigned int sliceNum_FT;
62 unsigned int nPiZero_FT;
63 unsigned int nPiPlus_FT;
64 unsigned int nPiMinus_FT;
65 unsigned int nNeutron_FT;
66 unsigned int nProton_FT;
88 double hitsPerPlane_FT;
91 TTree* outTree =
new TTree(
"outTree",
"CorrETree");
92 outTree->Branch(
"runNum" , &runNum_FT);
93 outTree->Branch(
"subrunNum" , &subrunNum_FT);
94 outTree->Branch(
"cycleNum" , &cycleNum_FT);
95 outTree->Branch(
"eventNum" , &eventNum_FT);
96 outTree->Branch(
"sliceNum" , &sliceNum_FT);
97 outTree->Branch(
"det" , &det_FT);
98 outTree->Branch(
"parentPDG" , &parentPDG_FT);
99 outTree->Branch(
"intType" , &intType_FT);
100 outTree->Branch(
"mode" , &mode_FT);
101 outTree->Branch(
"nPng2D" , &nPng2D_FT);
102 outTree->Branch(
"nPng3D" , &nPng3D_FT);
103 outTree->Branch(
"n1HitPng2D" , &n1HitPng2D_FT);
104 outTree->Branch(
"nPiZero" , &nPiZero_FT);
105 outTree->Branch(
"nPiPlus" , &nPiPlus_FT);
106 outTree->Branch(
"nPiMinus" , &nPiMinus_FT);
107 outTree->Branch(
"nNeutron" , &nNeutron_FT);
108 outTree->Branch(
"nProton" , &nProton_FT);
109 outTree->Branch(
"recoVtxX" , &recoVtxX_FT);
110 outTree->Branch(
"recoVtxY" , &recoVtxY_FT);
111 outTree->Branch(
"recoVtxZ" , &recoVtxZ_FT);
112 outTree->Branch(
"trueVtxX" , &trueVtxX_FT);
113 outTree->Branch(
"trueVtxY" , &trueVtxY_FT);
114 outTree->Branch(
"trueVtxZ" , &trueVtxZ_FT);
115 outTree->Branch(
"isNC" , &isNC_FT);
116 outTree->Branch(
"isInDet" , &isInDet_FT);
117 outTree->Branch(
"passedNDCuts", &passedNDCuts_FT);
118 outTree->Branch(
"passedFDCuts", &passedFDCuts_FT);
119 outTree->Branch(
"y" , &y_FT);
120 outTree->Branch(
"eNu" , &eNu_FT);
121 outTree->Branch(
"eDep" , &eDep_FT);
122 outTree->Branch(
"eCal" , &eCal_FT);
123 outTree->Branch(
"eCor" , &eCor_FT);
124 outTree->Branch(
"eHad" , &eHad_FT);
125 outTree->Branch(
"eEM" , &eEM_FT);
126 outTree->Branch(
"eVis" , &eVis_FT);
127 outTree->Branch(
"ePng2D" , &ePng2D_FT);
128 outTree->Branch(
"ePng3D" , &ePng3D_FT);
129 outTree->Branch(
"ppfxWgt" , &ppfxWgt_FT);
130 outTree->Branch(
"hitsPerPlane", &hitsPerPlane_FT);
131 outTree->Branch(
"nonPngE" , &nonPngE_FT);
134 for(
int i = 0;
i < recTree->GetEntries();
i++)
137 recTree->GetEntry(
i);
143 bool hasNu = srProxy.
mc.
nnu != 0;
146 runNum_FT = srProxy.
hdr.
run;
149 eventNum_FT = srProxy.
hdr.
evt;
152 parentPDG_FT = (hasNu) ? srProxy.
mc.
nu[0].beam.ptype.GetValue() : 0;
153 intType_FT = (hasNu) ? srProxy.
mc.
nu[0].inttype.GetValue() : 0;
154 mode_FT = (hasNu) ? srProxy.
mc.
nu[0].mode.GetValue() : -1;
159 for (
int i = 0;
i < nPng2D_FT;
i++)
164 nPiZero_FT = (hasNu) ? srProxy.
mc.
nu[0].npizero.GetValue() : 0;
165 nPiPlus_FT = (hasNu) ? srProxy.
mc.
nu[0].npiplus.GetValue() : 0;
166 nPiMinus_FT = (hasNu) ? srProxy.
mc.
nu[0].npiminus.GetValue() : 0;
167 nNeutron_FT = (hasNu) ? srProxy.
mc.
nu[0].nneutron.GetValue() : 0;
168 nProton_FT = (hasNu) ? srProxy.
mc.
nu[0].nproton.GetValue() : 0;
172 trueVtxX_FT = (hasNu) ? srProxy.
mc.
nu[0].vtx.X() : -5555;
173 trueVtxY_FT = (hasNu) ? srProxy.
mc.
nu[0].vtx.Y() : -5555;
174 trueVtxZ_FT = (hasNu) ? srProxy.
mc.
nu[0].vtx.Z() : -5555;
175 isNC_FT = (hasNu) ? not srProxy.
mc.
nu[0].iscc :
false;
176 isInDet_FT = (hasNu) ? srProxy.
mc.
nu[0].isvtxcont.GetValue() :
false;
179 y_FT = (hasNu) ? srProxy.
mc.
nu[0].y.GetValue() : -5.f;
180 eNu_FT = (hasNu) ? srProxy.
mc.
nu[0].E.GetValue() : -5.f;
181 eDep_FT = (hasNu) ? y_FT * eNu_FT : -5.
f;
186 eVis_FT = (hasNu) ? srProxy.
mc.
nu[0].visE.GetValue() : -5.f;
188 for (
int j = 0;
j < nPng3D_FT;
j++)
caf::Proxy< size_t > npng
caf::Proxy< caf::SRFuzzyK > fuzzyk
caf::Proxy< std::vector< caf::SRProng > > png2d
caf::Proxy< caf::SRHeader > hdr
Proxy for caf::StandardRecord.
const Var kNus20Energy([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kNEARDET &&sr->hdr.det!=caf::kFARDET) return(double) sr->slc.calE;double pars[4][6]={{1.049, 0.795, 0.8409, 0.17, 0.82,-1.00},{1.025, 0.797, 0.9162, 0.53,-0.26,-1.48},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00}};int detCur=(sr->hdr.det==caf::kFARDET)+((sr->spill.isRHC)<< 1);double e_EM=ana::kEME_2020(sr);double e_Had=ana::kHADE_2020(sr);double e_Cal=sr->slc.calE;return(e_EM/pars[detCur][0]+e_Had/pars[detCur][1])/(pars[detCur][2]+pars[detCur][3]*std::pow(e_Cal+pars[detCur][4], 2)*std::exp(pars[detCur][5]*e_Cal));})
caf::Proxy< std::vector< caf::SRNeutrino > > nu
caf::Proxy< short int > nnu
caf::Proxy< caf::SRElastic > elastic
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
const Var kHADE_2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const double CVNha_CalE=sr->slc.calE;return std::max(CVNha_CalE-kEME_2020(sr), 0.);})
caf::Proxy< caf::SRTruthBranch > mc
caf::Proxy< float > orphCalE
The StandardRecord is the primary top-level object in the Common Analysis File trees.
caf::Proxy< caf::SRVector3D > vtx
void NuSTreeMaker(const std::string file="")
caf::Proxy< size_t > npng2d
std::string to_string(ModuleType mt)
caf::Proxy< caf::SRVertexBranch > vtx
const Var kEME_2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNem_CalE=0.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+ (double) png.cvnpart.pizeroid+ (double) png.cvnpart.electronid);double haPID=((double) png.cvnpart.protonid+ (double) png.cvnpart.pionid+ (double) png.cvnpart.neutronid+ (double) png.cvnpart.otherid+ (double) png.cvnpart.muonid);if(emPID< 0) continue;if(emPID >=haPID ) CVNem_CalE+=png_CalE;else continue;}if(kLongestProng(sr) >=500) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE;})