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)
4 void tree_maker_simple()
5 {
6  std::cout << "Sorry, you must run in compiled mode" << std::endl;
7 }
8 
9 #else
10 
11 #include "CAFAna/Core/SpectrumLoaderBase.h"
12 #include "CAFAna/Core/Progress.h"
13 #include "CAFAna/Core/SAMQuerySource.h"
14 #include "CAFAna/Core/SAMProjectSource.h"
15 #include "CAFAna/Core/WildcardSource.h"
16 #include "CAFAna/Core/IFileSource.h"
17 #include "CAFAna/Core/FileListSource.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"
23 #include "CAFAna/Core/SpectrumLoader.h"
24 #include "CAFAna/Core/OscillatableSpectrum.h"
25 #include "CAFAna/Core/Var.h"
26 #include "CAFAna/Vars/Vars.h"
27 #include "CAFAna/Core/OscCurve.h"
28 #include "CAFAna/Core/ReweightableSpectrum.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 
54 #include "StandardRecord/StandardRecord.h"
55 
56 using namespace ana;
57 
58 //void tree_maker_simple(const std::string fGENIE,const std::string fileout)
59 void tree_maker_simple()
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
caf::SRCVNResult::ncid
float ncid
Likelihood Neutral Current.
Definition: SRCVNResult.h:23
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
caf::SRFuzzyK::png
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
caf::StandardRecord::me
SRMichelE me
Michel electron branch.
Definition: StandardRecord.h:45
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::SRMichelE::nkalman
size_t nkalman
Definition: SRMichelE.h:28
make_root_from_grid_output.remid
remid
Definition: make_root_from_grid_output.py:79
submit_syst.x2
x2
Definition: submit_syst.py:250
demo1.nhit
nhit
Definition: demo1.py:25
submit_syst.y2
y2
Definition: submit_syst.py:250
caf::SRMichelE::trkkalman
std::vector< SRTrkME > trkkalman
Definition: SRMichelE.h:22
caf::SRCVNResult::nutauid
float nutauid
Likelihood Charge Current NuTau.
Definition: SRCVNResult.h:22
caf::SRCVNResult::nueid
float nueid
Likelihood Charge Current NuE.
Definition: SRCVNResult.h:21
tree_maker_simple
void tree_maker_simple()
Definition: tree_maker_simple.C:59
caf::SRRemid::pid
float pid
PID value output by kNN.
Definition: SRRemid.h:25
caf::SRFuzzyK::nshwlid
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
caf::SRIDBranch::cvn
SRCVNResult cvn
Horrible hack to appease CAFAna.
Definition: SRIDBranch.h:53
caf::SRIDBranch::remid
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
caf::SRSlice::calE
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
release_diff.diff
diff
Definition: release_diff.py:80
caf::SRMichelE::slc
std::vector< SRSlcME > slc
Definition: SRMichelE.h:26
caf::SRElastic::vtx
SRVector3D vtx
Vertex position in detector coordinates. [cm].
Definition: SRElastic.h:40
demo0.loader
loader
Definition: demo0.py:10
caf::SRTruthBranch::nnu
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
caf::SRFuzzyK::npng
size_t npng
Definition: SRFuzzyK.h:26
caf::SRSlice::nhit
unsigned int nhit
number of hits
Definition: SRSlice.h:22
om::cout
OStream cout
Definition: OStream.cxx:6
make_template_knob_config.recTree
recTree
Definition: make_template_knob_config.py:87
caf::StandardRecord
The StandardRecord is the primary top-level object in the Common Analysis File trees.
Definition: StandardRecord.h:33
ana::Progress::SetProgress
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
caf::StandardRecord::sel
SRIDBranch sel
Selector (PID) branch.
Definition: StandardRecord.h:47
caf::SRVertexBranch::elastic
SRElastic elastic
Single vertex found by Elastic Arms.
Definition: SRVertexBranch.h:26
caf::SRCVNResult::numuid
float numuid
Likelihood Charge Current NuMu.
Definition: SRCVNResult.h:20
ana::Progress
A simple ascii-art progress bar.
Definition: Progress.h:9
caf::StandardRecord::slc
SRSlice slc
Slice branch: nhit, extents, time, etc.
Definition: StandardRecord.h:42
caf::StandardRecord::mc
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
Definition: StandardRecord.h:48
ana::SpectrumLoader::HandleFile
virtual void HandleFile(TFile *f, Progress *prog=0)
Definition: SpectrumLoader.cxx:128
caf::SRElastic::fuzzyk
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRElastic.h:44
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
e
Float_t e
Definition: plot.C:35
ana::Progress::Done
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
caf::SRTruthBranch::nu
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
fiducial
TCut fiducial(x1cut &&y1cut &&z1cut)
caf::StandardRecord::vtx
SRVertexBranch vtx
Vertex branch: location, time, etc.
Definition: StandardRecord.h:44
caf::SRMichelE::nslc
size_t nslc
Definition: SRMichelE.h:32
string
enum BeamMode string
Generated on Sun Feb 28 2021 14:13:14 for NOvA by   doxygen 1.8.11