AliPhysics  608b256 (608b256)
AliHFMassFitterVAR.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 // Class deriving from AliHFMassFitter for the fit of D0 invariant mass distribution
21 // including the possibility of using templates for reflections
22 //
23 // Main changes:
24 // - more flexibility in defining the number of back, signal, and reflection parameters (--> background function).
25 // - add possibility to fit with reflection template (more options available)
26 //
27 // Probably a shorter and more efficient version could have been obtained writing a totally independent class.
28 // Author: A. Rossi, andrea.rossi@cern.ch
30 
31 #include <Riostream.h>
32 #include <TMath.h>
33 #include <TNtuple.h>
34 #include <TH1F.h>
35 #include <TF1.h>
36 #include <TList.h>
37 #include <TFile.h>
38 #include <TCanvas.h>
39 #include <TVirtualPad.h>
40 #include <TGraph.h>
41 #include <TVirtualFitter.h>
42 #include <TMinuit.h>
43 #include <TStyle.h>
44 #include <TPaveText.h>
45 #include <TDatabasePDG.h>
46 #include "AliHFMassFitter.h"
47 #include "AliHFMassFitterVAR.h"
48 #include "AliVertexingHFUtils.h"
49 
50 
51 using std::cout;
52 using std::endl;
53 
55 ClassImp(AliHFMassFitterVAR);
57 
58 //************** constructors
61  fNparSignal(3),
62  fNparBack(2),
63  fNparRefl(0),
64  fhTemplRefl(0x0),
65  fSignParNames(0x0),
66  fBackParNames(0x0),
67  fReflParNames(0x0),
68  fSmoothRefl(kFALSE),
69  fReflInit(0),
70  fFixParSign(0x0),
71  fFixParBack(0x0),
72  fFixParRefl(0x0),
73  fFixParSignExternalValue(0x0),
74  fFixParBackExternalValue(0x0),
75  fFixParReflExternalValue(0x0),
76  fparSignFixExt(0x0),
77  fparBackFixExt(0x0),
78  fparReflFixExt(0x0),
79  fRawYieldHelp(0),
80  fpolbackdegreeTay(4),
81  fpolbackdegreeTayHelp(-1),
82  fMassParticle(1.864)
83 {
84  // default constructor
85 
86  cout<<"Default constructor:"<<endl;
87  cout<<"Remember to set the Histo, the Type, the FixPar"<<endl;
88 
89 }
90 
91 //___________________________________________________________________________
92 AliHFMassFitterVAR::AliHFMassFitterVAR (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
94  fNparSignal(3),
95  fNparBack(2),
96  fNparRefl(0),
97  fhTemplRefl(0x0),
98  fSignParNames(0x0),
99  fBackParNames(0x0),
100  fReflParNames(0x0),
101  fSmoothRefl(kFALSE),
102  fReflInit(0),
103  fFixParSign(0x0),
104  fFixParBack(0x0),
105  fFixParRefl(0x0),
109  fparSignFixExt(0x0),
110  fparBackFixExt(0x0),
111  fparReflFixExt(0x0),
112  fRawYieldHelp(0),
115  fMassParticle(1.864)
116 {
117  // standard constructor
118 
119  fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
120  fhistoInvMass->SetDirectory(0);
121  fminMass=minvalue;
122  fmaxMass=maxvalue;
123  if(rebin!=1) RebinMass(rebin);
124  else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
125  CheckRangeFit();
126  ftypeOfFit4Bkg=fittypeb;
127  ftypeOfFit4Sgn=fittypes;
129  else fWithBkg=kTRUE;
130  if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
131  else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
132 
133  ComputeParSize();
136 
138 
139  fContourGraph=new TList();
140  fContourGraph->SetOwner();
141 
142 }
143 
144 
146  AliHFMassFitter(mfit),
147  fNparSignal(mfit.fNparSignal),
148  fNparBack(mfit.fNparBack),
149  fNparRefl(mfit.fNparRefl),
150  fhTemplRefl(mfit.fhTemplRefl),
151  fSignParNames(0x0),
152  fBackParNames(0x0),
153  fReflParNames(0x0),
154  fSmoothRefl(mfit.fSmoothRefl),
155  fReflInit(mfit.fReflInit),
156  fFixParSign(0x0),
157  fFixParBack(0x0),
158  fFixParRefl(0x0),
162  fparSignFixExt(0x0),
163  fparBackFixExt(0x0),
164  fparReflFixExt(0x0),
169 {
170  //copy constructor
172  for(Int_t j=0;j<fNparSignal;j++){
173  fSignParNames[j]=mfit.fSignParNames[j];
174  }
175 
177  for(Int_t j=0;j<fNparBack;j++){
178  fBackParNames[j]=mfit.fBackParNames[j];
179  }
180 
182  for(Int_t j=0;j<fNparRefl;j++){
183  fReflParNames[j]=mfit.fReflParNames[j];
184  }
185 
189 
193 
197 
198  memcpy(fFixParSign,mfit.fFixParSign,mfit.fNparSignal*sizeof(Bool_t));
199  memcpy(fFixParBack,mfit.fFixParBack,mfit.fNparBack*sizeof(Bool_t));
200  memcpy(fFixParRefl,mfit.fFixParRefl,mfit.fNparRefl*sizeof(Double_t));
201 
205 
206  memcpy(fparSignFixExt,mfit.fparSignFixExt,mfit.fNparSignal*sizeof(Bool_t));
207  memcpy(fparBackFixExt,mfit.fparBackFixExt,mfit.fNparBack*sizeof(Bool_t));
208  memcpy(fparReflFixExt,mfit.fparReflFixExt,mfit.fNparRefl*sizeof(Double_t));
209 }
210 
211 //_________________________________________________________________________
213 {
214  //destructor
215 
216  if(fhTemplRefl) delete fhTemplRefl;
217  delete [] fSignParNames;
218  delete [] fBackParNames;
219  delete [] fReflParNames;
220 
221  delete [] fFixParSign;
222  delete [] fFixParBack;
223  delete [] fFixParRefl;
224 
225 
226  delete [] fFixParSignExternalValue;
227  delete [] fFixParBackExternalValue;
228  delete [] fFixParReflExternalValue;
229 
230  delete [] fparSignFixExt;
231  delete [] fparBackFixExt;
232  delete [] fparReflFixExt;
233 }
234 
235 //_________________________________________________________________________
236 
238 
239  //assignment operator
240 
241  if(&mfit == this) return *this;
243 
245  fNparBack=mfit.fNparBack;
246  fNparRefl=mfit.fNparRefl;
248 
253 
255 
256  delete [] fSignParNames;
257  delete [] fBackParNames;
258  delete [] fReflParNames;
259 
263 
264  delete [] fparSignFixExt;
265  delete [] fparBackFixExt;
266  delete [] fparReflFixExt;
267 
271 
272 
273  delete [] fFixParSign;
274  delete [] fFixParBack;
275  delete [] fFixParRefl;
276 
280 
281 
282  delete [] fFixParSignExternalValue;
283  delete [] fFixParBackExternalValue;
284  delete [] fFixParReflExternalValue;
285 
289 
290 
292  for(Int_t j=0;j<fNparSignal;j++){
293  fSignParNames[j]=mfit.fSignParNames[j];
294  }
295 
297  for(Int_t j=0;j<fNparBack;j++){
298  fBackParNames[j]=mfit.fBackParNames[j];
299  }
300 
302  for(Int_t j=0;j<fNparRefl;j++){
303  fReflParNames[j]=mfit.fReflParNames[j];
304  }
305 
306 
307 
308 
309  memcpy(fFixParSign,mfit.fFixParSign,mfit.fNparSignal*sizeof(Bool_t));
310  memcpy(fFixParBack,mfit.fFixParBack,mfit.fNparBack*sizeof(Bool_t));
311  memcpy(fFixParRefl,mfit.fFixParRefl,mfit.fNparRefl*sizeof(Bool_t));
312 
316 
317 
318  memcpy(fparSignFixExt,mfit.fparSignFixExt,mfit.fNparSignal*sizeof(Double_t));
319  memcpy(fparBackFixExt,mfit.fparBackFixExt,mfit.fNparBack*sizeof(Double_t));
320  memcpy(fparReflFixExt,mfit.fparReflFixExt,mfit.fNparRefl*sizeof(Double_t));
321 
322 
323  return *this;
324 }
325 
327  Printf("AliHFMassFitterVAR::SetBackHighPolDegree, remeber that this method has to be called after the constructor");
328  fpolbackdegreeTay=deg;
329  if(fFitPars)delete fFitPars;
330  ComputeParSize();
334 }
335 
336 
338  fhTemplRefl=(TH1F*)h->Clone("hTemplRefl");
339 
340  if(opt.EqualTo("templ")||opt.EqualTo("Templ")||opt.EqualTo("TEMPL")||opt.EqualTo("template")||opt.EqualTo("Template")||opt.EqualTo("TEMPLATE")){
341  return fhTemplRefl;
342  }
343  TF1 *f;
344  Bool_t isPoissErr=kTRUE;
345  if(opt.EqualTo("1gaus")||opt.EqualTo("singlegaus")){
346  if(minRange>=0&&maxRange>=0){
347  f=new TF1("mygaus","gaus",TMath::Max(minRange,h->GetBinLowEdge(1)),TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX())));
348  }
349  else f=new TF1("mygaus","gaus",h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
350  f->SetParameter(0,h->GetMaximum());
351  // f->SetParLimits(0,0,100.*h->Integral());
352  f->SetParameter(1,1.865);
353  f->SetParameter(2,0.050);
354 
355  fhTemplRefl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
356  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
357  fhTemplRefl->SetBinContent(j,f->Integral(fhTemplRefl->GetBinLowEdge(j),fhTemplRefl->GetXaxis()->GetBinUpEdge(j))/fhTemplRefl->GetBinWidth(j));
358  if(fhTemplRefl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
359  }
360  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
361  if(isPoissErr){
362  fhTemplRefl->SetBinError(j,TMath::Sqrt(fhTemplRefl->GetBinContent(j)));
363  }
364  else fhTemplRefl->SetBinError(j,0.001*fhTemplRefl->GetBinContent(j));
365  }
366  return fhTemplRefl;
367  }
368 
369 
370  if(opt.EqualTo("2gaus")||opt.EqualTo("doublegaus")){
371  if(minRange>=0&&maxRange>=0){
372  f=new TF1("my2gaus","[0]*([3]/( TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[3])/( TMath::Sqrt(2.*TMath::Pi())*[5])*TMath::Exp(-(x-[4])*(x-[4])/(2.*[5]*[5])))",TMath::Max(minRange,h->GetBinLowEdge(1)),TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX())));
373  }
374  else f=new TF1("my2gaus","[0]*([3]/( TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[3])/( TMath::Sqrt(2.*TMath::Pi())*[5])*TMath::Exp(-(x-[4])*(x-[4])/(2.*[5]*[5])))",h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
375 
376  f->SetParameter(0,h->GetMaximum());
377  // f->SetParLimits(0,0,100.*h->Integral());
378  f->SetParLimits(3,0.,1.);
379  f->SetParameter(3,0.5);
380 
381  f->SetParameter(1,1.84);
382  f->SetParameter(2,0.050);
383 
384  f->SetParameter(4,1.88);
385  f->SetParameter(5,0.050);
386 
387  fhTemplRefl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
388  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
389  fhTemplRefl->SetBinContent(j,f->Integral(fhTemplRefl->GetBinLowEdge(j),fhTemplRefl->GetXaxis()->GetBinUpEdge(j))/fhTemplRefl->GetBinWidth(j));
390  if(fhTemplRefl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
391  }
392  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
393  if(isPoissErr){
394  fhTemplRefl->SetBinError(j,TMath::Sqrt(fhTemplRefl->GetBinContent(j)));
395  }
396  else fhTemplRefl->SetBinError(j,0.001*fhTemplRefl->GetBinContent(j));
397  }
398  return fhTemplRefl;
399  }
400 
401 
402  if(opt.EqualTo("pol3")||opt.EqualTo("POL3")){
403  if(minRange>=0&&maxRange>=0){
404  f=new TF1("mypol3","pol3",TMath::Max(minRange,h->GetBinLowEdge(1)),TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX())));
405  }
406  else f=new TF1("mypol3","pol3",h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
407  f->SetParameter(0,h->GetMaximum());
408 
409  // f->SetParLimits(0,0,100.*h->Integral());
410  // Hard to initialize the other parameters...
411  for(Int_t nf=0;nf<10;nf++){
412  fhTemplRefl->Fit(f,"REM","");
413  //,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
414  }
415  // Printf("We USED %d POINTS in the Fit",f->GetNumberFitPoints());
416  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
417  fhTemplRefl->SetBinContent(j,f->Integral(fhTemplRefl->GetBinLowEdge(j),fhTemplRefl->GetXaxis()->GetBinUpEdge(j))/fhTemplRefl->GetBinWidth(j));
418  if(fhTemplRefl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
419  }
420  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
421  if(isPoissErr){
422  fhTemplRefl->SetBinError(j,TMath::Sqrt(fhTemplRefl->GetBinContent(j)));
423  }
424  else fhTemplRefl->SetBinError(j,0.001*fhTemplRefl->GetBinContent(j));
425  }
426  return fhTemplRefl;
427  }
428 
429 
430  if(opt.EqualTo("pol6")||opt.EqualTo("POL6")){
431 
432  if(minRange>=0&&maxRange>=0){
433  f=new TF1("mypol6","pol6",TMath::Max(minRange,h->GetBinLowEdge(1)),TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX())));
434  }
435  else f=new TF1("mypol6","pol6",h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
436  f->SetParameter(0,h->GetMaximum());
437  // f->SetParLimits(0,0,100.*h->Integral());
438  // Hard to initialize the other parameters...
439 
440  fhTemplRefl->Fit(f,"RLEMI","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
441 
442 
443  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
444  fhTemplRefl->SetBinContent(j,f->Integral(fhTemplRefl->GetBinLowEdge(j),fhTemplRefl->GetXaxis()->GetBinUpEdge(j))/fhTemplRefl->GetBinWidth(j));
445  if(fhTemplRefl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
446  }
447  for(Int_t j=1;j<=fhTemplRefl->GetNbinsX();j++){
448  if(isPoissErr){
449  fhTemplRefl->SetBinError(j,TMath::Sqrt(fhTemplRefl->GetBinContent(j)));
450  }
451  else fhTemplRefl->SetBinError(j,0.001*fhTemplRefl->GetBinContent(j));
452  }
453  return fhTemplRefl;
454  }
455 
456  // no good option passed
457  Printf("AlliHFMassFitterVAR::SetTemplateReflections: bad option");
458  return 0x0;
459 
460 
461 }
462 
463 
464 //___________________________________________________________________________
465 
467  // Fit function for signal+background
468 
469  //exponential or linear fit
470  //
471  // par[0] = tot integral
472  // par[1] = slope
473  // par[2] = gaussian integral
474  // par[3] = gaussian mean
475  // par[4] = gaussian sigma
476 
477  Double_t bkg=0,sgn=0,refl=0.;
478  Double_t parbkg[fNparBack];
479  sgn = FitFunction4Sgn(x,&par[fNparBack]);
480  if(ftypeOfFit4Sgn==2){
481  if(ftypeOfFit4Bkg==6){
482  parbkg[0]=par[0];
483  }
484  else parbkg[0]=par[0]-par[fNparBack]-par[fNparBack+fNparSignal]*par[fNparBack];
485  }
486  else {
487  if(ftypeOfFit4Bkg==6){
488  parbkg[0]=par[0];
489  }
490  else parbkg[0]=par[0]-par[fNparBack];
491  }
492  for(Int_t jb=1;jb<fNparBack;jb++){
493  parbkg[jb]=par[jb];
494  }
495  bkg = FitFunction4Bkg(x,parbkg);
496 
497 
498  if(ftypeOfFit4Sgn==2)refl=FitFunction4Refl(x,&par[fNparBack+fNparSignal]);
499 
500  return sgn+bkg+refl;
501 
502 
503 }
504 
505 //_________________________________________________________________________
507  // Fit function for the signal
508 
509  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
510  //Par:
511  // * [0] = integralSgn
512  // * [1] = mean
513  // * [2] = sigma
514  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
515  fRawYieldHelp=par[0]/fhistoInvMass->GetBinWidth(1);
516  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]);
517 
518 }
521  if(ftypeOfFit4Sgn!=2){
522  Printf("FitFunction4BkgAndRefl called but w/o reasons: this is just for drawing and requires parameters to be set from outside!!!");
523  }
524  // The following lines useful for debugging:
525  // Double_t parbkg[fNparBack];
526  // parbkg[0]=par[0]-fRawYieldHelp-par[fNparBack]*fRawYieldHelp;
527  // for(Int_t jb=1;jb<fNparBack;jb++){
528  // parbkg[jb]=par[jb];
529  // }
530  // Printf("Raw Yield: %f Func For Refl at x=%f: %f",fRawYieldHelp,x[0],FitFunction4Refl(x,&par[fNparBack]));
531  return FitFunction4Bkg(x,par)+FitFunction4Refl(x,&par[fNparBack]);
532 
533 }
536 
537  Int_t bin =fhTemplRefl->FindBin(x[0]);
538  Double_t value=fhTemplRefl->GetBinContent(bin);
539  Int_t binmin=fhTemplRefl->FindBin(fminMass*1.00001);
540  Int_t binmax=fhTemplRefl->FindBin(fmaxMass*0.99999);
541  Double_t norm=fhTemplRefl->Integral(binmin,binmax)*fhTemplRefl->GetBinWidth(bin);
542  if(TMath::Abs(value)<1.e-14&&fSmoothRefl){// very rough, assume a constant trend, much better would be a pol1 or pol2 over a broader range
543  value+=fhTemplRefl->GetBinContent(bin-1)+fhTemplRefl->GetBinContent(bin+1);
544  value/=3.;
545  }
546  return par[0]*value/norm*fRawYieldHelp*fhistoInvMass->GetBinWidth(1);
547 
548 }
549 
552  Double_t back=par[0];
553  for(Int_t it=1;it<=fpolbackdegreeTayHelp;it++){
554  back+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
555  }
556  return back;
557 
558 
559 }
560 
561 //__________________________________________________________________________
562 
564  // Fit function for the background
565 
566  Double_t maxDeltaM = 4.*fSigmaSgn;
567  if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
568  TF1::RejectPoint();
569  return 0;
570  }
571 
572  Int_t firstPar=0;
573  Double_t gaus2=0.,total=-1.;
574  if(ftypeOfFit4Sgn == 1){
575  firstPar=3;
576  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
577  //Par:
578  // * [0] = integralSgn
579  // * [1] = mean
580  // * [2] = sigma
581  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
582  gaus2 = FitFunction4Sgn(x,par);
583  }
584 
585  switch (ftypeOfFit4Bkg){
586  case 0:
587  //exponential
588  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
589  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
590  //as integralTot- integralGaus (=par [2])
591  //Par:
592  // * [0] = integralBkg;
593  // * [1] = B;
594  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
595  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]);
596  break;
597  case 1:
598  //linear
599  //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)
600  // * [0] = integralBkg;
601  // * [1] = b;
602  total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
603  break;
604  case 2:
605  //polynomial
606  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
607  //+ 1/3*c*(max^3-min^3) ->
608  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
609  // * [0] = integralBkg;
610  // * [1] = b;
611  // * [2] = c;
612  total = par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass))+par[2+firstPar]*(x[0]*x[0]-1/3.*(fmaxMass*fmaxMass*fmaxMass-fminMass*fminMass*fminMass)/(fmaxMass-fminMass));
613  break;
614  case 3:
615  total=par[0+firstPar];
616  break;
617  case 4:
618  //power function
619  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
620  //
621  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
622  // * [0] = integralBkg;
623  // * [1] = b;
624  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
625  {
626  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
627 
628  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]);
629  }
630  break;
631  case 5:
632  //power function wit exponential
633  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
634  {
635  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
636 
637  total = par[1+firstPar]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2+firstPar]*(x[0]-mpi));
638  }
639  break;
640  case 6:
641  // the following comment must be removed
642  // // pol 3, following convention for pol 2
643  // //y=a+b*x+c*x**2+d*x**3 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
644  // //+ 1/3*c*(max^3-min^3) + 1/4 d * (max^4-min^4) ->
645  // //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3) - 1/4 d * (max^4-min^4) )/(max-min)
646  // // * [0] = integralBkg;
647  // // * [1] = b;
648  // // * [2] = c;
649  // // * [3] = d;
650  {
651  total=par[0+firstPar];
652  for(Int_t it=1;it<=fpolbackdegreeTay;it++){
653  total+=par[it+firstPar]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
654  }
655 
656  }
657  break;
658  // default:
659  // Types of Fit Functions for Background:
660  // * 0 = exponential;
661  // * 1 = linear;
662  // * 2 = polynomial 2nd order
663  // * 3 = no background"<<endl;
664  // * 4 = Power function
665  // * 5 = Power function with exponential
666 
667  }
668  return total+gaus2;
669 }
670 
671 //__________________________________________________________________________
673 
674  //determines the ranges of the side bands
675 
676  if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
677  Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
678  Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin+1);
679 
680  Double_t sidebandldouble,sidebandrdouble;
681  Bool_t leftok=kFALSE, rightok=kFALSE;
682 
683  if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
684  cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
685  return kFALSE;
686  }
687 
688  //histo limit = fit function limit
689  if((TMath::Abs(fminMass-minHisto) < 1e6 || TMath::Abs(fmaxMass - maxHisto) < 1e6) && (fMass-4.*fSigmaSgn-fminMass) < 1e6){
690  Double_t coeff = (fMass-fminMass)/fSigmaSgn;
691  sidebandldouble=(fMass-0.5*coeff*fSigmaSgn);
692  sidebandrdouble=(fMass+0.5*coeff*fSigmaSgn);
693  cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
694  if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
695  if (coeff<2) {
696  cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
697  ftypeOfFit4Bkg=3;
698  //set binleft and right without considering SetRangeFit- anyway no bkg!
699  sidebandldouble=(fMass-4.*fSigmaSgn);
700  sidebandrdouble=(fMass+4.*fSigmaSgn);
701  }
702  }
703  else {
704  sidebandldouble=(fMass-4.*fSigmaSgn);
705  sidebandrdouble=(fMass+4.*fSigmaSgn);
706  }
707 
708  cout<<"Left side band ";
709  Double_t tmp=0.;
710  tmp=sidebandldouble;
711  //calculate bin corresponding to fSideBandl
712  fSideBandl=fhistoInvMass->FindBin(sidebandldouble);
713  if (sidebandldouble >= fhistoInvMass->GetBinCenter(fSideBandl)) fSideBandl++;
714  sidebandldouble=fhistoInvMass->GetBinLowEdge(fSideBandl);
715 
716  if(TMath::Abs(tmp-sidebandldouble) > 1e-6){
717  cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
718  leftok=kTRUE;
719  }
720  cout<<sidebandldouble<<" (bin "<<fSideBandl<<")"<<endl;
721 
722  cout<<"Right side band ";
723  tmp=sidebandrdouble;
724  //calculate bin corresponding to fSideBandr
725  fSideBandr=fhistoInvMass->FindBin(sidebandrdouble);
726  if (sidebandrdouble < fhistoInvMass->GetBinCenter(fSideBandr)) fSideBandr--;
727  sidebandrdouble=fhistoInvMass->GetBinLowEdge(fSideBandr+1);
728 
729  if(TMath::Abs(tmp-sidebandrdouble) > 1e-6){
730  cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
731  rightok=kTRUE;
732  }
733  cout<<sidebandrdouble<<" (bin "<<fSideBandr<<")"<<endl;
734  if (fSideBandl==0 || fSideBandr==fNbin) {
735  cout<<"Error! Range too little";
736  return kFALSE;
737  }
738  return kTRUE;
739 }
740 
741 //__________________________________________________________________________
743  //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
744 
745  if (!fhistoInvMass) {
746  cout<<"No histogram to fit! SetHisto(TH1F* h) before! "<<endl;
747  return kFALSE;
748  }
749  Bool_t leftok=kFALSE, rightok=kFALSE;
750  Int_t nbins=fhistoInvMass->GetNbinsX();
751  Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins+1);
752 
753  //check if limits are inside histogram range
754 
755  if( fminMass-minhisto < 0. ) {
756  cout<<"Out of histogram left bound! Setting to "<<minhisto<<endl;
757  fminMass=minhisto;
758  }
759  if( fmaxMass-maxhisto > 0. ) {
760  cout<<"Out of histogram right bound! Setting to"<<maxhisto<<endl;
761  fmaxMass=maxhisto;
762  }
763 
764  Double_t tmp=0.;
765  tmp=fminMass;
766  //calculate bin corresponding to fminMass
768  if (fminMass >= fhistoInvMass->GetBinCenter(fminBinMass)) fminBinMass++;
769  fminMass=fhistoInvMass->GetBinLowEdge(fminBinMass);
770  if(TMath::Abs(tmp-fminMass) > 1e-6){
771  cout<<"Left bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fminMass<<endl;
772  leftok=kTRUE;
773  }
774 
775  tmp=fmaxMass;
776  //calculate bin corresponding to fmaxMass
778  if (fmaxMass < fhistoInvMass->GetBinCenter(fmaxBinMass)) fmaxBinMass--;
779  fmaxMass=fhistoInvMass->GetBinLowEdge(fmaxBinMass+1);
780  if(TMath::Abs(tmp-fmaxMass) > 1e-6){
781  cout<<"Right bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fmaxMass<<endl;
782  rightok=kTRUE;
783  }
784 
785  return (leftok && rightok);
786 
787 }
788 
789 
790 //__________________________________________________________________________
791 
793  // Main method of the class: performs the fit of the histogram
794 
795  //Set default fitter Minuit in order to use gMinuit in the contour plots
796  TVirtualFitter::SetDefaultFitter("Minuit");
797 
798  Bool_t isBkgOnly=kFALSE;
799  Double_t slope1=-1,slope2=1,slope3=1;
800  Double_t diffUnderBands=0;
801  TF1* func;
802  Int_t fit1status=RefitWithBkgOnly(kFALSE);
803  if(fit1status){
804  Int_t checkinnsigma=4;
805  Double_t range[2]={fMass-checkinnsigma*fSigmaSgn,fMass+checkinnsigma*fSigmaSgn};
806  func=fhistoInvMass->GetFunction("funcbkgonly");
807  Double_t intUnderFunc=func->Integral(range[0],range[1]);
808  Double_t intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(range[0]),fhistoInvMass->FindBin(range[1]),"width");
809  cout<<"Pick zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
810  Double_t diffUnderPick=(intUnderHisto-intUnderFunc);
811  intUnderFunc=func->Integral(fminMass,fminMass+checkinnsigma*fSigmaSgn);
812  intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fminMass+checkinnsigma*fSigmaSgn),"width");
813  cout<<"Band (l) zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
814  diffUnderBands=(intUnderHisto-intUnderFunc);
815  Double_t relDiff=diffUnderPick/diffUnderBands;
816  cout<<"Relative difference = "<<relDiff<<endl;
817  if(TMath::Abs(relDiff) < 0.25) isBkgOnly=kTRUE;
818  else{
819  cout<<"Relative difference = "<<relDiff<<": I suppose there is some signal, continue with total fit!"<<endl;
820  }
821  }
822  if(isBkgOnly) {
823  cout<<"INFO!! The histogram contains only background"<<endl;
824  if(draw)DrawFit();
825 
826  //increase counter of number of fits done
827  fcounter++;
828 
829  return kTRUE;
830  }
831 
832 
833  //function names
834  TString bkgname="funcbkg";
835  TString bkg1name="funcbkg1";
836  TString massname="funcmass";
837 
838  //Total integral
839  Double_t totInt = fhistoInvMass->Integral(fminBinMass,fmaxBinMass,"width");
840  //cout<<"Here tot integral is = "<<totInt<<"; integral in whole range is "<<fhistoInvMass->Integral("width")<<endl;
841  fSideBands = kTRUE;
842  Double_t width=fhistoInvMass->GetBinWidth(8);
843  //cout<<"fNbin = "<<fNbin<<endl;
844  if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
845 
846  Bool_t ok=SideBandsBounds();
847  if(!ok) return kFALSE;
848 
849  //sidebands integral - first approx (from histo)
850  Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(fminBinMass,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fmaxBinMass,"width");
851  cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
852  cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
853  if (sideBandsInt<=0) {
854  cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
855  return kFALSE;
856  }
857 
858  /*Fit Bkg*/
859  TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitterVAR::FitFunction4Bkg,fminMass,fmaxMass,fNparBack,"AliHFMassFitterVAR","FitFunction4Bkg");
860  // TF1 *funcbkg = GetBackgroundFunc();
861  cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
862 
863  funcbkg->SetLineColor(2); //red
864 
865  // cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
866  Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
867  Double_t *parBackInit=new Double_t[fNparBack];
868  if(func){
869  for(Int_t j=0;j<fNparBack;j++){
870  parBackInit[j]=func->GetParameter(j);
871  }
872  }
873  else {
874  Print("First fit background function not present");
875  delete funcbkg;
876  delete[] parBackInit;
877  return kFALSE;
878  }
879 
880  if(ftypeOfFit4Bkg!=6){// THIS IS LIKELY UNNECESSARY
881  parBackInit[0]=totInt;// this and the following 2 lines are from previous AliHFMassFitter version... can be removed
882  }
883 
884  for(Int_t j=0;j<fNparBack;j++){
885  // Printf(" HERE back %d",j);
886  funcbkg->SetParName(j,fBackParNames[j]);
887  funcbkg->SetParameter(j,parBackInit[j]);
889  funcbkg->FixParameter(j,fparBackFixExt[j]);
890  }
891  else if(fFixParBack[j]){
892  funcbkg->FixParameter(j,parBackInit[j]);
893  }
894  }
895 
896 
897 
898  Double_t intbkg1=0,conc1=0;
899  //if only signal and reflection: skip
900  if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
901  // ftypeOfFit4Sgn=0;
902  fhistoInvMass->Fit(funcbkg,"R,E,0");
903 
904  fSideBands = kFALSE;
905  //intbkg1 = funcbkg->GetParameter(0);
906 
907  intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
908  if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
909  if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
910  if(ftypeOfFit4Bkg==5) conc1 = funcbkg->GetParameter(2);
911 
912 
913  //cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
914  }
915  else cout<<"\t\t//"<<endl;
916  if(ftypeOfFit4Bkg==6){
917  // THE FOLLOWING LINES ARE NEEDED FOR THE FOLLOWING REASON:
918  // when a histogram is fitted with the function TF1 *f, the function that is added to the list of function is a clone of f,
919  // so the pointers are different
920  // Conversely, if the function is added to the list as f, of course the pointers remain the same
921  TF1 *fh=fhistoInvMass->GetFunction(bkgname.Data());
922  PrepareHighPolFit(fh);
923  }
924 
925 
926 
927 
928  //sidebands integral - second approx (from fit)
929  fSideBands = kFALSE;
930 
931  Double_t bkgInt;
932  //cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
933 
934  if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg->Integral(fminMass,fmaxMass);
935  else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
936  //cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
937 
938  //Signal integral - first approx
939  Double_t sgnInt=diffUnderBands;
940  sgnInt = totInt-bkgInt;
941  Printf("Estimates: bkgInt: %f, totInt: %f, diffUndBan:%f",bkgInt,totInt,diffUnderBands);
942  //cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
943  /*Fit All Mass distribution with exponential + gaussian (+gaussian braodened) */
944 
945  // Double_t sgnInt=diffUnderBands;
946  if(sgnInt<0)sgnInt=0.1*totInt;
947  TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitterVAR::FitFunction4MassDistr,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitterVAR","FitFunction4MassDistr");
948  cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
949  funcmass->SetLineColor(4); //blue
950 
951  //Set parameters
952  cout<<"\nTOTAL FIT"<<endl;
953 
954  Double_t parSignInit[3]={sgnInt,fMass,fSigmaSgn};
955  Double_t parReflInit[2]={fReflInit,0.};
956  fRawYieldHelp=sgnInt;
957  Printf("Initial parameters: \n back: tot int: %f, slope1:%f, slope2:%f \n sign: sgnInt %f, mean %f, sigma %f",totInt,slope1,slope2,sgnInt,fMass,fSigmaSgn);
958  for(Int_t j=0;j<fNparBack;j++){
959 
960  funcmass->SetParName(j,fBackParNames[j]);
961  funcmass->SetParameter(j,parBackInit[j]);
963  funcmass->FixParameter(j,fparBackFixExt[j]);
964  }
965  else if(fFixParBack[j]){
966  funcmass->FixParameter(j,parBackInit[j]);
967  }
968  }
969  for(Int_t j=0;j<fNparSignal;j++){
970 
971  funcmass->SetParName(j+fNparBack,fSignParNames[j]);
972  funcmass->SetParameter(j+fNparBack,parSignInit[j]);
974  funcmass->FixParameter(j+fNparBack,fparSignFixExt[j]);
975  }
976  else if(fFixParSign[j]){
977  funcmass->FixParameter(j+fNparBack,parSignInit[j]);
978  }
979 
980  }
981  for(Int_t j=0;j<fNparRefl;j++){
982 
983  funcmass->SetParName(j+fNparBack+fNparSignal,fReflParNames[j]);
984  funcmass->SetParameter(j+fNparBack+fNparSignal,parReflInit[j]);
985  funcmass->SetParLimits(j+fNparBack+fNparSignal,0.,1.);
987  funcmass->FixParameter(j+fNparBack+fNparSignal,fparReflFixExt[j]);
988  }
989  else if(fFixParRefl[j]){
990  funcmass->FixParameter(j+fNparBack+fNparSignal,parReflInit[j]);
991  }
992 
993  }
994 
995  Int_t status;
996  Printf("Fitting");
997  status = fhistoInvMass->Fit(massname.Data(),Form("R,%s,+,0",fFitOption.Data()));
998  if (status != 0){
999  cout<<"Minuit returned "<<status<<endl;
1000  delete funcbkg;
1001  delete[] parBackInit;
1002  return kFALSE;
1003  }
1004  Printf("Fitted");
1005  Printf("MassFitter: Value of func at 1.800: %f",funcmass->Eval(1.8000));
1006  cout<<"fit done"<<endl;
1007  //reset value of fMass and fSigmaSgn to those found from fit
1008  fMass=funcmass->GetParameter(fNparBack+1);
1009  fMassErr=funcmass->GetParError(fNparBack+1);
1010  fSigmaSgn=funcmass->GetParameter(fNparBack+2);
1011  fSigmaSgnErr=funcmass->GetParError(fNparBack+2);
1012  fRawYield=funcmass->GetParameter(fNparBack)/fhistoInvMass->GetBinWidth(1);
1014  fRawYieldErr=funcmass->GetParError(fNparBack)/fhistoInvMass->GetBinWidth(1);
1015 
1016  // The following lines can be useful for debugging
1017  // for(Int_t i=0;i<fNFinalPars;i++){
1018  // fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1019  // fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1020  // //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1021  // }
1022  // /*
1023  // //check: cout parameters
1024  // for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
1025  // cout<<i<<"\t"<<fFitPars[i]<<endl;
1026  // }
1027  // */
1028 
1029  if(funcmass->GetParameter(fNparBack+2) <0 || funcmass->GetParameter(fNparBack+1) <0 || funcmass->GetParameter(fNparBack) <0 ) {
1030  cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1031  delete funcbkg;
1032  delete[] parBackInit;
1033  delete funcmass;
1034  return kFALSE;
1035  }
1036 
1037  //increase counter of number of fits done
1038  fcounter++;
1039 
1040 
1041  delete[] parBackInit;
1042  delete funcbkg;
1043  delete funcmass;
1044 
1045  // Printf("Now calling add functions to hist in MassFitter");
1047  // Printf("After calling add functions to hist in MassFitter");
1048  if (draw) DrawFit();
1049  // Printf("After Draw Fit");
1050 
1051  return kTRUE;
1052 }
1053 //________________________________________________________________________
1055  // Perform intermediate fit steps up to fpolbackdegreeTay-1
1056  TH1F *hCp=(TH1F*)fhistoInvMass->Clone("htemp");
1057  Double_t estimatecent=0.5*(hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn))+hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn)));// just a first rough estimate
1058  Double_t estimateslope=(hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn))-hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn)))/(7*fSigmaSgn);// first rough estimate
1059 
1060  for(Int_t j=1;j<=hCp->GetNbinsX();j++){
1061  // h->SetBinContent(j,f->Eval(h->GetBinCenter(j))+fSignal->Integral(h->GetBinLowEdge(j),h->GetBinLowEdge(j+1)));
1062  // h->SetBinError(j,TMath::Sqrt(h->GetBinContent(j)));
1063  Double_t binCenter=hCp->GetBinCenter(j);
1064  if(binCenter>(fMass-2.5*fSigmaSgn) && binCenter<(fMass+2.5*fSigmaSgn)){//ranges are ok for D0 up to 16 GeV/c
1065  hCp->SetBinContent(j,0);
1066  hCp->SetBinError(j,0);
1067  }
1068  }
1069 
1071  TF1 *funcbkg,*funcPrev=0x0;
1073  funcbkg = new TF1(Form("temp%d",fpolbackdegreeTayHelp),this,&AliHFMassFitterVAR::BackFitFuncPolHelper,fminMass,fmaxMass,fpolbackdegreeTayHelp+1,"AliHFMassFitterVAR","BackFitFuncPolHelper");
1074  if(funcPrev){
1075  for(Int_t j=0;j<fpolbackdegreeTayHelp;j++){// now is +1 degree w.r.t. previous fit funct
1076  funcbkg->SetParameter(j,funcPrev->GetParameter(j));
1077  }
1078  delete funcPrev;
1079  }
1080  else{
1081  funcbkg->SetParameter(0,estimatecent);
1082  funcbkg->SetParameter(1,estimateslope);
1083 
1084  }
1085  hCp->Fit(funcbkg,"REMN","");
1086  funcPrev=(TF1*)funcbkg->Clone("ftemp");
1087  delete funcbkg;
1089  }
1090 
1091  // TString name=fback->GetName();
1092  // funcbkg->Copy(*fback);
1093  // fback->SetName(name.Data());
1094  for(Int_t j=0;j<=fpolbackdegreeTay;j++){
1095  fback->SetParameter(j,funcPrev->GetParameter(j));
1096  fback->SetParError(j,funcPrev->GetParError(j));
1097  }
1098  hCp->Fit(fback,"REMN","");// THIS IS JUST TO SET NOT ONLY THE PARAMETERS BUT ALSO chi2, etc...
1099 
1100 
1101  // The following lines might be useful for debugging
1102  // TCanvas *cDebug=new TCanvas();
1103  // cDebug->cd();
1104  // hCp->Draw();
1105  // TString strout=Form("Test%d.root",(Int_t)fhistoInvMass->GetBinContent(fhistoInvMass->FindBin(fMass)));
1106  // cDebug->Print(strout.Data());
1107  // delete cDebug;
1108 
1109  delete funcPrev;
1110  delete hCp;
1111  return kTRUE;
1112 
1113 }
1114 //______________________________________________________________________________
1115 
1117 
1118  //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
1119  //If you want to change the backgroud function or range use SetType or SetRangeFit before
1120 
1121  TString bkgname="funcbkgonly";
1122  fSideBands = kFALSE;
1123  Int_t typesSave=ftypeOfFit4Sgn;
1125  TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitterVAR::FitFunction4Bkg,fminMass,fmaxMass,fNparBack,"AliHFMassFitterVAR","FitFunction4Bkg");
1126 
1127  funcbkg->SetLineColor(kBlue+3); //dark blue
1128 
1129  Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
1130 
1131  switch (ftypeOfFit4Bkg) {
1132  case 0: //gaus+expo
1133  funcbkg->SetParNames("BkgInt","Slope");
1134  funcbkg->SetParameters(integral,-2.);
1135  break;
1136  case 1:
1137  funcbkg->SetParNames("BkgInt","Slope");
1138  funcbkg->SetParameters(integral,-100.);
1139  break;
1140  case 2:
1141  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1142  funcbkg->SetParameters(integral,-10.,5);
1143  break;
1144  case 3:
1145  cout<<"Warning! This choice does not make a lot of sense..."<<endl;
1146  if(ftypeOfFit4Sgn==0){
1147  funcbkg->SetParNames("Const");
1148  funcbkg->SetParameter(0,0.);
1149  funcbkg->FixParameter(0,0.);
1150  }
1151  break;
1152  case 4:
1153  funcbkg->SetParNames("BkgInt","Coef1");
1154  funcbkg->SetParameters(integral,0.5);
1155  break;
1156  case 5:
1157  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1158  funcbkg->SetParameters(integral,-10.,5.);
1159  break;
1160  case 6:
1161  for(Int_t ipb=0;ipb<=fpolbackdegreeTay;ipb++){
1162  funcbkg->SetParName(ipb,Form("Coef%d",ipb));
1163  }
1164  break;
1165  default:
1166  cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1167  ftypeOfFit4Sgn=typesSave;
1168  delete funcbkg;
1169  return kFALSE;
1170  break;
1171  }
1172 
1173  Int_t status=0;
1174  if(ftypeOfFit4Bkg==6){
1175  if(PrepareHighPolFit(funcbkg)){
1176 
1177  fhistoInvMass->GetListOfFunctions()->Add(funcbkg);
1178  fhistoInvMass->GetFunction(funcbkg->GetName())->SetBit(1<<9,kTRUE);
1179  }
1180  }
1181  else status=fhistoInvMass->Fit(bkgname.Data(),"R,E,+,0");
1182  if (status != 0){
1183  ftypeOfFit4Sgn=typesSave;
1184  cout<<"Minuit returned "<<status<<endl;
1185  return kFALSE;
1186  }
1188 
1189  if(draw) DrawFit();
1190  ftypeOfFit4Sgn=typesSave;
1191  return kTRUE;
1192 
1193 }
1194 
1195 
1196 //_________________________________________________________________________
1197 void AliHFMassFitterVAR::IntS(Float_t *valuewitherror) const {
1198 
1199  //gives the integral of signal obtained from fit parameters
1200  if(!valuewitherror) {
1201  printf("AliHFMassFitterVAR::IntS: got a null pointer\n");
1202  return;
1203  }
1204 
1205  // Int_t index=fParsSize/2 - 3;
1206  valuewitherror[0]=fRawYield;
1207  valuewitherror[1]=fRawYieldErr;
1208 }
1209 
1210 
1211 //_________________________________________________________________________
1213 
1214  //Add the background function in the complete range to the list of functions attached to the histogram
1215 
1216  //cout<<"AddFunctionsToHisto called"<<endl;
1217  TString bkgname = "funcbkg";
1218 
1219  Bool_t done1=kFALSE,done2=kFALSE;
1220  Printf(" AddFunctionsToHisto ############# LISTING ALL FUNCTIONS #################");
1221  fhistoInvMass->GetListOfFunctions()->Print();
1222 
1223  TString bkgnamesave=bkgname;
1224  TString testname=bkgname;
1225  testname += "FullRange";
1226  TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1227  if(testfunc){
1228  done1=kTRUE;
1229  testfunc=0x0;
1230  }
1231  testname="funcbkgonly";
1232  testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1233  if(testfunc){
1234  done2=kTRUE;
1235  testfunc=0x0;
1236  }
1237 
1238  if(done1 && done2){
1239  cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1240  return;
1241  }
1242 
1243  TList *hlist=fhistoInvMass->GetListOfFunctions();
1244  hlist->ls();
1245 
1246  if(!done2){
1247  TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1248  if(!bonly){
1249  cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1250  }else{
1251  bonly->SetLineColor(kBlue+3);
1252  hlist->Add((TF1*)bonly->Clone());
1253  delete bonly;
1254  }
1255 
1256  }
1257 
1258  if(!done1){
1259  TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1260  if(!b){
1261  cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1262  return;
1263  }
1264 
1265  bkgname += "FullRange";
1266  TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitterVAR::FitFunction4Bkg,fminMass,fmaxMass,fNparBack,"AliHFMassFitterVAR","FitFunction4Bkg");
1267  //cout<<bfullrange->GetName()<<endl;
1268  for(Int_t i=0;i<fNparBack;i++){
1269  bfullrange->SetParName(i,b->GetParName(i));
1270  bfullrange->SetParameter(i,b->GetParameter(i));
1271  bfullrange->SetParError(i,b->GetParError(i));
1272  }
1273  bfullrange->SetLineStyle(4);
1274  bfullrange->SetLineColor(14);
1275 
1276 
1277  bkgnamesave += "Recalc";
1278 
1279  TF1 *blastpar;
1280  if(ftypeOfFit4Sgn<2)blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitterVAR::FitFunction4Bkg,fminMass,fmaxMass,fNparBack,"AliHFMassFitterVAR","FitFunction4Bkg");
1281  else blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitterVAR::FitFunction4BkgAndReflDraw,fminMass,fmaxMass,fNparBack+fNparRefl,"AliHFMassFitterVAR","FitFunction4BkgAndReflDraw");
1282 
1283  TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1284 
1285  if (!mass){
1286  cout<<"funcmass doesn't exist "<<endl;
1287  delete bfullrange;
1288  delete blastpar;
1289  return;
1290  }
1291 
1292 
1293 
1294  //intBkg=intTot-intS
1295  if(ftypeOfFit4Bkg==6){
1296  blastpar->SetParameter(0,mass->GetParameter(0));
1297  blastpar->SetParError(0,mass->GetParError(0));
1298  }
1299  else{
1300  blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNparBack));
1301  blastpar->SetParError(0,mass->GetParError(fNparBack));
1302  }
1303  for(Int_t jb=1;jb<fNparBack;jb++){
1304  blastpar->SetParameter(jb,mass->GetParameter(jb));
1305  blastpar->SetParError(jb,mass->GetParError(jb));
1306  }
1307  if(ftypeOfFit4Sgn==2){
1308  if(ftypeOfFit4Bkg!=6){
1309  blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNparBack)-mass->GetParameter(fNparBack+fNparSignal)*mass->GetParameter(fNparBack));
1310  blastpar->SetParError(0,mass->GetParError(fNparBack));
1311  }
1312  blastpar->SetParameter(fNparBack,mass->GetParameter(fNparBack+fNparSignal));
1313  for(Int_t jr=1;jr<fNparRefl;jr++){
1314  blastpar->SetParameter(jr+fNparBack,mass->GetParameter(fNparBack+fNparSignal+jr));
1315  blastpar->SetParError(jr+fNparBack,mass->GetParError(fNparBack+fNparSignal+jr));
1316  }
1317  }
1318  else for(Int_t jr=0;jr<fNparRefl;jr++){
1319  blastpar->SetParameter(jr+fNparBack,mass->GetParameter(fNparBack+fNparSignal+jr));
1320  blastpar->SetParError(jr+fNparBack,mass->GetParError(fNparBack+fNparSignal+jr));
1321  }
1322 
1323  blastpar->SetLineStyle(1);
1324  blastpar->SetLineColor(2);
1325 
1326  hlist->Add((TF1*)bfullrange->Clone());
1327  hlist->Add((TF1*)blastpar->Clone());
1328  hlist->ls();
1329 
1330  delete bfullrange;
1331  delete blastpar;
1332 
1333  }
1334 
1335 
1336 }
1337 //_________________________________________________________________________
1338 void AliHFMassFitterVAR::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1339 
1340  //write the canvas in a root file
1341 
1342  gStyle->SetOptStat(0);
1343  gStyle->SetCanvasColor(0);
1344  gStyle->SetFrameFillColor(0);
1345 
1346  TString type="";
1347 
1348  switch (ftypeOfFit4Bkg){
1349  case 0:
1350  type="Exp"; //3+2
1351  break;
1352  case 1:
1353  type="Lin"; //3+2
1354  break;
1355  case 2:
1356  type="Pl2"; //3+3
1357  break;
1358  case 3:
1359  type="noB"; //3+1
1360  break;
1361  case 4:
1362  type="Pow"; //3+3
1363  break;
1364  case 5:
1365  type="PowExp"; //3+3
1366  break;
1367  }
1368 
1369  TString filename=Form("%sMassFit.root",type.Data());
1370  filename.Prepend(userIDstring);
1371  path.Append(filename);
1372 
1373  TCanvas* c = static_cast<TCanvas*>(GetPad(nsigma,writeFitInfo));
1374  c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1375  if (draw) c->DrawClone();
1376 
1377  TFile outputcv(path.Data(),"update");
1378  outputcv.cd();
1379  c->Write();
1380  outputcv.Close();
1381 }
1382 
1383 
1384 //_________________________________________________________________________
1385 
1386 void AliHFMassFitterVAR::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo){
1387  //draws histogram together with fit functions with default nice colors in user canvas
1388  PlotFitVAR(pd,nsigma,writeFitInfo);
1389 
1390  pd->Draw();
1391 
1392 }
1393 //_________________________________________________________________________
1394 void AliHFMassFitterVAR::PlotFitVAR(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo){
1395  //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1396  gStyle->SetOptStat(0);
1397  gStyle->SetCanvasColor(0);
1398  gStyle->SetFrameFillColor(0);
1399 
1400  cout<<"nsigma = "<<nsigma<<endl;
1401  cout<<"Verbosity = "<<writeFitInfo<<endl;
1402 
1403  TH1F* hdraw=GetHistoClone();
1404 
1405  if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1406  cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1407  return;
1408  }
1409 
1410  if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1411  cout<<"Drawing background fit only"<<endl;
1412  hdraw->SetMinimum(0);
1413  hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1414  pd->cd();
1415  hdraw->SetMarkerStyle(20);
1416  hdraw->DrawClone("PE");
1417  hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1418 
1419  if(writeFitInfo > 0){
1420  TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1421  pinfo->SetBorderSize(0);
1422  pinfo->SetFillStyle(0);
1423  TF1* f=hdraw->GetFunction("funcbkgonly");
1424  for (Int_t i=0;i<fNparBack;i++){
1425  pinfo->SetTextColor(kBlue+3);
1426  TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1427  pinfo->AddText(str);
1428  }
1429 
1430  for (Int_t i=fNparBack;i<fNparBack+fNparSignal;i++){
1431  pinfo->SetTextColor(kBlue+3);
1432  TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1433  pinfo->AddText(str);
1434  }
1435 
1436  for (Int_t i=fNparBack+fNparSignal;i<fNparBack+fNparSignal+fNparRefl;i++){
1437  pinfo->SetTextColor(kBlue+3);
1438  TString str=Form("%s = %.3f #pm %.3f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1439  pinfo->AddText(str);
1440  }
1441 
1442  pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1443  pd->cd();
1444  pinfo->DrawClone();
1445 
1446 
1447  }
1448 
1449  return;
1450  }
1451 
1452  hdraw->SetMinimum(0);
1453  hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1454  pd->cd();
1455  hdraw->SetMarkerStyle(20);
1456  hdraw->DrawClone("PE");
1457  // if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1458  // if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1459  if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1460 
1461  if(writeFitInfo > 0){
1462  TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1463  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1464  pinfob->SetBorderSize(0);
1465  pinfob->SetFillStyle(0);
1466  pinfom->SetBorderSize(0);
1467  pinfom->SetFillStyle(0);
1468  TF1* ff=fhistoInvMass->GetFunction("funcmass");
1469 
1470  // the following was a previous choice
1471  // for (Int_t i=0;i<fNparBack;i++){
1472  // pinfom->SetTextColor(kBlue);
1473  // TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1474  // if(!(writeFitInfo==1 && i==fNparBack)) pinfom->AddText(str);
1475  // }
1476 
1477  for (Int_t i=fNparBack;i<fNparBack+fNparSignal;i++){
1478  pinfom->SetTextColor(kBlue);
1479  // TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1480  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1481  if (ff->GetParameter(fNparBack+1)<0.2) str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i)*1000,ff->GetParError(i)*1000);
1482  if(!(writeFitInfo==1 && i==fNparBack)) pinfom->AddText(str);
1483  }
1484  for (Int_t i=fNparBack+fNparSignal;i<fNparBack+fNparSignal+fNparRefl;i++){
1485  pinfom->SetTextColor(kBlue+3);
1486  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1487  pinfom->AddText(str);
1488  }
1489 
1490  pd->cd();
1491  pinfom->DrawClone();
1492 
1493  TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1494  pinfo2->SetBorderSize(0);
1495  pinfo2->SetFillStyle(0);
1496 
1497  Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1498 
1499  Signal(nsigma,signal,errsignal);
1500  Background(nsigma,bkg, errbkg);
1501  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,bkg,errbkg,signif,errsignif);
1502 
1503 
1504  TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1505  pinfo2->AddText(str);
1506  str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1507  pinfo2->AddText(str);
1508  str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1509  pinfo2->AddText(str);
1510  if(bkg>0) str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
1511  pinfo2->AddText(str);
1512 
1513  pd->cd();
1514  pinfo2->Draw();
1515  if(writeFitInfo==1){
1516  TF1 *fitBkgRec=hdraw->GetFunction("funcbkgRecalc");
1517  fitBkgRec->SetLineColor(14); // does not work
1518  fitBkgRec->SetLineStyle(2); // does not work
1519 
1520  TF1 *fitCp;
1521  if(ftypeOfFit4Sgn==2){// drawing background function w/o reflection contribution
1522  fitCp= new TF1("fbackcpfordrawing",this,&AliHFMassFitterVAR::FitFunction4BkgAndReflDraw,fminMass,fmaxMass,fNparBack+fNparRefl,"AliHFMassFitterVAR","FitFunction4BkgAndReflDraw");
1523  fitCp->SetParameter(fNparBack,0);// set to 0 reflection normalization
1524 
1525  for(Int_t ibk=0;ibk<fNparBack;ibk++){
1526  fitCp->SetParameter(ibk,fitBkgRec->GetParameter(ibk));
1527  }
1528  fitCp->SetLineColor(kMagenta);
1529  fitCp->SetLineStyle(7);
1530  fitCp->SetLineWidth(2);
1531  fitCp->DrawCopy("Lsame");
1532  //Printf("WHERE I SHOULD BE: npars func=%d; par 0=%f, par1=%f,par2=%f",fitCp->GetNpar(),fitCp->GetParameter(0),fitCp->GetParameter(1),fitCp->GetParameter(2));
1533  delete fitCp;
1534  }
1535  }
1536 
1537  if(writeFitInfo > 1){
1538  for (Int_t i=0;i<fNparBack+fNparSignal;i++){
1539  pinfob->SetTextColor(kRed);
1540  str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1541  pinfob->AddText(str);
1542  }
1543  pd->cd();
1544  pinfob->DrawClone();
1545  }
1546 
1547 
1548  }
1549  return;
1550 }
1551 
1552 
1553 //_________________________________________________________________________
1554 
1555 void AliHFMassFitterVAR::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1556  // Return signal integral in mean+- n sigma
1557 
1558  if(fcounter==0) {
1559  cout<<"Use MassFitter method before Signal"<<endl;
1560  return;
1561  }
1562 
1563  Double_t min=fMass-nOfSigma*fSigmaSgn;
1564  Double_t max=fMass+nOfSigma*fSigmaSgn;
1565 
1566  Signal(min,max,signal,errsignal);
1567 
1568 
1569  return;
1570 
1571 }
1572 //______________________________________________
1574  TString massname="funcmass";
1575 
1576 
1577  TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1578  if(!funcmass){
1579  cout<<"AliHFMassFitterVAR::Signal() ERROR -> Mass distr function not found!"<<endl;
1580  return -1;
1581  }
1582  if(ftypeOfFit4Sgn != 2){
1583  Printf("The fit was done w/o reflection template");
1584  return -1;
1585  }
1586  err=funcmass->GetParError(fNparBack+fNparSignal);
1587  return funcmass->GetParameter(fNparBack+fNparSignal);
1588 
1589 }
1590 //_________________________________________________________________________
1591 
1592 void AliHFMassFitterVAR::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1593 
1594  // Return signal integral in a range
1595 
1596  if(fcounter==0) {
1597  cout<<"Use MassFitter method before Signal"<<endl;
1598  return;
1599  }
1600 
1601  //functions names
1602  TString bkgname="funcbkgRecalc";
1603  TString bkg1name="funcbkg1Recalc";
1604  TString massname="funcmass";
1605 
1606 
1607  TF1 *funcbkg=0;
1608  TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1609  if(!funcmass){
1610  cout<<"AliHFMassFitterVAR::Signal() ERROR -> Mass distr function not found!"<<endl;
1611  return;
1612  }
1613 
1614  if(ftypeOfFit4Sgn == 1) funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1615  else funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1616 
1617  if(!funcbkg){
1618  cout<<"AliHFMassFitterVAR::Signal() ERROR -> Bkg function not found!"<<endl;
1619  return;
1620  }
1621 
1622  Int_t np=fNparBack;
1623 
1624  Double_t intS,intSerr;
1625 
1626  //relative error evaluation
1627  intS=funcmass->GetParameter(np);
1628  intSerr=funcmass->GetParError(np);
1629 
1630  cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1631  Double_t background,errbackground;
1632  Background(min,max,background,errbackground);
1633 
1634  //signal +/- error in the range
1635 
1636  Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1637  signal=mass - background;
1638  errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1639 
1640 }
1641 
1642 //_________________________________________________________________________
1643 
1644 void AliHFMassFitterVAR::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1645  // Return background integral in mean+- n sigma
1646 
1647  if(fcounter==0) {
1648  cout<<"Use MassFitter method before Background"<<endl;
1649  return;
1650  }
1651  Double_t min=fMass-nOfSigma*fSigmaSgn;
1652  Double_t max=fMass+nOfSigma*fSigmaSgn;
1653 
1654  Background(min,max,background,errbackground);
1655 
1656  return;
1657 
1658 }
1659 //___________________________________________________________________________
1660 
1661 void AliHFMassFitterVAR::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1662  // Return background integral in a range
1663 
1664  if(fcounter==0) {
1665  cout<<"Use MassFitter method before Background"<<endl;
1666  return;
1667  }
1668 
1669  //functions names
1670  TString bkgname="funcbkgRecalc";
1671  TString bkg1name="funcbkg1Recalc";
1672 
1673  TF1 *funcbkg=0;
1674  if(ftypeOfFit4Sgn == 1) funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1675  else funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1676  if(!funcbkg){
1677  cout<<"AliHFMassFitterVAR::Background() ERROR -> Bkg function not found!"<<endl;
1678  return;
1679  }
1680 
1681 
1682  Double_t intB,intBerr;
1683 
1684  //relative error evaluation: from final parameters of the fit
1685  if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1686  else{
1687  intB=funcbkg->GetParameter(0);
1688  intBerr=funcbkg->GetParError(0);
1689  cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1690  }
1691 
1692  //relative error evaluation: from histo
1693 
1694  intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1695  Double_t sum2=0;
1696  for(Int_t i=1;i<=fSideBandl;i++){
1697  sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1698  }
1699  for(Int_t i=fSideBandr;i<=fNbin;i++){
1700  sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1701  }
1702 
1703  intBerr=TMath::Sqrt(sum2);
1704  cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1705 
1706  cout<<"Last estimation of bkg error is used"<<endl;
1707 
1708  //backround +/- error in the range
1709  if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1710  background = 0;
1711  errbackground = 0;
1712  }
1713  else{
1714  background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1715  errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1716  //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1717  }
1718  return;
1719 
1720 }
1721 
1722 
1723 //__________________________________________________________________________
1724 
1725 void AliHFMassFitterVAR::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
1726  // Return significance in mean+- n sigma
1727 
1728  Double_t min=fMass-nOfSigma*fSigmaSgn;
1729  Double_t max=fMass+nOfSigma*fSigmaSgn;
1730  Significance(min, max, significance, errsignificance);
1731 
1732  return;
1733 }
1734 
1735 //__________________________________________________________________________
1736 
1737 void AliHFMassFitterVAR::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1738  // Return significance integral in a range
1739 
1740  Double_t signal,errsignal,background,errbackground;
1741  Signal(min, max,signal,errsignal);
1742  Background(min, max,background,errbackground);
1743 
1744  if (signal+background <= 0.){
1745  cout<<"Cannot calculate significance because of div by 0!"<<endl;
1746  significance=-1;
1747  errsignificance=0;
1748  return;
1749  }
1750 
1751  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,background,errbackground,significance,errsignificance);
1752 
1753  return;
1754 }
1755 
1757 
1758  //compute the size of the parameter array and set the data member
1759 
1760  switch (ftypeOfFit4Bkg) {//npar background func
1761  case 0:
1762  fParsSize = 2*3;
1763  break;
1764  case 1:
1765  fParsSize = 2*3;
1766  break;
1767  case 2:
1768  fParsSize = 3*3;
1769  break;
1770  case 3:
1771  fParsSize = 1*3;
1772  break;
1773  case 4:
1774  fParsSize = 2*3;
1775  break;
1776  case 5:
1777  fParsSize = 3*3;
1778  break;
1779  case 6:
1780  fParsSize = (fpolbackdegreeTay+1)*3;
1781  break;
1782 
1783  default:
1784  cout<<"AliHFMassFitterVAR, Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
1785  break;
1786  }
1787 
1788  fParsSize += 3; // npar refl
1789  fParsSize += 3; // npar signal gaus
1790 
1791  fParsSize*=2; // add errors
1792  cout<<"Parameters array size "<<fParsSize<<endl;
1793 }
1794 
1795 
1797 
1798  //compute the number of parameters of the total (signal+bgk) function
1799  cout<<"Info:ComputeNFinalPars... ";
1800  delete [] fBackParNames;
1801  delete [] fFixParBack;
1802  delete [] fFixParBackExternalValue;
1803  delete [] fparBackFixExt;
1804 
1805  delete [] fSignParNames;
1806  delete [] fFixParSign;
1807  delete [] fFixParSignExternalValue;
1808  delete [] fparSignFixExt;
1809 
1810  delete [] fReflParNames;
1811  delete [] fFixParRefl;
1812  delete [] fFixParReflExternalValue;
1813  delete [] fparReflFixExt;
1814 
1815 
1816  switch (ftypeOfFit4Bkg) {//npar background func
1817  case 0:
1818  fNparBack=2;
1820  fBackParNames[0]="BkgInt";
1821  fBackParNames[1]="Slope";
1822  break;
1823  case 1:
1824  fNparBack=2;
1826  fBackParNames[0]="BkgInt";
1827  fBackParNames[1]="Slope";
1828  break;
1829  case 2:
1830  fNparBack=3;
1832  fBackParNames[0]="BkgInt";
1833  fBackParNames[1]="Coef1";
1834  fBackParNames[2]="Coef2";
1835  break;
1836  case 3:
1837  fNparBack=1;
1839  fBackParNames[0]="Const";
1840  break;
1841  case 4:
1842  fNparBack=2;
1844  fBackParNames[0]="BkgInt";
1845  fBackParNames[1]="Coef1";
1846  break;
1847  case 5:
1848  fNparBack=3;
1850  fBackParNames[0]="BkgInt";
1851  fBackParNames[1]="Coef1";
1852  fBackParNames[2]="Coef2";
1853  break;
1854  case 6:
1857  for(Int_t j=0;j<fNparBack;j++){
1858  fBackParNames[j]=Form("Coef%d",j);
1859  }
1860  break;
1861  default:
1862  cout<<"AliHFMassFitterVAR, Error in computing fNparBack: check ftypeOfFit4Bkg"<<endl;
1863  break;
1864  }
1865 
1869 
1870  for(Int_t ib=0;ib<fNparBack;ib++){
1871  fFixParBack[ib]=kFALSE;
1872  fFixParBackExternalValue[ib]=kFALSE;
1873  fparBackFixExt[ib]=0.;
1874  }
1875 
1876  fNparSignal=3;
1877  fSignParNames=new TString[3];
1878  fSignParNames[0]="Signal";
1879  fSignParNames[1]="Mean";
1880  fSignParNames[2]="Sigma";
1881 
1885 
1886  for(Int_t ib=0;ib<fNparSignal;ib++){
1887  fFixParSign[ib]=kFALSE;
1888  fFixParSignExternalValue[ib]=kFALSE;
1889  fparSignFixExt[ib]=0.;
1890  }
1891 
1892 
1893  fNparRefl=0;
1894  if(ftypeOfFit4Sgn==2){
1895  fNparRefl=1;
1897  fReflParNames[0]="ReflOverS";
1898  }
1899 
1900 
1901  if(fNparRefl>0){
1905  }
1906  else {
1907  fFixParRefl=0x0;
1908  fparReflFixExt=0x0;
1909  fparReflFixExt=0x0;
1910  }
1911 
1912 
1913  for(Int_t ib=0;ib<fNparRefl;ib++){
1914  fFixParRefl[ib]=kFALSE;
1915  fFixParReflExternalValue[ib]=kFALSE;
1916  fparReflFixExt[ib]=0.;
1917  }
1918 
1919  fNFinalPars=fNparBack+fNparSignal+fNparRefl;
1920 
1921  cout<<": "<<fNFinalPars<<endl;
1922 
1923  Printf("Total number of par: %d, Back:%d, Sign:%d, Refl: %d",fNFinalPars,fNparBack,fNparSignal,fNparRefl);
1924 }
1925 
1926 
1928  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1929  pinfom->SetBorderSize(0);
1930  pinfom->SetFillStyle(0);
1931  TF1* ff=fhistoInvMass->GetFunction("funcmass");
1932 
1933 
1934  for (Int_t i=fNparBack;i<fNparBack+fNparSignal;i++){
1935  pinfom->SetTextColor(kBlue);
1936  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1937  if(!(mode==1 && i==fNparBack)) pinfom->AddText(str);
1938  }
1939  for (Int_t i=fNparBack+fNparSignal;i<fNparBack+fNparSignal+fNparRefl;i++){
1940  pinfom->SetTextColor(kBlue+3);
1941  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1942  pinfom->AddText(str);
1943  }
1944 
1945  if(mode>1){
1946  pinfom->SetTextColor(kBlue+5);
1947  for (Int_t i=0;i<fNparBack;i++){
1948  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1949  pinfom->AddText(str);
1950  }
1951  }
1952  pinfom->AddText(Form("#chi^{2}/NDF=%.2f/%d",ff->GetChisquare(),ff->GetNDF()));
1953 
1954  return pinfom;
1955 }
1956 
1957 
1959  TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1960  pinfo2->SetBorderSize(0);
1961  pinfo2->SetFillStyle(0);
1962 
1963  Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1964 
1965  Signal(nsigma,signal,errsignal);
1966  Background(nsigma,bkg, errbkg);
1967  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,bkg,errbkg,signif,errsignif);
1968 
1969 
1970  TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1971  pinfo2->AddText(str);
1972  str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1973  pinfo2->AddText(str);
1974  str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1975  pinfo2->AddText(str);
1976  if(bkg>0) str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
1977  pinfo2->AddText(str);
1978 
1979  return pinfo2;
1980 
1981 }
1982 
1983 
1984 
1985 TH1F* AliHFMassFitterVAR::GetOverBackgroundResidualsAndPulls(Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
1986 
1987 
1988  if(!fhistoInvMass){
1989  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, invariant mass histogram not avaialble!");
1990  return 0x0;
1991  }
1992 
1993  TF1 *fback=fhistoInvMass->GetFunction("funcbkgRecalc");
1994  if(!fback){
1995  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, funcbkgRecalc not available!");
1996  return 0x0;
1997  }
1998 
1999  // THIS LONG WAY TO CP THE FUNC IS NEEDED ONLY TO EXTEND THE RANGE OF THE FUNCTION: NOT POSSIBLE OTHERWISE (WHY??? REALLY UNCOMFORTABLE)
2000  TF1 *fbackCp;
2001  if(ftypeOfFit4Sgn<2)fbackCp=new TF1("ftmpback",this,&AliHFMassFitterVAR::FitFunction4Bkg,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1), fNparBack,"AliHFMassFitterVAR","FitFunction4Bkg");
2002  else fbackCp=new TF1("ftmpback",this,&AliHFMassFitterVAR::FitFunction4BkgAndReflDraw,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1),fNparBack+fNparRefl,"AliHFMassFitterVAR","FitFunction4BkgAndReflDraw");
2003 
2004  for(Int_t i=0;i<fback->GetNpar();i++){
2005  fbackCp->SetParameter(i,fback->GetParameter(i));
2006  }
2007 
2008 
2009  TH1F *h=GetResidualsAndPulls(fhistoInvMass,fbackCp,minrange,maxrange,hPulls,hResidualTrend,hPullsTrend);
2010  delete fbackCp;
2011 
2012  if(fSigmaSgn<0){
2013  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, negative sigma: fit not performed or something went wrong, cannto draw gaussian term on top of residuals");
2014  return h;
2015  }
2016 
2017  if(hResidualTrend){
2018  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));
2019  fgauss->SetParameter(0,fRawYield*fhistoInvMass->GetBinWidth(1));
2020  fgauss->SetParameter(1,fMass);
2021  fgauss->SetParameter(2,fSigmaSgn);
2022  fgauss->SetLineColor(kBlue);
2023  hResidualTrend->GetListOfFunctions()->Add(fgauss);
2024  }
2025  return h;
2026 }
2027 
2028 
2029 TH1F* AliHFMassFitterVAR::GetAllRangeResidualsAndPulls(Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
2030 
2031 
2032  if(!fhistoInvMass){
2033  Printf("AliHFMassFitter::GetAllRangeResidualsAndPulls, invariant mass histogram not avaialble!");
2034  return 0x0;
2035  }
2036 
2037  TF1 *f=fhistoInvMass->GetFunction("funcmass");
2038  if(!f){
2039  Printf("AliHFMassFitter::GetAllRangeResidualsAndPulls, funcmass not available!");
2040  return 0x0;
2041  }
2042 
2043  // THIS LONG WAY TO CP THE FUNC IS NEEDED ONLY TO EXTEND THE RANGE OF THE FUNCTION: NOT POSSIBLE OTHERWISE (WHY??? REALLY UNCOMFORTABLE)
2044  TF1 *fmassCp=new TF1("fmassCp",this,&AliHFMassFitterVAR::FitFunction4MassDistr,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1), fNFinalPars,"AliHFMassFitterVAR","FitFunction4MassDistr");
2045  for(Int_t i=0;i< f->GetNpar();i++){
2046  fmassCp->SetParameter(i,f->GetParameter(i));
2047  }
2048 
2049  TH1F *h=GetResidualsAndPulls(fhistoInvMass,fmassCp,minrange,maxrange,hPulls,hResidualTrend,hPullsTrend);
2050  delete fmassCp;
2051  return h;
2052 
2053 }
TVirtualPad * GetPad(Double_t nsigma=3, Int_t writeFitInfo=1) const
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
backgournd in (min, max) with error
TH1F * GetHistoClone() const
Double_t fReflInit
smoothing refl template
Float_t * fFitPars
number of parameters of the final function
const char * filename
Definition: TestFCM.C:1
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
signal in (min, max) with error
Double_t FitFunction4BkgAndReflDraw(Double_t *x, Double_t *par)
TPaveText * GetYieldBox(Double_t nsigma=3.)
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
TPaveText * GetFitParametersBox(Double_t nsigma=3., Int_t mode=0)
TString * fSignParNames
template of reflection contribution
Double_t * fparReflFixExt
external values to fix back parameters
Int_t fSideBandr
left side band limit (bin number)
double Double_t
Definition: External.C:58
void PlotFitVAR(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1)
Int_t fNparBack
number of signal parameters
Bool_t * fFixParReflExternalValue
fix signal parameter from ext value
Bool_t RefitWithBkgOnly(Bool_t draw=kTRUE)
void WriteCanvas(TString userIDstring="", TString path="./", Double_t nsigma=3, Int_t writeFitInfo=1, Bool_t draw=kFALSE) const
setters
Bool_t * fFixParBack
fix signal parameter from ext value
void DrawFit(Double_t nsigma=3) const
void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
return total integral of the histogram
Double_t FitFunction4Sgn(Double_t *x, Double_t *par)
Double_t mass
void DrawHere(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1)
void SetBackHighPolDegree(Int_t deg)
Int_t fParsSize
number of bins
Int_t ftypeOfFit4Sgn
0 = exponential; 1 = linear; 2 = pol2
Int_t fNparRefl
number of bkg parameters
Bool_t * fFixParBackExternalValue
fix signal parameter from ext value
Double_t fmaxMass
lower mass limit
Double_t GetReflOverSignal(Double_t &err) const
significance in (min, max) with error
Int_t ftypeOfFit4Bkg
signal+background (kTRUE) or signal only (kFALSE)
TCanvas * c
Definition: TestFitELoss.C:172
Bool_t MassFitter(Bool_t draw=kTRUE)
Bool_t * fFixParRefl
fix signal parameter from ext value
TString * fReflParNames
back parameter names
Double_t fMassErr
signal gaussian mean value
Int_t fmaxBinMass
bin corresponding to fminMass
Double_t fminMass
histogram to fit
AliHFMassFitterVAR for the fit of invariant mass distribution of charmed mesons.
Int_t fNbin
bin corresponding to fmaxMass
int Int_t
Definition: External.C:63
Double_t fRawYieldErr
signal gaussian integral
Double_t fSigmaSgnErr
signal gaussian sigma
float Float_t
Definition: External.C:68
Double_t fMass
contains fit parameters
Int_t fcounter
right side band limit (bin number)
Double_t fRawYieldHelp
external values to fix refl parameters
Int_t fpolbackdegreeTay
internal variable used when fitting with reflections
Int_t fNFinalPars
size of fFitPars array
Double_t nsigma
TH1F * SetTemplateReflections(const TH1F *h, TString option="templ", Double_t minRange=1.72, Double_t maxRange=2.05)
Int_t mode
Definition: anaM.C:41
Double_t BackFitFuncPolHelper(Double_t *x, Double_t *par)
Double_t fRawYield
err signal gaussian sigma
TString * fBackParNames
signal parameter names
Double_t FitFunction4Bkg(Double_t *x, Double_t *par)
void IntS(Float_t *valuewitherror) const
TH1F * GetResidualsAndPulls(TH1 *h, TF1 *f, Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Double_t fMassParticle
help variable
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.
Bool_t draw[nPtBins]
virtual void SetDefaultFixParam()
Double_t fSigmaSgn
err signal gaussian mean value
Double_t FitFunction4Refl(Double_t *x, Double_t *par)
TList * fContourGraph
L, LW or Chi2.
AliHFMassFitterVAR & operator=(const AliHFMassFitterVAR &mfit)
Double_t * fparBackFixExt
external values to fix signal parameters
Bool_t fSmoothRefl
refl parameter names
Int_t rebin
TH1F * GetOverBackgroundResidualsAndPulls(Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Bool_t * fFixParSign
initial value of Refl/Signal
Int_t fminBinMass
upper mass limit
Int_t fpolbackdegreeTayHelp
degree of polynomial expansion for back fit (option 6 for back)
const Int_t nbins
bool Bool_t
Definition: External.C:53
Bool_t * fFixParSignExternalValue
fix signal parameter from ext value
TH1F * GetAllRangeResidualsAndPulls(Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Double_t FitFunction4MassDistr(Double_t *x, Double_t *par)
significance in (min, max) with error
Bool_t fSideBands
err on signal gaussian integral
TString fFitOption
Number of points used in the fit.
AliHFMassFitter for the fit of invariant mass distribution of charmed mesons.
Definition: External.C:196
TH1F * fhTemplRefl
number of reflection parameters
Double_t * fparSignFixExt
fix signal parameter from ext value
Bool_t PrepareHighPolFit(TF1 *fback)
AliHFMassFitter & operator=(const AliHFMassFitter &mfit)