FairRoot/PandaRoot
cutfinder.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 
18 #define MAX 1000
19 #define BINS 500
20 #define NCAN 3
21 
22 std::map<int, int> evcnt;
23 
24 double signi[MAX];
25 int idx[MAX];
26 TString vars[MAX];
27 double cut[MAX];
28 
29 double N0_sig, N0_bg;
30 
31 bool mycompare(int i, int j)
32 {
33  return signi[i]>signi[j];
34 }
35 
36 double bestSignificance(TH1* h1, TH1* h2, double &bestcut, double fac2=1.)
37 {
38  int n=h1->GetNbinsX();
39  double sum1=0, sum2=0;
40  double bestr=0, bestl=0;
41  double int1 = h1->Integral();
42  double int2 = h2->Integral();
43  double bestcutl=-999.;
44  double bestcutr=999.;
45 
46  bestcut=-999.;
47 
48  for (int i=1;i<n;++i)
49  {
50  double cur1 = h1->GetBinContent(i)/int1;
51  double cur2 = h2->GetBinContent(i)/int2;
52  sum1 += cur1;
53  sum2 += cur2;
54 
55  double vall = 0, valr=0;
56 
57  if ((sum1+fac2*sum2)>0) vall=sum1/sqrt(fac2*sum2+sum1);
58  if ((1.-sum1+fac2*(1.-sum2))>0) valr=(1.-sum1)/sqrt(1.-sum1+fac2*(1.-sum2));
59 
60  if (valr>bestr) { bestr=valr; bestcutr = h1->GetBinCenter(i);}
61  if (vall>bestl) { bestl=vall; bestcutl = h1->GetBinCenter(i);}
62  }
63 
64  if (bestr>bestl)
65  {
66  bestcut = bestcutr;
67  return bestr;
68  }
69 
70  bestcut = bestcutl;
71  return bestl;
72 }
73 
74 int countEvents(TTree *t, TString ccut)
75 {
76  t->SetEventList(0);
77  t->SetBranchStatus("*",1);
78  t->Draw(">>el",ccut);
79  t->SetBranchStatus("*",0);
80  t->SetBranchStatus("evt",1);
81 
82  TEventList *el=(TEventList*)gDirectory->Get("el");
83 
84  Float_t ev;
85  t->SetBranchAddress("evt",&ev);
86 
87  evcnt.clear();
88 
89  for (int i=0;i<el->GetN();++i)
90  {
91  t->GetEntry(el->GetEntry(i));
92  evcnt[ev]+=1;
93  }
94  t->SetBranchStatus("*",1);
95 
96  return evcnt.size();
97 }
98 
99 double diffQA(TH1F &ht1, TH1F &ht2)
100 {
101  int n=ht1.GetNbinsX();
102  TH1F h1=ht1;
103  TH1F h2=ht2;
104  h1.Scale(1.0/h1.Integral());
105  h2.Scale(1.0/h2.Integral());
106  h1.Add(&h2,-1);
107 
108  double sum=0.;
109 
110  for (int i=1;i<=n;++i)
111  {
112  sum+=fabs(h1.GetBinContent(i));
113  }
114 
115  return sum;
116 }
117 
118 double bestEff(TH1* h1, TH1* h2, double &bestcut, double supr2=0.8, double fac2=1.)
119 {
120  int n=h1->GetNbinsX();
121  double sum1=0, sum2=0;
122  double bestr=0, bestl=0;
123  double int1 = h1->Integral();
124  double int2 = h2->Integral();
125  double bestcutl=-999.;
126  double bestcutr=999.;
127 
128  bestcut=-999.;
129  int i=1;
130 
131  while (sum2<supr2)
132  {
133  sum1 += h1->GetBinContent(i)/int1;
134  sum2 += h2->GetBinContent(i++)/int2;
135  }
136  bestr=(1.0-sum1);
137  bestcutr = h1->GetBinCenter(i-1);
138 
139  sum1=0;
140  sum2=0;
141  i=n;
142 
143  while (sum2<supr2)
144  {
145  sum1 += h1->GetBinContent(i)/int1;
146  sum2 += h2->GetBinContent(i--)/int2;
147  }
148 
149  bestl=(1.0-sum1);
150  bestcutl = h1->GetBinCenter(i-1);
151 
152  if (bestr>bestl)
153  {
154  bestcut = bestcutr;
155  return bestr;
156  }
157 
158  bestcut = bestcutl;
159  return bestl;
160 }
161 double bestSuppression(TH1* h1, TH1* h2, double &bestcut, double eff=0.9)
162 {
163  int n=h1->GetNbinsX();
164  double sum1=0, sum2=0;
165  double bestr=0, bestl=0;
166  double int1 = h1->Integral();
167  double int2 = h2->Integral();
168  double bestcutl=-999.;
169  double bestcutr=999.;
170 
171  bestcut=-999.;
172  int i=1;
173 
174  while (sum1<eff)
175  {
176  sum1 += h1->GetBinContent(i)/int1;
177  sum2 += h2->GetBinContent(i++)/int2;
178  }
179  bestr=(1.0-sum2);
180  bestcutr = h1->GetBinCenter(i-1);
181 
182  sum1=0;
183  sum2=0;
184  i=n;
185 
186  while (sum1<eff)
187  {
188  sum1 += h1->GetBinContent(i)/int1;
189  sum2 += h2->GetBinContent(i--)/int2;
190  }
191 
192  bestl=(1.0-sum2);
193  bestcutl = h1->GetBinCenter(i-1);
194 
195  if (bestr>bestl)
196  {
197  bestcut = bestcutr;
198  return bestr;
199  }
200 
201  bestcut = bestcutl;
202  return bestl;
203 }
204 
205 void findLimits(TTree *t, TString var, TString ccut, double &low, double &high, double frac = 0.98)
206 {
207  low = t->GetMinimum(var);
208  high = t->GetMaximum(var);
209  return;
210 
211  double miss = (1.-frac)/2;
212 
213  t->SetEventList(0);
214  TEventList el("el");
215  t->Draw(">>el",ccut);
216  t->SetEventList(&el);
217 
218  double llow = t->GetMinimum(var);
219  double lhigh = t->GetMaximum(var);
220 
221  TH1F htemp("htemp","",500,llow, lhigh);
222  t->Project("htemp",var);
223 
224  double sum=0.;
225  double integ = htemp.Integral();
226 
227  int i=1;
228 
229  while (sum<miss) sum+=htemp.GetBinContent(i++)/integ;
230  low = htemp.GetBinCenter(i-2);
231  sum=0.; i=500;
232 
233  while (sum<miss) sum+=htemp.GetBinContent(i--)/integ;
234  high = htemp.GetBinCenter(i+2);
235 
236  t->SetEventList(0);
237 }
238 
239 int cutfinder(TString pre, TString ntp, TString sigcut, TString precut="", double n0s=50., double n0b=200., double supr=0.9 )
240 {
241 
242  gStyle->SetTitleY(0.993);
243  gStyle->SetTitleH(0.07);
244  gStyle->SetPadTopMargin(0.08);
245  gStyle->SetTitleBorderSize(0);
246  gStyle->SetOptStat(0);
247 
248  N0_sig = n0s*1000.;
249  N0_bg = n0b*1000.;
250 
251  int i,j;
252  double sigl, sigh, bgl, bgh;
253 
254  TFile *f=new TFile("M"+pre+"_"+ntp+".root","READ");
255  //TFile *f=new TFile(fname,"READ");
256  TTree *t=(TTree*)f->Get(ntp);
257 
258  TCanvas *c1=new TCanvas("c1","c1",10,10,1800,600);
259  c1->Divide(6,2);
260 
261  TObjArray* branches = t->GetListOfBranches();
262 
263  TString bgcut = "!("+sigcut+")";
264  TString bstring = precut+" "+sigcut;
265 
266 /* if (precut!="")
267  {
268  sigcut+="&&"+precut;
269  bgcut +="&&"+precut;
270  }*/
271 
272  t->SetEventList(0);
273  TEventList el1("el1");
274  t->Draw(">>el1",precut);
275  t->SetEventList(&el1);
276 
277 
278  int cnt=1;
279 
280  for (i=0;i<MAX;++i) { idx[i]=i; signi[i]=0.;}
281 
282 /* bstring.ReplaceAll("("," ");bstring.ReplaceAll(")"," ");
283  bstring.ReplaceAll("&"," ");bstring.ReplaceAll("|"," ");bstring.ReplaceAll("!"," ");
284  bstring.ReplaceAll(">"," ");bstring.ReplaceAll("<"," ");bstring.ReplaceAll("="," ");
285 
286  bstring+=" ev";
287 
288  TObjArray* ar=bstring.Tokenize(" ");
289 
290  t->SetBranchStatus("*",0);
291 
292  for (Int_t i = 0; i <ar->GetLast()+1; i++)
293  {
294  TString subStr = ((TObjString *)ar->At(i))->GetString();
295  if (!subStr.IsFloat())
296  {
297  cout << "\"" << subStr << "\" ";
298  t->SetBranchStatus(subStr,1);
299  }
300  }
301  cout <<endl;*/
302  for(i=0; i<=branches->GetLast(); ++i)
303  {
304 
305  TBranch* branch = (TBranch*)branches->UncheckedAt(i);
306  vars[i]=branch->GetName();
307  //t->SetBranchStatus(vars[i],1);
308 
309  if (vars[i]=="evt") continue;
310  if (vars[i]=="mode") continue;
311  if (vars[i].Contains("dec")) continue;
312  if (vars[i].BeginsWith("t") && vars[i]!="thr") continue;
313 
314  findLimits(t,vars[i],sigcut, sigl, sigh);
315  findLimits(t,vars[i],bgcut, bgl, bgh);
316 
317  if (sigl>bgl) sigl=bgl;
318  if (sigh<bgh) sigh=bgh;
319 
320 /* cout <<sigl<<" - "<<sigh<<endl;*/
321  if (sigl==sigh) continue;
322  TH1F h1("h1","",200,sigl,sigh);
323  TH1F h2("h2","",200,sigl,sigh);
324  h2.SetLineColor(2);
325 
326  t->Project("h1",vars[i],sigcut);
327  t->Project("h2",vars[i], bgcut);
328 
329  if (h1.GetEntries()<3 || h2.GetEntries()<3) continue;
330 
331  signi[i] = bestSuppression(&h1,&h2, cut[i], supr);
332  //signi[i] = bestEff(&h1,&h2, cut[i], supr);
333 
334  cout <<vars[i]<<" "<<flush;//<<" : "<< vars[i]<<" : signi="<<signi[i]<<endl;
335 
336 /* c1->cd(cnt%20+1);
337  h2.DrawNormalized();
338  h1.DrawNormalized("same");
339  c1->Update();
340  cnt++;*/
341  //t->SetBranchStatus(vars[i],0);
342 
343  }
344 
345  //t->SetBranchStatus("*",1);
346 
347  cout <<"\n\nBEST 12 vars:"<<endl<<endl;
348 
349  std::vector<int> myidx (idx, idx+MAX);
350 
351  std::sort(myidx.begin(), myidx.end(), mycompare);
352 
353  TLine l;
354  l.SetLineColor(6);
355  l.SetLineStyle(2);
356  l.SetLineWidth(2);
357  TLatex lt;
358  lt.SetTextSize(0.06);
359 
360  for (j=0;j<12;++j)
361  {
362  i=myidx[j];
363  cout << vars[i]<<" : qa = "<<signi[i]<<" cut = "<<cut[i]<<endl;
364  if (j<12)
365  {
366  findLimits(t,vars[i],sigcut, sigl, sigh);
367  findLimits(t,vars[i],bgcut, bgl, bgh);
368 
369  if (sigl>bgl) sigl=bgl;
370  if (sigh<bgh) sigh=bgh;
371 
372  TH1F h1("h1",vars[i],200,sigl,sigh);
373  TH1F h2("h2",vars[i],200,sigl,sigh);
374  h2.SetLineColor(2);
375 
376  t->Project("h1",vars[i], sigcut);
377  t->Project("h2",vars[i], bgcut);
378 
379  h1.Scale(1.0/h1.Integral());
380  h2.Scale(1.0/h2.Integral());
381  h1.SetTitleSize(0.05);
382 
383  c1->cd(j+1);
384  double maxi = h1.GetMaximum();
385  if (h2.GetMaximum()>maxi ) maxi = h2.GetMaximum();
386  maxi*=1.1;
387  h1.SetMaximum(maxi);
388  h2.SetMaximum(maxi);
389 
390  h1.DrawNormalized();
391  h2.DrawNormalized("same");
392 
393  double axmin = h1.GetXaxis()->GetXmin(),axmax = h1.GetXaxis()->GetXmax();
394 
395  l.DrawLine(cut[i],0, cut[i], maxi*0.9);
396 
397  lt.DrawLatex(axmin+(axmax-axmin)*0.6,1.01*maxi,TString::Format("sig = %6.4f",signi[i]));
398 
399  c1->Update();
400  }
401 
402  }
403 
404  if (precut!="")
405  {
406  sigcut+="&&"+precut;
407  bgcut +="&&"+precut;
408  }
409 
410  t->SetEventList(0);
411 
412  float nsig = countEvents(t,sigcut);
413  float nbg = countEvents(t,bgcut);
414 
415  cout <<"SIG EVT : "<<nsig<<" "<<"BG EVT : "<<nbg<<endl;
416  cout <<"SIG EFF : "<<nsig/N0_sig<<" "<<"BG EFF : "<< nbg/N0_bg<<endl;
417  return 0;
418 }
N0_bg
double N0_bg
Definition: cutfinder.C:29
i
Int_t i
Definition: run_full.C:25
low
int low
Definition: anaMvdDigi.C:51
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
MAX
#define MAX
Definition: cutfinder.C:18
ev
int ev
Definition: GammaSpectraAnalysis.C:44
n
int n
Definition: HypStatDecay.cxx:21
h1
TH1F * h1
Definition: anal_hit_digi_cluster_fwendcap.C:121
evcnt
std::map< int, int > evcnt
Definition: cutfinder.C:22
bestSignificance
double bestSignificance(TH1 *h1, TH1 *h2, double &bestcut, double fac2=1.)
Definition: cutfinder.C:36
vars
TString vars[MAX]
Definition: cutfinder.C:26
f
TFile * f
Definition: bump_analys.C:12
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
c1
c1
Definition: plot_dirc.C:35
cut
double cut[MAX]
Definition: cutfinder.C:27
findLimits
void findLimits(TTree *t, TString var, TString ccut, double &low, double &high, double frac=0.98)
Definition: cutfinder.C:205
N0_sig
double N0_sig
Definition: cutfinder.C:29
idx
int idx[MAX]
Definition: cutfinder.C:25
nsig
int nsig
Definition: toy_core.C:46
cutfinder
int cutfinder(TString pre, TString ntp, TString sigcut, TString precut="", double n0s=50., double n0b=200., double supr=0.9)
Definition: cutfinder.C:239
h2
TH1F * h2
Definition: anal_hit_digi_cluster_fwendcap.C:122
bestSuppression
double bestSuppression(TH1 *h1, TH1 *h2, double &bestcut, double eff=0.9)
Definition: cutfinder.C:161
cnt
Int_t cnt
Definition: hist-t7.C:106
high
double high
Definition: full_core_ntp.C:35
bestEff
double bestEff(TH1 *h1, TH1 *h2, double &bestcut, double supr2=0.8, double fac2=1.)
Definition: cutfinder.C:118
t
TTree * t
Definition: bump_analys.C:13
signi
double signi[MAX]
Definition: cutfinder.C:24
diffQA
double diffQA(TH1F &ht1, TH1F &ht2)
Definition: cutfinder.C:99
countEvents
int countEvents(TTree *t, TString ccut)
Definition: cutfinder.C:74
mycompare
bool mycompare(int i, int j)
Definition: cutfinder.C:31
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