AliPhysics  d565ceb (d565ceb)
AnalyzeCharmFractionHists.C
Go to the documentation of this file.
3 const Int_t timeWait=500;
4 
6 
7  //
8  // Macro for analyzing the impact parameter distribution of the D0 candidates.
9  //
10  // andrea.rossi@ts.infn.it
11  //
12  //==========================================================================
13 
14 
15  DrawHistos("d0D0_PureBack.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin);
16  SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin);
17  return;
18 
19 }
20 
21 TH1F* CumulativeHist(const TH1F *h1,TString nameHist=0x0){
22  //ASSUME SAME BINNING
23  TString histname="hCumulative";
24  if(!nameHist.IsNull())histname.Append(nameHist.Data());
25  Int_t nbins;
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);
31 
32  TH1F *h1copy=new TH1F();
33  *h1copy=*h1;
34  //h1copy->Sumw2();
35 
36  TH1F *hCumul=new TH1F(histname.Data(),histname.Data(),nbins,minX-binwidth/2.,maxX+binwidth/2.);
37 
38 
39  for(Int_t j=1;j<=nbins;j++){
40  cumul+=h1copy->GetBinContent(j);
41  hCumul->SetBinContent(j,cumul);
42  }
43  hCumul->Sumw2();
44  hCumul->Scale(1./h1copy->Integral(1,nbins));
45 
46  delete h1copy;
47  return hCumul;
48 }
49 
50 
52 //
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
56 //
58 
60  TString hist="hd0D0";
61  DrawHistos("","",hist,pt,rebin);
62  return;
63 }
64 void DrawHistos(TString file1,TString file2,TString hist,Int_t pt,Int_t rebin=1){
65 
66  if(pt>=0){
67  hist.Append("pt");
68  hist+=pt;
69  }
70  TString nameH=hist;
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());
79  nameH=hist;
80  TH1F *h2=(TH1F*)fSide->Get(hist.Data());
81  nameH.Append("_SideBands");
82  h2->SetName(nameH.Data());
83 
84  Bool_t satisfied=kFALSE;
85 
86  while(!satisfied){
87  h1->Rebin(rebin);
88  h2->Rebin(rebin);
89  TH1F *hCumul1=CumulativeHist(h1,"_Gaus");
90  hCumul1->SetDrawOption("E4");
91  TH1F *hCumul2=CumulativeHist(h2,"_SideBands");
92  hCumul2->SetDrawOption("E4");
93 
94 
95  TH1F *hDivCumul=new TH1F();
96  *hDivCumul=*hCumul2;
97  // hDivCumul->Sumw2();
98  hDivCumul->Divide(hCumul1,hCumul2);
99  hDivCumul->SetName("RatioCumulGausSdBands");
100  hDivCumul->SetTitle("Ratio Cumulative Gaus Over SideBands");
101 
102  TH1F *hGausComp=new TH1F();
103  *hGausComp=*h1;
104  nameH=h1->GetName();
105  nameH.Append("Compare");
106  hGausComp->SetName(nameH.Data());
107  hGausComp->SetLineColor(kRed);
108  hGausComp->Sumw2();
109  integr1=hGausComp->Integral();
110 
111  TH1F *hSideComp=new TH1F();
112  *hSideComp=*h2;
113  nameH=h2->GetName();
114  nameH.Append("Compare");
115  hSideComp->SetName(nameH.Data());
116  hSideComp->SetLineColor(kBlue);
117  hSideComp->Sumw2();
118  integr2=hSideComp->Integral();
119  if(integr2>integr1)hSideComp->Scale(integr1/integr2);
120  else hGausComp->Scale(integr2/integr1);
121 
122  TH1F *hDivision=new TH1F();
123  *hDivision=*h1;
124  hDivision->SetName("RatioGausSdBands");
125  hDivision->SetTitle("Ratio Gaus Over Sd Bands");
126  hDivision->Sumw2();
127  hDivision->Divide(hGausComp,hSideComp);
128 
129 
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);
133  hSideComp->Draw();
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);
141  hDivision->Draw();
142  cCompareSimple->Update();
143 
144  cCompareCum=new TCanvas("cCompareCum","Comparison of cumulative function",700,700);
145 
146  cCompareCum->Divide(1,2);
147 
148  cCompareCum->cd(1);
149 
150  hCumul1->SetLineColor(kRed);
151  hCumul2->SetLineColor(kBlue);
152  hCumul1->SetFillColor(kRed);
153  hCumul2->SetFillColor(kBlue);
154  hCumul1->Draw();
155  hCumul2->Draw("Sames");
156  cCompareCum->cd(2);
157  hDivCumul->Draw();
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);
165  //hCumul2->Draw();
166  // cCompareCum->Update();
167 
168  for(Int_t timewait=0;timewait<timeWait;timewait++){
169  cCompareCum->Update();
170  cCompareSimple->Update();
171  }
172  cout<<"Are you satisfied?"<<endl;
173  cin>>satisfied;
174 
175  if(!satisfied){
176 
177 
178  delete hCumul1;
179  delete hCumul2;
180  delete hDivCumul;
181  delete cCompareCum;
182  delete hSideComp;
183  delete hGausComp;
184  delete cCompareSimple;
185  delete hDivision;
186 
187  cout<<"Rebin"<<endl;
188  cin>>rebin;
189  if(rebin<0)return;
190  }
191  }
192 
193  return;
194 
195 }
196 
198 //
199 //
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.
204 //
206 
208  SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",pt,rebin);
209 
210  return;
211 }
212 
213 void SubtractHist(TString fileSignal,TString fileNoMC,TString fileNoMCSB,TString hist,Int_t pt,Int_t rebin){
214 
215  if(pt>=0){
216  hist.Append("pt");
217  hist+=pt;
218  }
219 
220  TString nameH=hist;
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");
226  return;
227  }
228  if(gSystem->AccessPathName(fileNoMC.Data())){
229  Printf("Wrong d0distr under inv mass peak file! \n");
230  return;
231  }
232  if(gSystem->AccessPathName(fileNoMCSB.Data())){
233  Printf("Wrong d0distr Side band file! \n");
234  return;
235  }
236 
237  TFile *fSide=TFile::Open(fileNoMCSB.Data());
238  TFile *fNoMCSignal=TFile::Open(fileNoMC.Data());
239 
240  TH1F *h1=(TH1F*)fNoMCSignal->Get(hist.Data());
241  nameH.Append("_SelSignal");
242  h1->SetName(nameH.Data());
243  nameH=hist;
244  h1->Sumw2();
245 
246 
247  TH1F *h2=(TH1F*)fSide->Get(hist.Data());
248  nameH.Append("_SideBands");
249  h2->SetName(nameH.Data());
250  h2->Sumw2();
251 
252 
253  Double_t integrGaus,integrSB,integrSign;
254  integrGaus=h1->Integral();
255  integrSB=h2->Integral();
256 
257 
258  TFile *fSignal=TFile::Open(fileSignal.Data());
259  TH1F *hSign=(TH1F*)fSignal->Get(hist.Data());
260  nameH=hist;
261  nameH.Append("_MCSignal");
262  hSign->SetName(nameH.Data());
263  hSign->Sumw2();
264 
265  Double_t integrGaus,integrSB,integrSign;
266  integrGaus=h1->Integral();
267  integrSB=h2->Integral();
268  integrSign=hSign->Integral();
269 
270 
271  Bool_t satisfied=kFALSE;
272 
273  while(!satisfied){
274 
275  h1->Rebin(rebin);
276  h2->Rebin(rebin);
277  hSign->Rebin(rebin);
278 
279  TH1F *hGausSubtract=new TH1F();
280  *hGausSubtract=*h1;
281  nameH=hist;
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();
287  return;*/
288  cSubtraction=new TCanvas("cSubtraction","cSubtraction",600,700);
289  cSubtraction->cd();
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: ";
302  strText+=integrSB;
303  TLatex *lat=new TLatex(500.,500.,strText.Data());// *text=p->AddText(strText.Data());
304  lat->Draw();
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());
311  pvInfo->Draw();
312  hGausSubtract->GetListOfFunctions()->Add(pvInfo);
313  // pvInfo->SetParent(hGausSubtract->GetListOfFunctions());
314 
315  // strText="MCSignal integral: ";
316  //strText+=integrSign;
317  text=p->AddText(strText.Data());
318  strText="Sel Signal integral: ";
319  strText+=integrGaus;
320  text=p->AddText(strText.Data());
321  cSubtraction->Modified();*/
322 
323 
324  TH1F *hCumulGausSubtract=CumulativeHist(hGausSubtract,"_BackSubtr");
325  hCumulGausSubtract->SetDrawOption("E4");
326  TH1F *hCumulSign=CumulativeHist(hSign,"_MCSignal");
327  hCumulSign->SetDrawOption("E4");
328 
329 
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");
335 
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();
343 
344  TH1F *hSignComp=new TH1F();
345  *hSignComp=*hSign;
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);
353 
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);
359 
360 
361  cCompareSubtractSimple=new TCanvas("cCompareSubtractSimple","Comparison of BackSubtracted and MCSignal",700,700);
362  cCompareSubtractSimple->Divide(1,2);
363  cCompareSubtractSimple->cd(1);
364  hSignComp->Draw();
365  hGausSubComp->Draw("Sames");
366  cCompareSubtractSimple->cd(2);
367  hDivision->SetLineColor(1);
368  hDivision->Draw();
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();
375 
376 
377  cCompareSubtractCum=new TCanvas("cCompareSubtractCumulative","Comparison of cumulative functions",700,700);
378 
379  cCompareSubtractCum->Divide(1,2);
380 
381  cCompareSubtractCum->cd(1);
382 
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);
391  hDivCumul->Draw();
392 
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();
399 
400 
401  for(Int_t timewait=0;timewait<timeWait;timewait++){
402  cCompareSubtractCum->Update();
403  cCompareSubtractSimple->Update();
404  cSubtraction->Update();
405  }
406 
407  cout<<"Are you satisfied?"<<endl;
408  cin>>satisfied;
409 
410  if(!satisfied){
411  cout<<"Rebin"<<endl;
412  cin>>rebin;
413  delete hCumulGausSubtract;
414  delete hCumulSign;
415  delete hDivCumul;
416  delete cCompareSubtractCum;
417  delete hSignComp;
418  delete hGausSubComp;
419  delete cCompareSubtractSimple;
420  delete hDivision;
421  }
422  }
423 
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);
428 
429  return;
430 
431 }
432 
434 //
435 // Fit a istogram with a gaus(mean,sigma) + a*exp(mean2,lambda) function
436 //
438 void DoFit(TString filename,TString histo,TString gausSide,Int_t ptbin,Int_t rebin=-1){
439 
440  TString histname=histo;
441  TString fileout=histo;
442 
443  if(ptbin>=0){
444  histname.Append("pt");
445  histname+=ptbin;
446  fileout.Append("pt");
447  fileout+=ptbin;
448  }
449  else fileout.Append("ptAll");
450 
451  TFile *fIN=TFile::Open(filename.Data());
452 
453  TH1F *h=(TH1F*)fIN->Get(histname.Data());
454 
455  Int_t ok=0;
456  Int_t binL=1,binR=h->GetNbinsX();
457 
458  while(ok!=1){
459  if(rebin>0)h->Rebin(rebin);
460 
461  Double_t integr=h->Integral(binL,binR);
462  Double_t binwidth=h->GetBinWidth(10);
463  Double_t minX,maxX;
464  Int_t nbin=h->GetNbinsX();
465 
466 
467  TCanvas *c1=new TCanvas("c1","c1",700,600);
468  c1->cd();
469  // h->SetFillColor(kRed);
470  h->Draw("");
471  // c1->SetLogy();
472  c1->Update();
473 
474  TF1 *f=CreateFunc(integr,binwidth);
475  h->Fit(f->GetName(),"RI","",h->GetBinCenter(binL)-binwidth/2.,h->GetBinCenter(binR)+binwidth/2.);
476  c1->Update();
477 
478  cout<<"Is it ok?"<<endl;
479  cin>>ok;
480  if(ok==1){
481 
482  fileout.Append("_");
483  fileout.Append(gausSide.Data());
484  fileout.Append("Fit.root");
485 
486  TFile *fOUT=new TFile(fileout.Data(),"RECREATE");
487  fOUT->cd();
488  c1->Write();
489  h->Write();
490  fOUT->Close();
491  delete c1;
492  delete f;
493  }
494  else if(ok==0){
495  cout<<"rebin value= ";
496  cin>>rebin;
497  cout<<endl;
498  cout<<"Change Interval?"<<endl;
499  cin>>ok;
500  if(ok){
501  cout<<"Lower value= ";
502  cin>>minX;
503  cout<<endl;
504  cout<<"Upper value= ";
505  cin>>maxX;
506  cout<<endl;
507 
508  binL=h->FindBin(minX);
509  binR=h->FindBin(maxX);
510  }
511  ok=0;
512  delete f;
513  delete c1;
514  }
515  else {
516  ok=1;
517  delete f;
518  delete c1;
519  }
520  }
521 
522  return;
523 }
524 
525 
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]))");
528 
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");
543 
544  func->SetParameter(1,50.);
545  func->SetParameter(4,50.);
546 
547  return func;
548 
549 }
550 
552 // IN PROGRESS
553 //
555 Double_t GetChi2(const TH1F *h1,const TH1F *h2,TH1F *hchi2,Int_t &ndof){//ASSUME SAME BINNING
556 
557  Int_t int1=h1->Integral();
558  Int_t int2=h2->Integral();
559  Int_t nm,nM;
560  TH1F **h=new TH1F*[2];
561  TH1F *hchisq=new TH1F("hChi2","Chi2 histo",1000,0.,10.);
562 
563  if(int2>int1){
564  nM=int2;
565  nm=int1;
566  h[0]=h2;
567  h[1]=h1;
568  }
569  else {
570  nM=int1;
571  nm=int2;
572  h[0]=h1;
573  h[1]=h2;
574  }
575  ndof=0;
576  Int_t nbin=h[0]->GetNbinsX();//ASSUME SAME BINNING
577 
578  Double_t chi2=0.,chi2Sum=0.;
579  Double_t f1,fTrue;
580  for(Int_t j=1;j<=nbin;j++){//THE BINNING IS NOT TAKEN INTO ACCOUNT
581 
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);
585  hchisq->Fill(chi2);
586  chi2Sum+=chi2;
587  ndof++;
588  }
589 
590  *hchi2=*hchisq;
591  delete hchisq;
592 
593  return chi2Sum;
594 
595 
596 }
597 
598 void TestChi2(){
599 
600 
601  TH1F *h1=new TH1F("h1","h1",50,-10.,10.);
602  TH1F *h2=new TH1F("h2","h2",50,-10.,10.);
603 
604  h1->FillRandom("gaus",50000);
605  h2->FillRandom("gaus",80000);
606  TH1F *hchi2=new TH1F();
607  Int_t ndof;
608  Double_t chi2=GetChi2(h1,h2,hchi2,ndof);
609  hchi2->Draw();
610  cout<<"The chi2 is "<<chi2;
611 }
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
TCanvas * cCompareSubtractCum
TCanvas * cCompareSubtractSimple
TH1F * CumulativeHist(const TH1F *h1, TString nameHist=0x0)
TSystem * gSystem
TCanvas * cCompareSimple
Double_t GetChi2(const TH1F *h1, const TH1F *h2, TH1F *hchi2, Int_t &ndof)
TCanvas * cSubtraction
void DoFit(TString filename, TString histo, TString gausSide, Int_t ptbin, Int_t rebin=-1)
int Int_t
Definition: External.C:63
const Int_t timeWait
void DrawHistos(Int_t pt, Int_t rebin=1)
void TestChi2()
void SubtractHist(Int_t pt, Int_t rebin=1)
TCanvas * cCompareCum
void AnalyzeCharmFractionHists(Int_t ptbin, Int_t rebin)
Int_t rebin
TF1 * CreateFunc(Double_t integral, Double_t binw)
const Int_t nbins
bool Bool_t
Definition: External.C:53