PWG3/vertexingHF/macros/plotHFMassFitterOutput.C

   1 #include <Riostream.h>
   2 #include <TH1F.h>
   3 #include <TF1.h>
   4 #include <TAttFill.h>
   5 #include <TFile.h>
   6 #include <TCanvas.h>
   7 #include <AliHFMassFitter.h>
   8
   9 //root[0] .L plotHFMassFitterOutput.C+
  10 //root[1] plotHFMassFitterOutput(2)
  11
  12 //input: pt bin you want to plot; path of the HFMassFitterOutput.root file
  13
  14 //three lines of fit:
  15 //grey= background with parameters from the fit on the side-bands
  16 //red= background with parameters from the total fit
  17 //blue= total fit
  18
  19 void plotHFMassFitterOutput(Int_t ptbin, TString path="./"){
  20   TString stdfilename="HFMassFitterOutput.root";
  21   stdfilename.Prepend(path);
  22
  23   cout<<"Opening "<<stdfilename<<endl;
  24   TFile *fin=new TFile(stdfilename.Data());
  25   if(!fin->IsOpen()){
  26     cout<<"File "<<stdfilename<<" not found"<<endl;
  27     return;
  28   }
  29   TString histoname="histMass_";
  30   histoname+=ptbin;
  31
  32   TH1F *h=(TH1F*)fin->Get(histoname);
  33   if(!h) {
  34     cout<<histoname<<" is not here, sorry!"<<endl;
  35     histoname="fhistoInvMass";
  36     h=(TH1F*)fin->Get(histoname);
  37     if(!h){
  38       cout<<histoname<<" is not here, sorry!"<<endl;
  39       cout<<"Write the name of the histo required (e.g. histMassSum<factor>_<ptbin>): ";
  40       cin>> histoname;
  41       h=(TH1F*)fin->Get(histoname);
  42       if(!h){
  43         cout<<histoname<<" is not in "<<stdfilename<<" check it, please!"<<endl;
  44           return;
  45       }
  46     } else {
  47       histoname+=ptbin;
  48       h->SetName(histoname);
  49     }
  50   }
  51
  52
  53     if(h){
  54     cout<<histoname<<" will be drawn!"<<endl;
  55     TCanvas *c1=new TCanvas("c1",histoname);
  56      //funcbkgfullrange_faint
  57     cout<<h<<endl;
  58     (h->GetFunction("funcbkgFullRange"))->SetLineColor(14);
  59     (h->GetFunction("funcbkgFullRange"))->SetLineStyle(4);
  60
  61     Double_t xmin,xmax;
  62     (h->GetFunction("funcbkgFullRange"))->GetRange(xmin,xmax);
  63     Int_t nfreepar=h->GetFunction("funcbkgFullRange")->GetNumberFreeParameters();
  64     cout<<"nfreepar = "<<nfreepar<<endl;
  65
  66     cout<<"Range = ("<<xmin<<", "<<xmax<<")"<<endl;
  67     cout<<"Bin Width "<<h->GetBinWidth(3)<<endl;
  68     cout<<"Initial parameters:\n";
  69     cout<<"par0= "<<h->GetFunction("funcbkgFullRange")->GetParameter(0)<<endl;
  70     cout<<"par1= "<<h->GetFunction("funcbkgFullRange")->GetParameter(1)<<endl;
  71     cout<<"par2= "<<h->GetFunction("funcbkgFullRange")->GetParameter(2)<<endl;
  72     cout<<"Formula= "<<h->GetFunction("funcbkgFullRange")->GetExpFormula()<<endl;
  73
  74     TF1 *fmass=h->GetFunction("funcmass");
  75     cout<<"Parametri massa input:\n";
  76     cout<<"par0= "<<h->GetFunction("funcmass")->GetParameter(0)-h->GetFunction("funcmass")->GetParameter(nfreepar)<<endl;
  77     cout<<"par1= "<<h->GetFunction("funcmass")->GetParameter(1)<<endl;
  78     cout<<"par2= "<<h->GetFunction("funcmass")->GetParameter(2)<<endl;
  79
  80
  81     AliHFMassFitter *fitter=new AliHFMassFitter(h,xmin,xmax,1,2,0);
  82     fitter->SetSideBands(kFALSE);
  83
  84     //funcbkgfullrange_deep
  85
  86     TF1 *fbkgfullrange_deep=new TF1("fbkgfullrange_deep",fitter,&AliHFMassFitter::FitFunction4Bkg,xmin,xmax,nfreepar,"AliHFMassFitter","FitFunction4Bkg");
  87
  88
  89     fbkgfullrange_deep->SetParameter(0,(fmass->GetParameter(0)-fmass->GetParameter(nfreepar))); //lin exp pol2 no-bkg
  90     if(nfreepar>=2) fbkgfullrange_deep->SetParameter(1,fmass->GetParameter(1)); //lin exp pol2
  91     if(nfreepar==3) fbkgfullrange_deep->SetParameter(2,fmass->GetParameter(2)); //pol2
  92
  93     fbkgfullrange_deep->SetMinimum(0);
  94
  95     cout<<"Final parameters:\n";
  96     cout<<"par0= "<<fbkgfullrange_deep->GetParameter(0)<<"\tcompare with "<<fbkgfullrange_deep->Integral(xmin,xmax)<<endl;
  97     cout<<"par1= "<<fbkgfullrange_deep->GetParameter(1)<<endl;
  98     cout<<"par2= "<<fbkgfullrange_deep->GetParameter(2)<<endl;
  99
 100     c1->cd();
 101     fbkgfullrange_deep->SetLineStyle(1);
 102     fbkgfullrange_deep->SetLineColor(2);
 103     h->Draw();
 104     fbkgfullrange_deep->Draw("sames");
 105     fmass->Draw("sames");
 106
 107   } else {
 108     cout<<histoname<<" is not here, sorry!"<<endl;
 109     return;
 110   }
 111 }