FairRoot/PandaRoot
showvarexample.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 "TLegend.h"
18 
19 #define MAX 1000
20 #define BINS 500
21 #define NCAN 3
22 
23 std::map<int, int> evcnt;
24 
25 int idx[MAX];
26 TString vars[MAX];
27 
28 double N0_sig, N0_bg;
29 
30 // *** setup style
31 void setStyle()
32 {
33  gStyle->SetPadTopMargin(0.12);
34  gStyle->SetPadBottomMargin(0.12);
35  gStyle->SetPadLeftMargin(0.17);
36  gStyle->SetPadRightMargin(0.02);
37 
38  gStyle->SetLabelSize(0.075,"X");
39  gStyle->SetLabelSize(0.065,"Y");
40 
41  gStyle->SetStatY(0.87);
42  gStyle->SetStatX(0.98);
43  gStyle->SetStatW(0.4);
44  gStyle->SetStatH(0.12);
45  gStyle->SetOptStat(0);//"e"); // only entries
46 
47  gStyle->SetTitleH(0.1);
48  gStyle->SetTitleX(0.15);
49  gStyle->SetTitleOffset(1.1,"x");
50  gStyle->SetTitleXSize(0.076);
51  gStyle->SetNdivisions(505);
52 }
53 
54 
55 
56 int countEvents(TTree *t, TString ccut)
57 {
58  t->SetEventList(0);
59  t->SetBranchStatus("*",1);
60  t->Draw(">>el",ccut);
61  t->SetBranchStatus("*",0);
62  t->SetBranchStatus("evt",1);
63 
64  TEventList *el=(TEventList*)gDirectory->Get("el");
65 
66  Float_t ev;
67  t->SetBranchAddress("evt",&ev);
68 
69  evcnt.clear();
70 
71  for (int i=0;i<el->GetN();++i)
72  {
73  t->GetEntry(el->GetEntry(i));
74  evcnt[ev]+=1;
75  }
76  t->SetBranchStatus("*",1);
77 
78  return evcnt.size();
79 }
80 
81 double suppression(TH1* h1, TH1* h2, double cut)
82 {
83  int n=h1->GetNbinsX();
84  double sum1=0, sum2=0;
85  double int1 = h1->Integral();
86  double int2 = h2->Integral();
87 
88  int i=1;
89 
90  while (h1->GetBinCenter(i)<cut)
91  {
92  sum1 += h1->GetBinContent(i)/int1;
93  sum2 += h2->GetBinContent(i++)/int2;
94  }
95  return sum2;
96 }
97 void findLimits(TTree *t, TString var, TString ccut, double &low, double &high, double frac = 0.98)
98 {
99  low = t->GetMinimum(var);
100  high = t->GetMaximum(var);
101  return;
102 
103  double miss = (1.-frac)/2;
104 
105  t->SetEventList(0);
106  TEventList el("el");
107  t->Draw(">>el",ccut);
108  t->SetEventList(&el);
109 
110  double llow = t->GetMinimum(var);
111  double lhigh = t->GetMaximum(var);
112 
113  TH1F htemp("htemp","",500,llow, lhigh);
114  t->Project("htemp",var);
115 
116  double sum=0.;
117  double integ = htemp.Integral();
118 
119  int i=1;
120 
121  while (sum<miss) sum+=htemp.GetBinContent(i++)/integ;
122  low = htemp.GetBinCenter(i-2);
123  sum=0.; i=500;
124 
125  while (sum<miss) sum+=htemp.GetBinContent(i--)/integ;
126  high = htemp.GetBinCenter(i+2);
127 
128  t->SetEventList(0);
129 }
130 
131 void drawHistos(TTree *t, TString var, TString title, TString precut, TString newcut, double val, double val2=-999.)
132 {
133  double sigl, sigh, bgl, bgh;
134 
135  TLine l;
136  l.SetLineColor(1);
137  l.SetLineStyle(2);
138  l.SetLineWidth(2);
139 
140  TLatex lt;
141  lt.SetTextSize(0.06);
142  lt.SetTextColor(2);
143  TLatex lt2;
144  lt2.SetTextSize(0.06);
145  lt2.SetTextColor(kBlue+2);
146 
147 /* double preeffb = countEvents(t, "mode==900&&tag");
148  double preeffs = countEvents(t, "mode!=900&&tag");
149  preeffb/=500000.;
150  preeffs/=50000.;*/
151 
152 
153  double eff=countEvents(t, "mode==900&&tag&&"+precut+"&&"+newcut);
154  cout <<eff<<endl;
155  eff/=500000.;
156  double seff=countEvents(t, "mode!=900&&tag&&"+precut+"&&"+newcut);
157  cout <<seff<<endl;
158  seff/=50000.;
159 
160  TString sigcut = "mode!=900&&tag&&"+precut;
161  TString bgcut = "mode==900&&tag&&"+precut;
162 
163  findLimits(t,var,sigcut, sigl, sigh);
164  findLimits(t,var,bgcut, bgl, bgh);
165 
166  if (sigl>bgl) sigl=bgl;
167  if (sigh<bgh) sigh=bgh;
168 
169  if (var=="dtht") sigh=0.5;
170 
171  TH1F *h1=new TH1F("h1",title,200,sigl,sigh);
172  TH1F *h2=new TH1F("h2",title,200,sigl,sigh);
173  h2->SetLineColor(2);
174  h1->SetLineWidth(2);
175  h2->SetLineWidth(2);
176 
177  t->Project("h1",var, sigcut);
178  t->Project("h2",var, bgcut);
179 
180  double signi = suppression(h1,h2, val );
181 
182  h1->Scale(1.0/h1->Integral());
183  h2->Scale(1.0/h2->Integral());
184  h1->SetTitleSize(0.05);
185 
186  double maxi = h1->GetMaximum();
187  if (h2->GetMaximum()>maxi ) maxi = h2->GetMaximum();
188  maxi*=1.1;
189  h1->SetMaximum(maxi);
190  h2->SetMaximum(maxi);
191 
192  h1->DrawNormalized();
193  h2->DrawNormalized("same");
194 
195  double axmin = h1->GetXaxis()->GetXmin(),axmax = h1->GetXaxis()->GetXmax();
196 
197  l.DrawLine(val,0, val, maxi*0.9);
198  if (val2!=-999) l.DrawLine(val2,0, val2, maxi*0.9);
199 
200  //lt.DrawLatex(axmin+(axmax-axmin)*0.03,0.65*maxi,TString::Format("#epsilon_{BG,cut} = %5.2f%%",100*(1-signi)));
201  lt.DrawLatex(axmin+(axmax-axmin)*0.53,0.53*maxi,TString::Format("#epsilon_{bg} = %5.2f%%",100*eff));
202  lt2.DrawLatex(axmin+(axmax-axmin)*0.53,0.65*maxi,TString::Format("#epsilon_{sig} = %5.2f%%",100*seff));
203 
204  TLegend *leg1 = new TLegend(0.67,0.7,0.98,0.88);
205  leg1->AddEntry("h1","signal","l");
206  leg1->AddEntry("h2","BG","l");
207  leg1->Draw();
208 
209 }
210 
211 int showvarexample()
212 {
213  setStyle();
214 
215  int i,j;
216 
217  TFile *f=new TFile("data/M55_ndpm.root","READ");
218  //TFile *f=new TFile(fname,"READ");
219  TTree *t=(TTree*)f->Get("ndpm");
220 
221  TCanvas *c1=new TCanvas("c1","c1",10,10,1600,400);
222  c1->Divide(4,1);
223  TObjArray* branches = t->GetListOfBranches();
224 
225  TString bgcut = "mode==900 && tag";
226 
227 /* TString name1= "p(#Lambda) cms [GeV/c]";
228  TString var1 = "lampcm";
229  TString cut1 = "lampcm>1.8";
230  double val1 = 1.8;
231 
232  TString name2= "Thrust";
233  TString var2 = "thr";
234  TString cut2 = "thr>0.95";
235  double val2 = 0.95;
236 
237  TString name3= "p_{t}(#Lambda) [GeV/c]";
238  TString var3 = "lampt";
239  TString cut3 = "lampt>1";
240  double val3 = 1;*/
241 
242 
243 //dpcm>1.9&&dp>4&&ptmax>0.7&&dtht>0.09 *** sup 1000
244 
245  TString name1= "dpcm [GeV/c]";
246  TString var1 = "dpcm";
247  TString cut1 = "dpcm>1.9";
248  double val1 = 1.9;
249  double val11 =-999;
250 
251  TString name2= "dp [GeV/c]";
252  TString var2 = "dp";
253  TString cut2 = "dp>4";
254  double val2 = 4;
255 
256  TString name3= "ptmax [GeV/c]";
257  TString var3 = "ptmax";
258  TString cut3 = "ptmax>0.7";
259  double val3 = 0.7;
260 
261  TString name4= "dtht [rad]";
262  TString var4 = "dtht";
263  TString cut4 = "dtht>0.09";
264  double val4 = 0.09;
265 
266  //abs(dpcm-2.05)<0.2&&dp>2&&dpt>0.5&&ptmax>0.5 *** eff 90
267 
268 
269 // TString name1= "dpcm [GeV/c]";
270 // TString var1 = "dpcm";
271 // TString cut1 = "abs(dpcm-2.05)<0.2";
272 // double val1 = 1.85;
273 // double val11 =2.25;
274 //
275 // TString name2= "dp [GeV/c]";
276 // TString var2 = "dp";
277 // TString cut2 = "dp>2";
278 // double val2 = 2;
279 //
280 // TString name3= "ptmax [GeV/c]";
281 // TString var3 = "ptmax";
282 // TString cut3 = "ptmax>0.5";
283 // double val3 = 0.5;
284 //
285 // TString name4= "dpt [GeV/c]";
286 // TString var4 = "dpt";
287 // TString cut4 = "dpt>0.5";
288 // double val4 = 0.5;
289 
290  double preEff = countEvents(t, bgcut);
291  cout <<preEff<<endl;
292  preEff/=500000.;
293 
294 
295 /* t->SetEventList(0);
296  TEventList el1("el1");
297  t->Draw(">>el1",precut);
298  t->SetEventList(&el1);*/
299 
300  c1->cd(1);
301  drawHistos(t, var1,name1, "1",cut1,val1,val11);
302  c1->cd(2);
303  drawHistos(t, var2,name2, cut1,cut2,val2);
304  c1->cd(3);
305  drawHistos(t, var3,name3, cut1+"&&"+cut2,cut3,val3);
306  c1->cd(4);
307  drawHistos(t, var4,name4, cut1+"&&"+cut2+"&&"+cut3,cut4,val4);
308  c1->Update();
309 
310 
311 
312  cout<<100*preEff<<endl;
313 
314  return 0;
315 }
N0_bg
double N0_bg
Definition: showvarexample.C:28
name3
TString name3
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
name2
TString name2
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
countEvents
int countEvents(TTree *t, TString ccut)
Definition: showvarexample.C:56
i
Int_t i
Definition: run_full.C:25
low
int low
Definition: anaMvdDigi.C:51
drawHistos
void drawHistos(TTree *t, TString var, TString title, TString precut, TString newcut, double val, double val2=-999.)
Definition: showvarexample.C:131
val
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
MAX
#define MAX
Definition: showvarexample.C:19
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
cut
double cut[MAX]
Definition: autocutx.C:36
findLimits
void findLimits(TTree *t, TString var, TString ccut, double &low, double &high, double frac=0.98)
Definition: showvarexample.C:97
N0_sig
double N0_sig
Definition: showvarexample.C:28
name4
TString name4
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
name1
TString name1
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
f
TFile * f
Definition: bump_analys.C:12
c1
c1
Definition: plot_dirc.C:35
suppression
double suppression(TH1 *h1, TH1 *h2, double cut)
Definition: showvarexample.C:81
h2
TH1F * h2
Definition: anal_hit_digi_cluster_fwendcap.C:122
vars
TString vars[MAX]
Definition: showvarexample.C:26
evcnt
std::map< int, int > evcnt
Definition: showvarexample.C:23
high
double high
Definition: full_core_ntp.C:35
setStyle
void setStyle()
Definition: showvarexample.C:31
t
TTree * t
Definition: bump_analys.C:13
showvarexample
int showvarexample()
Definition: showvarexample.C:211
idx
int idx[MAX]
Definition: showvarexample.C:25
signi
double signi[MAX]
Definition: cutfinderx.C:35
htemp
TH1D * htemp[nsteps]
Definition: dedx_bands.C:61
Generated on Sun Aug 18 2019 03:02:02 for FairRoot/PandaRoot by   doxygen 1.8.8