PWG3/vertexingHF/macros/AnalyzeCharmFractionHists.C

   1 TCanvas *cCompareSimple,*cCompareCum;
   2 TCanvas *cSubtraction,*cCompareSubtractSimple,*cCompareSubtractCum;
   3 const Int_t timeWait=500;
   4
   5 void AnalyzeCharmFractionHists(Int_t ptbin,Int_t rebin){
   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
  51 ///////////////////// HISTOGRAMS COMPARISON   //////////////////////////////////
  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 //
  57 ///////////////////////////////////////////////////////////////////////////////////
  58
  59 void DrawHistos(Int_t pt,Int_t rebin=1){
  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
 197 ///////////////////////////////////////   BACKGROUND SUBTRACTION EXAMPLE       ///////////////////
 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 //
 205 //////////////////////////////////////////////////////////////////////////////////////////////
 206
 207 void SubtractHist(Int_t pt,Int_t rebin=1){
 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
 433 //////////////////////////////// FIT DISTRIBUTION WITH GAUSSIAN + EXP TAILS    ///////////////////
 434 //
 435 //            Fit a istogram with a gaus(mean,sigma) + a*exp(mean2,lambda) function
 436 //
 437 /////////////////////////////////////////////////////////////////////////////////////////////////
 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
 551 ////////////////////////////////// CHI2 TEST: ////////////////////////
 552 //                               IN PROGRESS
 553 //
 554 ///////////////////////////////////////////////////////////////////////////////
 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 }