make_dst_cosrejbdttrain.C
Go to the documentation of this file.
1 // Dan's P. heritage
2 
5 
7 
8 #include "TFile.h"
9 #include "TTree.h"
10 #include "TVector3.h"
11 #include "TH1.h"
12 
13 #include <vector>
14 #include <string>
15 #include <iostream>
16 
17 using namespace ana;
18 
19 // Macro for stripping only most interesting info out of a set of CAF files
20 // Intended for making a dst for training a CosRej MVA. This macro will
21 // make these files light enough to run a standard TMVA training script on
23 {
24 
25  // Needed for deteriming cosZ
26  TVector3 beamDir = NuMIBeamDirection(caf::kFARDET).Unit();
27 
28  // New test vars
29 
30  float inel, hpp, wshw, lidcos;
31  int nshw;
32  float starty, stopy, cosang, maxy, efrac, nplfr, lid;
33 
34  // Basic size variables
35  float nhit,calE;
36  // PID variables
37  float cvne, remid, cvnsse, cvnprelec, cvnprmuon, lidmulll, lidelll, lidemullt, lidemulll;
38  // Basic reco variables
39  float prongl,sparseassym,costhetaZ,ptp,pxp,pyp,vtxx,vtxy,vtxz, xyhit_assym;
40  // Save distances to all walls, and dist to any wall but the top
41  float disttop,distbottom,distfront,distback,disteast,distwest,distnottop;
42 
43  // Needed to specifiy signal / bkg and for making plots after training
44  // with osc weights so that we can reasonably gauge performance
45  float woscdumb;
46  int pdg,pdgorig;
47  bool isGENIE, isCC;
48 
49  // ROOT complains about needing the output file opened before makeing
50  // TTree
51  TFile* outFile = TFile::Open(out.c_str(), "recreate");
52 
53  // Need three things - a TTree for sig/bkg, POT, and livetime
54  TH1F *pot = new TH1F("pot","",1,0,1);
55  TH1F *livetime = new TH1F("livetime","",1,0,1);
56  TTree *sigTree = new TTree("sigTree","");
57  TTree *bkgTree = new TTree("bkgTree","");
58 
59  // Set branch address of all the variables we're saving
60  //general
61  sigTree->Branch("nhit", &nhit);
62  sigTree->Branch("hitsperpl", &hpp);
63  sigTree->Branch("calE", &calE);
64  sigTree->Branch("inelast", &inel);
65  sigTree->Branch("costheta", &lidcos);
66  sigTree->Branch("cosang", &cosang);
67  sigTree->Branch("costhetaZ", &costhetaZ);
68  sigTree->Branch("numofshow", &nshw);
69  sigTree->Branch("widthofshow", &wshw);
70  sigTree->Branch("efrac", &efrac);
71  sigTree->Branch("xyhitassym", &xyhit_assym);
72  sigTree->Branch("prongl", &prongl);
73  sigTree->Branch("costhetaZ", &costhetaZ);
74  // pid vars
75  sigTree->Branch("cvne", &cvne);
76  sigTree->Branch("remid", &remid);
77  sigTree->Branch("lid", &lid);
78  sigTree->Branch("cvnsse", &cvnsse);
79  sigTree->Branch("cvnprelec", &cvnprelec);
80  sigTree->Branch("cvnprmuon", &cvnprmuon);
81  sigTree->Branch("lidmulll", &lidmulll);
82  sigTree->Branch("lidelll", &lidelll);
83  sigTree->Branch("lidemullt", &lidemullt);
84  sigTree->Branch("lidemulll", &lidemulll);
85  //cosrej vars
86  sigTree->Branch("sparseassym", &sparseassym);
87  sigTree->Branch("ptp", &ptp);
88  sigTree->Branch("pxp", &pxp);
89  sigTree->Branch("pyp", &pyp);
90  //geometry
91  sigTree->Branch("starty", &starty);
92  sigTree->Branch("stopy", &stopy);
93  sigTree->Branch("maxy", &maxy);
94  sigTree->Branch("nplfr", &nplfr);
95  sigTree->Branch("vtxx", &vtxx);
96  sigTree->Branch("vtxy", &vtxy);
97  sigTree->Branch("vtxz", &vtxz);
98  sigTree->Branch("disttop", &disttop);
99  sigTree->Branch("distbottom", &distbottom);
100  sigTree->Branch("distfront", &distfront);
101  sigTree->Branch("distback", &distback);
102  sigTree->Branch("disteast", &disteast);
103  sigTree->Branch("distwest", &distwest);
104  sigTree->Branch("distnottop", &distnottop);
105 
106  sigTree->Branch("pdg", &pdg);
107  sigTree->Branch("pdgorig", &pdgorig);
108  sigTree->Branch("woscdumb", &woscdumb);
109  sigTree->Branch("isGENIE", &isGENIE);
110  sigTree->Branch("isCC", &isCC);
111 
112  // And for bkg
113  // general
114  bkgTree->Branch("nhit", &nhit);
115  bkgTree->Branch("hitsperpl", &hpp);
116  bkgTree->Branch("calE", &calE);
117  bkgTree->Branch("inelast", &inel);
118  bkgTree->Branch("costheta", &lidcos);
119  bkgTree->Branch("cosang", &cosang);
120  bkgTree->Branch("costhetaZ", &costhetaZ);
121  bkgTree->Branch("numofshow", &nshw);
122  bkgTree->Branch("widthofshow", &wshw);
123  bkgTree->Branch("efrac", &efrac);
124  bkgTree->Branch("xyhitassym", &xyhit_assym);
125  bkgTree->Branch("prongl", &prongl);
126  bkgTree->Branch("costhetaZ", &costhetaZ);
127  // pid vars
128  bkgTree->Branch("cvne", &cvne);
129  bkgTree->Branch("remid", &remid);
130  bkgTree->Branch("lid", &lid);
131  bkgTree->Branch("cvnsse", &cvnsse);
132  bkgTree->Branch("cvnprelec", &cvnprelec);
133  bkgTree->Branch("cvnprmuon", &cvnprmuon);
134  bkgTree->Branch("lidmulll", &lidmulll);
135  bkgTree->Branch("lidelll", &lidelll);
136  bkgTree->Branch("lidemullt", &lidemullt);
137  bkgTree->Branch("lidemulll", &lidemulll);
138  //cosrej vars
139  bkgTree->Branch("sparseassym", &sparseassym);
140  bkgTree->Branch("ptp", &ptp);
141  bkgTree->Branch("pxp", &pxp);
142  bkgTree->Branch("pyp", &pyp);
143  //geometry
144  bkgTree->Branch("starty", &starty);
145  bkgTree->Branch("stopy", &stopy);
146  bkgTree->Branch("maxy", &maxy);
147  bkgTree->Branch("nplfr", &nplfr);
148  bkgTree->Branch("vtxx", &vtxx);
149  bkgTree->Branch("vtxy", &vtxy);
150  bkgTree->Branch("vtxz", &vtxz);
151  bkgTree->Branch("disttop", &disttop);
152  bkgTree->Branch("distbottom", &distbottom);
153  bkgTree->Branch("distfront", &distfront);
154  bkgTree->Branch("distback", &distback);
155  bkgTree->Branch("disteast", &disteast);
156  bkgTree->Branch("distwest", &distwest);
157  bkgTree->Branch("distnottop", &distnottop);
158 
159  bkgTree->Branch("pdg", &pdg);
160  bkgTree->Branch("pdgorig", &pdgorig);
161  bkgTree->Branch("woscdumb", &woscdumb);
162  bkgTree->Branch("isGENIE", &isGENIE);
163  bkgTree->Branch("isCC", &isCC);
164 
165 
166  // Need a way to run over a set of files on the grid through SAM.
167  // Use CAFAna's way of dealing with this so that it is maintained
168  int stride = getenv("CAFANA_STRIDE") ? atoi(getenv("CAFANA_STRIDE")) : -1;
169  int offset = getenv("CAFANA_OFFSET") ? atoi(getenv("CAFANA_OFFSET")) : -1;
170  SAMQuerySource fSource(fname, stride, offset, -1);
171  std::cout << "Running over dataset with stride: " << stride << std::endl;
172  std::cout << "Running over dataset with offset: " << offset << std::endl;
173 
174  // Ask for the next file
175  while (TFile *inFile = fSource.GetNextFile()){
176 
177  // We have a file, let's get fill our TTree
178 
179  // Grab POT / livetime
180  TH1* curpot = (TH1*)inFile->Get("TotalPOT");
181  pot->Fill(0.5,curpot->Integral()); // Only nonzero for MC, where we don't
182  // care about spill cuts
183  isGENIE = curpot->Integral() > 0;
184  delete curpot;
185  TH1* curlivetime = (TH1*)inFile->Get("TotalLivetime");
186  livetime->Fill(0.5,curlivetime->Integral()); // Only nonzero for cosmics
187  delete curlivetime;
188 
189  // Grap recTree from current file
190  TTree *recTree = (TTree*)inFile->Get("recTree");
191  caf::StandardRecord* sr = 0;
192  recTree->SetBranchAddress("rec", &sr);
193 
194  // Run over the TTree's entries
195  for(int i = 0; i<recTree->GetEntries(); i++) {
196  // if (i > 10000) break;
197  recTree->GetEntry(i);
198 
199  // Apply some preselection - dominated by the CVN cut
200  if (sr->sel.cvnProd3Train.nueid < 0.5) // new cvn -ok
201  continue; // Quickly cut lost of stuff
202  if (!sr->sel.veto.keep)
203  continue; // Quickly cut lots of stuff
204  if ( !sr->vtx.elastic.IsValid )
205  continue; // Need Vertex for prongs
206  // if (sr->slc.nhit < 20 || sr->slc.nhit > 150)
207  // continue; // Require event reasonably sized
208 
209  // Grab the variables you want before storing them in a TTree
210 
211  inel = 1.f;
212  if (sr->vtx.elastic.IsValid ){
213  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
214  inel = (sr->slc.calE - sr->vtx.elastic.fuzzyk.png[0].shwlid.calE)/(sr->slc.calE);
215  }
216  }
217 
218  hpp = sr->sel.nuecosrej.hitsperplane;
219 
220  nshw = 0u;
221  if (sr->vtx.elastic.IsValid) { nshw = sr->vtx.elastic.fuzzyk.nshwlid; }
222 
223  wshw = -1000.;
224  if (sr->vtx.elastic.IsValid ) {
225  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
226  wshw = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
227  }
228  }
229 
230  lidcos = -1000.;
231  if (sr->vtx.elastic.IsValid) {
232  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
233  lidcos = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.costheta;
234  }
235  }
236 
237  starty = -1000.;
238  if (sr->vtx.elastic.IsValid ) {
239  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
240  starty = sr->vtx.elastic.fuzzyk.png[0].shwlid.start.y;
241  }
242  }
243 
244  stopy = -1000.;
245  if (sr->vtx.elastic.IsValid) {
246  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
247  stopy = sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.y;
248  }
249  }
250 
251  cosang = sr->sel.nuecosrej.cosdang;
252 
253  efrac = -1000.;
254  if (sr->vtx.elastic.IsValid) {
255  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
256  efrac = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.shwEFrac;
257  }
258  }
259 
260  nplfr = sr->sel.contain.nplanestofront;
261 
262  maxy = -1000.0f;
263  if (sr->vtx.elastic.IsValid) {
264  if (sr->vtx.elastic.fuzzyk.nshwlid != 0) {
265  float maxyall = -999.0;
266  for (unsigned int i=0;i<sr->vtx.elastic.fuzzyk.nshwlid;i++){
267  maxyall = std::max(std::max(sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Y(),sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Y()),maxyall);
268  maxy = maxyall;
269  }
270  }
271  }
272 
273  // Energy estimators
274  nhit = sr->slc.nhit;
275  calE = sr->slc.calE;
276 
277  // Various PID's
278  cvne = sr->sel.cvn2017.nueid;
279  cvnsse = sr->sel.cvnProd3Train.nueid;
280  remid = (sr->trk.kalman.ntracks>0) ? float(sr->sel.remid.pid) : float(-1);
281  lid = sr->sel.lid.ann;
282  cvnprelec = -1000.;
283  if (sr->vtx.elastic.IsValid ) {
284  if (sr->vtx.elastic.fuzzyk.npng >= 1) {
285  cvnprelec = sr->vtx.elastic.fuzzyk.png[0].cvnpart.electronid;
286  }
287  }
288  cvnprmuon = -1000.;
289  if (sr->vtx.elastic.IsValid ) {
290  if (sr->vtx.elastic.fuzzyk.npng >= 1) {
291  cvnprmuon = sr->vtx.elastic.fuzzyk.png[0].cvnpart.muonid;
292  }
293  }
294  lidmulll = -1000.;
295  if (sr->vtx.elastic.IsValid ) {
296  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
297  lidmulll = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.mulll;
298  }
299  }
300  lidelll = -1000.;
301  if (sr->vtx.elastic.IsValid ) {
302  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
303  lidelll = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.elll;
304  }
305  }
306  lidemullt = -1000.;
307  if (sr->vtx.elastic.IsValid ) {
308  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
309  lidemullt = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.emullt;
310  }
311  }
312  lidemulll = -1000.;
313  if (sr->vtx.elastic.IsValid ) {
314  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
315  lidemulll = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.emulll;
316  }
317  }
318  // Some reco vars that can help pid / reject back
319  prongl = -1;
320  if (sr->vtx.elastic.IsValid ){
321  if (sr->vtx.elastic.fuzzyk.npng > 0){
322  prongl = sr->vtx.elastic.fuzzyk.png[0].len;
323  }
324  }
325  sparseassym = sr->sel.nuecosrej.sparsenessasymm;
326  xyhit_assym = -1000.0;
327  if (sr->vtx.elastic.IsValid ) {
328  if (sr->vtx.elastic.fuzzyk.nshwlid >= 1) {
329  float xyhitdiff = abs(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx -
330  sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity);
331  float xyhitsum = (sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx +
332  sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity);
333  if(xyhitsum>0) xyhit_assym = xyhitdiff/xyhitsum;
334  }
335  }
336 
337  costhetaZ = -2;
338  if (sr->vtx.elastic.IsValid ){
339  if (sr->vtx.elastic.fuzzyk.npng > 0){
340  TVector3 dir(sr->vtx.elastic.fuzzyk.png[0].dir);
341  costhetaZ = dir.Dot(beamDir);
342  }
343  }
344  ptp = sr->sel.nuecosrej.photptp;
345  pxp = sr->sel.nuecosrej.photpxp;
346  pyp = sr->sel.nuecosrej.photpyp;
347 
348  vtxx = -1e3;
349  vtxy = -1e3;
350  vtxz = -1e3;
351  if (sr->vtx.elastic.IsValid ){
352  TVector3 vtx = TVector3(sr->vtx.elastic.vtx);
353  vtxx = vtx[0];
354  vtxy = vtx[1];
355  vtxz = vtx[2];
356  }
357 
358  // Containment variables
359  disttop = sr->sel.nuecosrej.distallpngtop;
360  distbottom = sr->sel.nuecosrej.distallpngbottom;
361  distwest = sr->sel.nuecosrej.distallpngwest;
362  disteast = sr->sel.nuecosrej.distallpngeast;
363  distfront = sr->sel.nuecosrej.distallpngfront;
364  distback = sr->sel.nuecosrej.distallpngback;
365  // A way to do this
366  distnottop = distbottom;
367  if (distwest < distnottop) distnottop = distwest;
368  if (disteast < distnottop) distnottop = disteast;
369  if (distfront < distnottop) distnottop = distfront;
370  if (distback < distnottop) distnottop = distback;
371 
372  // Require an MCNeutrino to fill truth
373  // Don't continue if no neutrino - cosmic slices have no neutrino
374  if (sr->mc.nnu != 0){
375  pdg = sr->mc.nu[0].pdg;
376  pdgorig = sr->mc.nu[0].pdgorig;
377  isCC = sr->mc.nu[0].iscc;
378  woscdumb = sr->mc.nu[0].woscdumb;
379  }
380 
381  if(isGENIE && isCC && std::abs(pdgorig)==14 && std::abs(pdg)==12) // it will train numu->nue CC (FHC mode) vs everything else
382  sigTree->Fill();
383  else
384  bkgTree->Fill();
385  } // run over recTree
386 
387  }
388 
389 
390  // Go back to your output file and save your TTree
391  outFile->cd();
392  sigTree->Write();
393  bkgTree->Write();
394  pot->Write();
395  livetime->Write();
396  outFile->Close();
397 }
T max(const caf::Proxy< T > &a, T b)
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
float distallpngbottom
Definition: SRNueCosRej.h:95
size_t ntracks
Definition: SRKalman.h:23
double maxy
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
nhit
Definition: demo1.py:25
void abs(TH1 *hist)
float vtxz
Definition: NusVarsTemp.cxx:32
float abs(float number)
Definition: d0nt_math.hpp:39
ifstream inFile
Definition: AnaPlotMaker.h:34
SRKalman kalman
Tracks produced by KalmanTrack.
Definition: SRTrackBranch.h:24
float pid
PID value output by kNN.
Definition: SRRemid.h:25
bool keep
Definition: SRVeto.h:34
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
TFile * outFile
Definition: PlotXSec.C:135
float distallpngfront
Definition: SRNueCosRej.h:99
std::string getenv(std::string const &name)
SRNueCosRej nuecosrej
Output from NueCosRej (Nue Cosmic Rejection)
Definition: SRIDBranch.h:48
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
#define pot
caf::StandardRecord * sr
File source based on a SAM query or dataset (definition)
SRVector3D vtx
Vertex position in detector coordinates. [cm].
Definition: SRElastic.h:40
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
size_t npng
Definition: SRFuzzyK.h:26
float ann
ann output for the slice, currently the same as most energetic shower
Definition: SRELid.h:24
unsigned int nhit
number of hits
Definition: SRSlice.h:22
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The StandardRecord is the primary top-level object in the Common Analysis File trees.
void make_dst_cosrejbdttrain(std::string fname, std::string out)
float vtxy
Definition: NusVarsTemp.cxx:31
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
SRELid lid
Output from LIDBuilder (LID) package.
Definition: SRIDBranch.h:42
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
virtual TFile * GetNextFile() override
Returns the next file in sequence, ready for reading.
int nplanestofront
number of planes between the front of the detector (configuration dependent) and hit with the smalles...
Definition: SRContain.h:23
float vtxx
Definition: NusVarsTemp.cxx:30
SRSlice slc
Slice branch: nhit, extents, time, etc.
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:333
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRElastic.h:44
SRVeto veto
Output from CosmicVeto (Preliminary preselection)
Definition: SRIDBranch.h:52
SRTrackBranch trk
Track branch: nhit, len, etc.
SRContain contain
Output from SRContain (containment related variables)
Definition: SRIDBranch.h:51
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
SRVertexBranch vtx
Vertex branch: location, time, etc.