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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // AliHFMassFitter for the fit of invariant mass distribution
23 // Author: C.Bianchin, chiara.bianchin@pd.infn.it
24 /////////////////////////////////////////////////////////////
26 #include <Riostream.h>
34 #include <TVirtualPad.h>
36 #include <TVirtualFitter.h>
39 #include <TPaveText.h>
40 #include <TDatabasePDG.h>
42 #include "AliHFMassFitter.h"
46 ClassImp(AliHFMassFitter)
48 //************** constructors
49 AliHFMassFitter::AliHFMassFitter() :
74 // default constructor
76 cout<<"Default constructor:"<<endl;
77 cout<<"Remember to set the Histo, the Type, the FixPar"<<endl;
81 //___________________________________________________________________________
83 AliHFMassFitter::AliHFMassFitter (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
108 // standard constructor
110 fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
111 fhistoInvMass->SetDirectory(0);
114 if(rebin!=1) RebinMass(rebin);
115 else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
117 ftypeOfFit4Bkg=fittypeb;
118 ftypeOfFit4Sgn=fittypes;
119 if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2 && ftypeOfFit4Bkg!=4 && ftypeOfFit4Bkg!=5) fWithBkg=kFALSE;
121 if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
122 else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
125 fFitPars=new Float_t[fParsSize];
127 SetDefaultFixParam();
129 fContourGraph=new TList();
130 fContourGraph->SetOwner();
135 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
137 fhistoInvMass(mfit.fhistoInvMass),
138 fminMass(mfit.fminMass),
139 fmaxMass(mfit.fmaxMass),
140 fminBinMass(mfit.fminBinMass),
141 fmaxBinMass(mfit.fmaxBinMass),
143 fParsSize(mfit.fParsSize),
144 fNFinalPars(mfit.fNFinalPars),
146 fWithBkg(mfit.fWithBkg),
147 ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
148 ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
149 ffactor(mfit.ffactor),
150 fntuParam(mfit.fntuParam),
152 fSigmaSgn(mfit.fSigmaSgn),
153 fSideBands(mfit.fSideBands),
155 fSideBandl(mfit.fSideBandl),
156 fSideBandr(mfit.fSideBandr),
157 fcounter(mfit.fcounter),
158 fContourGraph(mfit.fContourGraph)
162 if(mfit.fParsSize > 0){
163 fFitPars=new Float_t[fParsSize];
164 fFixPar=new Bool_t[fNFinalPars];
165 memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
166 memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Bool_t));
171 //_________________________________________________________________________
173 AliHFMassFitter::~AliHFMassFitter() {
177 cout<<"AliHFMassFitter destructor called"<<endl;
179 delete fhistoInvMass;
190 //_________________________________________________________________________
192 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
194 //assignment operator
196 if(&mfit == this) return *this;
197 fhistoInvMass= mfit.fhistoInvMass;
198 fminMass= mfit.fminMass;
199 fmaxMass= mfit.fmaxMass;
201 fParsSize= mfit.fParsSize;
202 fWithBkg= mfit.fWithBkg;
203 ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
204 ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
205 ffactor= mfit.ffactor;
206 fntuParam= mfit.fntuParam;
208 fSigmaSgn= mfit.fSigmaSgn;
209 fSideBands= mfit.fSideBands;
210 fSideBandl= mfit.fSideBandl;
211 fSideBandr= mfit.fSideBandr;
212 fcounter= mfit.fcounter;
213 fContourGraph= mfit.fContourGraph;
215 if(mfit.fParsSize > 0){
218 fFitPars=new Float_t[fParsSize];
219 memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
223 fFixPar=new Bool_t[fNFinalPars];
224 memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Float_t));
230 //************ tools & settings
232 //__________________________________________________________________________
234 void AliHFMassFitter::ComputeNFinalPars() {
236 //compute the number of parameters of the total (signal+bgk) function
237 cout<<"Info:ComputeNFinalPars... ";
238 switch (ftypeOfFit4Bkg) {//npar background func
258 cout<<"Error in computing fNFinalPars: check ftypeOfFit4Bkg"<<endl;
262 fNFinalPars+=3; //gaussian signal
263 cout<<": "<<fNFinalPars<<endl;
265 //__________________________________________________________________________
267 void AliHFMassFitter::ComputeParSize() {
269 //compute the size of the parameter array and set the data member
271 switch (ftypeOfFit4Bkg) {//npar background func
291 cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
295 fParsSize += 3; // npar refl
296 fParsSize += 3; // npar signal gaus
298 fParsSize*=2; // add errors
299 cout<<"Parameters array size "<<fParsSize<<endl;
302 //___________________________________________________________________________
303 void AliHFMassFitter::SetDefaultFixParam(){
305 //Set default values for fFixPar (only total integral fixed)
308 fFixPar=new Bool_t[fNFinalPars];
310 fFixPar[0]=kTRUE; //default: IntTot fixed
311 cout<<"Parameter 0 is fixed"<<endl;
312 for(Int_t i=1;i<fNFinalPars;i++){
318 //___________________________________________________________________________
319 Bool_t AliHFMassFitter::SetFixThisParam(Int_t thispar,Bool_t fixpar){
321 //set the value (kFALSE or kTRUE) of one element of fFixPar
322 //return kFALSE if something wrong
324 if(thispar>=fNFinalPars) {
325 cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
329 cout<<"Initializing fFixPar...";
330 SetDefaultFixParam();
331 cout<<" done."<<endl;
334 fFixPar[thispar]=fixpar;
335 if(fixpar)cout<<"Parameter "<<thispar<<" is now fixed"<<endl;
336 else cout<<"Parameter "<<thispar<<" is now free"<<endl;
340 //___________________________________________________________________________
341 Bool_t AliHFMassFitter::GetFixThisParam(Int_t thispar)const{
342 //return the value of fFixPar[thispar]
343 if(thispar>=fNFinalPars) {
344 cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
348 cout<<"Error! Parameters to be fixed still not set"<<endl;
352 return fFixPar[thispar];
356 //___________________________________________________________________________
357 void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
359 fhistoInvMass = new TH1F(*histoToFit);
360 fhistoInvMass->SetDirectory(0);
361 //cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
364 //___________________________________________________________________________
366 void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
368 //set the type of fit to perform for signal and background
370 ftypeOfFit4Bkg = fittypeb;
371 ftypeOfFit4Sgn = fittypes;
374 fFitPars = new Float_t[fParsSize];
376 SetDefaultFixParam();
379 //___________________________________________________________________________
381 void AliHFMassFitter::Reset() {
383 //delete the histogram and reset the mean and sigma to default
385 cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
388 cout<<"Reset "<<fhistoInvMass<<endl;
389 delete fhistoInvMass;
392 //_________________________________________________________________________
394 void AliHFMassFitter::InitNtuParam(TString ntuname) {
396 // Create ntuple to keep fit parameters
399 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");
403 //_________________________________________________________________________
405 void AliHFMassFitter::FillNtuParam() {
406 // Fill ntuple with fit parameters
410 if (ftypeOfFit4Bkg==2) {
411 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
412 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
413 fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
414 fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
415 fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
416 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
417 fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
418 fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
419 fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
420 fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
421 fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
422 fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
423 fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
424 fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
425 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
427 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
428 fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
429 fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
430 fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
431 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
432 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
433 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
434 fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
435 fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
436 fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
437 fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
438 fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
439 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
440 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
441 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
445 if(ftypeOfFit4Bkg==3){
446 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
447 fntuParam->SetBranchAddress("slope1",¬hing);
448 fntuParam->SetBranchAddress("conc1",¬hing);
449 fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
450 fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
451 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
452 fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
453 fntuParam->SetBranchAddress("slope2",¬hing);
454 fntuParam->SetBranchAddress("conc2",¬hing);
455 fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
456 fntuParam->SetBranchAddress("slope3",¬hing);
457 fntuParam->SetBranchAddress("conc3",¬hing);
458 fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
459 fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
460 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
462 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
463 fntuParam->SetBranchAddress("slope1Err",¬hing);
464 fntuParam->SetBranchAddress("conc1Err",¬hing);
465 fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
466 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
467 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
468 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
469 fntuParam->SetBranchAddress("slope2Err",¬hing);
470 fntuParam->SetBranchAddress("conc2Err",¬hing);
471 fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
472 fntuParam->SetBranchAddress("slope3Err",¬hing);
473 fntuParam->SetBranchAddress("conc3Err",¬hing);
474 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
475 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
476 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
480 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
481 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
482 fntuParam->SetBranchAddress("conc1",¬hing);
483 fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
484 fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
485 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
486 fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
487 fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
488 fntuParam->SetBranchAddress("conc2",¬hing);
489 fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
490 fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
491 fntuParam->SetBranchAddress("conc3",¬hing);
492 fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
493 fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
494 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
496 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
497 fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
498 fntuParam->SetBranchAddress("conc1Err",¬hing);
499 fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
500 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
501 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
502 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
503 fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
504 fntuParam->SetBranchAddress("conc2Err",¬hing);
505 fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
506 fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
507 fntuParam->SetBranchAddress("conc3Err",¬hing);
508 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
509 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
510 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
514 fntuParam->TTree::Fill();
517 //_________________________________________________________________________
519 TNtuple* AliHFMassFitter::NtuParamOneShot(TString ntuname){
520 // Create, fill and return ntuple with fit parameters
522 InitNtuParam(ntuname.Data());
526 //_________________________________________________________________________
528 void AliHFMassFitter::RebinMass(Int_t bingroup){
529 // Rebin invariant mass histogram
532 cout<<"Histogram not set"<<endl;
535 Int_t nbinshisto=fhistoInvMass->GetNbinsX();
537 cout<<"Error! Cannot group "<<bingroup<<" bins\n";
539 cout<<"Kept original number of bins: "<<fNbin<<endl;
542 while(nbinshisto%bingroup != 0) {
545 cout<<"Group "<<bingroup<<" bins"<<endl;
546 fhistoInvMass->Rebin(bingroup);
547 fNbin = fhistoInvMass->GetNbinsX();
548 cout<<"New number of bins: "<<fNbin<<endl;
555 //___________________________________________________________________________
557 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
558 // Fit function for signal+background
561 //exponential or linear fit
563 // par[0] = tot integral
565 // par[2] = gaussian integral
566 // par[3] = gaussian mean
567 // par[4] = gaussian sigma
569 Double_t total,bkg=0,sgn=0;
571 if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
572 if(ftypeOfFit4Sgn == 0) {
574 Double_t parbkg[2] = {par[0]-par[2], par[1]};
575 bkg = FitFunction4Bkg(x,parbkg);
577 if(ftypeOfFit4Sgn == 1) {
578 Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
579 bkg = FitFunction4Bkg(x,parbkg);
582 sgn = FitFunction4Sgn(x,&par[2]);
588 // par[0] = tot integral
591 // par[3] = gaussian integral
592 // par[4] = gaussian mean
593 // par[5] = gaussian sigma
595 if (ftypeOfFit4Bkg==2) {
597 if(ftypeOfFit4Sgn == 0) {
599 Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
600 bkg = FitFunction4Bkg(x,parbkg);
602 if(ftypeOfFit4Sgn == 1) {
604 Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
605 bkg = FitFunction4Bkg(x,parbkg);
608 sgn = FitFunction4Sgn(x,&par[3]);
611 if (ftypeOfFit4Bkg==3) {
613 if(ftypeOfFit4Sgn == 0) {
614 bkg=FitFunction4Bkg(x,par);
615 sgn=FitFunction4Sgn(x,&par[1]);
617 if(ftypeOfFit4Sgn == 1) {
618 Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
619 bkg=FitFunction4Bkg(x,parbkg);
620 sgn=FitFunction4Sgn(x,&par[1]);
626 // par[0] = tot integral
628 // par[2] = gaussian integral
629 // par[3] = gaussian mean
630 // par[4] = gaussian sigma
632 if (ftypeOfFit4Bkg==4) {
634 if(ftypeOfFit4Sgn == 0) {
636 Double_t parbkg[2] = {par[0]-par[2], par[1]};
637 bkg = FitFunction4Bkg(x,parbkg);
639 if(ftypeOfFit4Sgn == 1) {
641 Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-par[2], par[1]};
642 bkg = FitFunction4Bkg(x,parbkg);
644 sgn = FitFunction4Sgn(x,&par[2]);
648 //Power and exponential fit
650 // par[0] = tot integral
653 // par[3] = gaussian integral
654 // par[4] = gaussian mean
655 // par[5] = gaussian sigma
657 if (ftypeOfFit4Bkg==5) {
659 if(ftypeOfFit4Sgn == 0) {
660 Double_t parbkg[3] = {par[0]-par[3],par[1],par[2]};
661 bkg = FitFunction4Bkg(x,parbkg);
663 if(ftypeOfFit4Sgn == 1) {
664 Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-par[3], par[1], par[2]};
665 bkg = FitFunction4Bkg(x,parbkg);
667 sgn = FitFunction4Sgn(x,&par[3]);
675 //_________________________________________________________________________
676 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
677 // Fit function for the signal
679 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
681 // * [0] = integralSgn
684 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
686 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]);
690 //__________________________________________________________________________
692 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
693 // Fit function for the background
695 Double_t maxDeltaM = 4.*fSigmaSgn;
696 if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
701 Double_t gaus2=0,total=-1;
702 if(ftypeOfFit4Sgn == 1){
704 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
706 // * [0] = integralSgn
709 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
710 gaus2 = FitFunction4Sgn(x,par);
713 switch (ftypeOfFit4Bkg){
716 //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
717 //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
718 //as integralTot- integralGaus (=par [2])
720 // * [0] = integralBkg;
722 //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
723 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]);
727 //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)
728 // * [0] = integralBkg;
730 total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
734 //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
735 //+ 1/3*c*(max^3-min^3) ->
736 //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
737 // * [0] = integralBkg;
740 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));
743 total=par[0+firstPar];
747 //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
749 //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
750 // * [0] = integralBkg;
752 // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
754 Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
756 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]);
760 //power function wit exponential
761 //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
763 Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
765 total = par[1+firstPar]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2+firstPar]*(x[0]-mpi));
769 // Types of Fit Functions for Background:
770 // * 0 = exponential;
772 // * 2 = polynomial 2nd order
773 // * 3 = no background"<<endl;
774 // * 4 = Power function
775 // * 5 = Power function with exponential
781 //__________________________________________________________________________
782 Bool_t AliHFMassFitter::SideBandsBounds(){
784 //determines the ranges of the side bands
786 if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
787 Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
788 Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin+1);
790 Double_t sidebandldouble,sidebandrdouble;
791 Bool_t leftok=kFALSE, rightok=kFALSE;
793 if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
794 cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
798 //histo limit = fit function limit
799 if((TMath::Abs(fminMass-minHisto) < 1e6 || TMath::Abs(fmaxMass - maxHisto) < 1e6) && (fMass-4.*fSigmaSgn-fminMass) < 1e6){
800 Double_t coeff = (fMass-fminMass)/fSigmaSgn;
801 sidebandldouble=(fMass-0.5*coeff*fSigmaSgn);
802 sidebandrdouble=(fMass+0.5*coeff*fSigmaSgn);
803 cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
804 if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
806 cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
808 //set binleft and right without considering SetRangeFit- anyway no bkg!
809 sidebandldouble=(fMass-4.*fSigmaSgn);
810 sidebandrdouble=(fMass+4.*fSigmaSgn);
814 sidebandldouble=(fMass-4.*fSigmaSgn);
815 sidebandrdouble=(fMass+4.*fSigmaSgn);
818 cout<<"Left side band ";
821 //calculate bin corresponding to fSideBandl
822 fSideBandl=fhistoInvMass->FindBin(sidebandldouble);
823 if (sidebandldouble >= fhistoInvMass->GetBinCenter(fSideBandl)) fSideBandl++;
824 sidebandldouble=fhistoInvMass->GetBinLowEdge(fSideBandl);
826 if(TMath::Abs(tmp-sidebandldouble) > 1e-6){
827 cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
830 cout<<sidebandldouble<<" (bin "<<fSideBandl<<")"<<endl;
832 cout<<"Right side band ";
834 //calculate bin corresponding to fSideBandr
835 fSideBandr=fhistoInvMass->FindBin(sidebandrdouble);
836 if (sidebandrdouble < fhistoInvMass->GetBinCenter(fSideBandr)) fSideBandr--;
837 sidebandrdouble=fhistoInvMass->GetBinLowEdge(fSideBandr+1);
839 if(TMath::Abs(tmp-sidebandrdouble) > 1e-6){
840 cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
843 cout<<sidebandrdouble<<" (bin "<<fSideBandr<<")"<<endl;
844 if (fSideBandl==0 || fSideBandr==fNbin) {
845 cout<<"Error! Range too little";
851 //__________________________________________________________________________
853 void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
855 // get the range of the side bands
857 if (fSideBandl==0 && fSideBandr==0){
858 cout<<"Use MassFitter method first"<<endl;
865 //__________________________________________________________________________
866 Bool_t AliHFMassFitter::CheckRangeFit(){
867 //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
869 if (!fhistoInvMass) {
870 cout<<"No histogram to fit! SetHisto(TH1F* h) before! "<<endl;
873 Bool_t leftok=kFALSE, rightok=kFALSE;
874 Int_t nbins=fhistoInvMass->GetNbinsX();
875 Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins+1);
877 //check if limits are inside histogram range
879 if( fminMass-minhisto < 0. ) {
880 cout<<"Out of histogram left bound! Setting to "<<minhisto<<endl;
883 if( fmaxMass-maxhisto > 0. ) {
884 cout<<"Out of histogram right bound! Setting to"<<maxhisto<<endl;
890 //calculate bin corresponding to fminMass
891 fminBinMass=fhistoInvMass->FindBin(fminMass);
892 if (fminMass >= fhistoInvMass->GetBinCenter(fminBinMass)) fminBinMass++;
893 fminMass=fhistoInvMass->GetBinLowEdge(fminBinMass);
894 if(TMath::Abs(tmp-fminMass) > 1e-6){
895 cout<<"Left bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fminMass<<endl;
900 //calculate bin corresponding to fmaxMass
901 fmaxBinMass=fhistoInvMass->FindBin(fmaxMass);
902 if (fmaxMass < fhistoInvMass->GetBinCenter(fmaxBinMass)) fmaxBinMass--;
903 fmaxMass=fhistoInvMass->GetBinLowEdge(fmaxBinMass+1);
904 if(TMath::Abs(tmp-fmaxMass) > 1e-6){
905 cout<<"Right bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fmaxMass<<endl;
909 return (leftok && rightok);
913 //__________________________________________________________________________
915 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
916 // Main method of the class: performs the fit of the histogram
918 //Set default fitter Minuit in order to use gMinuit in the contour plots
919 TVirtualFitter::SetDefaultFitter("Minuit");
922 Bool_t isBkgOnly=kFALSE;
924 Int_t fit1status=RefitWithBkgOnly(kFALSE);
926 Int_t checkinnsigma=4;
927 Double_t range[2]={fMass-checkinnsigma*fSigmaSgn,fMass+checkinnsigma*fSigmaSgn};
928 TF1* func=GetHistoClone()->GetFunction("funcbkgonly");
929 Double_t intUnderFunc=func->Integral(range[0],range[1]);
930 Double_t intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(range[0]),fhistoInvMass->FindBin(range[1]),"width");
931 cout<<"Pick zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
932 Double_t diffUnderPick=(intUnderHisto-intUnderFunc);
933 intUnderFunc=func->Integral(fminMass,fminMass+checkinnsigma*fSigmaSgn);
934 intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fminMass+checkinnsigma*fSigmaSgn),"width");
935 cout<<"Band (l) zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
936 Double_t diffUnderBands=(intUnderHisto-intUnderFunc);
937 Double_t relDiff=diffUnderPick/diffUnderBands;
938 cout<<"Relative difference = "<<relDiff<<endl;
939 if(TMath::Abs(relDiff) < 1) isBkgOnly=kTRUE;
941 cout<<"Relative difference = "<<relDiff<<": I suppose there is some signal, continue with total fit!"<<endl;
945 cout<<"INFO!! The histogram contains only background"<<endl;
948 //increase counter of number of fits done
954 Int_t bkgPar = fNFinalPars-3; //background function's number of parameters
956 cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
959 TString listname="contourplot";
962 fContourGraph=new TList();
963 fContourGraph->SetOwner();
966 fContourGraph->SetName(listname);
970 TString bkgname="funcbkg";
971 TString bkg1name="funcbkg1";
972 TString massname="funcmass";
975 Double_t totInt = fhistoInvMass->Integral(fminBinMass,fmaxBinMass, "width");
976 //cout<<"Here tot integral is = "<<totInt<<"; integral in whole range is "<<fhistoInvMass->Integral("width")<<endl;
978 Double_t width=fhistoInvMass->GetBinWidth(8);
979 //cout<<"fNbin = "<<fNbin<<endl;
980 if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
982 Bool_t ok=SideBandsBounds();
983 if(!ok) return kFALSE;
985 //sidebands integral - first approx (from histo)
986 Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
987 cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
988 cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
989 if (sideBandsInt<=0) {
990 cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
997 TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
998 cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
1000 funcbkg->SetLineColor(2); //red
1002 //first fit for bkg: approx bkgint
1004 switch (ftypeOfFit4Bkg) {
1006 funcbkg->SetParNames("BkgInt","Slope");
1007 funcbkg->SetParameters(sideBandsInt,-2.);
1010 funcbkg->SetParNames("BkgInt","Slope");
1011 funcbkg->SetParameters(sideBandsInt,-100.);
1014 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1015 funcbkg->SetParameters(sideBandsInt,-10.,5);
1018 if(ftypeOfFit4Sgn==0){
1019 funcbkg->SetParNames("Const");
1020 funcbkg->SetParameter(0,0.);
1021 funcbkg->FixParameter(0,0.);
1025 funcbkg->SetParNames("BkgInt","Coef2");
1026 funcbkg->SetParameters(sideBandsInt,0.5);
1029 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1030 funcbkg->SetParameters(sideBandsInt, -10., 5.);
1033 cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1037 cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
1038 Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
1040 Double_t intbkg1=0,slope1=0,conc1=0;
1041 //if only signal and reflection: skip
1042 if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
1044 fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
1046 for(Int_t i=0;i<bkgPar;i++){
1047 fFitPars[i]=funcbkg->GetParameter(i);
1048 //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
1049 fFitPars[fNFinalPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
1050 //cout<<fNFinalPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
1052 fSideBands = kFALSE;
1053 //intbkg1 = funcbkg->GetParameter(0);
1055 intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
1056 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
1057 if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
1058 if(ftypeOfFit4Bkg==5) conc1 = funcbkg->GetParameter(2);
1061 //cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
1063 else cout<<"\t\t//"<<endl;
1065 ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
1067 if (ftypeOfFit4Sgn == 1) {
1068 cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
1071 //cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
1073 funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
1074 cout<<"Function name = "<<funcbkg1->GetName()<<endl;
1076 funcbkg1->SetLineColor(2); //red
1078 switch (ftypeOfFit4Bkg) {
1081 cout<<"*** Exponential Fit ***"<<endl;
1082 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
1083 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1088 cout<<"*** Linear Fit ***"<<endl;
1089 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
1090 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1095 cout<<"*** Polynomial Fit ***"<<endl;
1096 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
1097 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
1101 //no background: gaus sign+ gaus broadened
1103 cout<<"*** No background Fit ***"<<endl;
1104 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
1105 funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
1106 funcbkg1->FixParameter(3,0.);
1111 cout<<"*** Power function Fit ***"<<endl;
1112 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef2");
1113 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1118 cout<<"*** Power function conv. with exponential Fit ***"<<endl;
1119 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
1120 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
1124 //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
1125 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
1127 Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
1129 cout<<"Minuit returned "<<status<<endl;
1133 for(Int_t i=0;i<bkgPar;i++){
1134 fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
1135 //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
1136 fFitPars[fNFinalPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
1137 //cout<<"\t"<<fNFinalPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
1140 intbkg1=funcbkg1->GetParameter(3);
1141 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
1142 if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
1143 if(ftypeOfFit4Bkg==5) conc1 = funcbkg1->GetParameter(5);
1149 for(Int_t i=0;i<3;i++){
1150 fFitPars[bkgPar-3+i]=0.;
1151 cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
1152 fFitPars[fNFinalPars+3*bkgPar-6+i]= 0.;
1153 cout<<fNFinalPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
1156 for(Int_t i=0;i<bkgPar-3;i++){
1157 fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
1158 cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
1159 fFitPars[fNFinalPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
1160 cout<<fNFinalPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
1166 //sidebands integral - second approx (from fit)
1167 fSideBands = kFALSE;
1169 //cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
1170 if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
1171 else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
1172 //cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
1174 //Signal integral - first approx
1176 sgnInt = totInt-bkgInt;
1177 //cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
1179 cout<<"Setting sgnInt = - sgnInt"<<endl;
1182 /*Fit All Mass distribution with exponential + gaussian (+gaussian braodened) */
1183 TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4MassDistr");
1184 cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
1185 funcmass->SetLineColor(4); //blue
1188 cout<<"\nTOTAL FIT"<<endl;
1191 funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
1192 funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
1194 //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1195 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1197 funcmass->FixParameter(0,totInt);
1200 funcmass->FixParameter(1,slope1);
1203 funcmass->FixParameter(2,sgnInt);
1206 funcmass->FixParameter(3,fMass);
1209 funcmass->FixParameter(4,fSigmaSgn);
1212 if (fNFinalPars==6){
1213 funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
1214 funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
1216 //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1217 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1218 if(fFixPar[0])funcmass->FixParameter(0,totInt);
1219 if(fFixPar[1])funcmass->FixParameter(1,slope1);
1220 if(fFixPar[2])funcmass->FixParameter(2,conc1);
1221 if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
1222 if(fFixPar[4])funcmass->FixParameter(4,fMass);
1223 if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
1225 //funcmass->FixParameter(2,sgnInt);
1228 funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
1229 if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
1230 else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
1231 if(fFixPar[0]) funcmass->FixParameter(0,0.);
1232 if(fFixPar[1])funcmass->FixParameter(1,sgnInt);
1233 if(fFixPar[2])funcmass->FixParameter(2,fMass);
1234 if(fFixPar[3])funcmass->FixParameter(3,fSigmaSgn);
1235 //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1236 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1242 status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
1244 cout<<"Minuit returned "<<status<<endl;
1248 cout<<"fit done"<<endl;
1249 //reset value of fMass and fSigmaSgn to those found from fit
1250 fMass=funcmass->GetParameter(fNFinalPars-2);
1251 fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
1253 for(Int_t i=0;i<fNFinalPars;i++){
1254 fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1255 fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1256 //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1259 //check: cout parameters
1260 for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
1261 cout<<i<<"\t"<<fFitPars[i]<<endl;
1265 if(funcmass->GetParameter(fNFinalPars-1) <0 || funcmass->GetParameter(fNFinalPars-2) <0 || funcmass->GetParameter(fNFinalPars-3) <0 ) {
1266 cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1270 //increase counter of number of fits done
1276 for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
1278 for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
1279 cout<<"Par "<<kpar<<" and "<<jpar<<endl;
1281 // produce 2 contours per couple of parameters
1282 TGraph* cont[2] = {0x0, 0x0};
1283 const Double_t errDef[2] = {1., 4.};
1284 for (Int_t i=0; i<2; i++) {
1285 gMinuit->SetErrorDef(errDef[i]);
1286 cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
1287 cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
1290 if(!cont[0] || !cont[1]){
1291 cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
1295 // set graph titles and add them to the list
1296 TString title = "Contour plot";
1297 TString titleX = funcmass->GetParName(kpar);
1298 TString titleY = funcmass->GetParName(jpar);
1299 for (Int_t i=0; i<2; i++) {
1300 cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
1301 cont[i]->SetTitle(title);
1302 cont[i]->GetXaxis()->SetTitle(titleX);
1303 cont[i]->GetYaxis()->SetTitle(titleY);
1304 cont[i]->GetYaxis()->SetLabelSize(0.033);
1305 cont[i]->GetYaxis()->SetTitleSize(0.033);
1306 cont[i]->GetYaxis()->SetTitleOffset(1.67);
1308 fContourGraph->Add(cont[i]);
1312 TString cvname = Form("c%d%d", kpar, jpar);
1313 TCanvas *c4=new TCanvas(cvname,cvname,600,600);
1315 cont[1]->SetFillColor(38);
1316 cont[1]->Draw("alf");
1317 cont[0]->SetFillColor(9);
1318 cont[0]->Draw("lf");
1326 if (ftypeOfFit4Sgn == 1) {
1332 AddFunctionsToHisto();
1333 if (draw) DrawFit();
1339 //______________________________________________________________________________
1341 Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
1343 //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
1344 //If you want to change the backgroud function or range use SetType or SetRangeFit before
1346 TString bkgname="funcbkgonly";
1347 fSideBands = kFALSE;
1349 TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1351 funcbkg->SetLineColor(kBlue+3); //dark blue
1353 Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
1355 switch (ftypeOfFit4Bkg) {
1357 funcbkg->SetParNames("BkgInt","Slope");
1358 funcbkg->SetParameters(integral,-2.);
1361 funcbkg->SetParNames("BkgInt","Slope");
1362 funcbkg->SetParameters(integral,-100.);
1365 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1366 funcbkg->SetParameters(integral,-10.,5);
1369 cout<<"Warning! This choice does not make a lot of sense..."<<endl;
1370 if(ftypeOfFit4Sgn==0){
1371 funcbkg->SetParNames("Const");
1372 funcbkg->SetParameter(0,0.);
1373 funcbkg->FixParameter(0,0.);
1377 funcbkg->SetParNames("BkgInt","Coef1");
1378 funcbkg->SetParameters(integral,0.5);
1381 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1382 funcbkg->SetParameters(integral,-10.,5.);
1385 cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1391 Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
1393 cout<<"Minuit returned "<<status<<endl;
1396 AddFunctionsToHisto();
1403 //_________________________________________________________________________
1404 Double_t AliHFMassFitter::GetChiSquare() const{
1406 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1408 cout<<"funcmass not found"<<endl;
1411 return funcmass->GetChisquare();
1414 //_________________________________________________________________________
1415 Double_t AliHFMassFitter::GetReducedChiSquare() const{
1416 //Get reduced Chi^2 method
1417 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1419 cout<<"funcmass not found"<<endl;
1422 return funcmass->GetChisquare()/funcmass->GetNDF();
1427 //_________________________________________________________________________
1428 void AliHFMassFitter::GetFitPars(Float_t *vector) const {
1429 // Return fit parameters
1431 for(Int_t i=0;i<fParsSize;i++){
1432 vector[i]=fFitPars[i];
1437 //_________________________________________________________________________
1438 void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1440 //gives the integral of signal obtained from fit parameters
1441 if(!valuewitherror) {
1442 printf("AliHFMassFitter::IntS: got a null pointer\n");
1446 Int_t index=fParsSize/2 - 3;
1447 valuewitherror[0]=fFitPars[index];
1448 index=fParsSize - 3;
1449 valuewitherror[1]=fFitPars[index];
1453 //_________________________________________________________________________
1454 void AliHFMassFitter::AddFunctionsToHisto(){
1456 //Add the background function in the complete range to the list of functions attached to the histogram
1458 //cout<<"AddFunctionsToHisto called"<<endl;
1459 TString bkgname = "funcbkg";
1461 Bool_t done1=kFALSE,done2=kFALSE;
1463 TString bkgnamesave=bkgname;
1464 TString testname=bkgname;
1465 testname += "FullRange";
1466 TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1471 testname="funcbkgonly";
1472 testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1479 cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1483 TList *hlist=fhistoInvMass->GetListOfFunctions();
1487 TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1489 cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1491 bonly->SetLineColor(kBlue+3);
1492 hlist->Add((TF1*)bonly->Clone());
1499 TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1501 cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1505 bkgname += "FullRange";
1506 TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1507 //cout<<bfullrange->GetName()<<endl;
1508 for(Int_t i=0;i<fNFinalPars-3;i++){
1509 bfullrange->SetParName(i,b->GetParName(i));
1510 bfullrange->SetParameter(i,b->GetParameter(i));
1511 bfullrange->SetParError(i,b->GetParError(i));
1513 bfullrange->SetLineStyle(4);
1514 bfullrange->SetLineColor(14);
1516 bkgnamesave += "Recalc";
1518 TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1520 TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1523 cout<<"funcmass doesn't exist "<<endl;
1527 //intBkg=intTot-intS
1528 blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
1529 blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
1530 if (fNFinalPars>=5) {
1531 blastpar->SetParameter(1,mass->GetParameter(1));
1532 blastpar->SetParError(1,mass->GetParError(1));
1534 if (fNFinalPars==6) {
1535 blastpar->SetParameter(2,mass->GetParameter(2));
1536 blastpar->SetParError(2,mass->GetParError(2));
1539 blastpar->SetLineStyle(1);
1540 blastpar->SetLineColor(2);
1542 hlist->Add((TF1*)bfullrange->Clone());
1543 hlist->Add((TF1*)blastpar->Clone());
1554 //_________________________________________________________________________
1556 TH1F* AliHFMassFitter::GetHistoClone() const{
1558 TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1561 //_________________________________________________________________________
1563 void AliHFMassFitter::WriteHisto(TString path) const {
1565 //Write the histogram in the default file HFMassFitterOutput.root
1567 if (fcounter == 0) {
1568 cout<<"Use MassFitter method before WriteHisto"<<endl;
1571 TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1573 path += "HFMassFitterOutput.root";
1576 if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1577 else output = new TFile(path.Data(),"update");
1580 fContourGraph->Write();
1585 cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1591 //_________________________________________________________________________
1593 void AliHFMassFitter::WriteNtuple(TString path) const{
1594 //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1595 path += "HFMassFitterOutput.root";
1596 TFile *output = new TFile(path.Data(),"update");
1601 //cout<<nget->GetName()<<" written in "<<path<<endl;
1602 cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1613 //_________________________________________________________________________
1614 void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1616 //write the canvas in a root file
1618 gStyle->SetOptStat(0);
1619 gStyle->SetCanvasColor(0);
1620 gStyle->SetFrameFillColor(0);
1624 switch (ftypeOfFit4Bkg){
1641 type="PowExp"; //3+3
1645 TString filename=Form("%sMassFit.root",type.Data());
1646 filename.Prepend(userIDstring);
1647 path.Append(filename);
1649 TFile* outputcv=new TFile(path.Data(),"update");
1651 TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1652 c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1653 if(draw)c->DrawClone();
1659 //_________________________________________________________________________
1661 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1662 //return a TVirtualPad with the fitted histograms and info
1664 TString cvtitle="fit of ";
1665 cvtitle+=fhistoInvMass->GetName();
1669 TCanvas *c=new TCanvas(cvname,cvtitle);
1670 PlotFit(c->cd(),nsigma,writeFitInfo);
1673 //_________________________________________________________________________
1675 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1676 //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1677 gStyle->SetOptStat(0);
1678 gStyle->SetCanvasColor(0);
1679 gStyle->SetFrameFillColor(0);
1681 cout<<"nsigma = "<<nsigma<<endl;
1682 cout<<"Verbosity = "<<writeFitInfo<<endl;
1684 TH1F* hdraw=GetHistoClone();
1686 if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1687 cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1691 if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1692 cout<<"Drawing background fit only"<<endl;
1693 hdraw->SetMinimum(0);
1694 hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1696 hdraw->SetMarkerStyle(20);
1697 hdraw->DrawClone("PE");
1698 hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1700 if(writeFitInfo > 0){
1701 TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1702 pinfo->SetBorderSize(0);
1703 pinfo->SetFillStyle(0);
1704 TF1* f=hdraw->GetFunction("funcbkgonly");
1705 for (Int_t i=0;i<fNFinalPars-3;i++){
1706 pinfo->SetTextColor(kBlue+3);
1707 TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1708 pinfo->AddText(str);
1711 pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1721 hdraw->SetMinimum(0);
1722 hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1724 hdraw->SetMarkerStyle(20);
1725 hdraw->DrawClone("PE");
1726 // if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1727 // if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1728 if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1730 if(writeFitInfo > 0){
1731 TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1732 TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1733 pinfob->SetBorderSize(0);
1734 pinfob->SetFillStyle(0);
1735 pinfom->SetBorderSize(0);
1736 pinfom->SetFillStyle(0);
1737 TF1* ff=fhistoInvMass->GetFunction("funcmass");
1739 for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1740 pinfom->SetTextColor(kBlue);
1741 TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1742 if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1745 pinfom->DrawClone();
1747 TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1748 pinfo2->SetBorderSize(0);
1749 pinfo2->SetFillStyle(0);
1751 Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1753 Significance(nsigma,signif,errsignif);
1754 Signal(nsigma,signal,errsignal);
1755 Background(nsigma,bkg, errbkg);
1757 Significance(1.828,1.892,signif,errsignif);
1758 Signal(1.828,1.892,signal,errsignal);
1759 Background(1.828,1.892,bkg, errbkg);
1761 TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1762 pinfo2->AddText(str);
1763 str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1764 pinfo2->AddText(str);
1765 str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1766 pinfo2->AddText(str);
1767 str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
1768 pinfo2->AddText(str);
1773 if(writeFitInfo > 1){
1774 for (Int_t i=0;i<fNFinalPars-3;i++){
1775 pinfob->SetTextColor(kRed);
1776 str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1777 pinfob->AddText(str);
1780 pinfob->DrawClone();
1788 //_________________________________________________________________________
1790 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1791 //draws histogram together with fit functions with default nice colors in user canvas
1792 PlotFit(pd,nsigma,writeFitInfo);
1797 //_________________________________________________________________________
1798 void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1800 //draws histogram together with fit functions with default nice colors
1801 gStyle->SetOptStat(0);
1802 gStyle->SetCanvasColor(0);
1803 gStyle->SetFrameFillColor(0);
1806 TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1812 //_________________________________________________________________________
1814 void AliHFMassFitter::PrintParTitles() const{
1816 //prints to screen the parameters names
1818 TF1 *f=fhistoInvMass->GetFunction("funcmass");
1820 cout<<"Fit function not found"<<endl;
1824 cout<<"Parameter Titles \n";
1825 for(Int_t i=0;i<fNFinalPars;i++){
1826 cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1832 //************ significance
1834 //_________________________________________________________________________
1836 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1837 // Return signal integral in mean+- n sigma
1840 cout<<"Use MassFitter method before Signal"<<endl;
1844 Double_t min=fMass-nOfSigma*fSigmaSgn;
1845 Double_t max=fMass+nOfSigma*fSigmaSgn;
1847 Signal(min,max,signal,errsignal);
1854 //_________________________________________________________________________
1856 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1858 // Return signal integral in a range
1861 cout<<"Use MassFitter method before Signal"<<endl;
1866 TString bkgname="funcbkgRecalc";
1867 TString bkg1name="funcbkg1Recalc";
1868 TString massname="funcmass";
1872 TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1874 cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1878 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1879 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1882 cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1886 Int_t np=fNFinalPars-3;
1888 Double_t intS,intSerr;
1890 //relative error evaluation
1891 intS=funcmass->GetParameter(np);
1892 intSerr=funcmass->GetParError(np);
1894 cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1895 Double_t background,errbackground;
1896 Background(min,max,background,errbackground);
1898 //signal +/- error in the range
1900 Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1901 signal=mass - background;
1902 errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1906 //_________________________________________________________________________
1908 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1909 // Return background integral in mean+- n sigma
1912 cout<<"Use MassFitter method before Background"<<endl;
1915 Double_t min=fMass-nOfSigma*fSigmaSgn;
1916 Double_t max=fMass+nOfSigma*fSigmaSgn;
1918 Background(min,max,background,errbackground);
1923 //___________________________________________________________________________
1925 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1926 // Return background integral in a range
1929 cout<<"Use MassFitter method before Background"<<endl;
1934 TString bkgname="funcbkgRecalc";
1935 TString bkg1name="funcbkg1Recalc";
1938 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1939 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1941 cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1946 Double_t intB,intBerr;
1948 //relative error evaluation: from final parameters of the fit
1949 if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1951 intB=funcbkg->GetParameter(0);
1952 intBerr=funcbkg->GetParError(0);
1953 cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1956 //relative error evaluation: from histo
1958 intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1960 for(Int_t i=1;i<=fSideBandl;i++){
1961 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1963 for(Int_t i=fSideBandr;i<=fNbin;i++){
1964 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1967 intBerr=TMath::Sqrt(sum2);
1968 cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1970 cout<<"Last estimation of bkg error is used"<<endl;
1972 //backround +/- error in the range
1973 if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1978 background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1979 errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1980 //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1987 //__________________________________________________________________________
1989 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
1990 // Return significance in mean+- n sigma
1992 Double_t min=fMass-nOfSigma*fSigmaSgn;
1993 Double_t max=fMass+nOfSigma*fSigmaSgn;
1994 Significance(min, max, significance, errsignificance);
1999 //__________________________________________________________________________
2001 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
2002 // Return significance integral in a range
2004 Double_t signal,errsignal,background,errbackground;
2005 Signal(min, max,signal,errsignal);
2006 Background(min, max,background,errbackground);
2008 if (signal+background <= 0.){
2009 cout<<"Cannot calculate significance because of div by 0!"<<endl;
2015 significance = signal/TMath::Sqrt(signal+background);
2017 errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);