AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fitD0InvMass.C
Go to the documentation of this file.
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 
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 }
double Double_t
Definition: External.C:58
void fitD0(Int_t rebin=0, TString listname="coutputmassD0mycuts", Int_t nsigma=3, TString pathin="./", TString pathout="./", Int_t btype=2)
Definition: fitD0InvMass.C:22
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Double_t nsigma
TList * fitter
Definition: DrawAnaELoss.C:26
TFile * file
TList with histograms for a given trigger.
Int_t rebin
bool Bool_t
Definition: External.C:53
AliHFMassFitter for the fit of invariant mass distribution of charmed mesons.
TDirectoryFile * dir