PWG3/vertexingHF/AliHFMassFitter.cxx

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