testfom.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TLegend.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "TTree.h"
6 #include "TCanvas.h"
7 #include "TGraph.h"
8 #include "TRandom3.h"
9 
10 #include <vector>
11 #include <string>
12 #include <iostream>
13 
14 void testfom(){
15 
16  float NHit, CVNe, SparsenessAsymm, PtP, PxP, PyP, DistTop, DistNotTop, BDT;
17  float vtxx, vtxy, vtxz, prongl;
18  float DistFront, DistEast, DistWest, DistBottom, DistBack;
19  float lid, starty, stopy;
20  float inelast, hitsperpl, numofshow, widthofshow, costheta, remid, cosang, efrac, nplfr;
21 
22  float weight, isGENIE, pdg, pdgorig, isCC, woscdumb;
23 
24  float issig, isbb;
25 
26  TFile *in = TFile::Open("files/test_current_bdt.root");
27 
28  TTree *tree = (TTree*)in->Get("dataset/TestTree");
29 
30  tree->SetBranchAddress("nhit", &NHit);
31  tree->SetBranchAddress("cvnsse", &CVNe);
32  tree->SetBranchAddress("sparseassym", &SparsenessAsymm);
33  tree->SetBranchAddress("ptp", &PtP);
34 // tree->SetBranchAddress("pxp", &PxP);
35 // tree->SetBranchAddress("pyp", &PyP);
36  tree->SetBranchAddress("disttop", &DistTop);
37  tree->SetBranchAddress("distnottop", &DistNotTop);
38 // tree->SetBranchAddress("vtxx", &vtxx);
39 // tree->SetBranchAddress("vtxy", &vtxy);
40 // tree->SetBranchAddress("vtxz", &vtxz);
41 // tree->SetBranchAddress("prongl", &vtxz);
42 // tree->SetBranchAddress("inelast", &inelast);
43 // tree->SetBranchAddress("hitsperpl", &hitsperpl);
44 // tree->SetBranchAddress("numofshow", &numofshow);
45 // tree->SetBranchAddress("widthofshow", &widthofshow);
46 // tree->SetBranchAddress("costheta", &costheta);
47 // tree->SetBranchAddress("remid", &remid);
48 // tree->SetBranchAddress("distfront", &DistFront);
49 // tree->SetBranchAddress("disteast", &DistEast);
50 // tree->SetBranchAddress("distwest", &DistWest);
51 // tree->SetBranchAddress("distbottom", &DistBottom);
52 // tree->SetBranchAddress("distback", &DistBack);
53 // tree->SetBranchAddress("lid", &lid);
54 // tree->SetBranchAddress("stopy", &stopy);
55 // tree->SetBranchAddress("starty", &starty);
56 // tree->SetBranchAddress("cosang", &cosang);
57 // tree->SetBranchAddress("efrac", &efrac);
58 // tree->SetBranchAddress("nplfr", &nplfr);
59 
60  tree->SetBranchAddress("woscdumb", &woscdumb);
61  tree->SetBranchAddress("weight", &weight);
62  tree->SetBranchAddress("isGENIE", &isGENIE);
63  tree->SetBranchAddress("pdg", &pdg);
64  tree->SetBranchAddress("pdgorig", &pdgorig);
65  tree->SetBranchAddress("isCC", &isCC);
66 
67  tree->SetBranchAddress("test_current_bdt", &BDT);
68 
69  int i; double weightbkg =0;
70  TH1F* sig = new TH1F("sig", "BDT", 60, 0.2, 0.8);
71  TH1F* cosm = new TH1F("cosm", "BDT", 60, 0.2, 0.8);
72  TH1F* bkg = new TH1F("bkg", "BDT", 60, 0.2, 0.8);
73  TH1F* tot = new TH1F("tot", "BDT", 60, 0.2, 0.8);
74  for (i = 0; i < tree->GetEntries(); i++){
75  tree->GetEntry(i);
76  if (CVNe > 0.75 ) tot->Fill(BDT,weight*2);
77  if (CVNe > 0.75 && isGENIE && pdg==12 && pdgorig==14 && isCC==1) {
78  sig->Fill(BDT, weight*2);
79  }
80  if (CVNe > 0.75 && isGENIE && (pdg!=12 || pdgorig!=14 || isCC!=1)) {
81  if (isCC) weightbkg=woscdumb*0.000586;
82  if (!isCC) weightbkg=0.5*0.000586; // see docdb 17179 p.3
83  bkg->Fill(BDT, weightbkg*2);}
84  if (CVNe > 0.75 && !isGENIE) {
85  cosm->Fill(BDT, weight*2);
86  }
87 
88  }
89 
90  cout<<"num "<<i<<" "<<sig->Integral()<<" "<<cosm->Integral()<<" "<<bkg->Integral()<<" "<<tot->Integral()<<endl;
91 
92  tot->GetXaxis()->SetTitle("BDT value");
93  tot->GetXaxis()->CenterTitle();
94  tot->GetYaxis()->SetTitle("Events / 9e20");
95  tot->GetYaxis()->CenterTitle();
96  double xlow = 0.6;
97  tot->SetLineColor(kRed);
98  sig->SetLineColor(kBlue);
99  cosm->SetLineColor(kBlack);
100  bkg->SetLineColor(kGreen+2);
101  double max = std::max(tot->GetMaximum(), cosm->GetMaximum());
102  tot->SetMaximum(max*1.2);
103 
104  TCanvas *can = new TCanvas();
105  tot->Add(bkg);
106  tot->Draw("hist");
107  bkg->Draw("same,hist");
108  sig->Draw("same,hist");
109  cosm->Draw("same,hist");
110 
111  TLegend *leg = new TLegend(xlow,0.54,xlow+0.2,0.85);
112  leg->AddEntry(tot,"Total","l");
113  leg->AddEntry(sig,"Signal","l");
114  leg->AddEntry(cosm,"Cosmic bkg","l");
115  leg->AddEntry(bkg,"Beam bkg","l");
116  leg->Draw();
117  can->SaveAs("figures/Fig_current_bdt.pdf");
118 
119  double bestfom = 0; double bestcvn = 0; double bestbdt = 0;
120  double curcvn, curbdt, curfom;
121  double bestsig = 0; double bestbkg = 0; double bestcos = 0;
122  double localbestfom, localbestbdt, localbestsig, localbestbkg, localbestcos;
123 
124  int k = 0;
125  for (k=20; k<50; k++) {
126  curcvn=0.5+k*0.01;
127  //cout<<curcvn;
128 
129  sig->Reset();
130  bkg->Reset();
131  cosm->Reset();
132  tot->Reset();
133  for (i = 0; i < tree->GetEntries(); i++){
134  tree->GetEntry(i);
135  if (CVNe < curcvn) continue;
136  tot->Fill(BDT,weight*2); // *2 since set was divided into test and training trees which are ~equal in number of events
137  if (isGENIE && pdg==12 && pdgorig==14 && isCC==1) {
138  sig->Fill(BDT, weight*2);
139  }
140  if (isGENIE && (pdg!=12 || pdgorig!=14 || isCC!=1)) {
141  if (isCC) weightbkg=woscdumb*0.000586;
142  if (!isCC) weightbkg=0.5*0.000586; // see docdb 17179 p.3
143  bkg->Fill(BDT, weightbkg*2);}
144  if (!isGENIE) {
145  cosm->Fill(BDT, weight*2);
146  }
147  }
148 
149 
150  for (int j = 1; j <= 60; j++){
151  double sign = sig->Integral(j,60);
152  double bkgd= bkg->Integral(j,60) + cosm->Integral(j,60);
153  double curfom = sign/sqrt(sign+bkgd);
154  curbdt = j;
155  double cursig=sig->Integral(j, 60);
156  double curbkg = bkg->Integral(j, 60);
157  double curcos = cosm->Integral(j, 60);
158  if (curfom > bestfom){
159  localbestfom = curfom;
160  localbestbdt = j;
161  localbestsig = cursig;
162  localbestbkg = curbkg;
163  localbestcos = curcos;
164  }
165  }
166 
167 
168 cout<<"----local : bestfom "<<localbestfom<<" bestcvn "<<curcvn<<" bestbdt "<<localbestbdt<<" signal "<<localbestsig<<" bkg "<<localbestbkg<<" cosmic "<<localbestcos<<endl;
169 
170  if (localbestfom > bestfom){
171  bestfom = localbestfom;
172  bestbdt = localbestbdt;
173  bestcvn = curcvn;
174  bestsig = localbestsig;
175  bestbkg = localbestbkg;
176  bestcos = localbestcos;
177  }
178 
179 
180 }
181 cout<<"bestfom "<<bestfom<<" bestcvn "<<bestcvn<<" bestbdt "<<bestbdt<<endl;
182 cout<<" signal "<<bestsig<<" bkg "<<bestbkg<<" cosmic "<<bestcos<<endl;
183 
184 }
T max(const caf::Proxy< T > &a, T b)
enum BeamMode kRed
const Var weight
T sqrt(T number)
Definition: d0nt_math.hpp:156
float vtxx
float vtxy
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
float vtxz
ifstream in
Definition: comparison.C:7
void testfom()
Definition: testfom.C:14
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kGreen
enum BeamMode kBlue
def sign(x)
Definition: canMan.py:197