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