AliPhysics  master (3d17d9d)
AliHFInvMassFitter.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2019, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TH1F.h>
17 #include <TF1.h>
18 #include <TMath.h>
19 #include <TVirtualFitter.h>
20 #include <TDatabasePDG.h>
21 #include <TCanvas.h>
22 #include <TStyle.h>
23 #include <TPaveText.h>
24 
25 #include "AliVertexingHFUtils.h"
26 #include "AliHFInvMassFitter.h"
27 
29 ClassImp(AliHFInvMassFitter);
31 
47 
48 //__________________________________________________________________________
50  TNamed(),
51  fHistoInvMass(0x0),
52  fMinMass(0),
53  fMaxMass(5),
54  fTypeOfFit4Bkg(kExpo),
55  fPolDegreeBkg(4),
56  fCurPolDegreeBkg(-1),
57  fMassParticle(1.864),
58  fTypeOfFit4Sgn(kGaus),
59  fMass(1.865),
60  fMassErr(0.),
61  fSigmaSgn(0.012),
62  fSigmaSgnErr(0.),
63  fSigmaSgn2Gaus(0.012),
64  fFixedMean(kFALSE),
65  fFixedSigma(kFALSE),
66  fBoundSigma(kFALSE),
67  fSigmaVar(0.012),
68  fParSig(0.1),
69  fFixedSigma2Gaus(kFALSE),
70  fFixedRawYield(-1.),
71  fFrac2Gaus(0.2),
72  fFixedFrac2Gaus(kFALSE),
73  fRatio2GausSigma(0.),
74  fFixedRatio2GausSigma(kFALSE),
75  fNParsSig(3),
76  fNParsBkg(2),
77  fOnlySideBands(kFALSE),
78  fNSigma4SideBands(4.),
79  fCheckSignalCountsAfterFirstFit(kTRUE),
80  fFitOption("L,E"),
81  fRawYield(0.),
82  fRawYieldErr(0.),
83  fSigFunc(0x0),
84  fBkgFuncSb(0x0),
85  fBkgFunc(0x0),
86  fBkgFuncRefit(0x0),
87  fReflections(kFALSE),
88  fNParsRfl(0),
89  fRflOverSig(0),
90  fFixRflOverSig(kFALSE),
91  fHistoTemplRfl(0x0),
92  fSmoothRfl(kFALSE),
93  fRawYieldHelp(0),
94  fRflFunc(0x0),
95  fBkRFunc(0x0),
96  fSecondPeak(kFALSE),
97  fNParsSec(0),
98  fSecMass(-999.),
99  fSecWidth(9999.),
100  fFixSecMass(kFALSE),
101  fFixSecWidth(kFALSE),
102  fSecFunc(0x0),
103  fTotFunc(0x0)
104 {
106 }
107 
108 //__________________________________________________________________________
109 AliHFInvMassFitter::AliHFInvMassFitter(const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t fittypeb, Int_t fittypes):
110  TNamed(),
111  fHistoInvMass(0x0),
112  fMinMass(minvalue),
113  fMaxMass(maxvalue),
114  fTypeOfFit4Bkg(fittypeb),
115  fPolDegreeBkg(4),
116  fCurPolDegreeBkg(-1),
117  fMassParticle(1.864),
118  fTypeOfFit4Sgn(fittypes),
119  fMass(1.865),
120  fMassErr(0.),
121  fSigmaSgn(0.012),
122  fSigmaSgnErr(0.),
123  fSigmaSgn2Gaus(0.012),
124  fFixedMean(kFALSE),
125  fFixedSigma(kFALSE),
126  fBoundSigma(kFALSE),
127  fSigmaVar(0.012),
128  fParSig(0.1),
129  fFixedSigma2Gaus(kFALSE),
130  fFixedRawYield(-1.),
131  fFrac2Gaus(0.2),
132  fFixedFrac2Gaus(kFALSE),
133  fRatio2GausSigma(0.),
134  fFixedRatio2GausSigma(kFALSE),
135  fNParsSig(3),
136  fNParsBkg(2),
137  fOnlySideBands(kFALSE),
138  fNSigma4SideBands(4.),
140  fFitOption("L,E"),
141  fRawYield(0.),
142  fRawYieldErr(0.),
143  fSigFunc(0x0),
144  fBkgFuncSb(0x0),
145  fBkgFunc(0x0),
146  fBkgFuncRefit(0x0),
147  fReflections(kFALSE),
148  fNParsRfl(0),
149  fRflOverSig(0),
150  fFixRflOverSig(kFALSE),
151  fHistoTemplRfl(0x0),
152  fSmoothRfl(kFALSE),
153  fRawYieldHelp(0),
154  fRflFunc(0x0),
155  fBkRFunc(0x0),
156  fSecondPeak(kFALSE),
157  fNParsSec(0),
158  fSecMass(-999.),
159  fSecWidth(9999.),
160  fFixSecMass(kFALSE),
161  fFixSecWidth(kFALSE),
162  fSecFunc(0x0),
163  fTotFunc(0x0)
164 {
166  fHistoInvMass=(TH1F*)histoToFit->Clone("fHistoInvMass");
167  fHistoInvMass->SetDirectory(0);
169 }
170 //_________________________________________________________________________
172 
174 
175  delete fSigFunc;
176  delete fBkgFunc;
177  delete fBkgFuncSb;
178  delete fBkgFuncRefit;
179  delete fHistoInvMass;
180  delete fRflFunc;
181  delete fBkRFunc;
182  delete fSecFunc;
183  delete fTotFunc;
184  delete fHistoTemplRfl;
185 
186 }
187 //__________________________________________________________________________
191 
192  switch (fTypeOfFit4Bkg) {
193  case 0:
194  fNParsBkg=2;
195  break;
196  case 1:
197  fNParsBkg=2;
198  break;
199  case 2:
200  fNParsBkg=3;
201  break;
202  case 3:
203  fNParsBkg=1;
204  break;
205  case 4:
206  fNParsBkg=2;
207  break;
208  case 5:
209  fNParsBkg=2;
210  break;
211  case 6:
213  break;
214  default:
215  AliError("Error in computing fNParsBkg: check fTypeOfFit4Bkg");
216  break;
217  }
218 
219  switch (fTypeOfFit4Sgn) {
220  case 0:
221  fNParsSig=3;
222  break;
223  case 1:
224  fNParsSig=5;
225  break;
226  case 2:
227  fNParsSig=5;
228  break;
229  default:
230  AliError("Error in computing fNParsSig: check fTypeOfFit4Sgn");
231  break;
232  }
233 
234  if(fReflections) fNParsRfl=1;
235  else fNParsRfl=0;
236 
237  if(fSecondPeak) fNParsSec=3;
238  else fNParsSec=0;
239 
240 }
241 //__________________________________________________________________________
247 
248  TVirtualFitter::SetDefaultFitter("Minuit");
249 
250  Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");
251 
252  fOnlySideBands = kTRUE;
253  fBkgFuncSb = CreateBackgroundFitFunction("funcbkgsb",integralHisto);
254  Int_t status=-1;
255  printf("\n--- First fit with only background on the side bands - Exclusion region = %.2f sigma ---\n",fNSigma4SideBands);
256  if(fTypeOfFit4Bkg==6){
258  fHistoInvMass->GetListOfFunctions()->Add(fBkgFuncSb);
259  fHistoInvMass->GetFunction(fBkgFuncSb->GetName())->SetBit(1<<9,kTRUE);
260  status=0;
261  }
262  }
263  else status=fHistoInvMass->Fit("funcbkgsb",Form("R,%s,+,0",fFitOption.Data()));
264  fBkgFuncSb->SetLineColor(kGray+1);
265  if (status != 0){
266  printf(" ---> Failed first fit with only background, minuit status = %d\n",status);
267  return 0;
268  }
269 
270  fOnlySideBands = kFALSE;
271  if(!fBkgFunc){
272  fBkgFunc = CreateBackgroundFitFunction("funcbkg",integralHisto);
273  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFunc->SetParameter(ipar,fBkgFuncSb->GetParameter(ipar));
274  }
275  fBkgFunc->SetLineColor(kGray+1);
276 
277  printf("\n--- Estimate signal counts in the peak region ---\n");
278  Double_t estimSignal=CheckForSignal(fMass,fSigmaSgn);
279  Bool_t doFinalFit=kTRUE;
280  if(fCheckSignalCountsAfterFirstFit && estimSignal<0.){
281  if(draw) DrawFit();
282  estimSignal=0.;
283  doFinalFit=kFALSE;
284  printf("Abandon fit: no signal counts after first fit\n");
285  }
286 
287  fRawYieldHelp=estimSignal; // needed for reflection normalization
288  if(!fBkgFuncRefit){
289  fBkgFuncRefit = CreateBackgroundFitFunction("funcbkgrefit",integralHisto);
290  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFuncRefit->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
291  }
292  fBkgFuncRefit->SetLineColor(2);
293  fSigFunc = CreateSignalFitFunction("fsigfit",estimSignal);
294  if(fSecondPeak){
295  printf(" ---> Final fit includes a second inv. mass peak\n");
297  fSecFunc = CreateSecondPeakFunction("fsecpeak",estimSec);
298  }
299  if(fReflections){
300  printf(" ---> Final fit includes reflections\n");
301  fRflFunc = CreateReflectionFunction("freflect");
303  }
304  fTotFunc = CreateTotalFitFunction("funcmass");
305 
306  if(doFinalFit){
307  printf("\n--- Final fit with signal+background on the full range ---\n");
308  status=fHistoInvMass->Fit("funcmass",Form("R,%s,+,0",fFitOption.Data()));
309  if (status != 0){
310  printf(" ---> Failed fit with signal+background, minuit status = %d\n",status);
311  return 0;
312  }
313  }
314 
315  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
316  fBkgFuncRefit->SetParameter(ipar,fTotFunc->GetParameter(ipar));
317  fBkgFuncRefit->SetParError(ipar,fTotFunc->GetParError(ipar));
318  }
319  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
320  fSigFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg));
321  fSigFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg));
322  }
323  if(fSecondPeak){
324  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
325  fSecFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig));
326  fSecFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig));
327  }
328  fSecFunc->SetLineColor(kMagenta+1);
329  fSecFunc->SetLineStyle(3);
330  }
331  if(fReflections){
332  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
333  fRflFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
334  fRflFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
335  fBkRFunc->SetParameter(ipar+fNParsBkg,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
336  fBkRFunc->SetParError(ipar+fNParsBkg,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
337  }
338  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
339  fBkRFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar));
340  fBkRFunc->SetParError(ipar,fTotFunc->GetParError(ipar));
341  }
342  fRflFunc->SetLineColor(kGreen+1);
343  fBkRFunc->SetLineColor(kRed+1);
344  fBkRFunc->SetLineStyle(7);
345  }
346  fMass=fSigFunc->GetParameter(1);
347  fMassErr=fSigFunc->GetParError(1);
348  fSigmaSgn=fSigFunc->GetParameter(2);
349  fSigmaSgnErr=fSigFunc->GetParError(2);
350  fTotFunc->SetLineColor(4);
351  fRawYield=fTotFunc->GetParameter(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
352  fRawYieldErr=fTotFunc->GetParError(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
354  if(draw) DrawFit();
355  if(doFinalFit) return 1;
356  else return 2;
357 }
358 //______________________________________________________________________________
362 
363  TCanvas* c0=new TCanvas("c0");
364  DrawHere(c0);
365 }
366 //______________________________________________________________________________
367 void AliHFInvMassFitter::DrawHere(TVirtualPad* c, Double_t nsigma,Int_t writeFitInfo){
370 
371  gStyle->SetOptStat(0);
372  gStyle->SetCanvasColor(0);
373  gStyle->SetFrameFillColor(0);
374  c->cd();
375  fHistoInvMass->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
376  fHistoInvMass->SetMarkerStyle(20);
377  fHistoInvMass->SetMinimum(0.);
378  fHistoInvMass->Draw("PE");
379  if(fBkgFunc) fBkgFunc->Draw("same");
380  if(fBkgFuncRefit) fBkgFuncRefit->Draw("same");
381  if(fBkRFunc) fBkRFunc->Draw("same");
382  if(fRflFunc) fRflFunc->Draw("same");
383  if(fSecFunc) fSecFunc->Draw("same");
384  if(fTotFunc) fTotFunc->Draw("same");
385 
386  if(writeFitInfo > 0){
387  TPaveText *pinfos=new TPaveText(0.12,0.65,0.47,0.89,"NDC");
388  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
389  pinfos->SetBorderSize(0);
390  pinfos->SetFillStyle(0);
391  pinfom->SetBorderSize(0);
392  pinfom->SetFillStyle(0);
393  if(fTotFunc){
394  pinfom->SetTextColor(kBlue);
395  for(Int_t ipar=1; ipar<fNParsSig; ipar++){
396  pinfom->AddText(Form("%s = %.3f #pm %.3f",fTotFunc->GetParName(ipar+fNParsBkg),fTotFunc->GetParameter(ipar+fNParsBkg),fTotFunc->GetParError(ipar+fNParsBkg)));
397  }
398  if(writeFitInfo>=1) pinfom->Draw();
399 
400  Double_t bkg,errbkg;
401  Background(nsigma,bkg,errbkg);
402  Double_t signif,errsignif;
403  Significance(nsigma,signif,errsignif);
404 
405  pinfos->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
406  pinfos->AddText(Form("B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
407  pinfos->AddText(Form("S/B (%.0f#sigma) = %.4g ",nsigma,fRawYield/bkg));
408  if(fRflFunc) pinfos->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fRflFunc->GetParameter(0),fRflFunc->GetParError(0)));
409  pinfos->AddText(Form("Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
410  if(writeFitInfo>=2) pinfos->Draw();
411  }
412  }
413  c->Update();
414  return;
415 }
416 //______________________________________________________________________________
417 void AliHFInvMassFitter::DrawHistoMinusFit(TVirtualPad* c, Int_t writeFitInfo){
420 
421  gStyle->SetOptStat(0);
422  gStyle->SetCanvasColor(0);
423  gStyle->SetFrameFillColor(0);
424  c->cd();
425  TH1F* hInvMassMinusFit=(TH1F*)fHistoInvMass->Clone(Form("%sMinusFit",fHistoInvMass->GetName()));
426  hInvMassMinusFit->Reset("ICEM");
427  Double_t ymin=0;
428  Double_t ymax=1;
429  for(Int_t jst=1;jst<=fHistoInvMass->GetNbinsX();jst++){
430  Double_t integBkg=fBkgFuncRefit->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst);
431  hInvMassMinusFit->SetBinContent(jst,fHistoInvMass->GetBinContent(jst)-integBkg);
432  hInvMassMinusFit->SetBinError(jst,fHistoInvMass->GetBinError(jst));
433  if(hInvMassMinusFit->GetBinCenter(jst)>fMinMass && hInvMassMinusFit->GetBinCenter(jst)<fMaxMass){
434  if(hInvMassMinusFit->GetBinContent(jst)+hInvMassMinusFit->GetBinError(jst)>ymax) ymax=hInvMassMinusFit->GetBinContent(jst)+hInvMassMinusFit->GetBinError(jst);
435  if(hInvMassMinusFit->GetBinContent(jst)-hInvMassMinusFit->GetBinError(jst)<ymin) ymin=hInvMassMinusFit->GetBinContent(jst)-hInvMassMinusFit->GetBinError(jst);
436  }
437  }
438  hInvMassMinusFit->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
439  hInvMassMinusFit->SetMarkerStyle(20);
440  hInvMassMinusFit->SetMinimum(ymin-0.08*TMath::Abs(ymin));
441  hInvMassMinusFit->SetMaximum(ymax+0.08*TMath::Abs(ymax));
442  hInvMassMinusFit->DrawCopy();
443  if(fSigFunc) fSigFunc->Draw("same");
444  if(fRflFunc) fRflFunc->Draw("same");
445  if(fSecFunc) fSecFunc->Draw("same");
446  if(writeFitInfo > 0){
447  TPaveText *pinfom=new TPaveText(0.15,0.6,0.48,.87,"NDC");
448  pinfom->SetBorderSize(0);
449  pinfom->SetFillStyle(0);
450  if(fTotFunc){
451  for(Int_t ipar=1; ipar<fNParsSig; ipar++){
452  pinfom->AddText(Form("%s = %.3f #pm %.3f",fTotFunc->GetParName(ipar+fNParsBkg),fTotFunc->GetParameter(ipar+fNParsBkg),fTotFunc->GetParError(ipar+fNParsBkg)));
453  }
454  pinfom->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
455  if(fRflFunc) pinfom->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fRflFunc->GetParameter(0),fRflFunc->GetParError(0)));
456  pinfom->Draw();
457  }
458  }
459  c->Update();
460  return;
461 }
462 //______________________________________________________________________________
466 
467  Double_t minForSig=mean-4.*sigma;
468  Double_t maxForSig=mean+4.*sigma;
469  Int_t binForMinSig=fHistoInvMass->FindBin(minForSig);
470  Int_t binForMaxSig=fHistoInvMass->FindBin(maxForSig);
471  Double_t sum=0.;
472  Double_t sumback=0.;
473  fBkgFunc->Print();
474  for(Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
475  sum+=fHistoInvMass->GetBinContent(ibin);
476  sumback+=fBkgFunc->Eval(fHistoInvMass->GetBinCenter(ibin));
477  }
478  Double_t diffUnderPeak=(sum-sumback);
479  printf(" ---> IntegralUnderHisto=%f IntegralUnderBkgFunc=%f EstimatedSignal=%f\n",sum,sumback,diffUnderPeak);
480  if(diffUnderPeak/TMath::Sqrt(sum)<1.){
481  printf(" ---> (Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal\n",diffUnderPeak/TMath::Sqrt(sum));
482  return -1;
483  }
484  return diffUnderPeak*fHistoInvMass->GetBinWidth(1);
485 }
486 
487 //______________________________________________________________________________
491 
493  TF1* funcbkg = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Bkg,fMinMass,fMaxMass,fNParsBkg,"AliHFInvMassFitter","FitFunction4Bkg");
494  switch (fTypeOfFit4Bkg) {
495  case 0: //gaus+expo
496  funcbkg->SetParNames("BkgInt","Slope");
497  funcbkg->SetParameters(integral,-2.);
498  break;
499  case 1:
500  funcbkg->SetParNames("BkgInt","Slope");
501  funcbkg->SetParameters(integral,-100.);
502  break;
503  case 2:
504  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
505  funcbkg->SetParameters(integral,-10.,5);
506  break;
507  case 3:
508  funcbkg->SetParNames("Const");
509  funcbkg->SetParameter(0,0.);
510  funcbkg->FixParameter(0,0.);
511  break;
512  case 4:
513  funcbkg->SetParNames("BkgInt","Coef1");
514  funcbkg->SetParameters(integral,0.5);
515  break;
516  case 5:
517  funcbkg->SetParNames("Coef1","Coef2");
518  funcbkg->SetParameters(-10.,5.);
519  break;
520  case 6:
521  for(Int_t j=0;j<fNParsBkg; j++){
522  funcbkg->SetParName(j,Form("Coef%d",j));
523  funcbkg->SetParameter(j,0);
524  }
525  funcbkg->SetParameter(0,integral);
526  break;
527  default:
528  AliError(Form("Wrong choice of fTypeOfFit4Bkg (%d)",fTypeOfFit4Bkg));
529  delete funcbkg;
530  return 0x0;
531  break;
532  }
533  // if(fFixToHistoIntegral) funcbkg->FixParameter(0,integral);
534  funcbkg->SetLineColor(kBlue+3);
535  return funcbkg;
536 }
537 //______________________________________________________________________________
541 
542  TF1* funcsec = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4SecPeak,fMinMass,fMaxMass,3,"AliHFInvMassFitter","FitFunction4SecPeak");
543  funcsec->SetParameter(0,integsig);
544  funcsec->SetParameter(1,fSecMass);
545  if(fFixSecMass) funcsec->FixParameter(1,fSecMass);
546  funcsec->SetParameter(2,fSecWidth);
547  if(fFixSecWidth) funcsec->FixParameter(2,fSecWidth);
548  funcsec->SetParNames("SecPeakInt","SecMean","SecSigma");
549  return funcsec;
550 }
551 //______________________________________________________________________________
555  TF1* funcrfl = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Refl,fMinMass,fMaxMass,1,"AliHFInvMassFitter","FitFunction4Refl");
556  funcrfl->SetParameter(0,fRflOverSig);
557  funcrfl->SetParLimits(0,0.,1.);
558  if(fFixRflOverSig) funcrfl->FixParameter(0,fRflOverSig);
559  funcrfl->SetParNames("ReflOverS");
560  return funcrfl;
561 }
562 //______________________________________________________________________________
567  Int_t totParams=fNParsBkg+fNParsRfl;
568  TF1* fbr=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4BkgAndRefl,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4BkgAndRefl");
569  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
570  fbr->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
571  fbr->SetParName(ipar,fBkgFunc->GetParName(ipar));
572  }
573  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
574  fbr->SetParameter(ipar+fNParsBkg,fRflFunc->GetParameter(ipar));
575  fbr->SetParName(ipar+fNParsBkg,fRflFunc->GetParName(ipar));
576  // par limits not set because this function is not used for fitting but only for drawing
577  }
578  return fbr;
579 }
580 //______________________________________________________________________________
584 
586  TF1* funcsig = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Sgn,fMinMass,fMaxMass,fNParsSig,"AliHFInvMassFitter","FitFunction4Sgn");
587  if(fTypeOfFit4Sgn==kGaus){
588  funcsig->SetParameter(0,integsig);
589  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
590  funcsig->SetParameter(1,fMass);
591  if(fFixedMean) funcsig->FixParameter(1,fMass);
592  funcsig->SetParameter(2,fSigmaSgn);
593  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
594  if(fBoundSigma) funcsig->SetParLimits(2,fSigmaVar*(1-fParSig), fSigmaVar*(1+fParSig));
595  funcsig->SetParNames("SgnInt","Mean","Sigma");
596  }
597  if(fTypeOfFit4Sgn==k2Gaus){
598  funcsig->SetParameter(0,integsig);
599  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
600  funcsig->SetParameter(1,fMass);
601  if(fFixedMean) funcsig->FixParameter(1,fMass);
602  funcsig->SetParameter(2,fSigmaSgn);
603  funcsig->SetParLimits(2,0.004,0.05);
604  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
605  funcsig->SetParameter(3,fFrac2Gaus);
606  if(fFixedFrac2Gaus) funcsig->FixParameter(3,fFrac2Gaus);
607  else funcsig->SetParLimits(3,0.,1.);
608  funcsig->SetParameter(4,fSigmaSgn2Gaus);
609  if(fFixedSigma2Gaus) funcsig->FixParameter(4,fSigmaSgn2Gaus);
610  else funcsig->SetParLimits(4,0.004,0.05);
611  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","Sigma2");
612  }
614  funcsig->SetParameter(0,integsig);
615  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
616  funcsig->SetParameter(1,fMass);
617  if(fFixedMean) funcsig->FixParameter(1,fMass);
618  funcsig->SetParameter(2,fSigmaSgn);
619  funcsig->SetParLimits(2,0.004,0.05);
620  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
621  funcsig->SetParameter(3,fFrac2Gaus);
622  if(fFixedFrac2Gaus) funcsig->FixParameter(3,fFrac2Gaus);
623  else funcsig->SetParLimits(3,0.,1.);
624  funcsig->SetParameter(4,fSigmaSgn2Gaus);
625  if(fFixedRatio2GausSigma) funcsig->FixParameter(4,fRatio2GausSigma);
626  else funcsig->SetParLimits(4,0.,20.);
627  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","RatioSigma12");
628  }
629  return funcsig;
630 }
631 
632 //______________________________________________________________________________
636 
639  TF1* ftot=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4Mass");
640  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
641  ftot->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
642  ftot->SetParName(ipar,fBkgFunc->GetParName(ipar));
643  // Double_t parmin,parmax;
644  // fBkgFunc->GetParLimits(ipar,parmin,parmax);
645  // ftot->SetParLimits(ipar,parmin,parmax);
646  }
647  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
648  ftot->SetParameter(ipar+fNParsBkg,fSigFunc->GetParameter(ipar));
649  ftot->SetParName(ipar+fNParsBkg,fSigFunc->GetParName(ipar));
650  Double_t parmin,parmax;
651  fSigFunc->GetParLimits(ipar,parmin,parmax);
652  ftot->SetParLimits(ipar+fNParsBkg,parmin,parmax);
653  }
654  if(fSecondPeak && fSecFunc){
655  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
656  ftot->SetParameter(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParameter(ipar));
657  ftot->SetParName(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParName(ipar));
658  Double_t parmin,parmax;
659  fSecFunc->GetParLimits(ipar,parmin,parmax);
660  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig,parmin,parmax);
661  }
662  }
663  if(fReflections && fRflFunc){
664  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
665  ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParameter(ipar));
666  ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParName(ipar));
667  Double_t parmin,parmax;
668  fRflFunc->GetParLimits(ipar,parmin,parmax);
669  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec,parmin,parmax);
670  }
671  }
672  return ftot;
673 }
674 //__________________________________________________________________________
678 
680  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
681  TF1::RejectPoint();
682  return 0;
683  }
684  if(fOnlySideBands && fSecondPeak && TMath::Abs(x[0]-fSecMass) < (fNSigma4SideBands*fSecWidth)){
685  TF1::RejectPoint();
686  return 0;
687  }
688  Double_t total=0;
689 
690  switch (fTypeOfFit4Bkg){
691  case 0:
692  //exponential
693  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
694  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
695  //as integralTot- integralGaus (=par [2])
696  //Par:
697  // * [0] = integralBkg;
698  // * [1] = B;
699  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
700  total = par[0]*par[1]/(TMath::Exp(par[1]*fMaxMass)-TMath::Exp(par[1]*fMinMass))*TMath::Exp(par[1]*x[0]);
701  // AliInfo("Background function set to: exponential");
702  break;
703  case 1:
704  //linear
705  //y=a+b*x -> integral = a(max-min)+1/2*b*(max^2-min^2) -> a = (integral-1/2*b*(max^2-min^2))/(max-min)=integral/(max-min)-1/2*b*(max+min)
706  // * [0] = integralBkg;
707  // * [1] = b;
708  total= par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
709  // AliInfo("Background function set to: linear");
710  break;
711  case 2:
712  //parabola
713  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
714  //+ 1/3*c*(max^3-min^3) ->
715  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
716  // * [0] = integralBkg;
717  // * [1] = b;
718  // * [2] = c;
719  total = par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass))+par[2]*(x[0]*x[0]-1/3.*(fMaxMass*fMaxMass*fMaxMass-fMinMass*fMinMass*fMinMass)/(fMaxMass-fMinMass));
720  // AliInfo("Background function set to: polynomial");
721  break;
722  case 3:
723  total=par[0];
724  break;
725  case 4:
726  //power function
727  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
728  //
729  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
730  // * [0] = integralBkg;
731  // * [1] = b;
732  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
733  {
734  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
735 
736  total = par[0]*(par[1]+1.)/(TMath::Power(fMaxMass-mpi,par[1]+1.)-TMath::Power(fMinMass-mpi,par[1]+1.))*TMath::Power(x[0]-mpi,par[1]);
737  // AliInfo("Background function set to: powerlaw");
738  }
739  break;
740  case 5:
741  //power function wit exponential
742  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
743  {
744  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
745 
746  total = par[0]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[1]*(x[0]-mpi));
747  // AliInfo("Background function set to: wit exponential");
748  }
749  break;
750  case 6:
751  // the following comment must be removed
752  // // pol 3, following convention for pol 2
753  // //y=a+b*x+c*x**2+d*x**3 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
754  // //+ 1/3*c*(max^3-min^3) + 1/4 d * (max^4-min^4) ->
755  // //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3) - 1/4 d * (max^4-min^4) )/(max-min)
756  // // * [0] = integralBkg;
757  // // * [1] = b;
758  // // * [2] = c;
759  // // * [3] = d;
760  {
761  total=par[0];
762  for(Int_t it=1;it<=fPolDegreeBkg;it++){
763  total+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
764  }
765  }
766  break;
767  }
768  return total;
769 }
770 //_________________________________________________________________________
774 
775  // AliInfo("Signal function set to: Gaussian");
776  Double_t sigval=0;
777  Double_t g1=0;
778  Double_t g2=0;
779  switch (fTypeOfFit4Sgn){
780  case 0:
781  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
782  //Par:
783  // * [0] = integralSgn
784  // * [1] = mean
785  // * [2] = sigma
786  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
787  sigval=par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
788  break;
789  case 1:
790  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
791  //Par:
792  // * [0] = integralSgn
793  // * [1] = mean
794  // * [2] = sigma1
795  // * [3] = 2nd gaussian ratio
796  // * [4] = deltaSigma
797  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
798  g1=(1.-par[3])/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
799  g2=par[3]/TMath::Sqrt(2.*TMath::Pi())/par[4]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[4]/par[4]);
800  sigval=par[0]*(g1+g2);
801  break;
802  case 2:
803  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
804  //Par:
805  // * [0] = integralSgn
806  // * [1] = mean
807  // * [2] = sigma1
808  // * [3] = 2nd gaussian ratio
809  // * [4] = ratio sigma12
810  g1=(1.-par[3])/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
811  g2=par[3]/TMath::Sqrt(2.*TMath::Pi())/(par[4]*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./(par[4]*par[2])/(par[4]*par[2]));
812  sigval=par[0]*(g1+g2);
813  break;
814  }
815  fRawYieldHelp=par[0]/fHistoInvMass->GetBinWidth(1);
816  return sigval;
817 }
818 //_________________________________________________________________________
822  if(!fHistoTemplRfl) return 0;
823 
824  Int_t bin =fHistoTemplRfl->FindBin(x[0]);
825  Double_t value=fHistoTemplRfl->GetBinContent(bin);
826  Int_t binmin=fHistoTemplRfl->FindBin(fMinMass*1.00001);
827  Int_t binmax=fHistoTemplRfl->FindBin(fMaxMass*0.99999);
828  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
829  if(TMath::Abs(value)<1.e-14&&fSmoothRfl){// very rough, assume a constant trend, much better would be a pol1 or pol2 over a broader range
830  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
831  value/=3.;
832  }
833  return par[0]*value/norm*fRawYieldHelp*fHistoInvMass->GetBinWidth(1);
834 }
835 //_________________________________________________________________________
839  Double_t bkg=FitFunction4Bkg(x,par);
840  Double_t refl=0;
841  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg]);
842  return bkg+refl;
843 }
844 //_________________________________________________________________________
848 
849  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
850  //Par:
851  // * [0] = integralSgn
852  // * [1] = mean
853  // * [2] = sigma
854  Double_t secgaval=par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
855  return secgaval;
856 }
857 //_________________________________________________________________________
861 
862  Double_t bkg=FitFunction4Bkg(x,par);
863  Double_t sig=FitFunction4Sgn(x,&par[fNParsBkg]);
864  Double_t sec=0.;
865  if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNParsBkg+fNParsSig]);
866  Double_t refl=0;
867  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg+fNParsSig+fNParsSec]);
868  return bkg+sig+sec+refl;
869 }
870 
871 //_________________________________________________________________________
872 void AliHFInvMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
875 
876  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
877  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
878  Signal(minMass,maxMass,signal,errsignal);
879  return;
880 }
881 
882 //_________________________________________________________________________
883 void AliHFInvMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
886 
887  signal=fSigFunc->Integral(min, max)/(Double_t)fHistoInvMass->GetBinWidth(1);
888  errsignal=(fRawYieldErr/fRawYield)*signal;/*assume relative error is the same as for total integral*/
889  return;
890 }
891 
892 //___________________________________________________________________________
893 void AliHFInvMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
896 
897  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
898  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
899  Background(minMass,maxMass,background,errbackground);
900  return;
901 
902 }
903 //___________________________________________________________________________
904 void AliHFInvMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
907 
908  TF1 *funcbkg=0x0;
909  if(fBkgFuncRefit) funcbkg=fBkgFuncRefit;
910  else if(fBkgFunc) funcbkg=fBkgFunc;
911  if(!funcbkg){
912  AliError("Bkg function not found!");
913  return;
914  }
915 
916  Double_t intB=funcbkg->GetParameter(0);
917  Double_t intBerr=funcbkg->GetParError(0);
918 
919  //relative error evaluation: from histo
920 
922  Int_t rightBand=fHistoInvMass->FindBin(fMass+fNSigma4SideBands*fSigmaSgn);
923  intB=fHistoInvMass->Integral(1,leftBand)+fHistoInvMass->Integral(rightBand,fHistoInvMass->GetNbinsX());
924  Double_t sum2=0;
925  for(Int_t i=1;i<=leftBand;i++){
926  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
927  }
928  for(Int_t i=rightBand; i<=fHistoInvMass->GetNbinsX();i++){
929  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
930  }
931 
932  intBerr=TMath::Sqrt(sum2);
933 
934  background=funcbkg->Integral(min,max)/(Double_t)fHistoInvMass->GetBinWidth(1);
935  errbackground=intBerr/intB*background;
936 
937  return;
938 
939 }
940 //__________________________________________________________________________
941 
942 void AliHFInvMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
945 
946  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
947  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
948  Significance(minMass, maxMass, significance, errsignificance);
949 
950  return;
951 }
952 
953 //__________________________________________________________________________
954 
955 void AliHFInvMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
958 
959  Double_t background,errbackground;
960  Background(min,max,background,errbackground);
961 
962  if (fRawYield+background <= 0.){
963  significance=-1;
964  errsignificance=0;
965  return;
966  }
967 
968  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldErr,background,errbackground,significance,errsignificance);
969 
970  return;
971 }
972 //________________________________________________________________________
976 
977  TH1F *hCp=(TH1F*)fHistoInvMass->Clone("htemp");
978  Double_t estimatecent=0.5*(hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn))+hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn)));// just a first rough estimate
979  Double_t estimateslope=(hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn))-hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn)))/(7*fSigmaSgn);// first rough estimate
980 
982  TF1 *funcbkg,*funcPrev=0x0;
984  funcbkg = new TF1(Form("temp%d",fCurPolDegreeBkg),this,&AliHFInvMassFitter::BackFitFuncPolHelper,fMinMass,fMaxMass,fCurPolDegreeBkg+1,"AliHFInvMassFitter","BackFitFuncPolHelper");
985  if(funcPrev){
986  for(Int_t j=0;j<fCurPolDegreeBkg;j++){// now is +1 degree w.r.t. previous fit funct
987  funcbkg->SetParameter(j,funcPrev->GetParameter(j));
988  }
989  delete funcPrev;
990  }
991  else{
992  funcbkg->SetParameter(0,estimatecent);
993  funcbkg->SetParameter(1,estimateslope);
994  }
995  printf(" ---> Pre-fit of background with pol degree %d ---\n",fCurPolDegreeBkg);
996  hCp->Fit(funcbkg,"REMN","");
997  funcPrev=(TF1*)funcbkg->Clone("ftemp");
998  delete funcbkg;
1000  }
1001 
1002  for(Int_t j=0;j<=fPolDegreeBkg;j++){
1003  fback->SetParameter(j,funcPrev->GetParameter(j));
1004  fback->SetParError(j,funcPrev->GetParError(j));
1005  }
1006  printf(" ---> Final background fit with pol degree %d ---\n",fPolDegreeBkg);
1007  hCp->Fit(fback,"REMN","");// THIS IS JUST TO SET NOT ONLY THE PARAMETERS BUT ALSO chi2, etc...
1008 
1009 
1010  // The following lines might be useful for debugging
1011  // TCanvas *cDebug=new TCanvas();
1012  // cDebug->cd();
1013  // hCp->Draw();
1014  // TString strout=Form("Test%d.root",(Int_t)fhistoInvMass->GetBinContent(fhistoInvMass->FindBin(fMass)));
1015  // cDebug->Print(strout.Data());
1016  // delete cDebug;
1017 
1018  delete funcPrev;
1019  delete hCp;
1020  return kTRUE;
1021 
1022 }
1023 // _______________________________________________________________________
1027 
1028  Double_t maxDeltaM = fNSigma4SideBands*fSigmaSgn;
1029  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
1030  TF1::RejectPoint();
1031  return 0;
1032  }
1033  Double_t back=par[0];
1034  for(Int_t it=1;it<=fCurPolDegreeBkg;it++){
1035  back+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
1036  }
1037  return back;
1038 }
1039 // _______________________________________________________________________
1048 
1049  if(!h){
1050  fReflections=kFALSE;
1051  return 0x0;
1052  }
1053  fHistoTemplRfl=(TH1F*)h->Clone("hTemplRfl");
1054  opt.ToLower();
1055  fReflections=kTRUE;
1056  printf("\n--- Reflection templates from simulation ---\n");
1057  if(opt.Contains("templ")){
1058  printf(" ---> Reflection contribution using directly the histogram from simulation\n");
1059  return fHistoTemplRfl;
1060  }
1061 
1062  TF1 *f=0x0;
1063  Bool_t isPoissErr=kTRUE;
1064  Double_t xMinForFit=h->GetBinLowEdge(1);
1065  Double_t xMaxForFit=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
1066  if(minRange>=0 && maxRange>=0){
1067  xMinForFit=TMath::Max(minRange,h->GetBinLowEdge(1));
1068  xMaxForFit=TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1069  }
1070  if(opt.EqualTo("1gaus") || opt.EqualTo("singlegaus")){
1071  printf(" ---> Reflection contribution from single-Gaussian fit to histogram from simulation\n");
1072  f=new TF1("mygaus","gaus",xMinForFit,xMaxForFit);
1073  f->SetParameter(0,h->GetMaximum());
1074  // f->SetParLimits(0,0,100.*h->Integral());
1075  f->SetParameter(1,1.865);
1076  f->SetParameter(2,0.050);
1077  fHistoTemplRfl->Fit(f,"REM0","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1078  }
1079  else if(opt.EqualTo("2gaus") || opt.EqualTo("doublegaus")){
1080  printf(" ---> Reflection contribution from double-Gaussian fit to histogram from simulation\n");
1081  f=new TF1("my2gaus","[0]*([3]/( TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[3])/( TMath::Sqrt(2.*TMath::Pi())*[5])*TMath::Exp(-(x-[4])*(x-[4])/(2.*[5]*[5])))",xMinForFit,xMaxForFit);
1082  f->SetParameter(0,h->GetMaximum());
1083  // f->SetParLimits(0,0,100.*h->Integral());
1084  f->SetParLimits(3,0.,1.);
1085  f->SetParameter(3,0.5);
1086  f->SetParameter(1,1.84);
1087  f->SetParameter(2,0.050);
1088  f->SetParameter(4,1.88);
1089  f->SetParameter(5,0.050);
1090  fHistoTemplRfl->Fit(f,"REM0","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1091  }
1092  else if(opt.EqualTo("pol3")){
1093  printf(" ---> Reflection contribution from pol3 fit to histogram from simulation\n");
1094  f=new TF1("mypol3","pol3",xMinForFit,xMaxForFit);
1095  f->SetParameter(0,h->GetMaximum());
1096  // f->SetParLimits(0,0,100.*h->Integral());
1097  // Hard to initialize the other parameters...
1098  fHistoTemplRfl->Fit(f,"REM0","");
1099  // Printf("We USED %d POINTS in the Fit",f->GetNumberFitPoints());
1100  }
1101  else if(opt.EqualTo("pol6")){
1102  printf(" ---> Reflection contribution from pol6 fit to histogram from simulation\n");
1103  f=new TF1("mypol6","pol6",xMinForFit,xMaxForFit);
1104  f->SetParameter(0,h->GetMaximum());
1105  // f->SetParLimits(0,0,100.*h->Integral());
1106  // Hard to initialize the other parameters...
1107  fHistoTemplRfl->Fit(f,"RLEMI0","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1108  }
1109  else{
1110  // no good option passed
1111  printf(" ---> Bad option for reflection configuration -> reflections will not be included in the fit\n");
1112  fReflections=kFALSE;
1113  delete fHistoTemplRfl;
1114  fHistoTemplRfl=0x0;
1115  return 0x0;
1116  }
1117 
1118  // Fill fHistoTemplRfl with values of fit function
1119  if(f){
1120  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1121  fHistoTemplRfl->SetBinContent(j,f->Integral(fHistoTemplRfl->GetBinLowEdge(j),fHistoTemplRfl->GetXaxis()->GetBinUpEdge(j))/fHistoTemplRfl->GetBinWidth(j));
1122  if(fHistoTemplRfl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
1123  }
1124  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1125  if(isPoissErr){
1126  if(fHistoTemplRfl->GetBinContent(j)>0) fHistoTemplRfl->SetBinError(j,TMath::Sqrt(fHistoTemplRfl->GetBinContent(j)));
1127  else fHistoTemplRfl->SetBinError(j,0);
1128  }
1129  else fHistoTemplRfl->SetBinError(j,0.001*fHistoTemplRfl->GetBinContent(j));
1130  }
1131  fReflections=kTRUE;
1132  return fHistoTemplRfl;
1133  }else{
1134  printf(" ---> Fit to MC template for reflection failed -> reflections will not be included in the fit\n");
1135  fReflections=kFALSE;
1136  delete fHistoTemplRfl;
1137  fHistoTemplRfl=0x0;
1138  return 0x0;
1139  }
1140  return 0x0;
1141 }
1142 // _______________________________________________________________________
1147  // else (default) mean of gaussian fit
1148 
1149  Double_t massVal=fMass;
1150  switch (pdgCode) {
1151  case 411:
1152  case 421:
1153  case 431:
1154  case 4122:
1155  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1156  break;
1157  case 413:
1158  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1159  massVal-=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1160  break;
1161  default:
1162  massVal=fMass;
1163  break;
1164  }
1165 
1166  Double_t minMass=massVal-nOfSigma*fSigmaSgn;
1167  Double_t maxMass=massVal+nOfSigma*fSigmaSgn;
1168  return GetRawYieldBinCounting(errRyBC,minMass,maxMass,option);
1169 }
1170 // _______________________________________________________________________
1176 
1177  Int_t minBinSum=fHistoInvMass->FindBin(minMass);
1178  Int_t maxBinSum=fHistoInvMass->FindBin(maxMass);
1179  if(minBinSum<1){
1180  printf("Left range for bin counting smaller than allowed by histogram axis, setting it to the lower edge of the first histo bin\n");
1181  minBinSum=1;
1182  }
1183  if(maxBinSum>fHistoInvMass->GetNbinsX()){
1184  printf("Right range for bin counting larger than allowed by histogram axis, setting it to the upper edge of the last histo bin\n");
1185  maxBinSum=fHistoInvMass->GetNbinsX();
1186  }
1187  Double_t cntSig=0.;
1188  Double_t cntErr=0.;
1189  errRyBC=0;
1190  TF1* fbackground=fBkgFunc;
1191  if(option==1) fbackground=fBkgFuncRefit;
1192  if(!fbackground) return 0.;
1193 
1194  for(Int_t jb=minBinSum; jb<=maxBinSum; jb++){
1195  Double_t cntTot=fHistoInvMass->GetBinContent(jb);
1196  Double_t cntBkg=fbackground->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1197  Double_t cntRefl=0;
1198  if(option==1 && fRflFunc) cntRefl=fRflFunc->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1199  //Double_t cntBkg=fbackground->Eval(fHistoInvMass->GetBinCenter(jb));
1200  Double_t cntSecPeak=0;
1201  if(option==1 && fSecondPeak && fSecFunc) cntSecPeak=fSecFunc->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1202  cntSig+=(cntTot-cntBkg-cntRefl-cntSecPeak);
1203  cntErr+=(fHistoInvMass->GetBinError(jb)*fHistoInvMass->GetBinError(jb));
1204  }
1205  errRyBC=TMath::Sqrt(cntErr);
1206  return cntSig;
1207 }
1208 
1209 
1210 // _______________________________________________________________________
1211 TH1F* AliHFInvMassFitter::GetResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange, Int_t option){
1212 
1214 
1215  Int_t binmi=fHistoInvMass->FindBin(fMinMass*1.001);
1216  Int_t binma=fHistoInvMass->FindBin(fMaxMass*0.9999);
1217  if(maxrange>minrange){
1218  binmi=fHistoInvMass->FindBin(minrange*1.001);
1219  binma=fHistoInvMass->FindBin(maxrange*0.9999);
1220  }
1221  if(hResidualTrend){
1222  //fHistoInvMass->Copy(hResidualTrend);
1223  hResidualTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1224  hResidualTrend->SetName(Form("%s_residualTrend",fHistoInvMass->GetName()));
1225  hResidualTrend->SetTitle(Form("%s (Residuals)",fHistoInvMass->GetTitle()));
1226  hResidualTrend->SetMarkerStyle(20);
1227  hResidualTrend->SetMarkerSize(1.0);
1228  hResidualTrend->Reset();
1229  }
1230  if(hPullsTrend){
1231  hPullsTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1232  hPullsTrend->Reset();
1233  hPullsTrend->SetName(Form("%s_pullTrend",fHistoInvMass->GetName()));
1234  hPullsTrend->SetTitle(Form("%s (Pulls)",fHistoInvMass->GetTitle()));
1235  hPullsTrend->SetMarkerStyle(20);
1236  hPullsTrend->SetMarkerSize(1.0);
1237  }
1238  if(hPulls){
1239  hPulls->SetName(Form("%s_pulls",fHistoInvMass->GetName()));
1240  hPulls->SetTitle(Form("%s ; Pulls",fHistoInvMass->GetTitle()));
1241  hPulls->SetBins(40,-10,10);
1242  hPulls->Reset();
1243  }
1244 
1245  Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
1246  TArrayD *arval=new TArrayD(binma-binmi+1);
1247  for(Int_t jst=1;jst<=fHistoInvMass->GetNbinsX();jst++){
1248  Double_t integFit=0;
1249  if(option==0) integFit=fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst));
1250  else{
1251  integFit=fBkgFuncRefit->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst));
1252  if(option==2) integFit+=fRflFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst));
1253  }
1254  res=fHistoInvMass->GetBinContent(jst)-integFit/fHistoInvMass->GetBinWidth(jst);
1255  if(jst>=binmi&&jst<=binma){
1256  arval->AddAt(res,jst-binmi);
1257  if(res<min)min=res;
1258  if(res>max)max=res;
1259  }
1260  // Printf("Res = %f from %f - %f",res,fHistoInvMass->GetBinContent(jst),fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst));
1261  if(hResidualTrend){
1262  hResidualTrend->SetBinContent(jst,res);
1263  hResidualTrend->SetBinError(jst,fHistoInvMass->GetBinError(jst));
1264  }
1265  if(hPulls){
1266  if(jst>=binmi&&jst<=binma)hPulls->Fill(res/fHistoInvMass->GetBinError(jst));
1267  }
1268  if(hPullsTrend){
1269  hPullsTrend->SetBinContent(jst,res/fHistoInvMass->GetBinError(jst));
1270  hPullsTrend->SetBinError(jst,0.0001);
1271  }
1272  }
1273  if(hResidualTrend){
1274  hResidualTrend->GetXaxis()->SetRange(binmi,binma);
1275  if(option!=0){
1276  TF1 *fgauss=new TF1("signalTermForRes","[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./[2]/[2])",fHistoInvMass->GetBinLowEdge(1),fHistoInvMass->GetBinLowEdge(fHistoInvMass->GetNbinsX()+1));
1277  fgauss->SetParameter(0,fRawYield*fHistoInvMass->GetBinWidth(1));
1278  fgauss->SetParameter(1,fMass);
1279  fgauss->SetParameter(2,fSigmaSgn);
1280  fgauss->SetLineColor(kBlue);
1281  hResidualTrend->GetListOfFunctions()->Add(fgauss);
1282  }
1283  }
1284  if(hPullsTrend){
1285  hPullsTrend->GetXaxis()->SetRange(binmi,binma);
1286  hPullsTrend->SetMinimum(-7);
1287  hPullsTrend->SetMaximum(+7);
1288  }
1289  if(TMath::Abs(min)>TMath::Abs(max))max=min;
1290 
1291  TH1F *hout=new TH1F(Form("%s_residuals",fHistoInvMass->GetName()),Form("%s ; residuals",fHistoInvMass->GetTitle()),25,-TMath::Abs(max)*1.5,TMath::Abs(max)*1.5);
1292  for(Int_t j=0;j<binma-binmi+1;j++){
1293  hout->Fill(arval->At(j));
1294  }
1295  hout->Sumw2();
1296  hout->Fit("gaus","LEM","",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
1297 
1298  if(hPulls){
1299  hPulls->Sumw2();
1300  hPulls->Fit("gaus","LEM","",-3,3);
1301  }
1302  delete arval;
1303  return hout;
1304 }
1305 // _______________________________________________________________________
1306 TH1F* AliHFInvMassFitter::GetOverBackgroundResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1308  return GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend,minrange,maxrange,1);
1309 }
1310 // _______________________________________________________________________
1311 TH1F* AliHFInvMassFitter::GetOverBackgroundPlusReflResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1313  return GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend,minrange,maxrange,2);
1314 }
1315 
1316 // _______________________________________________________________________
1320  if(fBkgFunc){
1321  printf("--- Background function in 1st step fit to side bands ---\n");
1322  fBkgFunc->Print();
1323  }
1324  if(fTotFunc){
1325  printf("--- Total fit function from 2nd step fit ---\n");
1326  fTotFunc->Print();
1327  }
1328  if(fBkgFuncRefit){
1329  printf("--- Background function in 2nd step fit ---\n");
1330  fBkgFuncRefit->Print();
1331  }
1332  if(fRflFunc){
1333  printf("--- Reflections ---\n");
1334  fRflFunc->Print();
1335  }
1336  if(fBkRFunc){
1337  printf("--- Background + reflections ---\n");
1338  fBkRFunc->Print();
1339  }
1340  if(fSecFunc){
1341  printf("--- Additional Gaussian peak ---\n");
1342  fSecFunc->Print();
1343  }
1344 }
Int_t fTypeOfFit4Sgn
pdg value of particle mass
Double_t fSigmaSgnErr
signal gaussian sigma
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1, Int_t option=0)
Bool_t fFixSecWidth
flag to fix the position of the 2nd peak
Double_t fRawYieldHelp
switch for smoothing of reflection template
const Double_t ymax
Definition: AddTaskCFDStar.C:7
double Double_t
Definition: External.C:58
Double_t fRawYield
L, LW or Chi2.
TF1 * CreateReflectionFunction(TString fname)
Int_t fCurPolDegreeBkg
degree of polynomial expansion for back fit (option 6 for back)
Bool_t PrepareHighPolFit(TF1 *fback)
Bool_t fOnlySideBands
fit parameters in background fit function
virtual void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
Bool_t fFixedRatio2GausSigma
initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar ...
Double_t FitFunction4SecPeak(Double_t *x, Double_t *par)
Double_t fMass
signal fit func
TF1 * CreateTotalFitFunction(TString fname)
Double_t fRawYieldErr
signal gaussian integral
Int_t MassFitter(Bool_t draw=kTRUE)
Int_t fNParsBkg
fit parameters in signal fit function
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
Int_t fNParsRfl
flag use/not use reflections
TCanvas * c
Definition: TestFitELoss.C:172
TF1 * CreateSecondPeakFunction(TString fname, Double_t integral)
Bool_t fFixedFrac2Gaus
initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar ...
TH1F * GetOverBackgroundPlusReflResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
TF1 * fTotFunc
fit function for second peak
Double_t fFixedRawYield
switch for fix Sigma of second gaussian in case of k2Gaus
Double_t FitFunction4Bkg(Double_t *x, Double_t *par)
Double_t fSigmaSgn
unc on signal gaussian mean value
Double_t fSecMass
fit parameters in 2nd peak fit function
TH1F * GetOverBackgroundResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
Bool_t fFixSecMass
width of the 2nd peak
TF1 * fSecFunc
flag to fix the width of the 2nd peak
Double_t * sigma
int Int_t
Definition: External.C:63
Bool_t fFixedSigma2Gaus
+/- range variation of bound Sigma of gaussian in %
Double_t fSecWidth
position of the 2nd peak
TString fFitOption
switch for check after first fit
TF1 * fBkRFunc
fit function for reflections
Double_t fFrac2Gaus
initialization for wa yield
Int_t fNParsSec
switch off/on second peak (for D+->KKpi in Ds)
TF1 * CreateSignalFitFunction(TString fname, Double_t integral)
Bool_t fFixedMean
signal second gaussian sigma in case of k2Gaus
Double_t fMassParticle
help variable
Double_t fSigmaVar
switch for bound Sigma of gaussian
Double_t FitFunction4Refl(Double_t *x, Double_t *par)
Double_t GetRawYieldBinCounting(Double_t &errRyBC, Double_t nSigma=3., Int_t option=0, Int_t pdgCode=0) const
TH1F * SetTemplateReflections(const TH1 *h, TString opt, Double_t minRange, Double_t maxRange)
Double_t nsigma
Int_t fTypeOfFit4Bkg
upper mass limit
TH1F * fHistoTemplRfl
switch for fix refl/signal
TF1 * fBkgFunc
background fit function (1st step, side bands only)
Double_t fMassErr
signal gaussian mean value
Double_t fRflOverSig
fit parameters in reflection fit function
Double_t fMinMass
histogram to fit
Bool_t fFixRflOverSig
reflection/signal
void DrawHistoMinusFit(TVirtualPad *c, Int_t writeFitInfo=1)
static void ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance, Double_t &errsignificance)
Significance calculator.
Double_t FitFunction4Mass(Double_t *x, Double_t *par)
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
Double_t BackFitFuncPolHelper(Double_t *x, Double_t *par)
Double_t fMaxMass
lower mass limit
TF1 * CreateBackgroundFitFunction(TString fname, Double_t integral)
TF1 * fRflFunc
internal variable for fit with reflections
Double_t fNSigma4SideBands
kTRUE = only side bands considered
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
Double_t minMass
Double_t fRatio2GausSigma
switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar ...
Double_t FitFunction4Sgn(Double_t *x, Double_t *par)
Int_t fPolDegreeBkg
background fit func
const Double_t ymin
Definition: AddTaskCFDStar.C:6
Double_t fSigmaSgn2Gaus
unc on signal gaussian sigma
Double_t maxMass
bool Bool_t
Definition: External.C:53
TF1 * fBkgFuncRefit
background fit function (1st step, extended in peak region)
Double_t CheckForSignal(Double_t mean, Double_t sigma)
TF1 * fSigFunc
err on signal gaussian integral
TF1 * fBkgFuncSb
Signal fit function.
Bool_t fCheckSignalCountsAfterFirstFit
number of sigmas to veto the signal peak
Bool_t draw[maxPtBins]
Double_t fParSig
value of bound Sigma of gaussian
Bool_t fSecondPeak
fit function for reflections
Bool_t fSmoothRfl
histogram with reflection template
Definition: External.C:196
TF1 * CreateBackgroundPlusReflectionFunction(TString fname)
Double_t FitFunction4BkgAndRefl(Double_t *x, Double_t *par)
Int_t fNParsSig
switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar ...
void DrawHere(TVirtualPad *c, Double_t nsigma=3, Int_t writeFitInfo=2)
Bool_t fBoundSigma
switch for fix Sigma of gaussian
Bool_t fReflections
background fit function (2nd step)
Bool_t fFixedSigma
switch for fix mean of gaussian