.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / PWGPP / pid / doeffPi.C
CommitLineData
7484404e 1#include"TF1.h"
2#include"TH1D.h"
3#include"TH2F.h"
4#include"TMath.h"
5#include"TSystem.h"
6#include"TCanvas.h"
7#include"TFile.h"
8#include"TGraphErrors.h"
9#include"AliPIDperfContainer.h"
cd9f95f6 10#include"TRandom.h"
7484404e 11
12Int_t LoadLib();
13void doeffPi(Int_t pos=1,Float_t prob=0.1,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8);
14TH2F *GetHistoPip(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkp=0,Float_t pMinkn=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8);
15TH2F *GetHistoPin(Float_t pt=1,Float_t ptM=1.1,Float_t pMinkn=0,Float_t pMinkp=0.,Float_t etaminkp=-0.8,Float_t etamaxkp=0.8);
bf2e05d5 16void fit(TH1D *h,Float_t *a=NULL,char *opt="",char *opt2="",Float_t pt=1.5);
7484404e 17void AddHisto(TH2F *h1,TH2F *h2,Float_t w);
18
a8ad4709 19TObject* fContPid1;
20TObject* fContPid2;
21const Int_t nBinPid = 14; // pt,eta, ptPip, ptPin, PPip, PPin, TOF3sigmaPip, TOF3sigmaPin, isPhiTrue, nsigmaPip, nsigmaPin
22// 0.985 < mass < 1.045 (60) and 0 < centrality < 100 (10)
d82096df 23Int_t binPid[nBinPid] = {1/*ptPhi*/,1/*EtaPi*/,20/*pt+*/,20/*pt-*/,5/*P+*/,5/*P-*/,2/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
a8ad4709 24Float_t xmin[nBinPid] = {1,-0.8,0.3,0.3,0,0,-0.5,-0.5,-0.5,0,0,-TMath::Pi(),-TMath::Pi(),-TMath::Pi()/2};
25Float_t xmax[nBinPid] = {5,0.8,4.3,4.3,1,1,1.5,1.5,1.5,7.5,7.5,TMath::Pi(),TMath::Pi(),TMath::Pi()/2};
26
27TF1 *fsign;
28TF1 *fall;
29TF1 *fback;
30
31Int_t ifunc=0;
32
33Float_t fitmin = 0.3;
34Float_t fitmax = 0.7;
35
08b5b548 36Int_t cmin = 1; // min 1
37Int_t cmax = 10;//max 10
a8ad4709 38
6490e37e 39Float_t weightS = -1.;
a8ad4709 40
7484404e 41Int_t rebinsize = 1;
a8ad4709 42
43Int_t parplotted = 2;
44
45Bool_t isMC = kFALSE; // don't change this (is set automatically)
46Bool_t selectTrue = kTRUE; // put it to true to remove background (only for MC)
7ae044be 47Bool_t keepTrue = kTRUE; // put it to false to fit only background (only for MC)
a8ad4709 48
49Bool_t kGoodMatch = kFALSE; // to check good matching
50
51Bool_t kSigma2vs3 = kFALSE; // to check good matching
52
d82096df 53Bool_t require5sigma = kFALSE; // don't touch this flag
a8ad4709 54
55Bool_t bayesVsigma = kFALSE; // only to do checks
56
57Bool_t kTOFmatch = kFALSE; // for combined PID requires TOF matching
58
08b5b548 59Bool_t kOverAll = kFALSE;
e34b28fe 60Bool_t kOverAllTOFmatch = kFALSE;
08b5b548 61Bool_t kOverAll2Sigma = kFALSE;
e34b28fe 62
167b8b0a 63Bool_t kPid2Sigma = kFALSE;
64Bool_t kPid3Sigma = kFALSE;
65
08b5b548 66TH2F *hmatched;
67TH2F *htracked;
a8ad4709 68
69Bool_t kLoaded=kFALSE;
7484404e 70
71Int_t LoadLib(){
6490e37e 72 weightS = -1.;
a8ad4709 73
bf2e05d5 74 require5sigma = kFALSE;
75
a8ad4709 76 if(! kLoaded){
4070f709 77 gSystem->Load("libVMC");
78 gSystem->Load("libPhysics");
79 gSystem->Load("libTree");
80 gSystem->Load("libMinuit");
81 gSystem->Load("libSTEERBase");
82 gSystem->Load("libANALYSIS");
83 gSystem->Load("libAOD");
84 gSystem->Load("libESD");
85 gSystem->Load("libANALYSIS");
86 gSystem->Load("libANALYSISalice");
87 gSystem->Load("libCORRFW");
88 gSystem->Load("libNetx");
89 gSystem->Load("libPWGPPpid");
a8ad4709 90
91 TFile *f = new TFile("AnalysisResults.root");
a8ad4709 92 TList *l = (TList *) f->Get("contK0sBayes1");
08b5b548 93 TList *l2 = (TList *) f->Get("contK0sBayes2");
7484404e 94
95 if(! (l && l2)) return 0;
96
08b5b548 97 fContPid1 = (AliPIDperfContainer *) l->FindObject("contPID");
98 fContPid2 = (AliPIDperfContainer *) l->FindObject("contPID2");
99 hmatched = (TH2F *) l2->FindObject("hMatchPi");
100 htracked = (TH2F *) l2->FindObject("hTrackingPi");
a8ad4709 101 }
102 kLoaded = kTRUE;
103
104 // check if MC
476e3d92
FN
105 Float_t x[] = {xmin[0]+0.001,xmin[1]+0.001,xmin[2]+0.001,xmin[3]+0.001,xmin[4]+0.001,xmin[5]+0.001,xmin[6]+0.001,xmin[7]+0.001,1/*trueMC*/,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
106 Float_t x2[] = {xmax[0],xmax[1],xmax[2],xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],xmax[8],xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
a8ad4709 107
7484404e 108 AliPIDperfContainer *tmp = (AliPIDperfContainer *) fContPid1;
a8ad4709 109 TH1D *h = tmp->GetQA(0, x, x2)->ProjectionX("checkMC");
110
111 if(h->GetEntries()) isMC = kTRUE;
112 else isMC=kFALSE;
113
114 if(!isMC){
115 selectTrue = kFALSE;
116 keepTrue = kTRUE;
117 }
118 else{
119 printf("MC truth found!!!!!!\nIt is MC!!!!!!");
120 }
121
bf2e05d5 122 fsign = new TF1("fsign","gaus(0) +0.5*[0]*TMath::Exp(-[3]*TMath::Abs(x-[1]))",fitmin,fitmax);
123 fback = new TF1("fback","pol2",fitmin,fitmax);
124 fall = new TF1("fall","gaus(0) +0.5*[0]*TMath::Exp(-[3]*TMath::Abs(x-[1])) + pol2(4)",fitmin,fitmax);
a8ad4709 125
126 fsign->SetLineColor(2);
127 fback->SetLineColor(4);
128
129 if(kSigma2vs3){
130 kGoodMatch=kFALSE;
131 kOverAll = 0;
132 }
133
134 if(bayesVsigma){
135 kOverAll = 0;
136 kGoodMatch=kFALSE;
137 kSigma2vs3=kFALSE;
138 kTOFmatch=kTRUE;
139 weightS = -0.7;
140 }
7484404e 141
142 return 1;
a8ad4709 143}
144
7484404e 145void doeffPi(Int_t pos,Float_t prob,Float_t etaminkp,Float_t etamaxkp){
a8ad4709 146 LoadLib();
08b5b548 147 TH1D *hm = hmatched->ProjectionX("matchingPiEff",cmin,cmax);
148 TH1D *ht = htracked->ProjectionX("tracking",cmin,cmax);
149
150 hm->GetYaxis()->SetTitle("TOF matching eff.");
151 hm->SetTitle("Using probability as weights");
152
153 hm->Sumw2();
154 ht->Sumw2();
155
156 hm->Divide(hm,ht,1,1,"B");
a8ad4709 157
158 Int_t nptbin = binPid[2];
159 Float_t minptbin = xmin[2];
160 Float_t maxptbin = xmax[2];
161
162 if(pos == 0){
163 nptbin = binPid[3];
164 minptbin = xmin[3];
165 maxptbin = xmax[3];
166 }
167
167b8b0a 168 if(prob > 0.1999|| kPid3Sigma ||kPid2Sigma){
a8ad4709 169 kGoodMatch = kFALSE;
170 kSigma2vs3 = kFALSE;
80d8d480 171// if(! kOverAll) require5sigma = kTRUE;
6490e37e 172 if(!isMC) weightS = -0.95;
a8ad4709 173 }
174
7484404e 175 TCanvas *c1 = new TCanvas();
176 c1->Divide((nptbin+1)/2,2);
177 TH2F *hh,*hh2;
178 TH1D *h;
a8ad4709 179 char name[100];
180 Float_t b[50][3];
181
182 Double_t xx[50],yy[50];
183 Double_t exx[50],eyy[50];
184
185 for(Int_t i=0;i < nptbin;i++){
7484404e 186 c1->cd(i+1);//->SetLogy();
a8ad4709 187 Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
188 Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
189
190 xx[i] = (ptmin+ptmax)/2;
191 exx[i] = (-ptmin+ptmax)/2;
192
193 Float_t pp=0.1;
194 if(prob < 0.2) pp = 0.;
195 if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
196 else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
197 sprintf(name,"TOF matched: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
198 hh->SetTitle(name);
199 sprintf(name,"hNoPid%i",i);
200
201 pp=prob;
202 if(prob < 0.2) pp = 0.1;
203 if(pos) hh2=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
204 else hh2=GetHistoPin(ptmin,ptmax,pp,0.0);
205 AddHisto(hh,hh2,weightS);
206
207 h = hh->ProjectionX(name,cmin,cmax);
208 h->RebinX(rebinsize);
209 h->Draw("ERR");
210 h->SetMarkerStyle(24);
211 b[i][0]=-1;
212 Int_t ntrial = 0;
213 Float_t chi2 = 10000;
08b5b548 214 while(ntrial < 3 && (chi2 > 20 + 1000*selectTrue)){
bf2e05d5 215 fit(h,b[i],"WW","",xx[i]);
a8ad4709 216 c1->Update();
217// getchar();
bf2e05d5 218 fit(h,b[i],"","",xx[i]);
a8ad4709 219 ntrial++;
220 chi2 = b[i][2];
221 printf("chi2 = %f\n",chi2);
222 c1->Update();
223// getchar();
224
225 }
226
227 yy[i] = fall->GetParameter(parplotted);
228 eyy[i] = fall->GetParError(parplotted);
229 }
230
231 TGraphErrors *gpar = new TGraphErrors(nptbin,xx,yy,exx,eyy);
7484404e 232 c1->cd(8);
08b5b548 233// gpar->Draw("AP");
a8ad4709 234 gpar->SetMarkerStyle(20);
235
236 TCanvas *c2 = new TCanvas();
237 c2->Divide((nptbin+1)/2,2);
238 Float_t b2[50][3];
239
240 for(Int_t i=0;i < nptbin;i++){
241 c2->cd(i+1);
242 Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
243 Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
244
245 Float_t pp=prob;
246 if(prob < 0.2) pp = 0.1;
247 if(pos) hh=GetHistoPip(ptmin,ptmax,pp,0.0,etaminkp,etamaxkp);
248 else hh=GetHistoPin(ptmin,ptmax,pp,0.0);
249 sprintf(name,"P_{TOF} > 0.8: %f < p_{T} < %f GeV/#it{c}",ptmin,ptmax);
250 hh->SetTitle(name);
251 sprintf(name,"hPid60_%i",i);
252 h = hh->ProjectionX(name,cmin,cmax);
253 h->RebinX(rebinsize);
254 h->Draw("ERR");
255 h->SetMarkerStyle(24);
256 b2[i][0]=-1;
257 Int_t ntrial = 0;
258 Float_t chi2 = 10000;
08b5b548 259 while(ntrial < 3 && (chi2 > 20 + 1000*selectTrue)){
bf2e05d5 260 fit(h,b2[i],"WW","",xx[i]);
261 fit(h,b2[i],"","",xx[i]);
a8ad4709 262 ntrial++;
263 chi2 = b2[i][2];
264 printf("chi2 = %f\n",chi2);
265 }
266 yy[i] = fall->GetParameter(parplotted);
267 eyy[i] = fall->GetParError(parplotted);
268
269 }
270
271 TGraphErrors *gpar2 = new TGraphErrors(nptbin,xx,yy,exx,eyy);
272 c2->cd(8);
08b5b548 273// gpar2->Draw("AP");
a8ad4709 274 gpar2->SetMarkerStyle(20);
275
276 Double_t xpt[50],expt[50],eff[50],efferr[50];
277 for(Int_t i=0;i<nptbin;i++){
278 printf("%f +/- %f - %f +/- %f\n",b[i][0],b[i][1],b2[i][0],b2[i][1]);
279
280 Float_t ptmin = minptbin+(maxptbin-minptbin)/nptbin*(i);
281 Float_t ptmax = minptbin+(maxptbin-minptbin)/nptbin*(i+1);
282
283 xpt[i] = (ptmin+ptmax)/2;
284 expt[i] = (-ptmin+ptmax)/2;
285 eff[i] = b2[i][0]/(b[i][0]-b2[i][0]*weightS);
286
ca4825ae 287 // b[i][0] = b[i][0]-b2[i][0]*weightS;
a8ad4709 288
ca4825ae 289 // efferr[i] = TMath::Sqrt(b[i][1]*b[i][1]/b[i][0]/b[i][0] + b2[i][1]*b2[i][1]/b2[i][0]/b2[i][0])*(b2[i][0]+b2[i][1])*(1+weightS*(b2[i][0]-b2[i][1])/b[i][0])/b[i][0];//*(1-eff[i]);//der2*der2*(b[i][1]*b[i][1] - b2[i][1]*b2[i][1]));
290
291 efferr[i] = 1./(b[i][0]-b2[i][0]*weightS)/(b[i][0]-b2[i][0]*weightS)*TMath::Sqrt(b[i][0]*b[i][0]*b2[i][1]*b2[i][1] + b2[i][0]*b2[i][0]*b[i][1]*b[i][1]);
a8ad4709 292
293 if(TMath::Abs(efferr[i]) > 1)efferr[i]=1;
294 }
295 new TCanvas();
296 TGraphErrors *geff = new TGraphErrors(nptbin,xpt,eff,expt,efferr);
297 geff->Draw("AP");
298
299 char flag[100];
7484404e 300 flag[0] = '\0';
a8ad4709 301
302 if(isMC){
303 if(selectTrue) sprintf(flag,"true");
304 else if(!keepTrue) sprintf(flag,"back");
305 }
306
08b5b548 307 Bool_t kWriteME = kFALSE;
308
a8ad4709 309 char flag2[100];
7484404e 310 flag2[0] = '\0';
a8ad4709 311
312 char etarange[100];
313 sprintf(etarange,"_%.1f-%.1f_",etaminkp,etamaxkp);
314
315 if(kGoodMatch)
316 sprintf(flag2,"GM");
317
318 if(bayesVsigma)
319 sprintf(flag2,"BayesVsSigma");
320
321 if(kSigma2vs3)
322 sprintf(flag2,"Sigma2vs3");
323
324 if(kOverAll)
325 sprintf(flag2,"OverAll");
e34b28fe 326 if(kOverAllTOFmatch)
327 sprintf(flag2,"OverAllTOF");
328 if(kOverAll2Sigma)
329 sprintf(flag2,"OverAll2sigma");
a8ad4709 330
167b8b0a 331 if(kPid3Sigma)
332 sprintf(flag2,"pid3sigma");
333 if(kPid2Sigma)
334 sprintf(flag2,"pid2sigma");
335
336
a8ad4709 337 if(pos){
338 if(prob >=0.2) sprintf(name,"pionPos%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
08b5b548 339 else{
340 sprintf(name,"pionPos%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
6490e37e 341 if(!(kOverAll || bayesVsigma || kGoodMatch || kSigma2vs3)) kWriteME = kTRUE;
08b5b548 342 }
a8ad4709 343 }
344 else{
345 if(prob >=0.2) sprintf(name,"pionNeg%sP%iEff%i_%i%s%s.root",etarange,Int_t(prob*100),(cmin-1)*10,cmax*10,flag,flag2);
346 else sprintf(name,"pionNeg%sMatchEff%i_%i%s%s.root",etarange,(cmin-1)*10,cmax*10,flag,flag2);
347 }
348
08b5b548 349 geff->SetTitle("#pi efficiency (from K^{0}_{s});p_{T} (GeV/#it{c};efficiency");
a8ad4709 350 TFile *fout = new TFile(name,"RECREATE");
351 geff->Write();
08b5b548 352 if(kWriteME) hm->Write();
a8ad4709 353 fout->Close();
354
08b5b548 355 if(kWriteME) hm->Draw("SAME");
a8ad4709 356}
357
7484404e 358TH2F *GetHistoPip(Float_t pt,Float_t ptM,Float_t pMinkp,Float_t pMinkn,Float_t etaminkp,Float_t etamaxkp){
a8ad4709 359
167b8b0a 360 Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,pt+0.001,xmin[3]+0.001,pMinkp+0.001,pMinkn+0.001,(pMinkp>0.09 || kPid3Sigma||kPid2Sigma)+0.001,kTOFmatch+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
a8ad4709 361 Float_t x2[] = {xmax[0],etamaxkp-0.001,ptM-0.001,xmax[3],xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
362
e34b28fe 363 if(kOverAll){
364 x[6] = 0.0001;
365 x2[9] = 5.9;
366 if(pMinkp > 0.19) x2[9] = 4.9;
367 }
368
369 if(kOverAllTOFmatch && pMinkp > 0.19){
370 x[6] = 1.0001;
6490e37e 371 x2[9] = 4.9;
e34b28fe 372 }
373
374 if(kOverAll2Sigma && pMinkp > 0.09){
08b5b548 375 x2[9] = 2.;
e34b28fe 376 x[6] = 1.0001;
377 }
378
a8ad4709 379 if(kGoodMatch){
380 x[6] = 1.0001;
381 if(pMinkp > 0)
382 x2[9] = 4.9;
383
384 }
385
386 if(kTOFmatch){
387 x[6] = 1.0001;
388 }
389
390 if(kSigma2vs3){
391 x[6] = 1.0001;
392 x2[9] = 3;
393 if(pMinkp > 0)
394 x2[9] = 2;
395 }
396
397 if(bayesVsigma){
398 if(pMinkp > 0){
399 x[4] = 0.2001;
400 x2[9] = 5;
401 }
402 else{
403 x2[9] = 3;
404 }
405
406
407 }
408
409 if(require5sigma) x2[9] = 4.9;
167b8b0a 410 if(kPid3Sigma && pMinkp>0.09) x2[9] = 2.9;
411 if(kPid2Sigma && pMinkp>0.09) x2[9] = 1.9;
a8ad4709 412
d82096df 413 printf("max sigma = %f\n",x2[9]);
414
7484404e 415 AliPIDperfContainer *tmp = (AliPIDperfContainer *) fContPid1;
a8ad4709 416
417 TH2F *h = tmp->GetQA(0, x, x2);
418
419 h->GetXaxis()->SetTitle("M_{K^{0}_{s}} (GeV/#it{c}^{2})");
420 h->GetYaxis()->SetTitle("centrality [%]");
421
422 return h;
423}
424
7484404e 425TH2F *GetHistoPin(Float_t pt,Float_t ptM,Float_t pMinkn,Float_t pMinkp,Float_t etaminkp,Float_t etamaxkp){
a8ad4709 426
167b8b0a 427 Float_t x[] = {xmin[0]+0.001,etaminkp+0.001,xmin[2]+0.001,pt+0.001,pMinkp+0.001,pMinkn+0.001,kTOFmatch+0.001,(pMinkn>0.09 || kPid3Sigma|| kPid2Sigma)+0.001,selectTrue,xmin[9],xmin[10],xmin[11],xmin[12],xmin[13]};
a8ad4709 428 Float_t x2[] = {xmax[0],etamaxkp-0.001,xmax[2],ptM-0.001,xmax[4],xmax[5],xmax[6],xmax[7],keepTrue,xmax[9],xmax[10],xmax[11],xmax[12],xmax[13]};
429
e34b28fe 430 if(kOverAll){
431 x[7] = 0.0001;
432 x2[10] = 5.9;
433 if(pMinkn > 0.19) x2[10] = 4.9;
434 }
435
436 if(kOverAllTOFmatch && pMinkn > 0.19){
437 x[7] = 1.0001;
6490e37e 438 x2[10] = 4.9;
e34b28fe 439 }
440
441 if(kOverAll2Sigma && pMinkn > 0.09){
442 x2[10] = 2;
443 x[7] = 1.0001;
444 }
445
a8ad4709 446 if(kGoodMatch){
447 x[7] = 1.0001;
448 if(pMinkn > 0)
449 x2[10] = 4.9;
450
451 }
452
453 if(kTOFmatch){
454 x[7] = 1.0001;
455 }
456
457 if(kSigma2vs3){
458 x[7] = 1.0001;
459 x2[10] = 3;
460 if(pMinkn > 0)
461 x2[10] = 2;
462 }
463
464 if(bayesVsigma){
465 if(pMinkn > 0){
466 x[5] = 0.2001;
467 x2[10] = 5;
468 }
469 else{
470 x2[10] = 3;
471 }
472 }
473
474 if(require5sigma) x2[10] = 4.9;
475
167b8b0a 476 if(kPid3Sigma && pMinkn>0.09) x2[10] = 2.9;
477 if(kPid2Sigma && pMinkn>0.09) x2[10] = 1.9;
478
d82096df 479 printf("max sigma = %f\n",x2[10]);
480
7484404e 481 AliPIDperfContainer *tmp = (AliPIDperfContainer *) fContPid2;
a8ad4709 482
483 TH2F *h = tmp->GetQA(0, x, x2);
484
485 h->GetXaxis()->SetTitle("M_{K^{0}_{s}} (GeV/#it{c}^{2})");
486 h->GetYaxis()->SetTitle("centrality [%]");
487
488 return h;
489}
490
bf2e05d5 491void fit(TH1D *h,Float_t *a,char *opt,char *opt2,Float_t pt){
a8ad4709 492 if(h->Integral(1,h->GetNbinsX()) < 1){
493 if(a){
494 a[0]=0.001;
495 a[1]=1;
496 }
497 return;
498 }
499
500
501 fall->SetParameter(0,100);
502 fall->SetParameter(1,0.4971);
bf2e05d5 503 fall->SetParameter(2,2.89748e-03);
504 fall->FixParameter(3,230+30/pt);
a8ad4709 505
bf2e05d5 506 fall->SetParLimits(0,0.00001,1000000);
a8ad4709 507 fall->SetParLimits(1,0.4965,0.4985);
bf2e05d5 508 fall->SetParLimits(2,0.0025,0.005);
509 //fall->SetParLimits(3,200,350);
a8ad4709 510
511 fall->ReleaseParameter(4);
512 fall->ReleaseParameter(5);
513
514 if(selectTrue){
515 fall->FixParameter(4,0);
516 fall->FixParameter(5,0);
bf2e05d5 517 fall->FixParameter(6,0);
a8ad4709 518 }
7ae044be 519 fall->FixParameter(6,0);
520
a8ad4709 521
cd9f95f6 522 char namenew[100];
523 sprintf(namenew,"%s_%i",h->GetName(),Int_t(gRandom->Rndm()*10000));
524 TH1D *h2 = new TH1D(*h);
525 h2->SetName(namenew);
526
476e3d92 527// Float_t entries = h2->GetBinContent(h2->FindBin(0.497));
cd9f95f6 528// printf("entries under the peak = %f, pt = %f\n",entries,pt);
529// getchar();
530
531 if(pt > 2.5){
532 if(pt < 2.8) h2->RebinX(2);
533 else if(pt < 3) h2->RebinX(4);
9309a88b 534 else h2->RebinX(4);
cd9f95f6 535 }
536
537 h=h2;
538
a8ad4709 539 char name[100];
540 TF1 *ftmp=fall;
541
542 TF1 *ftmp2=new TF1(*fsign);
543 sprintf(name,"fsign%i",ifunc);
544 ftmp2->SetName(name);
545
546 TF1 *ftmp3=new TF1(*fback);
547 sprintf(name,"ftmp3%i",ifunc);
548 ftmp3->SetName(name);
549
550 ifunc++;
551
552 h->Fit(ftmp,opt,opt2,fitmin,fitmax);
553 h->Draw("ERR");
554
555 ftmp2->SetParameter(0,ftmp->GetParameter(0));
556 ftmp2->SetParameter(1,ftmp->GetParameter(1));
557 ftmp2->SetParameter(2,ftmp->GetParameter(2));
558 ftmp2->SetParameter(3,ftmp->GetParameter(3));
559 ftmp2->Draw("SAME");
560 ftmp3->SetParameter(0,ftmp->GetParameter(4));
561 ftmp3->SetParameter(1,ftmp->GetParameter(5));
bf2e05d5 562 ftmp3->SetParameter(2,ftmp->GetParameter(6));
a8ad4709 563 ftmp3->Draw("SAME");
564
565 Float_t mean = ftmp->GetParameter(1);
bf2e05d5 566 Float_t sigma = TMath::Abs(ftmp->GetParameter(2));
a8ad4709 567
568 Float_t signI = ftmp2->Integral(mean-10*sigma,mean+10*sigma)/h->GetBinWidth(1);
569 Float_t backI = ftmp3->Integral(mean-3*sigma,mean+3*sigma)/h->GetBinWidth(1);
570
d82096df 571 if(signI < 0) signI = 0;
572 if(backI < 1) backI = 1;
573
a8ad4709 574 Float_t errI = TMath::Sqrt(ftmp->GetParError(0)*ftmp->GetParError(0)/(0.001+ftmp->GetParameter(0))/(0.001+ftmp->GetParameter(0)));
575
576 printf("signal(5 sigma) = %f +/- %f(fit) +/- %f(stat)\n",signI,errI*signI,TMath::Sqrt(signI));
577 printf("backgr(3sigma) = %f\n",backI);
578 printf("significance(3 sigma) = %f\n",signI/sqrt(signI+backI));
579
580 if(a){
581 a[0]=signI;
582 a[1]=signI*errI*signI*errI + signI;
583 a[1] = TMath::Sqrt(a[1]);
584 if(ftmp->GetNDF()) a[2] = ftmp->GetChisquare()/ftmp->GetNDF();
585
586
587 if(selectTrue){
588 a[0] = h->Integral(1,h->GetNbinsX());
589 a[1] = TMath::Sqrt(a[0]);
590 }
591 }
592}
593
7484404e 594void AddHisto(TH2F *h1,TH2F *h2,Float_t w){
a8ad4709 595 Int_t nbinx = h1->GetNbinsX();
596 Int_t nbiny = h1->GetNbinsY();
597
598 for(Int_t i=1;i<=nbinx;i++){
599 for(Int_t j=1;j<=nbiny;j++){
7484404e 600 Double_t val = h1->GetBinContent(i,j) + h2->GetBinContent(i,j)*w;
a8ad4709 601 Float_t err = TMath::Min(TMath::Sqrt(val),val);
602 h1->SetBinContent(i,j,val);
603 h1->SetBinError(i,j,err);
604 }
605 }
606}