19 #include <TVirtualFitter.h> 20 #include <TDatabasePDG.h> 23 #include <TPaveText.h> 54 fTypeOfFit4Bkg(
kExpo),
58 fTypeOfFit4Sgn(
kGaus),
63 fSigmaSgn2Gaus(0.012),
66 fFixedSigma2Gaus(kFALSE),
69 fFixedFrac2Gaus(kFALSE),
71 fFixedRatio2GausSigma(kFALSE),
74 fOnlySideBands(kFALSE),
75 fNSigma4SideBands(4.),
86 fFixRflOverSig(kFALSE),
207 AliError(
"Error in computing fNParsBkg: check fTypeOfFit4Bkg");
222 AliError(
"Error in computing fNParsSig: check fTypeOfFit4Sgn");
240 TVirtualFitter::SetDefaultFitter(
"Minuit");
247 printf(
"\n--- First fit with only background on the side bands - Exclusion region = %.2f sigma ---\n",
fNSigma4SideBands);
258 printf(
" ---> Failed first fit with only background, minuit status = %d\n",status);
269 printf(
"\n--- Estimate signal counts in the peak region ---\n");
286 printf(
" ---> Final fit includes a second inv. mass peak\n");
291 printf(
" ---> Final fit includes reflections\n");
298 printf(
"\n--- Final fit with signal+background on the full range ---\n");
301 printf(
" ---> Failed fit with signal+background, minuit status = %d\n",status);
316 fSecFunc->SetParameter(ipar,
fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig));
317 fSecFunc->SetParError(ipar,
fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig));
346 if(doFinalFit)
return 1;
354 TCanvas* c0=
new TCanvas(
"c0");
362 gStyle->SetOptStat(0);
363 gStyle->SetCanvasColor(0);
364 gStyle->SetFrameFillColor(0);
377 if(writeFitInfo > 0){
378 TPaveText *pinfos=
new TPaveText(0.12,0.65,0.47,0.89,
"NDC");
379 TPaveText *pinfom=
new TPaveText(0.6,0.7,1.,.87,
"NDC");
380 pinfos->SetBorderSize(0);
381 pinfos->SetFillStyle(0);
382 pinfom->SetBorderSize(0);
383 pinfom->SetFillStyle(0);
385 pinfom->SetTextColor(kBlue);
397 pinfos->AddText(Form(
"B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
398 pinfos->AddText(Form(
"S/B (%.0f#sigma) = %.4f ",nsigma,
fRawYield/bkg));
400 pinfos->AddText(Form(
"Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
419 for(
Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
423 Double_t diffUnderPeak=(sum-sumback);
424 printf(
" ---> IntegralUnderFitFunc=%f IntegralUnderHisto=%f EstimatedSignal=%f\n",sum,sumback,diffUnderPeak);
425 if(diffUnderPeak/TMath::Sqrt(sum)<1.){
426 printf(
" ---> (Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal/\n",diffUnderPeak/TMath::Sqrt(sum));
441 funcbkg->SetParNames(
"BkgInt",
"Slope");
442 funcbkg->SetParameters(integral,-2.);
445 funcbkg->SetParNames(
"BkgInt",
"Slope");
446 funcbkg->SetParameters(integral,-100.);
449 funcbkg->SetParNames(
"BkgInt",
"Coef1",
"Coef2");
450 funcbkg->SetParameters(integral,-10.,5);
453 funcbkg->SetParNames(
"Const");
454 funcbkg->SetParameter(0,0.);
455 funcbkg->FixParameter(0,0.);
458 funcbkg->SetParNames(
"BkgInt",
"Coef1");
459 funcbkg->SetParameters(integral,0.5);
462 funcbkg->SetParNames(
"Coef1",
"Coef2");
463 funcbkg->SetParameters(-10.,5.);
467 funcbkg->SetParName(j,Form(
"Coef%d",j));
468 funcbkg->SetParameter(j,0);
470 funcbkg->SetParameter(0,integral);
473 AliError(Form(
"Wrong choice of fTypeOfFit4Bkg (%d)",
fTypeOfFit4Bkg));
479 funcbkg->SetLineColor(kBlue+3);
488 funcsec->SetParameter(0,integsig);
493 funcsec->SetParNames(
"SecPeakInt",
"SecMean",
"SecSigma");
502 funcrfl->SetParLimits(0,0.,1.);
504 funcrfl->SetParNames(
"ReflOverS");
515 fbr->SetParameter(ipar,
fBkgFunc->GetParameter(ipar));
516 fbr->SetParName(ipar,
fBkgFunc->GetParName(ipar));
519 fbr->SetParameter(ipar+fNParsBkg,
fRflFunc->GetParameter(ipar));
520 fbr->SetParName(ipar+fNParsBkg,
fRflFunc->GetParName(ipar));
533 funcsig->SetParameter(0,integsig);
535 funcsig->SetParameter(1,
fMass);
539 funcsig->SetParNames(
"SgnInt",
"Mean",
"Sigma");
542 funcsig->SetParameter(0,integsig);
544 funcsig->SetParameter(1,
fMass);
547 funcsig->SetParLimits(2,0.004,0.05);
551 else funcsig->SetParLimits(3,0.,1.);
554 else funcsig->SetParLimits(4,0.004,0.05);
555 funcsig->SetParNames(
"SgnInt",
"Mean",
"Sigma1",
"Frac",
"Sigma2");
558 funcsig->SetParameter(0,integsig);
560 funcsig->SetParameter(1,
fMass);
563 funcsig->SetParLimits(2,0.004,0.05);
567 else funcsig->SetParLimits(3,0.,1.);
570 else funcsig->SetParLimits(4,0.,20.);
571 funcsig->SetParNames(
"SgnInt",
"Mean",
"Sigma1",
"Frac",
"RatioSigma12");
585 ftot->SetParameter(ipar,
fBkgFunc->GetParameter(ipar));
586 ftot->SetParName(ipar,
fBkgFunc->GetParName(ipar));
592 ftot->SetParameter(ipar+fNParsBkg,
fSigFunc->GetParameter(ipar));
593 ftot->SetParName(ipar+fNParsBkg,
fSigFunc->GetParName(ipar));
595 fSigFunc->GetParLimits(ipar,parmin,parmax);
596 ftot->SetParLimits(ipar+fNParsBkg,parmin,parmax);
600 ftot->SetParameter(ipar+fNParsBkg+fNParsSig,
fSecFunc->GetParameter(ipar));
601 ftot->SetParName(ipar+fNParsBkg+fNParsSig,
fSecFunc->GetParName(ipar));
603 fSecFunc->GetParLimits(ipar,parmin,parmax);
604 ftot->SetParLimits(ipar+fNParsBkg+fNParsSig,parmin,parmax);
609 ftot->SetParameter(ipar+fNParsBkg+fNParsSig+
fNParsSec,
fRflFunc->GetParameter(ipar));
612 fRflFunc->GetParLimits(ipar,parmin,parmax);
613 ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+
fNParsSec,parmin,parmax);
644 total = par[0]*par[1]/(TMath::Exp(par[1]*
fMaxMass)-TMath::Exp(par[1]*
fMinMass))*TMath::Exp(par[1]*x[0]);
652 total= par[0]/(fMaxMass-
fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
663 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));
678 Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
680 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]);
688 Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
690 total = par[0]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[1]*(x[0]-mpi));
707 total+=par[it]*TMath::Power(x[0]-
fMassParticle,it)/TMath::Factorial(it);
731 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]);
742 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]);
743 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]);
744 sigval=par[0]*(g1+g2);
754 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]);
755 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]));
756 sigval=par[0]*(g1+g2);
798 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]);
812 return bkg+sig+sec+refl;
822 Signal(minMass,maxMass,signal,errsignal);
843 Background(minMass,maxMass,background,errbackground);
856 AliError(
"Bkg function not found!");
860 Double_t intB=funcbkg->GetParameter(0);
861 Double_t intBerr=funcbkg->GetParError(0);
869 for(
Int_t i=1;i<=leftBand;i++){
876 intBerr=TMath::Sqrt(sum2);
879 errbackground=intBerr/intB*background;
892 Significance(minMass, maxMass, significance, errsignificance);
926 TF1 *funcbkg,*funcPrev=0x0;
931 funcbkg->SetParameter(j,funcPrev->GetParameter(j));
936 funcbkg->SetParameter(0,estimatecent);
937 funcbkg->SetParameter(1,estimateslope);
939 printf(
" ---> Pre-fit of background with pol degree %d ---\n",
fCurPolDegreeBkg);
940 hCp->Fit(funcbkg,
"REMN",
"");
941 funcPrev=(TF1*)funcbkg->Clone(
"ftemp");
947 fback->SetParameter(j,funcPrev->GetParameter(j));
948 fback->SetParError(j,funcPrev->GetParError(j));
950 printf(
" ---> Final background fit with pol degree %d ---\n",
fPolDegreeBkg);
951 hCp->Fit(fback,
"REMN",
"");
979 back+=par[it]*TMath::Power(x[0]-
fMassParticle,it)/TMath::Factorial(it);
1000 printf(
"\n--- Reflection templates from simulation ---\n");
1001 if(opt.Contains(
"templ")){
1002 printf(
" ---> Reflection contribution using directly the histogram from simulation\n");
1008 Double_t xMinForFit=h->GetBinLowEdge(1);
1009 Double_t xMaxForFit=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
1010 if(minRange>=0 && maxRange>=0){
1011 xMinForFit=TMath::Max(minRange,h->GetBinLowEdge(1));
1012 xMaxForFit=TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1014 if(opt.EqualTo(
"1gaus") || opt.EqualTo(
"singlegaus")){
1015 printf(
" ---> Reflection contribution from single-Gaussian fit to histogram from simulation\n");
1016 f=
new TF1(
"mygaus",
"gaus",xMinForFit,xMaxForFit);
1017 f->SetParameter(0,h->GetMaximum());
1019 f->SetParameter(1,1.865);
1020 f->SetParameter(2,0.050);
1023 else if(opt.EqualTo(
"2gaus") || opt.EqualTo(
"doublegaus")){
1024 printf(
" ---> Reflection contribution from double-Gaussian fit to histogram from simulation\n");
1025 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);
1026 f->SetParameter(0,h->GetMaximum());
1028 f->SetParLimits(3,0.,1.);
1029 f->SetParameter(3,0.5);
1030 f->SetParameter(1,1.84);
1031 f->SetParameter(2,0.050);
1032 f->SetParameter(4,1.88);
1033 f->SetParameter(5,0.050);
1036 else if(opt.EqualTo(
"pol3")){
1037 printf(
" ---> Reflection contribution from pol3 fit to histogram from simulation\n");
1038 f=
new TF1(
"mypol3",
"pol3",xMinForFit,xMaxForFit);
1039 f->SetParameter(0,h->GetMaximum());
1045 else if(opt.EqualTo(
"pol6")){
1046 printf(
" ---> Reflection contribution from pol6 fit to histogram from simulation\n");
1047 f=
new TF1(
"mypol6",
"pol6",xMinForFit,xMaxForFit);
1048 f->SetParameter(0,h->GetMaximum());
1055 printf(
" ---> Bad option for reflection configuration -> reflections will not be included in the fit\n");
1066 if(
fHistoTemplRfl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
1077 printf(
" ---> Fit to MC template for reflection failed -> reflections will not be included in the fit\n");
1098 massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1101 massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1102 massVal-=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1123 printf(
"Left range for bin counting smaller than allowed by histogram axis, setting it to the lower edge of the first histo bin\n");
1127 printf(
"Right range for bin counting larger than allowed by histogram axis, setting it to the upper edge of the last histo bin\n");
1135 if(!fbackground)
return 0.;
1137 for(
Int_t jb=minBinSum; jb<=maxBinSum; jb++){
1145 cntSig+=(cntTot-cntBkg-cntRefl-cntSecPeak);
1148 errRyBC=TMath::Sqrt(cntErr);
1160 if(maxrange>minrange){
1167 hResidualTrend->SetName(Form(
"%s_residualTrend",
fHistoInvMass->GetName()));
1168 hResidualTrend->SetTitle(Form(
"%s (Residuals)",
fHistoInvMass->GetTitle()));
1169 hResidualTrend->SetMarkerStyle(20);
1170 hResidualTrend->SetMarkerSize(1.0);
1171 hResidualTrend->Reset();
1175 hPullsTrend->Reset();
1176 hPullsTrend->SetName(Form(
"%s_pullTrend",
fHistoInvMass->GetName()));
1177 hPullsTrend->SetTitle(Form(
"%s (Pulls)",
fHistoInvMass->GetTitle()));
1178 hPullsTrend->SetMarkerStyle(20);
1179 hPullsTrend->SetMarkerSize(1.0);
1183 hPulls->SetTitle(Form(
"%s ; Pulls",
fHistoInvMass->GetTitle()));
1184 hPulls->SetBins(40,-10,10);
1188 Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
1193 if(jst>=binmi&&jst<=binma){
1194 arval->AddAt(res,jst-binmi);
1200 hResidualTrend->SetBinContent(jst,res);
1201 hResidualTrend->SetBinError(jst,
fHistoInvMass->GetBinError(jst));
1204 if(jst>=binmi&&jst<=binma)hPulls->Fill(res/
fHistoInvMass->GetBinError(jst));
1207 hPullsTrend->SetBinContent(jst,res/
fHistoInvMass->GetBinError(jst));
1208 hPullsTrend->SetBinError(jst,0.0001);
1211 if(hResidualTrend)hResidualTrend->GetXaxis()->SetRange(binmi,binma);
1213 hPullsTrend->GetXaxis()->SetRange(binmi,binma);
1214 hPullsTrend->SetMinimum(-7);
1215 hPullsTrend->SetMaximum(+7);
1217 if(TMath::Abs(min)>TMath::Abs(max))max=min;
1219 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);
1220 for(
Int_t j=0;j<binma-binmi+1;j++){
1221 hout->Fill(arval->At(j));
1224 hout->Fit(
"gaus",
"LEM",
"",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
1228 hPulls->Fit(
"gaus",
"LEM",
"",-3,3);
1238 printf(
"--- Background function in 1st step fit to side bands ---\n");
1242 printf(
"--- Total fit function from 2nd step fit ---\n");
1246 printf(
"--- Background function in 2nd step fit ---\n");
1250 printf(
"--- Reflections ---\n");
1254 printf(
"--- Background + reflections ---\n");
1258 printf(
"--- Additional Gaussian peak ---\n");
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_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
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
TF1 * CreateSecondPeakFunction(TString fname, Double_t integral)
Bool_t fFixedFrac2Gaus
initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar ...
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
Bool_t fFixSecMass
width of the 2nd peak
TF1 * fSecFunc
flag to fix the width of the 2nd peak
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
Bool_t fFixedSigma2Gaus
switch for fix Sigma of gaussian
Double_t fSecWidth
position of the 2nd peak
TString fFitOption
number of sigmas to veto the signal peak
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 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)
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)
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 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
Double_t fSigmaSgn2Gaus
unc on signal gaussian sigma
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
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 ...
Bool_t fReflections
background fit function (2nd step)
Bool_t fFixedSigma
switch for fix mean of gaussian