PWG3/vertexingHF/macros/fitD0InvMass.C

   1 #include <Riostream.h>
   2 #include <TStyle.h>
   3 #include <TFile.h>
   4 #include <TList.h>
   5 #include <TDirectoryFile.h>
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TCanvas.h>
   9
  10 #include <AliHFMassFitter.h>
  11
  12 #include <fstream>
  13 #include <cmath>
  14
  15 //root[0] .x fitD0InvMass.C+
  16 //root[1] fitD0(...)
  17
  18 //input: nsigma is the number of sigma in which significance, signal and background are calculated
  19
  20 //fit of histograms in the directory PWG3_D2H_D0InvMass of AnalysisResults.root . Must specify the list name where the histo are stored. Produce a root file with the fitted histos and a text file with signal, background and significance
  21
  22 void fitD0(Int_t rebin=0,TString listname="coutputmassD0mycuts",Int_t nsigma=3, TString pathin="./",TString pathout="./",Int_t btype=2){
  23
  24   TString file="AnalysisResults.root"; //"D0InvMass.root";
  25
  26   file.Prepend(pathin);
  27   cout<<"Opening "<<file<<endl;
  28   TFile *fin=new TFile(file.Data());
  29   if(!fin->IsOpen()){
  30     cout<<"File "<<file.Data()<<" not found"<<endl;
  31     return;
  32   }
  33
  34   gStyle->SetOptFit(0111);
  35   gStyle->SetOptStat("nemrou");
  36
  37   const Int_t npt=5;
  38   TString dirname="PWG3_D2H_D0InvMass";
  39   TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname);
  40
  41   TList *lista = (TList*) dir->Get(listname.Data());
  42   if(!lista) {
  43     cout<<listname<<" doesn't exist"<<endl;
  44     return;
  45   }
  46
  47   TString fileout="table.dat";
  48   fileout.Prepend(pathout);
  49   ofstream out(fileout.Data());
  50
  51   AliHFMassFitter *fitter=new AliHFMassFitter();
  52
  53   cout<<"mean = "<<fitter->GetMean()<<"\tsigma= "<<fitter->GetSigma()<<endl;
  54   Bool_t init=kTRUE;
  55
  56   out<<"Bakground type = "<<btype<<";\trebin = "<<rebin<<endl;
  57   Int_t i=0;
  58   for(Int_t ipt=0;ipt<npt;ipt++){ //npt
  59
  60     out<<"************************************************************\n";
  61     out<<"* ptbin * entriesMass * entriesSgn * entriesBkg * entriesRfl *\n";
  62
  63     //ALL w cuts
  64     TString nameall="histMass_";
  65     nameall+=(ipt+1);
  66     TH1F *hMass=(TH1F*)lista->FindObject(nameall);
  67     if(!hMass){
  68       cout<<nameall<<" not found"<<endl;
  69       return;
  70     }
  71
  72     out<<"* "<< ipt+1 <<"  *  "<<hMass->GetEntries()<<"  *  ";
  73     //SIGNAL from MC
  74     TString namesgn="histSgn_";
  75     namesgn+=(ipt+1);
  76     TH1F *hSgn=(TH1F*)lista->FindObject(namesgn);
  77     if(!hSgn){
  78       cout<<namesgn<<" not found"<<endl;
  79       return;
  80     }
  81
  82     out<<hSgn->GetEntries()<<"  *  ";
  83     //BACKGROUND from MC
  84     TString namebkg="histBkg_";
  85     namebkg+=(ipt+1);
  86     TH1F *hBkg=(TH1F*)lista->FindObject(namebkg);
  87     if(!hBkg){
  88       cout<<namebkg<<" not found"<<endl;
  89       return;
  90     }
  91
  92     out<<hBkg->GetEntries()<<"  *  ";
  93     //REFLECTED SIGNAL from MC
  94     TString namerfl="histRfl_";
  95     namerfl+=(ipt+1);
  96     TH1F *hRfl=(TH1F*)lista->FindObject(namerfl);
  97     if(!hRfl){
  98       cout<<namerfl<<" not found"<<endl;
  99       return;
 100     }
 101
 102     out<<hRfl->GetEntries()<<"  *\n";
 103     out<<"*************************************************\n";
 104     if(rebin!=0){
 105       hSgn->Rebin(rebin);
 106       hBkg->Rebin(rebin);
 107       hRfl->Rebin(rebin);
 108     }
 109     Double_t min, max;
 110     Int_t    nbin;
 111     Double_t sgn,errsgn,errsgn2=0,bkg,errbkg,errbkg2=0,sgnf,errsgnf,sgnfMC;
 112     Double_t meanfrfit,sigmafrfit, meanMC=hSgn->GetMean(),sigmaMC=hSgn->GetRMS();
 113     Double_t width=0;
 114     nbin=hMass->GetNbinsX();
 115     min=hMass->GetBinLowEdge(1);
 116     max=min+nbin*hMass->GetBinWidth(nbin);
 117
 118     //FIT:
 119
 120     TString namentu="ntupt3bin";
 121     if(init) fitter->InitNtuParam((char*)namentu.Data());
 122     // - all
 123     fitter->SetHisto(hMass);
 124     fitter->SetRangeFit(min, max);
 125     //fitter->SetRangeFit(1.83,1.89);
 126     fitter->SetType(btype,0);//(b,s)
 127     if(rebin!=0) fitter->RebinMass(rebin);
 128     Bool_t fitOK=fitter->MassFitter(kFALSE); //kFALSE = do not draw
 129     if(!fitOK) {
 130       out<<"Fit return kFALSE, skip "<<hMass->GetName()<<endl;
 131       fitter->Reset(); //delete histogram set
 132       continue;
 133     }
 134     width=fitter->GetHistoClone()->GetBinWidth(3);
 135     cout<<"\nChi^2 = "<<fitter->GetChiSquare()<<"\t Reduced Chi^2 = "<<fitter->GetReducedChiSquare()<<endl;
 136     meanfrfit=fitter->GetMean();
 137     sigmafrfit=fitter->GetSigma();
 138     out<<"mean = "<<meanfrfit<<" (MC "<<meanMC<<")"<<"\tsigma= "<<sigmafrfit<<" (MC "<<sigmaMC<<")"<<endl;
 139     //cout<<"old nbin = "<<hMass->GetNbinsX()<<"\tnew nbin = "<<fitter->GetHistoClone()->GetNbinsX()<<endl;
 140     if(meanfrfit<0 || sigmafrfit<0 || meanfrfit<min || meanfrfit>max) {
 141       cout<<"Fit failed, check"<<endl;
 142       out<<hMass->GetName();
 143       out<<": \nSgn not available"<<endl;
 144       out<<"Bkg  not available"<<endl;
 145       out<<"Sgnf not available"<<endl;
 146       out<<"*************************************************\n";
 147     }
 148     else{
 149       //cout<<"Writing..."<<endl;
 150       fitter->WriteHisto(pathout);
 151       hMass= fitter->GetHistoClone();
 152       //TH1F *hc=fitter->GetHistoClone();
 153       //TF1 *fbtest=hc->GetFunction("funcbkgRecalc"); //new version of ALiHFMassFitter
 154
 155       fitter->FillNtuParam();
 156
 157       Double_t limsx,limdx;
 158       limsx=meanfrfit-nsigma*sigmafrfit;
 159       limdx=meanfrfit+nsigma*sigmafrfit;
 160       // limsx=meanMC-nsigma*sigmaMC;
 161       //limdx=meanMC+nsigma*sigmaMC;
 162
 163
 164       //determine limit of nsigma in bins
 165       Int_t binsx,bindx;
 166       binsx=hMass->FindBin(limsx);
 167       if (limsx > hMass->GetBinCenter(binsx)) binsx++;
 168       bindx=hMass->FindBin(limdx);
 169       if (limdx < hMass->GetBinCenter(bindx)) bindx--;
 170
 171       //reconvert bin in x
 172       Double_t sxr,dxr;
 173       sxr=hMass->GetBinLowEdge(binsx);
 174       dxr=hMass->GetBinLowEdge(bindx+1);
 175
 176       fitter->Signal(sxr,dxr,sgn,errsgn);
 177       fitter->Background(sxr,dxr,bkg,errbkg);
 178       fitter->Significance(sxr,dxr,sgnf,errsgnf);
 179
 180
 181       Float_t inttot,intsgn,intsgnerr;
 182
 183       Int_t np=-99;
 184       switch (btype){
 185       case 0: //expo
 186         np=2;
 187         break;
 188       case 1: //linear
 189         np=2;
 190         break;
 191       case 2: //pol2
 192         np=3;
 193         break;
 194       case 3: //no bkg
 195         np=1;
 196         break;
 197       }
 198
 199
 200       TF1 *fmass=hMass->GetFunction("funcmass");
 201       if (fmass){
 202
 203         inttot=fmass->GetParameter(0);
 204         intsgn=fmass->GetParameter(np);
 205         intsgnerr=fmass->GetParError(np);
 206         //cout<<"i = "<<i<<"inttot = "<<inttot<<"\tintsgn = "<<intsgn<<"\tintsgnErr = "<<intsgnerr<<"\twidth = "<<width<<endl;
 207         Double_t errbkg2rel=errbkg/bkg*TMath::Sqrt((inttot-intsgn)/width/bkg);//intsgnerr/(inttot-intsgn)*TMath::Sqrt((inttot-intsgn)/width/bkg);
 208         Double_t errsgn2rel=errsgn/sgn*TMath::Sqrt(intsgn/width/sgn);
 209         errbkg2=errbkg2rel*bkg;
 210         errsgn2=errsgn2rel*sgn;
 211       }
 212
 213
 214       cout<<"bin sx = "<<binsx<<"\t xsx = "<<sxr<<endl;
 215       cout<<"bin dx = "<<bindx<<"\t xdx = "<<dxr<<endl;
 216
 217       out<<hMass->GetName();
 218       Double_t sgnMC,bkgMC,rflMC;
 219       sgnMC=hSgn->Integral(binsx,bindx);
 220       bkgMC=hBkg->Integral(binsx,bindx);
 221       rflMC=hRfl->Integral(binsx,bindx);
 222       sgnfMC=sgnMC/TMath::Sqrt(sgnMC+bkgMC+rflMC);
 223       out<<": \nSgn  "<<sgn<<"+/-"<<errsgn<<" ("<<errsgn2<<")\tCompare with: "<<sgnMC<<endl;
 224       out<<"Bkg  "<<bkg<<"+/-"<<errbkg<<" ("<<errbkg2<<")\tCompare with: "<<bkgMC<<" + "<<rflMC<<" = "<<bkgMC+rflMC<<endl;
 225       out<<"Sgnf "<<sgnf<<"+/-"<<errsgnf<<"\tCompare with: "<<sgnfMC<<endl;
 226       if (sgn != 0 && sgnfMC != 0) out<<"sigmaS/S = "<<errsgn/sgn<<"\tCompare with 1/signif: "<<1./sgnfMC<<endl;
 227       else {
 228         out<<"sigmaS/S = "<<errsgn<<"/"<<sgn<<"\tCompare with 1/signif: 1./"<<sgnfMC<<endl;
 229       }
 230       out<<"nsigma considered for comparison = \ndx: "<<(dxr-meanfrfit)/sigmafrfit<<"\nsx: "<<(meanfrfit-sxr)/sigmafrfit<<endl;
 231       out<<"Mean = "<<meanfrfit<<"\tSigma = "<<sigmafrfit<<endl;
 232       out<<"*************************************************\n";
 233       i++;
 234     }
 235
 236     sgn=0; bkg=0; sgnf=0;
 237     errsgn=0; errbkg=0; errsgnf=0;
 238     if (ipt == npt-1) fitter->WriteNtuple(pathout);
 239
 240     out<<endl;
 241     fitter->Reset(); //delete histogram set
 242
 243     init=kFALSE;
 244
 245   }
 246
 247   out.close();
 248 }