1 TCanvas *cCompareSimple,*cCompareCum;
2 TCanvas *cSubtraction,*cCompareSubtractSimple,*cCompareSubtractCum;
3 const Int_t timeWait=500;
5 void AnalyzeCharmFractionHists(Int_t ptbin,Int_t rebin){
8 // Macro for analyzing the impact parameter distribution of the D0 candidates.
10 // andrea.rossi@ts.infn.it
12 //==========================================================================
15 DrawHistos("d0D0_PureBack.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin);
16 SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin);
21 TH1F* CumulativeHist(const TH1F *h1,TString nameHist=0x0){
23 TString histname="hCumulative";
24 if(!nameHist.IsNull())histname.Append(nameHist.Data());
26 Double_t minX,maxX,binwidth,cumul=0.,errcumul=0.;
27 nbins=h1->GetNbinsX();
28 minX=h1->GetBinCenter(1);
29 maxX=h1->GetBinCenter(nbins);
30 binwidth=h1->GetBinWidth(1);
32 TH1F *h1copy=new TH1F();
36 TH1F *hCumul=new TH1F(histname.Data(),histname.Data(),nbins,minX-binwidth/2.,maxX+binwidth/2.);
39 for(Int_t j=1;j<=nbins;j++){
40 cumul+=h1copy->GetBinContent(j);
41 hCumul->SetBinContent(j,cumul);
44 hCumul->Scale(1./h1copy->Integral(1,nbins));
51 ///////////////////// HISTOGRAMS COMPARISON //////////////////////////////////
53 // Compare the desired histograms. The default case is the comparison
54 // between the d0 distribution for background (from MC) D0 candidates with the
55 // d0 distribution for "side bands" d0 candidates
57 ///////////////////////////////////////////////////////////////////////////////////
59 void DrawHistos(Int_t pt,Int_t rebin=1){
61 DrawHistos("","",hist,pt,rebin);
64 void DrawHistos(TString file1,TString file2,TString hist,Int_t pt,Int_t rebin=1){
71 Double_t integr1,integr2;
72 if(file1.IsNull())file1="d0D0_PureBack.root";
73 if(file2.IsNull())file2="d0D0NoMCSel_SideBand.root";
74 TFile *fSide=TFile::Open(file2.Data());
75 TFile *fmer=TFile::Open(file1.Data());
76 TH1F *h1=(TH1F*)fmer->Get(hist.Data());
77 nameH.Append("_Gaus");
78 h1->SetName(nameH.Data());
80 TH1F *h2=(TH1F*)fSide->Get(hist.Data());
81 nameH.Append("_SideBands");
82 h2->SetName(nameH.Data());
84 Bool_t satisfied=kFALSE;
89 TH1F *hCumul1=CumulativeHist(h1,"_Gaus");
90 hCumul1->SetDrawOption("E4");
91 TH1F *hCumul2=CumulativeHist(h2,"_SideBands");
92 hCumul2->SetDrawOption("E4");
95 TH1F *hDivCumul=new TH1F();
97 // hDivCumul->Sumw2();
98 hDivCumul->Divide(hCumul1,hCumul2);
99 hDivCumul->SetName("RatioCumulGausSdBands");
100 hDivCumul->SetTitle("Ratio Cumulative Gaus Over SideBands");
102 TH1F *hGausComp=new TH1F();
105 nameH.Append("Compare");
106 hGausComp->SetName(nameH.Data());
107 hGausComp->SetLineColor(kRed);
109 integr1=hGausComp->Integral();
111 TH1F *hSideComp=new TH1F();
114 nameH.Append("Compare");
115 hSideComp->SetName(nameH.Data());
116 hSideComp->SetLineColor(kBlue);
118 integr2=hSideComp->Integral();
119 if(integr2>integr1)hSideComp->Scale(integr1/integr2);
120 else hGausComp->Scale(integr2/integr1);
122 TH1F *hDivision=new TH1F();
124 hDivision->SetName("RatioGausSdBands");
125 hDivision->SetTitle("Ratio Gaus Over Sd Bands");
127 hDivision->Divide(hGausComp,hSideComp);
130 cCompareSimple=new TCanvas("cCompareSimple","Comparison of Under Gaus and Under Side Band background",700,700);
131 cCompareSimple->Divide(1,2);
132 cCompareSimple->cd(1);
134 hGausComp->Draw("Sames");
135 cCompareSimple->Update();
136 TPaveStats *p=hGausComp->FindObject("stats");
137 p->SetTextColor(kRed);
138 p=(TPaveStats*)hSideComp->FindObject("stats");
139 p->SetTextColor(kBlue);
140 cCompareSimple->cd(2);
142 cCompareSimple->Update();
144 cCompareCum=new TCanvas("cCompareCum","Comparison of cumulative function",700,700);
146 cCompareCum->Divide(1,2);
150 hCumul1->SetLineColor(kRed);
151 hCumul2->SetLineColor(kBlue);
152 hCumul1->SetFillColor(kRed);
153 hCumul2->SetFillColor(kBlue);
155 hCumul2->Draw("Sames");
158 cCompareCum->Update();
159 p=(TPaveStats*)hCumul1->FindObject("stats");
160 p->SetTextColor(kRed);
161 p=(TPaveStats*)hCumul2->FindObject("stats");
162 p->SetTextColor(kBlue);
163 cCompareCum->Update();
164 // cCompareCum->cd(2);
166 // cCompareCum->Update();
168 for(Int_t timewait=0;timewait<timeWait;timewait++){
169 cCompareCum->Update();
170 cCompareSimple->Update();
172 cout<<"Are you satisfied?"<<endl;
184 delete cCompareSimple;
197 /////////////////////////////////////// BACKGROUND SUBTRACTION EXAMPLE ///////////////////
200 // Performs: Signal Imp.Par Distr - B* Normalized Background Imp.Par Distr.(from Side Bands)
201 // and compare it with the MC Signal Imp. Par distribution.
202 // The amount of Background B is taken as the true one from the difference between
203 // the integrals of the "Observed Signal" and the MC Signal d0 distributions.
205 //////////////////////////////////////////////////////////////////////////////////////////////
207 void SubtractHist(Int_t pt,Int_t rebin=1){
208 SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",pt,rebin);
213 void SubtractHist(TString fileSignal,TString fileNoMC,TString fileNoMCSB,TString hist,Int_t pt,Int_t rebin){
221 Double_t integr1,integr2;
222 /* if(fileNoMC.IsNull())fileNoMC="d0D0merged.root";
223 if(fileNoMCSB.IsNull())fileNoMCSB="d0D0SideBmerged.root";*/
224 if(gSystem->AccessPathName(fileSignal.Data())){
225 Printf("Wrong signal file! \n");
228 if(gSystem->AccessPathName(fileNoMC.Data())){
229 Printf("Wrong d0distr under inv mass peak file! \n");
232 if(gSystem->AccessPathName(fileNoMCSB.Data())){
233 Printf("Wrong d0distr Side band file! \n");
237 TFile *fSide=TFile::Open(fileNoMCSB.Data());
238 TFile *fNoMCSignal=TFile::Open(fileNoMC.Data());
240 TH1F *h1=(TH1F*)fNoMCSignal->Get(hist.Data());
241 nameH.Append("_SelSignal");
242 h1->SetName(nameH.Data());
247 TH1F *h2=(TH1F*)fSide->Get(hist.Data());
248 nameH.Append("_SideBands");
249 h2->SetName(nameH.Data());
253 Double_t integrGaus,integrSB,integrSign;
254 integrGaus=h1->Integral();
255 integrSB=h2->Integral();
258 TFile *fSignal=TFile::Open(fileSignal.Data());
259 TH1F *hSign=(TH1F*)fSignal->Get(hist.Data());
261 nameH.Append("_MCSignal");
262 hSign->SetName(nameH.Data());
265 Double_t integrGaus,integrSB,integrSign;
266 integrGaus=h1->Integral();
267 integrSB=h2->Integral();
268 integrSign=hSign->Integral();
271 Bool_t satisfied=kFALSE;
279 TH1F *hGausSubtract=new TH1F();
282 hist.Append("_Subtracted");
283 hGausSubtract->SetName(hist.Data());
284 //hGausSubtract->Sumw2();
285 hGausSubtract->Add(h1,h2,1.,-1./integrSB*(integrGaus-integrSign));
286 /* hGausSubtract->Draw();
288 cSubtraction=new TCanvas("cSubtraction","cSubtraction",600,700);
290 hGausSubtract->SetLineColor(kRed);
291 hGausSubtract->Draw("E0");
292 hSign->SetLineColor(kBlue);
293 hSign->Draw("sames");
294 cSubtraction->Update();
295 TPaveStats *p=hGausSubtract->FindObject("stats");
296 p->SetTextColor(kRed);
297 cSubtraction->Update();
298 p=(TPaveStats*)hSign->FindObject("stats");
299 p->SetTextColor(kBlue);
300 cSubtraction->Update();
301 /*TString strText="Side Band integral: ";
303 TLatex *lat=new TLatex(500.,500.,strText.Data());// *text=p->AddText(strText.Data());
305 TPaveText *pvInfo = new TPaveText(72.74276,79.18782,92.84497,95.43147,"brNDC");
306 pvInfo->SetName("pvInfo");
307 pvInfo->SetBorderSize(2);
308 pvInfo->SetFillColor(19);
309 pvInfo->SetTextAlign(12);
310 TText *text=p->AddText(strText.Data());
312 hGausSubtract->GetListOfFunctions()->Add(pvInfo);
313 // pvInfo->SetParent(hGausSubtract->GetListOfFunctions());
315 // strText="MCSignal integral: ";
316 //strText+=integrSign;
317 text=p->AddText(strText.Data());
318 strText="Sel Signal integral: ";
320 text=p->AddText(strText.Data());
321 cSubtraction->Modified();*/
324 TH1F *hCumulGausSubtract=CumulativeHist(hGausSubtract,"_BackSubtr");
325 hCumulGausSubtract->SetDrawOption("E4");
326 TH1F *hCumulSign=CumulativeHist(hSign,"_MCSignal");
327 hCumulSign->SetDrawOption("E4");
330 TH1F *hDivCumul=new TH1F();
331 *hDivCumul=*hCumulSign;
332 hDivCumul->Divide(hCumulGausSubtract,hCumulSign);
333 hDivCumul->SetName("RatioCumulBackSubtr_MCSignal");
334 hDivCumul->SetTitle("Ratio Cumulative BackSubtracted Over MCSignal");
336 TH1F *hGausSubComp=new TH1F();
337 *hGausSubComp=*hGausSubtract;
338 nameH=hGausSubtract->GetName();
339 nameH.Append("Compare");
340 hGausSubComp->SetName(nameH.Data());
341 hGausSubComp->SetLineColor(kRed);
342 integr1=hGausSubComp->Integral();
344 TH1F *hSignComp=new TH1F();
346 nameH=hSign->GetName();
347 nameH.Append("Compare");
348 hSignComp->SetName(nameH.Data());
349 hSignComp->SetLineColor(kBlue);
350 integr2=hSignComp->Integral();
351 if(integr2>integr1)hSignComp->Scale(integr1/integr2);
352 else hGausSubComp->Scale(integr2/integr1);
354 TH1F *hDivision=new TH1F();
355 *hDivision=*hGausSubtract;
356 hDivision->SetName("RatioBackSubtr_Signal");
357 hDivision->SetTitle("Ratio BackSubtracted Over MCSignal");
358 hDivision->Divide(hGausSubComp,hSignComp);
361 cCompareSubtractSimple=new TCanvas("cCompareSubtractSimple","Comparison of BackSubtracted and MCSignal",700,700);
362 cCompareSubtractSimple->Divide(1,2);
363 cCompareSubtractSimple->cd(1);
365 hGausSubComp->Draw("Sames");
366 cCompareSubtractSimple->cd(2);
367 hDivision->SetLineColor(1);
369 cCompareSubtractSimple->Update();
370 p=(TPaveStats*)hGausSubComp->FindObject("stats");
371 p->SetTextColor(kRed);
372 p=(TPaveStats*)hSignComp->FindObject("stats");
373 p->SetTextColor(kBlue);
374 cCompareSubtractSimple->Update();
377 cCompareSubtractCum=new TCanvas("cCompareSubtractCumulative","Comparison of cumulative functions",700,700);
379 cCompareSubtractCum->Divide(1,2);
381 cCompareSubtractCum->cd(1);
383 hCumulGausSubtract->SetLineColor(kRed);
384 hCumulSign->SetLineColor(kBlue);
385 hCumulGausSubtract->SetFillColor(kRed);
386 hCumulSign->SetFillColor(kBlue);
387 hCumulGausSubtract->Draw();
388 hCumulSign->Draw("Sames");
389 cCompareSubtractCum->cd(2);
390 hDivCumul->SetLineColor(1);
393 cCompareSubtractCum->Update();
394 p=(TPaveStats*)hCumulGausSubtract->FindObject("stats");
395 p->SetTextColor(kRed);
396 p=(TPaveStats*)hCumulSign->FindObject("stats");
397 p->SetTextColor(kBlue);
398 cCompareSubtractCum->Update();
401 for(Int_t timewait=0;timewait<timeWait;timewait++){
402 cCompareSubtractCum->Update();
403 cCompareSubtractSimple->Update();
404 cSubtraction->Update();
407 cout<<"Are you satisfied?"<<endl;
413 delete hCumulGausSubtract;
416 delete cCompareSubtractCum;
419 delete cCompareSubtractSimple;
424 printf("Side Band integral: %d \n",integrSB);
425 printf("Selected Signal integral: %d \n",integrGaus);
426 printf("MC Signal integral: %d \n",integrSign);
427 printf(" -> Background = %d \n",integrGaus-integrSign);
433 //////////////////////////////// FIT DISTRIBUTION WITH GAUSSIAN + EXP TAILS ///////////////////
435 // Fit a istogram with a gaus(mean,sigma) + a*exp(mean2,lambda) function
437 /////////////////////////////////////////////////////////////////////////////////////////////////
438 void DoFit(TString filename,TString histo,TString gausSide,Int_t ptbin,Int_t rebin=-1){
440 TString histname=histo;
441 TString fileout=histo;
444 histname.Append("pt");
446 fileout.Append("pt");
449 else fileout.Append("ptAll");
451 TFile *fIN=TFile::Open(filename.Data());
453 TH1F *h=(TH1F*)fIN->Get(histname.Data());
456 Int_t binL=1,binR=h->GetNbinsX();
459 if(rebin>0)h->Rebin(rebin);
461 Double_t integr=h->Integral(binL,binR);
462 Double_t binwidth=h->GetBinWidth(10);
464 Int_t nbin=h->GetNbinsX();
467 TCanvas *c1=new TCanvas("c1","c1",700,600);
469 // h->SetFillColor(kRed);
474 TF1 *f=CreateFunc(integr,binwidth);
475 h->Fit(f->GetName(),"RI","",h->GetBinCenter(binL)-binwidth/2.,h->GetBinCenter(binR)+binwidth/2.);
478 cout<<"Is it ok?"<<endl;
483 fileout.Append(gausSide.Data());
484 fileout.Append("Fit.root");
486 TFile *fOUT=new TFile(fileout.Data(),"RECREATE");
495 cout<<"rebin value= ";
498 cout<<"Change Interval?"<<endl;
501 cout<<"Lower value= ";
504 cout<<"Upper value= ";
508 binL=h->FindBin(minX);
509 binR=h->FindBin(maxX);
526 TF1* CreateFunc(Double_t integral,Double_t binw){
527 TF1 *func=new TF1("func","[5]*[6]*((1.-[2])*1./2./[1]*TMath::Exp(-TMath::Abs((x-[0])/[1]))+[2]/TMath::Sqrt(2*TMath::Pi())/[4]*TMath::Exp(-1*(x-[3])*(x-[3])/2./[4]/[4]))");
529 func->FixParameter(5,integral);
530 func->SetParName(5,"histInt");
531 func->FixParameter(6,binw);
532 func->SetParName(6,"binW");
533 func->SetParLimits(0,-100,100.);
534 func->SetParName(0,"expMean");
535 func->SetParLimits(1,0.001,1000);
536 func->SetParName(1,"expDecay");
537 func->SetParLimits(2,0.00001,1.);
538 func->SetParName(2,"fracGaus");
539 func->SetParLimits(3,-100,100.);
540 func->SetParName(3,"gausMean");
541 func->SetParLimits(4,5.,200.);
542 func->SetParName(4,"gausSigma");
544 func->SetParameter(1,50.);
545 func->SetParameter(4,50.);
551 ////////////////////////////////// CHI2 TEST: ////////////////////////
554 ///////////////////////////////////////////////////////////////////////////////
555 Double_t GetChi2(const TH1F *h1,const TH1F *h2,TH1F *hchi2,Int_t &ndof){//ASSUME SAME BINNING
557 Int_t int1=h1->Integral();
558 Int_t int2=h2->Integral();
560 TH1F **h=new TH1F*[2];
561 TH1F *hchisq=new TH1F("hChi2","Chi2 histo",1000,0.,10.);
576 Int_t nbin=h[0]->GetNbinsX();//ASSUME SAME BINNING
578 Double_t chi2=0.,chi2Sum=0.;
580 for(Int_t j=1;j<=nbin;j++){//THE BINNING IS NOT TAKEN INTO ACCOUNT
582 if(h[1]->GetBinContent(j)==0||h[0]->GetBinContent(j)==0)continue;
583 fTrue=h[0]->GetBinContent(j)/nM;
584 chi2=(h[1]->GetBinContent(j)-fTrue*nm)*(h[1]->GetBinContent(j)-fTrue*nm)/h[1]->GetBinContent(j);
601 TH1F *h1=new TH1F("h1","h1",50,-10.,10.);
602 TH1F *h2=new TH1F("h2","h2",50,-10.,10.);
604 h1->FillRandom("gaus",50000);
605 h2->FillRandom("gaus",80000);
606 TH1F *hchi2=new TH1F();
608 Double_t chi2=GetChi2(h1,h2,hchi2,ndof);
610 cout<<"The chi2 is "<<chi2;