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