AliPhysics  8d00e07 (8d00e07)
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  fFixedMean(kFALSE),
64  fFixedSigma(kFALSE),
65  fFixedRawYield(-1.),
66  fNParsSig(3),
67  fNParsBkg(2),
68  fOnlySideBands(kFALSE),
69  fNSigma4SideBands(4.),
70  fFitOption("L,E"),
71  fRawYield(0.),
72  fRawYieldErr(0.),
73  fSigFunc(0x0),
74  fBkgFuncSb(0x0),
75  fBkgFunc(0x0),
76  fBkgFuncRefit(0x0),
77  fReflections(kFALSE),
78  fNParsRfl(0),
79  fRflOverSig(0),
80  fFixRflOverSig(kFALSE),
81  fHistoTemplRfl(0x0),
82  fSmoothRfl(kFALSE),
83  fRawYieldHelp(0),
84  fRflFunc(0x0),
85  fBkRFunc(0x0),
86  fSecondPeak(kFALSE),
87  fNParsSec(0),
88  fSecMass(-999.),
89  fSecWidth(9999.),
90  fFixSecMass(kFALSE),
91  fFixSecWidth(kFALSE),
92  fSecFunc(0x0),
93  fTotFunc(0x0)
94 {
96 }
97 
98 //__________________________________________________________________________
99 AliHFInvMassFitter::AliHFInvMassFitter(const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t fittypeb, Int_t fittypes):
100  TNamed(),
101  fHistoInvMass(0x0),
102  fMinMass(minvalue),
103  fMaxMass(maxvalue),
104  fTypeOfFit4Bkg(fittypeb),
105  fPolDegreeBkg(4),
106  fCurPolDegreeBkg(-1),
107  fMassParticle(1.864),
108  fTypeOfFit4Sgn(fittypes),
109  fMass(1.865),
110  fMassErr(0.),
111  fSigmaSgn(0.012),
112  fSigmaSgnErr(0.),
113  fFixedMean(kFALSE),
114  fFixedSigma(kFALSE),
115  fFixedRawYield(-1.),
116  fNParsSig(3),
117  fNParsBkg(2),
118  fOnlySideBands(kFALSE),
119  fNSigma4SideBands(4.),
120  fFitOption("L,E"),
121  fRawYield(0.),
122  fRawYieldErr(0.),
123  fSigFunc(0x0),
124  fBkgFuncSb(0x0),
125  fBkgFunc(0x0),
126  fBkgFuncRefit(0x0),
127  fReflections(kFALSE),
128  fNParsRfl(0),
129  fRflOverSig(0),
130  fFixRflOverSig(kFALSE),
131  fHistoTemplRfl(0x0),
132  fSmoothRfl(kFALSE),
133  fRawYieldHelp(0),
134  fRflFunc(0x0),
135  fBkRFunc(0x0),
136  fSecondPeak(kFALSE),
137  fNParsSec(0),
138  fSecMass(-999.),
139  fSecWidth(9999.),
140  fFixSecMass(kFALSE),
141  fFixSecWidth(kFALSE),
142  fSecFunc(0x0),
143  fTotFunc(0x0)
144 {
146  fHistoInvMass=(TH1F*)histoToFit->Clone("fHistoInvMass");
147  fHistoInvMass->SetDirectory(0);
149 }
150 //_________________________________________________________________________
152 
154 
155  delete fSigFunc;
156  delete fBkgFunc;
157  delete fBkgFuncSb;
158  delete fBkgFuncRefit;
159  delete fHistoInvMass;
160  delete fRflFunc;
161  delete fBkRFunc;
162  delete fSecFunc;
163  delete fTotFunc;
164  delete fHistoTemplRfl;
165 
166 }
167 //__________________________________________________________________________
171 
172  switch (fTypeOfFit4Bkg) {
173  case 0:
174  fNParsBkg=2;
175  break;
176  case 1:
177  fNParsBkg=2;
178  break;
179  case 2:
180  fNParsBkg=3;
181  break;
182  case 3:
183  fNParsBkg=1;
184  break;
185  case 4:
186  fNParsBkg=2;
187  break;
188  case 5:
189  fNParsBkg=2;
190  break;
191  case 6:
193  break;
194  default:
195  AliError("Error in computing fNParsBkg: check fTypeOfFit4Bkg");
196  break;
197  }
198 
199  switch (fTypeOfFit4Sgn) {
200  case 0:
201  fNParsSig=3;
202  break;
203  case 1:
204  fNParsSig=5;
205  break;
206  default:
207  AliError("Error in computing fNParsSig: check fTypeOfFit4Sgn");
208  break;
209  }
210 
211  if(fReflections) fNParsRfl=1;
212  else fNParsRfl=0;
213 
214  if(fSecondPeak) fNParsSec=3;
215  else fNParsSec=0;
216 
217 }
218 //__________________________________________________________________________
224 
225  TVirtualFitter::SetDefaultFitter("Minuit");
226 
227  Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");
228 
229  fOnlySideBands = kTRUE;
230  fBkgFuncSb = CreateBackgroundFitFunction("funcbkgsb",integralHisto);
231  Int_t status=-1;
232  printf("\n--- First fit with only background on the side bands - Exclusion region = %.2f sigma ---\n",fNSigma4SideBands);
233  if(fTypeOfFit4Bkg==6){
235  fHistoInvMass->GetListOfFunctions()->Add(fBkgFuncSb);
236  fHistoInvMass->GetFunction(fBkgFuncSb->GetName())->SetBit(1<<9,kTRUE);
237  status=0;
238  }
239  }
240  else status=fHistoInvMass->Fit("funcbkgsb",Form("R,%s,+,0",fFitOption.Data()));
241  fBkgFuncSb->SetLineColor(kGray+1);
242  if (status != 0){
243  printf(" ---> Failed first fit with only background, minuit status = %d\n",status);
244  return 0;
245  }
246 
247  fOnlySideBands = kFALSE;
248  if(!fBkgFunc){
249  fBkgFunc = CreateBackgroundFitFunction("funcbkg",integralHisto);
250  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFunc->SetParameter(ipar,fBkgFuncSb->GetParameter(ipar));
251  }
252  fBkgFunc->SetLineColor(kGray+1);
253 
254  printf("\n--- Estimate signal counts in the peak region ---\n");
255  Double_t estimSignal=CheckForSignal(fMass,fSigmaSgn);
256  Bool_t doFinalFit=kTRUE;
257  if(estimSignal<0.){
258  if(draw) DrawFit();
259  estimSignal=0.;
260  doFinalFit=kFALSE;
261  }
262 
263  fRawYieldHelp=estimSignal; // needed for reflection normalization
264  if(!fBkgFuncRefit){
265  fBkgFuncRefit = CreateBackgroundFitFunction("funcbkgrefit",integralHisto);
266  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFuncRefit->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
267  }
268  fBkgFuncRefit->SetLineColor(2);
269  fSigFunc = CreateSignalFitFunction("fsigfit",estimSignal);
270  if(fSecondPeak){
271  printf(" ---> Final fit includes a second inv. mass peak\n");
273  fSecFunc = CreateSecondPeakFunction("fsecpeak",estimSec);
274  }
275  if(fReflections){
276  printf(" ---> Final fit includes reflections\n");
277  fRflFunc = CreateReflectionFunction("freflect");
279  }
280  fTotFunc = CreateTotalFitFunction("funcmass");
281 
282  if(doFinalFit){
283  printf("\n--- Final fit with signal+background on the full range ---\n");
284  status=fHistoInvMass->Fit("funcmass",Form("R,%s,+,0",fFitOption.Data()));
285  if (status != 0){
286  printf(" ---> Failed fit with signal+background, minuit status = %d\n",status);
287  return 0;
288  }
289  }
290 
291  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
292  fBkgFuncRefit->SetParameter(ipar,fTotFunc->GetParameter(ipar));
293  fBkgFuncRefit->SetParError(ipar,fTotFunc->GetParError(ipar));
294  }
295  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
296  fSigFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg));
297  fSigFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg));
298  }
299  if(fSecondPeak){
300  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
301  fSecFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig));
302  fSecFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig));
303  }
304  fSecFunc->SetLineColor(kMagenta+1);
305  fSecFunc->SetLineStyle(3);
306  }
307  if(fReflections){
308  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
309  fRflFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
310  fRflFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
311  fBkRFunc->SetParameter(ipar+fNParsBkg,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
312  fBkRFunc->SetParError(ipar+fNParsBkg,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
313  }
314  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
315  fBkRFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar));
316  fBkRFunc->SetParError(ipar,fTotFunc->GetParError(ipar));
317  }
318  fRflFunc->SetLineColor(kGreen+1);
319  fBkRFunc->SetLineColor(kRed+1);
320  fBkRFunc->SetLineStyle(7);
321  }
322  fMass=fSigFunc->GetParameter(1);
323  fMassErr=fSigFunc->GetParError(1);
324  fSigmaSgn=fSigFunc->GetParameter(2);
325  fSigmaSgnErr=fSigFunc->GetParError(2);
326  fTotFunc->SetLineColor(4);
327  fRawYield=fTotFunc->GetParameter(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
328  fRawYieldErr=fTotFunc->GetParError(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
330  if(draw) DrawFit();
331  if(doFinalFit) return 1;
332  else return 2;
333 }
334 //______________________________________________________________________________
338 
339  TCanvas* c0=new TCanvas("c0");
340  DrawHere(c0);
341 }
342 //______________________________________________________________________________
343 void AliHFInvMassFitter::DrawHere(TVirtualPad* c, Double_t nsigma,Int_t writeFitInfo){
346 
347  gStyle->SetOptStat(0);
348  gStyle->SetCanvasColor(0);
349  gStyle->SetFrameFillColor(0);
350  c->cd();
351  fHistoInvMass->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
352  fHistoInvMass->SetMarkerStyle(20);
353  fHistoInvMass->SetMinimum(0.);
354  fHistoInvMass->Draw("PE");
355  if(fBkgFunc) fBkgFunc->Draw("same");
356  if(fBkgFuncRefit) fBkgFuncRefit->Draw("same");
357  if(fBkRFunc) fBkRFunc->Draw("same");
358  if(fRflFunc) fRflFunc->Draw("same");
359  if(fSecFunc) fSecFunc->Draw("same");
360  if(fTotFunc) fTotFunc->Draw("same");
361 
362  if(writeFitInfo > 0){
363  TPaveText *pinfos=new TPaveText(0.12,0.65,0.47,0.89,"NDC");
364  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
365  pinfos->SetBorderSize(0);
366  pinfos->SetFillStyle(0);
367  pinfom->SetBorderSize(0);
368  pinfom->SetFillStyle(0);
369  if(fTotFunc){
370  pinfom->SetTextColor(kBlue);
371  for(Int_t ipar=1; ipar<fNParsSig; ipar++){
372  pinfom->AddText(Form("%s = %.3f #pm %.3f",fTotFunc->GetParName(ipar+fNParsBkg),fTotFunc->GetParameter(ipar+fNParsBkg),fTotFunc->GetParError(ipar+fNParsBkg)));
373  }
374  pinfom->Draw();
375 
376  Double_t bkg,errbkg;
377  Background(nsigma,bkg,errbkg);
378  Double_t signif,errsignif;
379  Significance(nsigma,signif,errsignif);
380 
381  pinfos->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
382  pinfos->AddText(Form("B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
383  pinfos->AddText(Form("S/B = (%.0f#sigma) %.4f ",nsigma,fRawYield/bkg));
384  if(fRflFunc) pinfos->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fRflFunc->GetParameter(0),fRflFunc->GetParError(0)));
385  pinfos->AddText(Form("Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
386  pinfos->Draw();
387  }
388  }
389  c->Update();
390  return;
391 }
392 //______________________________________________________________________________
396 
397  Double_t minForSig=mean-4.*sigma;
398  Double_t maxForSig=mean+4.*sigma;
399  Int_t binForMinSig=fHistoInvMass->FindBin(minForSig);
400  Int_t binForMaxSig=fHistoInvMass->FindBin(maxForSig);
401  Double_t sum=0.;
402  Double_t sumback=0.;
403  fBkgFunc->Print();
404  for(Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
405  sum+=fHistoInvMass->GetBinContent(ibin);
406  sumback+=fBkgFunc->Eval(fHistoInvMass->GetBinCenter(ibin));
407  }
408  Double_t diffUnderPeak=(sum-sumback);
409  printf(" ---> IntegralUnderFitFunc=%f IntegralUnderHisto=%f EstimatedSignal=%f\n",sum,sumback,diffUnderPeak);
410  if(diffUnderPeak/TMath::Sqrt(sum)<1.){
411  printf(" ---> (Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal/\n",diffUnderPeak/TMath::Sqrt(sum));
412  return -1;
413  }
414  return diffUnderPeak*fHistoInvMass->GetBinWidth(1);
415 }
416 
417 //______________________________________________________________________________
421 
423  TF1* funcbkg = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Bkg,fMinMass,fMaxMass,fNParsBkg,"AliHFInvMassFitter","FitFunction4Bkg");
424  switch (fTypeOfFit4Bkg) {
425  case 0: //gaus+expo
426  funcbkg->SetParNames("BkgInt","Slope");
427  funcbkg->SetParameters(integral,-2.);
428  break;
429  case 1:
430  funcbkg->SetParNames("BkgInt","Slope");
431  funcbkg->SetParameters(integral,-100.);
432  break;
433  case 2:
434  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
435  funcbkg->SetParameters(integral,-10.,5);
436  break;
437  case 3:
438  funcbkg->SetParNames("Const");
439  funcbkg->SetParameter(0,0.);
440  funcbkg->FixParameter(0,0.);
441  break;
442  case 4:
443  funcbkg->SetParNames("BkgInt","Coef1");
444  funcbkg->SetParameters(integral,0.5);
445  break;
446  case 5:
447  funcbkg->SetParNames("Coef1","Coef2");
448  funcbkg->SetParameters(-10.,5.);
449  break;
450  case 6:
451  for(Int_t j=0;j<fNParsBkg; j++){
452  funcbkg->SetParName(j,Form("Coef%d",j));
453  funcbkg->SetParameter(j,0);
454  }
455  funcbkg->SetParameter(0,integral);
456  break;
457  default:
458  AliError(Form("Wrong choice of fTypeOfFit4Bkg (%d)",fTypeOfFit4Bkg));
459  delete funcbkg;
460  return 0x0;
461  break;
462  }
463  // if(fFixToHistoIntegral) funcbkg->FixParameter(0,integral);
464  funcbkg->SetLineColor(kBlue+3);
465  return funcbkg;
466 }
467 //______________________________________________________________________________
471 
472  TF1* funcsec = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4SecPeak,fMinMass,fMaxMass,3,"AliHFInvMassFitter","FitFunction4SecPeak");
473  funcsec->SetParameter(0,integsig);
474  funcsec->SetParameter(1,fSecMass);
475  if(fFixSecMass) funcsec->FixParameter(1,fSecMass);
476  funcsec->SetParameter(2,fSecWidth);
477  if(fFixSecWidth) funcsec->FixParameter(2,fSecWidth);
478  funcsec->SetParNames("SecPeakInt","SecMean","SecSigma");
479  return funcsec;
480 }
481 //______________________________________________________________________________
485  TF1* funcrfl = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Refl,fMinMass,fMaxMass,1,"AliHFInvMassFitter","FitFunction4Refl");
486  funcrfl->SetParameter(0,fRflOverSig);
487  funcrfl->SetParLimits(0,0.,1.);
488  if(fFixRflOverSig) funcrfl->FixParameter(0,fRflOverSig);
489  funcrfl->SetParNames("ReflOverS");
490  return funcrfl;
491 }
492 //______________________________________________________________________________
497  Int_t totParams=fNParsBkg+fNParsRfl;
498  TF1* fbr=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4BkgAndRefl,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4BkgAndRefl");
499  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
500  fbr->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
501  fbr->SetParName(ipar,fBkgFunc->GetParName(ipar));
502  }
503  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
504  fbr->SetParameter(ipar+fNParsBkg,fRflFunc->GetParameter(ipar));
505  fbr->SetParName(ipar+fNParsBkg,fRflFunc->GetParName(ipar));
506  // par limits not set because this function is not used for fitting but only for drawing
507  }
508  return fbr;
509 }
510 //______________________________________________________________________________
514 
516  TF1* funcsig = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Sgn,fMinMass,fMaxMass,fNParsSig,"AliHFInvMassFitter","FitFunction4Sgn");
517  if(fTypeOfFit4Sgn==kGaus){
518  funcsig->SetParameter(0,integsig);
519  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
520  funcsig->SetParameter(1,fMass);
521  if(fFixedMean) funcsig->FixParameter(1,fMass);
522  funcsig->SetParameter(2,fSigmaSgn);
523  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
524  funcsig->SetParNames("SgnInt","Mean","Sigma");
525  }
526  if(fTypeOfFit4Sgn==k2Gaus){
527  funcsig->SetParameter(0,integsig);
528  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
529  funcsig->SetParameter(1,fMass);
530  if(fFixedMean) funcsig->FixParameter(1,fMass);
531  funcsig->SetParameter(2,fSigmaSgn);
532  funcsig->SetParLimits(2,0.004,0.05);
533  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
534  funcsig->SetParameter(3,0.2);
535  funcsig->SetParLimits(3,0.,1.);
536  funcsig->SetParameter(4,fSigmaSgn);
537  funcsig->SetParLimits(4,0.004,0.05);
538  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","Sigma2");
539  }
540  return funcsig;
541 }
542 
543 //______________________________________________________________________________
547 
550  TF1* ftot=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4Mass");
551  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
552  ftot->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
553  ftot->SetParName(ipar,fBkgFunc->GetParName(ipar));
554  // Double_t parmin,parmax;
555  // fBkgFunc->GetParLimits(ipar,parmin,parmax);
556  // ftot->SetParLimits(ipar,parmin,parmax);
557  }
558  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
559  ftot->SetParameter(ipar+fNParsBkg,fSigFunc->GetParameter(ipar));
560  ftot->SetParName(ipar+fNParsBkg,fSigFunc->GetParName(ipar));
561  Double_t parmin,parmax;
562  fSigFunc->GetParLimits(ipar,parmin,parmax);
563  ftot->SetParLimits(ipar+fNParsBkg,parmin,parmax);
564  }
565  if(fSecondPeak && fSecFunc){
566  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
567  ftot->SetParameter(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParameter(ipar));
568  ftot->SetParName(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParName(ipar));
569  Double_t parmin,parmax;
570  fSecFunc->GetParLimits(ipar,parmin,parmax);
571  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig,parmin,parmax);
572  }
573  }
574  if(fReflections && fRflFunc){
575  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
576  ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParameter(ipar));
577  ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParName(ipar));
578  Double_t parmin,parmax;
579  fRflFunc->GetParLimits(ipar,parmin,parmax);
580  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec,parmin,parmax);
581  }
582  }
583  return ftot;
584 }
585 //__________________________________________________________________________
589 
591  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
592  TF1::RejectPoint();
593  return 0;
594  }
595  if(fOnlySideBands && fSecondPeak && TMath::Abs(x[0]-fSecMass) < (fNSigma4SideBands*fSecWidth)){
596  TF1::RejectPoint();
597  return 0;
598  }
599  Double_t total=0;
600 
601  switch (fTypeOfFit4Bkg){
602  case 0:
603  //exponential
604  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
605  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
606  //as integralTot- integralGaus (=par [2])
607  //Par:
608  // * [0] = integralBkg;
609  // * [1] = B;
610  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
611  total = par[0]*par[1]/(TMath::Exp(par[1]*fMaxMass)-TMath::Exp(par[1]*fMinMass))*TMath::Exp(par[1]*x[0]);
612  // AliInfo("Background function set to: exponential");
613  break;
614  case 1:
615  //linear
616  //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)
617  // * [0] = integralBkg;
618  // * [1] = b;
619  total= par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
620  // AliInfo("Background function set to: linear");
621  break;
622  case 2:
623  //parabola
624  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
625  //+ 1/3*c*(max^3-min^3) ->
626  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
627  // * [0] = integralBkg;
628  // * [1] = b;
629  // * [2] = c;
630  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));
631  // AliInfo("Background function set to: polynomial");
632  break;
633  case 3:
634  total=par[0];
635  break;
636  case 4:
637  //power function
638  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
639  //
640  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
641  // * [0] = integralBkg;
642  // * [1] = b;
643  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
644  {
645  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
646 
647  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]);
648  // AliInfo("Background function set to: powerlaw");
649  }
650  break;
651  case 5:
652  //power function wit exponential
653  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
654  {
655  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
656 
657  total = par[0]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[1]*(x[0]-mpi));
658  // AliInfo("Background function set to: wit exponential");
659  }
660  break;
661  case 6:
662  // the following comment must be removed
663  // // pol 3, following convention for pol 2
664  // //y=a+b*x+c*x**2+d*x**3 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
665  // //+ 1/3*c*(max^3-min^3) + 1/4 d * (max^4-min^4) ->
666  // //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)
667  // // * [0] = integralBkg;
668  // // * [1] = b;
669  // // * [2] = c;
670  // // * [3] = d;
671  {
672  total=par[0];
673  for(Int_t it=1;it<=fPolDegreeBkg;it++){
674  total+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
675  }
676  }
677  break;
678  }
679  return total;
680 }
681 //_________________________________________________________________________
685 
686  // AliInfo("Signal function set to: Gaussian");
687  Double_t sigval=0;
688  switch (fTypeOfFit4Sgn){
689  case 0:
690  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
691  //Par:
692  // * [0] = integralSgn
693  // * [1] = mean
694  // * [2] = sigma
695  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
696  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]);
697  break;
698  case 1:
699  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
700  //Par:
701  // * [0] = integralSgn
702  // * [1] = mean
703  // * [2] = sigma1
704  // * [3] = 2nd gaussian ratio
705  // * [4] = deltaSigma
706  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
707  Double_t 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]);
708  Double_t 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]);
709  sigval=par[0]*(g1+g2);
710  break;
711  }
712  fRawYieldHelp=par[0]/fHistoInvMass->GetBinWidth(1);
713  return sigval;
714 }
715 //_________________________________________________________________________
719  if(!fHistoTemplRfl) return 0;
720 
721  Int_t bin =fHistoTemplRfl->FindBin(x[0]);
722  Double_t value=fHistoTemplRfl->GetBinContent(bin);
723  Int_t binmin=fHistoTemplRfl->FindBin(fMinMass*1.00001);
724  Int_t binmax=fHistoTemplRfl->FindBin(fMaxMass*0.99999);
725  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
726  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
727  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
728  value/=3.;
729  }
730  return par[0]*value/norm*fRawYieldHelp*fHistoInvMass->GetBinWidth(1);
731 }
732 //_________________________________________________________________________
736  Double_t bkg=FitFunction4Bkg(x,par);
737  Double_t refl=0;
738  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg]);
739  return bkg+refl;
740 }
741 //_________________________________________________________________________
745 
746  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
747  //Par:
748  // * [0] = integralSgn
749  // * [1] = mean
750  // * [2] = sigma
751  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]);
752  return secgaval;
753 }
754 //_________________________________________________________________________
758 
759  Double_t bkg=FitFunction4Bkg(x,par);
760  Double_t sig=FitFunction4Sgn(x,&par[fNParsBkg]);
761  Double_t sec=0.;
762  if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNParsBkg+fNParsSig]);
763  Double_t refl=0;
764  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg+fNParsSig+fNParsSec]);
765  return bkg+sig+sec+refl;
766 }
767 
768 //_________________________________________________________________________
769 void AliHFInvMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
772 
773  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
774  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
775  Signal(minMass,maxMass,signal,errsignal);
776  return;
777 }
778 
779 //_________________________________________________________________________
780 void AliHFInvMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
783 
784  signal=fSigFunc->Integral(min, max)/(Double_t)fHistoInvMass->GetBinWidth(1);
785  errsignal=(fRawYieldErr/fRawYield)*signal;/*assume relative error is the same as for total integral*/
786  return;
787 }
788 
789 //___________________________________________________________________________
790 void AliHFInvMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
793 
794  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
795  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
796  Background(minMass,maxMass,background,errbackground);
797  return;
798 
799 }
800 //___________________________________________________________________________
801 void AliHFInvMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
804 
805  TF1 *funcbkg=0x0;
806  if(fBkgFuncRefit) funcbkg=fBkgFuncRefit;
807  else if(fBkgFunc) funcbkg=fBkgFunc;
808  if(!funcbkg){
809  AliError("Bkg function not found!");
810  return;
811  }
812 
813  Double_t intB=funcbkg->GetParameter(0);
814  Double_t intBerr=funcbkg->GetParError(0);
815 
816  //relative error evaluation: from histo
817 
819  Int_t rightBand=fHistoInvMass->FindBin(fMass+fNSigma4SideBands*fSigmaSgn);
820  intB=fHistoInvMass->Integral(1,leftBand)+fHistoInvMass->Integral(rightBand,fHistoInvMass->GetNbinsX());
821  Double_t sum2=0;
822  for(Int_t i=1;i<=leftBand;i++){
823  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
824  }
825  for(Int_t i=rightBand; i<=fHistoInvMass->GetNbinsX();i++){
826  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
827  }
828 
829  intBerr=TMath::Sqrt(sum2);
830 
831  background=funcbkg->Integral(min,max)/(Double_t)fHistoInvMass->GetBinWidth(1);
832  errbackground=intBerr/intB*background;
833 
834  return;
835 
836 }
837 //__________________________________________________________________________
838 
839 void AliHFInvMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
842 
843  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
844  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
845  Significance(minMass, maxMass, significance, errsignificance);
846 
847  return;
848 }
849 
850 //__________________________________________________________________________
851 
852 void AliHFInvMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
855 
856  Double_t background,errbackground;
857  Background(min,max,background,errbackground);
858 
859  if (fRawYield+background <= 0.){
860  significance=-1;
861  errsignificance=0;
862  return;
863  }
864 
865  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldErr,background,errbackground,significance,errsignificance);
866 
867  return;
868 }
869 //________________________________________________________________________
873 
874  TH1F *hCp=(TH1F*)fHistoInvMass->Clone("htemp");
875  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
876  Double_t estimateslope=(hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn))-hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn)))/(7*fSigmaSgn);// first rough estimate
877 
879  TF1 *funcbkg,*funcPrev=0x0;
881  funcbkg = new TF1(Form("temp%d",fCurPolDegreeBkg),this,&AliHFInvMassFitter::BackFitFuncPolHelper,fMinMass,fMaxMass,fCurPolDegreeBkg+1,"AliHFInvMassFitter","BackFitFuncPolHelper");
882  if(funcPrev){
883  for(Int_t j=0;j<fCurPolDegreeBkg;j++){// now is +1 degree w.r.t. previous fit funct
884  funcbkg->SetParameter(j,funcPrev->GetParameter(j));
885  }
886  delete funcPrev;
887  }
888  else{
889  funcbkg->SetParameter(0,estimatecent);
890  funcbkg->SetParameter(1,estimateslope);
891  }
892  printf(" ---> Pre-fit of background with pol degree %d ---\n",fCurPolDegreeBkg);
893  hCp->Fit(funcbkg,"REMN","");
894  funcPrev=(TF1*)funcbkg->Clone("ftemp");
895  delete funcbkg;
897  }
898 
899  for(Int_t j=0;j<=fPolDegreeBkg;j++){
900  fback->SetParameter(j,funcPrev->GetParameter(j));
901  fback->SetParError(j,funcPrev->GetParError(j));
902  }
903  printf(" ---> Final background fit with pol degree %d ---\n",fPolDegreeBkg);
904  hCp->Fit(fback,"REMN","");// THIS IS JUST TO SET NOT ONLY THE PARAMETERS BUT ALSO chi2, etc...
905 
906 
907  // The following lines might be useful for debugging
908  // TCanvas *cDebug=new TCanvas();
909  // cDebug->cd();
910  // hCp->Draw();
911  // TString strout=Form("Test%d.root",(Int_t)fhistoInvMass->GetBinContent(fhistoInvMass->FindBin(fMass)));
912  // cDebug->Print(strout.Data());
913  // delete cDebug;
914 
915  delete funcPrev;
916  delete hCp;
917  return kTRUE;
918 
919 }
920 // _______________________________________________________________________
924 
926  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
927  TF1::RejectPoint();
928  return 0;
929  }
930  Double_t back=par[0];
931  for(Int_t it=1;it<=fCurPolDegreeBkg;it++){
932  back+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
933  }
934  return back;
935 }
936 // _______________________________________________________________________
945 
946  if(!h){
947  fReflections=kFALSE;
948  return 0x0;
949  }
950  fHistoTemplRfl=(TH1F*)h->Clone("hTemplRfl");
951  opt.ToLower();
952  fReflections=kTRUE;
953  printf("\n--- Reflection templates from simulation ---\n");
954  if(opt.Contains("templ")){
955  printf(" ---> Reflection contribution using directly the histogram from simulation\n");
956  return fHistoTemplRfl;
957  }
958 
959  TF1 *f=0x0;
960  Bool_t isPoissErr=kTRUE;
961  Double_t xMinForFit=h->GetBinLowEdge(1);
962  Double_t xMaxForFit=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
963  if(minRange>=0 && maxRange>=0){
964  xMinForFit=TMath::Max(minRange,h->GetBinLowEdge(1));
965  xMaxForFit=TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
966  }
967  if(opt.EqualTo("1gaus") || opt.EqualTo("singlegaus")){
968  printf(" ---> Reflection contribution from single-Gaussian fit to histogram from simulation\n");
969  f=new TF1("mygaus","gaus",xMinForFit,xMaxForFit);
970  f->SetParameter(0,h->GetMaximum());
971  // f->SetParLimits(0,0,100.*h->Integral());
972  f->SetParameter(1,1.865);
973  f->SetParameter(2,0.050);
974  fHistoTemplRfl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
975  }
976  else if(opt.EqualTo("2gaus") || opt.EqualTo("doublegaus")){
977  printf(" ---> Reflection contribution from double-Gaussian fit to histogram from simulation\n");
978  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);
979  f->SetParameter(0,h->GetMaximum());
980  // f->SetParLimits(0,0,100.*h->Integral());
981  f->SetParLimits(3,0.,1.);
982  f->SetParameter(3,0.5);
983  f->SetParameter(1,1.84);
984  f->SetParameter(2,0.050);
985  f->SetParameter(4,1.88);
986  f->SetParameter(5,0.050);
987  fHistoTemplRfl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
988  }
989  else if(opt.EqualTo("pol3")){
990  printf(" ---> Reflection contribution from pol3 fit to histogram from simulation\n");
991  f=new TF1("mypol3","pol3",xMinForFit,xMaxForFit);
992  f->SetParameter(0,h->GetMaximum());
993  // f->SetParLimits(0,0,100.*h->Integral());
994  // Hard to initialize the other parameters...
995  fHistoTemplRfl->Fit(f,"REM","");
996  // Printf("We USED %d POINTS in the Fit",f->GetNumberFitPoints());
997  }
998  else if(opt.EqualTo("pol6")){
999  printf(" ---> Reflection contribution from pol6 fit to histogram from simulation\n");
1000  f=new TF1("mypol6","pol6",xMinForFit,xMaxForFit);
1001  f->SetParameter(0,h->GetMaximum());
1002  // f->SetParLimits(0,0,100.*h->Integral());
1003  // Hard to initialize the other parameters...
1004  fHistoTemplRfl->Fit(f,"RLEMI","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1005  }
1006  else{
1007  // no good option passed
1008  printf(" ---> Bad option for reflection configuration -> reflections will not be included in the fit\n");
1009  fReflections=kFALSE;
1010  delete fHistoTemplRfl;
1011  fHistoTemplRfl=0x0;
1012  return 0x0;
1013  }
1014 
1015  // Fill fHistoTemplRfl with values of fit function
1016  if(f){
1017  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1018  fHistoTemplRfl->SetBinContent(j,f->Integral(fHistoTemplRfl->GetBinLowEdge(j),fHistoTemplRfl->GetXaxis()->GetBinUpEdge(j))/fHistoTemplRfl->GetBinWidth(j));
1019  if(fHistoTemplRfl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
1020  }
1021  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1022  if(isPoissErr){
1023  fHistoTemplRfl->SetBinError(j,TMath::Sqrt(fHistoTemplRfl->GetBinContent(j)));
1024  }
1025  else fHistoTemplRfl->SetBinError(j,0.001*fHistoTemplRfl->GetBinContent(j));
1026  }
1027  fReflections=kTRUE;
1028  return fHistoTemplRfl;
1029  }else{
1030  printf(" ---> Fit to MC template for reflection failed -> reflections will not be included in the fit\n");
1031  fReflections=kFALSE;
1032  delete fHistoTemplRfl;
1033  fHistoTemplRfl=0x0;
1034  return 0x0;
1035  }
1036  return 0x0;
1037 }
1038 // _______________________________________________________________________
1043  // else (default) mean of gaussian fit
1044 
1045  Double_t massVal=fMass;
1046  switch (pdgCode) {
1047  case 411:
1048  case 421:
1049  case 431:
1050  case 4122:
1051  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1052  break;
1053  case 413:
1054  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1055  massVal-=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1056  break;
1057  default:
1058  massVal=fMass;
1059  break;
1060  }
1061 
1062  Double_t minMass=massVal-nOfSigma*fSigmaSgn;
1063  Double_t maxMass=massVal+nOfSigma*fSigmaSgn;
1064  return GetRawYieldBinCounting(errRyBC,minMass,maxMass,option);
1065 }
1066 // _______________________________________________________________________
1072 
1073  Int_t minBinSum=fHistoInvMass->FindBin(minMass);
1074  Int_t maxBinSum=fHistoInvMass->FindBin(maxMass);
1075  if(minBinSum<1){
1076  printf("Left range for bin counting smaller than allowed by histogram axis, setting it to the lower edge of the first histo bin\n");
1077  minBinSum=1;
1078  }
1079  if(maxBinSum>fHistoInvMass->GetNbinsX()){
1080  printf("Right range for bin counting larger than allowed by histogram axis, setting it to the upper edge of the last histo bin\n");
1081  maxBinSum=fHistoInvMass->GetNbinsX();
1082  }
1083  Double_t cntSig=0.;
1084  Double_t cntErr=0.;
1085  errRyBC=0;
1086  TF1* fbackground=fBkgFunc;
1087  if(option==1) fbackground=fBkgFuncRefit;
1088  if(!fbackground) return 0.;
1089 
1090  for(Int_t jb=minBinSum; jb<=maxBinSum; jb++){
1091  Double_t cntTot=fHistoInvMass->GetBinContent(jb);
1092  Double_t cntBkg=fbackground->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1093  Double_t cntRefl=0;
1094  if(option==1 && fRflFunc) cntRefl=fRflFunc->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1095  //Double_t cntBkg=fbackground->Eval(fHistoInvMass->GetBinCenter(jb));
1096  cntSig+=(cntTot-cntBkg-cntRefl);
1097  cntErr+=(fHistoInvMass->GetBinError(jb)*fHistoInvMass->GetBinError(jb));
1098  }
1099  errRyBC=TMath::Sqrt(cntErr);
1100  return cntSig;
1101 }
1102 
1103 
1104 // _______________________________________________________________________
1105 TH1F* AliHFInvMassFitter::GetResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1106 
1108 
1109  Int_t binmi=fHistoInvMass->FindBin(fMinMass*1.001);
1110  Int_t binma=fHistoInvMass->FindBin(fMaxMass*0.9999);
1111  if(maxrange>minrange){
1112  binmi=fHistoInvMass->FindBin(minrange*1.001);
1113  binma=fHistoInvMass->FindBin(maxrange*0.9999);
1114  }
1115  if(hResidualTrend){
1116  //fHistoInvMass->Copy(hResidualTrend);
1117  hResidualTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1118  hResidualTrend->SetName(Form("%s_residualTrend",fHistoInvMass->GetName()));
1119  hResidualTrend->SetTitle(Form("%s (Residuals)",fHistoInvMass->GetTitle()));
1120  hResidualTrend->SetMarkerStyle(20);
1121  hResidualTrend->SetMarkerSize(1.0);
1122  hResidualTrend->Reset();
1123  }
1124  if(hPullsTrend){
1125  hPullsTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1126  hPullsTrend->Reset();
1127  hPullsTrend->SetName(Form("%s_pullTrend",fHistoInvMass->GetName()));
1128  hPullsTrend->SetTitle(Form("%s (Pulls)",fHistoInvMass->GetTitle()));
1129  hPullsTrend->SetMarkerStyle(20);
1130  hPullsTrend->SetMarkerSize(1.0);
1131  }
1132  if(hPulls){
1133  hPulls->SetName(Form("%s_pulls",fHistoInvMass->GetName()));
1134  hPulls->SetTitle(Form("%s ; Pulls",fHistoInvMass->GetTitle()));
1135  hPulls->SetBins(40,-10,10);
1136  hPulls->Reset();
1137  }
1138 
1139  Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
1140  TArrayD *arval=new TArrayD(binma-binmi+1);
1141  for(Int_t jst=1;jst<=fHistoInvMass->GetNbinsX();jst++){
1142 
1143  res=fHistoInvMass->GetBinContent(jst)-fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst);
1144  if(jst>=binmi&&jst<=binma){
1145  arval->AddAt(res,jst-binmi);
1146  if(res<min)min=res;
1147  if(res>max)max=res;
1148  }
1149  // Printf("Res = %f from %f - %f",res,fHistoInvMass->GetBinContent(jst),fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst));
1150  if(hResidualTrend){
1151  hResidualTrend->SetBinContent(jst,res);
1152  hResidualTrend->SetBinError(jst,fHistoInvMass->GetBinError(jst));
1153  }
1154  if(hPulls){
1155  if(jst>=binmi&&jst<=binma)hPulls->Fill(res/fHistoInvMass->GetBinError(jst));
1156  }
1157  if(hPullsTrend){
1158  hPullsTrend->SetBinContent(jst,res/fHistoInvMass->GetBinError(jst));
1159  hPullsTrend->SetBinError(jst,0.0001);
1160  }
1161  }
1162  if(hResidualTrend)hResidualTrend->GetXaxis()->SetRange(binmi,binma);
1163  if(hPullsTrend){
1164  hPullsTrend->GetXaxis()->SetRange(binmi,binma);
1165  hPullsTrend->SetMinimum(-7);
1166  hPullsTrend->SetMaximum(+7);
1167  }
1168  if(TMath::Abs(min)>TMath::Abs(max))max=min;
1169 
1170  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);
1171  for(Int_t j=0;j<binma-binmi+1;j++){
1172  hout->Fill(arval->At(j));
1173  }
1174  hout->Sumw2();
1175  hout->Fit("gaus","LEM","",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
1176 
1177  if(hPulls){
1178  hPulls->Sumw2();
1179  hPulls->Fit("gaus","LEM","",-3,3);
1180  }
1181  delete arval;
1182  return hout;
1183 }
1184 // _______________________________________________________________________
1188  if(fBkgFunc){
1189  printf("--- Background function in 1st step fit to side bands ---\n");
1190  fBkgFunc->Print();
1191  }
1192  if(fTotFunc){
1193  printf("--- Total fit function from 2nd step fit ---\n");
1194  fTotFunc->Print();
1195  }
1196  if(fBkgFuncRefit){
1197  printf("--- Background function in 2nd step fit ---\n");
1198  fBkgFuncRefit->Print();
1199  }
1200  if(fRflFunc){
1201  printf("--- Reflections ---\n");
1202  fRflFunc->Print();
1203  }
1204  if(fBkRFunc){
1205  printf("--- Background + reflections ---\n");
1206  fBkRFunc->Print();
1207  }
1208  if(fSecFunc){
1209  printf("--- Additional Gaussian peak ---\n");
1210  fSecFunc->Print();
1211  }
1212 }
Int_t fTypeOfFit4Sgn
pdg value of particle mass
Double_t fSigmaSgnErr
signal gaussian sigma
Bool_t fFixSecWidth
flag to fix the position of the 2nd peak
Double_t fRawYieldHelp
switch for smoothing of reflection template
double Double_t
Definition: External.C:58
Double_t fRawYield
L, LW or Chi2.
void DrawHere(TVirtualPad *c, Double_t nsigma=3, Int_t writeFitInfo=1)
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
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)
TF1 * fTotFunc
fit function for second peak
Double_t fFixedRawYield
switch for fix Sigma of gaussian
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
Bool_t fFixSecMass
width of the 2nd peak
TF1 * fSecFunc
flag to fix the width of the 2nd peak
Double_t * sigma
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
int Int_t
Definition: External.C:63
Double_t fSecWidth
position of the 2nd peak
TString fFitOption
number of sigmas to veto the signal peak
TF1 * fBkRFunc
fit function for reflections
Int_t fNParsSec
switch off/on second peak (for D+->KKpi in Ds)
TF1 * CreateSignalFitFunction(TString fname, Double_t integral)
Bool_t fFixedMean
unc on signal gaussian sigma
Double_t fMassParticle
help variable
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
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)
Bool_t draw[nPtBins]
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 FitFunction4Sgn(Double_t *x, Double_t *par)
Int_t fPolDegreeBkg
background fit func
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 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
initialization for wa yield
Bool_t fReflections
background fit function (2nd step)
Bool_t fFixedSigma
switch for fix mean of gaussian