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