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