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