--- /dev/null
+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;
+}
\ No newline at end of file
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff