AliPhysics  41af4b0 (41af4b0)
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
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),
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;
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();
143 
145 
146  fContourGraph=new TList();
147  fContourGraph->SetOwner();
148 
149 }
150 
151 
153  TNamed(),
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),
166  ffactor(mfit.ffactor),
167  fntuParam(mfit.fntuParam),
168  fMass(mfit.fMass),
169  fMassErr(mfit.fMassErr),
170  fSigmaSgn(mfit.fSigmaSgn),
172  fRawYield(mfit.fRawYield),
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),
182 {
183  //copy constructor
184 
185  if(mfit.fParsSize > 0){
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 
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 
215 
216  //assignment operator
217 
218  if(&mfit == this) return *this;
220  fminMass= mfit.fminMass;
221  fmaxMass= mfit.fmaxMass;
222  fNbin= mfit.fNbin;
223  fParsSize= mfit.fParsSize;
224  fWithBkg= mfit.fWithBkg;
227  ffactor= mfit.ffactor;
228  fntuParam= mfit.fntuParam;
229  fMass= mfit.fMass;
230  fMassErr= mfit.fMassErr;
231  fSigmaSgn= mfit.fSigmaSgn;
233  fRawYield= mfit.fRawYield;
235  fSideBands= mfit.fSideBands;
236  fSideBandl= mfit.fSideBandl;
237  fSideBandr= mfit.fSideBandr;
238  fFitOption= mfit.fFitOption;
239  fcounter= mfit.fcounter;
240  fNpfits = mfit.fNpfits;
242 
243  if(mfit.fParsSize > 0){
244  delete[] fFitPars;
245 
247  memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
248 
249  delete[] fFixPar;
250 
252  memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Float_t));
253  }
254 
255  return *this;
256 }
257 
258 //************ tools & settings
259 
260 //__________________________________________________________________________
261 
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 
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 //___________________________________________________________________________
332 
333  //Set default values for fFixPar (only total integral fixed)
334 
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 //___________________________________________________________________________
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...";
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 //___________________________________________________________________________
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 
405 }
406 
407 //___________________________________________________________________________
408 
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 
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 
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 
548  // Create, fill and return ntuple with fit parameters
549 
550  InitNtuParam(ntuname.Data());
551  FillNtuParam();
552  return fntuParam;
553 }
554 //_________________________________________________________________________
555 
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 
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 //_________________________________________________________________________
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 
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 //__________________________________________________________________________
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 
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 //__________________________________________________________________________
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
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
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 
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(fminBinMass,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fmaxBinMass,"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 
1368  if (draw) DrawFit();
1369 
1370 
1371  return kTRUE;
1372 }
1373 
1374 //______________________________________________________________________________
1375 
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  }
1433 
1434  if(draw) DrawFit();
1435 
1436  return kTRUE;
1437 
1438 }
1439 //_________________________________________________________________________
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 //_________________________________________________________________________
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 //_________________________________________________________________________
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 //_________________________________________________________________________
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 //_________________________________________________________________________
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 //_________________________________________________________________________
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 //_________________________________________________________________________
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 
1659 
1660  TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1661  return hout;
1662 }
1663 //_________________________________________________________________________
1664 
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 
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  TCanvas* c = static_cast<TCanvas*>(GetPad(nsigma,writeFitInfo));
1752  c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1753  if (draw) c->DrawClone();
1754 
1755  TFile outputcv(path.Data(),"update");
1756  outputcv.cd();
1757  c->Write();
1758  outputcv.Close();
1759 }
1760 
1761 //_________________________________________________________________________
1762 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo) const
1763 {
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, nsigma, writeFitInfo);
1773  return c;
1774 }
1775 //_________________________________________________________________________
1776 
1777 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const
1778 {
1779  //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1780  gStyle->SetOptStat(0);
1781  gStyle->SetCanvasColor(0);
1782  gStyle->SetFrameFillColor(0);
1783 
1784  cout<<"nsigma = "<<nsigma<<endl;
1785  cout<<"Verbosity = "<<writeFitInfo<<endl;
1786 
1787  TH1F* hdraw=GetHistoClone();
1788 
1789  if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1790  cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1791  return;
1792  }
1793 
1794  if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1795  cout<<"Drawing background fit only"<<endl;
1796  hdraw->SetMinimum(0);
1797  hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1798  pd->cd();
1799  hdraw->SetMarkerStyle(20);
1800  hdraw->DrawClone("PE");
1801  hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1802 
1803  if(writeFitInfo > 0){
1804  TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1805  pinfo->SetBorderSize(0);
1806  pinfo->SetFillStyle(0);
1807  TF1* f=hdraw->GetFunction("funcbkgonly");
1808  for (Int_t i=0;i<fNFinalPars-3;i++){
1809  pinfo->SetTextColor(kBlue+3);
1810  TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1811  pinfo->AddText(str);
1812  }
1813 
1814  pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1815  pd->cd();
1816  pinfo->Draw();
1817  }
1818 
1819  return;
1820  }
1821 
1822  hdraw->SetMinimum(0);
1823  hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1824  pd->cd();
1825  hdraw->SetMarkerStyle(20);
1826  hdraw->DrawClone("PE");
1827  // if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1828  // if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1829  if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1830 
1831  if(writeFitInfo > 0){
1832  TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1833  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1834  pinfob->SetBorderSize(0);
1835  pinfob->SetFillStyle(0);
1836  pinfom->SetBorderSize(0);
1837  pinfom->SetFillStyle(0);
1838  TF1* ff=fhistoInvMass->GetFunction("funcmass");
1839 
1840  for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1841  pinfom->SetTextColor(kBlue);
1842  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1843  if (ff->GetParameter(fNFinalPars-2)<0.2) str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i)*1000,ff->GetParError(i)*1000);
1844  if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1845  }
1846  pd->cd();
1847  pinfom->Draw();
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->Draw();
1883  }
1884  }
1885  return;
1886 }
1887 
1888 //_________________________________________________________________________
1889 
1890 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1891  //draws histogram together with fit functions with default nice colors in user canvas
1892  PlotFit(pd,nsigma,writeFitInfo);
1893 
1894  pd->Draw();
1895 
1896 }
1897 //_________________________________________________________________________
1899 
1900  //draws histogram together with fit functions with default nice colors
1901  gStyle->SetOptStat(0);
1902  gStyle->SetCanvasColor(0);
1903  gStyle->SetFrameFillColor(0);
1904 
1905 
1906  TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1907  c->Draw();
1908 
1909 }
1910 
1911 
1912 //_________________________________________________________________________
1913 
1915 
1916  //prints to screen the parameters names
1917 
1918  TF1 *f=fhistoInvMass->GetFunction("funcmass");
1919  if(!f) {
1920  cout<<"Fit function not found"<<endl;
1921  return;
1922  }
1923 
1924  cout<<"Parameter Titles \n";
1925  for(Int_t i=0;i<fNFinalPars;i++){
1926  cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1927  }
1928  cout<<endl;
1929 
1930 }
1931 
1932 //************ significance
1933 
1934 //_________________________________________________________________________
1935 
1936 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1937  // Return signal integral in mean+- n sigma
1938 
1939  if(fcounter==0) {
1940  cout<<"Use MassFitter method before Signal"<<endl;
1941  return;
1942  }
1943 
1944  Double_t min=fMass-nOfSigma*fSigmaSgn;
1945  Double_t max=fMass+nOfSigma*fSigmaSgn;
1946 
1947  Signal(min,max,signal,errsignal);
1948 
1949 
1950  return;
1951 
1952 }
1953 
1954 //_________________________________________________________________________
1955 
1956 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1957 
1958  // Return signal integral in a range
1959 
1960  if(fcounter==0) {
1961  cout<<"Use MassFitter method before Signal"<<endl;
1962  return;
1963  }
1964 
1965  //functions names
1966  TString bkgname="funcbkgRecalc";
1967  TString bkg1name="funcbkg1Recalc";
1968  TString massname="funcmass";
1969 
1970 
1971  TF1 *funcbkg=0;
1972  TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1973  if(!funcmass){
1974  cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1975  return;
1976  }
1977 
1978  if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1979  else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1980 
1981  if(!funcbkg){
1982  cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1983  return;
1984  }
1985 
1986  Int_t np=fNFinalPars-3;
1987 
1988  Double_t intS,intSerr;
1989 
1990  //relative error evaluation
1991  intS=funcmass->GetParameter(np);
1992  intSerr=funcmass->GetParError(np);
1993 
1994  cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1995  Double_t background,errbackground;
1996  Background(min,max,background,errbackground);
1997 
1998  //signal +/- error in the range
1999 
2000  Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
2001  signal=mass - background;
2002  errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
2003 
2004 }
2005 
2006 //_________________________________________________________________________
2007 
2008 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
2009  // Return background integral in mean+- n sigma
2010 
2011  if(fcounter==0) {
2012  cout<<"Use MassFitter method before Background"<<endl;
2013  return;
2014  }
2015  Double_t min=fMass-nOfSigma*fSigmaSgn;
2016  Double_t max=fMass+nOfSigma*fSigmaSgn;
2017 
2018  Background(min,max,background,errbackground);
2019 
2020  return;
2021 
2022 }
2023 //___________________________________________________________________________
2024 
2025 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
2026  // Return background integral in a range
2027 
2028  if(fcounter==0) {
2029  cout<<"Use MassFitter method before Background"<<endl;
2030  return;
2031  }
2032 
2033  //functions names
2034  TString bkgname="funcbkgRecalc";
2035  TString bkg1name="funcbkg1Recalc";
2036 
2037  TF1 *funcbkg=0;
2038  if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
2039  else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
2040  if(!funcbkg){
2041  cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
2042  return;
2043  }
2044 
2045 
2046  Double_t intB,intBerr;
2047 
2048  //relative error evaluation: from final parameters of the fit
2049  if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
2050  else{
2051  intB=funcbkg->GetParameter(0);
2052  intBerr=funcbkg->GetParError(0);
2053  cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
2054  }
2055 
2056  //relative error evaluation: from histo
2057 
2058  intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
2059  Double_t sum2=0;
2060  for(Int_t i=1;i<=fSideBandl;i++){
2061  sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
2062  }
2063  for(Int_t i=fSideBandr;i<=fNbin;i++){
2064  sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
2065  }
2066 
2067  intBerr=TMath::Sqrt(sum2);
2068  cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
2069 
2070  cout<<"Last estimation of bkg error is used"<<endl;
2071 
2072  //backround +/- error in the range
2073  if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
2074  background = 0;
2075  errbackground = 0;
2076  }
2077  else{
2078  background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
2079  errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
2080  //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
2081  }
2082  return;
2083 
2084 }
2085 
2086 
2087 //__________________________________________________________________________
2088 
2089 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
2090  // Return significance in mean+- n sigma
2091 
2092  Double_t min=fMass-nOfSigma*fSigmaSgn;
2093  Double_t max=fMass+nOfSigma*fSigmaSgn;
2094  Significance(min, max, significance, errsignificance);
2095 
2096  return;
2097 }
2098 
2099 //__________________________________________________________________________
2100 
2101 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
2102  // Return significance integral in a range
2103 
2104  if(fcounter==0){
2105  AliError("Number of fits is zero, check whether you made the fit before computing the significance!\n");
2106  return;
2107  }
2108 
2109  Double_t signal,errsignal,background,errbackground;
2110  Signal(min, max,signal,errsignal);
2111  Background(min, max,background,errbackground);
2112 
2113  if (signal+background <= 0.){
2114  cout<<"Cannot calculate significance because of div by 0!"<<endl;
2115  significance=-1;
2116  errsignificance=0;
2117  return;
2118  }
2119 
2120  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,background,errbackground,significance,errsignificance);
2121 
2122  return;
2123 }
2124 
2125 
2126 TH1F* AliHFMassFitter::GetResidualsAndPulls(TH1 *h,TF1 *f,Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
2127  Int_t binmi=1,binma=h->GetNbinsX();
2128 
2129  if(maxrange>minrange){
2130  binmi=h->FindBin(minrange*1.001);
2131  binma=h->FindBin(maxrange*0.9999);
2132  }
2133  if(hResidualTrend){
2134  //h->Copy(hResidualTrend);
2135  hResidualTrend->SetBins(h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax());
2136  hResidualTrend->SetName(Form("%s_residualTrend",h->GetName()));
2137  hResidualTrend->SetTitle(Form("%s (Residuals)",h->GetTitle()));
2138  hResidualTrend->SetMarkerStyle(20);
2139  hResidualTrend->SetMarkerSize(1.0);
2140  hResidualTrend->Reset();
2141  }
2142  if(hPullsTrend){
2143  hPullsTrend->SetBins(h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax());
2144  hPullsTrend->Reset();
2145  hPullsTrend->SetName(Form("%s_pullTrend",h->GetName()));
2146  hPullsTrend->SetTitle(Form("%s (Pulls)",h->GetTitle()));
2147  hPullsTrend->SetMarkerStyle(20);
2148  hPullsTrend->SetMarkerSize(1.0);
2149  }
2150  if(hPulls){
2151  hPulls->SetName(Form("%s_pulls",h->GetName()));
2152  hPulls->SetTitle(Form("%s ; Pulls",h->GetTitle()));
2153  hPulls->SetBins(40,-10,10);
2154  hPulls->Reset();
2155  }
2156 
2157  Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
2158  TArrayD *arval=new TArrayD(binma-binmi+1);
2159  for(Int_t jst=1;jst<=h->GetNbinsX();jst++){
2160 
2161  res=h->GetBinContent(jst)-f->Integral(h->GetBinLowEdge(jst),h->GetBinLowEdge(jst)+h->GetBinWidth(jst))/h->GetBinWidth(jst);
2162  if(jst>=binmi&&jst<=binma){
2163  arval->AddAt(res,jst-binmi);
2164  if(res<min)min=res;
2165  if(res>max)max=res;
2166  }
2167  // Printf("Res = %f from %f - %f",res,h->GetBinContent(jst),f->Integral(h->GetBinLowEdge(jst),h->GetBinLowEdge(jst)+h->GetBinWidth(jst))/h->GetBinWidth(jst));
2168  if(hResidualTrend){
2169  hResidualTrend->SetBinContent(jst,res);
2170  hResidualTrend->SetBinError(jst,h->GetBinError(jst));
2171  }
2172  if(hPulls){
2173  if(jst>=binmi&&jst<=binma)hPulls->Fill(res/h->GetBinError(jst));
2174  }
2175  if(hPullsTrend){
2176  hPullsTrend->SetBinContent(jst,res/h->GetBinError(jst));
2177  hPullsTrend->SetBinError(jst,0.0001);
2178  }
2179  }
2180  if(hResidualTrend)hResidualTrend->GetXaxis()->SetRange(binmi,binma);
2181  if(hPullsTrend){
2182  hPullsTrend->GetXaxis()->SetRange(binmi,binma);
2183  hPullsTrend->SetMinimum(-7);
2184  hPullsTrend->SetMaximum(+7);
2185  }
2186  if(TMath::Abs(min)>TMath::Abs(max))max=min;
2187 
2188  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);
2189  for(Int_t j=0;j<binma-binmi+1;j++){
2190  hout->Fill(arval->At(j));
2191  }
2192  hout->Sumw2();
2193  hout->Fit("gaus","LEM","",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
2194 
2195  if(hPulls){
2196  hPulls->Sumw2();
2197  hPulls->Fit("gaus","LEM","",-3,3);
2198  }
2199  delete arval;
2200  return hout;
2201 }
2202 
2203 
2204 TH1F* AliHFMassFitter::GetOverBackgroundResidualsAndPulls(Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
2205 
2206 
2207  if(!fhistoInvMass){
2208  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, invariant mass histogram not avaialble!");
2209  return 0x0;
2210  }
2211 
2212  TF1 *fback=fhistoInvMass->GetFunction("funcbkgRecalc");
2213  if(!fback){
2214  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, funcbkgRecalc not available!");
2215  return 0x0;
2216  }
2217 
2218  // THIS LONG WAY TO CP THE FUNC IS NEEDED ONLY TO EXTEND THE RANGE OF THE FUNCTION: NOT POSSIBLE OTHERWISE (WHY??? REALLY UNCOMFORTABLE)
2219 
2220  // TF1 *fbackCp=(TF1*)fback->Clone("ftmpback");
2221  TF1 *fbackCp=new TF1("ftmpback",this,&AliHFMassFitter::FitFunction4Bkg,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1), fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");;
2222  for(Int_t i=0;i< fNFinalPars-3;i++){
2223  fbackCp->SetParameter(i,fback->GetParameter(i));
2224  }
2225  // fbackCp->SetName("ftmpback");
2226  // fbackCp->SetRange(minrange*0.5,maxrange*1.5);
2227  // fback->Copy(*fbackCp);
2228  // fbackCp->SetRange(minrange*0.5,maxrange*1.5);
2229  // fbackCp->Paint();
2230  // fbackCp->Update();
2231 
2232  TH1F *h=GetResidualsAndPulls(fhistoInvMass,fbackCp,minrange,maxrange,hPulls,hResidualTrend,hPullsTrend);
2233  delete fbackCp;
2234 
2235  if(fSigmaSgn<0){
2236  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, negative sigma: fit not performed or something went wrong, cannto draw gaussian term on top of residuals");
2237  return h;
2238  }
2239 
2240  if(hResidualTrend){
2241  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));
2242  fgauss->SetParameter(0,fRawYield*fhistoInvMass->GetBinWidth(1));
2243  fgauss->SetParameter(1,fMass);
2244  fgauss->SetParameter(2,fSigmaSgn);
2245  fgauss->SetLineColor(kBlue);
2246  hResidualTrend->GetListOfFunctions()->Add(fgauss);
2247  }
2248  return h;
2249 }
2250 
2251 
2252 TH1F* AliHFMassFitter::GetAllRangeResidualsAndPulls(Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
2253 
2254 
2255  if(!fhistoInvMass){
2256  Printf("AliHFMassFitter::GetAllRangeResidualsAndPulls, invariant mass histogram not avaialble!");
2257  return 0x0;
2258  }
2259 
2260  TF1 *f=fhistoInvMass->GetFunction("funcmass");
2261  if(!f){
2262  Printf("AliHFMassFitter::GetAllRangeResidualsAndPulls, funcmass not available!");
2263  return 0x0;
2264  }
2265 
2266  // THIS LONG WAY TO CP THE FUNC IS NEEDED ONLY TO EXTEND THE RANGE OF THE FUNCTION: NOT POSSIBLE OTHERWISE (WHY??? REALLY UNCOMFORTABLE)
2267  TF1 *fmassCp=new TF1("fmassCp",this,&AliHFMassFitter::FitFunction4MassDistr,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1), fNFinalPars,"AliHFMassFitter","FitFunction4MassDistr");
2268  for(Int_t i=0;i< fNFinalPars;i++){
2269  fmassCp->SetParameter(i,f->GetParameter(i));
2270  }
2271 
2272  TH1F *h=GetResidualsAndPulls(fhistoInvMass,fmassCp,minrange,maxrange,hPulls,hResidualTrend,hPullsTrend);
2273  delete fmassCp;
2274  return h;
2275 
2276 }
2277 
2279 
2280  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
2281  pinfom->SetBorderSize(0);
2282  pinfom->SetFillStyle(0);
2283  TF1* ff=fhistoInvMass->GetFunction("funcmass");
2284 
2285  for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
2286  pinfom->SetTextColor(kBlue);
2287  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
2288  if(!(mode==1 && i==fNFinalPars-3)) pinfom->AddText(str);
2289  }
2290 
2291 
2292  if(mode>1){
2293  pinfom->SetTextColor(kBlue+5);
2294  for (Int_t i=0;i<fNFinalPars-3;i++){
2295  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
2296  pinfom->AddText(str);
2297  }
2298  }
2299  pinfom->AddText(Form("#chi^{2}/NDF=%.2f/%d",ff->GetChisquare(),ff->GetNDF()));
2300  return pinfom;
2301 }
2302 
2303 
2305  TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
2306  pinfo2->SetBorderSize(0);
2307  pinfo2->SetFillStyle(0);
2308 
2309  Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
2310 
2311  Signal(nsigma,signal,errsignal);
2312  Background(nsigma,bkg, errbkg);
2313  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,bkg,errbkg,signif,errsignif);
2314 
2315 
2316  TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
2317  pinfo2->AddText(str);
2318  str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
2319  pinfo2->AddText(str);
2320  str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
2321  pinfo2->AddText(str);
2322  if(bkg>0) str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
2323  pinfo2->AddText(str);
2324 
2325  return pinfo2;
2326 
2327 }
virtual TH1F * GetOverBackgroundResidualsAndPulls(Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
TVirtualPad * GetPad(Double_t nsigma=3, Int_t writeFitInfo=1) const
TH1F * GetHistoClone() const
Float_t * fFitPars
number of parameters of the final function
const char * filename
Definition: TestFCM.C:1
void DrawHere(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1) const
write the canvas in a root file
Bool_t * fFixPar
kTRUE = only side bands considered
virtual Bool_t CheckRangeFit()
virtual TPaveText * GetFitParametersBox(Double_t nsigma=3., Int_t mode=0)
void SetHisto(const TH1F *histoToFit)
setters
void WriteNtuple(TString path="./") const
write the histogram
Int_t fSideBandr
left side band limit (bin number)
double Double_t
Definition: External.C:58
void GetSideBandsBounds(Int_t &lb, Int_t &hb) const
const char * title
Definition: MakeQAPdf.C:27
virtual void PlotFit(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1) const
void DrawFit(Double_t nsigma=3) const
Double_t mass
Int_t fParsSize
number of bins
virtual void IntS(Float_t *valuewitherror) const
Int_t ftypeOfFit4Sgn
0 = exponential; 1 = linear; 2 = pol2
virtual void ComputeParSize()
Double_t fmaxMass
lower mass limit
Int_t ftypeOfFit4Bkg
signal+background (kTRUE) or signal only (kFALSE)
TCanvas * c
Definition: TestFitELoss.C:172
Double_t GetFitProbability() const
virtual ~AliHFMassFitter()
Double_t fMassErr
signal gaussian mean value
virtual Bool_t RefitWithBkgOnly(Bool_t draw=kTRUE)
virtual void AddFunctionsToHisto()
void PrintParTitles() const
Int_t fmaxBinMass
bin corresponding to fminMass
virtual TPaveText * GetYieldBox(Double_t nsigma=3.)
Double_t fminMass
histogram to fit
Double_t * sigma
virtual void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
return total integral of the histogram
Int_t fNbin
bin corresponding to fmaxMass
int Int_t
Definition: External.C:63
virtual Double_t FitFunction4Bkg(Double_t *x, Double_t *par)
Double_t fRawYieldErr
signal gaussian integral
Double_t fSigmaSgnErr
signal gaussian sigma
float Float_t
Definition: External.C:68
Double_t GetChiSquare() const
Double_t GetMean() const
Double_t GetReducedChiSquare() const
Double_t fMass
contains fit parameters
Bool_t GetFixThisParam(Int_t thispar) const
Int_t fcounter
right side band limit (bin number)
void FillNtuParam()
initialize TNtuple to store the parameters
void SetType(Int_t fittypeb, Int_t fittypes)
Int_t fNFinalPars
size of fFitPars array
void InitNtuParam(TString ntuname="ntupar")
void GetFitPars(Float_t *pars) const
virtual void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
signal in (min, max) with error
Double_t nsigma
Int_t mode
Definition: anaM.C:41
virtual Double_t FitFunction4MassDistr(Double_t *x, Double_t *par)
significance in (min, max) with error
Double_t fRawYield
err signal gaussian sigma
TH1F * GetResidualsAndPulls(TH1 *h, TF1 *f, Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
void RebinMass(Int_t bingroup=1)
static void ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance, Double_t &errsignificance)
Significance calculator.
virtual Bool_t MassFitter(Bool_t draw=kTRUE)
Bool_t draw[nPtBins]
TNtuple * fntuParam
number to multiply to the sigma of the signal to obtain the reflected gaussian
virtual void ComputeNFinalPars()
virtual void SetDefaultFixParam()
Double_t fSigmaSgn
err signal gaussian mean value
TList * fContourGraph
L, LW or Chi2.
Double_t GetBkgChiSquare()
Double_t GetBkgReducedChiSquare()
virtual TH1F * GetAllRangeResidualsAndPulls(Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Int_t rebin
Int_t ffactor
0 = gaus; 1 = gaus+gaus broadened
void WriteHisto(TString path="./") const
the three functions above all together
Int_t fminBinMass
upper mass limit
const Int_t nbins
Double_t GetSigma() const
bool Bool_t
Definition: External.C:53
virtual Bool_t SetFixThisParam(Int_t thispar, Bool_t fixpar)
Bool_t fSideBands
err on signal gaussian integral
Int_t fNpfits
internal counter
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
backgournd in (min, max) with error
TString fFitOption
Number of points used in the fit.
virtual Double_t FitFunction4Sgn(Double_t *x, Double_t *par)
AliHFMassFitter for the fit of invariant mass distribution of charmed mesons.
Definition: External.C:196
virtual void WriteCanvas(TString userIDstring="", TString path="./", Double_t nsigma=3, Int_t writeFitInfo=1, Bool_t draw=kFALSE) const
write the TNtuple
AliHFMassFitter & operator=(const AliHFMassFitter &mfit)
TNtuple * NtuParamOneShot(TString ntuname="ntupar")
return the TNtuple