[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / AnalyzeCharmFractionHists.C

TCanvas *cCompareSimple,*cCompareCum;
TCanvas *cSubtraction,*cCompareSubtractSimple,*cCompareSubtractCum;
const Int_t timeWait=500;

void AnalyzeCharmFractionHists(Int_t ptbin,Int_t rebin){

  //
  // Macro for analyzing the impact parameter distribution of the D0 candidates.
  //  
  // andrea.rossi@ts.infn.it
  //
  //==========================================================================


  DrawHistos("d0D0_PureBack.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin);
  SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",ptbin,rebin);
  return;

}

TH1F* CumulativeHist(const TH1F *h1,TString nameHist=0x0){
  //ASSUME SAME BINNING
  TString histname="hCumulative";
  if(!nameHist.IsNull())histname.Append(nameHist.Data());
  Int_t nbins;
  Double_t minX,maxX,binwidth,cumul=0.,errcumul=0.;
  nbins=h1->GetNbinsX();
  minX=h1->GetBinCenter(1);
  maxX=h1->GetBinCenter(nbins);
  binwidth=h1->GetBinWidth(1);

  TH1F *h1copy=new TH1F();
  *h1copy=*h1;
  //h1copy->Sumw2();

  TH1F *hCumul=new TH1F(histname.Data(),histname.Data(),nbins,minX-binwidth/2.,maxX+binwidth/2.);


  for(Int_t j=1;j<=nbins;j++){
    cumul+=h1copy->GetBinContent(j);
    hCumul->SetBinContent(j,cumul);
  }
  hCumul->Sumw2();
  hCumul->Scale(1./h1copy->Integral(1,nbins));

  delete h1copy;
  return hCumul;
}


///////////////////// HISTOGRAMS COMPARISON   //////////////////////////////////
//              
//    Compare the desired histograms. The default case is the comparison
//    between the d0 distribution for background (from MC) D0 candidates with the
//    d0 distribution for "side bands" d0 candidates
//
///////////////////////////////////////////////////////////////////////////////////
             
void DrawHistos(Int_t pt,Int_t rebin=1){
  TString hist="hd0D0";
  DrawHistos("","",hist,pt,rebin);
  return;
}
void DrawHistos(TString file1,TString file2,TString hist,Int_t pt,Int_t rebin=1){
  
  if(pt>=0){
    hist.Append("pt");
    hist+=pt;
  }
  TString nameH=hist;
  Double_t integr1,integr2;
  if(file1.IsNull())file1="d0D0_PureBack.root";
  if(file2.IsNull())file2="d0D0NoMCSel_SideBand.root";
  TFile *fSide=TFile::Open(file2.Data());
  TFile *fmer=TFile::Open(file1.Data());
  TH1F *h1=(TH1F*)fmer->Get(hist.Data());
  nameH.Append("_Gaus");
  h1->SetName(nameH.Data());
  nameH=hist;
  TH1F *h2=(TH1F*)fSide->Get(hist.Data());
  nameH.Append("_SideBands");
  h2->SetName(nameH.Data());

  Bool_t satisfied=kFALSE;
  
  while(!satisfied){
    h1->Rebin(rebin);
    h2->Rebin(rebin);
    TH1F *hCumul1=CumulativeHist(h1,"_Gaus");
    hCumul1->SetDrawOption("E4");
    TH1F *hCumul2=CumulativeHist(h2,"_SideBands");
    hCumul2->SetDrawOption("E4");
    
    
    TH1F *hDivCumul=new TH1F();
    *hDivCumul=*hCumul2;
    //    hDivCumul->Sumw2();
    hDivCumul->Divide(hCumul1,hCumul2);
    hDivCumul->SetName("RatioCumulGausSdBands");
    hDivCumul->SetTitle("Ratio Cumulative Gaus Over SideBands");
    
    TH1F *hGausComp=new TH1F();
    *hGausComp=*h1;
    nameH=h1->GetName();
    nameH.Append("Compare");
    hGausComp->SetName(nameH.Data());
    hGausComp->SetLineColor(kRed);
    hGausComp->Sumw2();
    integr1=hGausComp->Integral();
    
    TH1F *hSideComp=new TH1F();
    *hSideComp=*h2;
    nameH=h2->GetName();
    nameH.Append("Compare");
    hSideComp->SetName(nameH.Data());
    hSideComp->SetLineColor(kBlue);
    hSideComp->Sumw2();
    integr2=hSideComp->Integral();
    if(integr2>integr1)hSideComp->Scale(integr1/integr2);
    else hGausComp->Scale(integr2/integr1);
    
    TH1F *hDivision=new TH1F();
    *hDivision=*h1;
    hDivision->SetName("RatioGausSdBands");
    hDivision->SetTitle("Ratio Gaus Over Sd Bands");
    hDivision->Sumw2();
    hDivision->Divide(hGausComp,hSideComp);
    
    
    cCompareSimple=new TCanvas("cCompareSimple","Comparison of Under Gaus and Under Side Band background",700,700);
    cCompareSimple->Divide(1,2);
    cCompareSimple->cd(1);
    hSideComp->Draw();
    hGausComp->Draw("Sames");
    cCompareSimple->Update();
    TPaveStats *p=hGausComp->FindObject("stats");
    p->SetTextColor(kRed);
    p=(TPaveStats*)hSideComp->FindObject("stats");
    p->SetTextColor(kBlue);
    cCompareSimple->cd(2);
    hDivision->Draw();
    cCompareSimple->Update();
    
    cCompareCum=new TCanvas("cCompareCum","Comparison of cumulative function",700,700);
    
    cCompareCum->Divide(1,2);
    
    cCompareCum->cd(1);
    
    hCumul1->SetLineColor(kRed);
    hCumul2->SetLineColor(kBlue);
    hCumul1->SetFillColor(kRed);
    hCumul2->SetFillColor(kBlue);
    hCumul1->Draw();
    hCumul2->Draw("Sames");
    cCompareCum->cd(2);
    hDivCumul->Draw();
    cCompareCum->Update();
    p=(TPaveStats*)hCumul1->FindObject("stats");
    p->SetTextColor(kRed);
    p=(TPaveStats*)hCumul2->FindObject("stats");
    p->SetTextColor(kBlue);
    cCompareCum->Update();
    // cCompareCum->cd(2);
    //hCumul2->Draw();
    //  cCompareCum->Update();
  
    for(Int_t timewait=0;timewait<timeWait;timewait++){
      cCompareCum->Update();
      cCompareSimple->Update();    
    }
    cout<<"Are you satisfied?"<<endl;
    cin>>satisfied;
    
    if(!satisfied){
     
    
      delete hCumul1;
      delete hCumul2;
      delete hDivCumul;
      delete cCompareCum;
      delete hSideComp;
      delete hGausComp;
      delete cCompareSimple;
      delete hDivision;
   
      cout<<"Rebin"<<endl;
      cin>>rebin;
      if(rebin<0)return;
    }
  }
  
  return;
  
}

///////////////////////////////////////   BACKGROUND SUBTRACTION EXAMPLE       ///////////////////
//
//
//   Performs: Signal Imp.Par Distr - B* Normalized Background Imp.Par Distr.(from Side Bands) 
//             and compare it with the MC Signal Imp. Par distribution. 
//             The amount of Background B is taken as the true one from the difference between 
//             the integrals of the "Observed Signal" and the MC Signal d0 distributions.
//           
//////////////////////////////////////////////////////////////////////////////////////////////

void SubtractHist(Int_t pt,Int_t rebin=1){
  SubtractHist("d0D0_Signal.root","d0D0NoMCSel.root","d0D0NoMCSel_SideBand.root","hd0D0",pt,rebin);
  
  return;
}

void SubtractHist(TString fileSignal,TString fileNoMC,TString fileNoMCSB,TString hist,Int_t pt,Int_t rebin){
  
  if(pt>=0){
    hist.Append("pt");
    hist+=pt;
  }
  
  TString nameH=hist;
  Double_t integr1,integr2;
  /*  if(fileNoMC.IsNull())fileNoMC="d0D0merged.root";
      if(fileNoMCSB.IsNull())fileNoMCSB="d0D0SideBmerged.root";*/
  if(gSystem->AccessPathName(fileSignal.Data())){
    Printf("Wrong signal file! \n");
    return;
  }
  if(gSystem->AccessPathName(fileNoMC.Data())){
    Printf("Wrong d0distr under inv mass peak file! \n");
    return;
  }
  if(gSystem->AccessPathName(fileNoMCSB.Data())){
    Printf("Wrong  d0distr Side band file! \n");
    return;
  }

  TFile *fSide=TFile::Open(fileNoMCSB.Data());
  TFile *fNoMCSignal=TFile::Open(fileNoMC.Data());
  
  TH1F *h1=(TH1F*)fNoMCSignal->Get(hist.Data());
  nameH.Append("_SelSignal");
  h1->SetName(nameH.Data());
  nameH=hist;
  h1->Sumw2();


  TH1F *h2=(TH1F*)fSide->Get(hist.Data());
  nameH.Append("_SideBands");
  h2->SetName(nameH.Data());
  h2->Sumw2();


  Double_t integrGaus,integrSB,integrSign;
  integrGaus=h1->Integral();
  integrSB=h2->Integral();
  

  TFile *fSignal=TFile::Open(fileSignal.Data());
  TH1F *hSign=(TH1F*)fSignal->Get(hist.Data());
  nameH=hist;
  nameH.Append("_MCSignal");
  hSign->SetName(nameH.Data());
  hSign->Sumw2();

  Double_t integrGaus,integrSB,integrSign;
  integrGaus=h1->Integral();
  integrSB=h2->Integral();
  integrSign=hSign->Integral();


   Bool_t satisfied=kFALSE;
   
   while(!satisfied){
     
     h1->Rebin(rebin);
     h2->Rebin(rebin);
     hSign->Rebin(rebin);
     
     TH1F *hGausSubtract=new TH1F();
     *hGausSubtract=*h1;
     nameH=hist;
     hist.Append("_Subtracted");
     hGausSubtract->SetName(hist.Data());
     //hGausSubtract->Sumw2();
     hGausSubtract->Add(h1,h2,1.,-1./integrSB*(integrGaus-integrSign));
     /*     hGausSubtract->Draw();
	    return;*/
     cSubtraction=new TCanvas("cSubtraction","cSubtraction",600,700);
     cSubtraction->cd();
     hGausSubtract->SetLineColor(kRed);
     hGausSubtract->Draw("E0"); 
     hSign->SetLineColor(kBlue);
     hSign->Draw("sames");
     cSubtraction->Update();
     TPaveStats *p=hGausSubtract->FindObject("stats");
     p->SetTextColor(kRed);
     cSubtraction->Update();
     p=(TPaveStats*)hSign->FindObject("stats");
     p->SetTextColor(kBlue);
     cSubtraction->Update();
     /*TString strText="Side Band integral: ";
       strText+=integrSB;
       TLatex *lat=new TLatex(500.,500.,strText.Data());// *text=p->AddText(strText.Data());
       lat->Draw();
       TPaveText *pvInfo = new TPaveText(72.74276,79.18782,92.84497,95.43147,"brNDC");
       pvInfo->SetName("pvInfo");
       pvInfo->SetBorderSize(2);
       pvInfo->SetFillColor(19);
       pvInfo->SetTextAlign(12);
       TText *text=p->AddText(strText.Data());
       pvInfo->Draw();
       hGausSubtract->GetListOfFunctions()->Add(pvInfo);
       //     pvInfo->SetParent(hGausSubtract->GetListOfFunctions());
       
       //     strText="MCSignal integral: ";
       //strText+=integrSign;
       text=p->AddText(strText.Data());
       strText="Sel Signal integral: ";
       strText+=integrGaus;
       text=p->AddText(strText.Data());
       cSubtraction->Modified();*/
  
     
     TH1F *hCumulGausSubtract=CumulativeHist(hGausSubtract,"_BackSubtr");
     hCumulGausSubtract->SetDrawOption("E4");
     TH1F *hCumulSign=CumulativeHist(hSign,"_MCSignal");
     hCumulSign->SetDrawOption("E4");
     
     
     TH1F *hDivCumul=new TH1F();
     *hDivCumul=*hCumulSign;
     hDivCumul->Divide(hCumulGausSubtract,hCumulSign);
     hDivCumul->SetName("RatioCumulBackSubtr_MCSignal");
     hDivCumul->SetTitle("Ratio Cumulative BackSubtracted Over MCSignal");
    
     TH1F *hGausSubComp=new TH1F();
     *hGausSubComp=*hGausSubtract;
     nameH=hGausSubtract->GetName();
     nameH.Append("Compare");
     hGausSubComp->SetName(nameH.Data());
     hGausSubComp->SetLineColor(kRed);
     integr1=hGausSubComp->Integral();
     
     TH1F *hSignComp=new TH1F();
     *hSignComp=*hSign;
     nameH=hSign->GetName();
     nameH.Append("Compare");
     hSignComp->SetName(nameH.Data());
     hSignComp->SetLineColor(kBlue);
     integr2=hSignComp->Integral();
     if(integr2>integr1)hSignComp->Scale(integr1/integr2);
     else hGausSubComp->Scale(integr2/integr1);
     
     TH1F *hDivision=new TH1F();
     *hDivision=*hGausSubtract;
     hDivision->SetName("RatioBackSubtr_Signal");
     hDivision->SetTitle("Ratio BackSubtracted Over MCSignal");
     hDivision->Divide(hGausSubComp,hSignComp);
     
     
     cCompareSubtractSimple=new TCanvas("cCompareSubtractSimple","Comparison of BackSubtracted and MCSignal",700,700);
     cCompareSubtractSimple->Divide(1,2);
     cCompareSubtractSimple->cd(1);
     hSignComp->Draw();
     hGausSubComp->Draw("Sames");
     cCompareSubtractSimple->cd(2);
     hDivision->SetLineColor(1);
     hDivision->Draw();
     cCompareSubtractSimple->Update();
     p=(TPaveStats*)hGausSubComp->FindObject("stats");
     p->SetTextColor(kRed);
     p=(TPaveStats*)hSignComp->FindObject("stats");
     p->SetTextColor(kBlue);
     cCompareSubtractSimple->Update();

     
     cCompareSubtractCum=new TCanvas("cCompareSubtractCumulative","Comparison of cumulative functions",700,700);
    
    cCompareSubtractCum->Divide(1,2);
    
    cCompareSubtractCum->cd(1);
    
    hCumulGausSubtract->SetLineColor(kRed);
    hCumulSign->SetLineColor(kBlue);
    hCumulGausSubtract->SetFillColor(kRed);
    hCumulSign->SetFillColor(kBlue);
    hCumulGausSubtract->Draw();
    hCumulSign->Draw("Sames");
    cCompareSubtractCum->cd(2);
    hDivCumul->SetLineColor(1);
    hDivCumul->Draw();

    cCompareSubtractCum->Update();
    p=(TPaveStats*)hCumulGausSubtract->FindObject("stats");
    p->SetTextColor(kRed);
    p=(TPaveStats*)hCumulSign->FindObject("stats");
    p->SetTextColor(kBlue);
    cCompareSubtractCum->Update();

    
    for(Int_t timewait=0;timewait<timeWait;timewait++){
      cCompareSubtractCum->Update();
      cCompareSubtractSimple->Update();  
      cSubtraction->Update();  
    }

    cout<<"Are you satisfied?"<<endl;
    cin>>satisfied;
    
    if(!satisfied){
      cout<<"Rebin"<<endl;
      cin>>rebin;
      delete hCumulGausSubtract;
      delete hCumulSign;
      delete hDivCumul;
      delete cCompareSubtractCum;
      delete hSignComp;
      delete hGausSubComp;
      delete cCompareSubtractSimple;
      delete hDivision;
    }
  }
  
   printf("Side Band integral: %d \n",integrSB);
   printf("Selected Signal integral: %d \n",integrGaus);
   printf("MC Signal integral: %d \n",integrSign);
   printf(" -> Background = %d \n",integrGaus-integrSign);
   
  return;
  
}

//////////////////////////////// FIT DISTRIBUTION WITH GAUSSIAN + EXP TAILS    ///////////////////
//
//            Fit a istogram with a gaus(mean,sigma) + a*exp(mean2,lambda) function
//
/////////////////////////////////////////////////////////////////////////////////////////////////
void DoFit(TString filename,TString histo,TString gausSide,Int_t ptbin,Int_t rebin=-1){

  TString  histname=histo;
  TString fileout=histo;

  if(ptbin>=0){
    histname.Append("pt");
    histname+=ptbin;
    fileout.Append("pt");
    fileout+=ptbin;
  }
  else fileout.Append("ptAll");
  
  TFile *fIN=TFile::Open(filename.Data());
  
  TH1F *h=(TH1F*)fIN->Get(histname.Data());

  Int_t ok=0;  
  Int_t binL=1,binR=h->GetNbinsX();
  
  while(ok!=1){
    if(rebin>0)h->Rebin(rebin);

    Double_t integr=h->Integral(binL,binR);
    Double_t binwidth=h->GetBinWidth(10);
    Double_t minX,maxX;
    Int_t nbin=h->GetNbinsX();
 
    
    TCanvas *c1=new TCanvas("c1","c1",700,600);
    c1->cd();
    //    h->SetFillColor(kRed);
    h->Draw("");
    //    c1->SetLogy();
    c1->Update();
  
    TF1 *f=CreateFunc(integr,binwidth);
    h->Fit(f->GetName(),"RI","",h->GetBinCenter(binL)-binwidth/2.,h->GetBinCenter(binR)+binwidth/2.);
    c1->Update();
    
    cout<<"Is it ok?"<<endl;
    cin>>ok;
    if(ok==1){
 
      fileout.Append("_");
      fileout.Append(gausSide.Data());
      fileout.Append("Fit.root");
      
      TFile *fOUT=new TFile(fileout.Data(),"RECREATE");
      fOUT->cd();
      c1->Write();
      h->Write();
      fOUT->Close();
      delete c1;
      delete f;
    }
    else if(ok==0){
      cout<<"rebin value= ";
      cin>>rebin;
      cout<<endl;
      cout<<"Change Interval?"<<endl;
      cin>>ok;
      if(ok){
	cout<<"Lower value= ";
	cin>>minX;
	cout<<endl;
	cout<<"Upper value= ";
	cin>>maxX;
	cout<<endl;
      
	binL=h->FindBin(minX);
	binR=h->FindBin(maxX);
      }
      ok=0;
      delete f;
      delete c1;
    }
    else {
      ok=1;
      delete f;
      delete c1;
    }
  }

  return;
}


TF1* CreateFunc(Double_t integral,Double_t binw){
  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]))");

  func->FixParameter(5,integral);
  func->SetParName(5,"histInt");
  func->FixParameter(6,binw);
  func->SetParName(6,"binW");
  func->SetParLimits(0,-100,100.);
  func->SetParName(0,"expMean");
  func->SetParLimits(1,0.001,1000);
  func->SetParName(1,"expDecay");
  func->SetParLimits(2,0.00001,1.);
  func->SetParName(2,"fracGaus");
  func->SetParLimits(3,-100,100.);
  func->SetParName(3,"gausMean");
  func->SetParLimits(4,5.,200.);
  func->SetParName(4,"gausSigma");

  func->SetParameter(1,50.);
  func->SetParameter(4,50.);
  
  return func;

}

