fNFinalPars(1),
fFitPars(0),
fWithBkg(0),
- ftypeOfFit4Bkg(0),
- ftypeOfFit4Sgn(0),
+ ftypeOfFit4Bkg(kExpo),
+ ftypeOfFit4Sgn(kGaus),
ffactor(1),
fntuParam(0),
fMass(1.865),
+ fMassErr(0.01),
fSigmaSgn(0.012),
+ fSigmaSgnErr(0.001),
+ fRawYield(0.),
+ fRawYieldErr(0.),
fSideBands(0),
fFixPar(0),
fSideBandl(0),
fSideBandr(0),
fcounter(0),
+ fFitOption("L,"),
fContourGraph(0)
{
// default constructor
//___________________________________________________________________________
-AliHFMassFitter::AliHFMassFitter (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
+AliHFMassFitter::AliHFMassFitter (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin, Int_t fittypeb, Int_t fittypes):
TNamed(),
fhistoInvMass(0),
fminMass(0),
fNFinalPars(1),
fFitPars(0),
fWithBkg(0),
- ftypeOfFit4Bkg(0),
- ftypeOfFit4Sgn(0),
+ ftypeOfFit4Bkg(kExpo),
+ ftypeOfFit4Sgn(kGaus),
ffactor(1),
fntuParam(0),
fMass(1.865),
+ fMassErr(0.01),
fSigmaSgn(0.012),
+ fSigmaSgnErr(0.001),
+ fRawYield(0.),
+ fRawYieldErr(0.),
fSideBands(0),
fFixPar(0),
fSideBandl(0),
fSideBandr(0),
fcounter(0),
+ fFitOption("L,"),
fContourGraph(0)
{
// standard constructor
ffactor(mfit.ffactor),
fntuParam(mfit.fntuParam),
fMass(mfit.fMass),
+ fMassErr(mfit.fMassErr),
fSigmaSgn(mfit.fSigmaSgn),
+ fSigmaSgnErr(mfit.fSigmaSgnErr),
+ fRawYield(mfit.fRawYield),
+ fRawYieldErr(mfit.fRawYieldErr),
fSideBands(mfit.fSideBands),
fFixPar(0),
fSideBandl(mfit.fSideBandl),
fSideBandr(mfit.fSideBandr),
fcounter(mfit.fcounter),
+ fFitOption(mfit.fFitOption),
fContourGraph(mfit.fContourGraph)
{
//copy constructor
ffactor= mfit.ffactor;
fntuParam= mfit.fntuParam;
fMass= mfit.fMass;
+ fMassErr= mfit.fMassErr;
fSigmaSgn= mfit.fSigmaSgn;
+ fSigmaSgnErr= mfit.fSigmaSgnErr;
+ fRawYield= mfit.fRawYield;
+ fRawYieldErr= mfit.fRawYieldErr;
fSideBands= mfit.fSideBands;
fSideBandl= mfit.fSideBandl;
fSideBandr= mfit.fSideBandr;
+ fFitOption= mfit.fFitOption;
fcounter= mfit.fcounter;
fContourGraph= mfit.fContourGraph;
//return kFALSE if something wrong
if(thispar>=fNFinalPars) {
- cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
+ AliError(Form("Error! Parameter out of bounds! Max is %d\n",fNFinalPars-1));
return kFALSE;
}
if(!fFixPar){
Bool_t AliHFMassFitter::GetFixThisParam(Int_t thispar)const{
//return the value of fFixPar[thispar]
if(thispar>=fNFinalPars) {
- cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
+ AliError(Form("Error! Parameter out of bounds! Max is %d\n",fNFinalPars-1));
return kFALSE;
}
if(!fFixPar) {
- cout<<"Error! Parameters to be fixed still not set"<<endl;
+ AliError("Error! Parameters to be fixed still not set");
return kFALSE;
}
// Rebin invariant mass histogram
if(!fhistoInvMass){
- cout<<"Histogram not set"<<endl;
+ AliError("Histogram not set!!");
return;
}
Int_t nbinshisto=fhistoInvMass->GetNbinsX();
// * [2] = sigma
//gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
+ // AliInfo("Signal function set to: Gaussian");
return par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
}
// * [1] = B;
//exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
total = par[0+firstPar]*par[1+firstPar]/(TMath::Exp(par[1+firstPar]*fmaxMass)-TMath::Exp(par[1+firstPar]*fminMass))*TMath::Exp(par[1+firstPar]*x[0]);
+ // AliInfo("Background function set to: exponential");
break;
case 1:
//linear
// * [0] = integralBkg;
// * [1] = b;
total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
+ // AliInfo("Background function set to: linear");
break;
case 2:
//polynomial
// * [1] = b;
// * [2] = c;
total = par[0+firstPar]/(fmaxMass-fminMass)+par[1]*(x[0]-0.5*(fmaxMass+fminMass))+par[2+firstPar]*(x[0]*x[0]-1/3.*(fmaxMass*fmaxMass*fmaxMass-fminMass*fminMass*fminMass)/(fmaxMass-fminMass));
+ // AliInfo("Background function set to: polynomial");
break;
case 3:
total=par[0+firstPar];
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]);
+ // AliInfo("Background function set to: powerlaw");
}
break;
case 5:
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));
+ // AliInfo("Background function set to: wit exponential");
}
break;
// default:
Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin+1);
- Double_t sidebandldouble,sidebandrdouble;
- Bool_t leftok=kFALSE, rightok=kFALSE;
+ Double_t sidebandldouble=0.,sidebandrdouble=0.;
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;
+ AliError("Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value");
return kFALSE;
}
//histo limit = fit function limit
- if((TMath::Abs(fminMass-minHisto) < 1e6 || TMath::Abs(fmaxMass - maxHisto) < 1e6) && (fMass-4.*fSigmaSgn-fminMass) < 1e6){
+ if((TMath::Abs(fminMass-minHisto) < 1e-6 || TMath::Abs(fmaxMass - maxHisto) < 1e-6) && (fMass-4.*fSigmaSgn-fminMass) < 1e-6){
Double_t coeff = (fMass-fminMass)/fSigmaSgn;
sidebandldouble=(fMass-0.5*coeff*fSigmaSgn);
sidebandrdouble=(fMass+0.5*coeff*fSigmaSgn);
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;
sidebandrdouble=fhistoInvMass->GetBinLowEdge(fSideBandr+1);
if(TMath::Abs(tmp-sidebandrdouble) > 1e-6){
- cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
- rightok=kTRUE;
+ AliWarning(Form("%f is not allowed, changing it to the nearest value allowed: \n",tmp));
}
cout<<sidebandrdouble<<" (bin "<<fSideBandr<<")"<<endl;
if (fSideBandl==0 || fSideBandr==fNbin) {
- cout<<"Error! Range too little";
+ AliError("Error! Range too little");
return kFALSE;
}
return kTRUE;
Double_t diffUnderBands=(intUnderHisto-intUnderFunc);
Double_t relDiff=diffUnderPick/diffUnderBands;
cout<<"Relative difference = "<<relDiff<<endl;
- if(TMath::Abs(relDiff) < 1) isBkgOnly=kTRUE;
+ if(TMath::Abs(relDiff) < 0.25) isBkgOnly=kTRUE;
else{
cout<<"Relative difference = "<<relDiff<<": I suppose there is some signal, continue with total fit!"<<endl;
}
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);
+ Double_t width=fhistoInvMass->GetBinWidth(1);
//cout<<"fNbin = "<<fNbin<<endl;
if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
//if only signal and reflection: skip
if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
ftypeOfFit4Sgn=0;
- fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
+ fhistoInvMass->Fit(bkgname.Data(),Form("R,%sE,0",fFitOption.Data()));
for(Int_t i=0;i<bkgPar;i++){
fFitPars[i]=funcbkg->GetParameter(i);
funcbkg1->SetLineColor(2); //red
- switch (ftypeOfFit4Bkg) {
+ switch (ftypeOfFit4Bkg) {
case 0:
{
cout<<"*** Exponential Fit ***"<<endl;
}
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;
+ //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");
+ Int_t status=fhistoInvMass->Fit(bkg1name.Data(),Form("R,%sE,+,0",fFitOption.Data()));
if (status != 0){
cout<<"Minuit returned "<<status<<endl;
return kFALSE;
Int_t status;
- status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
+ status = fhistoInvMass->Fit(massname.Data(),Form("R,%sE,+,0",fFitOption.Data()));
if (status != 0){
cout<<"Minuit returned "<<status<<endl;
return kFALSE;
cout<<"fit done"<<endl;
//reset value of fMass and fSigmaSgn to those found from fit
fMass=funcmass->GetParameter(fNFinalPars-2);
+ fMassErr=funcmass->GetParError(fNFinalPars-2);
fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
-
+ fSigmaSgnErr=funcmass->GetParError(fNFinalPars-1);
+ fRawYield=funcmass->GetParameter(fNFinalPars-3)/fhistoInvMass->GetBinWidth(1);
+ fRawYieldErr=funcmass->GetParError(fNFinalPars-3)/fhistoInvMass->GetBinWidth(1);
+
for(Int_t i=0;i<fNFinalPars;i++){
fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
}
- Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
+ Int_t status=fhistoInvMass->Fit(bkgname.Data(),Form("R,%sE,+,0",fFitOption.Data()));
if (status != 0){
cout<<"Minuit returned "<<status<<endl;
return kFALSE;
pinfo2->AddText(str);
str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
pinfo2->AddText(str);
- str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
+ if(bkg>0) str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
pinfo2->AddText(str);
pd->cd();
void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
// Return significance integral in a range
+ if(fcounter==0){
+ AliError("Number of fits is zero, check whether you made the fit before computing the significance!\n");
+ return;
+ }
+
Double_t signal,errsignal,background,errbackground;
Signal(min, max,signal,errsignal);
Background(min, max,background,errbackground);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff