1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // AliHFMassFitter for the fit of invariant mass distribution
21 // Author: C.Bianchin, chiara.bianchin@pd.infn.it
22 /////////////////////////////////////////////////////////////
26 #include "AliHFMassFitter.h"
29 ClassImp(AliHFMassFitter)
31 //************** constructors
32 AliHFMassFitter::AliHFMassFitter() :
50 // default constructor
52 cout<<"Default constructor"<<endl;
55 //___________________________________________________________________________
57 AliHFMassFitter::AliHFMassFitter (TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
75 // standard constructor
77 fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
80 if(rebin!=1) RebinMass(rebin);
81 else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
82 ftypeOfFit4Bkg=fittypeb;
83 ftypeOfFit4Sgn=fittypes;
84 if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
86 if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
87 else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
90 fFitPars=new Float_t[fParsSize];
94 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
96 fhistoInvMass(mfit.fhistoInvMass),
97 fminMass(mfit.fminMass),
98 fmaxMass(mfit.fmaxMass),
100 fParsSize(mfit.fParsSize),
102 fWithBkg(mfit.fWithBkg),
103 ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
104 ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
105 ffactor(mfit.ffactor),
106 fntuParam(mfit.fntuParam),
108 fSigmaSgn(mfit.fSigmaSgn),
109 fSideBands(mfit.fSideBands),
110 fcounter(mfit.fcounter)
115 if(mfit.fParsSize > 0){
116 fFitPars=new Float_t[fParsSize];
117 memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
119 //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
122 //_________________________________________________________________________
124 AliHFMassFitter::~AliHFMassFitter() {
125 cout<<"AliHFMassFitter destructor called"<<endl;
126 if(fhistoInvMass) delete fhistoInvMass;
127 if(fntuParam) delete fntuParam;
128 if(fFitPars) delete[] fFitPars;
132 //_________________________________________________________________________
134 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
136 //assignment operator
138 if(&mfit == this) return *this;
139 fhistoInvMass= mfit.fhistoInvMass;
140 fminMass= mfit.fminMass;
141 fmaxMass= mfit.fmaxMass;
143 fParsSize= mfit.fParsSize;
144 fWithBkg= mfit.fWithBkg;
145 ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
146 ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
147 ffactor= mfit.ffactor;
148 fntuParam= mfit.fntuParam;
150 fSigmaSgn= mfit.fSigmaSgn;
151 fSideBands= mfit.fSideBands;
152 fcounter= mfit.fcounter;
153 if(mfit.fParsSize > 0){
154 if(fFitPars) delete[] fFitPars;
155 fFitPars=new Float_t[fParsSize];
156 memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
158 // fFitPars=new Float_t[fParsSize];
159 // for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
164 //************ tools & settings
166 //__________________________________________________________________________
168 void AliHFMassFitter::ComputeParSize() {
169 switch (ftypeOfFit4Bkg) {//npar backround func
183 cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
187 fParsSize += 3; // npar refl
188 fParsSize += 3; // npar signal gaus
190 fParsSize*=2; // add errors
191 cout<<"Parameters array size "<<fParsSize<<endl;
194 //___________________________________________________________________________
195 void AliHFMassFitter::SetHisto(TH1F *histoToFit){
196 fhistoInvMass = (TH1F*)histoToFit->Clone();
199 //___________________________________________________________________________
201 void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
202 ftypeOfFit4Bkg = fittypeb;
203 ftypeOfFit4Sgn = fittypes;
205 fFitPars = new Float_t[fParsSize];
209 //___________________________________________________________________________
211 void AliHFMassFitter::Reset() {
214 delete fhistoInvMass;
217 //_________________________________________________________________________
219 void AliHFMassFitter::InitNtuParam(char *ntuname) {
220 // Create ntuple to keep fit parameters
222 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");
226 //_________________________________________________________________________
228 void AliHFMassFitter::FillNtuParam() {
229 // Fill ntuple with fit parameters
233 if (ftypeOfFit4Bkg==2) {
234 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
235 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
236 fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
237 fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
238 fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
239 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
240 fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
241 fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
242 fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
243 fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
244 fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
245 fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
246 fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
247 fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
248 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
250 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
251 fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
252 fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
253 fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
254 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
255 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
256 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
257 fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
258 fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
259 fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
260 fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
261 fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
262 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
263 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
264 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
268 if(ftypeOfFit4Bkg==3){
269 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
270 fntuParam->SetBranchAddress("slope1",¬hing);
271 fntuParam->SetBranchAddress("conc1",¬hing);
272 fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
273 fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
274 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
275 fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
276 fntuParam->SetBranchAddress("slope2",¬hing);
277 fntuParam->SetBranchAddress("conc2",¬hing);
278 fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
279 fntuParam->SetBranchAddress("slope3",¬hing);
280 fntuParam->SetBranchAddress("conc3",¬hing);
281 fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
282 fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
283 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
285 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
286 fntuParam->SetBranchAddress("slope1Err",¬hing);
287 fntuParam->SetBranchAddress("conc1Err",¬hing);
288 fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
289 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
290 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
291 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
292 fntuParam->SetBranchAddress("slope2Err",¬hing);
293 fntuParam->SetBranchAddress("conc2Err",¬hing);
294 fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
295 fntuParam->SetBranchAddress("slope3Err",¬hing);
296 fntuParam->SetBranchAddress("conc3Err",¬hing);
297 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
298 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
299 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
303 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
304 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
305 fntuParam->SetBranchAddress("conc1",¬hing);
306 fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
307 fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
308 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
309 fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
310 fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
311 fntuParam->SetBranchAddress("conc2",¬hing);
312 fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
313 fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
314 fntuParam->SetBranchAddress("conc3",¬hing);
315 fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
316 fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
317 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
319 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
320 fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
321 fntuParam->SetBranchAddress("conc1Err",¬hing);
322 fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
323 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
324 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
325 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
326 fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
327 fntuParam->SetBranchAddress("conc2Err",¬hing);
328 fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
329 fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
330 fntuParam->SetBranchAddress("conc3Err",¬hing);
331 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
332 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
333 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
337 fntuParam->TTree::Fill();
340 //_________________________________________________________________________
342 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
343 // Create, fill and return ntuple with fit parameters
345 InitNtuParam(ntuname);
349 //_________________________________________________________________________
351 void AliHFMassFitter::RebinMass(Int_t bingroup){
352 // Rebin invariant mass histogram
355 cout<<"Error! Cannot group "<<bingroup<<" bins\n";
356 fNbin=fhistoInvMass->GetNbinsX();
357 cout<<"Kept original number of bins: "<<fNbin<<endl;
359 fhistoInvMass->Rebin(bingroup);
360 fNbin = fhistoInvMass->GetNbinsX();
361 cout<<"New number of bins: "<<fNbin<<endl;
369 //___________________________________________________________________________
371 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
372 // Fit function for signal+background
375 //exponential or linear fit
377 // par[0] = tot integral
379 // par[2] = gaussian integral
380 // par[3] = gaussian mean
381 // par[4] = gaussian sigma
383 Double_t total,bkg=0,sgn=0;
385 if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
386 if(ftypeOfFit4Sgn == 0) {
388 Double_t parbkg[2] = {par[0]-par[2], par[1]};
389 bkg = FitFunction4Bkg(x,parbkg);
391 if(ftypeOfFit4Sgn == 1) {
392 Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
393 bkg = FitFunction4Bkg(x,parbkg);
396 sgn = FitFunction4Sgn(x,&par[2]);
402 // par[0] = tot integral
405 // par[3] = gaussian integral
406 // par[4] = gaussian mean
407 // par[5] = gaussian sigma
409 if (ftypeOfFit4Bkg==2) {
411 if(ftypeOfFit4Sgn == 0) {
413 Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
414 bkg = FitFunction4Bkg(x,parbkg);
416 if(ftypeOfFit4Sgn == 1) {
418 Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
419 bkg = FitFunction4Bkg(x,parbkg);
422 sgn = FitFunction4Sgn(x,&par[3]);
425 if (ftypeOfFit4Bkg==3) {
427 if(ftypeOfFit4Sgn == 0) {
428 bkg=FitFunction4Bkg(x,par);
429 sgn=FitFunction4Sgn(x,&par[1]);
431 if(ftypeOfFit4Sgn == 1) {
432 Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
433 bkg=FitFunction4Bkg(x,parbkg);
434 sgn=FitFunction4Sgn(x,&par[1]);
443 //_________________________________________________________________________
444 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
445 // Fit function for the signal
447 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
449 // * [0] = integralSgn
452 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
454 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]);
458 //__________________________________________________________________________
460 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
461 // Fit function for the background
463 Double_t maxDeltaM = 4.*fSigmaSgn;
464 if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
469 Double_t gaus2=0,total=-1;
470 if(ftypeOfFit4Sgn == 1){
472 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
474 // * [0] = integralSgn
477 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
478 gaus2 = FitFunction4Sgn(x,par);
481 switch (ftypeOfFit4Bkg){
484 //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
485 //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
486 //as integralTot- integralGaus (=par [2])
488 // * [0] = integralBkg;
490 //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
491 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]);
495 //y=a+b*x -> integral = a(max-min)+1/2*b*(max^2-min^2) -> a = (integral-1/2*b*(max^2-min^2))/(max-min)=integral/(max-min)-1/2*b*(max+min)
496 // * [0] = integralBkg;
498 total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
502 //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
503 //+ 1/3*c*(max^3-min^3) ->
504 //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
505 // * [0] = integralBkg;
508 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));
511 total=par[0+firstPar];
514 // Types of Fit Functions for Background:
515 // * 0 = exponential;
517 // * 2 = polynomial 2nd order
518 // * 3 = no background"<<endl;
524 //__________________________________________________________________________
526 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
527 // Main method of the class: performs the fit of the histogram
529 Int_t nFitPars=0; //total function's number of parameters
530 switch (ftypeOfFit4Bkg){
545 Int_t bkgPar = nFitPars-3; //background function's number of parameters
547 cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
550 TString bkgname="funcbkg";
551 TString bkg1name="funcbkg1";
552 TString massname="funcmass";
555 Double_t totInt = fhistoInvMass->Integral("width");
558 Double_t width=fhistoInvMass->GetBinWidth(8);
559 Int_t binleft,binright;
560 if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
561 Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
562 Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
564 if(fMass-fminMass < 0) {
565 cout<<"Left limit of range > mean: change left limit or initial mean value"<<endl;
569 binleft=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
570 binright=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
573 if((fminMass!=minHisto || fmaxMass!= maxHisto) && (fMass-4.*fSigmaSgn-fminMass) < 0){
574 Double_t coeff = (fMass-fminMass)/fSigmaSgn;
575 binleft=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
576 binright=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
577 cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
581 cout<<"Error! Range too little";
584 Int_t fbin=(Int_t)((fminMass-minHisto)/width)+1, lbin=(Int_t)((fmaxMass-minHisto)/width)+1;
585 cout<<"Range: "<<fminMass<<", "<<fmaxMass<<endl;
586 cout<<"Range in bin = "<<fbin<<" --> "<<lbin<<endl;
588 //sidebands integral - first approx (from histo)
589 Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,binleft,"width") + (Double_t)fhistoInvMass->Integral(binright,fNbin,"width");
590 cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<binleft<<"\tbinright = "<<binright<<endl;
591 cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
593 //side bands usando il range utente
594 Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(fbin,binleft,"width") + (Double_t)fhistoInvMass->Integral(binright,lbin,"width");
595 cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<binleft<<"\tbinright = "<<binright<<endl;
596 cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
600 //TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
601 TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
602 cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
604 funcbkg->SetLineColor(2); //red
606 //first fit for bkg: approx bkgint
608 switch (ftypeOfFit4Bkg) {
610 funcbkg->SetParNames("BkgInt","Slope");
611 funcbkg->SetParameters(sideBandsInt,-2.);
614 funcbkg->SetParNames("BkgInt","Slope");
615 funcbkg->SetParameters(sideBandsInt,-100.);
618 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
619 funcbkg->SetParameters(sideBandsInt,-10.,5);
622 if(ftypeOfFit4Sgn==0){
623 funcbkg->SetParNames("Const");
624 funcbkg->SetParameter(0,0.);
625 funcbkg->FixParameter(0,0.);
629 cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
633 cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
634 Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
636 Double_t intbkg1=0,slope1=0,conc1=0;
637 //if only signal and reflection: skip
638 if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
640 fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
642 for(Int_t i=0;i<bkgPar;i++){
643 fFitPars[i]=funcbkg->GetParameter(i);
644 //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
645 fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
646 //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
649 //intbkg1 = funcbkg->GetParameter(0);
650 funcbkg->SetRange(fminMass,fmaxMass);
651 intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
652 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
653 if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
654 cout<<"Primo fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
656 else cout<<"\t\t//"<<endl;
658 ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
660 if (ftypeOfFit4Sgn == 1) {
661 cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
664 cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
666 funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
667 cout<<"Function name = "<<funcbkg1->GetName()<<endl;
669 funcbkg1->SetLineColor(2); //red
671 if(ftypeOfFit4Bkg==2){
672 cout<<"*** Polynomial Fit ***"<<endl;
673 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
674 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
675 // cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
676 // cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
678 if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
680 cout<<"*** No background Fit ***"<<endl;
681 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
682 funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
683 funcbkg1->FixParameter(3,0.);
684 } else{ //expo or linear
685 if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
686 if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
687 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
688 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
691 fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
693 for(Int_t i=0;i<bkgPar;i++){
694 fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
695 //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
696 fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
697 //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
700 intbkg1=funcbkg1->GetParameter(3);
701 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
702 if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
707 for(Int_t i=0;i<3;i++){
708 fFitPars[bkgPar-3+i]=0.;
709 cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
710 fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
711 cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
714 for(Int_t i=0;i<bkgPar-3;i++){
715 fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
716 cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
717 fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
718 cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
724 //sidebands integral - second approx (from fit)
728 if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
729 else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
730 cout<<"------BkgInt(Fit) = "<<bkgInt<<endl;
732 //Signal integral - first approx
734 sgnInt = totInt-bkgInt;
735 cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
737 /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
738 TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
739 cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
740 funcmass->SetLineColor(4); //blue
743 cout<<"\nTOTAL FIT"<<endl;
746 funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
747 funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
748 //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
749 // cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
750 funcmass->FixParameter(0,totInt);
753 funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
754 funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
755 // cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
756 // cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
757 funcmass->FixParameter(0,totInt);
760 funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
761 if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
762 else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
763 funcmass->FixParameter(0,0.);
764 //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
765 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
769 fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
770 cout<<"fit done"<<endl;
771 //reset value of fMass and fSigmaSgn to those found from fit
772 fMass=funcmass->GetParameter(nFitPars-2);
773 fSigmaSgn=funcmass->GetParameter(nFitPars-1);
775 for(Int_t i=0;i<nFitPars;i++){
776 fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
777 fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
778 //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
781 //check: cout parameters
782 for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
783 cout<<i<<"\t"<<fFitPars[i]<<endl;
788 TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
789 TH1F *fhistocopy=new TH1F(*fhistoInvMass);
791 fhistocopy->DrawClone();
792 if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
793 else funcbkg->DrawClone("sames");
794 funcmass->DrawClone("sames");
799 //increase counter of number of fits done
802 if (ftypeOfFit4Sgn == 1) delete funcbkg1;
807 //_________________________________________________________________________
808 Double_t AliHFMassFitter::GetChiSquare() const{
809 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
810 return funcmass->GetChisquare();
815 //_________________________________________________________________________
816 void AliHFMassFitter::GetFitPars(Float_t *vector) const {
817 // Return fit parameters
819 for(Int_t i=0;i<fParsSize;i++){
820 vector[i]=fFitPars[i];
823 //_________________________________________________________________________
825 TH1F* AliHFMassFitter::GetHistoClone() const{
826 TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
829 //_________________________________________________________________________
831 void AliHFMassFitter::WriteHisto(TString path) {
833 cout<<"Use MassFitter method before WriteHisto"<<endl;
837 TH1F* hget=(TH1F*)fhistoInvMass->Clone();
839 TString bkgname = "funcbkg";
841 switch (ftypeOfFit4Bkg){
855 if (ftypeOfFit4Sgn == 1) {
859 TList *hlist=hget->GetListOfFunctions();
861 TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
862 cout<< b->GetParameter(0)<<endl;
863 //cout<<"WRITEHISTO np = "<<np<<"\n"<<bkgname<<endl;
864 bkgname += "FullRange";
865 TF1 *bwrite=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
866 cout<<bwrite->GetName()<<endl;
867 for(Int_t i=0;i<np;i++){
868 //cout<<i<<" di "<<np<<endl;
869 bwrite->SetParName(i,b->GetParName(i));
870 bwrite->SetParameter(i,b->GetParameter(i));
874 path += "HFMassFitterOutput.root";
877 if (fcounter == 1) output = new TFile(path.Data(),"recreate");
878 else output = new TFile(path.Data(),"update");
883 cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
890 //_________________________________________________________________________
892 void AliHFMassFitter::WriteNtuple(TString path) const{
893 TNtuple* nget=(TNtuple*)fntuParam->Clone();
894 path += "HFMassFitterOutput.root";
895 TFile *output = new TFile(path.Data(),"update");
899 cout<<nget->GetName()<<" written in "<<path<<endl;
905 //************ significance
907 //_________________________________________________________________________
909 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
910 // Return signal integral in mean+- n sigma
914 TString bkgname="funcbkg";
915 TString bkg1name="funcbkg1";
916 TString massname="funcmass";
918 TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
919 Double_t mean=0,sigma=1;
920 Double_t intS,intSerr;
922 if(ftypeOfFit4Bkg == 2) { //polynomial
923 mean=funcmass->GetParameter(4); //mean
924 sigma=funcmass->GetParameter(5); //sigma
926 intSerr=fFitPars[27];
927 } else if(ftypeOfFit4Bkg == 3){ //no background
928 mean=funcmass->GetParameter(2); //mean
929 sigma=funcmass->GetParameter(3); //sigma
931 intSerr=fFitPars[15];
932 } else { //expo or linear
933 mean=funcmass->GetParameter(3); //mean
934 sigma=funcmass->GetParameter(4); //sigma
936 intSerr=fFitPars[21];
940 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
941 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
942 Double_t min=mean-nOfSigma*sigma;
943 Double_t max=mean+nOfSigma*sigma;
944 if(ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) signal=funcmass->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
945 else signal=(funcmass->Integral(min,max)-funcbkg->Integral(min,max))/(Double_t)fhistoInvMass->GetBinWidth(2);
946 errsignal=intSerr/intS*signal; // assume relative error is the same as for total integral
950 //_________________________________________________________________________
952 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
953 // Return background integral in mean+- n sigma
956 TString bkgname="funcbkg";
957 TString bkg1name="funcbkg1";
958 TString massname="funcmass";
960 TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
961 Double_t mean=0,sigma=1;
962 Double_t intB,intBerr,err;
964 if(ftypeOfFit4Bkg == 2) { //polynomial
965 mean=funcmass->GetParameter(4); //mean
966 sigma=funcmass->GetParameter(5); //sigma
967 intB=fFitPars[9]-fFitPars[12];
968 intBerr=fFitPars[27];
969 } else if(ftypeOfFit4Bkg == 3){ //no background
970 mean=funcmass->GetParameter(2); //mean
971 sigma=funcmass->GetParameter(3); //sigma
972 if (ftypeOfFit4Sgn == 1){ //reflection
979 } else { //expo or linear
980 mean=funcmass->GetParameter(3); //mean
981 sigma=funcmass->GetParameter(4); //sigma
982 intB=fFitPars[7]-fFitPars[9];
983 intBerr=fFitPars[21];
988 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
989 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
990 Double_t min=mean-nOfSigma*sigma;
991 Double_t max=mean+nOfSigma*sigma;
992 if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
997 background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
998 errbackground=intBerr/intB*background;
1004 //__________________________________________________________________________
1006 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
1007 // Return significance in mean+- n sigma
1009 Double_t signal,errsignal,background,errbackground;
1010 Signal(nOfSigma,signal,errsignal);
1011 Background(nOfSigma,background,errbackground);
1013 significance = signal/TMath::Sqrt(signal+background);
1014 errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1019 //__________________________________________________________________________