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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliHFMassFitter for the fit of invariant mass distribution
  21 // of charmed mesons
  22 //
  23 // Author: C.Bianchin, chiara.bianchin@pd.infn.it
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <Riostream.h>
  27 #include <TMath.h>
  28 #include <TNtuple.h>
  29 #include <TH1F.h>
  30 #include <TF1.h>
  31 #include <TList.h>
  32 #include <TFile.h>
  33 #include <TCanvas.h>
  34 #include <TVirtualPad.h>
  35 #include <TGraph.h>
  36 #include <TVirtualFitter.h>
  37 #include <TMinuit.h>
  38 #include <TStyle.h>
  39 #include <TPaveText.h>
  40 #include <TDatabasePDG.h>
  41
  42 #include "AliHFMassFitter.h"
  43
  44
  45
  46 ClassImp(AliHFMassFitter)
  47
  48  //************** constructors
  49 AliHFMassFitter::AliHFMassFitter() :
  50   TNamed(),
  51   fhistoInvMass(0),
  52   fminMass(0),
  53   fmaxMass(0),
  54   fminBinMass(0),
  55   fmaxBinMass(0),
  56   fNbin(0),
  57   fParsSize(1),
  58   fNFinalPars(1),
  59   fFitPars(0),
  60   fWithBkg(0),
  61   ftypeOfFit4Bkg(0),
  62   ftypeOfFit4Sgn(0),
  63   ffactor(1),
  64   fntuParam(0),
  65   fMass(1.865),
  66   fSigmaSgn(0.012),
  67   fSideBands(0),
  68   fFixPar(0),
  69   fSideBandl(0),
  70   fSideBandr(0),
  71   fcounter(0),
  72   fContourGraph(0)
  73 {
  74   // default constructor
  75
  76   cout<<"Default constructor:"<<endl;
  77   cout<<"Remember to set the Histo, the Type, the FixPar"<<endl;
  78
  79 }
  80
  81 //___________________________________________________________________________
  82
  83 AliHFMassFitter::AliHFMassFitter (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
  84  TNamed(),
  85  fhistoInvMass(0),
  86  fminMass(0),
  87  fmaxMass(0),
  88  fminBinMass(0),
  89  fmaxBinMass(0),
  90  fNbin(0),
  91  fParsSize(1),
  92  fNFinalPars(1),
  93  fFitPars(0),
  94  fWithBkg(0),
  95  ftypeOfFit4Bkg(0),
  96  ftypeOfFit4Sgn(0),
  97  ffactor(1),
  98  fntuParam(0),
  99  fMass(1.865),
 100  fSigmaSgn(0.012),
 101  fSideBands(0),
 102  fFixPar(0),
 103  fSideBandl(0),
 104  fSideBandr(0),
 105  fcounter(0),
 106  fContourGraph(0)
 107 {
 108   // standard constructor
 109
 110   fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
 111   fhistoInvMass->SetDirectory(0);
 112   fminMass=minvalue;
 113   fmaxMass=maxvalue;
 114   if(rebin!=1) RebinMass(rebin);
 115   else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
 116   CheckRangeFit();
 117   ftypeOfFit4Bkg=fittypeb;
 118   ftypeOfFit4Sgn=fittypes;
 119   if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
 120   else fWithBkg=kTRUE;
 121   if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
 122   else  cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
 123
 124   ComputeParSize();
 125   fFitPars=new Float_t[fParsSize];
 126
 127   SetDefaultFixParam();
 128
 129   fContourGraph=new TList();
 130   fContourGraph->SetOwner();
 131
 132 }
 133
 134
 135 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
 136   TNamed(),
 137   fhistoInvMass(mfit.fhistoInvMass),
 138   fminMass(mfit.fminMass),
 139   fmaxMass(mfit.fmaxMass),
 140   fminBinMass(mfit.fminBinMass),
 141   fmaxBinMass(mfit.fmaxBinMass),
 142   fNbin(mfit.fNbin),
 143   fParsSize(mfit.fParsSize),
 144   fNFinalPars(mfit.fNFinalPars),
 145   fFitPars(0),
 146   fWithBkg(mfit.fWithBkg),
 147   ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
 148   ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
 149   ffactor(mfit.ffactor),
 150   fntuParam(mfit.fntuParam),
 151   fMass(mfit.fMass),
 152   fSigmaSgn(mfit.fSigmaSgn),
 153   fSideBands(mfit.fSideBands),
 154   fFixPar(0),
 155   fSideBandl(mfit.fSideBandl),
 156   fSideBandr(mfit.fSideBandr),
 157   fcounter(mfit.fcounter),
 158   fContourGraph(mfit.fContourGraph)
 159 {
 160   //copy constructor
 161
 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));
 167   }
 168
 169 }
 170
 171 //_________________________________________________________________________
 172
 173 AliHFMassFitter::~AliHFMassFitter() {
 174
 175   //destructor
 176
 177   cout<<"AliHFMassFitter destructor called"<<endl;
 178
 179   delete fhistoInvMass;
 180
 181   delete fntuParam;
 182
 183   delete[] fFitPars;
 184
 185   delete[] fFixPar;
 186
 187   fcounter = 0;
 188 }
 189
 190 //_________________________________________________________________________
 191
 192 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
 193
 194   //assignment operator
 195
 196   if(&mfit == this) return *this;
 197   fhistoInvMass= mfit.fhistoInvMass;
 198   fminMass= mfit.fminMass;
 199   fmaxMass= mfit.fmaxMass;
 200   fNbin= mfit.fNbin;
 201   fParsSize= mfit.fParsSize;
 202   fWithBkg= mfit.fWithBkg;
 203   ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
 204   ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
 205   ffactor= mfit.ffactor;
 206   fntuParam= mfit.fntuParam;
 207   fMass= mfit.fMass;
 208   fSigmaSgn= mfit.fSigmaSgn;
 209   fSideBands= mfit.fSideBands;
 210   fSideBandl= mfit.fSideBandl;
 211   fSideBandr= mfit.fSideBandr;
 212   fcounter= mfit.fcounter;
 213   fContourGraph= mfit.fContourGraph;
 214
 215   if(mfit.fParsSize > 0){
 216     delete[] fFitPars;
 217
 218     fFitPars=new Float_t[fParsSize];
 219     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 220
 221     delete[] fFixPar;
 222
 223     fFixPar=new Bool_t[fNFinalPars];
 224     memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Float_t));
 225   }
 226
 227   return *this;
 228 }
 229
 230 //************ tools & settings
 231
 232 //__________________________________________________________________________
 233
 234 void AliHFMassFitter::ComputeNFinalPars() {
 235
 236   //compute the number of parameters of the total (signal+bgk) function
 237   cout<<"Info:ComputeNFinalPars... ";
 238   switch (ftypeOfFit4Bkg) {//npar background func
 239   case 0:
 240     fNFinalPars=2;
 241     break;
 242   case 1:
 243     fNFinalPars=2;
 244     break;
 245   case 2:
 246     fNFinalPars=3;
 247     break;
 248   case 3:
 249     fNFinalPars=1;
 250   case 4:
 251     fNFinalPars=2;
 252     break;
 253   case 5:
 254     fNFinalPars=3;
 255     break;
 256   default:
 257     cout<<"Error in computing fNFinalPars: check ftypeOfFit4Bkg"<<endl;
 258     break;
 259   }
 260
 261   fNFinalPars+=3; //gaussian signal
 262   cout<<": "<<fNFinalPars<<endl;
 263 }
 264 //__________________________________________________________________________
 265
 266 void AliHFMassFitter::ComputeParSize() {
 267
 268   //compute the size of the parameter array and set the data member
 269
 270   switch (ftypeOfFit4Bkg) {//npar background func
 271   case 0:
 272     fParsSize = 2*3;
 273     break;
 274   case 1:
 275     fParsSize = 2*3;
 276     break;
 277   case 2:
 278     fParsSize = 3*3;
 279     break;
 280   case 3:
 281     fParsSize = 1*3;
 282     break;
 283   case 4:
 284     fParsSize = 2*3;
 285     break;
 286   case 5:
 287     fParsSize = 3*3;
 288     break;
 289   default:
 290     cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
 291     break;
 292   }
 293
 294   fParsSize += 3; // npar refl
 295   fParsSize += 3; // npar signal gaus
 296
 297   fParsSize*=2;   // add errors
 298   cout<<"Parameters array size "<<fParsSize<<endl;
 299 }
 300
 301 //___________________________________________________________________________
 302 void AliHFMassFitter::SetDefaultFixParam(){
 303
 304   //Set default values for fFixPar (only total integral fixed)
 305
 306   ComputeNFinalPars();
 307   fFixPar=new Bool_t[fNFinalPars];
 308
 309   fFixPar[0]=kTRUE; //default: IntTot fixed
 310   cout<<"Parameter 0 is fixed"<<endl;
 311   for(Int_t i=1;i<fNFinalPars;i++){
 312     fFixPar[i]=kFALSE;
 313   }
 314
 315 }
 316
 317 //___________________________________________________________________________
 318 Bool_t AliHFMassFitter::SetFixThisParam(Int_t thispar,Bool_t fixpar){
 319
 320   //set the value (kFALSE or kTRUE) of one element of fFixPar
 321   //return kFALSE if something wrong
 322
 323   if(thispar>=fNFinalPars) {
 324     cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
 325     return kFALSE;
 326   }
 327   if(!fFixPar){
 328     cout<<"Initializing fFixPar...";
 329     SetDefaultFixParam();
 330     cout<<" done."<<endl;
 331   }
 332
 333   fFixPar[thispar]=fixpar;
 334   if(fixpar)cout<<"Parameter "<<thispar<<" is now fixed"<<endl;
 335   else cout<<"Parameter "<<thispar<<" is now free"<<endl;
 336   return kTRUE;
 337 }
 338
 339 //___________________________________________________________________________
 340 Bool_t AliHFMassFitter::GetFixThisParam(Int_t thispar)const{
 341   //return the value of fFixPar[thispar]
 342   if(thispar>=fNFinalPars) {
 343     cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
 344     return kFALSE;
 345   }
 346   if(!fFixPar) {
 347     cout<<"Error! Parameters to be fixed still not set"<<endl;
 348     return kFALSE;
 349
 350   }
 351   return fFixPar[thispar];
 352
 353 }
 354
 355 //___________________________________________________________________________
 356 void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
 357
 358   fhistoInvMass = new TH1F(*histoToFit);
 359   fhistoInvMass->SetDirectory(0);
 360   //cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
 361 }
 362
 363 //___________________________________________________________________________
 364
 365 void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
 366
 367   //set the type of fit to perform for signal and background
 368
 369   ftypeOfFit4Bkg = fittypeb;
 370   ftypeOfFit4Sgn = fittypes;
 371
 372   ComputeParSize();
 373   fFitPars = new Float_t[fParsSize];
 374
 375   SetDefaultFixParam();
 376 }
 377
 378 //___________________________________________________________________________
 379
 380 void AliHFMassFitter::Reset() {
 381
 382   //delete the histogram and reset the mean and sigma to default
 383
 384   cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
 385   fMass=1.85;
 386   fSigmaSgn=0.012;
 387   cout<<"Reset "<<fhistoInvMass<<endl;
 388   delete fhistoInvMass;
 389 }
 390
 391 //_________________________________________________________________________
 392
 393 void AliHFMassFitter::InitNtuParam(char *ntuname) {
 394
 395   // Create ntuple to keep fit parameters
 396
 397   fntuParam=0;
 398   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");
 399
 400 }
 401
 402 //_________________________________________________________________________
 403
 404 void AliHFMassFitter::FillNtuParam() {
 405   // Fill ntuple with fit parameters
 406
 407   Float_t nothing=0.;
 408
 409   if (ftypeOfFit4Bkg==2) {
 410       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 411       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 412       fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
 413       fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
 414       fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
 415       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
 416       fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
 417       fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
 418       fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
 419       fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
 420       fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
 421       fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
 422       fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
 423       fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
 424       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
 425
 426       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
 427       fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
 428       fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
 429       fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
 430       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
 431       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
 432       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
 433       fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
 434       fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
 435       fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
 436       fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
 437       fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
 438       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
 439       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
 440       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
 441
 442   } else {
 443
 444     if(ftypeOfFit4Bkg==3){
 445       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 446       fntuParam->SetBranchAddress("slope1",&nothing);
 447       fntuParam->SetBranchAddress("conc1",&nothing);
 448       fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
 449       fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
 450       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
 451       fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
 452       fntuParam->SetBranchAddress("slope2",&nothing);
 453       fntuParam->SetBranchAddress("conc2",&nothing);
 454       fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
 455       fntuParam->SetBranchAddress("slope3",&nothing);
 456       fntuParam->SetBranchAddress("conc3",&nothing);
 457       fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
 458       fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
 459       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
 460
 461       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
 462       fntuParam->SetBranchAddress("slope1Err",&nothing);
 463       fntuParam->SetBranchAddress("conc1Err",&nothing);
 464       fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
 465       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
 466       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
 467       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
 468       fntuParam->SetBranchAddress("slope2Err",&nothing);
 469       fntuParam->SetBranchAddress("conc2Err",&nothing);
 470       fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
 471       fntuParam->SetBranchAddress("slope3Err",&nothing);
 472       fntuParam->SetBranchAddress("conc3Err",&nothing);
 473       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
 474       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
 475       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
 476
 477     }
 478     else{
 479       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 480       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 481       fntuParam->SetBranchAddress("conc1",&nothing);
 482       fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
 483       fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
 484       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
 485       fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
 486       fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
 487       fntuParam->SetBranchAddress("conc2",&nothing);
 488       fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
 489       fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
 490       fntuParam->SetBranchAddress("conc3",&nothing);
 491       fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
 492       fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
 493       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
 494
 495       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
 496       fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
 497       fntuParam->SetBranchAddress("conc1Err",&nothing);
 498       fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
 499       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
 500       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
 501       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
 502       fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
 503       fntuParam->SetBranchAddress("conc2Err",&nothing);
 504       fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
 505       fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
 506       fntuParam->SetBranchAddress("conc3Err",&nothing);
 507       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
 508       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
 509       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
 510     }
 511
 512   }
 513   fntuParam->TTree::Fill();
 514 }
 515
 516 //_________________________________________________________________________
 517
 518 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
 519   // Create, fill and return ntuple with fit parameters
 520
 521   InitNtuParam(ntuname);
 522   FillNtuParam();
 523   return fntuParam;
 524 }
 525 //_________________________________________________________________________
 526
 527 void AliHFMassFitter::RebinMass(Int_t bingroup){
 528   // Rebin invariant mass histogram
 529
 530   if(!fhistoInvMass){
 531     cout<<"Histogram not set"<<endl;
 532     return;
 533   }
 534   Int_t nbinshisto=fhistoInvMass->GetNbinsX();
 535   if(bingroup<1){
 536     cout<<"Error! Cannot group "<<bingroup<<" bins\n";
 537     fNbin=nbinshisto;
 538     cout<<"Kept original number of bins: "<<fNbin<<endl;
 539   } else{
 540
 541     while(nbinshisto%bingroup != 0) {
 542       bingroup--;
 543     }
 544     cout<<"Group "<<bingroup<<" bins"<<endl;
 545     fhistoInvMass->Rebin(bingroup);
 546     fNbin = fhistoInvMass->GetNbinsX();
 547     cout<<"New number of bins: "<<fNbin<<endl;
 548   }
 549
 550 }
 551
 552 //************* fit
 553
 554 //___________________________________________________________________________
 555
 556 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
 557   // Fit function for signal+background
 558
 559
 560   //exponential or linear fit
 561   //
 562   // par[0] = tot integral
 563   // par[1] = slope
 564   // par[2] = gaussian integral
 565   // par[3] = gaussian mean
 566   // par[4] = gaussian sigma
 567
 568   Double_t total,bkg=0,sgn=0;
 569
 570   if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
 571     if(ftypeOfFit4Sgn == 0) {
 572
 573       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 574       bkg = FitFunction4Bkg(x,parbkg);
 575     }
 576     if(ftypeOfFit4Sgn == 1) {
 577       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
 578       bkg = FitFunction4Bkg(x,parbkg);
 579     }
 580
 581     sgn = FitFunction4Sgn(x,&par[2]);
 582
 583   }
 584
 585   //polynomial fit
 586
 587     // par[0] = tot integral
 588     // par[1] = coef1
 589     // par[2] = coef2
 590     // par[3] = gaussian integral
 591     // par[4] = gaussian mean
 592     // par[5] = gaussian sigma
 593
 594   if (ftypeOfFit4Bkg==2) {
 595
 596     if(ftypeOfFit4Sgn == 0) {
 597
 598       Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
 599       bkg = FitFunction4Bkg(x,parbkg);
 600     }
 601     if(ftypeOfFit4Sgn == 1) {
 602
 603       Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
 604      bkg = FitFunction4Bkg(x,parbkg);
 605     }
 606
 607     sgn = FitFunction4Sgn(x,&par[3]);
 608   }
 609
 610   if (ftypeOfFit4Bkg==3) {
 611
 612     if(ftypeOfFit4Sgn == 0) {
 613         bkg=FitFunction4Bkg(x,par);
 614         sgn=FitFunction4Sgn(x,&par[1]);
 615     }
 616     if(ftypeOfFit4Sgn == 1) {
 617       Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
 618       bkg=FitFunction4Bkg(x,parbkg);
 619       sgn=FitFunction4Sgn(x,&par[1]);
 620     }
 621   }
 622
 623   //Power fit
 624
 625     // par[0] = tot integral
 626     // par[1] = coef1
 627     // par[2] = gaussian integral
 628     // par[3] = gaussian mean
 629     // par[4] = gaussian sigma
 630
 631   if (ftypeOfFit4Bkg==4) {
 632
 633     if(ftypeOfFit4Sgn == 0) {
 634
 635       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 636       bkg = FitFunction4Bkg(x,parbkg);
 637     }
 638     if(ftypeOfFit4Sgn == 1) {
 639
 640       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-par[2], par[1]};
 641       bkg = FitFunction4Bkg(x,parbkg);
 642     }
 643     sgn = FitFunction4Sgn(x,&par[2]);
 644   }
 645
 646
 647   //Power and exponential fit
 648
 649     // par[0] = tot integral
 650     // par[1] = coef1
 651     // par[2] = coef2
 652     // par[3] = gaussian integral
 653     // par[4] = gaussian mean
 654     // par[5] = gaussian sigma
 655
 656   if (ftypeOfFit4Bkg==5) {
 657
 658     if(ftypeOfFit4Sgn == 0) {
 659         Double_t parbkg[3] = {par[0]-par[3],par[1],par[2]};
 660         bkg = FitFunction4Bkg(x,parbkg);
 661     }
 662     if(ftypeOfFit4Sgn == 1) {
 663      Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-par[3], par[1], par[2]};
 664      bkg = FitFunction4Bkg(x,parbkg);
 665     }
 666     sgn = FitFunction4Sgn(x,&par[3]);
 667   }
 668
 669   total = bkg + sgn;
 670
 671   return  total;
 672 }
 673
 674 //_________________________________________________________________________
 675 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
 676   // Fit function for the signal
 677
 678   //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 679   //Par:
 680   // * [0] = integralSgn
 681   // * [1] = mean
 682   // * [2] = sigma
 683   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 684
 685   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]);
 686
 687 }
 688
 689 //__________________________________________________________________________
 690
 691 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
 692   // Fit function for the background
 693
 694   Double_t maxDeltaM = 4.*fSigmaSgn;
 695   if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
 696     TF1::RejectPoint();
 697     return 0;
 698   }
 699   Int_t firstPar=0;
 700   Double_t gaus2=0,total=-1;
 701   if(ftypeOfFit4Sgn == 1){
 702     firstPar=3;
 703     //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 704     //Par:
 705     // * [0] = integralSgn
 706     // * [1] = mean
 707     // * [2] = sigma
 708     //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 709     gaus2 = FitFunction4Sgn(x,par);
 710   }
 711
 712   switch (ftypeOfFit4Bkg){
 713   case 0:
 714     //exponential
 715     //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
 716     //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
 717     //as integralTot- integralGaus (=par [2])
 718     //Par:
 719     // * [0] = integralBkg;
 720     // * [1] = B;
 721     //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
 722     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]);
 723     break;
 724   case 1:
 725     //linear
 726     //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)
 727     // * [0] = integralBkg;
 728     // * [1] = b;
 729     total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
 730     break;
 731   case 2:
 732     //polynomial
 733     //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
 734     //+ 1/3*c*(max^3-min^3) ->
 735     //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
 736     // * [0] = integralBkg;
 737     // * [1] = b;
 738     // * [2] = c;
 739     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));
 740     break;
 741   case 3:
 742     total=par[0+firstPar];
 743     break;
 744   case 4:
 745     //power function
 746     //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
 747     //
 748     //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
 749     // * [0] = integralBkg;
 750     // * [1] = b;
 751     // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
 752     {
 753     Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 754
 755     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]);
 756     }
 757     break;
 758   case 5:
 759    //power function wit exponential
 760     //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
 761     {
 762     Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 763
 764     total = par[1+firstPar]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2+firstPar]*(x[0]-mpi));
 765     }
 766     break;
 767 //   default:
 768 //     Types of Fit Functions for Background:
 769 //     * 0 = exponential;
 770 //     * 1 = linear;
 771 //     * 2 = polynomial 2nd order
 772 //     * 3 = no background"<<endl;
 773 //     * 4 = Power function
 774 //     * 5 = Power function with exponential
 775
 776   }
 777   return total+gaus2;
 778 }
 779
 780 //__________________________________________________________________________
 781 Bool_t AliHFMassFitter::SideBandsBounds(){
 782
 783   //determines the ranges of the side bands
 784
 785   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 786   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 787   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin+1);
 788
 789   Double_t sidebandldouble,sidebandrdouble;
 790   Bool_t leftok=kFALSE, rightok=kFALSE;
 791
 792   if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
 793     cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
 794     return kFALSE;
 795   }
 796
 797   //histo limit = fit function limit
 798   if((TMath::Abs(fminMass-minHisto) < 1e6 || TMath::Abs(fmaxMass - maxHisto) < 1e6) && (fMass-4.*fSigmaSgn-fminMass) < 1e6){
 799     Double_t coeff = (fMass-fminMass)/fSigmaSgn;
 800     sidebandldouble=(fMass-0.5*coeff*fSigmaSgn);
 801     sidebandrdouble=(fMass+0.5*coeff*fSigmaSgn);
 802     cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
 803     if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
 804     if (coeff<2) {
 805       cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
 806       ftypeOfFit4Bkg=3;
 807       //set binleft and right without considering SetRangeFit- anyway no bkg!
 808       sidebandldouble=(fMass-4.*fSigmaSgn);
 809       sidebandrdouble=(fMass+4.*fSigmaSgn);
 810     }
 811   }
 812   else {
 813     sidebandldouble=(fMass-4.*fSigmaSgn);
 814     sidebandrdouble=(fMass+4.*fSigmaSgn);
 815   }
 816
 817   cout<<"Left side band ";
 818   Double_t tmp=0.;
 819   tmp=sidebandldouble;
 820   //calculate bin corresponding to fSideBandl
 821   fSideBandl=fhistoInvMass->FindBin(sidebandldouble);
 822   if (sidebandldouble >= fhistoInvMass->GetBinCenter(fSideBandl)) fSideBandl++;
 823   sidebandldouble=fhistoInvMass->GetBinLowEdge(fSideBandl);
 824
 825   if(TMath::Abs(tmp-sidebandldouble) > 1e-6){
 826     cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
 827     leftok=kTRUE;
 828   }
 829   cout<<sidebandldouble<<" (bin "<<fSideBandl<<")"<<endl;
 830
 831   cout<<"Right side band ";
 832   tmp=sidebandrdouble;
 833   //calculate bin corresponding to fSideBandr
 834   fSideBandr=fhistoInvMass->FindBin(sidebandrdouble);
 835   if (sidebandrdouble < fhistoInvMass->GetBinCenter(fSideBandr)) fSideBandr--;
 836   sidebandrdouble=fhistoInvMass->GetBinLowEdge(fSideBandr+1);
 837
 838   if(TMath::Abs(tmp-sidebandrdouble) > 1e-6){
 839     cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
 840     rightok=kTRUE;
 841   }
 842   cout<<sidebandrdouble<<" (bin "<<fSideBandr<<")"<<endl;
 843   if (fSideBandl==0 || fSideBandr==fNbin) {
 844     cout<<"Error! Range too little";
 845     return kFALSE;
 846   }
 847   return kTRUE;
 848 }
 849
 850 //__________________________________________________________________________
 851
 852 void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
 853
 854   // get the range of the side bands
 855
 856   if (fSideBandl==0 && fSideBandr==0){
 857     cout<<"Use MassFitter method first"<<endl;
 858     return;
 859   }
 860   left=fSideBandl;
 861   right=fSideBandr;
 862 }
 863
 864 //__________________________________________________________________________
 865 Bool_t AliHFMassFitter::CheckRangeFit(){
 866   //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
 867
 868   if (!fhistoInvMass) {
 869     cout<<"No histogram to fit! SetHisto(TH1F* h) before! "<<endl;
 870     return kFALSE;
 871   }
 872   Bool_t leftok=kFALSE, rightok=kFALSE;
 873   Int_t nbins=fhistoInvMass->GetNbinsX();
 874   Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins+1);
 875
 876   //check if limits are inside histogram range
 877
 878   if( fminMass-minhisto < 0. ) {
 879     cout<<"Out of histogram left bound! Setting to "<<minhisto<<endl;
 880     fminMass=minhisto;
 881   }
 882   if( fmaxMass-maxhisto > 0. ) {
 883     cout<<"Out of histogram right bound! Setting to"<<maxhisto<<endl;
 884     fmaxMass=maxhisto;
 885   }
 886
 887   Double_t tmp=0.;
 888   tmp=fminMass;
 889   //calculate bin corresponding to fminMass
 890   fminBinMass=fhistoInvMass->FindBin(fminMass);
 891   if (fminMass >= fhistoInvMass->GetBinCenter(fminBinMass)) fminBinMass++;
 892   fminMass=fhistoInvMass->GetBinLowEdge(fminBinMass);
 893   if(TMath::Abs(tmp-fminMass) > 1e-6){
 894     cout<<"Left bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fminMass<<endl;
 895     leftok=kTRUE;
 896   }
 897
 898   tmp=fmaxMass;
 899   //calculate bin corresponding to fmaxMass
 900   fmaxBinMass=fhistoInvMass->FindBin(fmaxMass);
 901   if (fmaxMass < fhistoInvMass->GetBinCenter(fmaxBinMass)) fmaxBinMass--;
 902   fmaxMass=fhistoInvMass->GetBinLowEdge(fmaxBinMass+1);
 903   if(TMath::Abs(tmp-fmaxMass) > 1e-6){
 904     cout<<"Right bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fmaxMass<<endl;
 905     rightok=kTRUE;
 906   }
 907
 908   return (leftok && rightok);
 909
 910 }
 911
 912 //__________________________________________________________________________
 913
 914 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 915   // Main method of the class: performs the fit of the histogram
 916
 917   //Set default fitter Minuit in order to use gMinuit in the contour plots
 918   TVirtualFitter::SetDefaultFitter("Minuit");
 919
 920
 921   Bool_t isBkgOnly=kFALSE;
 922
 923   Int_t fit1status=RefitWithBkgOnly(kFALSE);
 924   if(fit1status){
 925     Int_t checkinnsigma=4;
 926     Double_t range[2]={fMass-checkinnsigma*fSigmaSgn,fMass+checkinnsigma*fSigmaSgn};
 927     TF1* func=GetHistoClone()->GetFunction("funcbkgonly");
 928     Double_t intUnderFunc=func->Integral(range[0],range[1]);
 929     Double_t intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(range[0]),fhistoInvMass->FindBin(range[1]),"width");
 930     cout<<"Pick zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
 931     Double_t diffUnderPick=(intUnderHisto-intUnderFunc);
 932     intUnderFunc=func->Integral(fminMass,fminMass+checkinnsigma*fSigmaSgn);
 933     intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fminMass+checkinnsigma*fSigmaSgn),"width");
 934     cout<<"Band (l) zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
 935     Double_t diffUnderBands=(intUnderHisto-intUnderFunc);
 936     Double_t relDiff=diffUnderPick/diffUnderBands;
 937     cout<<"Relative difference = "<<relDiff<<endl;
 938     if(TMath::Abs(relDiff) < 1) isBkgOnly=kTRUE;
 939     else{
 940       cout<<"Relative difference = "<<relDiff<<": I suppose there is some signal, continue with total fit!"<<endl;
 941     }
 942   }
 943   if(isBkgOnly) {
 944     cout<<"INFO!! The histogram contains only background"<<endl;
 945     if(draw)DrawFit();
 946
 947     //increase counter of number of fits done
 948     fcounter++;
 949
 950     return kTRUE;
 951   }
 952
 953   Int_t bkgPar = fNFinalPars-3; //background function's number of parameters
 954
 955   cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
 956
 957
 958   TString listname="contourplot";
 959   listname+=fcounter;
 960   if(!fContourGraph){
 961     fContourGraph=new TList();
 962     fContourGraph->SetOwner();
 963   }
 964
 965   fContourGraph->SetName(listname);
 966
 967
 968   //function names
 969   TString bkgname="funcbkg";
 970   TString bkg1name="funcbkg1";
 971   TString massname="funcmass";
 972
 973   //Total integral
 974   Double_t totInt = fhistoInvMass->Integral(fminBinMass,fmaxBinMass, "width");
 975   //cout<<"Here tot integral is = "<<totInt<<"; integral in whole range is "<<fhistoInvMass->Integral("width")<<endl;
 976   fSideBands = kTRUE;
 977   Double_t width=fhistoInvMass->GetBinWidth(8);
 978   //cout<<"fNbin = "<<fNbin<<endl;
 979   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 980
 981   Bool_t ok=SideBandsBounds();
 982   if(!ok) return kFALSE;
 983
 984   //sidebands integral - first approx (from histo)
 985   Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
 986   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 987   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 988   if (sideBandsInt<=0) {
 989     cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
 990     return kFALSE;
 991   }
 992
 993   /*Fit Bkg*/
 994
 995
 996   TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 997   cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
 998
 999   funcbkg->SetLineColor(2); //red
