]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliHFMassFitter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Use AddTaskXXX macros for analysis
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliHFMassFitter.cxx
CommitLineData
fabf4d8e 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/////////////////////////////////////////////////////////////
17//
18// AliHFMassFitter for the fit of invariant mass distribution
19// of charmed mesons
20//
21// Author: C.Bianchin, chiara.bianchin@pd.infn.it
22/////////////////////////////////////////////////////////////
23
24#include <TCanvas.h>
25
26#include "AliHFMassFitter.h"
27
28
29ClassImp(AliHFMassFitter)
30
31 //constructors
32AliHFMassFitter::AliHFMassFitter() : TNamed(),
33 fhistoInvMass(0),
34 fminMass(0),
35 fmaxMass(0),
36 fNbin(1),
37 fWithBkg(0),
38 ftypeOfFit4Bkg(0),
39 ftypeOfFit4Sgn(0),
40 ffactor(1),
41 fntuParam(0),
42 fMass(1.85),
43 fSigmaSgn(0.012),
44 fSideBands(0)
45{
46 // default constructor
47
48 cout<<"Default constructor"<<endl;
49}
50
51//___________________________________________________________________________
52
53AliHFMassFitter::AliHFMassFitter (TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
54 TNamed(),
55 fhistoInvMass(0),
56 fminMass(0),
57 fmaxMass(0),
58 fNbin(1),
59 fWithBkg(0),
60 ftypeOfFit4Bkg(0),
61 ftypeOfFit4Sgn(0),
62 ffactor(1),
63 fntuParam(0),
64 fMass(1.85),
65 fSigmaSgn(0.012),
66 fSideBands(0)
67{
68 // standard constructor
69
70 fhistoInvMass=histoToFit;
71 fminMass=minvalue;
72 fmaxMass=maxvalue;
73 if(rebin!=1) RebinMass(rebin);
74 else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
75 ftypeOfFit4Bkg=fittypeb;
76 ftypeOfFit4Sgn=fittypes;
77 if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
78 else fWithBkg=kTRUE;
79 if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
80 else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
81}
82
83//___________________________________________________________________________
84
85Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
86 // Fit function for signal+background
87
88
89 //exponential or linear fit
90 //
91 // par[0] = tot integral
92 // par[1] = slope
93 // par[2] = gaussian integral
94 // par[3] = gaussian mean
95 // par[4] = gaussian sigma
96
97 Double_t total,bkg=0,sgn=0;
98
99 if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
100 if(ftypeOfFit4Sgn == 0) {
101
102 Double_t parbkg[2] = {par[0]-par[2], par[1]};
103 bkg = FitFunction4Bkg(x,parbkg);
104 }
105 if(ftypeOfFit4Sgn == 1) {
106 Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
107 bkg = FitFunction4Bkg(x,parbkg);
108 }
109
110 sgn = FitFunction4Sgn(x,&par[2]);
111
112 }
113
114 //polynomial fit
115
116 // par[0] = tot integral
117 // par[1] = coef1
118 // par[2] = coef2
119 // par[3] = gaussian integral
120 // par[4] = gaussian mean
121 // par[5] = gaussian sigma
122
123 if (ftypeOfFit4Bkg==2) {
124
125 if(ftypeOfFit4Sgn == 0) {
126 //parbkg = new Double_t[2];
127 Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
128 bkg = FitFunction4Bkg(x,parbkg);
129 }
130 if(ftypeOfFit4Sgn == 1) {
131 //parbkg = new Double_t[6];
132 Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
133 bkg = FitFunction4Bkg(x,parbkg);
134 }
135
136 sgn = FitFunction4Sgn(x,&par[3]);
137 }
138
139 if (ftypeOfFit4Bkg==3) {
140
141 if(ftypeOfFit4Sgn == 0) {
142 bkg=FitFunction4Bkg(x,par);
143 sgn=FitFunction4Sgn(x,&par[1]);
144 }
145 if(ftypeOfFit4Sgn == 1) {
146 Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
147 bkg=FitFunction4Bkg(x,parbkg);
148 sgn=FitFunction4Sgn(x,&par[1]);
149 }
150 }
151
152 total = bkg + sgn;
153
154 return total;
155}
156
157//_________________________________________________________________________
158Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
159 // Fit function for the signal
160
161 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
162 //Par:
163 // * [0] = integralSgn
164 // * [1] = mean
165 // * [2] = sigma
166 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
167
168 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]);
169
170}
171
172//__________________________________________________________________________
173
174Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
175 // Fit function for the background
176
177 Double_t maxDeltaM = 4.*fSigmaSgn;
178 if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
179 TF1::RejectPoint();
180 return 0;
181 }
182 Int_t firstPar=0;
183 Double_t gaus2=0,total=-1;
184 if(ftypeOfFit4Sgn == 1){
185 firstPar=3;
186 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
187 //Par:
188 // * [0] = integralSgn
189 // * [1] = mean
190 // * [2] = sigma
191 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
192 gaus2 = FitFunction4Sgn(x,par);
193 }
194
195 switch (ftypeOfFit4Bkg){
196 case 0:
197 //exponential
198 //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
199 //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
200 //as integralTot- integralGaus (=par [2])
201 //Par:
202 // * [0] = integralBkg;
203 // * [1] = B;
204 //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
205 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]);
206 break;
207 case 1:
208 //linear
209 //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)
210 // * [0] = integralBkg;
211 // * [1] = b;
212 total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
213 break;
214 case 2:
215 //polynomial
216 //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
217 //+ 1/3*c*(max^3-min^3) ->
218 //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
219 // * [0] = integralBkg;
220 // * [1] = b;
221 // * [2] = c;
222 total = par[0+firstPar]/(fmaxMass-fminMass)+par[1]*(x[0]-0.5*(fmaxMass+fminMass))+par[2+firstPar]*(x[0]*x[0]-1/3.*(fmaxMass*fmaxMass*fmaxMass-fminMass*fminMass*fminMass)/(fmaxMass-fminMass));
223 break;
224 case 3:
225 total=par[0+firstPar];
226 break;
227// default:
228// Types of Fit Functions for Background:
229// * 0 = exponential;
230// * 1 = linear;
231// * 2 = polynomial 2nd order
232// * 3 = no background"<<endl;
233
234 }
235 return total+gaus2;
236}
237
238//__________________________________________________________________________
239
240void AliHFMassFitter::MassFitter(Bool_t draw){
241 // Main method of the class: performs the fit of the histogram
242
243 Int_t nFitPars=0; //total function's number of parameters
244 switch (ftypeOfFit4Bkg){
245 case 0:
246 nFitPars=5; //3+2
247 break;
248 case 1:
249 nFitPars=5; //3+2
250 break;
251 case 2:
252 nFitPars=6; //3+3
253 break;
254 case 3:
255 nFitPars=4; //3+1
256 break;
257 }
258
259 Int_t bkgPar = nFitPars-3; //background function's number of parameters
260
261 cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
262
263 //Total integral
264 Double_t totInt = fhistoInvMass->Integral("width");
265
266 fSideBands = kTRUE;
267 Double_t width;
268 Int_t binleft,binright;
269 fNbin=fhistoInvMass->GetNbinsX();
270 width=fhistoInvMass->GetBinWidth(8);
271 //width=(fmaxMass-fminMass)/(Double_t)fNbin;
272 binleft=(Int_t)((fMass-4.*fSigmaSgn-fminMass)/width);
273 binright=(Int_t)((fMass+4.*fSigmaSgn-fminMass)/width);
274
275 //sidebands integral - first approx (from histo)
276 Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,binleft,"width") + (Double_t)fhistoInvMass->Integral(binright,fNbin,"width");
277 cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<binleft<<"\tbinright = "<<binright<<endl;
278 cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
279
280 /*Fit Bkg*/
281
282 TF1 *funcbkg = new TF1("funcbkg",this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
283
284 funcbkg->SetLineColor(2); //red
285
286 //first fit for bkg: approx bkgint
287
288 switch (ftypeOfFit4Bkg) {
289 case 0: //gaus+expo
290 funcbkg->SetParNames("BkgInt","Slope");
291 funcbkg->SetParameters(sideBandsInt,-2.);
292 break;
293 case 1:
294 funcbkg->SetParNames("BkgInt","Slope");
295 funcbkg->SetParameters(sideBandsInt,-100.);
296 break;
297 case 2:
298 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
299 funcbkg->SetParameters(sideBandsInt,-10.,5);
300 break;
301 case 3:
302 if(ftypeOfFit4Sgn==0){
303 //in principle it doesn't have effects
304 funcbkg->SetParNames("Const");
305 funcbkg->SetParameter(0,0.);
306 funcbkg->FixParameter(0,0.);
307 }
308 break;
309 default:
310 cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
311 return;
312 break;
313 }
314 cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
315 Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
316
317 Double_t intbkg1=0,slope1=0,conc1=0;
318 //if only signal and reflection skip
319 if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
320 ftypeOfFit4Sgn=0;
321 fhistoInvMass->Fit("funcbkg","R,L,E,0");
322
323 for(Int_t i=0;i<bkgPar;i++){
324 fFitPars[i]=funcbkg->GetParameter(i);
325 //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
326 fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
327 //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
328 }
329
330 intbkg1 = funcbkg->GetParameter(0);
331 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
332 if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
333 cout<<"Primo fit: \nintbkg1 = "<<intbkg1<<"\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
334 } else cout<<"\t\t//"<<endl;
335
336 ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
337 TF1 *funcbkg1=0;
338 if (ftypeOfFit4Sgn == 1) {
339 cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
340 bkgPar+=3;
341
342 cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
343 funcbkg1 = new TF1("funcbkg1",this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
344 funcbkg1->SetLineColor(2); //red
345
346 if(ftypeOfFit4Bkg==2){
347 cout<<"*** Polynomial Fit ***"<<endl;
348 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
349 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
350// cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
351// cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
352 } else{
353 if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
354 {
355 cout<<"*** No background Fit ***"<<endl;
356 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
357 funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
358 funcbkg1->FixParameter(3,0.);
359 } else{ //expo or linear
360 if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
361 if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
362 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
363 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
364 }
365 }
366 fhistoInvMass->Fit("funcbkg1","R,L,E,+,0");
367
368 for(Int_t i=0;i<bkgPar;i++){
369 fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
370 //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
371 fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
372 //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
373
374 }
375
376 intbkg1=funcbkg1->GetParameter(3);
377 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
378 if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
379
380 } else {
381 bkgPar+=3;
382
383 for(Int_t i=0;i<3;i++){
384 fFitPars[bkgPar-3+i]=0.;
385 //cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
386 fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
387 //cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
388 }
389
390 for(Int_t i=0;i<bkgPar-3;i++){
391 fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
392 //cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
393 fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
394 //cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
395 }
396
397
398 }
399
400 //sidebands integral - second approx (from fit)
401 fSideBands = kFALSE;
402 Double_t bkgInt;
403
404 if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
405 else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
406 cout<<"------BkgInt(Fit) = "<<bkgInt<<endl;
407
408 //Signal integral - first approx
409 Double_t sgnInt;
410 sgnInt = totInt-bkgInt;
411 cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
412
413 /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
414 TF1 *funcmass = new TF1("funcmass",this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
415
416 funcmass->SetLineColor(4); //blue
417
418 //Set parameters
419 cout<<"\nTOTAL FIT"<<endl;
420
421 if(nFitPars==5){
422 funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
423 funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
424 // cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
425// cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
426 funcmass->FixParameter(0,totInt);
427 }
428 if (nFitPars==6){
429 funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
430 funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
431// cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
432// cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
433 funcmass->FixParameter(0,totInt);
434 }
435 if(nFitPars==4){
436 funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
437 funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
438 funcmass->FixParameter(0,0.);
439 //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
440 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
441
442 }
443
444 fhistoInvMass->Fit("funcmass","R,L,E,+,0");
445 cout<<"fit done"<<endl;
446
447 for(Int_t i=0;i<nFitPars;i++){
448 fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
449 fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
450 cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
451 }
452 /*
453 //check: cout parameters
454 for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
455 cout<<i<<"\t"<<fFitPars[i]<<endl;
456 }
457 */
458 /*
459 if(draw){
460 TCanvas *canvas=new TCanvas("canvas",fhistoInvMass->GetName());
461 TH1F *fhistocopy=new TH1F(*fhistoInvMass);
462 canvas->cd();
463 fhistocopy->Draw();
464 if(ftypeOfFit4Sgn == 1) {
465 cout<<"funcbkg1 "<<funcbkg1<<endl;
466 funcbkg1->Draw("sames");
467 }
468 else {
469 cout<<"funcbkg "<<funcbkg<<endl;
470 funcbkg->Draw("sames");
471 }
472 cout<<"funcmass "<<funcmass<<endl;
473 funcmass->Draw("sames");
474
475 }
476 */
477}
478
479//_________________________________________________________________________
480void AliHFMassFitter::GetFitPars(Float_t *vector25) const {
481 // Return fit parameters
482
483 for(Int_t i=0;i<25;i++){
484 vector25[i]=fFitPars[i];
485 }
486}
487
488//_________________________________________________________________________
489
490void AliHFMassFitter::InitNtuParam(char *ntuname) {
491 // Create ntuple to keep fit parameters
492
493 fntuParam=new TNtuple(ntuname,"Contains fit parameters","intbkg1:slope1:conc1:intGB:meanGB:sigmaGB:intbkg2:slope2:conc2:inttot:slope3:conc3:intsgn:meansgn:sigmasgn:intbkg1Err:slope1Err:conc1Err:intGBErr:meanGBErr:sigmaGBErr:intbkg2Err:slope2Err:conc2Err:inttotErr:slope3Err:conc3Err:intsgnErr:meansgnErr:sigmasgnErr");
494
495}
496
497//_________________________________________________________________________
498
499void AliHFMassFitter::FillNtuParam() {
500 // Fill ntuple with fit parameters
501
502 Float_t nothing=0.;
503
504 if (ftypeOfFit4Bkg==2) {
505 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
506 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
507 fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
508 fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
509 fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
510 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
511 fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
512 fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
513 fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
514 fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
515 fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
516 fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
517 fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
518 fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
519 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
520
521 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
522 fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
523 fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
524 fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
525 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
526 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
527 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
528 fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
529 fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
530 fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
531 fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
532 fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
533 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
534 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
535 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
536
537 } else {
538
539 if(ftypeOfFit4Bkg==3){
540 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
541 fntuParam->SetBranchAddress("slope1",&nothing);
542 fntuParam->SetBranchAddress("conc1",&nothing);
543 fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
544 fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
545 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
546 fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
547 fntuParam->SetBranchAddress("slope2",&nothing);
548 fntuParam->SetBranchAddress("conc2",&nothing);
549 fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
550 fntuParam->SetBranchAddress("slope3",&nothing);
551 fntuParam->SetBranchAddress("conc3",&nothing);
552 fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
553 fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
554 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
555
556 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
557 fntuParam->SetBranchAddress("slope1Err",&nothing);
558 fntuParam->SetBranchAddress("conc1Err",&nothing);
559 fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
560 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
561 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
562 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
563 fntuParam->SetBranchAddress("slope2Err",&nothing);
564 fntuParam->SetBranchAddress("conc2Err",&nothing);
565 fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
566 fntuParam->SetBranchAddress("slope3Err",&nothing);
567 fntuParam->SetBranchAddress("conc3Err",&nothing);
568 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
569 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
570 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
571
572 }
573 else{
574 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
575 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
576 fntuParam->SetBranchAddress("conc1",&nothing);
577 fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
578 fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
579 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
580 fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
581 fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
582 fntuParam->SetBranchAddress("conc2",&nothing);
583 fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
584 fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
585 fntuParam->SetBranchAddress("conc3",&nothing);
586 fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
587 fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
588 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
589
590 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
591 fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
592 fntuParam->SetBranchAddress("conc1Err",&nothing);
593 fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
594 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
595 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
596 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
597 fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
598 fntuParam->SetBranchAddress("conc2Err",&nothing);
599 fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
600 fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
601 fntuParam->SetBranchAddress("conc3Err",&nothing);
602 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
603 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
604 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
605 }
606
607 }
608 fntuParam->TTree::Fill();
609}
610
611//_________________________________________________________________________
612
613TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
614 // Create, fill and return ntuple with fit parameters
615
616 InitNtuParam(ntuname);
617 FillNtuParam();
618 return fntuParam;
619}
620//_________________________________________________________________________
621
622void AliHFMassFitter::RebinMass(Int_t bingroup){
623 // Rebin invariant mass histogram
624
625 if(bingroup<1){
626 cout<<"Error! Cannot group "<<bingroup<<" bins\n";
627 fNbin=fhistoInvMass->GetNbinsX();
628 cout<<"Kept original number of bins: "<<fNbin<<endl;
629 } else{
630 fhistoInvMass->Rebin(bingroup);
631 fNbin = fhistoInvMass->GetNbinsX();
632 cout<<"New number of bins: "<<fNbin<<endl;
633 }
634
635
636}
637
638//_________________________________________________________________________
639
640void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
641 // Return signal integral in mean+- n sigma
642
643 TF1 *funcmass=fhistoInvMass->GetFunction("funcmass");
644 Double_t mean=0,sigma=1;
645 Double_t intS,intSerr;
646
647 if(ftypeOfFit4Bkg == 2) { //polynomial
648 mean=funcmass->GetParameter(4); //mean
649 sigma=funcmass->GetParameter(5); //sigma
650 intS=fFitPars[12];
651 intSerr=fFitPars[27];
652 } else if(ftypeOfFit4Bkg == 3){ //no background
653 mean=funcmass->GetParameter(2); //mean
654 sigma=funcmass->GetParameter(3); //sigma
655 intS=fFitPars[6];
656 intSerr=fFitPars[15];
657 } else { //expo or linear
658 mean=funcmass->GetParameter(3); //mean
659 sigma=funcmass->GetParameter(4); //sigma
660 intS=fFitPars[9];
661 intSerr=fFitPars[21];
662 }
663
664 TF1 *funcbkg=0;
665 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction("funcbkg");
666 else funcbkg=fhistoInvMass->GetFunction("funcbkg1");
667 Double_t min=mean-nOfSigma*sigma;
668 Double_t max=mean+nOfSigma*sigma;
669 signal=funcmass->Integral(min,max)-funcbkg->Integral(min,max);
670 errsignal=intSerr/intS*signal; // assume relative error is the same as for total integral
671
672 return;
673}
674//_________________________________________________________________________
675
676void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
677 // Return background integral in mean+- n sigma
678
679 TF1 *funcmass=fhistoInvMass->GetFunction("funcmass");
680 Double_t mean=0,sigma=1;
681 Double_t intB,intBerr,err;
682
683 if(ftypeOfFit4Bkg == 2) { //polynomial
684 mean=funcmass->GetParameter(4); //mean
685 sigma=funcmass->GetParameter(5); //sigma
686 intB=fFitPars[9]-fFitPars[12];
687 intBerr=fFitPars[27];
688 } else if(ftypeOfFit4Bkg == 3){ //no background
689 mean=funcmass->GetParameter(2); //mean
690 sigma=funcmass->GetParameter(3); //sigma
691 if (ftypeOfFit4Sgn == 1){ //reflection
692 intB=fFitPars[1];
693 intBerr=fFitPars[9];
694 } else {
695 intB=-1;intBerr=-1;
696 err=0;
697 }
698 } else { //expo or linear
699 mean=funcmass->GetParameter(3); //mean
700 sigma=funcmass->GetParameter(4); //sigma
701 intB=fFitPars[7]-fFitPars[9];
702 intBerr=fFitPars[21];
703 }
704
705 TF1 *funcbkg=0;
706 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction("funcbkg");
707 else funcbkg=fhistoInvMass->GetFunction("funcbkg1");
708 Double_t min=mean-nOfSigma*sigma;
709 Double_t max=mean+nOfSigma*sigma;
710 background=funcbkg->Integral(min,max);
711 errbackground=intBerr/intB*background;
712
713 return;
714
715}
716
717//__________________________________________________________________________
718
719void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
720 // Return significance in mean+- n sigma
721
722 Double_t signal,errsignal,background,errbackground;
723 Signal(nOfSigma,signal,errsignal);
724 Background(nOfSigma,background,errbackground);
725
726 significance = signal/TMath::Sqrt(signal+background);
727 errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
728
729 return;
730}
731
732//__________________________________________________________________________