* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
+
/////////////////////////////////////////////////////////////
//
// AliHFMassFitter for the fit of invariant mass distribution
// Author: C.Bianchin, chiara.bianchin@pd.infn.it
/////////////////////////////////////////////////////////////
+#include <Riostream.h>
+#include <TMath.h>
+#include <TNtuple.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TList.h>
+#include <TFile.h>
#include <TCanvas.h>
+#include <TVirtualPad.h>
+#include <TGraph.h>
+#include <TVirtualFitter.h>
+#include <TMinuit.h>
+#include <TStyle.h>
+#include <TPaveText.h>
+#include <TDatabasePDG.h>
#include "AliHFMassFitter.h"
+
ClassImp(AliHFMassFitter)
//************** constructors
fhistoInvMass(0),
fminMass(0),
fmaxMass(0),
+ fminBinMass(0),
+ fmaxBinMass(0),
fNbin(0),
fParsSize(1),
+ fNFinalPars(1),
fFitPars(0),
fWithBkg(0),
ftypeOfFit4Bkg(0),
ftypeOfFit4Sgn(0),
ffactor(1),
fntuParam(0),
- fMass(1.85),
+ fMass(1.865),
fSigmaSgn(0.012),
fSideBands(0),
+ fFixPar(0),
fSideBandl(0),
fSideBandr(0),
- fcounter(0)
+ fcounter(0),
+ fContourGraph(0)
{
// default constructor
+
+ cout<<"Default constructor:"<<endl;
+ cout<<"Remember to set the Histo, the Type, the FixPar"<<endl;
- cout<<"Default constructor"<<endl;
}
//___________________________________________________________________________
-AliHFMassFitter::AliHFMassFitter (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),
fmaxMass(0),
+ fminBinMass(0),
+ fmaxBinMass(0),
fNbin(0),
fParsSize(1),
+ fNFinalPars(1),
fFitPars(0),
fWithBkg(0),
ftypeOfFit4Bkg(0),
ftypeOfFit4Sgn(0),
ffactor(1),
fntuParam(0),
- fMass(1.85),
+ fMass(1.865),
fSigmaSgn(0.012),
fSideBands(0),
+ fFixPar(0),
fSideBandl(0),
fSideBandr(0),
- fcounter(0)
+ fcounter(0),
+ fContourGraph(0)
{
// standard constructor
fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
+ fhistoInvMass->SetDirectory(0);
fminMass=minvalue;
fmaxMass=maxvalue;
if(rebin!=1) RebinMass(rebin);
else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
+ 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;
ComputeParSize();
fFitPars=new Float_t[fParsSize];
+
+ SetDefaultFixParam();
+
+ fContourGraph=new TList();
+ fContourGraph->SetOwner();
+
}
fhistoInvMass(mfit.fhistoInvMass),
fminMass(mfit.fminMass),
fmaxMass(mfit.fmaxMass),
+ fminBinMass(mfit.fminBinMass),
+ fmaxBinMass(mfit.fmaxBinMass),
fNbin(mfit.fNbin),
fParsSize(mfit.fParsSize),
+ fNFinalPars(mfit.fNFinalPars),
fFitPars(0),
fWithBkg(mfit.fWithBkg),
ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
fMass(mfit.fMass),
fSigmaSgn(mfit.fSigmaSgn),
fSideBands(mfit.fSideBands),
+ fFixPar(0),
fSideBandl(mfit.fSideBandl),
fSideBandr(mfit.fSideBandr),
- fcounter(mfit.fcounter)
-
+ fcounter(mfit.fcounter),
+ fContourGraph(mfit.fContourGraph)
{
//copy constructor
if(mfit.fParsSize > 0){
fFitPars=new Float_t[fParsSize];
+ fFixPar=new Bool_t[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];
+
}
//_________________________________________________________________________
AliHFMassFitter::~AliHFMassFitter() {
+
+ //destructor
+
cout<<"AliHFMassFitter destructor called"<<endl;
- if(fhistoInvMass) delete fhistoInvMass;
- if(fntuParam) delete fntuParam;
- if(fFitPars) delete[] fFitPars;
+
+ delete fhistoInvMass;
+
+ delete fntuParam;
+
+ delete[] fFitPars;
+
+ delete[] fFixPar;
+
fcounter = 0;
}
fSideBandl= mfit.fSideBandl;
fSideBandr= mfit.fSideBandr;
fcounter= mfit.fcounter;
+ fContourGraph= mfit.fContourGraph;
+
if(mfit.fParsSize > 0){
- if(fFitPars) delete[] fFitPars;
+ delete[] fFitPars;
+
fFitPars=new Float_t[fParsSize];
memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
+
+ 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;
}
//__________________________________________________________________________
+void AliHFMassFitter::ComputeNFinalPars() {
+
+ //compute the number of parameters of the total (signal+bgk) function
+ cout<<"Info:ComputeNFinalPars... ";
+ switch (ftypeOfFit4Bkg) {//npar background func
+ case 0:
+ fNFinalPars=2;
+ break;
+ case 1:
+ fNFinalPars=2;
+ break;
+ case 2:
+ fNFinalPars=3;
+ break;
+ 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;
+ }
+
+ fNFinalPars+=3; //gaussian signal
+ cout<<": "<<fNFinalPars<<endl;
+}
+//__________________________________________________________________________
+
void AliHFMassFitter::ComputeParSize() {
- switch (ftypeOfFit4Bkg) {//npar backround func
+
+ //compute the size of the parameter array and set the data member
+
+ switch (ftypeOfFit4Bkg) {//npar background func
case 0:
fParsSize = 2*3;
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(TH1F *histoToFit){
- fhistoInvMass = (TH1F*)histoToFit->Clone();
+void AliHFMassFitter::SetDefaultFixParam(){
+
+ //Set default values for fFixPar (only total integral fixed)
+
+ ComputeNFinalPars();
+ fFixPar=new Bool_t[fNFinalPars];
+
+ fFixPar[0]=kTRUE; //default: IntTot fixed
+ cout<<"Parameter 0 is fixed"<<endl;
+ for(Int_t i=1;i<fNFinalPars;i++){
+ fFixPar[i]=kFALSE;
+ }
+
+}
+
+//___________________________________________________________________________
+Bool_t AliHFMassFitter::SetFixThisParam(Int_t thispar,Bool_t fixpar){
+
+ //set the value (kFALSE or kTRUE) of one element of fFixPar
+ //return kFALSE if something wrong
+
+ if(thispar>=fNFinalPars) {
+ cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
+ return kFALSE;
+ }
+ if(!fFixPar){
+ cout<<"Initializing fFixPar...";
+ SetDefaultFixParam();
+ cout<<" done."<<endl;
+ }
+
+ fFixPar[thispar]=fixpar;
+ if(fixpar)cout<<"Parameter "<<thispar<<" is now fixed"<<endl;
+ else cout<<"Parameter "<<thispar<<" is now free"<<endl;
+ return kTRUE;
+}
+
+//___________________________________________________________________________
+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;
+ return kFALSE;
+ }
+ if(!fFixPar) {
+ cout<<"Error! Parameters to be fixed still not set"<<endl;
+ return kFALSE;
+
+ }
+ return fFixPar[thispar];
+
}
//___________________________________________________________________________
+void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
+
+ fhistoInvMass = new TH1F(*histoToFit);
+ fhistoInvMass->SetDirectory(0);
+ //cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
+}
- void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
- ftypeOfFit4Bkg = fittypeb;
- ftypeOfFit4Sgn = fittypes;
- ComputeParSize();
- fFitPars = new Float_t[fParsSize];
+//___________________________________________________________________________
+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();
}
//___________________________________________________________________________
void AliHFMassFitter::Reset() {
+
+ //delete the histogram and reset the mean and sigma to default
+
+ cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
fMass=1.85;
fSigmaSgn=0.012;
+ cout<<"Reset "<<fhistoInvMass<<endl;
delete fhistoInvMass;
}
//_________________________________________________________________________
-void AliHFMassFitter::InitNtuParam(char *ntuname) {
+void AliHFMassFitter::InitNtuParam(TString ntuname) {
+
// Create ntuple to keep fit parameters
+
fntuParam=0;
- fntuParam=new TNtuple(ntuname,"Contains fit parameters","intbkg1:slope1:conc1:intGB:meanGB:sigmaGB:intbkg2:slope2:conc2:inttot:slope3:conc3:intsgn:meansgn:sigmasgn:intbkg1Err:slope1Err:conc1Err:intGBErr:meanGBErr:sigmaGBErr:intbkg2Err:slope2Err:conc2Err:inttotErr:slope3Err:conc3Err:intsgnErr:meansgnErr:sigmasgnErr");
+ fntuParam=new TNtuple(ntuname.Data(),"Contains fit parameters","intbkg1:slope1:conc1:intGB:meanGB:sigmaGB:intbkg2:slope2:conc2:inttot:slope3:conc3:intsgn:meansgn:sigmasgn:intbkg1Err:slope1Err:conc1Err:intGBErr:meanGBErr:sigmaGBErr:intbkg2Err:slope2Err:conc2Err:inttotErr:slope3Err:conc3Err:intsgnErr:meansgnErr:sigmasgnErr");
}
//_________________________________________________________________________
-TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
+TNtuple* AliHFMassFitter::NtuParamOneShot(TString ntuname){
// Create, fill and return ntuple with fit parameters
- InitNtuParam(ntuname);
+ InitNtuParam(ntuname.Data());
FillNtuParam();
return fntuParam;
}
void AliHFMassFitter::RebinMass(Int_t bingroup){
// Rebin invariant mass histogram
+ if(!fhistoInvMass){
+ cout<<"Histogram not set"<<endl;
+ return;
+ }
+ Int_t nbinshisto=fhistoInvMass->GetNbinsX();
if(bingroup<1){
cout<<"Error! Cannot group "<<bingroup<<" bins\n";
- fNbin=fhistoInvMass->GetNbinsX();
+ fNbin=nbinshisto;
cout<<"Kept original number of bins: "<<fNbin<<endl;
} else{
+
+ while(nbinshisto%bingroup != 0) {
+ bingroup--;
+ }
+ cout<<"Group "<<bingroup<<" bins"<<endl;
fhistoInvMass->Rebin(bingroup);
fNbin = fhistoInvMass->GetNbinsX();
cout<<"New number of bins: "<<fNbin<<endl;
}
-
}
if(ftypeOfFit4Sgn == 1) {
Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
- bkg = FitFunction4Bkg(x,parbkg);
+ bkg = FitFunction4Bkg(x,parbkg);
}
sgn = FitFunction4Sgn(x,&par[3]);
}
}
+ //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;
//__________________________________________________________________________
Bool_t AliHFMassFitter::SideBandsBounds(){
- Double_t width=fhistoInvMass->GetBinWidth(8);
+ //determines the ranges of the side bands
+
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((fminMass!=minHisto || fmaxMass!= maxHisto) && (fMass-4.*fSigmaSgn-fminMass) <= 0){
+ //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);
+ sidebandldouble=(fMass-4.*fSigmaSgn);
+ sidebandrdouble=(fMass+4.*fSigmaSgn);
}
- if (fSideBandl==0) {
+ 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(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;
}
- /*
- Int_t fbin=(Int_t)((fminMass-minHisto)/width)+1, lbin=(Int_t)((fmaxMass-minHisto)/width)+1;
- cout<<"Range: "<<fminMass<<", "<<fmaxMass<<endl;
- cout<<"Range in bin = "<<fbin<<" --> "<<lbin<<endl;
- */
return kTRUE;
}
+//__________________________________________________________________________
+
void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
+
+ // get the range of the side bands
+
if (fSideBandl==0 && fSideBandr==0){
cout<<"Use MassFitter method first"<<endl;
return;
right=fSideBandr;
}
+//__________________________________________________________________________
+Bool_t AliHFMassFitter::CheckRangeFit(){
+ //check if the limit of the range correspond to the limit of bins. If not reset the limit to the nearer value which satisfy this condition
+
+ if (!fhistoInvMass) {
+ cout<<"No histogram to fit! SetHisto(TH1F* h) before! "<<endl;
+ return kFALSE;
+ }
+ Bool_t leftok=kFALSE, rightok=kFALSE;
+ Int_t nbins=fhistoInvMass->GetNbinsX();
+ Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins+1);
+
+ //check if limits are inside histogram range
+
+ if( fminMass-minhisto < 0. ) {
+ cout<<"Out of histogram left bound! Setting to "<<minhisto<<endl;
+ fminMass=minhisto;
+ }
+ if( fmaxMass-maxhisto > 0. ) {
+ cout<<"Out of histogram right bound! Setting to"<<maxhisto<<endl;
+ fmaxMass=maxhisto;
+ }
+
+ Double_t tmp=0.;
+ tmp=fminMass;
+ //calculate bin corresponding to fminMass
+ 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
+ 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;
+ }
+
+ return (leftok && rightok);
+
+}
+
//__________________________________________________________________________
Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
// Main method of the class: performs the fit of the histogram
-
- 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;
+
+ //Set default fitter Minuit in order to use gMinuit in the contour plots
+ TVirtualFitter::SetDefaultFitter("Minuit");
+
+
+ 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++;
+
+ return kTRUE;
+ }
+
+ Int_t bkgPar = fNFinalPars-3; //background function's number of parameters
+
+ cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
+
+
+ TString listname="contourplot";
+ listname+=fcounter;
+ if(!fContourGraph){
+ fContourGraph=new TList();
+ fContourGraph->SetOwner();
}
- Int_t bkgPar = nFitPars-3; //background function's number of parameters
+ fContourGraph->SetName(listname);
- cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
//function names
TString bkgname="funcbkg";
TString massname="funcmass";
//Total integral
- Double_t totInt = fhistoInvMass->Integral("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;
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;
return kFALSE;
}
- /*
- //side bands usando il range utente
-Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(fbin,binleft,"width") + (Double_t)fhistoInvMass->Integral(binright,lbin,"width");
- cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<binleft<<"\tbinright = "<<binright<<endl;
- cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
-
- */
/*Fit Bkg*/
- //TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
- TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
+
+ TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
funcbkg->SetLineColor(2); //red
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<<"Primo 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<<"*** 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.);
- } 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<<"*** 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;
}
- fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
-
+ //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;
+ return kFALSE;
+ }
+
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;
-
- /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
- TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
+ //cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
+ if (sgnInt <= 0){
+ cout<<"Setting sgnInt = - sgnInt"<<endl;
+ sgnInt=(-1)*sgnInt;
+ }
+ /*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;
- funcmass->FixParameter(0,totInt);
+ //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,sgnInt);
+ }
+ if(fFixPar[3]){
+ funcmass->FixParameter(3,fMass);
+ }
+ if(fFixPar[4]){
+ 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;
- funcmass->FixParameter(0,totInt);
+
+ //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<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);
+ //
+ //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);
- 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[0]) funcmass->FixParameter(0,0.);
+ 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;
}
-
- fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
+
+ Int_t status;
+
+ status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
+ 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(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 ) {
- cout<<"IntS or mean or sigma negative"<<endl;
+ 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;
}
//increase counter of number of fits done
fcounter++;
- if (ftypeOfFit4Sgn == 1) delete funcbkg1;
+ //contour plots
+ if(draw){
+
+ for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
+
+ for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
+ cout<<"Par "<<kpar<<" and "<<jpar<<endl;
+
+ // produce 2 contours per couple of parameters
+ TGraph* cont[2] = {0x0, 0x0};
+ const Double_t errDef[2] = {1., 4.};
+ for (Int_t i=0; i<2; i++) {
+ gMinuit->SetErrorDef(errDef[i]);
+ cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
+ cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
+ }
+
+ if(!cont[0] || !cont[1]){
+ cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
+ continue;
+ }
+
+ // set graph titles and add them to the list
+ TString title = "Contour plot";
+ TString titleX = funcmass->GetParName(kpar);
+ TString titleY = funcmass->GetParName(jpar);
+ for (Int_t i=0; i<2; i++) {
+ cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
+ cont[i]->SetTitle(title);
+ cont[i]->GetXaxis()->SetTitle(titleX);
+ cont[i]->GetYaxis()->SetTitle(titleY);
+ cont[i]->GetYaxis()->SetLabelSize(0.033);
+ cont[i]->GetYaxis()->SetTitleSize(0.033);
+ cont[i]->GetYaxis()->SetTitleOffset(1.67);
+
+ fContourGraph->Add(cont[i]);
+ }
+
+ // plot them
+ TString cvname = Form("c%d%d", kpar, jpar);
+ TCanvas *c4=new TCanvas(cvname,cvname,600,600);
+ c4->cd();
+ cont[1]->SetFillColor(38);
+ cont[1]->Draw("alf");
+ cont[0]->SetFillColor(9);
+ cont[0]->Draw("lf");
+
+ }
+
+ }
+
+ }
+
+ if (ftypeOfFit4Sgn == 1) {
+ delete funcbkg1;
+ }
delete funcbkg;
delete funcmass;
+
+ AddFunctionsToHisto();
+ if (draw) DrawFit();
+
+
return kTRUE;
}
-//_________________________________________________________________________
-Double_t AliHFMassFitter::GetChiSquare() const{
- TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
- return funcmass->GetChisquare();
-}
-//_________________________________________________________________________
-Double_t AliHFMassFitter::GetReducedChiSquare() const{
- TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
- return funcmass->GetChisquare()/funcmass->GetNDF();
-}
+//______________________________________________________________________________
-//*********output
+Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
-//_________________________________________________________________________
-void AliHFMassFitter::GetFitPars(Float_t *vector) const {
- // Return fit parameters
+ //perform a fit with background function only. Can be useful to try when fit fails to understand if it is because there's no signal
+ //If you want to change the backgroud function or range use SetType or SetRangeFit before
- for(Int_t i=0;i<fParsSize;i++){
- vector[i]=fFitPars[i];
- }
-}
-//_________________________________________________________________________
+ TString bkgname="funcbkgonly";
+ fSideBands = kFALSE;
-TH1F* AliHFMassFitter::GetHistoClone() const{
- TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
- return hout;
-}
-//_________________________________________________________________________
+ TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
-void AliHFMassFitter::WriteHisto(TString path) {
- if (fcounter == 0) {
- cout<<"Use MassFitter method before WriteHisto"<<endl;
- return;
- }
+ funcbkg->SetLineColor(kBlue+3); //dark blue
- TH1F* hget=(TH1F*)fhistoInvMass->Clone();
+ Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
- TString bkgname = "funcbkg";
- Int_t np=-99;
- switch (ftypeOfFit4Bkg){
- case 0:
- np=2;
+ switch (ftypeOfFit4Bkg) {
+ case 0: //gaus+expo
+ funcbkg->SetParNames("BkgInt","Slope");
+ funcbkg->SetParameters(integral,-2.);
break;
case 1:
- np=2;
+ funcbkg->SetParNames("BkgInt","Slope");
+ funcbkg->SetParameters(integral,-100.);
break;
case 2:
- np=3;
+ funcbkg->SetParNames("BkgInt","Coef1","Coef2");
+ funcbkg->SetParameters(integral,-10.,5);
break;
case 3:
- np=1;
+ 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;
break;
}
- if (ftypeOfFit4Sgn == 1) {
- bkgname += 1;
- np+=3;
+
+
+ Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
+ if (status != 0){
+ cout<<"Minuit returned "<<status<<endl;
+ return kFALSE;
+ }
+ AddFunctionsToHisto();
+
+ if(draw) DrawFit();
+
+ return kTRUE;
+
+}
+//_________________________________________________________________________
+Double_t AliHFMassFitter::GetChiSquare() const{
+ //Get Chi^2 method
+ TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
+ if(!funcmass) {
+ cout<<"funcmass not found"<<endl;
+ return -1;
}
- TList *hlist=hget->GetListOfFunctions();
+ 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;
+ return -1;
+ }
+ return funcmass->GetChisquare()/funcmass->GetNDF();
+}
+
+//*********output
+
+//_________________________________________________________________________
+void AliHFMassFitter::GetFitPars(Float_t *vector) const {
+ // Return fit parameters
+
+ for(Int_t i=0;i<fParsSize;i++){
+ vector[i]=fFitPars[i];
+ }
+}
+
+
+//_________________________________________________________________________
+void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
+
+ //gives the integral of signal obtained from fit parameters
+ 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];
+}
+
+
+//_________________________________________________________________________
+void AliHFMassFitter::AddFunctionsToHisto(){
+
+ //Add the background function in the complete range to the list of functions attached to the histogram
+
+ //cout<<"AddFunctionsToHisto called"<<endl;
+ TString bkgname = "funcbkg";
+
+ Bool_t done1=kFALSE,done2=kFALSE;
+
+ TString bkgnamesave=bkgname;
+ TString testname=bkgname;
+ testname += "FullRange";
+ TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
+ if(testfunc){
+ done1=kTRUE;
+ testfunc=0x0;
+ }
+ testname="funcbkgonly";
+ testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
+ if(testfunc){
+ done2=kTRUE;
+ testfunc=0x0;
+ }
+
+ if(done1 && done2){
+ cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
+ return;
+ }
+
+ TList *hlist=fhistoInvMass->GetListOfFunctions();
hlist->ls();
- TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
- cout<< b->GetParameter(0)<<endl;
- cout<<"WRITEHISTO np = "<<np<<"\n"<<bkgname<<endl;
- bkgname += "FullRange";
- TF1 *bwrite=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
- cout<<bwrite->GetName()<<endl;
- for(Int_t i=0;i<np;i++){
- cout<<i<<" di "<<np<<endl;
- bwrite->SetParName(i,b->GetParName(i));
- bwrite->SetParameter(i,b->GetParameter(i));
+
+ if(!done2){
+ TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
+ if(!bonly){
+ cout<<testname.Data()<<" not found looking for complete fit"<<endl;
+ }else{
+ bonly->SetLineColor(kBlue+3);
+ hlist->Add((TF1*)bonly->Clone());
+ delete bonly;
+ }
+
}
+
+ if(!done1){
+ TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
+ if(!b){
+ cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
+ return;
+ }
+
+ bkgname += "FullRange";
+ 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<fNFinalPars-3;i++){
+ bfullrange->SetParName(i,b->GetParName(i));
+ bfullrange->SetParameter(i,b->GetParameter(i));
+ bfullrange->SetParError(i,b->GetParError(i));
+ }
+ bfullrange->SetLineStyle(4);
+ bfullrange->SetLineColor(14);
+
+ bkgnamesave += "Recalc";
+
+ TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
+
+ TF1 *mass=fhistoInvMass->GetFunction("funcmass");
+
+ if (!mass){
+ cout<<"funcmass doesn't exist "<<endl;
+ return;
+ }
+
+ //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 (fNFinalPars==6) {
+ blastpar->SetParameter(2,mass->GetParameter(2));
+ blastpar->SetParError(2,mass->GetParError(2));
+ }
+
+ blastpar->SetLineStyle(1);
+ blastpar->SetLineColor(2);
+
+ hlist->Add((TF1*)bfullrange->Clone());
+ hlist->Add((TF1*)blastpar->Clone());
+ hlist->ls();
- hlist->Add(bwrite);
+ delete bfullrange;
+ delete blastpar;
+
+ }
+
+
+}
+
+//_________________________________________________________________________
+
+TH1F* AliHFMassFitter::GetHistoClone() const{
+
+ TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
+ return hout;
+}
+//_________________________________________________________________________
+
+void AliHFMassFitter::WriteHisto(TString path) const {
+
+ //Write the histogram in the default file HFMassFitterOutput.root
+
+ if (fcounter == 0) {
+ cout<<"Use MassFitter method before WriteHisto"<<endl;
+ return;
+ }
+ TH1F* hget=(TH1F*)fhistoInvMass->Clone();
+
path += "HFMassFitterOutput.root";
TFile *output;
else output = new TFile(path.Data(),"update");
output->cd();
hget->Write();
+ fContourGraph->Write();
+
+
output->Close();
cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
- delete bwrite;
delete output;
-
+
}
//_________________________________________________________________________
void AliHFMassFitter::WriteNtuple(TString path) const{
- TNtuple* nget=(TNtuple*)fntuParam->Clone();
+ //TNtuple* nget=(TNtuple*)fntuParam->Clone();
path += "HFMassFitterOutput.root";
TFile *output = new TFile(path.Data(),"update");
output->cd();
- nget->Write();
+ fntuParam->Write();
+ //nget->Write();
output->Close();
- cout<<nget->GetName()<<" written in "<<path<<endl;
- delete nget;
+ //cout<<nget->GetName()<<" written in "<<path<<endl;
+ cout<<fntuParam->GetName()<<" written in "<<path<<endl;
+ /*
+ if(nget) {
+ //delete nget;
+ nget=NULL;
+ }
+ */
+
delete output;
}
+//_________________________________________________________________________
+void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
+
+ //write the canvas in a root file
-//************ significance
+ gStyle->SetOptStat(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetFrameFillColor(0);
+
+ TString type="";
+
+ switch (ftypeOfFit4Bkg){
+ case 0:
+ type="Exp"; //3+2
+ break;
+ case 1:
+ type="Lin"; //3+2
+ break;
+ case 2:
+ type="Pl2"; //3+3
+ break;
+ 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());
+ filename.Prepend(userIDstring);
+ path.Append(filename);
+
+ TFile* outputcv=new TFile(path.Data(),"update");
+
+ TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
+ c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
+ if(draw)c->DrawClone();
+ outputcv->cd();
+ c->Write();
+ outputcv->Close();
+}
//_________________________________________________________________________
-void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
- // Return signal integral in mean+- n sigma
+TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
+ //return a TVirtualPad with the fitted histograms and info
+ TString cvtitle="fit of ";
+ cvtitle+=fhistoInvMass->GetName();
+ TString cvname="c";
+ cvname+=fcounter;
- //functions names
- TString bkgname="funcbkg";
- TString bkg1name="funcbkg1";
- TString massname="funcmass";
+ TCanvas *c=new TCanvas(cvname,cvtitle);
+ PlotFit(c->cd(),nsigma,writeFitInfo);
+ return c->cd();
+}
+//_________________________________________________________________________
+void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
+ //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
+ gStyle->SetOptStat(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetFrameFillColor(0);
- TF1 *funcbkg=0;
- TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
- if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
- else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
+ cout<<"nsigma = "<<nsigma<<endl;
+ cout<<"Verbosity = "<<writeFitInfo<<endl;
+
+ TH1F* hdraw=GetHistoClone();
- //put final parameter in background function
- if (ftypeOfFit4Bkg != 3) {//signal + background
+ if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
+ cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
+ return;
+ }
+
+ if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
+ cout<<"Drawing background fit only"<<endl;
+ hdraw->SetMinimum(0);
+ hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
+ pd->cd();
+ hdraw->SetMarkerStyle(20);
+ hdraw->DrawClone("PE");
+ hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
+
+ if(writeFitInfo > 0){
+ TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
+ pinfo->SetBorderSize(0);
+ pinfo->SetFillStyle(0);
+ TF1* f=hdraw->GetFunction("funcbkgonly");
+ for (Int_t i=0;i<fNFinalPars-3;i++){
+ pinfo->SetTextColor(kBlue+3);
+ TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
+ pinfo->AddText(str);
+ }
+
+ pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
+ pd->cd();
+ pinfo->DrawClone();
- funcbkg->SetParameter(1,funcmass->GetParameter(1)); //slope
+
+ }
- if(ftypeOfFit4Bkg == 2) { //polynomial
- funcbkg->SetParameter(0,funcmass->GetParameter(0)-funcmass->GetParameter(3));//intbkg
- funcbkg->SetParameter(2,funcmass->GetParameter(2)); //coef
+ return;
+ }
+
+ hdraw->SetMinimum(0);
+ hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
+ pd->cd();
+ hdraw->SetMarkerStyle(20);
+ hdraw->DrawClone("PE");
+// 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");
+ TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
+ pinfob->SetBorderSize(0);
+ pinfob->SetFillStyle(0);
+ pinfom->SetBorderSize(0);
+ pinfom->SetFillStyle(0);
+ TF1* ff=fhistoInvMass->GetFunction("funcmass");
+
+ for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
+ pinfom->SetTextColor(kBlue);
+ 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);
}
- else {//expo or linear
- funcbkg->SetParameter(0,funcmass->GetParameter(0)-funcmass->GetParameter(2));//intbkg
+ pd->cd();
+ pinfom->DrawClone();
+
+ TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
+ pinfo2->SetBorderSize(0);
+ pinfo2->SetFillStyle(0);
+
+ Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
+
+ Significance(nsigma,signif,errsignif);
+ Signal(nsigma,signal,errsignal);
+ Background(nsigma,bkg, errbkg);
+ /*
+ Significance(1.828,1.892,signif,errsignif);
+ Signal(1.828,1.892,signal,errsignal);
+ Background(1.828,1.892,bkg, errbkg);
+ */
+ TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
+ pinfo2->AddText(str);
+ str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
+ pinfo2->AddText(str);
+ str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
+ pinfo2->AddText(str);
+
+ pd->cd();
+ pinfo2->Draw();
+
+ if(writeFitInfo > 1){
+ for (Int_t i=0;i<fNFinalPars-3;i++){
+ pinfob->SetTextColor(kRed);
+ str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
+ pinfob->AddText(str);
+ }
+ pd->cd();
+ pinfob->DrawClone();
}
+
+
}
+ return;
+}
- //Double_t mean=0,sigma=1;
- Double_t intS,intSerr;
+//_________________________________________________________________________
- //relative error evaluation
- if(ftypeOfFit4Bkg == 2) { //polynomial
- //mean=funcmass->GetParameter(4); //mean
- //sigma=funcmass->GetParameter(5); //sigma
- intS=fFitPars[12];
- intSerr=fFitPars[27];
- } else if(ftypeOfFit4Bkg == 3){ //no background
- //mean=funcmass->GetParameter(2); //mean
- //sigma=funcmass->GetParameter(3); //sigma
- intS=fFitPars[6];
- intSerr=fFitPars[15];
- } else { //expo or linear
- //mean=funcmass->GetParameter(3); //mean
- //sigma=funcmass->GetParameter(4); //sigma
- intS=fFitPars[9];
- intSerr=fFitPars[21];
+void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
+ //draws histogram together with fit functions with default nice colors in user canvas
+ PlotFit(pd,nsigma,writeFitInfo);
+
+ pd->Draw();
+
+}
+//_________________________________________________________________________
+void AliHFMassFitter::DrawFit(Double_t nsigma) const{
+
+ //draws histogram together with fit functions with default nice colors
+ gStyle->SetOptStat(0);
+ gStyle->SetCanvasColor(0);
+ gStyle->SetFrameFillColor(0);
+
+
+ TCanvas* c=(TCanvas*)GetPad(nsigma,1);
+ c->Draw();
+
+}
+
+
+//_________________________________________________________________________
+
+void AliHFMassFitter::PrintParTitles() const{
+
+ //prints to screen the parameters names
+
+ TF1 *f=fhistoInvMass->GetFunction("funcmass");
+ if(!f) {
+ cout<<"Fit function not found"<<endl;
+ return;
}
- cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
- //signal +/- error in nsigma
+
+ cout<<"Parameter Titles \n";
+ for(Int_t i=0;i<fNFinalPars;i++){
+ cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
+ }
+ cout<<endl;
+
+}
+
+//************ significance
+
+//_________________________________________________________________________
+
+void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
+ // Return signal integral in mean+- n sigma
+
+ if(fcounter==0) {
+ cout<<"Use MassFitter method before Signal"<<endl;
+ return;
+ }
+
Double_t min=fMass-nOfSigma*fSigmaSgn;
Double_t max=fMass+nOfSigma*fSigmaSgn;
- if(ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) signal=funcmass->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
- else signal=(funcmass->Integral(min,max)-funcbkg->Integral(min,max))/(Double_t)fhistoInvMass->GetBinWidth(2);
- errsignal=intSerr/intS*signal; // assume relative error is the same as for total integral
-
+
+ Signal(min,max,signal,errsignal);
+
+
return;
+
}
+
//_________________________________________________________________________
-void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
- // Return background integral in mean+- n sigma
-
+void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
+
+ // Return signal integral in a range
+
+ if(fcounter==0) {
+ cout<<"Use MassFitter method before Signal"<<endl;
+ return;
+ }
+
//functions names
- TString bkgname="funcbkg";
- TString bkg1name="funcbkg1";
+ 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());
- //relative error evaluation: from first fit
- Double_t intB,intBerr;
- intB=funcbkg->GetParameter(0);
- intBerr=funcbkg->GetParError(0);
- cout<<"Bkg relative error evaluation: from first fit: "<<intBerr/intB<<endl;
+ if(!funcbkg){
+ cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
+ return;
+ }
+ Int_t np=fNFinalPars-3;
- //put final parameter in background function
- if (ftypeOfFit4Bkg != 3) {//signal + background
+ Double_t intS,intSerr;
- funcbkg->SetParameter(1,funcmass->GetParameter(1)); //slope
+ //relative error evaluation
+ intS=funcmass->GetParameter(np);
+ intSerr=funcmass->GetParError(np);
- if(ftypeOfFit4Bkg == 2) { //polynomial
- funcbkg->SetParameter(0,funcmass->GetParameter(0)-funcmass->GetParameter(3));//intbkg
- funcbkg->SetParameter(2,funcmass->GetParameter(2)); //coef
- }
- else {//expo or linear
- funcbkg->SetParameter(0,funcmass->GetParameter(0)-funcmass->GetParameter(2));//intbkg
- }
- }
+ cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
+ Double_t background,errbackground;
+ Background(min,max,background,errbackground);
- //Double_t mean=0,sigma=1;
-
- /*
- if(ftypeOfFit4Bkg == 2) { //polynomial
- //mean=funcmass->GetParameter(4); //mean
- //sigma=funcmass->GetParameter(5); //sigma
- intB=fFitPars[9]-fFitPars[12];
- intBerr=fFitPars[27];
- } else if(ftypeOfFit4Bkg == 3){ //no background
- //mean=funcmass->GetParameter(2); //mean
- //sigma=funcmass->GetParameter(3); //sigma
- if (ftypeOfFit4Sgn == 1){ //reflection
- intB=fFitPars[1];
- intBerr=fFitPars[9];
- } else {
- intB=-1;intBerr=-1;
- err=0;
- }
- } else { //expo or linear
- //mean=funcmass->GetParameter(3); //mean
- //sigma=funcmass->GetParameter(4); //sigma
- intB=fFitPars[7]-fFitPars[9];
- intBerr=fFitPars[21];
- }
- */
-
+ //signal +/- error in the range
+
+ Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
+ signal=mass - background;
+ errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
+
+}
+
+//_________________________________________________________________________
+void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
+ // Return background integral in mean+- n sigma
+
+ if(fcounter==0) {
+ cout<<"Use MassFitter method before Background"<<endl;
+ return;
+ }
Double_t min=fMass-nOfSigma*fSigmaSgn;
Double_t max=fMass+nOfSigma*fSigmaSgn;
+
+ Background(min,max,background,errbackground);
+
+ return;
+}
+//___________________________________________________________________________
+
+void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
+ // Return background integral in a range
+
+ if(fcounter==0) {
+ cout<<"Use MassFitter method before Background"<<endl;
+ return;
+ }
+
+ //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;
+ }
+
+
+ Double_t intB,intBerr;
+
+ //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;
+ }
+
//relative error evaluation: from histo
intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
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
+ //backround +/- error in the range
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;
}
+
//__________________________________________________________________________
void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
// Return significance in mean+- n sigma
+
+ Double_t min=fMass-nOfSigma*fSigmaSgn;
+ Double_t max=fMass+nOfSigma*fSigmaSgn;
+ Significance(min, max, significance, errsignificance);
+
+ return;
+}
+
+//__________________________________________________________________________
+
+void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
+ // Return significance integral in a range
Double_t signal,errsignal,background,errbackground;
- Signal(nOfSigma,signal,errsignal);
- Background(nOfSigma,background,errbackground);
+ Signal(min, max,signal,errsignal);
+ Background(min, max,background,errbackground);
+
+ if (signal+background <= 0.){
+ cout<<"Cannot calculate significance because of div by 0!"<<endl;
+ significance=-1;
+ errsignificance=0;
+ return;
+ }
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;
}
-//__________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff