AliPhysics  958ad07 (958ad07)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
59 AliHFMassFitterVAR::AliHFMassFitterVAR() :
60 AliHFMassFitter(),
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):
93  AliHFMassFitter(),
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),
106  fFixParSignExternalValue(0x0),
107  fFixParBackExternalValue(0x0),
108  fFixParReflExternalValue(0x0),
109  fparSignFixExt(0x0),
110  fparBackFixExt(0x0),
111  fparReflFixExt(0x0),
112  fRawYieldHelp(0),
113  fpolbackdegreeTay(4),
114  fpolbackdegreeTayHelp(-1),
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;
128  if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2 && ftypeOfFit4Bkg!=4 && ftypeOfFit4Bkg!=5&& ftypeOfFit4Bkg!=6) fWithBkg=kFALSE;
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();
134  ComputeNFinalPars();
135  fFitPars=new Float_t[fParsSize];
136 
137  SetDefaultFixParam();
138 
139  fContourGraph=new TList();
140  fContourGraph->SetOwner();
141 
142 }
143 
144 
145 AliHFMassFitterVAR::AliHFMassFitterVAR(const AliHFMassFitterVAR &mfit):
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),
159  fFixParSignExternalValue(0x0),
160  fFixParBackExternalValue(0x0),
161  fFixParReflExternalValue(0x0),
162  fparSignFixExt(0x0),
163  fparBackFixExt(0x0),
164  fparReflFixExt(0x0),
165  fRawYieldHelp(mfit.fRawYieldHelp),
166  fpolbackdegreeTay(mfit.fpolbackdegreeTay),
167  fpolbackdegreeTayHelp(mfit.fpolbackdegreeTayHelp),
168  fMassParticle(mfit.fMassParticle)
169 {
170  //copy constructor
171  fSignParNames=new TString[fNparSignal];
172  for(Int_t j=0;j<fNparSignal;j++){
173  fSignParNames[j]=mfit.fSignParNames[j];
174  }
175 
176  fBackParNames=new TString[fNparBack];
177  for(Int_t j=0;j<fNparBack;j++){
178  fBackParNames[j]=mfit.fBackParNames[j];
179  }
180 
181  fReflParNames=new TString[fNparRefl];
182  for(Int_t j=0;j<fNparRefl;j++){
183  fReflParNames[j]=mfit.fReflParNames[j];
184  }
185 
186  fFixParSign=new Bool_t[fNparSignal];
187  fFixParSignExternalValue=new Bool_t[fNparSignal];
188  fparSignFixExt=new Double_t[fNparSignal];
189 
190  fFixParBack=new Bool_t[fNparBack];
191  fFixParBackExternalValue=new Bool_t[fNparBack];
192  fparBackFixExt=new Double_t[fNparBack];
193 
194  fFixParRefl=new Bool_t[fNparRefl];
195  fFixParReflExternalValue=new Bool_t[fNparRefl];
196  fparReflFixExt=new Double_t[fNparRefl];
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 
202  memcpy(fFixParSignExternalValue,mfit.fFixParSignExternalValue,mfit.fNparSignal*sizeof(Bool_t));
203  memcpy(fFixParBackExternalValue,mfit.fFixParBackExternalValue,mfit.fNparBack*sizeof(Bool_t));
204  memcpy(fFixParReflExternalValue,mfit.fFixParReflExternalValue,mfit.fNparRefl*sizeof(Double_t));
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 //_________________________________________________________________________
212 AliHFMassFitterVAR::~AliHFMassFitterVAR()
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 
237 AliHFMassFitterVAR& AliHFMassFitterVAR::operator=(const AliHFMassFitterVAR &mfit){
238 
239  //assignment operator
240 
241  if(&mfit == this) return *this;
242  AliHFMassFitter::operator=(mfit);
243 
244  fNparSignal=mfit.fNparSignal;
245  fNparBack=mfit.fNparBack;
246  fNparRefl=mfit.fNparRefl;
247  fhTemplRefl=mfit.fhTemplRefl;
248 
249  fSmoothRefl=mfit.fSmoothRefl;
250  fRawYieldHelp=mfit.fRawYieldHelp;
251  fpolbackdegreeTay=mfit.fpolbackdegreeTay;
252  fpolbackdegreeTayHelp=mfit.fpolbackdegreeTayHelp;
253 
254  fMassParticle=mfit.fMassParticle;
255 
256  delete [] fSignParNames;
257  delete [] fBackParNames;
258  delete [] fReflParNames;
259 
260  fSignParNames=new TString[fNparSignal];
261  fBackParNames=new TString[fNparBack];
262  fReflParNames=new TString[fNparRefl];
263 
264  delete [] fparSignFixExt;
265  delete [] fparBackFixExt;
266  delete [] fparReflFixExt;
267 
268  fparSignFixExt=new Double_t[fNparSignal];
269  fparBackFixExt=new Double_t[fNparBack];
270  fparReflFixExt=new Double_t[fNparRefl];
271 
272 
273  delete [] fFixParSign;
274  delete [] fFixParBack;
275  delete [] fFixParRefl;
276 
277  fFixParSign=new Bool_t[fNparSignal];
278  fFixParBack=new Bool_t[fNparBack];
279  fFixParRefl=new Bool_t[fNparRefl];
280 
281 
282  delete [] fFixParSignExternalValue;
283  delete [] fFixParBackExternalValue;
284  delete [] fFixParReflExternalValue;
285 
286  fFixParSignExternalValue=new Bool_t[fNparSignal];
287  fFixParBackExternalValue=new Bool_t[fNparBack];
288  fFixParReflExternalValue=new Bool_t[fNparRefl];
289 
290 
291  fSignParNames=new TString[fNparSignal];
292  for(Int_t j=0;j<fNparSignal;j++){
293  fSignParNames[j]=mfit.fSignParNames[j];
294  }
295 
296  fBackParNames=new TString[fNparBack];
297  for(Int_t j=0;j<fNparBack;j++){
298  fBackParNames[j]=mfit.fBackParNames[j];
299  }
300 
301  fReflParNames=new TString[fNparRefl];
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 
313  memcpy(fFixParSignExternalValue,mfit.fFixParSignExternalValue,mfit.fNparSignal*sizeof(Bool_t));
314  memcpy(fFixParBackExternalValue,mfit.fFixParBackExternalValue,mfit.fNparBack*sizeof(Bool_t));
315  memcpy(fFixParReflExternalValue,mfit.fFixParReflExternalValue,mfit.fNparRefl*sizeof(Bool_t));
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 
326 void AliHFMassFitterVAR::SetBackHighPolDegree(Int_t deg){
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();
331  ComputeNFinalPars();
332  fFitPars=new Float_t[fParsSize];
333  SetDefaultFixParam();
334 }
335 
336 
337 TH1F* AliHFMassFitterVAR::SetTemplateReflections(const TH1F *h, TString opt,Double_t minRange,Double_t maxRange){
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 
466 Double_t AliHFMassFitterVAR::FitFunction4MassDistr (Double_t *x, Double_t *par){
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 //_________________________________________________________________________
506 Double_t AliHFMassFitterVAR::FitFunction4Sgn (Double_t *x, Double_t *par){
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 }
520 Double_t AliHFMassFitterVAR::FitFunction4BkgAndReflDraw(Double_t *x, Double_t *par){
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 }
535 Double_t AliHFMassFitterVAR::FitFunction4Refl(Double_t *x,Double_t *par){
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 
551 Double_t AliHFMassFitterVAR::BackFitFuncPolHelper(Double_t *x,Double_t *par){
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 
563 Double_t AliHFMassFitterVAR::FitFunction4Bkg(Double_t *x, Double_t *par){
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 //__________________________________________________________________________
672 Bool_t AliHFMassFitterVAR::SideBandsBounds(){
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 //__________________________________________________________________________
742 Bool_t AliHFMassFitterVAR::CheckRangeFit(){
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
767  fminBinMass=fhistoInvMass->FindBin(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
777  fmaxBinMass=fhistoInvMass->FindBin(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 
792 Bool_t AliHFMassFitterVAR::MassFitter(Bool_t draw){
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]);
888  if(fFixParBackExternalValue[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]);
962  if(fFixParBackExternalValue[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]);
973  if(fFixParSignExternalValue[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.);
986  if(fFixParReflExternalValue[j]){
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);
1013  fRawYieldHelp=fRawYield;
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");
1046  AddFunctionsToHisto();
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 //________________________________________________________________________
1054 Bool_t AliHFMassFitterVAR::PrepareHighPolFit(TF1 *fback){
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 
1070  fpolbackdegreeTayHelp=2;
1071  TF1 *funcbkg,*funcPrev=0x0;
1072  while(fpolbackdegreeTayHelp<=fpolbackdegreeTay){
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;
1088  fpolbackdegreeTayHelp++;
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 
1116 Bool_t AliHFMassFitterVAR::RefitWithBkgOnly(Bool_t draw){
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;
1124  if(ftypeOfFit4Sgn==2)ftypeOfFit4Sgn=0;
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  }
1187  AddFunctionsToHisto();
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 //_________________________________________________________________________
1212 void AliHFMassFitterVAR::AddFunctionsToHisto(){
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  if(!(writeFitInfo==1 && i==fNparBack)) pinfom->AddText(str);
1481  }
1482  for (Int_t i=fNparBack+fNparSignal;i<fNparBack+fNparSignal+fNparRefl;i++){
1483  pinfom->SetTextColor(kBlue+3);
1484  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1485  pinfom->AddText(str);
1486  }
1487 
1488  pd->cd();
1489  pinfom->DrawClone();
1490 
1491  TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1492  pinfo2->SetBorderSize(0);
1493  pinfo2->SetFillStyle(0);
1494 
1495  Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1496 
1497  Signal(nsigma,signal,errsignal);
1498  Background(nsigma,bkg, errbkg);
1499  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,bkg,errbkg,signif,errsignif);
1500 
1501 
1502  TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1503  pinfo2->AddText(str);
1504  str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1505  pinfo2->AddText(str);
1506  str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1507  pinfo2->AddText(str);
1508  if(bkg>0) str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
1509  pinfo2->AddText(str);
1510 
1511  pd->cd();
1512  pinfo2->Draw();
1513  if(writeFitInfo==1){
1514  TF1 *fitBkgRec=hdraw->GetFunction("funcbkgRecalc");
1515  fitBkgRec->SetLineColor(14); // does not work
1516  fitBkgRec->SetLineStyle(2); // does not work
1517 
1518  TF1 *fitCp;
1519  if(ftypeOfFit4Sgn==2){// drawing background function w/o reflection contribution
1520  fitCp= new TF1("fbackcpfordrawing",this,&AliHFMassFitterVAR::FitFunction4BkgAndReflDraw,fminMass,fmaxMass,fNparBack+fNparRefl,"AliHFMassFitterVAR","FitFunction4BkgAndReflDraw");
1521  fitCp->SetParameter(fNparBack,0);// set to 0 reflection normalization
1522 
1523  for(Int_t ibk=0;ibk<fNparBack;ibk++){
1524  fitCp->SetParameter(ibk,fitBkgRec->GetParameter(ibk));
1525  }
1526  fitCp->SetLineColor(kMagenta);
1527  fitCp->SetLineStyle(7);
1528  fitCp->SetLineWidth(2);
1529  fitCp->DrawCopy("Lsame");
1530  //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));
1531  delete fitCp;
1532  }
1533  }
1534 
1535  if(writeFitInfo > 1){
1536  for (Int_t i=0;i<fNparBack+fNparSignal;i++){
1537  pinfob->SetTextColor(kRed);
1538  str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1539  pinfob->AddText(str);
1540  }
1541  pd->cd();
1542  pinfob->DrawClone();
1543  }
1544 
1545 
1546  }
1547  return;
1548 }
1549 
1550 
1551 //_________________________________________________________________________
1552 
1553 void AliHFMassFitterVAR::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1554  // Return signal integral in mean+- n sigma
1555 
1556  if(fcounter==0) {
1557  cout<<"Use MassFitter method before Signal"<<endl;
1558  return;
1559  }
1560 
1561  Double_t min=fMass-nOfSigma*fSigmaSgn;
1562  Double_t max=fMass+nOfSigma*fSigmaSgn;
1563 
1564  Signal(min,max,signal,errsignal);
1565 
1566 
1567  return;
1568 
1569 }
1570 //______________________________________________
1571 Double_t AliHFMassFitterVAR::GetReflOverSignal(Double_t &err)const{
1572  TString massname="funcmass";
1573 
1574 
1575  TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1576  if(!funcmass){
1577  cout<<"AliHFMassFitterVAR::Signal() ERROR -> Mass distr function not found!"<<endl;
1578  return -1;
1579  }
1580  if(ftypeOfFit4Sgn != 2){
1581  Printf("The fit was done w/o reflection template");
1582  return -1;
1583  }
1584  err=funcmass->GetParError(fNparBack+fNparSignal);
1585  return funcmass->GetParameter(fNparBack+fNparSignal);
1586 
1587 }
1588 //_________________________________________________________________________
1589 
1590 void AliHFMassFitterVAR::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1591 
1592  // Return signal integral in a range
1593 
1594  if(fcounter==0) {
1595  cout<<"Use MassFitter method before Signal"<<endl;
1596  return;
1597  }
1598 
1599  //functions names
1600  TString bkgname="funcbkgRecalc";
1601  TString bkg1name="funcbkg1Recalc";
1602  TString massname="funcmass";
1603 
1604 
1605  TF1 *funcbkg=0;
1606  TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1607  if(!funcmass){
1608  cout<<"AliHFMassFitterVAR::Signal() ERROR -> Mass distr function not found!"<<endl;
1609  return;
1610  }
1611 
1612  if(ftypeOfFit4Sgn == 1) funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1613  else funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1614 
1615  if(!funcbkg){
1616  cout<<"AliHFMassFitterVAR::Signal() ERROR -> Bkg function not found!"<<endl;
1617  return;
1618  }
1619 
1620  Int_t np=fNparBack;
1621 
1622  Double_t intS,intSerr;
1623 
1624  //relative error evaluation
1625  intS=funcmass->GetParameter(np);
1626  intSerr=funcmass->GetParError(np);
1627 
1628  cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1629  Double_t background,errbackground;
1630  Background(min,max,background,errbackground);
1631 
1632  //signal +/- error in the range
1633 
1634  Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1635  signal=mass - background;
1636  errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1637 
1638 }
1639 
1640 //_________________________________________________________________________
1641 
1642 void AliHFMassFitterVAR::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1643  // Return background integral in mean+- n sigma
1644 
1645  if(fcounter==0) {
1646  cout<<"Use MassFitter method before Background"<<endl;
1647  return;
1648  }
1649  Double_t min=fMass-nOfSigma*fSigmaSgn;
1650  Double_t max=fMass+nOfSigma*fSigmaSgn;
1651 
1652  Background(min,max,background,errbackground);
1653 
1654  return;
1655 
1656 }
1657 //___________________________________________________________________________
1658 
1659 void AliHFMassFitterVAR::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1660  // Return background integral in a range
1661 
1662  if(fcounter==0) {
1663  cout<<"Use MassFitter method before Background"<<endl;
1664  return;
1665  }
1666 
1667  //functions names
1668  TString bkgname="funcbkgRecalc";
1669  TString bkg1name="funcbkg1Recalc";
1670 
1671  TF1 *funcbkg=0;
1672  if(ftypeOfFit4Sgn == 1) funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1673  else funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1674  if(!funcbkg){
1675  cout<<"AliHFMassFitterVAR::Background() ERROR -> Bkg function not found!"<<endl;
1676  return;
1677  }
1678 
1679 
1680  Double_t intB,intBerr;
1681 
1682  //relative error evaluation: from final parameters of the fit
1683  if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1684  else{
1685  intB=funcbkg->GetParameter(0);
1686  intBerr=funcbkg->GetParError(0);
1687  cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1688  }
1689 
1690  //relative error evaluation: from histo
1691 
1692  intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1693  Double_t sum2=0;
1694  for(Int_t i=1;i<=fSideBandl;i++){
1695  sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1696  }
1697  for(Int_t i=fSideBandr;i<=fNbin;i++){
1698  sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1699  }
1700 
1701  intBerr=TMath::Sqrt(sum2);
1702  cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1703 
1704  cout<<"Last estimation of bkg error is used"<<endl;
1705 
1706  //backround +/- error in the range
1707  if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1708  background = 0;
1709  errbackground = 0;
1710  }
1711  else{
1712  background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1713  errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1714  //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1715  }
1716  return;
1717 
1718 }
1719 
1720 
1721 //__________________________________________________________________________
1722 
1723 void AliHFMassFitterVAR::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
1724  // Return significance in mean+- n sigma
1725 
1726  Double_t min=fMass-nOfSigma*fSigmaSgn;
1727  Double_t max=fMass+nOfSigma*fSigmaSgn;
1728  Significance(min, max, significance, errsignificance);
1729 
1730  return;
1731 }
1732 
1733 //__________________________________________________________________________
1734 
1735 void AliHFMassFitterVAR::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1736  // Return significance integral in a range
1737 
1738  Double_t signal,errsignal,background,errbackground;
1739  Signal(min, max,signal,errsignal);
1740  Background(min, max,background,errbackground);
1741 
1742  if (signal+background <= 0.){
1743  cout<<"Cannot calculate significance because of div by 0!"<<endl;
1744  significance=-1;
1745  errsignificance=0;
1746  return;
1747  }
1748 
1749  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,background,errbackground,significance,errsignificance);
1750 
1751  return;
1752 }
1753 
1754 void AliHFMassFitterVAR::ComputeParSize() {
1755 
1756  //compute the size of the parameter array and set the data member
1757 
1758  switch (ftypeOfFit4Bkg) {//npar background func
1759  case 0:
1760  fParsSize = 2*3;
1761  break;
1762  case 1:
1763  fParsSize = 2*3;
1764  break;
1765  case 2:
1766  fParsSize = 3*3;
1767  break;
1768  case 3:
1769  fParsSize = 1*3;
1770  break;
1771  case 4:
1772  fParsSize = 2*3;
1773  break;
1774  case 5:
1775  fParsSize = 3*3;
1776  break;
1777  case 6:
1778  fParsSize = (fpolbackdegreeTay+1)*3;
1779  break;
1780 
1781  default:
1782  cout<<"AliHFMassFitterVAR, Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
1783  break;
1784  }
1785 
1786  fParsSize += 3; // npar refl
1787  fParsSize += 3; // npar signal gaus
1788 
1789  fParsSize*=2; // add errors
1790  cout<<"Parameters array size "<<fParsSize<<endl;
1791 }
1792 
1793 
1794 void AliHFMassFitterVAR::ComputeNFinalPars() {
1795 
1796  //compute the number of parameters of the total (signal+bgk) function
1797  cout<<"Info:ComputeNFinalPars... ";
1798  delete [] fBackParNames;
1799  delete [] fFixParBack;
1800  delete [] fFixParBackExternalValue;
1801  delete [] fparBackFixExt;
1802 
1803  delete [] fSignParNames;
1804  delete [] fFixParSign;
1805  delete [] fFixParSignExternalValue;
1806  delete [] fparSignFixExt;
1807 
1808  delete [] fReflParNames;
1809  delete [] fFixParRefl;
1810  delete [] fFixParReflExternalValue;
1811  delete [] fparReflFixExt;
1812 
1813 
1814  switch (ftypeOfFit4Bkg) {//npar background func
1815  case 0:
1816  fNparBack=2;
1817  fBackParNames=new TString[fNparBack];
1818  fBackParNames[0]="BkgInt";
1819  fBackParNames[1]="Slope";
1820  break;
1821  case 1:
1822  fNparBack=2;
1823  fBackParNames=new TString[fNparBack];
1824  fBackParNames[0]="BkgInt";
1825  fBackParNames[1]="Slope";
1826  break;
1827  case 2:
1828  fNparBack=3;
1829  fBackParNames=new TString[fNparBack];
1830  fBackParNames[0]="BkgInt";
1831  fBackParNames[1]="Coef1";
1832  fBackParNames[2]="Coef2";
1833  break;
1834  case 3:
1835  fNparBack=1;
1836  fBackParNames=new TString[fNparBack];
1837  fBackParNames[0]="Const";
1838  break;
1839  case 4:
1840  fNparBack=2;
1841  fBackParNames=new TString[fNparBack];
1842  fBackParNames[0]="BkgInt";
1843  fBackParNames[1]="Coef1";
1844  break;
1845  case 5:
1846  fNparBack=3;
1847  fBackParNames=new TString[fNparBack];
1848  fBackParNames[0]="BkgInt";
1849  fBackParNames[1]="Coef1";
1850  fBackParNames[2]="Coef2";
1851  break;
1852  case 6:
1853  fNparBack=fpolbackdegreeTay+1;
1854  fBackParNames=new TString[fNparBack];
1855  for(Int_t j=0;j<fNparBack;j++){
1856  fBackParNames[j]=Form("Coef%d",j);
1857  }
1858  break;
1859  default:
1860  cout<<"AliHFMassFitterVAR, Error in computing fNparBack: check ftypeOfFit4Bkg"<<endl;
1861  break;
1862  }
1863 
1864  fFixParBack=new Bool_t[fNparBack];
1865  fFixParBackExternalValue=new Bool_t[fNparBack];
1866  fparBackFixExt=new Double_t[fNparBack];
1867 
1868  for(Int_t ib=0;ib<fNparBack;ib++){
1869  fFixParBack[ib]=kFALSE;
1870  fFixParBackExternalValue[ib]=kFALSE;
1871  fparBackFixExt[ib]=0.;
1872  }
1873 
1874  fNparSignal=3;
1875  fSignParNames=new TString[3];
1876  fSignParNames[0]="Signal";
1877  fSignParNames[1]="Mean";
1878  fSignParNames[2]="Sigma";
1879 
1880  fFixParSign=new Bool_t[fNparSignal];
1881  fFixParSignExternalValue=new Bool_t[fNparSignal];
1882  fparSignFixExt=new Double_t[fNparSignal];
1883 
1884  for(Int_t ib=0;ib<fNparSignal;ib++){
1885  fFixParSign[ib]=kFALSE;
1886  fFixParSignExternalValue[ib]=kFALSE;
1887  fparSignFixExt[ib]=0.;
1888  }
1889 
1890 
1891  fNparRefl=0;
1892  if(ftypeOfFit4Sgn==2){
1893  fNparRefl=1;
1894  fReflParNames=new TString[fNparRefl];
1895  fReflParNames[0]="ReflOverS";
1896  }
1897 
1898 
1899  if(fNparRefl>0){
1900  fFixParRefl=new Bool_t[fNparRefl];
1901  fFixParReflExternalValue=new Bool_t[fNparRefl];
1902  fparReflFixExt=new Double_t[fNparRefl];
1903  }
1904  else {
1905  fFixParRefl=0x0;
1906  fparReflFixExt=0x0;
1907  fparReflFixExt=0x0;
1908  }
1909 
1910 
1911  for(Int_t ib=0;ib<fNparRefl;ib++){
1912  fFixParRefl[ib]=kFALSE;
1913  fFixParReflExternalValue[ib]=kFALSE;
1914  fparReflFixExt[ib]=0.;
1915  }
1916 
1917  fNFinalPars=fNparBack+fNparSignal+fNparRefl;
1918 
1919  cout<<": "<<fNFinalPars<<endl;
1920 
1921  Printf("Total number of par: %d, Back:%d, Sign:%d, Refl: %d",fNFinalPars,fNparBack,fNparSignal,fNparRefl);
1922 }
1923 
1924 
1925 TPaveText* AliHFMassFitterVAR::GetFitParametersBox(Double_t nsigma,Int_t mode){
1926  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1927  pinfom->SetBorderSize(0);
1928  pinfom->SetFillStyle(0);
1929  TF1* ff=fhistoInvMass->GetFunction("funcmass");
1930 
1931 
1932  for (Int_t i=fNparBack;i<fNparBack+fNparSignal;i++){
1933  pinfom->SetTextColor(kBlue);
1934  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1935  if(!(mode==1 && i==fNparBack)) pinfom->AddText(str);
1936  }
1937  for (Int_t i=fNparBack+fNparSignal;i<fNparBack+fNparSignal+fNparRefl;i++){
1938  pinfom->SetTextColor(kBlue+3);
1939  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1940  pinfom->AddText(str);
1941  }
1942 
1943  if(mode>1){
1944  pinfom->SetTextColor(kBlue+5);
1945  for (Int_t i=0;i<fNparBack;i++){
1946  TString str=Form("%s = %.3f #pm %.3f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1947  pinfom->AddText(str);
1948  }
1949  }
1950  pinfom->AddText(Form("#chi^{2}/NDF=%.2f/%d",ff->GetChisquare(),ff->GetNDF()));
1951 
1952  return pinfom;
1953 }
1954 
1955 
1956 TPaveText* AliHFMassFitterVAR::GetYieldBox(Double_t nsigma){
1957  TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1958  pinfo2->SetBorderSize(0);
1959  pinfo2->SetFillStyle(0);
1960 
1961  Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1962 
1963  Signal(nsigma,signal,errsignal);
1964  Background(nsigma,bkg, errbkg);
1965  AliVertexingHFUtils::ComputeSignificance(signal,errsignal,bkg,errbkg,signif,errsignif);
1966 
1967 
1968  TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1969  pinfo2->AddText(str);
1970  str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1971  pinfo2->AddText(str);
1972  str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1973  pinfo2->AddText(str);
1974  if(bkg>0) str=Form("S/B (%.0f#sigma) %.4f ",nsigma,signal/bkg);
1975  pinfo2->AddText(str);
1976 
1977  return pinfo2;
1978 
1979 }
1980 
1981 
1982 
1983 TH1F* AliHFMassFitterVAR::GetOverBackgroundResidualsAndPulls(Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
1984 
1985 
1986  if(!fhistoInvMass){
1987  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, invariant mass histogram not avaialble!");
1988  return 0x0;
1989  }
1990 
1991  TF1 *fback=fhistoInvMass->GetFunction("funcbkgRecalc");
1992  if(!fback){
1993  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, funcbkgRecalc not available!");
1994  return 0x0;
1995  }
1996 
1997  // THIS LONG WAY TO CP THE FUNC IS NEEDED ONLY TO EXTEND THE RANGE OF THE FUNCTION: NOT POSSIBLE OTHERWISE (WHY??? REALLY UNCOMFORTABLE)
1998  TF1 *fbackCp;
1999  if(ftypeOfFit4Sgn<2)fbackCp=new TF1("ftmpback",this,&AliHFMassFitterVAR::FitFunction4Bkg,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1), fNparBack,"AliHFMassFitterVAR","FitFunction4Bkg");
2000  else fbackCp=new TF1("ftmpback",this,&AliHFMassFitterVAR::FitFunction4BkgAndReflDraw,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1),fNparBack+fNparRefl,"AliHFMassFitterVAR","FitFunction4BkgAndReflDraw");
2001 
2002  for(Int_t i=0;i<fback->GetNpar();i++){
2003  fbackCp->SetParameter(i,fback->GetParameter(i));
2004  }
2005 
2006 
2007  TH1F *h=GetResidualsAndPulls(fhistoInvMass,fbackCp,minrange,maxrange,hPulls,hResidualTrend,hPullsTrend);
2008  delete fbackCp;
2009 
2010  if(fSigmaSgn<0){
2011  Printf("AliHFMassFitter::GetOverBackgroundResidualsAndPulls, negative sigma: fit not performed or something went wrong, cannto draw gaussian term on top of residuals");
2012  return h;
2013  }
2014 
2015  if(hResidualTrend){
2016  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));
2017  fgauss->SetParameter(0,fRawYield*fhistoInvMass->GetBinWidth(1));
2018  fgauss->SetParameter(1,fMass);
2019  fgauss->SetParameter(2,fSigmaSgn);
2020  fgauss->SetLineColor(kBlue);
2021  hResidualTrend->GetListOfFunctions()->Add(fgauss);
2022  }
2023  return h;
2024 }
2025 
2026 
2027 TH1F* AliHFMassFitterVAR::GetAllRangeResidualsAndPulls(Double_t minrange,Double_t maxrange,TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend){
2028 
2029 
2030  if(!fhistoInvMass){
2031  Printf("AliHFMassFitter::GetAllRangeResidualsAndPulls, invariant mass histogram not avaialble!");
2032  return 0x0;
2033  }
2034 
2035  TF1 *f=fhistoInvMass->GetFunction("funcmass");
2036  if(!f){
2037  Printf("AliHFMassFitter::GetAllRangeResidualsAndPulls, funcmass not available!");
2038  return 0x0;
2039  }
2040 
2041  // THIS LONG WAY TO CP THE FUNC IS NEEDED ONLY TO EXTEND THE RANGE OF THE FUNCTION: NOT POSSIBLE OTHERWISE (WHY??? REALLY UNCOMFORTABLE)
2042  TF1 *fmassCp=new TF1("fmassCp",this,&AliHFMassFitterVAR::FitFunction4MassDistr,fhistoInvMass->GetBinLowEdge(1),fhistoInvMass->GetBinLowEdge(fhistoInvMass->GetNbinsX()+1), fNFinalPars,"AliHFMassFitterVAR","FitFunction4MassDistr");
2043  for(Int_t i=0;i< f->GetNpar();i++){
2044  fmassCp->SetParameter(i,f->GetParameter(i));
2045  }
2046 
2047  TH1F *h=GetResidualsAndPulls(fhistoInvMass,fmassCp,minrange,maxrange,hPulls,hResidualTrend,hPullsTrend);
2048  delete fmassCp;
2049  return h;
2050 
2051 }
AliHFMassFitter::GetPad
TVirtualPad * GetPad(Double_t nsigma=3, Int_t writeFitInfo=1) const
Definition: AliHFMassFitter.cxx:1762
AliHFMassFitterVAR::Significance
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
backgournd in (min, max) with error
Definition: AliHFMassFitterVAR.cxx:1723
AliHFMassFitter::GetHistoClone
TH1F * GetHistoClone() const
Definition: AliHFMassFitter.cxx:1658
AliHFMassFitterVAR::fReflInit
Double_t fReflInit
smoothing refl template
Definition: AliHFMassFitterVAR.h:152
AliHFMassFitter::fFitPars
Float_t * fFitPars
number of parameters of the final function
Definition: AliHFMassFitter.h:151
filename
const char * filename
Definition: TestFCM.C:1
AliHFMassFitterVAR::Background
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
signal in (min, max) with error
Definition: AliHFMassFitterVAR.cxx:1642
str
return jsonbuilder str().c_str()
AliHFMassFitterVAR::FitFunction4BkgAndReflDraw
Double_t FitFunction4BkgAndReflDraw(Double_t *x, Double_t *par)
Definition: AliHFMassFitterVAR.cxx:520
AliHFMassFitterVAR::GetYieldBox
TPaveText * GetYieldBox(Double_t nsigma=3.)
Definition: AliHFMassFitterVAR.cxx:1956
Print
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
AliHFMassFitterVAR::GetFitParametersBox
TPaveText * GetFitParametersBox(Double_t nsigma=3., Int_t mode=0)
Definition: AliHFMassFitterVAR.cxx:1925
AliHFMassFitterVAR::fSignParNames
TString * fSignParNames
template of reflection contribution
Definition: AliHFMassFitterVAR.h:148
AliHFMassFitterVAR::fparReflFixExt
Double_t * fparReflFixExt
external values to fix back parameters
Definition: AliHFMassFitterVAR.h:161
AliHFMassFitter::fSideBandr
Int_t fSideBandr
left side band limit (bin number)
Definition: AliHFMassFitter.h:166
Double_t
double Double_t
Definition: External.C:58
AliHFMassFitterVAR::PlotFitVAR
void PlotFitVAR(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1)
Definition: AliHFMassFitterVAR.cxx:1394
AliHFMassFitterVAR::fNparBack
Int_t fNparBack
number of signal parameters
Definition: AliHFMassFitterVAR.h:145
AliHFMassFitterVAR::fFixParReflExternalValue
Bool_t * fFixParReflExternalValue
fix signal parameter from ext value
Definition: AliHFMassFitterVAR.h:158
AliHFMassFitterVAR::RefitWithBkgOnly
Bool_t RefitWithBkgOnly(Bool_t draw=kTRUE)
Definition: AliHFMassFitterVAR.cxx:1116
AliHFMassFitterVAR::WriteCanvas
void WriteCanvas(TString userIDstring="", TString path="./", Double_t nsigma=3, Int_t writeFitInfo=1, Bool_t draw=kFALSE) const
setters
Definition: AliHFMassFitterVAR.cxx:1338
AliHFMassFitterVAR::fFixParBack
Bool_t * fFixParBack
fix signal parameter from ext value
Definition: AliHFMassFitterVAR.h:154
AliHFMassFitter::DrawFit
void DrawFit(Double_t nsigma=3) const
Definition: AliHFMassFitter.cxx:1898
AliHFMassFitterVAR::Signal
void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
return total integral of the histogram
Definition: AliHFMassFitterVAR.cxx:1553
AliHFMassFitterVAR::FitFunction4Sgn
Double_t FitFunction4Sgn(Double_t *x, Double_t *par)
Definition: AliHFMassFitterVAR.cxx:506
mass
Double_t mass
Definition: charmCutsOptimization.C:35
AliHFMassFitterVAR::DrawHere
void DrawHere(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1)
Definition: AliHFMassFitterVAR.cxx:1386
AliHFMassFitterVAR::SetBackHighPolDegree
void SetBackHighPolDegree(Int_t deg)
Definition: AliHFMassFitterVAR.cxx:326
AliVertexingHFUtils.h
AliHFMassFitter::fParsSize
Int_t fParsSize
number of bins
Definition: AliHFMassFitter.h:149
AliHFMassFitter::ftypeOfFit4Sgn
Int_t ftypeOfFit4Sgn
0 = exponential; 1 = linear; 2 = pol2
Definition: AliHFMassFitter.h:154
AliHFMassFitterVAR::fNparRefl
Int_t fNparRefl
number of bkg parameters
Definition: AliHFMassFitterVAR.h:146
AliHFMassFitterVAR::fFixParBackExternalValue
Bool_t * fFixParBackExternalValue
fix signal parameter from ext value
Definition: AliHFMassFitterVAR.h:157
AliHFMassFitter::fmaxMass
Double_t fmaxMass
lower mass limit
Definition: AliHFMassFitter.h:145
AliHFMassFitterVAR::GetReflOverSignal
Double_t GetReflOverSignal(Double_t &err) const
significance in (min, max) with error
Definition: AliHFMassFitterVAR.cxx:1571
AliHFMassFitter::ftypeOfFit4Bkg
Int_t ftypeOfFit4Bkg
signal+background (kTRUE) or signal only (kFALSE)
Definition: AliHFMassFitter.h:153
c
TCanvas * c
Definition: TestFitELoss.C:172
AliHFMassFitterVAR::MassFitter
Bool_t MassFitter(Bool_t draw=kTRUE)
Definition: AliHFMassFitterVAR.cxx:792
AliHFMassFitterVAR::fFixParRefl
Bool_t * fFixParRefl
fix signal parameter from ext value
Definition: AliHFMassFitterVAR.h:155
AliHFMassFitterVAR::fReflParNames
TString * fReflParNames
back parameter names
Definition: AliHFMassFitterVAR.h:150
AliHFMassFitter::fMassErr
Double_t fMassErr
signal gaussian mean value
Definition: AliHFMassFitter.h:158
AliHFMassFitter::fmaxBinMass
Int_t fmaxBinMass
bin corresponding to fminMass
Definition: AliHFMassFitter.h:147
AliHFMassFitter::fminMass
Double_t fminMass
histogram to fit
Definition: AliHFMassFitter.h:144
AliHFMassFitterVAR
AliHFMassFitterVAR for the fit of invariant mass distribution of charmed mesons.
Definition: AliHFMassFitterVAR.h:28
AliHFMassFitter::fNbin
Int_t fNbin
bin corresponding to fmaxMass
Definition: AliHFMassFitter.h:148
Int_t
int Int_t
Definition: External.C:63
AliHFMassFitter::fRawYieldErr
Double_t fRawYieldErr
signal gaussian integral
Definition: AliHFMassFitter.h:162
AliHFMassFitter::fSigmaSgnErr
Double_t fSigmaSgnErr
signal gaussian sigma
Definition: AliHFMassFitter.h:160
Float_t
float Float_t
Definition: External.C:68
AliHFMassFitter::fMass
Double_t fMass
contains fit parameters
Definition: AliHFMassFitter.h:157
AliHFMassFitter::fcounter
Int_t fcounter
right side band limit (bin number)
Definition: AliHFMassFitter.h:167
AliHFMassFitterVAR::fRawYieldHelp
Double_t fRawYieldHelp
external values to fix refl parameters
Definition: AliHFMassFitterVAR.h:162
AliHFMassFitterVAR::fpolbackdegreeTay
Int_t fpolbackdegreeTay
internal variable used when fitting with reflections
Definition: AliHFMassFitterVAR.h:163
AliHFMassFitter::fNFinalPars
Int_t fNFinalPars
size of fFitPars array
Definition: AliHFMassFitter.h:150
nsigma
Double_t nsigma
Definition: charmCutsOptimization.C:29
AliHFMassFitterVAR::SetTemplateReflections
TH1F * SetTemplateReflections(const TH1F *h, TString option="templ", Double_t minRange=1.72, Double_t maxRange=2.05)
Definition: AliHFMassFitterVAR.cxx:337
mode
Int_t mode
Definition: anaM.C:41
AliHFMassFitterVAR::BackFitFuncPolHelper
Double_t BackFitFuncPolHelper(Double_t *x, Double_t *par)
Definition: AliHFMassFitterVAR.cxx:551
AliHFMassFitter::fRawYield
Double_t fRawYield
err signal gaussian sigma
Definition: AliHFMassFitter.h:161
AliHFMassFitterVAR::fBackParNames
TString * fBackParNames
signal parameter names
Definition: AliHFMassFitterVAR.h:149
AliHFMassFitterVAR::FitFunction4Bkg
Double_t FitFunction4Bkg(Double_t *x, Double_t *par)
Definition: AliHFMassFitterVAR.cxx:563
AliHFMassFitterVAR::IntS
void IntS(Float_t *valuewitherror) const
Definition: AliHFMassFitterVAR.cxx:1197
AliHFMassFitter::GetResidualsAndPulls
TH1F * GetResidualsAndPulls(TH1 *h, TF1 *f, Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Definition: AliHFMassFitter.cxx:2126
AliHFMassFitterVAR::fMassParticle
Double_t fMassParticle
help variable
Definition: AliHFMassFitterVAR.h:165
AliHFMassFitter::RebinMass
void RebinMass(Int_t bingroup=1)
Definition: AliHFMassFitter.cxx:556
ClassImp
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
Definition: AliAnalysisTaskDeltaPt.cxx:27
AliVertexingHFUtils::ComputeSignificance
static void ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance, Double_t &errsignificance)
Significance calculator.
Definition: AliVertexingHFUtils.cxx:73
draw
Bool_t draw[nPtBins]
Definition: ComputeDmesonYield.C:45
AliHFMassFitter::SetDefaultFixParam
virtual void SetDefaultFixParam()
Definition: AliHFMassFitter.cxx:331
AliHFMassFitter::fSigmaSgn
Double_t fSigmaSgn
err signal gaussian mean value
Definition: AliHFMassFitter.h:159
AliHFMassFitterVAR::FitFunction4Refl
Double_t FitFunction4Refl(Double_t *x, Double_t *par)
Definition: AliHFMassFitterVAR.cxx:535
AliHFMassFitter::fContourGraph
TList * fContourGraph
L, LW or Chi2.
Definition: AliHFMassFitter.h:170
AliHFMassFitterVAR::operator=
AliHFMassFitterVAR & operator=(const AliHFMassFitterVAR &mfit)
Definition: AliHFMassFitterVAR.cxx:237
AliHFMassFitterVAR::fparBackFixExt
Double_t * fparBackFixExt
external values to fix signal parameters
Definition: AliHFMassFitterVAR.h:160
AliHFMassFitterVAR::fSmoothRefl
Bool_t fSmoothRefl
refl parameter names
Definition: AliHFMassFitterVAR.h:151
rebin
Int_t rebin
Definition: charmCutsOptimization.C:32
AliHFMassFitterVAR::GetOverBackgroundResidualsAndPulls
TH1F * GetOverBackgroundResidualsAndPulls(Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Definition: AliHFMassFitterVAR.cxx:1983
AliHFMassFitterVAR::fFixParSign
Bool_t * fFixParSign
initial value of Refl/Signal
Definition: AliHFMassFitterVAR.h:153
AliHFMassFitter::fminBinMass
Int_t fminBinMass
upper mass limit
Definition: AliHFMassFitter.h:146
AliHFMassFitterVAR::fpolbackdegreeTayHelp
Int_t fpolbackdegreeTayHelp
degree of polynomial expansion for back fit (option 6 for back)
Definition: AliHFMassFitterVAR.h:164
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
Bool_t
bool Bool_t
Definition: External.C:53
AliHFMassFitterVAR::fFixParSignExternalValue
Bool_t * fFixParSignExternalValue
fix signal parameter from ext value
Definition: AliHFMassFitterVAR.h:156
AliHFMassFitterVAR::GetAllRangeResidualsAndPulls
TH1F * GetAllRangeResidualsAndPulls(Double_t minrange=0, Double_t maxrange=-1, TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0)
Definition: AliHFMassFitterVAR.cxx:2027
AliHFMassFitterVAR::FitFunction4MassDistr
Double_t FitFunction4MassDistr(Double_t *x, Double_t *par)
significance in (min, max) with error
Definition: AliHFMassFitterVAR.cxx:466
AliHFMassFitter::fSideBands
Bool_t fSideBands
err on signal gaussian integral
Definition: AliHFMassFitter.h:163
AliHFMassFitter.h
AliHFMassFitter::fFitOption
TString fFitOption
Number of points used in the fit.
Definition: AliHFMassFitter.h:169
AliHFMassFitter
AliHFMassFitter for the fit of invariant mass distribution of charmed mesons.
Definition: AliHFMassFitter.h:28
TH1
Definition: External.C:196
AliHFMassFitterVAR.h
AliHFMassFitterVAR::fhTemplRefl
TH1F * fhTemplRefl
number of reflection parameters
Definition: AliHFMassFitterVAR.h:147
AliHFMassFitterVAR::fparSignFixExt
Double_t * fparSignFixExt
fix signal parameter from ext value
Definition: AliHFMassFitterVAR.h:159
AliHFMassFitterVAR::PrepareHighPolFit
Bool_t PrepareHighPolFit(TF1 *fback)
Definition: AliHFMassFitterVAR.cxx:1054
AliHFMassFitter::operator=
AliHFMassFitter & operator=(const AliHFMassFitter &mfit)
Definition: AliHFMassFitter.cxx:214