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