* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
+
/////////////////////////////////////////////////////////////
//
// AliHFMassFitter for the fit of invariant mass distribution
#include <TMinuit.h>
#include <TStyle.h>
#include <TPaveText.h>
+#include <TDatabasePDG.h>
#include "AliHFMassFitter.h"
fhistoInvMass(0),
fminMass(0),
fmaxMass(0),
+ fminBinMass(0),
+ fmaxBinMass(0),
fNbin(0),
fParsSize(1),
fNFinalPars(1),
ftypeOfFit4Sgn(0),
ffactor(1),
fntuParam(0),
- fMass(1.85),
+ fMass(1.865),
fSigmaSgn(0.012),
fSideBands(0),
fFixPar(0),
fhistoInvMass(0),
fminMass(0),
fmaxMass(0),
+ fminBinMass(0),
+ fmaxBinMass(0),
fNbin(0),
fParsSize(1),
fNFinalPars(1),
ftypeOfFit4Sgn(0),
ffactor(1),
fntuParam(0),
- fMass(1.85),
+ fMass(1.865),
fSigmaSgn(0.012),
fSideBands(0),
fFixPar(0),
CheckRangeFit();
ftypeOfFit4Bkg=fittypeb;
ftypeOfFit4Sgn=fittypes;
- if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
+ if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2 && ftypeOfFit4Bkg!=4 && ftypeOfFit4Bkg!=5) fWithBkg=kFALSE;
else fWithBkg=kTRUE;
if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
fhistoInvMass(mfit.fhistoInvMass),
fminMass(mfit.fminMass),
fmaxMass(mfit.fmaxMass),
+ fminBinMass(mfit.fminBinMass),
+ fmaxBinMass(mfit.fmaxBinMass),
fNbin(mfit.fNbin),
fParsSize(mfit.fParsSize),
fNFinalPars(mfit.fNFinalPars),
memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Bool_t));
}
- //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
+
}
//_________________________________________________________________________
//destructor
cout<<"AliHFMassFitter destructor called"<<endl;
- if(fhistoInvMass) {
- cout<<"deleting histogram..."<<endl;
- delete fhistoInvMass;
- fhistoInvMass=NULL;
- }
- if(fntuParam){
- cout<<"deleting ntuple..."<<endl;
- delete fntuParam;
- fntuParam=NULL;
- }
+ delete fhistoInvMass;
- if(fFitPars) {
- delete[] fFitPars;
- cout<<"deleting parameter array..."<<endl;
- fFitPars=NULL;
- }
-
- if(fFixPar) {
- delete[] fFixPar;
- cout<<"deleting bool array..."<<endl;
- fFixPar=NULL;
- }
+ delete fntuParam;
+ delete[] fFitPars;
+
+ delete[] fFixPar;
+
fcounter = 0;
}
fContourGraph= mfit.fContourGraph;
if(mfit.fParsSize > 0){
- if(fFitPars) {
- delete[] fFitPars;
- fFitPars=NULL;
- }
+ delete[] fFitPars;
+
fFitPars=new Float_t[fParsSize];
memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
- if(fFixPar) {
- delete[] fFixPar;
- fFixPar=NULL;
- }
+ delete[] fFixPar;
+
fFixPar=new Bool_t[fNFinalPars];
memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Float_t));
}
-// fFitPars=new Float_t[fParsSize];
-// for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
return *this;
}
case 3:
fNFinalPars=1;
break;
+ case 4:
+ fNFinalPars=2;
+ break;
+ case 5:
+ fNFinalPars=3;
+ break;
default:
cout<<"Error in computing fNFinalPars: check ftypeOfFit4Bkg"<<endl;
break;
case 3:
fParsSize = 1*3;
break;
+ case 4:
+ fParsSize = 2*3;
+ break;
+ case 5:
+ fParsSize = 3*3;
+ break;
default:
cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
break;
//___________________________________________________________________________
void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
- //fhistoInvMass = (TH1F*)histoToFit->Clone();
+
fhistoInvMass = new TH1F(*histoToFit);
fhistoInvMass->SetDirectory(0);
- cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
+ //cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
}
//___________________________________________________________________________
- void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
-
- //set the type of fit to perform for signal and background
-
- ftypeOfFit4Bkg = fittypeb;
- ftypeOfFit4Sgn = fittypes;
-
- /*
- if(fFitPars) {
- delete[] fFitPars;
- fFitPars=NULL;
- }
- */
- ComputeParSize();
- fFitPars = new Float_t[fParsSize];
- /*
- if(fFixPar){
- delete[] fFixPar;
- fFixPar=NULL;
- }
- */
- SetDefaultFixParam();
-
-
+void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
+
+ //set the type of fit to perform for signal and background
+
+ ftypeOfFit4Bkg = fittypeb;
+ ftypeOfFit4Sgn = fittypes;
+
+ ComputeParSize();
+ fFitPars = new Float_t[fParsSize];
+
+ SetDefaultFixParam();
}
//___________________________________________________________________________
fMass=1.85;
fSigmaSgn=0.012;
cout<<"Reset "<<fhistoInvMass<<endl;
- if(fhistoInvMass) {
- //cout<<"esiste"<<endl;
- delete fhistoInvMass;
- fhistoInvMass=NULL;
- cout<<fhistoInvMass<<endl;
- }
- else cout<<"histogram doesn't exist, do not delete"<<endl;
-
-
+ delete fhistoInvMass;
}
//_________________________________________________________________________
}
}
+ //Power fit
+
+ // par[0] = tot integral
+ // par[1] = coef1
+ // par[2] = gaussian integral
+ // par[3] = gaussian mean
+ // par[4] = gaussian sigma
+
+ if (ftypeOfFit4Bkg==4) {
+
+ if(ftypeOfFit4Sgn == 0) {
+
+ Double_t parbkg[2] = {par[0]-par[2], par[1]};
+ bkg = FitFunction4Bkg(x,parbkg);
+ }
+ if(ftypeOfFit4Sgn == 1) {
+
+ Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-par[2], par[1]};
+ bkg = FitFunction4Bkg(x,parbkg);
+ }
+ sgn = FitFunction4Sgn(x,&par[2]);
+ }
+
+
+ //Power and exponential fit
+
+ // par[0] = tot integral
+ // par[1] = coef1
+ // par[2] = coef2
+ // par[3] = gaussian integral
+ // par[4] = gaussian mean
+ // par[5] = gaussian sigma
+
+ if (ftypeOfFit4Bkg==5) {
+
+ if(ftypeOfFit4Sgn == 0) {
+ Double_t parbkg[3] = {par[0]-par[3],par[1],par[2]};
+ bkg = FitFunction4Bkg(x,parbkg);
+ }
+ if(ftypeOfFit4Sgn == 1) {
+ Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-par[3], par[1], par[2]};
+ bkg = FitFunction4Bkg(x,parbkg);
+ }
+ sgn = FitFunction4Sgn(x,&par[3]);
+ }
+
total = bkg + sgn;
return total;
case 3:
total=par[0+firstPar];
break;
+ case 4:
+ //power function
+ //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
+ //
+ //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
+ // * [0] = integralBkg;
+ // * [1] = b;
+ // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
+ {
+ Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
+
+ total = par[0+firstPar]*(par[1+firstPar]+1.)/(TMath::Power(fmaxMass-mpi,par[1+firstPar]+1.)-TMath::Power(fminMass-mpi,par[1+firstPar]+1.))*TMath::Power(x[0]-mpi,par[1+firstPar]);
+ }
+ break;
+ case 5:
+ //power function wit exponential
+ //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
+ {
+ Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
+
+ total = par[1+firstPar]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2+firstPar]*(x[0]-mpi));
+ }
+ break;
// default:
// Types of Fit Functions for Background:
// * 0 = exponential;
// * 1 = linear;
// * 2 = polynomial 2nd order
// * 3 = no background"<<endl;
+// * 4 = Power function
+// * 5 = Power function with exponential
}
return total+gaus2;
//determines the ranges of the side bands
- Double_t width=fhistoInvMass->GetBinWidth(8);
if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
- Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
+ Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin+1);
+
+ Double_t sidebandldouble,sidebandrdouble;
+ Bool_t leftok=kFALSE, rightok=kFALSE;
if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
return kFALSE;
}
- if((TMath::Abs(fminMass-minHisto) < 10e6 || TMath::Abs(fmaxMass - maxHisto) < 10e6) && (fMass-4.*fSigmaSgn-fminMass) < 10e6){
+ //histo limit = fit function limit
+ if((TMath::Abs(fminMass-minHisto) < 1e6 || TMath::Abs(fmaxMass - maxHisto) < 1e6) && (fMass-4.*fSigmaSgn-fminMass) < 1e6){
Double_t coeff = (fMass-fminMass)/fSigmaSgn;
-
- fSideBandl=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
- fSideBandr=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
+ sidebandldouble=(fMass-0.5*coeff*fSigmaSgn);
+ sidebandrdouble=(fMass+0.5*coeff*fSigmaSgn);
cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
if (coeff<2) {
cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
ftypeOfFit4Bkg=3;
//set binleft and right without considering SetRangeFit- anyway no bkg!
- fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
- fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
+ sidebandldouble=(fMass-4.*fSigmaSgn);
+ sidebandrdouble=(fMass+4.*fSigmaSgn);
}
}
else {
- fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
- fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
-// cout<<"\tfMass = "<<fMass<<"\tfSigmaSgn = "<<fSigmaSgn<<"\tminHisto = "<<minHisto<<endl;
-// cout<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
+ sidebandldouble=(fMass-4.*fSigmaSgn);
+ sidebandrdouble=(fMass+4.*fSigmaSgn);
+ }
+
+ cout<<"Left side band ";
+ Double_t tmp=0.;
+ tmp=sidebandldouble;
+ //calculate bin corresponding to fSideBandl
+ fSideBandl=fhistoInvMass->FindBin(sidebandldouble);
+ if (sidebandldouble >= fhistoInvMass->GetBinCenter(fSideBandl)) fSideBandl++;
+ sidebandldouble=fhistoInvMass->GetBinLowEdge(fSideBandl);
+
+ if(TMath::Abs(tmp-sidebandldouble) > 1e-6){
+ cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
+ leftok=kTRUE;
}
+ cout<<sidebandldouble<<" (bin "<<fSideBandl<<")"<<endl;
+
+ cout<<"Right side band ";
+ tmp=sidebandrdouble;
+ //calculate bin corresponding to fSideBandr
+ fSideBandr=fhistoInvMass->FindBin(sidebandrdouble);
+ if (sidebandrdouble < fhistoInvMass->GetBinCenter(fSideBandr)) fSideBandr--;
+ sidebandrdouble=fhistoInvMass->GetBinLowEdge(fSideBandr+1);
- if (fSideBandl==0) {
+ if(TMath::Abs(tmp-sidebandrdouble) > 1e-6){
+ cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
+ rightok=kTRUE;
+ }
+ cout<<sidebandrdouble<<" (bin "<<fSideBandr<<")"<<endl;
+ if (fSideBandl==0 || fSideBandr==fNbin) {
cout<<"Error! Range too little";
return kFALSE;
}
}
Bool_t leftok=kFALSE, rightok=kFALSE;
Int_t nbins=fhistoInvMass->GetNbinsX();
- Double_t width=fhistoInvMass->GetBinWidth(1);
- Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins)+width;
+ Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins+1);
//check if limits are inside histogram range
fmaxMass=maxhisto;
}
- Int_t binl,binr;
Double_t tmp=0.;
tmp=fminMass;
//calculate bin corresponding to fminMass
- binl=fhistoInvMass->FindBin(fminMass);
- if (fminMass > fhistoInvMass->GetBinCenter(binl)) binl++;
- fminMass=fhistoInvMass->GetBinLowEdge(binl);
+ fminBinMass=fhistoInvMass->FindBin(fminMass);
+ if (fminMass >= fhistoInvMass->GetBinCenter(fminBinMass)) fminBinMass++;
+ fminMass=fhistoInvMass->GetBinLowEdge(fminBinMass);
if(TMath::Abs(tmp-fminMass) > 1e-6){
cout<<"Left bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fminMass<<endl;
leftok=kTRUE;
tmp=fmaxMass;
//calculate bin corresponding to fmaxMass
- binr=fhistoInvMass->FindBin(fmaxMass);
- if (fmaxMass < fhistoInvMass->GetBinCenter(binr)) binr--;
- fmaxMass=fhistoInvMass->GetBinLowEdge(binr)+width;
+ fmaxBinMass=fhistoInvMass->FindBin(fmaxMass);
+ if (fmaxMass < fhistoInvMass->GetBinCenter(fmaxBinMass)) fmaxBinMass--;
+ fmaxMass=fhistoInvMass->GetBinLowEdge(fmaxBinMass+1);
if(TMath::Abs(tmp-fmaxMass) > 1e-6){
cout<<"Right bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fmaxMass<<endl;
rightok=kTRUE;
//Set default fitter Minuit in order to use gMinuit in the contour plots
TVirtualFitter::SetDefaultFitter("Minuit");
- Int_t nFitPars=0; //total function's number of parameters
- switch (ftypeOfFit4Bkg){
- case 0:
- nFitPars=5; //3+2
- break;
- case 1:
- nFitPars=5; //3+2
- break;
- case 2:
- nFitPars=6; //3+3
- break;
- case 3:
- nFitPars=4; //3+1
- break;
+
+ Bool_t isBkgOnly=kFALSE;
+
+ Int_t fit1status=RefitWithBkgOnly(kFALSE);
+ if(fit1status){
+ Int_t checkinnsigma=4;
+ Double_t range[2]={fMass-checkinnsigma*fSigmaSgn,fMass+checkinnsigma*fSigmaSgn};
+ TF1* func=GetHistoClone()->GetFunction("funcbkgonly");
+ Double_t intUnderFunc=func->Integral(range[0],range[1]);
+ Double_t intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(range[0]),fhistoInvMass->FindBin(range[1]),"width");
+ cout<<"Pick zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
+ Double_t diffUnderPick=(intUnderHisto-intUnderFunc);
+ intUnderFunc=func->Integral(fminMass,fminMass+checkinnsigma*fSigmaSgn);
+ intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fminMass+checkinnsigma*fSigmaSgn),"width");
+ cout<<"Band (l) zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
+ Double_t diffUnderBands=(intUnderHisto-intUnderFunc);
+ Double_t relDiff=diffUnderPick/diffUnderBands;
+ cout<<"Relative difference = "<<relDiff<<endl;
+ if(TMath::Abs(relDiff) < 1) isBkgOnly=kTRUE;
+ else{
+ cout<<"Relative difference = "<<relDiff<<": I suppose there is some signal, continue with total fit!"<<endl;
+ }
}
+ if(isBkgOnly) {
+ cout<<"INFO!! The histogram contains only background"<<endl;
+ if(draw)DrawFit();
+
+ //increase counter of number of fits done
+ fcounter++;
- Int_t bkgPar = nFitPars-3; //background function's number of parameters
+ return kTRUE;
+ }
- cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
+ Int_t bkgPar = fNFinalPars-3; //background function's number of parameters
+
+ cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
TString listname="contourplot";
TString massname="funcmass";
//Total integral
- Double_t totInt = fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass), fhistoInvMass->FindBin(fmaxMass), "width");
-
+ Double_t totInt = fhistoInvMass->Integral(fminBinMass,fmaxBinMass, "width");
+ //cout<<"Here tot integral is = "<<totInt<<"; integral in whole range is "<<fhistoInvMass->Integral("width")<<endl;
fSideBands = kTRUE;
Double_t width=fhistoInvMass->GetBinWidth(8);
- cout<<"fNbin"<<fNbin<<endl;
+ //cout<<"fNbin = "<<fNbin<<endl;
if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
- //Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
- //Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
+
Bool_t ok=SideBandsBounds();
if(!ok) return kFALSE;
funcbkg->FixParameter(0,0.);
}
break;
+ case 4:
+ funcbkg->SetParNames("BkgInt","Coef2");
+ funcbkg->SetParameters(sideBandsInt,0.5);
+ break;
+ case 5:
+ funcbkg->SetParNames("BkgInt","Coef1","Coef2");
+ funcbkg->SetParameters(sideBandsInt, -10., 5.);
+ break;
default:
cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
return kFALSE;
for(Int_t i=0;i<bkgPar;i++){
fFitPars[i]=funcbkg->GetParameter(i);
//cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
- fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
- //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
+ fFitPars[fNFinalPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
+ //cout<<fNFinalPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
}
fSideBands = kFALSE;
//intbkg1 = funcbkg->GetParameter(0);
- funcbkg->SetRange(fminMass,fmaxMass);
+
intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
- cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
+ if(ftypeOfFit4Bkg==5) conc1 = funcbkg->GetParameter(2);
+
+
+ //cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
}
else cout<<"\t\t//"<<endl;
cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
bkgPar+=3;
- cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
+ //cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
cout<<"Function name = "<<funcbkg1->GetName()<<endl;
funcbkg1->SetLineColor(2); //red
- if(ftypeOfFit4Bkg==2){
- cout<<"*** Polynomial Fit ***"<<endl;
- funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
- funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
-
- //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
- //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
- } else{
- if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
+ switch (ftypeOfFit4Bkg) {
+ case 0:
{
- cout<<"*** No background Fit ***"<<endl;
- funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
- funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
- funcbkg1->FixParameter(3,0.);
- } else{ //expo or linear
- if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
- if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
- funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
- funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
+ cout<<"*** Exponential Fit ***"<<endl;
+ funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
+ funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
}
+ break;
+ case 1:
+ {
+ cout<<"*** Linear Fit ***"<<endl;
+ funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
+ funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
+ }
+ break;
+ case 2:
+ {
+ cout<<"*** Polynomial Fit ***"<<endl;
+ funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
+ funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
+ }
+ break;
+ case 3:
+ //no background: gaus sign+ gaus broadened
+ {
+ cout<<"*** No background Fit ***"<<endl;
+ funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
+ funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
+ funcbkg1->FixParameter(3,0.);
+ }
+ break;
+ case 4:
+ {
+ cout<<"*** Power function Fit ***"<<endl;
+ funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef2");
+ funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
+ }
+ break;
+ case 5:
+ {
+ cout<<"*** Power function conv. with exponential Fit ***"<<endl;
+ funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
+ funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
+ }
+ break;
}
+ //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
+ //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
+
Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
if (status != 0){
cout<<"Minuit returned "<<status<<endl;
for(Int_t i=0;i<bkgPar;i++){
fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
//cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
- fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
- //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
+ fFitPars[fNFinalPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
+ //cout<<"\t"<<fNFinalPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
}
intbkg1=funcbkg1->GetParameter(3);
if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
+ if(ftypeOfFit4Bkg==5) conc1 = funcbkg1->GetParameter(5);
+
} else {
bkgPar+=3;
for(Int_t i=0;i<3;i++){
fFitPars[bkgPar-3+i]=0.;
cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
- fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
- cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
+ fFitPars[fNFinalPars+3*bkgPar-6+i]= 0.;
+ cout<<fNFinalPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
}
for(Int_t i=0;i<bkgPar-3;i++){
fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
- fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
- cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
+ fFitPars[fNFinalPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
+ cout<<fNFinalPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
}
//sidebands integral - second approx (from fit)
fSideBands = kFALSE;
Double_t bkgInt;
- cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
+ //cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
- cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
+ //cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
//Signal integral - first approx
Double_t sgnInt;
sgnInt = totInt-bkgInt;
- cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
+ //cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
if (sgnInt <= 0){
cout<<"Setting sgnInt = - sgnInt"<<endl;
- sgnInt=- sgnInt;
+ sgnInt=(-1)*sgnInt;
}
- /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
- TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
+ /*Fit All Mass distribution with exponential + gaussian (+gaussian braodened) */
+ TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4MassDistr");
cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
funcmass->SetLineColor(4); //blue
//Set parameters
cout<<"\nTOTAL FIT"<<endl;
- if(nFitPars==5){
+ if(fNFinalPars==5){
funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
//cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
- //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
+ //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
if(fFixPar[0]){
- cout<<"fix1"<<endl;
funcmass->FixParameter(0,totInt);
}
if(fFixPar[1]){
- cout<<"fix2"<<endl;
funcmass->FixParameter(1,slope1);
}
if(fFixPar[2]){
- cout<<"fix3"<<endl;
funcmass->FixParameter(2,sgnInt);
}
if(fFixPar[3]){
- cout<<"fix4"<<endl;
funcmass->FixParameter(3,fMass);
}
if(fFixPar[4]){
- cout<<"fix5"<<endl;
funcmass->FixParameter(4,fSigmaSgn);
}
}
- if (nFitPars==6){
+ if (fNFinalPars==6){
funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
//cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
- //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
+ //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
if(fFixPar[0])funcmass->FixParameter(0,totInt);
if(fFixPar[1])funcmass->FixParameter(1,slope1);
if(fFixPar[2])funcmass->FixParameter(2,conc1);
if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
if(fFixPar[4])funcmass->FixParameter(4,fMass);
- if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
+ if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
//
//funcmass->FixParameter(2,sgnInt);
}
- if(nFitPars==4){
+ if(fNFinalPars==4){
funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
if(fFixPar[0]) funcmass->FixParameter(0,0.);
- //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
- //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
+ if(fFixPar[1])funcmass->FixParameter(1,sgnInt);
+ if(fFixPar[2])funcmass->FixParameter(2,fMass);
+ if(fFixPar[3])funcmass->FixParameter(3,fSigmaSgn);
+ //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
+ //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
}
- //funcmass->FixParameter(nFitPars-2,fMass);
- //funcmass->SetParLimits(nFitPars-1,0.007,0.05);
- //funcmass->SetParLimits(nFitPars-2,fMass-0.01,fMass+0.01);
+
Int_t status;
status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
cout<<"fit done"<<endl;
//reset value of fMass and fSigmaSgn to those found from fit
- fMass=funcmass->GetParameter(nFitPars-2);
- fSigmaSgn=funcmass->GetParameter(nFitPars-1);
+ fMass=funcmass->GetParameter(fNFinalPars-2);
+ fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
- for(Int_t i=0;i<nFitPars;i++){
+ for(Int_t i=0;i<fNFinalPars;i++){
fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
- fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
- //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
+ fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
+ //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
}
/*
//check: cout parameters
- for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
+ for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
cout<<i<<"\t"<<fFitPars[i]<<endl;
}
*/
-// if(draw){
-// TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
-// TH1F *fhistocopy=new TH1F(*fhistoInvMass);
-// canvas->cd();
-// fhistocopy->DrawClone();
-// if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
-// else funcbkg->DrawClone("sames");
-// funcmass->DrawClone("sames");
-// cout<<"Drawn"<<endl;
-// }
-
- if(funcmass->GetParameter(nFitPars-1) <0 || funcmass->GetParameter(nFitPars-2) <0 || funcmass->GetParameter(nFitPars-3) <0 ) {
+ if(funcmass->GetParameter(fNFinalPars-1) <0 || funcmass->GetParameter(fNFinalPars-2) <0 || funcmass->GetParameter(fNFinalPars-3) <0 ) {
cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
return kFALSE;
}
//contour plots
if(draw){
- for (Int_t kpar=1; kpar<nFitPars;kpar++){
+ for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
- for(Int_t jpar=kpar+1;jpar<nFitPars;jpar++){
+ for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
cout<<"Par "<<kpar<<" and "<<jpar<<endl;
// produce 2 contours per couple of parameters
}
- if (ftypeOfFit4Sgn == 1 && funcbkg1) {
+ if (ftypeOfFit4Sgn == 1) {
delete funcbkg1;
- funcbkg1=NULL;
- }
- if (funcbkg) {
- delete funcbkg;
- funcbkg=NULL;
- }
- if (funcmass) {
- delete funcmass;
- funcmass=NULL;
}
-
+ delete funcbkg;
+ delete funcmass;
+
AddFunctionsToHisto();
if (draw) DrawFit();
funcbkg->SetParameters(integral,-10.,5);
break;
case 3:
- cout<<"Warning! This choice does not have a lot of sense..."<<endl;
+ cout<<"Warning! This choice does not make a lot of sense..."<<endl;
if(ftypeOfFit4Sgn==0){
funcbkg->SetParNames("Const");
funcbkg->SetParameter(0,0.);
funcbkg->FixParameter(0,0.);
}
break;
+ case 4:
+ funcbkg->SetParNames("BkgInt","Coef1");
+ funcbkg->SetParameters(integral,0.5);
+ break;
+ case 5:
+ funcbkg->SetParNames("BkgInt","Coef1","Coef2");
+ funcbkg->SetParameters(integral,-10.,5.);
+ break;
default:
cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
return kFALSE;
}
//_________________________________________________________________________
Double_t AliHFMassFitter::GetChiSquare() const{
+ //Get Chi^2 method
TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
+ if(!funcmass) {
+ cout<<"funcmass not found"<<endl;
+ return -1;
+ }
return funcmass->GetChisquare();
}
//_________________________________________________________________________
Double_t AliHFMassFitter::GetReducedChiSquare() const{
+ //Get reduced Chi^2 method
TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
if(!funcmass) {
cout<<"funcmass not found"<<endl;
void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
//gives the integral of signal obtained from fit parameters
- if(!valuewitherror)valuewitherror=new Float_t[2];
+ if(!valuewitherror) {
+ printf("AliHFMassFitter::IntS: got a null pointer\n");
+ return;
+ }
Int_t index=fParsSize/2 - 3;
valuewitherror[0]=fFitPars[index];
index=fParsSize - 3;
valuewitherror[1]=fFitPars[index];
- }
+}
//_________________________________________________________________________
//Add the background function in the complete range to the list of functions attached to the histogram
- cout<<"AddFunctionsToHisto called"<<endl;
+ //cout<<"AddFunctionsToHisto called"<<endl;
TString bkgname = "funcbkg";
- Int_t np=-99;
- switch (ftypeOfFit4Bkg){
- case 0: //expo
- np=2;
- break;
- case 1: //linear
- np=2;
- break;
- case 2: //pol2
- np=3;
- break;
- case 3: //no bkg
- np=1;
- break;
- }
- if (ftypeOfFit4Sgn == 1) {
- bkgname += 1;
- np+=3;
- }
Bool_t done1=kFALSE,done2=kFALSE;
}else{
bonly->SetLineColor(kBlue+3);
hlist->Add((TF1*)bonly->Clone());
- if(bonly) {
- delete bonly;
- bonly=NULL;
- }
+ delete bonly;
}
}
}
bkgname += "FullRange";
- TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
+ TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
//cout<<bfullrange->GetName()<<endl;
- for(Int_t i=0;i<np;i++){
- //cout<<i<<" di "<<np<<endl;
+ for(Int_t i=0;i<fNFinalPars-3;i++){
bfullrange->SetParName(i,b->GetParName(i));
bfullrange->SetParameter(i,b->GetParameter(i));
bfullrange->SetParError(i,b->GetParError(i));
bkgnamesave += "Recalc";
- TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
+ TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
TF1 *mass=fhistoInvMass->GetFunction("funcmass");
return;
}
- blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(np));
- blastpar->SetParError(0,mass->GetParError(np));
- if (np>=2) {
+ //intBkg=intTot-intS
+ blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
+ blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
+ if (fNFinalPars>=5) {
blastpar->SetParameter(1,mass->GetParameter(1));
blastpar->SetParError(1,mass->GetParError(1));
}
- if (np==3) {
+ if (fNFinalPars==6) {
blastpar->SetParameter(2,mass->GetParameter(2));
blastpar->SetParError(2,mass->GetParError(2));
}
hlist->Add((TF1*)blastpar->Clone());
hlist->ls();
- if(bfullrange) {
- delete bfullrange;
- bfullrange=NULL;
- }
- if(blastpar) {
- delete blastpar;
- blastpar=NULL;
- }
+ delete bfullrange;
+ delete blastpar;
+
}
cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
- if(output) {
- delete output;
- output=NULL;
- }
-
+ delete output;
+
}
//_________________________________________________________________________
nget=NULL;
}
*/
- if(output) {
- delete output;
- output=NULL;
- }
+
+ delete output;
}
//_________________________________________________________________________
case 3:
type="noB"; //3+1
break;
+ case 4:
+ type="Pow"; //3+3
+ break;
+ case 5:
+ type="PowExp"; //3+3
+ break;
}
TString filename=Form("%sMassFit.root",type.Data());
TF1* f=hdraw->GetFunction("funcbkgonly");
for (Int_t i=0;i<fNFinalPars-3;i++){
pinfo->SetTextColor(kBlue+3);
- TString str=Form("%s = %f #pm %f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
+ TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
pinfo->AddText(str);
}
pd->cd();
hdraw->SetMarkerStyle(20);
hdraw->DrawClone("PE");
- hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
- hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
- hdraw->GetFunction("funcmass")->DrawClone("same");
+// if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
+// if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
+ if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
if(writeFitInfo > 0){
TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
pinfom->SetTextColor(kBlue);
- TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
+ TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
}
pd->cd();
if(writeFitInfo > 1){
for (Int_t i=0;i<fNFinalPars-3;i++){
pinfob->SetTextColor(kRed);
- TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
+ str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
pinfob->AddText(str);
}
pd->cd();
return;
}
- Int_t np=0;
- switch (ftypeOfFit4Bkg){
- case 0: //expo
- np=2;
- break;
- case 1: //linear
- np=2;
- break;
- case 2: //pol2
- np=3;
- break;
- case 3: //no bkg
- np=1;
- break;
- }
-
- np+=3; //3 parameter for signal
- if (ftypeOfFit4Sgn == 1) np+=3;
-
cout<<"Parameter Titles \n";
- for(Int_t i=0;i<np;i++){
+ for(Int_t i=0;i<fNFinalPars;i++){
cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
}
cout<<endl;
Signal(min,max,signal,errsignal);
-// //functions names
-// TString bkgname= "funcbkgRecalc";
-// TString bkg1name="funcbkg1Recalc";
-// TString massname="funcmass";
-
-
-// TF1 *funcbkg=0;
-// TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
-// if(!funcmass){
-// cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
-// return;
-// }
-
-// if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
-// else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
-
-// if(!funcbkg){
-// cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
-// return;
-// }
-
-// Int_t np=-99;
-// switch (ftypeOfFit4Bkg){
-// case 0: //expo
-// np=2;
-// break;
-// case 1: //linear
-// np=2;
-// break;
-// case 2: //pol2
-// np=3;
-// break;
-// case 3: //no bkg
-// np=1;
-// break;
-// }
-
-// Float_t intS,intSerr;
-
-// //relative error evaluation
-
-// intS=funcmass->GetParameter(np);
-// intSerr=funcmass->GetParError(np);
-
-// cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
-// Double_t background,errbackground;
-// Background(nOfSigma,background,errbackground);
-
-// //signal +/- error in nsigma
-// Double_t min=fMass-nOfSigma*fSigmaSgn;
-// Double_t max=fMass+nOfSigma*fSigmaSgn;
-// Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
-// signal=mass - background;
-// errsignal=TMath::Sqrt((intSerr/intS*mass)*(intSerr/intS*mass)/*assume relative error is the same as for total integral*/ + errbackground*errbackground);
return;
}
return;
}
- Int_t np=-99;
- switch (ftypeOfFit4Bkg){
- case 0: //expo
- np=2;
- break;
- case 1: //linear
- np=2;
- break;
- case 2: //pol2
- np=3;
- break;
- case 3: //no bkg
- np=1;
- break;
- }
+ Int_t np=fNFinalPars-3;
Double_t intS,intSerr;
Background(min,max,background,errbackground);
-// //functions names
-// TString bkgname="funcbkgRecalc";
-// TString bkg1name="funcbkg1Recalc";
-
-// TF1 *funcbkg=0;
-// if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
-// else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
-// if(!funcbkg){
-// cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
-// return;
-// }
-
-// Float_t intB,intBerr;//, intT,intTerr,intS,intSerr;
-
-// //relative error evaluation: from final parameters of the fit
-// if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
-// else{
-
-// intB=funcbkg->GetParameter(0);
-// intBerr=funcbkg->GetParError(0);
-
-// cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
-// }
-
-// Double_t min=fMass-nOfSigma*fSigmaSgn;
-// Double_t max=fMass+nOfSigma*fSigmaSgn;
-
-// //relative error evaluation: from histo
-
-// intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
-// Double_t sum2=0;
-
-// for(Int_t i=1;i<=fSideBandl;i++){
-// sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
-// }
-// for(Int_t i=fSideBandr;i<=fNbin;i++){
-// sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
-// }
-
-// intBerr=TMath::Sqrt(sum2);
-// cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
-
-// cout<<"Last estimation of bkg error is used"<<endl;
-
-// //backround +/- error in nsigma
-// if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
-// background = 0;
-// errbackground = 0;
-// }
-// else{
-// background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
-// errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
-// //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
-// }
return;
-
+
}
//___________________________________________________________________________
Double_t min=fMass-nOfSigma*fSigmaSgn;
Double_t max=fMass+nOfSigma*fSigmaSgn;
Significance(min, max, significance, errsignificance);
- /*
- Double_t signal,errsignal,background,errbackground;
- Signal(nOfSigma,signal,errsignal);
- Background(nOfSigma,background,errbackground);
- significance = signal/TMath::Sqrt(signal+background);
-
- errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
- */
return;
}
cout<<"Cannot calculate significance because of div by 0!"<<endl;
significance=-1;
errsignificance=0;
+ return;
}
significance = signal/TMath::Sqrt(signal+background);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff