18 #include "CAFAna/Core/Binning.h" 49 #include "TStopwatch.h" 53 #if not defined(__CINT__) || defined(__MAKECINT__) 54 #include "TMVA/Tools.h" 55 #include "TMVA/Reader.h" 56 #include "TMVA/MethodCuts.h" 65 const std::string fin =
"defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
83 float cutbdtvalhigh=0.45;
85 int countBDTselected=0;
87 int countTrueSignal=0;
90 int countBDTselected_sb=0;
92 int countTrueSignal_sb=0;
93 int countTrueBkg_sb=0;
95 int countSignalBeforeRemid=0;
96 int countBkgBeforeRemid=0;
97 int countCCBeforeRemid=0;
98 int countnoCCBeforeRemid=0;
99 int countNCBeforeRemid=0;
101 int countSignalBeforeReco=0;
102 int countBkgBeforeReco=0;
103 int countCCBeforeReco=0;
104 int countnoCCBeforeReco=0;
105 int countNCBeforeReco=0;
108 int countSignalBeforeME=0;
109 int countBkgBeforeME=0;
110 int countCCBeforeME=0;
111 int countnoCCBeforeME=0;
112 int countNCBeforeME=0;
115 int countSignalBeforeBdt=0;
116 int countBkgBeforeBdt=0;
117 int countCCBeforeBdt=0;
118 int countnoCCBeforeBdt=0;
119 int countNCBeforeBdt=0;
121 int countSelCCpi0bkg_tot=0;
122 int countSelCCnopi0bkg_tot=0;
123 int countSelNCbkg_tot=0;
125 int countSelCCpi0bkg_rej=0;
126 int countSelCCnopi0bkg_rej=0;
127 int countSelNCbkg_rej=0;
128 int countSelNCnopi0bkg_rej=0;
130 int countSelCCpi0bkg_tot_sb=0;
131 int countSelCCnopi0bkg_tot_sb=0;
132 int countSelNCbkg_tot_sb=0;
134 int countSelCCpi0bkg_rej_sb=0;
135 int countSelCCnopi0bkg_rej_sb=0;
136 int countSelNCbkg_rej_sb=0;
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);
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);
185 float nue, nue_Vis, nue_Vis_slice;
190 float cvnncid, cvnnumu,cvnnue, cvnnutau, cvnphoton;
194 float prong1MissingPl;
195 float prong2MissingPl;
205 float HadronicEnergy;
243 float thetabeam_true, costhetabeam_true;
253 float distPngStartPos;
255 gROOT->LoadMacro(
"weights/MVAnalysis_BDTG.class.C++");
257 std::vector<string> inputVariableNames;
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");
266 inputVariableNames.push_back(
"prong1Width");
267 inputVariableNames.push_back(
"prong1dedx");
268 inputVariableNames.push_back(
"prong1epLLL");
271 inputVariableNames.push_back(
"prong1MissingPl");
281 int Nfiles = filesrc->
NFiles();
286 while(TFile* file1 =(TFile*)filesrc->
GetNextFile()){
291 loader.
HandleFile(file1, Nfiles == 1 ? prog : 0);
293 if(Nfiles > 1 && prog) prog->
SetProgress((fileIdx+1.)/Nfiles);
299 if( !(file1) || (file1->IsZombie())){
306 TH1D *h_pot = (TH1D*)file1->Get(
"TotalPOT");
307 count_pot=h_pot->Integral();
308 pot_initial=pot_initial+count_pot;
315 TTree *
recTree = (TTree*)file1->Get(
"recTree");
319 recTree->SetBranchAddress(
"rec", &recTreeObject);
322 recTree->SetBranchStatus(
"*",0);
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);
333 recTree->SetBranchStatus(
"mc.nnu",1);
334 recTree->SetBranchStatus(
"mc.nu",1);
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);
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");
351 recTree->SetBranchStatus(
"mc.nu.prim.p.pz",1);
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);
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);
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);
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);
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);
406 recTree->SetBranchStatus(
"vtx.elastic.fuzzyk.png.shwlid.lid.vtxgev",1);
407 recTree->SetBranchStatus(
"vtx.elastic.fuzzyk.png.shwlid.lid.pi0mass",1);
409 recTree->SetBranchStatus(
"sand.nue.dedxpng1",1);
410 recTree->SetBranchStatus(
"sand.nue.dedxpng2",1);
413 recTree->SetBranchStatus(
"slc.calE",1);
414 recTree->SetBranchStatus(
"slc.nhit",1);
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);
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);
449 Int_t numberOfEntries = recTree->GetEntriesFast();
450 for (Int_t event = 0;
event < numberOfEntries; ++event) {
452 recTree->GetEntry(event);
453 cerr <<
"\r-- Processing event " <<
event <<
" of " << numberOfEntries;
460 if(recTreeObject->
mc.
nnu == 0)
continue;
461 if(recTreeObject->
vtx.nelastic == 0)
continue;
466 Int_t nbOfVtx = recTreeObject->
vtx.nelastic;
469 if(ProngNumber2d==0)
continue;
470 if(ProngNumber==0)
continue;
471 if(ShowerNumber<=0)
continue;
473 if(ProngNumber !=1)
continue;
474 if(nbOfVtx !=1)
continue;
480 Bool_t contained =
true;
492 if(!IsFiducial)
continue;
494 for(
int i=0;
i<ShowerNumber;
i++)
505 IsContained = contained;
506 if(!IsContained)
continue;
510 SlcHits = recTreeObject->
slc.
nhit;
516 Bool_t ncnumu =
false;
517 if(recTreeObject->
mc.
nu[0].pdg==14 && !recTreeObject->
mc.
nu[0].iscc){ncnumu =
true;}
521 if(!recTreeObject->
mc.
nu[0].iscc)
530 Bool_t
ispi0 =
false;
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){
541 e=recTreeObject->
mc.
nu[0].prim[
i].p.E;
544 if(t>maximum) maximum=
t;
548 if(nbOfPi>0) ispi0=
true;
552 IsCC = recTreeObject->
mc.
nu[0].iscc;
555 if(ProngNumber !=1)
continue;
556 if(nbOfVtx !=1)
continue;
560 if (IsNC && IsPi0all){
561 countSignalBeforeRemid++;
565 countBkgBeforeRemid++;
568 if(IsCC && IsPi0all) {
569 countCCBeforeRemid++;
571 if(IsCC && !IsPi0all){
572 countnoCCBeforeRemid++;
575 countNCBeforeRemid++;
585 if (remid>0.375)
continue;
591 if (IsNC && IsPi0all){
592 countSignalBeforeReco++;
594 countBkgBeforeReco++;
596 if(IsCC && IsPi0all) countCCBeforeReco++;
597 if(IsCC && !IsPi0all) countnoCCBeforeReco++;
598 if(!IsCC) countNCBeforeReco++;
608 float RecoEnergy_truepimass=-5.0;
609 float RecoEnergy_recopimass=-5.0;
610 float pi0mass1=0.135;
618 SliceEnergy = recTreeObject->
slc.
calE;
620 RecoEnergy_truepimass=SliceEnergy;
623 if (RecoEnergy_truepimass<0.0)
continue;
628 if (IsNC && IsPi0all){
629 countSignalBeforeBdt++;
637 if(IsCC && IsPi0all) {
641 if(IsCC && !IsPi0all){
642 countnoCCBeforeBdt++;
666 Short_t mc_nnu = recTreeObject->
mc.
nnu;
668 Short_t mc_nu_pdg = recTreeObject->
mc.
nu[0].pdg;
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;
687 prong1dedx=recTreeObject->sand.nue.dedxpng1;
688 prong2dedx=recTreeObject->sand.nue.dedxpng2;
706 Run = recTreeObject->
hdr.
run;
708 Event = recTreeObject->
hdr.
evt;
721 recopimass=recTreeObject->
mc.
nu[0].E;
741 RecoEnergy_recopimass= recopimass;
743 sliceEn=recTreeObject->
slc.
calE;
746 HadronicEnergy=sliceEn - totEn;
752 std::vector<double> inputVariableValues;
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);
761 inputVariableValues.push_back(prong1Width);
762 inputVariableValues.push_back(prong1dedx);
763 inputVariableValues.push_back(prong1epLLL);
766 inputVariableValues.push_back(prong1MissingPl);
774 Double_t bdtVal= classReader->
GetMvaValue(inputVariableValues);
775 totBdt->Fill(bdtVal);
779 float selSignalPimass=-100;
780 float selBkgPimass=-100;
781 float trueSignalPimass=-100;
782 float trueBkgPimass=-100;
795 selPiMass->Fill(recopimass);
797 if(IsNC && IsPi0all){
805 selSigPiMass->Fill(recopimass);
812 selBkgPiMass->Fill(recopimass);
815 if(IsCC && IsPi0all){
817 countSelCCpi0bkg_rej++;
820 selBkgPiMass_ccpi0->Fill(recopimass);
825 if(IsCC && !IsPi0all)
827 countSelCCnopi0bkg_rej++;
829 selBkgPiMass_ccnopi0->Fill(recopimass);
833 if(!IsCC && IsPi0all){
838 selBkgPiMass_ncpi0->Fill(recopimass);
842 if(!IsCC && !IsPi0all){
844 countSelNCnopi0bkg_rej++;
846 selBkgPiMass_ncnopi0->Fill(recopimass);
858 if(IsNC && IsPi0all){
861 selSigBDT->Fill(bdtVal);
867 selBkgBDT->Fill(bdtVal);
871 if(IsCC && IsPi0all){
873 countSelCCpi0bkg_tot++;
875 selBkgBDT_ccpi0->Fill(bdtVal);
880 if(IsCC && !IsPi0all) {
882 countSelCCnopi0bkg_tot++;
884 selBkgBDT_ccnopi0->Fill(bdtVal);
891 selBkgBDT_nc->Fill(bdtVal);
909 TFile* outputFile = TFile::Open(TString(fout),
"RECREATE" );
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;
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;
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;
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;
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;
945 selBkgBDT_ccpi0->Write();
946 selBkgBDT_ccnopi0->Write();
947 selBkgBDT_nc->Write();
952 selSigPiMass->Write();
953 selBkgPiMass->Write();
954 selBkgPiMass_ccpi0->Write();
955 selBkgPiMass_ccnopi0->Write();
956 selBkgPiMass_ncpi0->Write();
957 selBkgPiMass_ncnopi0->Write();
967 std::cout <<
"--- Created root file: \"BDTTest.root\" containing the MVA output histograms" <<
std::endl;
float ncid
Likelihood Neutral Current.
Cuts and Vars for the 2020 FD DiF Study.
SRHeader hdr
Header branch: run, subrun, etc.
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
float nutauid
Likelihood Charge Current NuTau.
float nueid
Likelihood Charge Current NuE.
float pid
PID value output by kNN.
virtual int NFiles() const
May return -1 indicating the number of files is not known.
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
virtual double GetMvaValue(const std::vector< double > &inputValues) const =0
SRCVNResult cvn
Horrible hack to appease CAFAna.
bool IsFiducial(const TVector3 &nuVtx, const TVector3 &min, const TVector3 &max)
SRRemid remid
Output from RecoMuonID (ReMId) package.
float calE
Calorimetric energy of the cluster [GeV].
SRVector3D vtx
Vertex position in detector coordinates. [cm].
short nnu
Number of neutrinos in nu vector (0 or 1)
unsigned int nhit
number of hits
virtual TFile * GetNextFile()=0
Returns the next file in sequence, ready for reading.
The StandardRecord is the primary top-level object in the Common Analysis File trees.
Interface class for accessing ROOT files in sequence.
void SetProgress(double frac)
Update the progress fraction between zero and one.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
IFileSource * WildcardOrSAMQuery(const std::string &str) const
Figure out if str is a wildcard or SAM query and return a source.
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
float numuid
Likelihood Charge Current NuMu.
A simple ascii-art progress bar.
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.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
void Done()
Call this when action is completed.
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
TCut fiducial(x1cut &&y1cut &&z1cut)
SRVertexBranch vtx
Vertex branch: location, time, etc.