tree_maker_simple.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 
3 //void tree_maker_simple(const std::string fGENIE,const std::string fileout)
5 {
6  std::cout << "Sorry, you must run in compiled mode" << std::endl;
7 }
8 
9 #else
10 
12 #include "CAFAna/Core/Progress.h"
18 #include "CAFAna/Core/Utilities.h"
19 #include "CAFAna/Core/Binning.h"
20 #include "CAFAna/Cuts/Cuts.h"
21 #include "CAFAna/Cuts/TruthCuts.h"
22 #include "CAFAna/Core/Spectrum.h"
25 #include "CAFAna/Core/Var.h"
26 #include "CAFAna/Vars/Vars.h"
27 #include "CAFAna/Core/OscCurve.h"
29 #include "CAFAna/Core/EventList.h"
30 #include "CAFAna/Analysis/Plots.h"
31 #include "CAFAna/Cuts/Cuts.h"
32 #include "CAFAna/Cuts/SpillCuts.h"
33 
34 #include "TLorentzVector.h"
35 #include <cmath>
36 #include <iostream>
37 #include <functional>
38 #include <list>
39 #include <memory>
40 #include <set>
41 #include <string>
42 #include <vector>
43 
44 //From Root
45 #include "TStyle.h"
46 #include "TFile.h"
47 #include "TROOT.h"
48 #include "TTree.h"
49 //#include "TLeaf.h"
50 //#include "TLeafElement.h"
51 #include "TVector3.h"
52 #include "TH1.h"
53 
55 
56 using namespace ana;
57 
58 //void tree_maker_simple(const std::string fGENIE,const std::string fileout)
60 {
61 
62  const std::string fGENIE = "defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
63 
64  SpectrumLoader loader(fGENIE);
65 
66  TFile outputFile("BdtTree9VarsOptimizedVolumes.root","RECREATE");
67 
68  TTree *fPreSelectionTree = new TTree("fPreSelectionTree", "Storing some useful variables for MVA training");
69 
70  /*******************************************************************/
71  // Define variables
72  int Run;
73  int SRun;
74  int Evt;
75  int IsCC;
76  int IsNumuNC;
77  int IsNC;
78  int IsContained;
79  int IsFid;
80  int TrueFid;
81  int IsPi0all;
82  int IsPi0020;
83  int IsPi0025;
84  int IsPi0030;
85  int IsPi0035;
86  int IsPi0040;
87  int IsPi0045;
88  int IsPi0050;
89  int Is2png;
90  int ShowerNumber;
91  int ProngNumber;
92  int nbOfShw;
93  int nbOfPng;
94  int nbOfVtx;
95  int nhit;
96 
97  float remid;
98  float michel;
99  float slcCalE;
100  float cvnncid, cvnnumu, cvnnue, cvnnutau, cvnphoton ;
101  float prong1MissingPl;
102  float prong1ContPl;
103  float prong1Length;
104  float prong1Energy;
105  float prong2Energy;
106  float prong1Width;
107  float prong1VtxGev;
108  float prong1dedx;
109  float prong2dedx;
110  float prong1epLLL;
111  float prong2epLLL;
112  float prong1epLLT;
113  float prong2epLLT;
114  float prong1epiLLL;
115  float prong2epiLLL;
116  float prong1epiLLT;
117  float prong2epiLLT;
118  float prong1epi0LLL;
119  float prong2epi0LLL;
120  float prong1epi0LLT;
121  float prong2epi0LLT;
122  float shwlidpimass;
123  float recopimass;
124  float distPngStartPos;
125  float png1dirx, png2dirx;
126  float png1diry, png2diry;
127  float png1dirz, png2dirz;
128  float Png1Pur, Png1Eff, Png1Pdg, Png1MotherPdg;
129  float Png2Pur, Png2Eff, Png2Pdg, Png2MotherPdg;
130 
131  /*****************************************************************/
132  //Define the branches on the Preselection Tree
133  fPreSelectionTree->Branch("IsCC",&IsCC,"IsCC/I");
134  fPreSelectionTree->Branch("IsNumuNC",&IsNumuNC,"IsNumuNC/I");
135  fPreSelectionTree->Branch("IsNC",&IsNC,"IsNC/I");
136  fPreSelectionTree->Branch("IsPi0all",&IsPi0all,"IsPi0all/I");
137  fPreSelectionTree->Branch("IsPi0020",&IsPi0020,"IsPi0020/I");
138  fPreSelectionTree->Branch("IsPi0025",&IsPi0025,"IsPi0025/I");
139  fPreSelectionTree->Branch("IsPi0030",&IsPi0030,"IsPi0030/I");
140  fPreSelectionTree->Branch("IsPi0035",&IsPi0035,"IsPi0035/I");
141  fPreSelectionTree->Branch("IsPi0040",&IsPi0040,"IsPi0040/I");
142  fPreSelectionTree->Branch("IsPi0045",&IsPi0045,"IsPi0045/I");
143  fPreSelectionTree->Branch("IsPi0050",&IsPi0050,"IsPi0050/I");
144  fPreSelectionTree->Branch("remid",&remid,"remid/F");
145  fPreSelectionTree->Branch("michel",&michel,"michel/F");
146 
147  //Main Variables
148  fPreSelectionTree->Branch("nhit",&nhit,"nhit/I");
149  fPreSelectionTree->Branch("prong1Length",&prong1Length,"prong1Length/F");
150  fPreSelectionTree->Branch("prong2Energy",&prong2Energy,"prong2Energy/F");
151  fPreSelectionTree->Branch("png2diry",&png2diry,"png2diry/F");
152  fPreSelectionTree->Branch("png1dirx",&png1dirx,"png1dirx/F");
153  fPreSelectionTree->Branch("png1diry",&png1diry,"png1diry/F");
154  fPreSelectionTree->Branch("prong1VtxGev",&prong1VtxGev,"prong1VtxGev/F");
155  fPreSelectionTree->Branch("prong1Width",&prong1Width,"prong1Width/F");
156  fPreSelectionTree->Branch("recopimass",&recopimass,"recopimass/F");
157  fPreSelectionTree->Branch("shwlidpimass",&shwlidpimass,"shwlidpimass/F");
158  fPreSelectionTree->Branch("prong2dedx",&prong2dedx,"prong2dedx/F");
159  fPreSelectionTree->Branch("png2dirx",&png2dirx,"png2dirx/F");
160  fPreSelectionTree->Branch("prong1Energy",&prong1Energy,"prong1Energy/F");
161  fPreSelectionTree->Branch("png2dirz",&png2dirz,"png2dirz/F");
162  fPreSelectionTree->Branch("png1dirz",&png1dirz,"png1dirz/F");
163  fPreSelectionTree->Branch("remid",&remid,"remid/F");
164  fPreSelectionTree->Branch("prong1MissingPl",&prong1MissingPl,"prong1MissingPl/F");
165  fPreSelectionTree->Branch("prong1ContPl",&prong1ContPl,"prong1ContPl/F");
166  fPreSelectionTree->Branch("prong1dedx",&prong1dedx,"prong1dedx/F");
167  fPreSelectionTree->Branch("prong1epLLL",&prong1epLLL,"prong1epLLL/F");
168  fPreSelectionTree->Branch("prong1epLLT",&prong1epLLT,"prong1epLLT/F");
169  fPreSelectionTree->Branch("prong1epiLLL",&prong1epiLLL,"prong1epiLLL/F");
170  fPreSelectionTree->Branch("prong1epiLLT",&prong1epiLLT,"prong1epiLLT/F");
171  fPreSelectionTree->Branch("prong1epi0LLL",&prong1epi0LLL,"prong1epi0LLL/F");
172  fPreSelectionTree->Branch("prong1epi0LLT",&prong1epi0LLT,"prong1epi0LLT/F");
173  fPreSelectionTree->Branch("cvnncid",&cvnncid,"cvnncid/F");
174  fPreSelectionTree->Branch("cvnnumu",&cvnnumu,"cvnnumu/F");
175  fPreSelectionTree->Branch("cvnnue",&cvnnue,"cvnnue/F");
176  fPreSelectionTree->Branch("cvnnutau",&cvnnutau,"cvnnutau/F");
177  fPreSelectionTree->Branch("cvnphoton",&cvnphoton,"cvnphoton/F");
178  fPreSelectionTree->Branch("slcCalE",&slcCalE,"slcCalE/F");
179  fPreSelectionTree->Branch("distPngStartPos",&distPngStartPos,"distPngStartPos/F");
180 
181  //Truth Matching
182  fPreSelectionTree->Branch("Png1Pdg",&Png1Pdg,"Png1Pdg/F");
183  fPreSelectionTree->Branch("Png1MotherPdg",&Png1MotherPdg,"Png1MotherPdg/F");
184  fPreSelectionTree->Branch("Png1Eff",&Png1Eff,"Png1Eff/F");
185  fPreSelectionTree->Branch("Png1Pur",&Png1Pur,"Png1Pur/F");
186  fPreSelectionTree->Branch("Png2Pdg",&Png2Pdg,"Png2Pdg/F");
187  fPreSelectionTree->Branch("Png2MotherPdg",&Png2MotherPdg,"Png2MotherPdg/F");
188  fPreSelectionTree->Branch("Png2Eff",&Png2Eff,"Png2Eff/F");
189  fPreSelectionTree->Branch("Png2Pur",&Png2Pur,"Png2Pur/F");
190 
191  /********************************************************************/
192  //Begin looping over the selected files from the SAM Database.
193  SAMQueryStatus *filesrc=new SAMQueryStatus(fGENIE);
194 
195  int Nfiles = filesrc->NFiles();
196  Progress* prog = 0;
197  TFile* file1;
198 
199  int fileIdx = -1;
200 
201  while(TFile* file1 =(TFile*)filesrc->GetNextFile())
202  {
203  //Start with some file management.
204  ++fileIdx;
205 
206  if(Nfiles >= 0 && !prog) prog = new Progress(TString::Format("Filling tree from %d files", Nfiles).Data());
207 
208  loader.HandleFile(file1, Nfiles == 1 ? prog : 0);
209 
210  if(Nfiles > 1 && prog) prog->SetProgress((fileIdx+1.)/Nfiles);
211 
212  /********************************************************************/
213  //Define the recTree object
214  TTree *recTree = (TTree*)file1->Get("recTree");
215 
216  // Every time a GetEntry(i) is done, the data in the tree is copied to this object
217  caf::StandardRecord* recTreeObject = 0;
218  recTree->SetBranchAddress("rec", &recTreeObject);
219 
220  // Do not load any branches into the memory yet
221  recTree->SetBranchStatus("*",0);
222 
223  /********************************************************************/
224  //Start loading branches into the memory.
225 
226  // Header chain
227  recTree->SetBranchStatus("hdr.run",1);
228  recTree->SetBranchStatus("hdr.subrun",1);
229  recTree->SetBranchStatus("hdr.evt",1);
230  recTree->SetBranchStatus("hdr.subevt",1);
231 
232  // MC Chain
233  recTree->SetBranchStatus("mc.nnu",1);
234  recTree->SetBranchStatus("mc.nu",1);
235  recTree->SetBranchStatus("mc.nu.iscc",1);
236  recTree->SetBranchStatus("mc.nu.mode",1);
237  recTree->SetBranchStatus("mc.nu.pdg",1);
238  recTree->SetBranchStatus("mc.nu.prim",1);
239  recTree->SetBranchStatus("mc.nu.prim.pdg",1);
240  recTree->SetBranchStatus("mc.nu.prim.p",1);
241  recTree->SetBranchStatus("mc.nu.prim.p.E",1);
242  recTree->SetBranchStatus("mc.nu.vtx.x",1);
243  recTree->SetBranchStatus("mc.nu.vtx.y",1);
244  recTree->SetBranchStatus("mc.nu.vtx.z",1);
245 
246  // Vtx Chain
247  recTree->SetBranchStatus("vtx.elastic",1);
248  recTree->SetBranchStatus("vtx.nelastic",1);
249  recTree->SetBranchStatus("vtx.elastic.vtx",1);
250  recTree->SetBranchStatus("vtx.elastic.vtx.x",1);
251  recTree->SetBranchStatus("vtx.elastic.vtx.y",1);
252  recTree->SetBranchStatus("vtx.elastic.vtx.z",1);
253 
254  // ShwLid Chain
255  recTree->SetBranchStatus("vtx.elastic.fuzzyk.nshwlid",1);
256  recTree->SetBranchStatus("vtx.elastic.fuzzyk.npng",1);
257  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png",1);
258  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid",1);
259  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop",1);
260  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop.x",1);
261  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop.y",1);
262  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop.z",1);
263  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start",1);
264  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start.x",1);
265  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start.y",1);
266  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start.z",1);
267  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.calE",1);
268  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.width",1);
269  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.len",1);
270  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.maxplanecont",1);
271  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.maxplanegap",1);
272  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir",1);
273  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir.x",1);
274  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir.y",1);
275  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir.z",1);
276 
277  //Longitudinal and Transverse Likelihoods
278  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid",1);
279  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.eplll",1);
280  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epllt",1);
281  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epilll",1);
282  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epillt",1);
283  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epi0lll",1);
284  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epi0llt",1);
285  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.vtxgev",1);
286 
287  recTree->SetBranchStatus("sand.nue.dedxpng1",1);
288  recTree->SetBranchStatus("sand.nue.dedxpng2",1);
289 
290  // Slice Chain
291  recTree->SetBranchStatus("slc.calE",1);
292  recTree->SetBranchStatus("slc.nhit",1);
293 
294  // Selection Chain
295  recTree->SetBranchStatus("sel.remid.pid",1);
296  recTree->SetBranchStatus("sel.cvn.ncid",1);
297  recTree->SetBranchStatus("sel.cvn.numuid",1);
298  recTree->SetBranchStatus("sel.cvn.nueid",1);
299  recTree->SetBranchStatus("sel.cvn.nutauid",1);
300  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.cvnpart.photonid",1);
301 
302  //Truth Matching
303  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth",1);
304  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.pdg",1);
305  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.motherpdg",1);
306  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.eff",1);
307  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.pur",1);
308 
309  //Tagging ME
310  recTree->SetBranchStatus("me.nslc",1);
311  recTree->SetBranchStatus("me.slc.mid",1);
312  recTree->SetBranchStatus("me.slc.deltat",1);
313  recTree->SetBranchStatus("me.nkalman",1);
314  recTree->SetBranchStatus("me.trkkalman.mid",1);
315  recTree->SetBranchStatus("me.trkkalman.deltat",1);
316 
317  /****************************************************************/
318  //Looping over entries
319  Int_t numberOfEntries = recTree->GetEntriesFast();
320 
321  //Begin filling the branches with variables and cuts applied.
322  for (Int_t event = 0; event < numberOfEntries; ++event)
323  {
324  recTree->GetEntry(event);
325 
326  if(recTreeObject->mc.nnu == 0) continue;
327  if(recTreeObject->vtx.nelastic == 0) continue;
328 
329  //Ensure at least one prong and shower.
330  nbOfShw = recTreeObject->vtx.elastic[0].fuzzyk.nshwlid;
331  nbOfPng = recTreeObject->vtx.elastic[0].fuzzyk.npng;
332  nbOfVtx = recTreeObject->vtx.nelastic;
333  if(nbOfPng==0) continue;
334  if(nbOfShw<=0) continue;
335  if(nbOfVtx<=0) continue;
336 
337  //Predefine containment, fiducial, and true fiducial cuts.
338  Bool_t contained = true;
339  Bool_t fiducial = false;
340  Bool_t truefid = false;
341 
342  //Check that the true vertex is inside the fiducial volume.
343  if((recTreeObject->mc.nu[0].vtx.x) < 170.0
344  && (recTreeObject->mc.nu[0].vtx.x) > -170.0
345  && (recTreeObject->mc.nu[0].vtx.y) < 170.0
346  && (recTreeObject->mc.nu[0].vtx.y) > -160.0
347  && (recTreeObject->mc.nu[0].vtx.z) < 1100.0
348  && (recTreeObject->mc.nu[0].vtx.z) > 50.0) truefid=true;
349 
350  //Check that the reconstructed vertex is within the fiducial volume.
351  if((recTreeObject->vtx.elastic[0].vtx.x) < 170.0
352  && (recTreeObject->mc.nu[0].vtx.x) > -170.0
353  && (recTreeObject->vtx.elastic[0].vtx.y) < 170.0
354  && (recTreeObject->mc.nu[0].vtx.y) > -160.0
355  && recTreeObject->vtx.elastic[0].vtx.z < 1100.0
356  && recTreeObject->vtx.elastic[0].vtx.z > 50.0) fiducial=true;
357 
358  //Check that the reconstructed shower is entirely within the containment volume.
359  for(int i=0;i<nbOfShw;i++)
360  {
361  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) > 125) contained=false;
362  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) < -100) contained=false;
363  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) > 100) contained=false;
364  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) < -120) contained=false;
365  if(recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z > 1200) contained=false;
366  if(recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z < 300) contained=false;
367  }
368 
369  IsContained = contained;
370  IsFid = fiducial;
371  TrueFid = truefid;
372 
373  if(!TrueFid) continue;
374  if(!IsContained) continue;
375  if(!IsFid) continue;
376 
377  //if(nbOfShw !=2) continue;
378  //if(nbOfShw != 1) continue;
379  if(nbOfPng !=1) continue;
380  if(nbOfVtx !=1) continue;
381 
382  /*************************************************************/
383  //Define selection cuts
384 
385  /*
386  //Slice Energy Cut
387  slcCalE = recTreeObject->slc.calE;
388  if(slcCalE<.2) continue;
389 
390  // Numu NC cut
391  Bool_t numunc = false;
392  if(recTreeObject->mc.nu[0].pdg==14 && !recTreeObject->mc.nu[0].iscc)
393  {
394  numunc = true;
395  }
396  IsNumuNC = numunc;
397  */
398 
399  //NC definition, cut is made in other script.
400  Bool_t isnc = false;
401  if(!recTreeObject->mc.nu[0].iscc)
402  {
403  isnc = true;
404  }
405  IsNC = isnc;
406 
407  //Primay Pi0 Cut
408  Bool_t ispi0all = false;
409  Bool_t ispi0020 = false;
410  Bool_t ispi0025 = false;
411  Bool_t ispi0030 = false;
412  Bool_t ispi0035 = false;
413  Bool_t ispi0040 = false;
414  Bool_t ispi0045 = false;
415  Bool_t ispi0050 = false;
416  float t=-5.0;
417  float e=-10.0;
418  float pi0mass=0.135;
419  int nbOfPiall=0;
420  int nbOfPi020=0;
421  int nbOfPi025=0;
422  int nbOfPi030=0;
423  int nbOfPi035=0;
424  int nbOfPi040=0;
425  int nbOfPi045=0;
426  int nbOfPi050=0;
427 
428  int nprim=recTreeObject->mc.nu[0].prim.size();
429  // std::cout<<"nb of Primaries: "<<nprim<<" E set to "<<e<<std::endl;
430  for(int i=0;i<nprim;i++)
431  {
432  //std::cout<<"Primaries: "<<recTreeObject->mc.nu[0].prim[i].pdg<<std::endl;
433  // std::cout<<"Primaries: "<<recTreeObject->mc.nu[0].prim[i].p.E<<std::endl;
434  if(recTreeObject->mc.nu[0].prim[i].pdg==111)
435  {
436  e=recTreeObject->mc.nu[0].prim[i].p.E;
437  // std::cout<<"Pi0 energy: "<<e<<std::endl;
438  t=e - pi0mass;
439  //std::cout<<"Primary pi0 kin en: "<<t<<std::endl;
440  if(t>0.0) nbOfPiall++;
441  if(t>0.20) nbOfPi020++;
442  if(t>0.25) nbOfPi025++;
443  if(t>0.30) nbOfPi030++;
444  if(t>0.35) nbOfPi035++;
445  if(t>0.40) nbOfPi040++;
446  if(t>0.45) nbOfPi045++;
447  if(t>0.50) nbOfPi050++;
448 
449  }
450  }
451  if(nbOfPiall>0) ispi0all=true;
452  if(nbOfPi020>0) ispi0020=true;
453  if(nbOfPi025>0) ispi0025=true;
454  if(nbOfPi030>0) ispi0030=true;
455  if(nbOfPi035>0) ispi0035=true;
456  if(nbOfPi040>0) ispi0040=true;
457  if(nbOfPi045>0) ispi0045=true;
458  if(nbOfPi050>0) ispi0050=true;
459 
460 
461  IsPi0all=ispi0all;
462  IsPi0020=ispi0020;
463  IsPi0025=ispi0025;
464  IsPi0030=ispi0030;
465  IsPi0035=ispi0035;
466  IsPi0040=ispi0040;
467  IsPi0045=ispi0045;
468  IsPi0050=ispi0050;
469 
470  /*
471  //Primary Pi0 cut
472  Bool_t ispi0 = false;
473  float t=-5.0;
474  float e=-10.0;
475  float pi0mass=0.135;
476  int nbOfPi=0;
477  int nprim=recTreeObject->mc.nu[0].prim.size();
478  for(int i=0;i<nprim;i++)
479  {
480  if(recTreeObject->mc.nu[0].prim[i].pdg==111)
481  {
482  e=recTreeObject->mc.nu[0].prim[i].p.E;
483  t=e - pi0mass;
484  if(t>0.0) nbOfPiall++;
485  if(t>0.5) nbOf
486  }
487  }
488  if(nbOfPi>0)
489  ispi0=true;
490  IsPi0=ispi0;
491  */
492 
493 
494 
495 
496  /***************************************************************/
497  //Define the main variables
498 
499  //Defining nature of samples
500  Short_t mc_nnu = recTreeObject->mc.nnu;
501  Bool_t mc_nu_iscc = recTreeObject->mc.nu[0].iscc;
502  Short_t mc_nu_pdg = recTreeObject->mc.nu[0].pdg;
503 
504  remid=recTreeObject->sel.remid.pid;
505  nhit=recTreeObject->slc.nhit;
506  prong1Length=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.len;
507  prong1Energy=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.calE;
508  prong1Width=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.width;
509  prong1VtxGev=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.vtxgev;
510  png1dirx=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.x;
511  png1diry=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.y;
512  png1dirz=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.z;
513  png2dirx=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.x;
514  png2diry=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.y;
515  png2dirz=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.z;
516  prong1dedx=recTreeObject->sand.nue.dedxpng1;
517  prong2dedx=recTreeObject->sand.nue.dedxpng2;
518  prong1epLLL=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.eplll;
519  prong1epLLT=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epllt;
520  prong1epiLLL=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epilll;
521  prong1epiLLT=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epillt;
522  prong1epi0LLL=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epi0lll;
523  prong1epi0LLT=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epi0llt;
524  prong1MissingPl=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.maxplanegap;
525  prong1ContPl=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.maxplanecont;
526 
527 
528  //Define the secondary variables
529  cvnncid=recTreeObject->sel.cvn.ncid;
530  cvnnumu=recTreeObject->sel.cvn.numuid;
531  cvnnue=recTreeObject->sel.cvn.nueid;
532  cvnnutau=recTreeObject->sel.cvn.nutauid;
533  cvnphoton=recTreeObject->vtx.elastic[0].fuzzyk.png[0].cvnpart.photonid;
534  slcCalE=recTreeObject->slc.calE;
535 
536 
537  int nslc = recTreeObject->me.nslc;
538  int ntrkklmn = recTreeObject->me.nkalman;
539  float n_me = 0;
540 
541  for (int i = 0; i < nslc; i++)
542  if (recTreeObject->me.slc[i].mid > 2. && recTreeObject->me.slc[i].deltat > 800.)
543  n_me++;
544  // Add MID selected TrkMEs
545 
546  for (int i = 0; i < ntrkklmn; i++)
547  if (recTreeObject->me.trkkalman[i].mid > 2. && recTreeObject->me.trkkalman[i].deltat > 800.)
548  n_me++;
549 
550  // checked all n_me also....but then as suggested by Dan, need to use > 2 candidates
551  // michel information is orthogonal to everything CVN (it can be CVN features or CVN output: need to understand this statement fully) knows but need to check if it helps in our case or not ..so using this new variable
552 
553  if (n_me > 2.0) n_me = 2.0;
554 
555  michel = n_me;
556 
557 
558 
559 
560  float thispi0mass1;
561  float mass1;
562  float dotproduct;
563 
564  dotproduct= recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.Unit().Dot(recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.Unit());
565  thispi0mass1=sqrt(recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.calE*recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.calE*(1-dotproduct));
566 
567  recopimass=thispi0mass1;
568  // ncvnout=recTreeObject->training.cvnfeatures.noutput;
569 
570  float thispi0mass[2];
571  float mass=1000;
572  float mindiff=1000;
573  float diff=1000;
574  for(int i=0;i<nbOfShw;i++)
575  {
576  thispi0mass[i]=recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.lid.pi0mass;
577  diff=fabs(thispi0mass[i]-pi0mass);
578  if(diff<mindiff)
579  {
580  mindiff=diff;
581  mass=thispi0mass[i];
582  }
583  }//end for
584 
585  shwlidpimass=mass;
586  //std::cout<<"Shwlid mass :" << shwlidpimass << std::endl;
587 
588  Png1Pdg=recTreeObject->vtx.elastic[0].fuzzyk.png[0].truth.pdg;
589  Png1MotherPdg=recTreeObject->vtx.elastic[0].fuzzyk.png[0].truth.motherpdg;
590  Png1Pur=recTreeObject->vtx.elastic[0].fuzzyk.png[0].truth.pur;
591  Png1Eff=recTreeObject->vtx.elastic[0].fuzzyk.png[0].truth.eff;
592 
593 
594  Png2Pdg=recTreeObject->vtx.elastic[0].fuzzyk.png[1].truth.pdg;
595  Png2MotherPdg=recTreeObject->vtx.elastic[0].fuzzyk.png[1].truth.motherpdg;
596  Png2Pur=recTreeObject->vtx.elastic[0].fuzzyk.png[1].truth.pur;
597  Png2Eff=recTreeObject->vtx.elastic[0].fuzzyk.png[1].truth.eff;
598 
599  float x1 = recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.start.x;
600  float x2 = recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.start.x;
601  float y1 = recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.start.y;
602  float y2 = recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.start.y;
603  float z1 = recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.start.z;
604  float z2 = recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.start.z;
605 
606  float distanceProngStart=(z2-z1)*(z2-z1)+(x2-x1)*(x2-x1) +(y2-y1)*(y2-y1);
607  distPngStartPos=sqrt(distanceProngStart);
608 
609 
610 
611  fPreSelectionTree->Fill();
612 
613  }// end loop over events
614 
615  recTree->Delete();
616 
617 
618  }// end for fileIdx
619 
620  if(prog)
621  {
622  prog->Done();
623  delete prog;
624  }
625  //loader.ReportExposures();
626 
627  outputFile.Write();
628  outputFile.Close();
629  std::cout<<"**** End **** "<<std::endl;
630  //return 0;
631 }
632 
633 #endif
float ncid
Likelihood Neutral Current.
Definition: SRCVNResult.h:23
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
SRMichelE me
Michel electron branch.
T sqrt(T number)
Definition: d0nt_math.hpp:156
size_t nkalman
Definition: SRMichelE.h:28
nhit
Definition: demo1.py:25
std::vector< SRTrkME > trkkalman
Definition: SRMichelE.h:22
float nutauid
Likelihood Charge Current NuTau.
Definition: SRCVNResult.h:22
float nueid
Likelihood Charge Current NuE.
Definition: SRCVNResult.h:21
void tree_maker_simple()
float pid
PID value output by kNN.
Definition: SRRemid.h:25
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
SRCVNResult cvn
Horrible hack to appease CAFAna.
Definition: SRIDBranch.h:53
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
std::vector< SRSlcME > slc
Definition: SRMichelE.h:26
SRVector3D vtx
Vertex position in detector coordinates. [cm].
Definition: SRElastic.h:40
loader
Definition: demo0.py:10
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
size_t npng
Definition: SRFuzzyK.h:26
unsigned int nhit
number of hits
Definition: SRSlice.h:22
OStream cout
Definition: OStream.cxx:6
The StandardRecord is the primary top-level object in the Common Analysis File trees.
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
float numuid
Likelihood Charge Current NuMu.
Definition: SRCVNResult.h:20
A simple ascii-art progress bar.
Definition: Progress.h:9
SRSlice slc
Slice branch: nhit, extents, time, etc.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
virtual void HandleFile(TFile *f, Progress *prog=0)
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRElastic.h:44
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Float_t e
Definition: plot.C:35
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
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
TCut fiducial(x1cut &&y1cut &&z1cut)
SRVertexBranch vtx
Vertex branch: location, time, etc.
size_t nslc
Definition: SRMichelE.h:32
enum BeamMode string