FairRoot/PandaRoot
cutfinder_toyev.C
Go to the documentation of this file.
1 #include <algorithm>
2 #include "TFile.h"
3 #include "TTree.h"
4 #include "TString.h"
5 #include "TH1F.h"
6 #include "TCanvas.h"
7 #include "TROOT.h"
8 #include "TEventList.h"
9 #include "TDirectory.h"
10 #include <iostream>
11 #include "TLine.h"
12 #include "TLatex.h"
13 #include "TStyle.h"
14 #include "TObjArray.h"
15 #include "TPRegexp.h"
16 #include <map>
17 #include <utility>
18 #include <algorithm>
19 
20 #define MAX 1000
21 #define BINS 500
22 #define NCAN 3
23 
24 std::map<int, int> evcnt;
25 
26 double signi[MAX];
27 int idx[MAX];
28 TString vars[MAX];
29 double cut[MAX];
30 
31 double minh[MAX], maxh[MAX];
32 
33 double N0_sig, N0_bg;
34 double Nsigev, Nbgev;
35 
36 bool mycompare(int i, int j)
37 {
38  return signi[i]>signi[j];
39 }
40 
41 int countEvents(TTree *t, TString ccut)
42 {
43  t->SetEventList(0);
44  t->SetBranchStatus("*",1);
45  t->Draw(">>el",ccut);
46  t->SetBranchStatus("*",0);
47  t->SetBranchStatus("evt",1);
48 
49  TEventList *el=(TEventList*)gDirectory->Get("el");
50 
51  Float_t ev;
52  t->SetBranchAddress("evt",&ev);
53 
54  evcnt.clear();
55 
56  for (int i=0;i<el->GetN();++i)
57  {
58  t->GetEntry(el->GetEntry(i));
59  evcnt[ev]+=1;
60  }
61  t->SetBranchStatus("*",1);
62 
63  return evcnt.size();
64 }
65 int countEvents(TTree *t, TEventList &el)
66 {
67  t->SetEventList(&el);
68  t->SetBranchStatus("*",0);
69  t->SetBranchStatus("evt",1);
70 
71  Float_t ev;
72  t->SetBranchAddress("evt",&ev);
73 
74  evcnt.clear();
75 
76  for (int i=0;i<el.GetN();++i)
77  {
78  t->GetEntry(el.GetEntry(i));
79  evcnt[ev]+=1;
80  }
81  t->SetBranchStatus("*",1);
82 
83  return evcnt.size();
84 }
85 
86 double bestSepEvt(TTree *t, TString varname, TEventList &els, TEventList &elb, double &bestcut, double eff, int idx)
87 {
88  int i,j;
89  t->SetBranchStatus("*",0);
90  t->SetBranchStatus("evt",1);
91  t->SetBranchStatus(varname,1);
92 
93  Float_t ev, var;
94  t->SetBranchAddress("evt",&ev);
95  t->SetBranchAddress(varname,&var);
96 
97  std::map<int, int> sigcnt;
98  std::map<int, int> bgcnt;
99  std::map<int, int> bgcnt2;
100 
101  std::vector<pair<double, int> > sigvals;
102  std::vector<pair<double, int> > bgvals;
103 
104  int Nsigval = els.GetN(), Nbgval = elb.GetN();
105  double min = 1e9, max=-1e9;
106 
107  // prepare the pairs of variable value, eventnumber for signal and count signals
108  for (i=0;i<els.GetN();++i)
109  {
110  t->GetEntry(els.GetEntry(i));
111  sigvals.push_back(std::make_pair(var, ev));
112  if (var>max) max = var;
113  if (var<min) min = var;
114  }
115 
116  // prepare the pairs of variable value, eventnumber for background and count background
117  for (i=0;i<elb.GetN();++i)
118  {
119  t->GetEntry(elb.GetEntry(i));
120  bgvals.push_back(std::make_pair(var, ev));
121  if (var>max) max = var;
122  if (var<min) min = var;
123  }
124 
125  minh[idx]=min; maxh[idx]=max;
126 
127  // sort signals by first element = variable value
128  sort(sigvals.begin(), sigvals.end());
129  sort(bgvals.begin(), bgvals.end());
130 
131  double pos=min, step = (max-min)/500., bestqal=0, bestcutl=0, bestqar=0, bestcutr=0;
132  int bestsig=0, bestbg=0;
133 
134  i=0, j=0;
135  sigcnt.clear();
136  bgcnt.clear();
137 
138  for (pos=min+step; pos<max; pos+=step)
139  {
140  while ( sigvals[i].first<=pos && i<Nsigval) sigcnt[sigvals[i++].second]+=1;
141  while ( bgvals[j].first<=pos && j<Nbgval) bgcnt[bgvals[j++].second]+=1;
142 
143  double qa = sigcnt.size()/Nsigev * (Nbgev-bgcnt.size())/Nbgev;
144  if (qa>bestqal)
145  {
146  bestqal = qa;
147  bestcutl = pos;
148  bestsig = sigcnt.size();
149  bestbg = bgcnt.size();
150  }
151  }
152 
153  cout <<varname<<" qa(l)="<<bestqal<<" S="<<bestsig <<" B="<<bestbg;
154 
155  i=sigvals.size()-1, j=bgvals.size()-1;
156  sigcnt.clear();
157  bgcnt.clear();
158 
159  for (pos=max-step; pos>min; pos-=step)
160  {
161  while ( sigvals[i].first>=pos && i>=0) sigcnt[sigvals[i--].second]+=1;
162  while ( bgvals[j].first>=pos && j>=0) bgcnt[bgvals[j--].second]+=1;
163 
164  double qa = sigcnt.size()/Nsigev * (Nbgev-bgcnt.size())/Nbgev;
165  if (qa>bestqar)
166  {
167  bestqar = qa;
168  bestcutr = pos;
169  bestsig = sigcnt.size();
170  bestbg = bgcnt.size();
171  }
172  }
173  cout <<" qa(r)="<<bestqar<<" S="<<bestsig <<" B="<<bestbg<<endl;
174 
175 
176  t->SetBranchStatus("*",1);
177 
178  bestcut=bestcutr;
179  if (bestqal>bestqar)
180  {
181  bestcut = bestcut;
182  return bestqal;
183  }
184 
185  return bestqar;
186 }
187 
188 
189 double bestSuppressionEvt(TTree *t, TString varname, TEventList &els, TEventList &elb, double &bestcut, double eff, int idx)
190 {
191  int i;
192  t->SetBranchStatus("*",0);
193  t->SetBranchStatus("evt",1);
194  t->SetBranchStatus(varname,1);
195 
196  Float_t ev, var;
197  t->SetBranchAddress("evt",&ev);
198  t->SetBranchAddress(varname,&var);
199 
200  std::map<int, int> sigcnt;
201  std::map<int, int> bgcnt;
202  std::map<int, int> bgcnt2;
203 
204  std::vector<pair<double, int> > sigvals;
205  std::vector<pair<double, int> > bgvals;
206 
207  int Nsigval = els.GetN(), Nbgval = elb.GetN();
208  double min = 1e9, max=-1e9;
209 
210  // prepare the pairs of variable value, eventnumber for signal and count signals
211  for (i=0;i<els.GetN();++i)
212  {
213  t->GetEntry(els.GetEntry(i));
214  sigvals.push_back(std::make_pair(var, ev));
215  if (var>max) max = var;
216  if (var<min) min = var;
217  }
218 
219  // prepare the pairs of variable value, eventnumber for background and count background
220  for (i=0;i<elb.GetN();++i)
221  {
222  t->GetEntry(elb.GetEntry(i));
223  bgvals.push_back(std::make_pair(var, ev));
224  if (var>max) max = var;
225  if (var<min) min = var;
226  }
227 
228  // sort signals by first element = variable value
229  sort(sigvals.begin(), sigvals.end());
230 
231  i=0;
232  sigcnt.clear();
233  // find cut for current efficiency requirement left cut
234  while ( (sigcnt.size()/(double)Nsigev)<eff && i<Nsigval ) sigcnt[sigvals[i++].second]+=1;
235  double leftcut = sigvals[i].first;
236 
237  sigcnt.clear();
238  i=sigvals.size()-1;
239  // find cut for current efficiency requirement right cut
240  while ( (sigcnt.size()/(double)Nsigev)<eff && i>=0) sigcnt[sigvals[i--].second]+=1;
241  double rightcut = sigvals[i].first;
242 
243  bgcnt.clear();
244  bgcnt2.clear();
245 
246  for (i=0;i<bgvals.size();++i)
247  {
248  if (bgvals[i].first<=leftcut) bgcnt[bgvals[i].second];
249  if (bgvals[i].first>=rightcut) bgcnt2[bgvals[i].second];
250  }
251 
252  minh[idx]=min; maxh[idx]=max;
253 
254  t->SetBranchStatus("*",1);
255 
256 
257  double leftsupr = (Nbgev-bgcnt.size())/Nbgev;
258  double rightsupr = (Nbgev-bgcnt2.size())/Nbgev;
259 
260  cout <<varname<<"(S="<<sigcnt.size()<<") l:"<<leftcut<<"("<<bgcnt.size()<<"/"<<leftsupr<<") r:"<<rightcut<<"("<<bgcnt2.size()<<"/"<<rightsupr<<")"<<endl;
261 
262  bestcut = rightcut;
263  if (leftsupr>rightsupr)
264  {
265  bestcut = leftcut;
266  return leftsupr;
267  }
268  return rightsupr;
269 }
270 
271 double bestEffEvt(TTree *t, TString varname, TEventList &els, TEventList &elb, double &bestcut, double supr, int idx)
272 {
273  int i;
274  t->SetBranchStatus("*",0);
275  t->SetBranchStatus("evt",1);
276  t->SetBranchStatus(varname,1);
277 
278  Float_t ev, var;
279  t->SetBranchAddress("evt",&ev);
280  t->SetBranchAddress(varname,&var);
281 
282  std::map<int, int> sigcnt;
283  std::map<int, int> bgcnt;
284  std::map<int, int> sigcnt2;
285 
286  std::vector<pair<double, int> > sigvals;
287  std::vector<pair<double, int> > bgvals;
288 
289  int Nsigval = els.GetN(), Nbgval = elb.GetN();
290  double min = 1e9, max=-1e9;
291 
292  // prepare the pairs of variable value, eventnumber for signal and count signals
293  for (i=0;i<els.GetN();++i)
294  {
295  t->GetEntry(els.GetEntry(i));
296  sigvals.push_back(std::make_pair(var, ev));
297  if (var>max) max = var;
298  if (var<min) min = var;
299  }
300 
301  // prepare the pairs of variable value, eventnumber for background and count background
302  for (i=0;i<elb.GetN();++i)
303  {
304  t->GetEntry(elb.GetEntry(i));
305  bgvals.push_back(std::make_pair(var, ev));
306  if (var>max) max = var;
307  if (var<min) min = var;
308  }
309 
310  // sort signals by first element = variable value
311  sort(bgvals.begin(), bgvals.end());
312 
313  i=0;
314  bgcnt.clear();
315  // find cut for current efficiency requirement left cut
316  while ( (bgcnt.size()/(double)Nbgev)<(1-supr) && i<Nbgval ) bgcnt[bgvals[i++].second]+=1;
317  double leftcut = bgvals[i].first;
318 
319  bgcnt.clear();
320  i=bgvals.size()-1;
321  // find cut for current efficiency requirement right cut
322  while ( (bgcnt.size()/(double)Nbgev)<(1-supr) && i>=0) bgcnt[bgvals[i--].second]+=1;
323  double rightcut = bgvals[i].first;
324 
325  sigcnt.clear();
326  sigcnt2.clear();
327 
328  for (i=0;i<sigvals.size();++i)
329  {
330  if (sigvals[i].first<leftcut) sigcnt[sigvals[i].second]+=1;
331  if (sigvals[i].first>rightcut) sigcnt2[sigvals[i].second]+=1;
332  }
333 
334  minh[idx]=min; maxh[idx]=max;
335 
336  t->SetBranchStatus("*",1);
337 
338  double lefteff = sigcnt.size()/Nsigev;
339  double righteff = sigcnt2.size()/Nsigev;
340 
341  cout <<varname<<"(B="<<bgcnt.size()<<") l:"<<leftcut<<"("<<sigcnt.size()<<"/"<<lefteff<<") r:"<<rightcut<<"("<<sigcnt2.size()<<"/"<<righteff<<")"<<endl;
342 
343  bestcut = rightcut;
344  if (lefteff>righteff)
345  {
346  bestcut = leftcut;
347  return lefteff;
348  }
349  return righteff;
350 }
351 
352 void findLimits(TTree *t, TString var, TString ccut, double &low, double &high, double frac = 0.98)
353 {
354  low = t->GetMinimum(var);
355  high = t->GetMaximum(var);
356  return;
357 
358  double miss = (1.-frac)/2;
359 
360  t->SetEventList(0);
361  TEventList el("el");
362  t->Draw(">>el",ccut);
363  t->SetEventList(&el);
364 
365  double llow = t->GetMinimum(var);
366  double lhigh = t->GetMaximum(var);
367 
368  TH1F htemp("htemp","",500,llow, lhigh);
369  t->Project("htemp",var);
370 
371  double sum=0.;
372  double integ = htemp.Integral();
373 
374  int i=1;
375 
376  while (sum<miss) sum+=htemp.GetBinContent(i++)/integ;
377  low = htemp.GetBinCenter(i-2);
378  sum=0.; i=500;
379 
380  while (sum<miss) sum+=htemp.GetBinContent(i--)/integ;
381  high = htemp.GetBinCenter(i+2);
382 
383  t->SetEventList(0);
384 }
385 
386 int cutfinder_toyev(TString pre, TString ntp, TString sigcut, TString precut="", double n0s=50., double n0b=500., double supr=0.9 )
387 {
388  gStyle->SetTitleY(0.993);
389  gStyle->SetTitleH(0.07);
390  gStyle->SetPadTopMargin(0.08);
391  gStyle->SetTitleBorderSize(0);
392  gStyle->SetOptStat(0);
393 
394  N0_sig = n0s*1000.;
395  N0_bg = n0b*1000.;
396 
397 
398  int i,j;
399  double sigl, sigh, bgl, bgh;
400 
401  TFile *f=new TFile("data/M"+pre+"_"+ntp+".root","READ");
402  //TFile *f=new TFile(fname,"READ");
403  TTree *t=(TTree*)f->Get(ntp);
404 
405  TEventList els("els");
406  TEventList elsall("elsall");
407  TEventList elb("elb");
408 
409  TCanvas *c1=new TCanvas("c1","c1",10,10,1500,900);
410  c1->Divide(5,3);
411 
412  TObjArray* branches = t->GetListOfBranches();
413 
414  TString bgcut = "!("+sigcut+")";
415  //sigcut+="&&mct";
416 
417  if (precut!="")
418  {
419  sigcut+="&&"+precut;
420  bgcut +="&&"+precut;
421  }
422  cout <<sigcut <<" "<<bgcut<<endl;
423  t->Draw(">>elsall",sigcut);
424  t->Draw(">>els",sigcut+"&&mct");
425  t->Draw(">>elb",bgcut);
426 
427  Nsigev = countEvents(t, els);
428  Nbgev = countEvents(t, elb);
429 
430  cout <<"SIG:"<<Nsigev<<" ev "<<els.GetN()<<" cn BG:"<<Nbgev<<" ev "<<elb.GetN()<<" cn"<<endl;
431 
432  int cnt=1;
433 
434  for (i=0;i<MAX;++i) { idx[i]=i; signi[i]=0.;}
435 
436  for(i=0; i<=branches->GetLast(); ++i)
437  {
438 
439  TBranch* branch = (TBranch*)branches->UncheckedAt(i);
440  vars[i]=branch->GetName();
441  //t->SetBranchStatus(vars[i],1);
442 
443  if (vars[i]=="evt") continue;
444  if (vars[i]=="mode") continue;
445  if (vars[i].Contains("pdg")) continue;
446  if (vars[i].Contains("mct")) continue;
447  if (vars[i].Contains("dec")) continue;
448  if (vars[i].BeginsWith("t") && vars[i]!="thr") continue;
449 
450  if (supr<0)
451  signi[i] = bestEffEvt(t, vars[i], els, elb, cut[i], -supr, i);
452  else
453  signi[i] = bestSuppressionEvt(t, vars[i], els, elb, cut[i], supr, i);
454 // signi[i] = bestEffEvt(t, vars[i], els, elb, cut[i], supr, i);
455 
456  }
457 
458  //t->SetBranchStatus("*",1);
459 
460  cout <<"\n\nBEST 12 vars:"<<endl<<endl;
461 
462  std::vector<int> myidx (idx, idx+MAX);
463 
464  std::sort(myidx.begin(), myidx.end(), mycompare);
465 
466  TLine l;
467  l.SetLineColor(6);
468  l.SetLineStyle(2);
469  l.SetLineWidth(2);
470  TLatex lt;
471  lt.SetTextSize(0.06);
472  TString target="supr";
473  if (supr<0) target="eff";
474 
475  for (j=0;j<15;++j)
476  {
477  i=myidx[j];
478  cout << vars[i]<<" : qa = "<<signi[i]<<" cut = "<<cut[i]<<endl;
479  if (j<15)
480  {
481 /* findLimits(t,vars[i],sigcut, sigl, sigh);
482  findLimits(t,vars[i],bgcut, bgl, bgh);
483 
484  if (sigl>bgl) sigl=bgl;
485  if (sigh<bgh) sigh=bgh;
486  */
487  TH1F h1("h1",vars[i],200,minh[i],maxh[i]);
488  TH1F h2("h2",vars[i],200,minh[i],maxh[i]);
489  h2.SetLineColor(2);
490 
491  t->SetEventList(&els);
492  t->Project("h1",vars[i]);
493  t->SetEventList(&elb);
494  t->Project("h2",vars[i]);
495 
496  h1.Scale(1.0/h1.Integral());
497  h2.Scale(1.0/h2.Integral());
498  h1.SetTitleSize(0.05);
499 
500  c1->cd(j+1);
501  double maxi = h1.GetMaximum();
502  if (h2.GetMaximum()>maxi ) maxi = h2.GetMaximum();
503  maxi*=1.1;
504  h1.SetMaximum(maxi);
505  h2.SetMaximum(maxi);
506 
507  h1.DrawNormalized();
508  h2.DrawNormalized("same");
509 
510  double axmin = h1.GetXaxis()->GetXmin(),axmax = h1.GetXaxis()->GetXmax();
511 
512  l.DrawLine(cut[i],0, cut[i], maxi*0.9);
513 
514  lt.DrawLatex(axmin+(axmax-axmin)*0.6,1.01*maxi,TString::Format("%s = %6.4f",target.Data(),signi[i]));
515 
516  //c1->Update();
517  }
518 
519  }
520 
521  c1->Update();
522 
523  float nsig = countEvents(t,elsall);
524  float nbg = countEvents(t,elb);
525 
526  cout <<"SIG EVT : "<<nsig<<" "<<"BG EVT : "<<nbg<<endl;
527  cout <<"SIG EFF : "<<nsig/N0_sig<<" "<<"BG EFF : "<< nbg/N0_bg<<endl;
528  return 0;
529 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
cut
double cut[MAX]
Definition: cutfinder_toyev.C:29
els
TEventList * els
Definition: findcuts.C:52
Nsigev
double Nsigev
Definition: cutfinder_toyev.C:34
cutfinder_toyev
int cutfinder_toyev(TString pre, TString ntp, TString sigcut, TString precut="", double n0s=50., double n0b=500., double supr=0.9)
Definition: cutfinder_toyev.C:386
i
Int_t i
Definition: run_full.C:25
low
int low
Definition: anaMvdDigi.C:51
bestSepEvt
double bestSepEvt(TTree *t, TString varname, TEventList &els, TEventList &elb, double &bestcut, double eff, int idx)
Definition: cutfinder_toyev.C:86
ev
int ev
Definition: GammaSpectraAnalysis.C:44
vars
TString vars[MAX]
Definition: cutfinder_toyev.C:28
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
evcnt
std::map< int, int > evcnt
Definition: cutfinder_toyev.C:24
N0_sig
double N0_sig
Definition: cutfinder_toyev.C:33
h1
TH1F * h1
Definition: anal_hit_digi_cluster_fwendcap.C:121
findLimits
void findLimits(TTree *t, TString var, TString ccut, double &low, double &high, double frac=0.98)
Definition: cutfinder_toyev.C:352
bestSuppressionEvt
double bestSuppressionEvt(TTree *t, TString varname, TEventList &els, TEventList &elb, double &bestcut, double eff, int idx)
Definition: cutfinder_toyev.C:189
MAX
#define MAX
Definition: cutfinder_toyev.C:20
Nbgev
double Nbgev
Definition: cutfinder_toyev.C:34
minh
double minh[MAX]
Definition: cutfinder_toyev.C:31
bestEffEvt
double bestEffEvt(TTree *t, TString varname, TEventList &els, TEventList &elb, double &bestcut, double supr, int idx)
Definition: cutfinder_toyev.C:271
idx
int idx[MAX]
Definition: cutfinder_toyev.C:27
countEvents
int countEvents(TTree *t, TString ccut)
Definition: cutfinder_toyev.C:41
f
TFile * f
Definition: bump_analys.C:12
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
maxh
double maxh[MAX]
Definition: cutfinder_toyev.C:31
mycompare
bool mycompare(int i, int j)
Definition: cutfinder_toyev.C:36
signi
double signi[MAX]
Definition: cutfinder_toyev.C:26
c1
c1
Definition: plot_dirc.C:35
nsig
int nsig
Definition: toy_core.C:46
N0_bg
double N0_bg
Definition: cutfinder_toyev.C:33
h2
TH1F * h2
Definition: anal_hit_digi_cluster_fwendcap.C:122
elb
TEventList * elb
Definition: findcuts.C:52
cnt
Int_t cnt
Definition: hist-t7.C:106
high
double high
Definition: full_core_ntp.C:35
t
TTree * t
Definition: bump_analys.C:13
htemp
TH1D * htemp[nsteps]
Definition: dedx_bands.C:61
Generated on Sun Aug 18 2019 03:01:09 for FairRoot/PandaRoot by   doxygen 1.8.8