PWG3/vertexingHF/AliHFMassFitter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   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  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliHFMassFitter for the fit of invariant mass distribution
  19 // of charmed mesons
  20 //
  21 // Author: C.Bianchin, chiara.bianchin@pd.infn.it
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TCanvas.h>
  25 #include <TMinuit.h>
  26 #include <TGraph.h>
  27 #include <TVirtualFitter.h>
  28 #include <TClonesArray.h>
  29
  30 #include "AliHFMassFitter.h"
  31
  32
  33 ClassImp(AliHFMassFitter)
  34
  35  //************** constructors
  36 AliHFMassFitter::AliHFMassFitter() :
  37   TNamed(),
  38   fhistoInvMass(0),
  39   fminMass(0),
  40   fmaxMass(0),
  41   fNbin(0),
  42   fParsSize(1),
  43   fFitPars(0),
  44   fWithBkg(0),
  45   ftypeOfFit4Bkg(0),
  46   ftypeOfFit4Sgn(0),
  47   ffactor(1),
  48   fntuParam(0),
  49   fMass(1.85),
  50   fSigmaSgn(0.012),
  51   fSideBands(0),
  52   fSideBandl(0),
  53   fSideBandr(0),
  54   fcounter(0),
  55   fContourGraph(0)
  56 {
  57   // default constructor
  58
  59   cout<<"Default constructor"<<endl;
  60
  61 }
  62
  63 //___________________________________________________________________________
  64
  65 AliHFMassFitter::AliHFMassFitter (TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
  66  TNamed(),
  67  fhistoInvMass(0),
  68  fminMass(0),
  69  fmaxMass(0),
  70  fNbin(0),
  71  fParsSize(1),
  72  fFitPars(0),
  73  fWithBkg(0),
  74  ftypeOfFit4Bkg(0),
  75  ftypeOfFit4Sgn(0),
  76  ffactor(1),
  77  fntuParam(0),
  78  fMass(1.85),
  79  fSigmaSgn(0.012),
  80  fSideBands(0),
  81  fSideBandl(0),
  82  fSideBandr(0),
  83  fcounter(0),
  84  fContourGraph(0)
  85 {
  86   // standard constructor
  87
  88   fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
  89   fminMass=minvalue;
  90   fmaxMass=maxvalue;
  91   if(rebin!=1) RebinMass(rebin);
  92   else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
  93   ftypeOfFit4Bkg=fittypeb;
  94   ftypeOfFit4Sgn=fittypes;
  95   if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
  96   else fWithBkg=kTRUE;
  97   if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
  98   else  cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
  99
 100   ComputeParSize();
 101   fFitPars=new Float_t[fParsSize];
 102
 103   fContourGraph=new TList();
 104   fContourGraph->SetOwner();
 105
 106 }
 107
 108
 109 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
 110   TNamed(),
 111   fhistoInvMass(mfit.fhistoInvMass),
 112   fminMass(mfit.fminMass),
 113   fmaxMass(mfit.fmaxMass),
 114   fNbin(mfit.fNbin),
 115   fParsSize(mfit.fParsSize),
 116   fFitPars(0),
 117   fWithBkg(mfit.fWithBkg),
 118   ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
 119   ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
 120   ffactor(mfit.ffactor),
 121   fntuParam(mfit.fntuParam),
 122   fMass(mfit.fMass),
 123   fSigmaSgn(mfit.fSigmaSgn),
 124   fSideBands(mfit.fSideBands),
 125   fSideBandl(mfit.fSideBandl),
 126   fSideBandr(mfit.fSideBandr),
 127   fcounter(mfit.fcounter),
 128   fContourGraph(mfit.fContourGraph)
 129 {
 130   //copy constructor
 131
 132   if(mfit.fParsSize > 0){
 133     fFitPars=new Float_t[fParsSize];
 134     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 135   }
 136   //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
 137 }
 138
 139 //_________________________________________________________________________
 140
 141 AliHFMassFitter::~AliHFMassFitter() {
 142   cout<<"AliHFMassFitter destructor called"<<endl;
 143   if(fhistoInvMass) {
 144     cout<<"deleting histogram..."<<endl;
 145     delete fhistoInvMass;
 146     fhistoInvMass=NULL;
 147   }
 148   if(fntuParam){
 149     cout<<"deleting ntuple..."<<endl;
 150     delete fntuParam;
 151
 152     fntuParam=NULL;
 153   }
 154
 155   if(fFitPars) {
 156     delete[] fFitPars;
 157     cout<<"deleting parameter array..."<<endl;
 158     fFitPars=NULL;
 159   }
 160
 161   fcounter = 0;
 162 }
 163
 164 //_________________________________________________________________________
 165
 166 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
 167
 168   //assignment operator
 169
 170   if(&mfit == this) return *this;
 171   fhistoInvMass= mfit.fhistoInvMass;
 172   fminMass= mfit.fminMass;
 173   fmaxMass= mfit.fmaxMass;
 174   fNbin= mfit.fNbin;
 175   fParsSize= mfit.fParsSize;
 176   fWithBkg= mfit.fWithBkg;
 177   ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
 178   ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
 179   ffactor= mfit.ffactor;
 180   fntuParam= mfit.fntuParam;
 181   fMass= mfit.fMass;
 182   fSigmaSgn= mfit.fSigmaSgn;
 183   fSideBands= mfit.fSideBands;
 184   fSideBandl= mfit.fSideBandl;
 185   fSideBandr= mfit.fSideBandr;
 186   fcounter= mfit.fcounter;
 187   fContourGraph= mfit.fContourGraph;
 188
 189   if(mfit.fParsSize > 0){
 190     if(fFitPars) {
 191       delete[] fFitPars;
 192       fFitPars=NULL;
 193     }
 194     fFitPars=new Float_t[fParsSize];
 195     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 196   }
 197 // fFitPars=new Float_t[fParsSize];
 198 //   for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
 199
 200   return *this;
 201 }
 202
 203 //************ tools & settings
 204
 205 //__________________________________________________________________________
 206
 207 void AliHFMassFitter::ComputeParSize() {
 208   switch (ftypeOfFit4Bkg) {//npar backround func
 209   case 0:
 210     fParsSize = 2*3;
 211     break;
 212   case 1:
 213     fParsSize = 2*3;
 214     break;
 215   case 2:
 216     fParsSize = 3*3;
 217     break;
 218   case 3:
 219     fParsSize = 1*3;
 220     break;
 221   default:
 222     cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
 223     break;
 224   }
 225
 226   fParsSize += 3; // npar refl
 227   fParsSize += 3; // npar signal gaus
 228
 229   fParsSize*=2;   // add errors
 230   cout<<"Parameters array size "<<fParsSize<<endl;
 231 }
 232
 233 //___________________________________________________________________________
 234 void AliHFMassFitter::SetHisto(TH1F *histoToFit){
 235   //fhistoInvMass = (TH1F*)histoToFit->Clone();
 236   fhistoInvMass = new TH1F(*histoToFit);
 237   cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
 238 }
 239
 240 //___________________________________________________________________________
 241
 242   void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
 243     ftypeOfFit4Bkg = fittypeb;
 244     ftypeOfFit4Sgn = fittypes;
 245     ComputeParSize();
 246     fFitPars = new Float_t[fParsSize];
 247
 248 }
 249
 250 //___________________________________________________________________________
 251
 252 void AliHFMassFitter::Reset() {
 253   cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
 254   fMass=1.85;
 255   fSigmaSgn=0.012;
 256   cout<<"Reset "<<fhistoInvMass<<endl;
 257   if(fhistoInvMass) {
 258     //cout<<"esiste"<<endl;
 259     //delete fhistoInvMass;
 260     fhistoInvMass=NULL;
 261     cout<<fhistoInvMass<<endl;
 262   }
 263   else cout<<"histogram doesn't exist, do not delete"<<endl;
 264
 265 }
 266
 267 //_________________________________________________________________________
 268
 269 void AliHFMassFitter::InitNtuParam(char *ntuname) {
 270   // Create ntuple to keep fit parameters
 271   fntuParam=0;
 272   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");
 273
 274 }
 275
 276 //_________________________________________________________________________
 277
 278 void AliHFMassFitter::FillNtuParam() {
 279   // Fill ntuple with fit parameters
 280
 281   Float_t nothing=0.;
 282
 283   if (ftypeOfFit4Bkg==2) {
 284       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 285       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 286       fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
 287       fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
 288       fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
 289       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
 290       fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
 291       fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
 292       fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
 293       fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
 294       fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
 295       fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
 296       fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
 297       fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
 298       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
 299
 300       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
 301       fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
 302       fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
 303       fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
 304       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
 305       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
 306       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
 307       fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
 308       fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
 309       fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
 310       fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
 311       fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
 312       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
 313       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
 314       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
 315
 316   } else {
 317
 318     if(ftypeOfFit4Bkg==3){
 319       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 320       fntuParam->SetBranchAddress("slope1",&nothing);
 321       fntuParam->SetBranchAddress("conc1",&nothing);
 322       fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
 323       fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
 324       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
 325       fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
 326       fntuParam->SetBranchAddress("slope2",&nothing);
 327       fntuParam->SetBranchAddress("conc2",&nothing);
 328       fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
 329       fntuParam->SetBranchAddress("slope3",&nothing);
 330       fntuParam->SetBranchAddress("conc3",&nothing);
 331       fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
 332       fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
 333       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
 334
 335       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
 336       fntuParam->SetBranchAddress("slope1Err",&nothing);
 337       fntuParam->SetBranchAddress("conc1Err",&nothing);
 338       fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
 339       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
 340       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
 341       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
 342       fntuParam->SetBranchAddress("slope2Err",&nothing);
 343       fntuParam->SetBranchAddress("conc2Err",&nothing);
 344       fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
 345       fntuParam->SetBranchAddress("slope3Err",&nothing);
 346       fntuParam->SetBranchAddress("conc3Err",&nothing);
 347       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
 348       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
 349       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
 350
 351     }
 352     else{
 353       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 354       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 355       fntuParam->SetBranchAddress("conc1",&nothing);
 356       fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
 357       fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
 358       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
 359       fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
 360       fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
 361       fntuParam->SetBranchAddress("conc2",&nothing);
 362       fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
 363       fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
 364       fntuParam->SetBranchAddress("conc3",&nothing);
 365       fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
 366       fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
 367       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
 368
 369       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
 370       fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
 371       fntuParam->SetBranchAddress("conc1Err",&nothing);
 372       fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
 373       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
 374       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
 375       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
 376       fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
 377       fntuParam->SetBranchAddress("conc2Err",&nothing);
 378       fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
 379       fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
 380       fntuParam->SetBranchAddress("conc3Err",&nothing);
 381       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
 382       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
 383       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
 384     }
 385
 386   }
 387   fntuParam->TTree::Fill();
 388 }
 389
 390 //_________________________________________________________________________
 391
 392 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
 393   // Create, fill and return ntuple with fit parameters
 394
 395   InitNtuParam(ntuname);
 396   FillNtuParam();
 397   return fntuParam;
 398 }
 399 //_________________________________________________________________________
 400
 401 void AliHFMassFitter::RebinMass(Int_t bingroup){
 402   // Rebin invariant mass histogram
 403
 404   if(bingroup<1){
 405     cout<<"Error! Cannot group "<<bingroup<<" bins\n";
 406     fNbin=fhistoInvMass->GetNbinsX();
 407     cout<<"Kept original number of bins: "<<fNbin<<endl;
 408   } else{
 409     fhistoInvMass->Rebin(bingroup);
 410     fNbin = fhistoInvMass->GetNbinsX();
 411     cout<<"New number of bins: "<<fNbin<<endl;
 412   }
 413
 414
 415 }
 416
 417 //************* fit
 418
 419 //___________________________________________________________________________
 420
 421 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
 422   // Fit function for signal+background
 423
 424
 425   //exponential or linear fit
 426   //
 427   // par[0] = tot integral
 428   // par[1] = slope
 429   // par[2] = gaussian integral
 430   // par[3] = gaussian mean
 431   // par[4] = gaussian sigma
 432
 433   Double_t total,bkg=0,sgn=0;
 434
 435   if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
 436     if(ftypeOfFit4Sgn == 0) {
 437
 438       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 439       bkg = FitFunction4Bkg(x,parbkg);
 440     }
 441     if(ftypeOfFit4Sgn == 1) {
 442       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
 443       bkg = FitFunction4Bkg(x,parbkg);
 444     }
 445
 446     sgn = FitFunction4Sgn(x,&par[2]);
 447
 448   }
 449
 450   //polynomial fit
 451
 452     // par[0] = tot integral
 453     // par[1] = coef1
 454     // par[2] = coef2
 455     // par[3] = gaussian integral
 456     // par[4] = gaussian mean
 457     // par[5] = gaussian sigma
 458
 459   if (ftypeOfFit4Bkg==2) {
 460
 461     if(ftypeOfFit4Sgn == 0) {
 462
 463       Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
 464       bkg = FitFunction4Bkg(x,parbkg);
 465     }
 466     if(ftypeOfFit4Sgn == 1) {
 467
 468       Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
 469      bkg = FitFunction4Bkg(x,parbkg);
 470     }
 471
 472     sgn = FitFunction4Sgn(x,&par[3]);
 473   }
 474
 475   if (ftypeOfFit4Bkg==3) {
 476
 477     if(ftypeOfFit4Sgn == 0) {
 478         bkg=FitFunction4Bkg(x,par);
 479         sgn=FitFunction4Sgn(x,&par[1]);
 480     }
 481     if(ftypeOfFit4Sgn == 1) {
 482       Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
 483       bkg=FitFunction4Bkg(x,parbkg);
 484       sgn=FitFunction4Sgn(x,&par[1]);
 485     }
 486   }
 487
 488   total = bkg + sgn;
 489
 490   return  total;
 491 }
 492
 493 //_________________________________________________________________________
 494 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
 495   // Fit function for the signal
 496
 497   //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 498   //Par:
 499   // * [0] = integralSgn
 500   // * [1] = mean
 501   // * [2] = sigma
 502   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 503
 504   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]);
 505
 506 }
 507
 508 //__________________________________________________________________________
 509
 510 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
 511   // Fit function for the background
 512
 513   Double_t maxDeltaM = 4.*fSigmaSgn;
 514   if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
 515     TF1::RejectPoint();
 516     return 0;
 517   }
 518   Int_t firstPar=0;
 519   Double_t gaus2=0,total=-1;
 520   if(ftypeOfFit4Sgn == 1){
 521     firstPar=3;
 522     //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 523     //Par:
 524     // * [0] = integralSgn
 525     // * [1] = mean
 526     // * [2] = sigma
 527     //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 528     gaus2 = FitFunction4Sgn(x,par);
 529   }
 530
 531   switch (ftypeOfFit4Bkg){
 532   case 0:
 533     //exponential
 534     //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
 535     //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
 536     //as integralTot- integralGaus (=par [2])
 537     //Par:
 538     // * [0] = integralBkg;
 539     // * [1] = B;
 540     //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
 541     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]);
 542     break;
 543   case 1:
 544     //linear
 545     //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)
 546     // * [0] = integralBkg;
 547     // * [1] = b;
 548     total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
 549     break;
 550   case 2:
 551     //polynomial
 552     //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
 553     //+ 1/3*c*(max^3-min^3) ->
 554     //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
 555     // * [0] = integralBkg;
 556     // * [1] = b;
 557     // * [2] = c;
 558     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));
 559     break;
 560   case 3:
 561     total=par[0+firstPar];
 562     break;
 563 //   default:
 564 //     Types of Fit Functions for Background:
 565 //     * 0 = exponential;
 566 //     * 1 = linear;
 567 //     * 2 = polynomial 2nd order
 568 //     * 3 = no background"<<endl;
 569
 570   }
 571   return total+gaus2;
 572 }
 573
 574 //__________________________________________________________________________
 575 Bool_t AliHFMassFitter::SideBandsBounds(){
 576
 577   Double_t width=fhistoInvMass->GetBinWidth(8);
 578   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 579   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 580   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 581
 582   if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
 583     cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
 584     return kFALSE;
 585   }
 586
 587   if((fminMass!=minHisto || fmaxMass!= maxHisto) && (fMass-4.*fSigmaSgn-fminMass) <= 0){
 588     Double_t coeff = (fMass-fminMass)/fSigmaSgn;
 589
 590     fSideBandl=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
 591     fSideBandr=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
 592     cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
 593     if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
 594
 595     if (coeff<2) {
 596       cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
 597       ftypeOfFit4Bkg=3;
 598       //set binleft and right without considering SetRangeFit- anyway no bkg!
 599       fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 600       fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 601     }
 602   }
 603   else {
 604   fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 605   fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 606 //   cout<<"\tfMass = "<<fMass<<"\tfSigmaSgn = "<<fSigmaSgn<<"\tminHisto = "<<minHisto<<endl;
 607 //   cout<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 608   }
 609
 610   if (fSideBandl==0) {
 611     cout<<"Error! Range too little";
 612     return kFALSE;
 613   }
 614   return kTRUE;
 615 }
 616
 617 //__________________________________________________________________________
 618
 619 void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
 620   if (fSideBandl==0 && fSideBandr==0){
 621     cout<<"Use MassFitter method first"<<endl;
 622     return;
 623   }
 624   left=fSideBandl;
 625   right=fSideBandr;
 626 }
 627
 628 //__________________________________________________________________________
 629
 630 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 631   // Main method of the class: performs the fit of the histogram
 632
 633   //Set default fitter Minuit in order to use gMinuit in the contour plots
 634   TVirtualFitter::SetDefaultFitter("Minuit");
 635
 636   Int_t nFitPars=0; //total function's number of parameters
 637   switch (ftypeOfFit4Bkg){
 638   case 0:
 639     nFitPars=5; //3+2
 640     break;
 641   case 1:
 642     nFitPars=5; //3+2
 643     break;
 644   case 2:
 645     nFitPars=6; //3+3
 646     break;
 647   case 3:
 648     nFitPars=4; //3+1
 649     break;
 650   }
 651
 652   Int_t bkgPar = nFitPars-3; //background function's number of parameters
 653
 654   cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 655
 656
 657   TString listname="contourplot";
 658   listname+=fcounter;
 659   if(!fContourGraph){
 660     fContourGraph=new TList();
 661     fContourGraph->SetOwner();
 662   }
 663
 664   fContourGraph->SetName(listname);
 665
 666
 667   //function names
 668   TString bkgname="funcbkg";
 669   TString bkg1name="funcbkg1";
 670   TString massname="funcmass";
 671
 672   //Total integral
 673   Double_t totInt = fhistoInvMass->Integral("width");
 674   cout<<"Integral"<<endl;
 675
 676   fSideBands = kTRUE;
 677   Double_t width=fhistoInvMass->GetBinWidth(8);
 678   cout<<"fNbin"<<fNbin<<endl;
 679   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 680   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 681   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 682   Bool_t ok=SideBandsBounds();
 683   if(!ok) return kFALSE;
 684
 685   //sidebands integral - first approx (from histo)
 686   Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
 687   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 688   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 689   if (sideBandsInt<=0) {
 690     cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
 691     return kFALSE;
 692   }
 693
 694   /*Fit Bkg*/
 695
 696
 697   TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 698   cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
 699
 700   funcbkg->SetLineColor(2); //red
 701
 702   //first fit for bkg: approx bkgint
 703
 704   switch (ftypeOfFit4Bkg) {
 705   case 0: //gaus+expo
 706     funcbkg->SetParNames("BkgInt","Slope");
 707     funcbkg->SetParameters(sideBandsInt,-2.);
 708     break;
 709   case 1:
 710     funcbkg->SetParNames("BkgInt","Slope");
 711     funcbkg->SetParameters(sideBandsInt,-100.);
 712     break;
 713   case 2:
 714     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
 715     funcbkg->SetParameters(sideBandsInt,-10.,5);
 716     break;
 717   case 3:
 718     if(ftypeOfFit4Sgn==0){
 719       funcbkg->SetParNames("Const");
 720       funcbkg->SetParameter(0,0.);
 721       funcbkg->FixParameter(0,0.);
 722     }
 723     break;
 724   default:
 725     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
 726     return kFALSE;
 727     break;
 728   }
 729   cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
 730   Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
 731
 732   Double_t intbkg1=0,slope1=0,conc1=0;
 733   //if only signal and reflection: skip
 734   if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
 735     ftypeOfFit4Sgn=0;
 736     fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
 737
 738     for(Int_t i=0;i<bkgPar;i++){
 739       fFitPars[i]=funcbkg->GetParameter(i);
 740       //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 741       fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
 742       //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 743     }
 744     fSideBands = kFALSE;
 745     //intbkg1 = funcbkg->GetParameter(0);
 746     funcbkg->SetRange(fminMass,fmaxMass);
 747     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
 748     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
 749     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
 750     cout<<"Primo fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
 751   }
 752   else cout<<"\t\t//"<<endl;
 753
 754   ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
 755   TF1 *funcbkg1=0;
 756   if (ftypeOfFit4Sgn == 1) {
 757     cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
 758     bkgPar+=3;
 759
 760     cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 761
 762     funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 763     cout<<"Function name = "<<funcbkg1->GetName()<<endl;
 764
 765     funcbkg1->SetLineColor(2); //red
 766
 767     if(ftypeOfFit4Bkg==2){
 768       cout<<"*** Polynomial Fit ***"<<endl;
 769       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
 770       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
 771 //     cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
 772 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
 773     } else{
 774       if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
 775         {
 776           cout<<"*** No background Fit ***"<<endl;
 777           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
 778           funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
 779           funcbkg1->FixParameter(3,0.);
 780         } else{ //expo or linear
 781           if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
 782           if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
 783           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
 784           funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
 785         }
 786     }
 787     fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
 788
 789     for(Int_t i=0;i<bkgPar;i++){
 790       fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
 791       //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
 792       fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
 793       //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
 794     }
 795
 796     intbkg1=funcbkg1->GetParameter(3);
 797     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
 798     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
 799
 800   } else {
 801     bkgPar+=3;
 802
 803     for(Int_t i=0;i<3;i++){
 804       fFitPars[bkgPar-3+i]=0.;
 805       cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
 806       fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
 807       cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
 808     }
 809
 810     for(Int_t i=0;i<bkgPar-3;i++){
 811       fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
 812       cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 813       fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
 814       cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 815     }
 816
 817
 818   }
 819
 820   //sidebands integral - second approx (from fit)
 821   fSideBands = kFALSE;
 822   Double_t bkgInt;
 823   cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
 824   if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
 825   else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
 826   cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
 827
 828   //Signal integral - first approx
 829   Double_t sgnInt;
 830   sgnInt = totInt-bkgInt;
 831   cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
 832
 833   /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
 834   TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
 835   cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
 836   funcmass->SetLineColor(4); //blue
 837
 838   //Set parameters
 839   cout<<"\nTOTAL FIT"<<endl;
 840
 841   if(nFitPars==5){
 842     funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
 843     funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
 844     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 845 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 846     funcmass->FixParameter(0,totInt);
 847   }
 848   if (nFitPars==6){
 849     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
 850     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
 851 //     cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 852 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 853     funcmass->FixParameter(0,totInt);
 854   }
 855   if(nFitPars==4){
 856     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
 857     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
 858     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
 859     funcmass->FixParameter(0,0.);
 860     //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 861     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 862
 863   }
 864
 865   fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
 866   cout<<"fit done"<<endl;
 867   //reset value of fMass and fSigmaSgn to those found from fit
 868   fMass=funcmass->GetParameter(nFitPars-2);
 869   fSigmaSgn=funcmass->GetParameter(nFitPars-1);
 870
 871   for(Int_t i=0;i<nFitPars;i++){
 872     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
 873     fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
 874     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
 875   }
 876   /*
 877   //check: cout parameters
 878   for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
 879     cout<<i<<"\t"<<fFitPars[i]<<endl;
 880     }
 881   */
 882
 883 //   if(draw){
 884 //     TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
 885 //     TH1F *fhistocopy=new TH1F(*fhistoInvMass);
 886 //     canvas->cd();
 887 //     fhistocopy->DrawClone();
 888 //     if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
 889 //     else funcbkg->DrawClone("sames");
 890 //     funcmass->DrawClone("sames");
 891 //     cout<<"Drawn"<<endl;
 892 //   }
 893
 894   if(funcmass->GetParameter(nFitPars-1) <0 || funcmass->GetParameter(nFitPars-2) <0 || funcmass->GetParameter(nFitPars-3) <0 ) {
 895     cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
 896     return kFALSE;
 897   }
 898
 899   //increase counter of number of fits done
 900   fcounter++;
 901
 902   //contour plots
 903   if(draw){
 904
 905     for (Int_t kpar=1; kpar<nFitPars;kpar++){
 906
 907       for(Int_t jpar=kpar+1;jpar<nFitPars;jpar++){
 908         cout<<"Par "<<kpar<<" and "<<jpar<<endl;
 909
 910         // produce 2 contours per couple of parameters
 911         TGraph* cont[2] = {0x0, 0x0};
 912         const Double_t errDef[2] = {1., 4.};
 913         for (Int_t i=0; i<2; i++) {
 914           gMinuit->SetErrorDef(errDef[i]);
 915           cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
 916         }
 917
 918         if(!cont[0] || !cont[1]){
 919           cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
 920           continue;
 921         }
 922
 923         // set graph titles and add them to the list
 924         TString title = "Contour plot";
 925         TString titleX = funcmass->GetParName(kpar);
 926         TString titleY = funcmass->GetParName(jpar);
 927         for (Int_t i=0; i<2; i++) {
 928           cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
 929           cont[i]->SetTitle(title);
 930           cont[i]->GetXaxis()->SetTitle(titleX);
 931           cont[i]->GetYaxis()->SetTitle(titleY);
 932           cont[i]->GetYaxis()->SetLabelSize(0.033);
 933           cont[i]->GetYaxis()->SetTitleSize(0.033);
 934           cont[i]->GetYaxis()->SetTitleOffset(1.67);
 935
 936           fContourGraph->Add(cont[i]);
 937         }
 938
 939         // plot them
 940         TString cvname = Form("c%d%d", kpar, jpar);
 941         TCanvas *c4=new TCanvas(cvname,cvname,600,600);
 942         c4->cd();
 943         cont[1]->SetFillColor(38);
 944         cont[1]->Draw("alf");
 945         cont[0]->SetFillColor(9);
 946         cont[0]->Draw("lf");
 947
 948       }
 949
 950     }
 951
 952   }
 953
 954   if (ftypeOfFit4Sgn == 1 && funcbkg1) {
 955     delete funcbkg1;
 956     funcbkg1=NULL;
 957   }
 958   if (funcbkg) {
 959     delete funcbkg;
 960     funcbkg=NULL;
 961   }
 962   if (funcmass) {
 963     delete funcmass;
 964     funcmass=NULL;
 965   }
 966
 967   AddFunctionsToHisto();
 968   if (draw) DrawFit();
 969
 970
 971   return kTRUE;
 972 }
 973 //_________________________________________________________________________
 974 Double_t AliHFMassFitter::GetChiSquare() const{
 975   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
 976   return funcmass->GetChisquare();
 977 }
 978
 979 //_________________________________________________________________________
 980 Double_t AliHFMassFitter::GetReducedChiSquare() const{
 981   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
 982   return funcmass->GetChisquare()/funcmass->GetNDF();
 983 }
 984
 985 //*********output
 986
 987 //_________________________________________________________________________
 988 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
 989   // Return fit parameters
 990
 991   for(Int_t i=0;i<fParsSize;i++){
 992     vector[i]=fFitPars[i];
 993   }
 994 }
 995
 996
 997 //_________________________________________________________________________
 998 void AliHFMassFitter::AddFunctionsToHisto(){
 999
1000   cout<<"AddFunctionsToHisto called"<<endl;
1001   TString bkgname = "funcbkg";
1002   Int_t np=-99;
1003   switch (ftypeOfFit4Bkg){
1004   case 0: //expo
1005     np=2;
1006     break;
1007   case 1: //linear
1008     np=2;
1009     break;
1010   case 2: //pol2
1011     np=3;
1012     break;
1013   case 3: //no bkg
1014     np=1;
1015     break;
1016   }
1017   if (ftypeOfFit4Sgn == 1) {
1018     bkgname += 1;
1019     np+=3;
1020   }
1021
1022   TString bkgnamesave=bkgname;
1023   TString testname=bkgname;
1024   testname += "FullRange";
1025   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1026   if(testfunc){
1027     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1028     return;
1029   }
1030
1031   TList *hlist=fhistoInvMass->GetListOfFunctions();
1032   hlist->ls();
1033
1034   TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1035   if(!b){
1036     cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1037     return;
1038   }
1039
1040   bkgname += "FullRange";
1041   TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
1042   //cout<<bfullrange->GetName()<<endl;
1043   for(Int_t i=0;i<np;i++){
1044     //cout<<i<<" di "<<np<<endl;
1045     bfullrange->SetParName(i,b->GetParName(i));
1046     bfullrange->SetParameter(i,b->GetParameter(i));
1047     bfullrange->SetParError(i,b->GetParError(i));
1048   }
1049   bfullrange->SetLineStyle(4);
1050   bfullrange->SetLineColor(14);
1051
1052   bkgnamesave += "Recalc";
1053
1054   TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
1055
1056   TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1057
1058   if (!mass){
1059     cout<<"funcmass doesn't exist "<<endl;
1060     return;
1061   }
1062
1063   blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(np));
1064   blastpar->SetParError(0,mass->GetParError(np));
1065   if (np>=2) {
1066     blastpar->SetParameter(1,mass->GetParameter(1));
1067     blastpar->SetParError(1,mass->GetParError(1));
1068   }
1069   if (np==3) {
1070     blastpar->SetParameter(2,mass->GetParameter(2));
1071     blastpar->SetParError(2,mass->GetParError(2));
1072   }
1073
1074   blastpar->SetLineStyle(1);
1075   blastpar->SetLineColor(2);
1076
1077   hlist->Add(bfullrange->Clone());
1078   hlist->Add(blastpar->Clone());
1079   hlist->ls();
1080
1081   if(bfullrange) {
1082     delete bfullrange;
1083     bfullrange=NULL;
1084   }
1085   if(blastpar) {
1086     delete blastpar;
1087     blastpar=NULL;
1088   }
1089 }
1090
1091 //_________________________________________________________________________
1092
1093 TH1F* AliHFMassFitter::GetHistoClone() const{
1094
1095   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1096   return hout;
1097 }
1098 //_________________________________________________________________________
1099
1100 void AliHFMassFitter::WriteHisto(TString path) {
1101   if (fcounter == 0) {
1102     cout<<"Use MassFitter method before WriteHisto"<<endl;
1103     return;
1104   }
1105   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1106
1107   path += "HFMassFitterOutput.root";
1108   TFile *output;
1109
1110   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1111   else output = new TFile(path.Data(),"update");
1112   output->cd();
1113   hget->Write();
1114   fContourGraph->Write();
1115   output->Close();
1116
1117   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1118
1119   if(output) {
1120     delete output;
1121     output=NULL;
1122   }
1123
1124 }
1125
1126 //_________________________________________________________________________
1127
1128 void AliHFMassFitter::WriteNtuple(TString path) const{
1129   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1130   path += "HFMassFitterOutput.root";
1131   TFile *output = new TFile(path.Data(),"update");
1132   output->cd();
1133   fntuParam->Write();
1134   //nget->Write();
1135   output->Close();
1136   //cout<<nget->GetName()<<" written in "<<path<<endl;
1137   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1138   /*
1139   if(nget) {
1140     //delete nget;
1141     nget=NULL;
1142   }
1143   */
1144   if(output) {
1145     delete output;
1146     output=NULL;
1147   }
1148 }
1149
1150 //_________________________________________________________________________
1151
1152 void AliHFMassFitter::DrawFit() const{
1153
1154   TString cvtitle="fit of ";
1155   cvtitle+=fhistoInvMass->GetName();
1156   TString cvname="c";
1157   cvname+=fcounter;
1158
1159   TCanvas *c=new TCanvas(cvname,cvtitle);
1160   c->cd();
1161   GetHistoClone()->Draw("hist");
1162   GetHistoClone()->GetFunction("funcbkgFullRange")->Draw("same");
1163   GetHistoClone()->GetFunction("funcbkgRecalc")->Draw("same");
1164   GetHistoClone()->GetFunction("funcmass")->Draw("same");
1165
1166
1167 }
1168
1169
1170 //_________________________________________________________________________
1171
1172 void AliHFMassFitter::PrintParTitles() const{
1173   TF1 *f=fhistoInvMass->GetFunction("funcmass");
1174   if(!f) {
1175     cout<<"Fit function not found"<<endl;
1176     return;
1177   }
1178
1179   Int_t np=0;
1180   switch (ftypeOfFit4Bkg){
1181   case 0: //expo
1182     np=2;
1183     break;
1184   case 1: //linear
1185     np=2;
1186     break;
1187   case 2: //pol2
1188     np=3;
1189     break;
1190   case 3: //no bkg
1191     np=1;
1192     break;
1193   }
1194
1195   np+=3; //3 parameter for signal
1196   if (ftypeOfFit4Sgn == 1) np+=3;
1197
1198   cout<<"Parameter Titles \n";
1199   for(Int_t i=0;i<np;i++){
1200     cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1201   }
1202   cout<<endl;
1203
1204 }
1205
1206 //************ significance
1207
1208 //_________________________________________________________________________
1209
1210 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1211   // Return signal integral in mean+- n sigma
1212
1213   if(fcounter==0) {
1214     cout<<"Use MassFitter method before Signal"<<endl;
1215     return;
1216   }
1217
1218   //functions names
1219   TString bkgname= "funcbkgRecalc";
1220   TString bkg1name="funcbkg1Recalc";
1221   TString massname="funcmass";
1222
1223
1224   TF1 *funcbkg=0;
1225   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1226   if(!funcmass){
1227     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1228     return;
1229   }
1230
1231   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1232   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1233
1234   if(!funcbkg){
1235     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1236     return;
1237   }
1238
1239   Int_t np=-99;
1240   switch (ftypeOfFit4Bkg){
1241   case 0: //expo
1242     np=2;
1243     break;
1244   case 1: //linear
1245     np=2;
1246     break;
1247   case 2: //pol2
1248     np=3;
1249     break;
1250   case 3: //no bkg
1251     np=1;
1252     break;
1253   }
1254
1255   Double_t intS,intSerr;
1256
1257  //relative error evaluation
1258   intS=funcmass->GetParameter(np);
1259   intSerr=funcmass->GetParError(np);
1260   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1261   Double_t background,errbackground;
1262   Background(nOfSigma,background,errbackground);
1263
1264   //signal +/- error in nsigma
1265   Double_t min=fMass-nOfSigma*fSigmaSgn;
1266   Double_t max=fMass+nOfSigma*fSigmaSgn;
1267
1268   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1269   signal=mass - background;
1270   errsignal=TMath::Sqrt((intSerr/intS*mass)*(intSerr/intS*mass)/*assume relative error is the same as for total integral*/ + errbackground*errbackground);
1271   return;
1272 }
1273
1274 //_________________________________________________________________________
1275
1276 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1277
1278   // Return signal integral in a range
1279
1280   if(fcounter==0) {
1281     cout<<"Use MassFitter method before Signal"<<endl;
1282     return;
1283   }
1284
1285   //functions names
1286   TString bkgname="funcbkgRecalc";
1287   TString bkg1name="funcbkg1Recalc";
1288   TString massname="funcmass";
1289
1290
1291   TF1 *funcbkg=0;
1292   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1293   if(!funcmass){
1294     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1295     return;
1296   }
1297
1298   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1299   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1300
1301   if(!funcbkg){
1302     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1303     return;
1304   }
1305
1306   Int_t np=-99;
1307   switch (ftypeOfFit4Bkg){
1308   case 0: //expo
1309     np=2;
1310     break;
1311   case 1: //linear
1312     np=2;
1313     break;
1314   case 2: //pol2
1315     np=3;
1316     break;
1317   case 3: //no bkg
1318     np=1;
1319     break;
1320   }
1321
1322   Double_t intS,intSerr;
1323
1324  //relative error evaluation
1325   intS=funcmass->GetParameter(np);
1326   intSerr=funcmass->GetParError(np);
1327
1328   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1329   Double_t background,errbackground;
1330   Background(min,max,background,errbackground);
1331
1332   //signal +/- error in the range
1333
1334   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1335   signal=mass - background;
1336   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1337
1338 }
1339
1340 //_________________________________________________________________________
1341
1342 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1343   // Return background integral in mean+- n sigma
1344
1345   if(fcounter==0) {
1346     cout<<"Use MassFitter method before Background"<<endl;
1347     return;
1348   }
1349
1350   //functions names
1351   TString bkgname="funcbkgRecalc";
1352   TString bkg1name="funcbkg1Recalc";
1353
1354   TF1 *funcbkg=0;
1355   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1356   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1357   if(!funcbkg){
1358     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1359     return;
1360   }
1361
1362   Double_t intB,intBerr;
1363
1364   //relative error evaluation: from final parameters of the fit
1365   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1366   else{
1367     intB=funcbkg->GetParameter(0);
1368     intBerr=funcbkg->GetParError(0);
1369     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1370   }
1371
1372   Double_t min=fMass-nOfSigma*fSigmaSgn;
1373   Double_t max=fMass+nOfSigma*fSigmaSgn;
1374
1375   //relative error evaluation: from histo
1376
1377   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1378   Double_t sum2=0;
1379   for(Int_t i=1;i<=fSideBandl;i++){
1380     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1381   }
1382   for(Int_t i=fSideBandr;i<=fNbin;i++){
1383     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1384   }
1385
1386   intBerr=TMath::Sqrt(sum2);
1387   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1388
1389   cout<<"Last estimation of bkg error is used"<<endl;
1390
1391   //backround +/- error in nsigma
1392   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1393     background = 0;
1394     errbackground = 0;
1395   }
1396   else{
1397     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1398     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1399     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1400   }
1401   return;
1402
1403 }
1404 //___________________________________________________________________________
1405
1406 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1407   // Return background integral in a range
1408
1409   if(fcounter==0) {
1410     cout<<"Use MassFitter method before Background"<<endl;
1411     return;
1412   }
1413
1414   //functions names
1415   TString bkgname="funcbkgRecalc";
1416   TString bkg1name="funcbkg1Recalc";
1417
1418   TF1 *funcbkg=0;
1419   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1420   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1421   if(!funcbkg){
1422     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1423     return;
1424   }
1425
1426
1427   Double_t intB,intBerr;
1428
1429   //relative error evaluation: from final parameters of the fit
1430   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1431   else{
1432     intB=funcbkg->GetParameter(0);
1433     intBerr=funcbkg->GetParError(0);
1434     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1435   }
1436
1437   //relative error evaluation: from histo
1438
1439   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1440   Double_t sum2=0;
1441   for(Int_t i=1;i<=fSideBandl;i++){
1442     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1443   }
1444   for(Int_t i=fSideBandr;i<=fNbin;i++){
1445     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1446   }
1447
1448   intBerr=TMath::Sqrt(sum2);
1449   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1450
1451   cout<<"Last estimation of bkg error is used"<<endl;
1452
1453   //backround +/- error in the range
1454   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1455     background = 0;
1456     errbackground = 0;
1457   }
1458   else{
1459     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1460     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1461     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1462   }
1463   return;
1464
1465 }
1466
1467
1468 //__________________________________________________________________________
1469
1470 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1471   // Return significance in mean+- n sigma
1472
1473   Double_t signal,errsignal,background,errbackground;
1474   Signal(nOfSigma,signal,errsignal);
1475   Background(nOfSigma,background,errbackground);
1476
1477   significance =  signal/TMath::Sqrt(signal+background);
1478
1479   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1480
1481   return;
1482 }
1483
1484 //__________________________________________________________________________
1485
1486 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1487   // Return significance integral in a range
1488
1489   Double_t signal,errsignal,background,errbackground;
1490   Signal(min, max,signal,errsignal);
1491   Background(min, max,background,errbackground);
1492
1493   significance =  signal/TMath::Sqrt(signal+background);
1494
1495   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1496
1497   return;
1498 }
1499
1500