Functions
NuSTreeMaker.C File Reference
#include "TFile.h"
#include "TH1.h"
#include "TTree.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include "StandardRecord/StandardRecord.h"
#include "StandardRecord/SRTrueParticle.h"
#include "StandardRecord/SRVector3D.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueEnergy2020.h"
#include "NuXAna/Cuts/NusCuts18.h"
#include "NuXAna/Cuts/NusCuts20.h"
#include "NuXAna/Analysis/NusPlotLists.h"
#include "NuXAna/Vars/NusVars.h"
#include <string>
#include <cmath>
#include <iostream>

Go to the source code of this file.

Functions

void NuSTreeMaker (const std::string file="")
 

Function Documentation

void NuSTreeMaker ( const std::string  file = "")

Definition at line 27 of file NuSTreeMaker.C.

References caf::Proxy< caf::SRHeader >::cycle, caf::Proxy< caf::SRHeader >::det, caf::Proxy< caf::SRVertexBranch >::elastic, caf::Proxy< caf::SRHeader >::evt, MakeMiniprodValidationCuts::f, file, submit_syst::fout, caf::Proxy< caf::SRElastic >::fuzzyk, caf::Proxy< caf::StandardRecord >::hdr, MECModelEnuComparisons::i, inFile, calib::j, ana::kCaloE, ana::kEME_2020, ana::kHADE_2020, ana::kHitsPerPlane, ana::kNus20Energy, ana::kNus20FDCuts, ana::kNus20NDCuts, ana::kPPFXFluxCVWgt, caf::Proxy< caf::StandardRecord >::mc, caf::Proxy< caf::SRTruthBranch >::nnu, caf::Proxy< caf::SRFuzzyK >::npng, caf::Proxy< caf::SRFuzzyK >::npng2d, caf::Proxy< caf::SRTruthBranch >::nu, caf::Proxy< caf::SRFuzzyK >::orphCalE, outFile, caf::Proxy< caf::SRFuzzyK >::png, caf::Proxy< caf::SRFuzzyK >::png2d, make_template_knob_config::recTree, caf::Proxy< caf::SRHeader >::run, sr, string, caf::Proxy< caf::SRHeader >::subevt, caf::Proxy< caf::SRHeader >::subrun, art::to_string(), caf::Proxy< caf::SRElastic >::vtx, caf::Proxy< caf::StandardRecord >::vtx, caf::Proxy< caf::SRVector3D >::X(), caf::Proxy< caf::SRVector3D >::Y(), and caf::Proxy< caf::SRVector3D >::Z().

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
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
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
const Cut kNus20NDCuts
Definition: NusCuts20.h:102
const Var kHitsPerPlane
ifstream inFile
Definition: AnaPlotMaker.h:34
TFile * outFile
Definition: PlotXSec.C:135
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
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
The StandardRecord is the primary top-level object in the Common Analysis File trees.
TFile * file
Definition: cellShifts.C:17
const Cut kNus20FDCuts
Definition: NusCuts20.h:174
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
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
enum BeamMode string