////////////////////////////////// CHI2 TEST: ////////////////////////
//                               IN PROGRESS
//
///////////////////////////////////////////////////////////////////////////////
Double_t GetChi2(const TH1F *h1,const TH1F *h2,TH1F *hchi2,Int_t &ndof){//ASSUME SAME BINNING
  
  Int_t int1=h1->Integral();
  Int_t int2=h2->Integral();
  Int_t nm,nM;
  TH1F **h=new TH1F*[2];
  TH1F *hchisq=new TH1F("hChi2","Chi2 histo",1000,0.,10.);
  
  if(int2>int1){
    nM=int2;
    nm=int1;
    h[0]=h2;
    h[1]=h1;
  }
  else {
    nM=int1;
    nm=int2;
    h[0]=h1;
    h[1]=h2;
  }
  ndof=0;
  Int_t nbin=h[0]->GetNbinsX();//ASSUME SAME BINNING
  
  Double_t chi2=0.,chi2Sum=0.;
  Double_t f1,fTrue;
  for(Int_t j=1;j<=nbin;j++){//THE BINNING IS NOT TAKEN INTO ACCOUNT
    
    if(h[1]->GetBinContent(j)==0||h[0]->GetBinContent(j)==0)continue;
    fTrue=h[0]->GetBinContent(j)/nM;
    chi2=(h[1]->GetBinContent(j)-fTrue*nm)*(h[1]->GetBinContent(j)-fTrue*nm)/h[1]->GetBinContent(j);
    hchisq->Fill(chi2);
    chi2Sum+=chi2;
    ndof++;
  }
  
  *hchi2=*hchisq;
  delete hchisq;

  return chi2Sum;
  

}

void TestChi2(){


  TH1F *h1=new TH1F("h1","h1",50,-10.,10.);
  TH1F *h2=new TH1F("h2","h2",50,-10.,10.);
  
  h1->FillRandom("gaus",50000);
  h2->FillRandom("gaus",80000);
  TH1F *hchi2=new TH1F();
  Int_t ndof;
  Double_t chi2=GetChi2(h1,h2,hchi2,ndof);
  hchi2->Draw();
  cout<<"The chi2 is "<<chi2;
}
Commit	Line	Data
1282b9e5	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(2TMath::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)-fTruenm)(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	}