NuSTreeMaker.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TH1.h"
3 #include "TTree.h"
4 #include "TTreeReader.h"
5 #include "TTreeReaderValue.h"
6 
7 #include "StandardRecord/Proxy/SRProxy.h"
8 #include "StandardRecord/StandardRecord.h"
9 #include "StandardRecord/SRTrueParticle.h"
10 #include "StandardRecord/SRVector3D.h"
11 
12 #include "CAFAna/Core/Var.h"
13 #include "CAFAna/Vars/PPFXWeights.h"
14 #include "CAFAna/Vars/Vars.h"
15 
16 #include "3FlavorAna/Vars/NueEnergy2020.h"
17 
18 #include "NuXAna/Cuts/NusCuts18.h"
19 #include "NuXAna/Cuts/NusCuts20.h"
20 #include "NuXAna/Analysis/NusPlotLists.h"
21 #include "NuXAna/Vars/NusVars.h"
22 
23 #include <string>
24 #include <cmath>
25 #include <iostream>
26 
27 void NuSTreeMaker(const std::string file = "")
28 {
29 
30  // Load input file
31  TFile* inFile = TFile::Open(file.c_str(), "READ");
32 
33  // Create tree reader
34  TTree* recTree = (TTree*) inFile->Get("recTree");
35  caf::StandardRecord* sr = nullptr;
36  recTree->SetBranchAddress("rec", &sr);
37  recTree->GetEntry(0);
38  static caf::SRProxy srProxy(nullptr, nullptr, "", 0, 0);
39  srProxy = *sr;
40 
41  // Set output file
42  // just name it after the run, subrun, and cycle of the first event
43  std::string fout = "CorrTree_r" + std::to_string(srProxy.hdr.run) +
44  "_s" + std::to_string(srProxy.hdr.subrun) +
45  "_c" + std::to_string(srProxy.hdr.cycle) + ".root";
46  TFile* outFile = new TFile(fout.c_str(),"RECREATE");
47  outFile->cd();
48 
49  // Setup output tree
50  unsigned int runNum_FT;
51  unsigned int subrunNum_FT;
52  int cycleNum_FT;
53  unsigned int eventNum_FT;
54  unsigned int sliceNum_FT;
55  int det_FT;
56  int parentPDG_FT;
57  int intType_FT;
58  int mode_FT;
59  int nPng2D_FT;
60  int nPng3D_FT;
61  int n1HitPng2D_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;
67  double recoVtxX_FT;
68  double recoVtxY_FT;
69  double recoVtxZ_FT;
70  double trueVtxX_FT;
71  double trueVtxY_FT;
72  double trueVtxZ_FT;
73  bool isNC_FT;
74  bool isInDet_FT;
75  bool passedNDCuts_FT;
76  bool passedFDCuts_FT;
77  double y_FT;
78  double eNu_FT;
79  double eDep_FT;
80  double eCal_FT;
81  double eCor_FT;
82  double eHad_FT;
83  double eEM_FT;
84  double eVis_FT;
85  double ePng2D_FT;
86  double ePng3D_FT;
87  double ppfxWgt_FT;
88  double hitsPerPlane_FT;
89  double nonPngE_FT;
90 
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);
132 
133  // Fill tree
134  for(int i = 0; i < recTree->GetEntries(); i++)
135  {
136  // Set entry
137  recTree->GetEntry(i);
138 
139  // Set proxy
140  srProxy = *sr;
141 
142  // Check if there is a mc neutrino
143  bool hasNu = srProxy.mc.nnu != 0;
144 
145  // Set tree values
146  runNum_FT = srProxy.hdr.run;
147  subrunNum_FT = srProxy.hdr.subrun;
148  cycleNum_FT = srProxy.hdr.cycle;
149  eventNum_FT = srProxy.hdr.evt;
150  sliceNum_FT = srProxy.hdr.subevt;
151  det_FT = srProxy.hdr.det;
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;
155  nPng2D_FT = srProxy.vtx.elastic.fuzzyk.npng2d;
156  nPng3D_FT = srProxy.vtx.elastic.fuzzyk.npng;
157  n1HitPng2D_FT = 0;
158  ePng2D_FT = 0;
159  for (int i = 0; i < nPng2D_FT; i++)
160  {
161  if (srProxy.vtx.elastic.fuzzyk.png2d[i].nhit == 1) n1HitPng2D_FT += 1;
162  ePng2D_FT += srProxy.vtx.elastic.fuzzyk.png2d[i].calE;
163  }
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;
169  recoVtxX_FT = srProxy.vtx.elastic.vtx.X();
170  recoVtxY_FT = srProxy.vtx.elastic.vtx.Y();
171  recoVtxZ_FT = srProxy.vtx.elastic.vtx.Z();
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;
177  passedNDCuts_FT = ana::kNus20NDCuts(&srProxy);
178  passedFDCuts_FT = ana::kNus20FDCuts(&srProxy);
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;
182  eCal_FT = ana::kCaloE(&srProxy);
183  eCor_FT = ana::kNus20Energy(&srProxy);
184  eHad_FT = ana::kHADE_2020(&srProxy);
185  eEM_FT = ana::kEME_2020 (&srProxy);
186  eVis_FT = (hasNu) ? srProxy.mc.nu[0].visE.GetValue() : -5.f;
187  ePng3D_FT = 0;
188  for (int j = 0; j < nPng3D_FT; j++)
189  {
190  ePng3D_FT += srProxy.vtx.elastic.fuzzyk.png[j].calE;
191  }
192  ppfxWgt_FT = ana::kPPFXFluxCVWgt(&srProxy);
193  hitsPerPlane_FT = ana::kHitsPerPlane(&srProxy);
194  nonPngE_FT = srProxy.vtx.elastic.fuzzyk.orphCalE;
195 
196  // Fill tree
197  outTree->Fill();
198  }// End fill loop
199 
200  // Write and close
201  outTree->Write();
202  outFile->Close();
203 }// End makeCorrTree
caf::Proxy< caf::SRFuzzyK >::npng
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< caf::SRElastic >::fuzzyk
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
caf::Proxy< caf::SRFuzzyK >::png2d
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2044
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
caf::Proxy< caf::SRHeader >::subrun
caf::Proxy< unsigned int > subrun
Definition: SRProxy.h:254
caf::Proxy< caf::StandardRecord >::hdr
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
submit_syst.fout
fout
Definition: submit_syst.py:244
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
ana::kNus20Energy
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));})
Definition: NusVars.h:64
caf::Proxy< caf::SRTruthBranch >::nu
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
caf::Proxy< caf::SRTruthBranch >::nnu
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
caf::Proxy< caf::SRHeader >::run
caf::Proxy< unsigned int > run
Definition: SRProxy.h:248
ana::kNus20NDCuts
const Cut kNus20NDCuts
Definition: NusCuts20.h:102
caf::Proxy< caf::SRVertexBranch >::elastic
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
ana::kHitsPerPlane
const Var kHitsPerPlane
Definition: NueVarsExtra.cxx:449
caf::Proxy< caf::SRFuzzyK >::png
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
outFile
TFile * outFile
Definition: PlotXSec.C:135
caf::Proxy< caf::SRHeader >::subevt
caf::Proxy< short unsigned int > subevt
Definition: SRProxy.h:250
ana::kCaloE
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
calib::j
const double j
Definition: BetheBloch.cxx:29
caf::Proxy< caf::SRHeader >::cycle
caf::Proxy< int > cycle
Definition: SRProxy.h:230
ana::kHADE_2020
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.);})
Definition: NueEnergy2020.h:15
caf::Proxy< caf::SRHeader >::evt
caf::Proxy< unsigned int > evt
Definition: SRProxy.h:237
make_template_knob_config.recTree
recTree
Definition: make_template_knob_config.py:87
caf::Proxy< caf::StandardRecord >::mc
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< caf::SRFuzzyK >::orphCalE
caf::Proxy< float > orphCalE
Definition: SRProxy.h:2042
caf::StandardRecord
The StandardRecord is the primary top-level object in the Common Analysis File trees.
Definition: StandardRecord.h:33
caf::Proxy< caf::SRElastic >::vtx
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
file
TFile * file
Definition: cellShifts.C:17
NuSTreeMaker
void NuSTreeMaker(const std::string file="")
Definition: NuSTreeMaker.C:27
ana::kNus20FDCuts
const Cut kNus20FDCuts
Definition: NusCuts20.h:174
caf::Proxy< caf::SRFuzzyK >::npng2d
caf::Proxy< size_t > npng2d
Definition: SRProxy.h:2039
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
caf::Proxy< caf::StandardRecord >::vtx
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
ana::kEME_2020
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;})
Definition: NueEnergy2020.h:14
caf::Proxy< caf::SRHeader >::det
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
string
enum BeamMode string
Generated on Sun Mar 7 2021 14:13:24 for NOvA by   doxygen 1.8.11