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