useBDTG.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 
3 void useBDTG()
4 {
5  std::cout << "Sorry, you must run in compiled mode" << std::endl;
6 }
7 
8 #else
9 
10 #include "CAFAna/Core/SpectrumLoaderBase.h"
11 #include "CAFAna/Core/Progress.h"
12 #include "CAFAna/Core/SAMQuerySource.h"
13 #include "CAFAna/Core/SAMProjectSource.h"
14 #include "CAFAna/Core/WildcardSource.h"
15 #include "CAFAna/Core/IFileSource.h"
16 #include "CAFAna/Core/FileListSource.h"
17 #include "CAFAna/Core/Utilities.h"
18 #include "CAFAna/Core/Binning.h"
19 #include "CAFAna/Cuts/Cuts.h"
20 #include "CAFAna/Core/Spectrum.h"
21 #include "CAFAna/Core/SpectrumLoader.h"
22 #include "CAFAna/Core/OscillatableSpectrum.h"
23 #include "CAFAna/Vars/Vars.h"
24 #include "CAFAna/Core/OscCurve.h"
25 #include "CAFAna/Core/ReweightableSpectrum.h"
26 #include "CAFAna/Core/EventList.h"
27 #include "CAFAna/Analysis/Plots.h"
28 #include "CAFAna/Cuts/Cuts.h"
29 #include "CAFAna/Cuts/SpillCuts.h"
30 
31 
32 #include <iostream>
33 #include <fstream>
34 #include <functional>
35 #include <list>
36 #include <memory>
37 #include <set>
38 #include <string>
39 #include <vector>
40 // ROOT includes
41 #include "TStyle.h"
42 #include "TFile.h"
43 #include "TROOT.h"
44 #include "TH2F.h"
45 #include "TTree.h"
46 #include "TVector3.h"
47 
48 #include "StandardRecord/StandardRecord.h"
49 #include "TStopwatch.h"
50 
51 #include "weights/MVAnalysis_BDTG.class.C"
52 
53 #if not defined(__CINT__) || defined(__MAKECINT__)
54 #include "TMVA/Tools.h"
55 #include "TMVA/Reader.h"
56 #include "TMVA/MethodCuts.h"
57 #endif
58 
59 using namespace ana;
60 //using namespace TMVA;
61 
62 //void useBDTG(const std::string fin,const std::string fout){
63 void useBDTG(){
64 
65  const std::string fin = "defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
66 
67  const std::string fout = "BDTGTest.root";
68 
69  /* ofstream filea;
70  filea.open("signal.txt");
71 
72  ofstream file1a;
73  file1a.open("bkg.txt");
74 
75 
76 ofstream file2a;
77  file2a.open("ncbkg.txt");
78  */
79  double pot_initial=0;
80  double count_pot=0;
81 
82  float cutbdtvallow=0;
83  float cutbdtvalhigh=0.45;
84 
85  int countBDTselected=0;
86 
87  int countTrueSignal=0;
88  int countTrueBkg=0;
89 
90  int countBDTselected_sb=0;
91 
92  int countTrueSignal_sb=0;
93  int countTrueBkg_sb=0;
94 
95  int countSignalBeforeRemid=0;
96  int countBkgBeforeRemid=0;
97  int countCCBeforeRemid=0;
98  int countnoCCBeforeRemid=0;
99  int countNCBeforeRemid=0;
100 
101  int countSignalBeforeReco=0;
102  int countBkgBeforeReco=0;
103  int countCCBeforeReco=0;
104  int countnoCCBeforeReco=0;
105  int countNCBeforeReco=0;
106 
107 
108  int countSignalBeforeME=0;
109  int countBkgBeforeME=0;
110  int countCCBeforeME=0;
111  int countnoCCBeforeME=0;
112  int countNCBeforeME=0;
113 
114 
115  int countSignalBeforeBdt=0;
116  int countBkgBeforeBdt=0;
117  int countCCBeforeBdt=0;
118  int countnoCCBeforeBdt=0;
119  int countNCBeforeBdt=0;
120 
121  int countSelCCpi0bkg_tot=0;
122  int countSelCCnopi0bkg_tot=0;
123  int countSelNCbkg_tot=0;
124 
125  int countSelCCpi0bkg_rej=0;
126  int countSelCCnopi0bkg_rej=0;
127  int countSelNCbkg_rej=0;
128  int countSelNCnopi0bkg_rej=0;
129 
130  int countSelCCpi0bkg_tot_sb=0;
131  int countSelCCnopi0bkg_tot_sb=0;
132  int countSelNCbkg_tot_sb=0;
133 
134  int countSelCCpi0bkg_rej_sb=0;
135  int countSelCCnopi0bkg_rej_sb=0;
136  int countSelNCbkg_rej_sb=0;
137 
138 
139  //-------------------------------
140 
141 
142  TH1F *selPiMass=new TH1F("selPimass","selPiMass",100,0.0,0.3);
143  TH1F *selSigPiMass=new TH1F("selSigPiMass","selSigPiMass",100,0.0,0.3);
144  TH1F *selBkgPiMass=new TH1F("selBkgPiMass","selBkgPiMass",100,0.0,0.3);
145  TH1F *selBkgPiMass_ccpi0=new TH1F("selBkgPimass_ccpi0","selBkgPiMass_ccpi0",100,0.0,0.3);
146  TH1F *selBkgPiMass_ccnopi0=new TH1F("selBkgPimass_ccnopi0","selBkgPiMass_ccnopi0",100,0.0,0.3);
147  TH1F *selBkgPiMass_ncpi0=new TH1F("selBkgPimass_ncpi0","selBkgPiMass_ncpi0",100,0.0,0.3);
148  TH1F *selBkgPiMass_ncnopi0=new TH1F("selBkgPimass_ncnopi0","selBkgPiMass_ncnopi0",100,0.0,0.3);
149 
150 
151 
152 
153  TH1F *totBdt=new TH1F("totBdt","totBdt",100,-1.0,1.0);
154  TH1F *selSigBDT=new TH1F("selSigBDT","selSigBDT",100,-1.0,1.0);
155  TH1F *selBkgBDT=new TH1F("selBkgBDT","selBkgBDT",100,-1.0,1.0);
156  TH1F *selBkgBDT_ccpi0=new TH1F("selBkgBDT_ccpi0","selBkgBDT_ccpi0",100,-1.0,1.0);
157  TH1F *selBkgBDT_ccnopi0=new TH1F("selBkgBDT_ccnopi0","selBkgBDT_ccnopi0",100,-1.0,1.0);
158  TH1F *selBkgBDT_nc=new TH1F("selBkgBDT_nc","selBkgBDT_nc",100,-1.0,1.0);
159 
160 
161 
162  SpectrumLoader loader(fin);
163 
164  // TFile outputFile(TString(fileout),"RECREATE");
165 
166  /// TTree *inTree = new TTree("inTree", "getting variables for BDT evaluation");
167 
168 // Debuging
169  // int Run;
170  // int SRun;
171  //int Evt;
172 
173  // Main variables
174  int IsNCNumu;
175  int IsNC;
176  int IsContained;
177  int IsFiducial;
178  int IsPi0all;
179  int IsTwoProngs;
180  int IsCC;
181  int inttype;
182  float michel;
183  int mode;
184 
185  float nue, nue_Vis, nue_Vis_slice;
186 
187  float IsRem;
188  float remid;
189  // float cvn;
190  float cvnncid, cvnnumu,cvnnue, cvnnutau, cvnphoton;
191  float sliceEn;
192  float TrueMuEn;
193 
194  float prong1MissingPl;
195  float prong2MissingPl;
196 
197  float prong1ContPl;
198  float prong2ContPl;
199 
200  float prong1Length;
201  float prong2Length;
202 
203  float prong1Energy;
204  float prong2Energy;
205  float HadronicEnergy;
206 
207  float prong1Width;
208  float prong2Width;
209 
210  float prong1VtxGev;
211 
212  float prong1dedx;
213  float prong2dedx;
214 
215  float prong1epLLL;
216  //float prong2epLLL;
217 
218  float prong1epLLT;
219  //float prong2epLLT;
220 
221  float prong1epiLLL;
222  //float prong2epiLLL;
223 
224  float prong1epiLLT;
225  //float prong2epiLLT;
226 
227 float prong1epi0LLL;
228 // float prong2epi0LLL;
229 
230  float prong1epi0LLT;
231  // float prong2epi0LLT;
232 
233  float png1dirx;
234  float png1diry;
235  float png1dirz;
236  float png2dirx;
237  float png2diry;
238  float png2dirz;
239 
240  TLorentzVector pi0p;
241  float pi0pz, pi0pt;
242 
243  float thetabeam_true, costhetabeam_true;
244  //not for training
245  float shwlidpi0mass;
246 float recopimass;
247 
248  int Run;
249  int Subrun;
250  int Event;
251  int Cycle;
252 int SliceNumber;
253  float distPngStartPos;
254 
255  gROOT->LoadMacro("weights/MVAnalysis_BDTG.class.C++");
256 
257  std::vector<string> inputVariableNames;
258 
259 
260  inputVariableNames.push_back("cvnphoton");
261  inputVariableNames.push_back("prong1epi0LLL");
262  inputVariableNames.push_back("prong1epiLLL");
263  inputVariableNames.push_back("prong1ContPl");
264  inputVariableNames.push_back("prong1epLLT");
265  // inputVariableNames.push_back("prong1epiLLT");
266  inputVariableNames.push_back("prong1Width");
267  inputVariableNames.push_back("prong1dedx");
268  inputVariableNames.push_back("prong1epLLL");
269  //inputVariableNames.push_back("prong1epi0LLT");
270  // inputVariableNames.push_back("recopimass");
271  inputVariableNames.push_back("prong1MissingPl");
272 
273 
274 
275 
276  IClassifierReader* classReader = new ReadBDTG(inputVariableNames);
277 
278 
279  IFileSource* filesrc=loader.WildcardOrSAMQuery(fin);
280 
281  int Nfiles = filesrc->NFiles();
282  std::cout<<"NUMBER OF FILES: "<<Nfiles<<std::endl;
283  Progress* prog = 0;
284  TFile* file1;
285  int fileIdx = -1;
286  while(TFile* file1 =(TFile*)filesrc->GetNextFile()){
287  ++fileIdx;
288  std::cout<<"*************** File NUMBER: "<<fileIdx<<std::endl;
289  if(Nfiles >= 0 && !prog) prog = new Progress(TString::Format("Filling tree from %d files", Nfiles).Data());
290 
291  loader.HandleFile(file1, Nfiles == 1 ? prog : 0);
292 
293  if(Nfiles > 1 && prog) prog->SetProgress((fileIdx+1.)/Nfiles);
294 
295  //******************************************
296 
297  // TFile *file1 = new TFile(f,"READ");
298 
299  if( !(file1) || (file1->IsZombie())){
300  std::cout << " --- INPUT file doesn't exists OR Zombie!!" << std::endl;
301  //return -1;
302  }
303 
304 
305 
306  TH1D *h_pot = (TH1D*)file1->Get("TotalPOT");
307  count_pot=h_pot->Integral();
308  pot_initial=pot_initial+count_pot;
309 
310 
311  //-------------------------------
312  // Defining variables
313 
314  // Defining the TTree object
315  TTree *recTree = (TTree*)file1->Get("recTree");
316 
317  // Every time a GetEntry(i) is done, the data in the tree is copied to this object
318  caf::StandardRecord* recTreeObject = 0;
319  recTree->SetBranchAddress("rec", &recTreeObject);
320 
321  // Do not load any branches into the memory yet
322  recTree->SetBranchStatus("*",0);
323 
324  // Header chain
325  recTree->SetBranchStatus("hdr.run",1);
326  recTree->SetBranchStatus("hdr.subrun",1);
327  recTree->SetBranchStatus("hdr.evt",1);
328  recTree->SetBranchStatus("hdr.subevt",1);
329  recTree->SetBranchStatus("hdr.cycle",1);
330 
331 
332  // MC Chain
333  recTree->SetBranchStatus("mc.nnu",1);
334  recTree->SetBranchStatus("mc.nu",1);
335 
336  recTree->SetBranchStatus("mc.nu.E",1);
337  recTree->SetBranchStatus("mc.nu.visE",1);
338  recTree->SetBranchStatus("mc.nu.visEinslc",1);
339  recTree->SetBranchStatus("mc.nu.iscc",1);
340  recTree->SetBranchStatus("mc.nu.inttype",1);
341  recTree->SetBranchStatus("mc.nu.mode",1);
342  recTree->SetBranchStatus("mc.nu.pdg",1);
343  // recTree->SetBranchStatus("mc.nu.pdgorig",1 );
344  recTree->SetBranchStatus("mc.nu.prim",1);
345  recTree->SetBranchStatus("mc.nu.prim.pdg",1);
346  recTree->SetBranchStatus("mc.nu.prim.p",1);
347  recTree->SetBranchStatus("mc.nu.prim.p.E");
348  // recTree->SetBranchStatus("mc.nu.prim.p.fE");
349  // recTree->SetBranchStatus("mc.nu.prim.p.P",1);
350  // recTree->SetBranchStatus("mc.nu.prim.p.Pt()",1);
351  recTree->SetBranchStatus("mc.nu.prim.p.pz",1);
352 
353  // Vtx Chain
354  recTree->SetBranchStatus("vtx.elastic",1);
355  recTree->SetBranchStatus("vtx.nelastic",1);
356  recTree->SetBranchStatus("vtx.elastic.vtx",1);
357  recTree->SetBranchStatus("vtx.elastic.vtx.x",1);
358  recTree->SetBranchStatus("vtx.elastic.vtx.y",1);
359  recTree->SetBranchStatus("vtx.elastic.vtx.z",1);
360 
361  // ShwLid Chain
362  recTree->SetBranchStatus("vtx.elastic.fuzzyk.npng2d",1);
363  recTree->SetBranchStatus("vtx.elastic.fuzzyk.nshwlid",1);
364  recTree->SetBranchStatus("vtx.elastic.fuzzyk.npng",1);
365  recTree->SetBranchStatus("vtx.elastic.fuzzyk.npng2d",1);
366  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png",1);
367  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d",1);
368  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid",1);
369  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop",1);
370  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop.x",1);
371  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop.y",1);
372  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.stop.z",1);
373  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start",1);
374  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start.x",1);
375  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start.y",1);
376  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.start.z",1);
377  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d.start",1);
378  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d.start.x",1);
379  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d.start.y",1);
380  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d.start.z",1);
381  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d.len",1);
382  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png2d.calE",1);
383 
384 
385  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.calE",1);
386  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.width",1);
387  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.len",1);
388  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.maxplanecont",1);
389  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.maxplanegap",1);
390 
391  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir",1);
392  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir.x",1);
393  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir.y",1);
394  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.dir.z",1);
395 
396 
397 //Longitudinal and Transverse Likelihoods
398  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid",1);
399  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.eplll",1);
400  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epllt",1);
401  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epilll",1);
402  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epillt",1);
403  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epi0lll",1);
404  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.epi0llt",1);
405 
406  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.vtxgev",1);
407  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.shwlid.lid.pi0mass",1);
408 
409  recTree->SetBranchStatus("sand.nue.dedxpng1",1);
410  recTree->SetBranchStatus("sand.nue.dedxpng2",1);
411 
412  // Slice Chain
413  recTree->SetBranchStatus("slc.calE",1);
414  recTree->SetBranchStatus("slc.nhit",1);
415 
416  // Selection Chain
417  recTree->SetBranchStatus("sel.remid.pid",1);
418  recTree->SetBranchStatus("sel.cvn.ncid",1);
419  recTree->SetBranchStatus("sel.cvn.numuid",1);
420  recTree->SetBranchStatus("sel.cvn.nueid",1);
421  recTree->SetBranchStatus("sel.cvn.nutauid",1);
422  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.cvnpart.photonid",1);
423 
424 
425 
426  // recTree->SetBranchStatus("vtx.elastic.vtx.X()",1);
427  //recTree->SetBranchStatus("vtx.elastic.vtx.Y()",1);
428  //recTree->SetBranchStatus("vtx.elastic.vtx.Z()",1);
429  // recTree->SetBranchStatus("vtx.elastic.npng3d",1);
430 
431 
432 
433  //recTree->SetBranchStatus("shw.shwlid.stop.X()",1);
434  //recTree->SetBranchStatus("shw.shwlid.stop.Y()",1);
435  //recTree->SetBranchStatus("shw.shwlid.stop.Z()",1);
436 
437 //Truth Matching
438  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth",1);
439  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.pdg",1);
440  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.motherpdg",1);
441  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.eff",1);
442  recTree->SetBranchStatus("vtx.elastic.fuzzyk.png.truth.pur",1);
443 
444 
445  //-------------------------------
446  // Looping over entries
447 
448  //Int_t numberOfEntries = recTree->GetEntries();
449  Int_t numberOfEntries = recTree->GetEntriesFast();
450  for (Int_t event = 0; event < numberOfEntries; ++event) {
451 
452  recTree->GetEntry(event);
453  cerr << "\r-- Processing event " << event << " of " << numberOfEntries;
454 
455  //-------------------------------
456  // 2 prongs cut
457  // Only slices - entries with 2 prongs
458  //std::cout<<"MC: "<<recTreeObject->mc.nnu<<std::endl;
459 
460  if(recTreeObject->mc.nnu == 0) continue;
461  if(recTreeObject->vtx.nelastic == 0) continue;
462 
463  Int_t ProngNumber = recTreeObject->vtx.elastic[0].fuzzyk.npng;
464  Int_t ShowerNumber = recTreeObject->vtx.elastic[0].fuzzyk.nshwlid;
465  Int_t ProngNumber2d = recTreeObject->vtx.elastic[0].fuzzyk.npng2d;
466  Int_t nbOfVtx = recTreeObject->vtx.nelastic;
467 
468  //std::cout<<"Number of Shw: "<<nbOfShw<<std::endl;
469  if(ProngNumber2d==0) continue;
470  if(ProngNumber==0) continue;
471  if(ShowerNumber<=0) continue;
472 
473  if(ProngNumber !=1) continue;
474  if(nbOfVtx !=1) continue;
475 
476  // std::cout << "Nb Of Prongs" << ProngNumber << std::endl;
477 
478  //-------------------------------
479  // Containment and fiducial cut
480  Bool_t contained = true;
481  Bool_t fiducial = false;
482 
483 if((recTreeObject->vtx.elastic[0].vtx.x) < 170.0
484  && (recTreeObject->vtx.elastic[0].vtx.x) > -170.0
485  && (recTreeObject->vtx.elastic[0].vtx.y) < 170.0
486  && (recTreeObject->vtx.elastic[0].vtx.y) < -160.0
487  && recTreeObject->vtx.elastic[0].vtx.z < 1100.0
488  && recTreeObject->vtx.elastic[0].vtx.z > 50.0) fiducial=true;
489 
490 
491  IsFiducial = fiducial;
492  if(!IsFiducial) continue;
493 
494  for(int i=0;i<ShowerNumber;i++)
495  {
496  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) > 125) contained=false;
497  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) < -100) contained=false;
498  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) > 100) contained=false;
499  if((recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) < -120) contained=false;
500  if(recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z > 1200) contained=false;
501  if(recTreeObject->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z < 300) contained=false;
502  }
503 
504 
505  IsContained = contained;
506  if(!IsContained) continue;
507 
508 
509  float SlcHits;
510 SlcHits = recTreeObject->slc.nhit;
511 
512 
513  //-------------------------------
514  // Numu NC and remid cut to select signal:
515 
516  Bool_t ncnumu = false;
517  if(recTreeObject->mc.nu[0].pdg==14 && !recTreeObject->mc.nu[0].iscc){ncnumu = true;}
518  IsNCNumu = ncnumu;
519 
520  Bool_t isnc = false;
521  if(!recTreeObject->mc.nu[0].iscc)
522  {
523  isnc = true;
524  }
525  IsNC = isnc;
526 
527 
528  //-------------------------------
529  // primaryPi0 cut
530  Bool_t ispi0 = false;
531  float t=-5.0;
532  float e=-10.0;
533  float pi0mass=0.135;
534  int nbOfPi=0;
535  float PiPrimary=0;
536  float maximum=0;
537  int nprim=recTreeObject->mc.nu[0].prim.size();
538  for(int i=0;i<nprim;i++){
539  if(recTreeObject->mc.nu[0].prim[i].pdg==111){
540  PiPrimary++;
541  e=recTreeObject->mc.nu[0].prim[i].p.E;
542  //std::cout<<"Pi0 energy: "<<e<<std::endl;
543  t=e - pi0mass;
544  if(t>maximum) maximum=t;
545  if(t>0.0) nbOfPi++;
546  }
547  }
548  if(nbOfPi>0) ispi0=true;
549 
550  IsPi0all=ispi0;
551 
552  IsCC = recTreeObject->mc.nu[0].iscc;
553 
554 
555  if(ProngNumber !=1) continue;
556  if(nbOfVtx !=1) continue;
557 
558 
559 
560  if (IsNC && IsPi0all){
561  countSignalBeforeRemid++;
562 
563 
564  }else{
565  countBkgBeforeRemid++;
566 
567 
568  if(IsCC && IsPi0all) {
569  countCCBeforeRemid++;
570  }
571  if(IsCC && !IsPi0all){
572  countnoCCBeforeRemid++;
573  }
574  if(!IsCC) {
575  countNCBeforeRemid++;
576  }
577 
578  }
579 
580 
581  //RemID cut
582 
583  remid=recTreeObject->sel.remid.pid;
584 
585  if (remid>0.375) continue;
586 
587 
588 
589 
590 
591  if (IsNC && IsPi0all){
592  countSignalBeforeReco++;
593  }else{
594  countBkgBeforeReco++;
595 
596  if(IsCC && IsPi0all) countCCBeforeReco++;
597  if(IsCC && !IsPi0all) countnoCCBeforeReco++;
598  if(!IsCC) countNCBeforeReco++;
599 
600 
601  }
602 
603 
604 
605  // For Reconstructed K.E.:
606  //I'm changing this to use slice energy but leaving the name as reco energy because I don't want to change the names all over the script.
607 
608  float RecoEnergy_truepimass=-5.0;
609  float RecoEnergy_recopimass=-5.0;
610  float pi0mass1=0.135;
611  float totEn=0;
612  float SliceEnergy=0;
613  float OtherEnergy=0;
614 
615  //totEn= recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.calE;// + recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.calE;
616  //totEn= prong1Energy+prong2Energy;
617 
618  SliceEnergy = recTreeObject->slc.calE;
619 
620  RecoEnergy_truepimass=SliceEnergy;
621 
622 
623  if (RecoEnergy_truepimass<0.0) continue;
624 
625 
626 
627 
628  if (IsNC && IsPi0all){
629  countSignalBeforeBdt++;
630 
631 
632  }else{
633  countBkgBeforeBdt++;
634 
635 
636 
637  if(IsCC && IsPi0all) {
638  countCCBeforeBdt++;
639 
640  }
641  if(IsCC && !IsPi0all){
642  countnoCCBeforeBdt++;
643 
644 
645  }
646  if(!IsCC) {
647  countNCBeforeBdt++;
648 
649 
650  }
651 
652  }
653 
654 
655 
656  //----------
657  // Run = recTreeObject->hdr.run;
658  // SRun = recTreeObject->hdr.subrun;
659  // Evt = recTreeObject->hdr.evt;
660  // SEvt = recTreeObject->hdr.subevt;
661 
662  //----------
663  // Main variables
664 
665  // Defining nature of samples
666  Short_t mc_nnu = recTreeObject->mc.nnu;
667  // Bool_t mc_nu_iscc = recTreeObject->mc.nu[0].iscc;
668  Short_t mc_nu_pdg = recTreeObject->mc.nu[0].pdg;
669 
670 
671  inttype= recTreeObject->mc.nu[0].inttype;
672  mode= recTreeObject->mc.nu[0].mode;
673  nue=recTreeObject->mc.nu[0].E;
674  nue_Vis=recTreeObject->mc.nu[0].visE;
675  nue_Vis_slice=recTreeObject->mc.nu[0].visEinslc;
676 
677 
678 
679  remid=recTreeObject->sel.remid.pid;
680  prong1MissingPl=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.maxplanegap;
681  prong1ContPl=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.maxplanecont;
682  prong1Length=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.len;
683  prong1Energy=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.calE;
684  prong2Energy=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.calE;
685  prong1Width=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.width;
686  prong1VtxGev=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.vtxgev;
687  prong1dedx=recTreeObject->sand.nue.dedxpng1;
688  prong2dedx=recTreeObject->sand.nue.dedxpng2;
689  prong1epLLL=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.eplll;
690  prong1epLLT=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epllt;
691  prong1epiLLL=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epilll;
692  prong1epiLLT=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epillt;
693  prong1epi0LLL=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epi0lll;
694  prong1epi0LLT=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.epi0llt;
695  png1dirx=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.x;
696  png1diry=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.y;
697  png1dirz=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.z;
698  png2dirx=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.x;
699  png2diry=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.y;
700  png2dirz=recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.z;
701  cvnncid=recTreeObject->sel.cvn.ncid;
702  cvnnumu=recTreeObject->sel.cvn.numuid;
703  cvnnue=recTreeObject->sel.cvn.nueid;
704  cvnnutau=recTreeObject->sel.cvn.nutauid;
705  cvnphoton=recTreeObject->vtx.elastic[0].fuzzyk.png[0].cvnpart.photonid;
706  Run = recTreeObject->hdr.run;
707  Subrun = recTreeObject->hdr.subrun;
708  Event = recTreeObject->hdr.evt;
709  Cycle = recTreeObject->hdr.cycle;
710  SliceNumber = recTreeObject->hdr.subevt;
711 
712 
713 
714 
715  // For Shw lid mass
716 
717  float masstrue=1000;
718 
719  //For now I'm going to make this into true neutrino energy. Again I don't want to have to change everything from here to the end so I'm going to leave the name the same.
720  // recopimass=recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.lid.pi0mass;
721  recopimass=recTreeObject->mc.nu[0].E;
722 
723 
724 
725 
726  /*
727  // For Reconstructed Mass:
728  float thispi0mass1;
729  float mass1;
730  float dotproduct; // here dotproduct is Cos theta as it is dot product of unit vectors and it is opening angle between two prongs
731 
732  dotproduct= recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.dir.Unit().Dot(recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.dir.Unit());
733  thispi0mass1=sqrt(recTreeObject->vtx.elastic[0].fuzzyk.png[0].shwlid.calE*recTreeObject->vtx.elastic[0].fuzzyk.png[1].shwlid.calE*(1-dotproduct));
734 
735 
736 
737  recopimass=thispi0mass1;
738 
739  */
740  //For RecoKE with reco mass
741  RecoEnergy_recopimass= recopimass;
742 
743  sliceEn=recTreeObject->slc.calE;
744 
745  // OtherEnergy=sliceEn- prong1Energy;
746  HadronicEnergy=sliceEn - totEn;
747 
748 
749 
750 
751 
752  std::vector<double> inputVariableValues;
753 
754 
755  inputVariableValues.push_back(cvnphoton);
756  inputVariableValues.push_back(prong1epi0LLL);
757  inputVariableValues.push_back(prong1epiLLL);
758  inputVariableValues.push_back(prong1ContPl);
759  inputVariableValues.push_back(prong1epLLT);
760  // inputVariableValues.push_back(prong1epiLLT);
761  inputVariableValues.push_back(prong1Width);
762  inputVariableValues.push_back(prong1dedx);
763  inputVariableValues.push_back(prong1epLLL);
764  //inputVariableValues.push_back(prong1epi0LLT);
765  // inputVariableValues.push_back(recopimass);
766  inputVariableValues.push_back(prong1MissingPl);
767 
768 
769 
770 
771  //std::cout << std::endl;
772  // std::cout <<" ==>Start TMVA BDT for NCpi0 selection... " << std::endl;
773 
774  Double_t bdtVal= classReader->GetMvaValue(inputVariableValues);
775  totBdt->Fill(bdtVal);
776 
777 
778 
779  float selSignalPimass=-100;
780  float selBkgPimass=-100;
781  float trueSignalPimass=-100;
782  float trueBkgPimass=-100;
783  //cut value from FOM max
784  // std::cout<<"The most probable pimass is: "<<shwlidpi0mass<<std::endl;
785  // if(bdtVal>0.5457){
786 
787 
788  if(bdtVal > 0.600){
789 
790  countBDTselected++;
791  // std::cout<<"The BDT value for this slice is: "<<bdtVal<<std::endl;
792 
793 
794 
795  selPiMass->Fill(recopimass);
796 
797  if(IsNC && IsPi0all){
798 
799  countTrueSignal++;
800 
801 
802  // filea << Run << "\t" << Subrun << "\t" << Cycle << "\t" << Event << "\t" << SliceNumber << "\n";
803 
804 
805  selSigPiMass->Fill(recopimass);
806 
807 
808  }else{
809  countTrueBkg++;
810 
811 
812  selBkgPiMass->Fill(recopimass);
813 
814 
815  if(IsCC && IsPi0all){
816 
817  countSelCCpi0bkg_rej++;
818 
819 
820  selBkgPiMass_ccpi0->Fill(recopimass);
821 
822 
823  }
824 
825  if(IsCC && !IsPi0all)
826  {
827 countSelCCnopi0bkg_rej++;
828 
829  selBkgPiMass_ccnopi0->Fill(recopimass);
830 
831 
832  }
833  if(!IsCC && IsPi0all){
834 
835  countSelNCbkg_rej++;
836 
837 
838  selBkgPiMass_ncpi0->Fill(recopimass);
839 
840  }
841 
842  if(!IsCC && !IsPi0all){
843 
844  countSelNCnopi0bkg_rej++;
845 
846  selBkgPiMass_ncnopi0->Fill(recopimass);
847 
848 
849 }
850 
851  }// end for bkg cut
852  } // end for bdt cut
853 
854 
855  // bkgPiMass->Fill(selBkgPimass);
856 
857  //Totals from MC truth
858  if(IsNC && IsPi0all){
859  // std::cout<<"This is a true signal event"<<std::endl;
860 
861  selSigBDT->Fill(bdtVal);
862 
863  }else{
864 
865 
866 
867  selBkgBDT->Fill(bdtVal);
868  //bkgNue->Fill(nue);
869 
870 
871  if(IsCC && IsPi0all){
872  // selBkgBDT_ccpi0->Fill(bdtVal);
873  countSelCCpi0bkg_tot++;
874 
875  selBkgBDT_ccpi0->Fill(bdtVal);
876 
877 
878  }
879 
880  if(IsCC && !IsPi0all) {
881  // selBkgBDT_ccnopi0->Fill(bdtVal);
882  countSelCCnopi0bkg_tot++;
883 
884  selBkgBDT_ccnopi0->Fill(bdtVal);
885 
886 
887 }
888  if(!IsCC){
889  // selBkgBDT_nc->Fill(bdtVal);
890  countSelNCbkg_tot++;
891  selBkgBDT_nc->Fill(bdtVal);
892 
893 }
894  }
895 
896 
897  } // end loop over events
898  recTree->Delete();
899  }// end for fileIdx
900 
901 
902  if(prog){
903  prog->Done();
904  delete prog;
905  }
906 
907 
908 
909  TFile* outputFile = TFile::Open(TString(fout), "RECREATE" );
910  std::cout<<std::endl;
911 
912  std::cout<<"***************Counting Events with Fiducial and Containment: " <<std::endl;
913  std::cout<<"Signal Before ReMID: "<< countSignalBeforeRemid << " Bkg: "<<countBkgBeforeRemid << " CC: "<< countCCBeforeRemid<< " NoCC: "<< countnoCCBeforeRemid<<" NC: "<< countNCBeforeRemid<<std::endl;
914 
915  std::cout<<"\n***************Counting Events with Fiducial and Containment and ReMID: " <<std::endl;
916  std::cout<<"Signal Before Slice Energy: "<< countSignalBeforeReco << "Bkg: "<<countBkgBeforeReco << "CC: "<< countCCBeforeReco<< "NoCC: "<< countnoCCBeforeReco<<"NC: "<< countNCBeforeReco<<std::endl;
917 
918  //std::cout<<"***************Counting Events with Fiducial and Containment, ReMID and Reco KE: " <<std::endl;
919  //std::cout<<"Signal Before ME: "<< countSignalBeforeME << "Bkg: "<<countBkgBeforeME << "CC: "<< countCCBeforeME<< "NoCC: "<< countnoCCBeforeME<<"NC: "<< countNCBeforeME<<std::endl;
920 
921 
922  std::cout<<"\n***************Counting Events with Fiducial and Containment, ReMID, and Slice Energy: " <<std::endl;
923  std::cout<<"Signal Before BDT: "<< countSignalBeforeBdt << "Bkg: "<<countBkgBeforeBdt << "CC: "<< countCCBeforeBdt<< "NoCC: "<< countnoCCBeforeBdt<<"NC: "<< countNCBeforeBdt<<std::endl;
924 
925 
926  std::cout<<"\n*********** Counting events: "<<std::endl;
927  std::cout<<"Selected signals with bdt: "<<countBDTselected<<std::endl;
928  std::cout<<"True signal: "<< countTrueSignal<<" it was bkg: "<< countTrueBkg<<std::endl;
929  std::cout<<"Bkg decomposition: CCpi0: "<< countSelCCpi0bkg_rej << " CCnopi0: "<< countSelCCnopi0bkg_rej <<" NC: "<<countSelNCbkg_rej<<std::endl;
930 
931  std::cout<<"\n*********** Number of events from MC truth: "<<std::endl;
932  std::cout<<"Signal: "<< selSigBDT->Integral()<<" Bkg: "<< selBkgBDT->Integral()<<std::endl;
933  std::cout<<"Bkg decomposition: CCpi0: "<< countSelCCpi0bkg_tot << " CCnopi0: "<< countSelCCnopi0bkg_tot <<" NC: "<<countSelNCbkg_tot<<std::endl;
934 
935 
936 
937 
938 
939  std::cout << "POT count: " << pot_initial << std::endl;
940 
941 
942 totBdt->Write();
943 selSigBDT->Write();
944 selBkgBDT->Write();
945 selBkgBDT_ccpi0->Write();
946 selBkgBDT_ccnopi0->Write();
947 selBkgBDT_nc->Write();
948 
949 
950 
951 selPiMass->Write();
952 selSigPiMass->Write();
953 selBkgPiMass->Write();
954 selBkgPiMass_ccpi0->Write();
955 selBkgPiMass_ccnopi0->Write();
956 selBkgPiMass_ncpi0->Write();
957 selBkgPiMass_ncnopi0->Write();
958 
959 
960 
961 //signalPiMass_sb->Write();
962 //truebkgPiMass_sb->Write();
963 
964  outputFile->Write();
965  outputFile->Close();
966 
967  std::cout << "--- Created root file: \"BDTTest.root\" containing the MVA output histograms" << std::endl;
968 
969  // delete reader;
970 
971  std::cout << "==> TMVAClassification is done!" << endl << std::endl;
972 
973 
974  std::cout<<"**** End **** "<<std::endl;
975  //return 0;
976 }
977 
978 #endif
caf::SRCVNResult::ncid
float ncid
Likelihood Neutral Current.
Definition: SRCVNResult.h:23
fin
TString fin
Definition: Style.C:24
make_root_from_grid_output.ispi0
ispi0
Definition: make_root_from_grid_output.py:96
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
caf::StandardRecord::hdr
SRHeader hdr
Header branch: run, subrun, etc.
Definition: StandardRecord.h:40
caf::SRHeader::subrun
unsigned int subrun
subrun number
Definition: SRHeader.h:22
caf::SRFuzzyK::png
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
submit_syst.fout
fout
Definition: submit_syst.py:244
make_root_from_grid_output.remid
remid
Definition: make_root_from_grid_output.py:79
om::cerr
OStream cerr
Definition: OStream.cxx:7
caf::SRHeader::run
unsigned int run
run number
Definition: SRHeader.h:21
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
caf::SRHeader::evt
unsigned int evt
ART event number, indexes trigger windows.
Definition: SRHeader.h:25
submit_nova_art.mode
mode
Definition: submit_nova_art.py:725
useBDTG
void useBDTG()
Definition: useBDTG.C:63
caf::SRRemid::pid
float pid
PID value output by kNN.
Definition: SRRemid.h:25
ana::IFileSource::NFiles
virtual int NFiles() const
May return -1 indicating the number of files is not known.
Definition: IFileSource.h:21
caf::SRFuzzyK::nshwlid
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
IClassifierReader::GetMvaValue
virtual double GetMvaValue(const std::vector< double > &inputValues) const =0
caf::SRIDBranch::cvn
SRCVNResult cvn
Horrible hack to appease CAFAna.
Definition: SRIDBranch.h:53
MVAnalysis_BDTG.class.C
ana::IsFiducial
bool IsFiducial(const TVector3 &nuVtx, const TVector3 &min, const TVector3 &max)
Definition: Flux.cxx:46
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
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
caf::SRFuzzyK::npng2d
size_t npng2d
Definition: SRFuzzyK.h:27
make_template_knob_config.recTree
recTree
Definition: make_template_knob_config.py:87
ana::IFileSource::GetNextFile
virtual TFile * GetNextFile()=0
Returns the next file in sequence, ready for reading.
caf::StandardRecord
The StandardRecord is the primary top-level object in the Common Analysis File trees.
Definition: StandardRecord.h:33
ana::IFileSource
Interface class for accessing ROOT files in sequence.
Definition: IFileSource.h:10
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
ana::SpectrumLoaderBase::WildcardOrSAMQuery
IFileSource * WildcardOrSAMQuery(const std::string &str) const
Figure out if str is a wildcard or SAM query and return a source.
Definition: SpectrumLoaderBase.cxx:72
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
caf::SRHeader::subevt
unsigned short subevt
slice number within spill
Definition: SRHeader.h:26
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::SRHeader::cycle
int cycle
MC simulation cycle number.
Definition: SRHeader.h:23
caf::StandardRecord::vtx
SRVertexBranch vtx
Vertex branch: location, time, etc.
Definition: StandardRecord.h:44
string
enum BeamMode string
Generated on Sun Feb 28 2021 14:13:14 for NOvA by   doxygen 1.8.11