1000
1001   //first fit for bkg: approx bkgint
1002
1003   switch (ftypeOfFit4Bkg) {
1004   case 0: //gaus+expo
1005     funcbkg->SetParNames("BkgInt","Slope");
1006     funcbkg->SetParameters(sideBandsInt,-2.);
1007     break;
1008   case 1:
1009     funcbkg->SetParNames("BkgInt","Slope");
1010     funcbkg->SetParameters(sideBandsInt,-100.);
1011     break;
1012   case 2:
1013     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1014     funcbkg->SetParameters(sideBandsInt,-10.,5);
1015     break;
1016   case 3:
1017     if(ftypeOfFit4Sgn==0){
1018       funcbkg->SetParNames("Const");
1019       funcbkg->SetParameter(0,0.);
1020       funcbkg->FixParameter(0,0.);
1021     }
1022     break;
1023   case 4:
1024     funcbkg->SetParNames("BkgInt","Coef2");
1025     funcbkg->SetParameters(sideBandsInt,0.5);
1026     break;
1027  case 5:
1028     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1029     funcbkg->SetParameters(sideBandsInt, 0.5, 50.);
1030     break;
1031   default:
1032     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1033     return kFALSE;
1034     break;
1035   }
1036   cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
1037   Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
1038
1039   Double_t intbkg1=0,slope1=0,conc1=0;
1040   //if only signal and reflection: skip
1041   if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
1042     ftypeOfFit4Sgn=0;
1043     fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
1044
1045     for(Int_t i=0;i<bkgPar;i++){
1046       fFitPars[i]=funcbkg->GetParameter(i);
1047       //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
1048       fFitPars[fNFinalPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
1049       //cout<<fNFinalPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
1050     }
1051     fSideBands = kFALSE;
1052     //intbkg1 = funcbkg->GetParameter(0);
1053
1054     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
1055     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
1056     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
1057     if(ftypeOfFit4Bkg==5) conc1 = funcbkg->GetParameter(2);
1058
1059
1060     //cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
1061   }
1062   else cout<<"\t\t//"<<endl;
1063
1064   ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
1065   TF1 *funcbkg1=0;
1066   if (ftypeOfFit4Sgn == 1) {
1067     cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
1068     bkgPar+=3;
1069
1070     //cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
1071
1072     funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
1073     cout<<"Function name = "<<funcbkg1->GetName()<<endl;
1074
1075     funcbkg1->SetLineColor(2); //red
1076
1077  switch (ftypeOfFit4Bkg) {
1078     case 0:
1079         {
1080         cout<<"*** Exponential Fit ***"<<endl;
1081         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
1082         funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1083         }
1084         break;
1085      case 1:
1086         {
1087         cout<<"*** Linear Fit ***"<<endl;
1088         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
1089         funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1090         }
1091         break;
1092      case 2:
1093         {
1094         cout<<"*** Polynomial Fit ***"<<endl;
1095         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
1096         funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
1097         }
1098         break;
1099     case 3:
1100        //no background: gaus sign+ gaus broadened
1101         {
1102         cout<<"*** No background Fit ***"<<endl;
1103         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
1104         funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
1105         funcbkg1->FixParameter(3,0.);
1106         }
1107         break;
1108      case 4:
1109         {
1110         cout<<"*** Power function Fit ***"<<endl;
1111         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef2");
1112         funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1113         }
1114         break;
1115       case 5:
1116         {
1117         cout<<"*** Power function conv. with exponential Fit ***"<<endl;
1118         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
1119         funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
1120         }
1121         break;
1122     }
1123       //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
1124       //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
1125
1126     Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
1127     if (status != 0){
1128       cout<<"Minuit returned "<<status<<endl;
1129       return kFALSE;
1130     }
1131
1132     for(Int_t i=0;i<bkgPar;i++){
1133       fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
1134       //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
1135       fFitPars[fNFinalPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
1136       //cout<<"\t"<<fNFinalPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
1137     }
1138
1139     intbkg1=funcbkg1->GetParameter(3);
1140     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
1141     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
1142     if(ftypeOfFit4Bkg==5) conc1 = funcbkg1->GetParameter(5);
1143
1144
1145   } else {
1146     bkgPar+=3;
1147
1148     for(Int_t i=0;i<3;i++){
1149       fFitPars[bkgPar-3+i]=0.;
1150       cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
1151       fFitPars[fNFinalPars+3*bkgPar-6+i]= 0.;
1152       cout<<fNFinalPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
1153     }
1154
1155     for(Int_t i=0;i<bkgPar-3;i++){
1156       fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
1157       cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
1158       fFitPars[fNFinalPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
1159       cout<<fNFinalPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
1160     }
1161
1162
1163   }
1164
1165   //sidebands integral - second approx (from fit)
1166   fSideBands = kFALSE;
1167   Double_t bkgInt;
1168   //cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
1169   if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
1170   else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
1171   //cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
1172
1173   //Signal integral - first approx
1174   Double_t sgnInt;
1175   sgnInt = totInt-bkgInt;
1176   //cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
1177   if (sgnInt <= 0){
1178     cout<<"Setting sgnInt = - sgnInt"<<endl;
1179     sgnInt=(-1)*sgnInt;
1180   }
1181   /*Fit All Mass distribution with exponential + gaussian (+gaussian braodened) */
1182   TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4MassDistr");
1183   cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
1184   funcmass->SetLineColor(4); //blue
1185
1186   //Set parameters
1187   cout<<"\nTOTAL FIT"<<endl;
1188
1189   if(fNFinalPars==5){
1190     funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
1191     funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
1192
1193     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1194     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1195     if(fFixPar[0]){
1196       funcmass->FixParameter(0,totInt);
1197     }
1198     if(fFixPar[1]){
1199       funcmass->FixParameter(1,slope1);
1200     }
1201     if(fFixPar[2]){
1202       funcmass->FixParameter(2,sgnInt);
1203     }
1204     if(fFixPar[3]){
1205       funcmass->FixParameter(3,fMass);
1206     }
1207     if(fFixPar[4]){
1208       funcmass->FixParameter(4,fSigmaSgn);
1209     }
1210   }
1211   if (fNFinalPars==6){
1212     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
1213     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
1214
1215     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1216     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1217     if(fFixPar[0])funcmass->FixParameter(0,totInt);
1218     if(fFixPar[1])funcmass->FixParameter(1,slope1);
1219     if(fFixPar[2])funcmass->FixParameter(2,conc1);
1220     if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
1221     if(fFixPar[4])funcmass->FixParameter(4,fMass);
1222     if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
1223     //
1224     //funcmass->FixParameter(2,sgnInt);
1225   }
1226   if(fNFinalPars==4){
1227     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
1228     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
1229     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
1230     if(fFixPar[0]) funcmass->FixParameter(0,0.);
1231     if(fFixPar[1])funcmass->FixParameter(1,sgnInt);
1232     if(fFixPar[2])funcmass->FixParameter(2,fMass);
1233     if(fFixPar[3])funcmass->FixParameter(3,fSigmaSgn);
1234    //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1235     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1236
1237   }
1238
1239   Int_t status;
1240
1241   status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
1242   if (status != 0){
1243     cout<<"Minuit returned "<<status<<endl;
1244     return kFALSE;
1245   }
1246
1247   cout<<"fit done"<<endl;
1248   //reset value of fMass and fSigmaSgn to those found from fit
1249   fMass=funcmass->GetParameter(fNFinalPars-2);
1250   fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
1251
1252   for(Int_t i=0;i<fNFinalPars;i++){
1253     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1254     fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1255     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1256   }
1257   /*
1258   //check: cout parameters
1259   for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
1260     cout<<i<<"\t"<<fFitPars[i]<<endl;
1261     }
1262   */
1263
1264   if(funcmass->GetParameter(fNFinalPars-1) <0 || funcmass->GetParameter(fNFinalPars-2) <0 || funcmass->GetParameter(fNFinalPars-3) <0 ) {
1265     cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1266     return kFALSE;
1267   }
1268
1269   //increase counter of number of fits done
1270   fcounter++;
1271
1272   //contour plots
1273   if(draw){
1274
1275     for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
1276
1277       for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
1278         cout<<"Par "<<kpar<<" and "<<jpar<<endl;
1279
1280         // produce 2 contours per couple of parameters
1281         TGraph* cont[2] = {0x0, 0x0};
1282         const Double_t errDef[2] = {1., 4.};
1283         for (Int_t i=0; i<2; i++) {
1284           gMinuit->SetErrorDef(errDef[i]);
1285           cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
1286           cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
1287         }
1288
1289         if(!cont[0] || !cont[1]){
1290           cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
1291           continue;
1292         }
1293
1294         // set graph titles and add them to the list
1295         TString title = "Contour plot";
1296         TString titleX = funcmass->GetParName(kpar);
1297         TString titleY = funcmass->GetParName(jpar);
1298         for (Int_t i=0; i<2; i++) {
1299           cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
1300           cont[i]->SetTitle(title);
1301           cont[i]->GetXaxis()->SetTitle(titleX);
1302           cont[i]->GetYaxis()->SetTitle(titleY);
1303           cont[i]->GetYaxis()->SetLabelSize(0.033);
1304           cont[i]->GetYaxis()->SetTitleSize(0.033);
1305           cont[i]->GetYaxis()->SetTitleOffset(1.67);
1306
1307           fContourGraph->Add(cont[i]);
1308         }
1309
1310         // plot them
1311         TString cvname = Form("c%d%d", kpar, jpar);
1312         TCanvas *c4=new TCanvas(cvname,cvname,600,600);
1313         c4->cd();
1314         cont[1]->SetFillColor(38);
1315         cont[1]->Draw("alf");
1316         cont[0]->SetFillColor(9);
1317         cont[0]->Draw("lf");
1318
1319       }
1320
1321     }
1322
1323   }
1324
1325   if (ftypeOfFit4Sgn == 1) {
1326     delete funcbkg1;
1327   }
1328   delete funcbkg;
1329   delete funcmass;
1330
1331   AddFunctionsToHisto();
1332   if (draw) DrawFit();
1333
1334
1335   return kTRUE;
1336 }
1337
1338 //______________________________________________________________________________
1339
1340 Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
1341
1342   //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
1343   //If you want to change the backgroud function or range use SetType or SetRangeFit before
1344
1345   TString bkgname="funcbkgonly";
1346   fSideBands = kFALSE;
1347
1348   TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1349
1350   funcbkg->SetLineColor(kBlue+3); //dark blue
1351
1352   Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
1353
1354   switch (ftypeOfFit4Bkg) {
1355   case 0: //gaus+expo
1356     funcbkg->SetParNames("BkgInt","Slope");
1357     funcbkg->SetParameters(integral,-2.);
1358     break;
1359   case 1:
1360     funcbkg->SetParNames("BkgInt","Slope");
1361     funcbkg->SetParameters(integral,-100.);
1362     break;
1363   case 2:
1364     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1365     funcbkg->SetParameters(integral,-10.,5);
1366     break;
1367   case 3:
1368     cout<<"Warning! This choice does not make a lot of sense..."<<endl;
1369     if(ftypeOfFit4Sgn==0){
1370       funcbkg->SetParNames("Const");
1371       funcbkg->SetParameter(0,0.);
1372       funcbkg->FixParameter(0,0.);
1373     }
1374     break;
1375   case 4:
1376     funcbkg->SetParNames("BkgInt","Coef1");
1377     funcbkg->SetParameters(integral,0.5);
1378     break;
1379   case 5:
1380     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1381     funcbkg->SetParameters(integral,-10.,5.);
1382     break;
1383   default:
1384     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1385     return kFALSE;
1386     break;
1387   }
1388
1389
1390   Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
1391   if (status != 0){
1392     cout<<"Minuit returned "<<status<<endl;
1393     return kFALSE;
1394   }
1395   AddFunctionsToHisto();
1396
1397   if(draw) DrawFit();
1398
1399   return kTRUE;
1400
1401 }
1402 //_________________________________________________________________________
1403 Double_t AliHFMassFitter::GetChiSquare() const{
1404   //Get Chi^2 method
1405   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1406   if(!funcmass) {
1407     cout<<"funcmass not found"<<endl;
1408     return -1;
1409   }
1410   return funcmass->GetChisquare();
1411 }
1412
1413 //_________________________________________________________________________
1414 Double_t AliHFMassFitter::GetReducedChiSquare() const{
1415   //Get reduced Chi^2 method
1416   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1417   if(!funcmass) {
1418     cout<<"funcmass not found"<<endl;
1419     return -1;
1420   }
1421   return funcmass->GetChisquare()/funcmass->GetNDF();
1422 }
1423
1424 //*********output
1425
1426 //_________________________________________________________________________
1427 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
1428   // Return fit parameters
1429
1430   for(Int_t i=0;i<fParsSize;i++){
1431     vector[i]=fFitPars[i];
1432   }
1433 }
1434
1435
1436 //_________________________________________________________________________
1437 void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1438
1439   //gives the integral of signal obtained from fit parameters
1440   if(!valuewitherror) {
1441     printf("AliHFMassFitter::IntS: got a null pointer\n");
1442     return;
1443   }
1444
1445   Int_t index=fParsSize/2 - 3;
1446   valuewitherror[0]=fFitPars[index];
1447   index=fParsSize - 3;
1448   valuewitherror[1]=fFitPars[index];
1449 }
1450
1451
1452 //_________________________________________________________________________
1453 void AliHFMassFitter::AddFunctionsToHisto(){
1454
1455   //Add the background function in the complete range to the list of functions attached to the histogram
1456
1457   //cout<<"AddFunctionsToHisto called"<<endl;
1458   TString bkgname = "funcbkg";
1459
1460   Bool_t done1=kFALSE,done2=kFALSE;
1461
1462   TString bkgnamesave=bkgname;
1463   TString testname=bkgname;
1464   testname += "FullRange";
1465   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1466   if(testfunc){
1467     done1=kTRUE;
1468     testfunc=0x0;
1469   }
1470   testname="funcbkgonly";
1471   testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1472   if(testfunc){
1473     done2=kTRUE;
1474     testfunc=0x0;
1475   }
1476
1477   if(done1 && done2){
1478     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1479     return;
1480   }
1481
1482   TList *hlist=fhistoInvMass->GetListOfFunctions();
1483   hlist->ls();
1484
1485   if(!done2){
1486     TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1487     if(!bonly){
1488       cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1489     }else{
1490       bonly->SetLineColor(kBlue+3);
1491       hlist->Add((TF1*)bonly->Clone());
1492       delete bonly;
1493     }
1494
1495   }
1496
1497   if(!done1){
1498     TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1499     if(!b){
1500       cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1501       return;
1502     }
1503
1504     bkgname += "FullRange";
1505     TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1506     //cout<<bfullrange->GetName()<<endl;
1507     for(Int_t i=0;i<fNFinalPars-3;i++){
1508       bfullrange->SetParName(i,b->GetParName(i));
1509       bfullrange->SetParameter(i,b->GetParameter(i));
1510       bfullrange->SetParError(i,b->GetParError(i));
1511     }
1512     bfullrange->SetLineStyle(4);
1513     bfullrange->SetLineColor(14);
1514
1515     bkgnamesave += "Recalc";
1516
1517     TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1518
1519     TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1520
1521     if (!mass){
1522       cout<<"funcmass doesn't exist "<<endl;
1523       return;
1524     }
1525
1526     //intBkg=intTot-intS
1527     blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
1528     blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
1529     if (fNFinalPars>=5) {
1530       blastpar->SetParameter(1,mass->GetParameter(1));
1531       blastpar->SetParError(1,mass->GetParError(1));
1532     }
1533     if (fNFinalPars==6) {
1534       blastpar->SetParameter(2,mass->GetParameter(2));
1535       blastpar->SetParError(2,mass->GetParError(2));
1536     }
1537
1538     blastpar->SetLineStyle(1);
1539     blastpar->SetLineColor(2);
1540
1541     hlist->Add((TF1*)bfullrange->Clone());
1542     hlist->Add((TF1*)blastpar->Clone());
1543     hlist->ls();
1544
1545     delete bfullrange;
1546     delete blastpar;
1547
1548   }
1549
1550
1551 }
1552
1553 //_________________________________________________________________________
1554
1555 TH1F* AliHFMassFitter::GetHistoClone() const{
1556
1557   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1558   return hout;
1559 }
1560 //_________________________________________________________________________
1561
1562 void AliHFMassFitter::WriteHisto(TString path) const {
1563
1564   //Write the histogram in the default file HFMassFitterOutput.root
1565
1566   if (fcounter == 0) {
1567     cout<<"Use MassFitter method before WriteHisto"<<endl;
1568     return;
1569   }
1570   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1571
1572   path += "HFMassFitterOutput.root";
1573   TFile *output;
1574
1575   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1576   else output = new TFile(path.Data(),"update");
1577   output->cd();
1578   hget->Write();
1579   fContourGraph->Write();
1580
1581
1582   output->Close();
1583
1584   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1585
1586   delete output;
1587
1588 }
1589
1590 //_________________________________________________________________________
1591
1592 void AliHFMassFitter::WriteNtuple(TString path) const{
1593   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1594   path += "HFMassFitterOutput.root";
1595   TFile *output = new TFile(path.Data(),"update");
1596   output->cd();
1597   fntuParam->Write();
1598   //nget->Write();
1599   output->Close();
1600   //cout<<nget->GetName()<<" written in "<<path<<endl;
1601   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1602   /*
1603   if(nget) {
1604     //delete nget;
1605     nget=NULL;
1606   }
1607   */
1608
1609   delete output;
1610 }
1611
1612 //_________________________________________________________________________
1613 void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1614
1615  //write the canvas in a root file
1616
1617   gStyle->SetOptStat(0);
1618   gStyle->SetCanvasColor(0);
1619   gStyle->SetFrameFillColor(0);
1620
1621   TString type="";
1622
1623   switch (ftypeOfFit4Bkg){
1624   case 0:
1625     type="Exp"; //3+2
1626     break;
1627   case 1:
1628     type="Lin"; //3+2
1629     break;
1630   case 2:
1631     type="Pl2"; //3+3
1632     break;
1633   case 3:
1634     type="noB"; //3+1
1635     break;
1636   case 4:
1637     type="Pow"; //3+3
1638     break;
1639   case 5:
1640     type="PowExp"; //3+3
1641     break;
1642   }
1643
1644   TString filename=Form("%sMassFit.root",type.Data());
1645   filename.Prepend(userIDstring);
1646   path.Append(filename);
1647
1648   TFile* outputcv=new TFile(path.Data(),"update");
1649
1650   TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1651   c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1652   if(draw)c->DrawClone();
1653   outputcv->cd();
1654   c->Write();
1655   outputcv->Close();
1656 }
1657
1658 //_________________________________________________________________________
1659
1660 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1661   //return a TVirtualPad with the fitted histograms and info
1662
1663   TString cvtitle="fit of ";
1664   cvtitle+=fhistoInvMass->GetName();
1665   TString cvname="c";
1666   cvname+=fcounter;
1667
1668   TCanvas *c=new TCanvas(cvname,cvtitle);
1669   PlotFit(c->cd(),nsigma,writeFitInfo);
1670   return c->cd();
1671 }
1672 //_________________________________________________________________________
1673
1674 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1675   //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1676   gStyle->SetOptStat(0);
1677   gStyle->SetCanvasColor(0);
1678   gStyle->SetFrameFillColor(0);
1679
1680   cout<<"nsigma = "<<nsigma<<endl;
1681   cout<<"Verbosity = "<<writeFitInfo<<endl;
1682
1683   TH1F* hdraw=GetHistoClone();
1684
1685   if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1686     cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1687     return;
1688   }
1689
1690   if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1691     cout<<"Drawing background fit only"<<endl;
1692     hdraw->SetMinimum(0);
1693     hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1694     pd->cd();
1695     hdraw->SetMarkerStyle(20);
1696     hdraw->DrawClone("PE");
1697     hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1698
1699     if(writeFitInfo > 0){
1700       TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1701       pinfo->SetBorderSize(0);
1702       pinfo->SetFillStyle(0);
1703       TF1* f=hdraw->GetFunction("funcbkgonly");
1704       for (Int_t i=0;i<fNFinalPars-3;i++){
1705         pinfo->SetTextColor(kBlue+3);
1706         TString str=Form("%s = %f #pm %f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1707         pinfo->AddText(str);
1708       }
1709
1710       pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1711       pd->cd();
1712       pinfo->DrawClone();
1713
1714
1715     }
1716
1717     return;
1718   }
1719
1720   hdraw->SetMinimum(0);
1721   hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1722   pd->cd();
1723   hdraw->SetMarkerStyle(20);
1724   hdraw->DrawClone("PE");
1725 //   if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1726 //   if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1727   if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1728
1729   if(writeFitInfo > 0){
1730     TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1731     TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1732     pinfob->SetBorderSize(0);
1733     pinfob->SetFillStyle(0);
1734     pinfom->SetBorderSize(0);
1735     pinfom->SetFillStyle(0);
1736     TF1* ff=fhistoInvMass->GetFunction("funcmass");
1737
1738     for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1739       pinfom->SetTextColor(kBlue);
1740       TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1741       if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1742     }
1743     pd->cd();
1744     pinfom->DrawClone();
1745
1746     TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1747     pinfo2->SetBorderSize(0);
1748     pinfo2->SetFillStyle(0);
1749
1750     Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1751
1752     Significance(nsigma,signif,errsignif);
1753     Signal(nsigma,signal,errsignal);
1754     Background(nsigma,bkg, errbkg);
1755     /*
1756       Significance(1.828,1.892,signif,errsignif);
1757       Signal(1.828,1.892,signal,errsignal);
1758       Background(1.828,1.892,bkg, errbkg);
1759     */
1760     TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1761     pinfo2->AddText(str);
1762     str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1763     pinfo2->AddText(str);
1764     str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1765     pinfo2->AddText(str);
1766
1767     pd->cd();
1768     pinfo2->Draw();
1769
1770     if(writeFitInfo > 1){
1771       for (Int_t i=0;i<fNFinalPars-3;i++){
1772         pinfob->SetTextColor(kRed);
1773         str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1774         pinfob->AddText(str);
1775       }
1776       pd->cd();
1777       pinfob->DrawClone();
1778     }
1779
1780
1781   }
1782   return;
1783 }
1784
1785 //_________________________________________________________________________
1786
1787 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1788   //draws histogram together with fit functions with default nice colors in user canvas
1789   PlotFit(pd,nsigma,writeFitInfo);
1790
1791   pd->Draw();
1792
1793 }
1794 //_________________________________________________________________________
1795 void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1796
1797   //draws histogram together with fit functions with default nice colors
1798   gStyle->SetOptStat(0);
1799   gStyle->SetCanvasColor(0);
1800   gStyle->SetFrameFillColor(0);
1801
1802
1803   TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1804   c->Draw();
1805
1806 }
1807
1808
1809 //_________________________________________________________________________
1810
1811 void AliHFMassFitter::PrintParTitles() const{
1812
1813   //prints to screen the parameters names
1814
1815   TF1 *f=fhistoInvMass->GetFunction("funcmass");
1816   if(!f) {
1817     cout<<"Fit function not found"<<endl;
1818     return;
1819   }
1820
1821   cout<<"Parameter Titles \n";
1822   for(Int_t i=0;i<fNFinalPars;i++){
1823     cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1824   }
1825   cout<<endl;
1826
1827 }
1828
1829 //************ significance
1830
1831 //_________________________________________________________________________
1832
1833 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1834   // Return signal integral in mean+- n sigma
1835
1836   if(fcounter==0) {
1837     cout<<"Use MassFitter method before Signal"<<endl;
1838     return;
1839   }
1840
1841   Double_t min=fMass-nOfSigma*fSigmaSgn;
1842   Double_t max=fMass+nOfSigma*fSigmaSgn;
1843
1844   Signal(min,max,signal,errsignal);
1845
1846
1847   return;
1848
1849 }
1850
1851 //_________________________________________________________________________
1852
1853 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1854
1855   // Return signal integral in a range
1856
1857   if(fcounter==0) {
1858     cout<<"Use MassFitter method before Signal"<<endl;
1859     return;
1860   }
1861
1862   //functions names
1863   TString bkgname="funcbkgRecalc";
1864   TString bkg1name="funcbkg1Recalc";
1865   TString massname="funcmass";
1866
1867
1868   TF1 *funcbkg=0;
1869   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1870   if(!funcmass){
1871     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1872     return;
1873   }
1874
1875   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1876   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1877
1878   if(!funcbkg){
1879     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1880     return;
1881   }
1882
1883   Int_t np=fNFinalPars-3;
1884
1885   Double_t intS,intSerr;
1886
1887  //relative error evaluation
1888   intS=funcmass->GetParameter(np);
1889   intSerr=funcmass->GetParError(np);
1890
1891   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1892   Double_t background,errbackground;
1893   Background(min,max,background,errbackground);
1894
1895   //signal +/- error in the range
1896
1897   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1898   signal=mass - background;
1899   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1900
1901 }
1902
1903 //_________________________________________________________________________
1904
1905 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1906   // Return background integral in mean+- n sigma
1907
1908   if(fcounter==0) {
1909     cout<<"Use MassFitter method before Background"<<endl;
1910     return;
1911   }
1912   Double_t min=fMass-nOfSigma*fSigmaSgn;
1913   Double_t max=fMass+nOfSigma*fSigmaSgn;
1914
1915   Background(min,max,background,errbackground);
1916
1917   return;
1918
1919 }
1920 //___________________________________________________________________________
1921
1922 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1923   // Return background integral in a range
1924
1925   if(fcounter==0) {
1926     cout<<"Use MassFitter method before Background"<<endl;
1927     return;
1928   }
1929
1930   //functions names
1931   TString bkgname="funcbkgRecalc";
1932   TString bkg1name="funcbkg1Recalc";
1933
1934   TF1 *funcbkg=0;
1935   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1936   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1937   if(!funcbkg){
1938     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1939     return;
1940   }
1941
1942
1943   Double_t intB,intBerr;
1944
1945   //relative error evaluation: from final parameters of the fit
1946   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1947   else{
1948     intB=funcbkg->GetParameter(0);
1949     intBerr=funcbkg->GetParError(0);
1950     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1951   }
1952
1953   //relative error evaluation: from histo
1954
1955   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1956   Double_t sum2=0;
1957   for(Int_t i=1;i<=fSideBandl;i++){
1958     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1959   }
1960   for(Int_t i=fSideBandr;i<=fNbin;i++){
1961     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1962   }
1963
1964   intBerr=TMath::Sqrt(sum2);
1965   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1966
1967   cout<<"Last estimation of bkg error is used"<<endl;
1968
1969   //backround +/- error in the range
1970   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1971     background = 0;
1972     errbackground = 0;
1973   }
1974   else{
1975     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1976     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1977     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1978   }
1979   return;
1980
1981 }
1982
1983
1984 //__________________________________________________________________________
1985
1986 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1987   // Return significance in mean+- n sigma
1988
1989   Double_t min=fMass-nOfSigma*fSigmaSgn;
1990   Double_t max=fMass+nOfSigma*fSigmaSgn;
1991   Significance(min, max, significance, errsignificance);
1992
1993   return;
1994 }
1995
1996 //__________________________________________________________________________
1997
1998 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1999   // Return significance integral in a range
2000
2001   Double_t signal,errsignal,background,errbackground;
2002   Signal(min, max,signal,errsignal);
2003   Background(min, max,background,errbackground);
2004
2005   if (signal+background <= 0.){
2006     cout<<"Cannot calculate significance because of div by 0!"<<endl;
2007     significance=-1;
2008     errsignificance=0;
2009     return;
2010   }
2011
2012   significance =  signal/TMath::Sqrt(signal+background);
2013
2014   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
2015
2016   return;
2017 }
2018
2019