]> 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
New class for PID of HF candidates (R. Romita)
[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
dc262918 24#include <Riostream.h>
25#include <TMath.h>
26#include <TNtuple.h>
27#include <TH1F.h>
28#include <TF1.h>
29#include <TList.h>
30#include <TFile.h>
fabf4d8e 31#include <TCanvas.h>
ad73ca8f 32#include <TGraph.h>
33#include <TVirtualFitter.h>
dc262918 34#include <TMinuit.h>
35
fabf4d8e 36
37#include "AliHFMassFitter.h"
38
39
dc262918 40
fabf4d8e 41ClassImp(AliHFMassFitter)
42
56cbeefb 43 //************** constructors
2f328d65 44AliHFMassFitter::AliHFMassFitter() :
45 TNamed(),
fabf4d8e 46 fhistoInvMass(0),
47 fminMass(0),
48 fmaxMass(0),
34c79b83 49 fNbin(0),
2f328d65 50 fParsSize(1),
51 fFitPars(0),
fabf4d8e 52 fWithBkg(0),
53 ftypeOfFit4Bkg(0),
54 ftypeOfFit4Sgn(0),
55 ffactor(1),
56 fntuParam(0),
57 fMass(1.85),
58 fSigmaSgn(0.012),
56cbeefb 59 fSideBands(0),
b7d4bc49 60 fSideBandl(0),
61 fSideBandr(0),
ad73ca8f 62 fcounter(0),
63 fContourGraph(0)
fabf4d8e 64{
65 // default constructor
66
67 cout<<"Default constructor"<<endl;
ad73ca8f 68
fabf4d8e 69}
70
71//___________________________________________________________________________
72
dc262918 73AliHFMassFitter::AliHFMassFitter (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
fabf4d8e 74 TNamed(),
75 fhistoInvMass(0),
76 fminMass(0),
77 fmaxMass(0),
34c79b83 78 fNbin(0),
2f328d65 79 fParsSize(1),
80 fFitPars(0),
fabf4d8e 81 fWithBkg(0),
82 ftypeOfFit4Bkg(0),
83 ftypeOfFit4Sgn(0),
84 ffactor(1),
85 fntuParam(0),
86 fMass(1.85),
87 fSigmaSgn(0.012),
56cbeefb 88 fSideBands(0),
b7d4bc49 89 fSideBandl(0),
90 fSideBandr(0),
ad73ca8f 91 fcounter(0),
92 fContourGraph(0)
fabf4d8e 93{
94 // standard constructor
95
56cbeefb 96 fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
6321ee46 97 fhistoInvMass->SetDirectory(0);
fabf4d8e 98 fminMass=minvalue;
99 fmaxMass=maxvalue;
100 if(rebin!=1) RebinMass(rebin);
101 else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
102 ftypeOfFit4Bkg=fittypeb;
103 ftypeOfFit4Sgn=fittypes;
104 if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
105 else fWithBkg=kTRUE;
106 if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
107 else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
2f328d65 108
109 ComputeParSize();
110 fFitPars=new Float_t[fParsSize];
ad73ca8f 111
112 fContourGraph=new TList();
113 fContourGraph->SetOwner();
114
fabf4d8e 115}
116
74179930 117
118AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
119 TNamed(),
120 fhistoInvMass(mfit.fhistoInvMass),
121 fminMass(mfit.fminMass),
122 fmaxMass(mfit.fmaxMass),
123 fNbin(mfit.fNbin),
2f328d65 124 fParsSize(mfit.fParsSize),
125 fFitPars(0),
74179930 126 fWithBkg(mfit.fWithBkg),
127 ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
128 ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
129 ffactor(mfit.ffactor),
130 fntuParam(mfit.fntuParam),
131 fMass(mfit.fMass),
132 fSigmaSgn(mfit.fSigmaSgn),
56cbeefb 133 fSideBands(mfit.fSideBands),
b7d4bc49 134 fSideBandl(mfit.fSideBandl),
135 fSideBandr(mfit.fSideBandr),
ad73ca8f 136 fcounter(mfit.fcounter),
137 fContourGraph(mfit.fContourGraph)
74179930 138{
139 //copy constructor
2f328d65 140
141 if(mfit.fParsSize > 0){
142 fFitPars=new Float_t[fParsSize];
143 memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
144 }
145 //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
146}
147
148//_________________________________________________________________________
149
150AliHFMassFitter::~AliHFMassFitter() {
dc262918 151
152 //destructor
153
56cbeefb 154 cout<<"AliHFMassFitter destructor called"<<endl;
a60f573d 155 if(fhistoInvMass) {
156 cout<<"deleting histogram..."<<endl;
157 delete fhistoInvMass;
158 fhistoInvMass=NULL;
159 }
160 if(fntuParam){
161 cout<<"deleting ntuple..."<<endl;
162 delete fntuParam;
163
164 fntuParam=NULL;
165 }
166
167 if(fFitPars) {
168 delete[] fFitPars;
169 cout<<"deleting parameter array..."<<endl;
170 fFitPars=NULL;
171 }
172
56cbeefb 173 fcounter = 0;
74179930 174}
175
176//_________________________________________________________________________
177
178AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
179
180 //assignment operator
181
182 if(&mfit == this) return *this;
183 fhistoInvMass= mfit.fhistoInvMass;
184 fminMass= mfit.fminMass;
185 fmaxMass= mfit.fmaxMass;
186 fNbin= mfit.fNbin;
2f328d65 187 fParsSize= mfit.fParsSize;
74179930 188 fWithBkg= mfit.fWithBkg;
189 ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
190 ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
191 ffactor= mfit.ffactor;
192 fntuParam= mfit.fntuParam;
193 fMass= mfit.fMass;
194 fSigmaSgn= mfit.fSigmaSgn;
195 fSideBands= mfit.fSideBands;
b7d4bc49 196 fSideBandl= mfit.fSideBandl;
197 fSideBandr= mfit.fSideBandr;
56cbeefb 198 fcounter= mfit.fcounter;
ad73ca8f 199 fContourGraph= mfit.fContourGraph;
200
2f328d65 201 if(mfit.fParsSize > 0){
a60f573d 202 if(fFitPars) {
203 delete[] fFitPars;
204 fFitPars=NULL;
205 }
2f328d65 206 fFitPars=new Float_t[fParsSize];
207 memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
208 }
209// fFitPars=new Float_t[fParsSize];
210// for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
74179930 211
212 return *this;
213}
56cbeefb 214
215//************ tools & settings
216
2f328d65 217//__________________________________________________________________________
218
219void AliHFMassFitter::ComputeParSize() {
dc262918 220
221 //compute the size of the parameter array and set the data member
222
2f328d65 223 switch (ftypeOfFit4Bkg) {//npar backround func
224 case 0:
225 fParsSize = 2*3;
226 break;
227 case 1:
228 fParsSize = 2*3;
229 break;
230 case 2:
231 fParsSize = 3*3;
232 break;
233 case 3:
56cbeefb 234 fParsSize = 1*3;
2f328d65 235 break;
236 default:
56cbeefb 237 cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
2f328d65 238 break;
239 }
240
56cbeefb 241 fParsSize += 3; // npar refl
242 fParsSize += 3; // npar signal gaus
2f328d65 243
56cbeefb 244 fParsSize*=2; // add errors
2f328d65 245 cout<<"Parameters array size "<<fParsSize<<endl;
246}
247
56cbeefb 248//___________________________________________________________________________
dc262918 249void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
a60f573d 250 //fhistoInvMass = (TH1F*)histoToFit->Clone();
251 fhistoInvMass = new TH1F(*histoToFit);
6321ee46 252 fhistoInvMass->SetDirectory(0);
a60f573d 253 cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
56cbeefb 254}
255
2f328d65 256//___________________________________________________________________________
257
258 void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
dc262918 259
260 //set the type of fit to perform for signal and background
261
2f328d65 262 ftypeOfFit4Bkg = fittypeb;
263 ftypeOfFit4Sgn = fittypes;
264 ComputeParSize();
265 fFitPars = new Float_t[fParsSize];
266
267}
74179930 268
fabf4d8e 269//___________________________________________________________________________
270
56cbeefb 271void AliHFMassFitter::Reset() {
dc262918 272
273 //delete the histogram and reset the mean and sigma to default
274
a60f573d 275 cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
0566386c 276 fMass=1.85;
277 fSigmaSgn=0.012;
a60f573d 278 cout<<"Reset "<<fhistoInvMass<<endl;
279 if(fhistoInvMass) {
ad73ca8f 280 //cout<<"esiste"<<endl;
6321ee46 281 delete fhistoInvMass;
a60f573d 282 fhistoInvMass=NULL;
283 cout<<fhistoInvMass<<endl;
284 }
285 else cout<<"histogram doesn't exist, do not delete"<<endl;
286
56cbeefb 287}
288
289//_________________________________________________________________________
290
291void AliHFMassFitter::InitNtuParam(char *ntuname) {
dc262918 292
56cbeefb 293 // Create ntuple to keep fit parameters
dc262918 294
56cbeefb 295 fntuParam=0;
296 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");
297
298}
299
300//_________________________________________________________________________
301
302void AliHFMassFitter::FillNtuParam() {
303 // Fill ntuple with fit parameters
304
305 Float_t nothing=0.;
306
307 if (ftypeOfFit4Bkg==2) {
308 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
309 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
310 fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
311 fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
312 fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
313 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
314 fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
315 fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
316 fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
317 fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
318 fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
319 fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
320 fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
321 fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
322 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
323
324 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
325 fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
326 fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
327 fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
328 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
329 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
330 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
331 fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
332 fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
333 fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
334 fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
335 fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
336 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
337 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
338 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
339
340 } else {
341
342 if(ftypeOfFit4Bkg==3){
343 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
344 fntuParam->SetBranchAddress("slope1",&nothing);
345 fntuParam->SetBranchAddress("conc1",&nothing);
346 fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
347 fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
348 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
349 fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
350 fntuParam->SetBranchAddress("slope2",&nothing);
351 fntuParam->SetBranchAddress("conc2",&nothing);
352 fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
353 fntuParam->SetBranchAddress("slope3",&nothing);
354 fntuParam->SetBranchAddress("conc3",&nothing);
355 fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
356 fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
357 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
358
359 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
360 fntuParam->SetBranchAddress("slope1Err",&nothing);
361 fntuParam->SetBranchAddress("conc1Err",&nothing);
362 fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
363 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
364 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
365 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
366 fntuParam->SetBranchAddress("slope2Err",&nothing);
367 fntuParam->SetBranchAddress("conc2Err",&nothing);
368 fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
369 fntuParam->SetBranchAddress("slope3Err",&nothing);
370 fntuParam->SetBranchAddress("conc3Err",&nothing);
371 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
372 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
373 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
374
375 }
376 else{
377 fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
378 fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
379 fntuParam->SetBranchAddress("conc1",&nothing);
380 fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
381 fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
382 fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
383 fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
384 fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
385 fntuParam->SetBranchAddress("conc2",&nothing);
386 fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
387 fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
388 fntuParam->SetBranchAddress("conc3",&nothing);
389 fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
390 fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
391 fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
392
393 fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
394 fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
395 fntuParam->SetBranchAddress("conc1Err",&nothing);
396 fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
397 fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
398 fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
399 fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
400 fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
401 fntuParam->SetBranchAddress("conc2Err",&nothing);
402 fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
403 fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
404 fntuParam->SetBranchAddress("conc3Err",&nothing);
405 fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
406 fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
407 fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
408 }
409
410 }
411 fntuParam->TTree::Fill();
412}
413
414//_________________________________________________________________________
415
416TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
417 // Create, fill and return ntuple with fit parameters
418
419 InitNtuParam(ntuname);
420 FillNtuParam();
421 return fntuParam;
422}
423//_________________________________________________________________________
424
425void AliHFMassFitter::RebinMass(Int_t bingroup){
426 // Rebin invariant mass histogram
427
428 if(bingroup<1){
429 cout<<"Error! Cannot group "<<bingroup<<" bins\n";
430 fNbin=fhistoInvMass->GetNbinsX();
431 cout<<"Kept original number of bins: "<<fNbin<<endl;
432 } else{
433 fhistoInvMass->Rebin(bingroup);
434 fNbin = fhistoInvMass->GetNbinsX();
435 cout<<"New number of bins: "<<fNbin<<endl;
436 }
437
438
439}
440
441//************* fit
442
443//___________________________________________________________________________
444
fabf4d8e 445Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
446 // Fit function for signal+background
447
448
449 //exponential or linear fit
450 //
451 // par[0] = tot integral
452 // par[1] = slope
453 // par[2] = gaussian integral
454 // par[3] = gaussian mean
455 // par[4] = gaussian sigma
456
457 Double_t total,bkg=0,sgn=0;
458
459 if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
460 if(ftypeOfFit4Sgn == 0) {
461
462 Double_t parbkg[2] = {par[0]-par[2], par[1]};
463 bkg = FitFunction4Bkg(x,parbkg);
464 }
465 if(ftypeOfFit4Sgn == 1) {
466 Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
467 bkg = FitFunction4Bkg(x,parbkg);
468 }
469
470 sgn = FitFunction4Sgn(x,&par[2]);
471
472 }
473
474 //polynomial fit
475
476 // par[0] = tot integral
477 // par[1] = coef1
478 // par[2] = coef2
479 // par[3] = gaussian integral
480 // par[4] = gaussian mean
481 // par[5] = gaussian sigma
482
483 if (ftypeOfFit4Bkg==2) {
484
485 if(ftypeOfFit4Sgn == 0) {
56cbeefb 486
fabf4d8e 487 Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
488 bkg = FitFunction4Bkg(x,parbkg);
489 }
490 if(ftypeOfFit4Sgn == 1) {
56cbeefb 491
fabf4d8e 492 Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
16856d6e 493 bkg = FitFunction4Bkg(x,parbkg);
fabf4d8e 494 }
495
496 sgn = FitFunction4Sgn(x,&par[3]);
497 }
498
499 if (ftypeOfFit4Bkg==3) {
500
501 if(ftypeOfFit4Sgn == 0) {
502 bkg=FitFunction4Bkg(x,par);
503 sgn=FitFunction4Sgn(x,&par[1]);
504 }
505 if(ftypeOfFit4Sgn == 1) {
506 Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
507 bkg=FitFunction4Bkg(x,parbkg);
508 sgn=FitFunction4Sgn(x,&par[1]);
509 }
510 }
511
512 total = bkg + sgn;
513
514 return total;
515}
516
517//_________________________________________________________________________
518Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
519 // Fit function for the signal
520
521 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
522 //Par:
523 // * [0] = integralSgn
524 // * [1] = mean
525 // * [2] = sigma
526 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
527
528 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]);
529
530}
531
532//__________________________________________________________________________
533
534Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
535 // Fit function for the background
536
537 Double_t maxDeltaM = 4.*fSigmaSgn;
538 if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
539 TF1::RejectPoint();
540 return 0;
541 }
542 Int_t firstPar=0;
543 Double_t gaus2=0,total=-1;
544 if(ftypeOfFit4Sgn == 1){
545 firstPar=3;
546 //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
547 //Par:
548 // * [0] = integralSgn
549 // * [1] = mean
550 // * [2] = sigma
551 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
552 gaus2 = FitFunction4Sgn(x,par);
553 }
554
555 switch (ftypeOfFit4Bkg){
556 case 0:
557 //exponential
558 //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
559 //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
560 //as integralTot- integralGaus (=par [2])
561 //Par:
562 // * [0] = integralBkg;
563 // * [1] = B;
564 //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
565 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]);
566 break;
567 case 1:
568 //linear
569 //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)
570 // * [0] = integralBkg;
571 // * [1] = b;
572 total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
573 break;
574 case 2:
575 //polynomial
576 //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
577 //+ 1/3*c*(max^3-min^3) ->
578 //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
579 // * [0] = integralBkg;
580 // * [1] = b;
581 // * [2] = c;
582 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));
583 break;
584 case 3:
585 total=par[0+firstPar];
586 break;
587// default:
588// Types of Fit Functions for Background:
589// * 0 = exponential;
590// * 1 = linear;
591// * 2 = polynomial 2nd order
592// * 3 = no background"<<endl;
593
594 }
595 return total+gaus2;
596}
597
b7d4bc49 598//__________________________________________________________________________
599Bool_t AliHFMassFitter::SideBandsBounds(){
600
dc262918 601 //determines the ranges of the side bands
602
b7d4bc49 603 Double_t width=fhistoInvMass->GetBinWidth(8);
604 if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
605 Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
606 Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
607
608 if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
609 cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
610 return kFALSE;
611 }
612
dc262918 613 if((TMath::Abs(fminMass-minHisto) < 10e6 || TMath::Abs(fmaxMass - maxHisto) < 10e6) && (fMass-4.*fSigmaSgn-fminMass) < 10e6){
b7d4bc49 614 Double_t coeff = (fMass-fminMass)/fSigmaSgn;
615
616 fSideBandl=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
617 fSideBandr=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
618 cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
619 if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
620
621 if (coeff<2) {
622 cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
623 ftypeOfFit4Bkg=3;
624 //set binleft and right without considering SetRangeFit- anyway no bkg!
625 fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
626 fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
627 }
628 }
629 else {
630 fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
631 fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
a60f573d 632// cout<<"\tfMass = "<<fMass<<"\tfSigmaSgn = "<<fSigmaSgn<<"\tminHisto = "<<minHisto<<endl;
633// cout<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
b7d4bc49 634 }
635
636 if (fSideBandl==0) {
637 cout<<"Error! Range too little";
638 return kFALSE;
639 }
b7d4bc49 640 return kTRUE;
641}
642
16856d6e 643//__________________________________________________________________________
644
b7d4bc49 645void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
dc262918 646
647 // get the range of the side bands
648
b7d4bc49 649 if (fSideBandl==0 && fSideBandr==0){
650 cout<<"Use MassFitter method first"<<endl;
651 return;
652 }
653 left=fSideBandl;
654 right=fSideBandr;
655}
656
fabf4d8e 657//__________________________________________________________________________
658
34c79b83 659Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
fabf4d8e 660 // Main method of the class: performs the fit of the histogram
ad73ca8f 661
662 //Set default fitter Minuit in order to use gMinuit in the contour plots
663 TVirtualFitter::SetDefaultFitter("Minuit");
664
fabf4d8e 665 Int_t nFitPars=0; //total function's number of parameters
666 switch (ftypeOfFit4Bkg){
667 case 0:
668 nFitPars=5; //3+2
669 break;
670 case 1:
671 nFitPars=5; //3+2
672 break;
673 case 2:
674 nFitPars=6; //3+3
675 break;
676 case 3:
677 nFitPars=4; //3+1
678 break;
679 }
680
681 Int_t bkgPar = nFitPars-3; //background function's number of parameters
682
683 cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
684
ad73ca8f 685
686 TString listname="contourplot";
687 listname+=fcounter;
688 if(!fContourGraph){
689 fContourGraph=new TList();
690 fContourGraph->SetOwner();
691 }
692
693 fContourGraph->SetName(listname);
694
695
56cbeefb 696 //function names
697 TString bkgname="funcbkg";
698 TString bkg1name="funcbkg1";
699 TString massname="funcmass";
700
fabf4d8e 701 //Total integral
702 Double_t totInt = fhistoInvMass->Integral("width");
ad73ca8f 703 cout<<"Integral"<<endl;
fabf4d8e 704
705 fSideBands = kTRUE;
34c79b83 706 Double_t width=fhistoInvMass->GetBinWidth(8);
ad73ca8f 707 cout<<"fNbin"<<fNbin<<endl;
34c79b83 708 if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
709 Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
710 Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
b7d4bc49 711 Bool_t ok=SideBandsBounds();
712 if(!ok) return kFALSE;
34c79b83 713
fabf4d8e 714 //sidebands integral - first approx (from histo)
b7d4bc49 715 Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
716 cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
fabf4d8e 717 cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
b7d4bc49 718 if (sideBandsInt<=0) {
719 cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
720 return kFALSE;
721 }
722
fabf4d8e 723 /*Fit Bkg*/
724
16856d6e 725
34c79b83 726 TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
56cbeefb 727 cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
fabf4d8e 728
729 funcbkg->SetLineColor(2); //red
730
731 //first fit for bkg: approx bkgint
732
733 switch (ftypeOfFit4Bkg) {
734 case 0: //gaus+expo
735 funcbkg->SetParNames("BkgInt","Slope");
736 funcbkg->SetParameters(sideBandsInt,-2.);
737 break;
738 case 1:
739 funcbkg->SetParNames("BkgInt","Slope");
740 funcbkg->SetParameters(sideBandsInt,-100.);
741 break;
742 case 2:
743 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
744 funcbkg->SetParameters(sideBandsInt,-10.,5);
745 break;
746 case 3:
747 if(ftypeOfFit4Sgn==0){
fabf4d8e 748 funcbkg->SetParNames("Const");
749 funcbkg->SetParameter(0,0.);
750 funcbkg->FixParameter(0,0.);
751 }
752 break;
753 default:
754 cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
34c79b83 755 return kFALSE;
fabf4d8e 756 break;
757 }
758 cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
759 Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
760
761 Double_t intbkg1=0,slope1=0,conc1=0;
56cbeefb 762 //if only signal and reflection: skip
fabf4d8e 763 if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
764 ftypeOfFit4Sgn=0;
56cbeefb 765 fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
fabf4d8e 766
767 for(Int_t i=0;i<bkgPar;i++){
768 fFitPars[i]=funcbkg->GetParameter(i);
769 //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
770 fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
771 //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
772 }
b7d4bc49 773 fSideBands = kFALSE;
34c79b83 774 //intbkg1 = funcbkg->GetParameter(0);
775 funcbkg->SetRange(fminMass,fmaxMass);
776 intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
fabf4d8e 777 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
778 if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
6321ee46 779 cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
56cbeefb 780 }
781 else cout<<"\t\t//"<<endl;
fabf4d8e 782
783 ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
784 TF1 *funcbkg1=0;
785 if (ftypeOfFit4Sgn == 1) {
786 cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
787 bkgPar+=3;
788
56cbeefb 789 cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
790
791 funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
792 cout<<"Function name = "<<funcbkg1->GetName()<<endl;
793
fabf4d8e 794 funcbkg1->SetLineColor(2); //red
795
796 if(ftypeOfFit4Bkg==2){
797 cout<<"*** Polynomial Fit ***"<<endl;
798 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
799 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
6321ee46 800
801 //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
802 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
fabf4d8e 803 } else{
804 if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
805 {
806 cout<<"*** No background Fit ***"<<endl;
807 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
808 funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
809 funcbkg1->FixParameter(3,0.);
810 } else{ //expo or linear
811 if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
812 if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
813 funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
814 funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
815 }
816 }
6321ee46 817 Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
818 if (status != 0){
819 cout<<"Minuit returned "<<status<<endl;
820 return kFALSE;
821 }
822
fabf4d8e 823 for(Int_t i=0;i<bkgPar;i++){
824 fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
825 //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
826 fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
56cbeefb 827 //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
fabf4d8e 828 }
829
830 intbkg1=funcbkg1->GetParameter(3);
831 if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
832 if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
833
834 } else {
835 bkgPar+=3;
836
837 for(Int_t i=0;i<3;i++){
838 fFitPars[bkgPar-3+i]=0.;
56cbeefb 839 cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
fabf4d8e 840 fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
56cbeefb 841 cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
fabf4d8e 842 }
843
844 for(Int_t i=0;i<bkgPar-3;i++){
845 fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
56cbeefb 846 cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
fabf4d8e 847 fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
56cbeefb 848 cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
fabf4d8e 849 }
850
851
852 }
853
854 //sidebands integral - second approx (from fit)
855 fSideBands = kFALSE;
856 Double_t bkgInt;
b7d4bc49 857 cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
fabf4d8e 858 if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
859 else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
b7d4bc49 860 cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
fabf4d8e 861
862 //Signal integral - first approx
863 Double_t sgnInt;
864 sgnInt = totInt-bkgInt;
865 cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
6321ee46 866 if (sgnInt <= 0){
867 cout<<"Setting sgnInt = - sgnInt"<<endl;
868 sgnInt=- sgnInt;
869 }
fabf4d8e 870 /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
56cbeefb 871 TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
872 cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
fabf4d8e 873 funcmass->SetLineColor(4); //blue
874
875 //Set parameters
876 cout<<"\nTOTAL FIT"<<endl;
877
878 if(nFitPars==5){
879 funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
880 funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
56cbeefb 881 //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
6321ee46 882 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
fabf4d8e 883 funcmass->FixParameter(0,totInt);
884 }
885 if (nFitPars==6){
886 funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
887 funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
6321ee46 888 //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
889 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
fabf4d8e 890 funcmass->FixParameter(0,totInt);
891 }
892 if(nFitPars==4){
893 funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
56cbeefb 894 if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
895 else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
fabf4d8e 896 funcmass->FixParameter(0,0.);
897 //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
898 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
899
900 }
6321ee46 901
902 Int_t status;
903
904 status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
905 if (status != 0){
906 cout<<"Minuit returned "<<status<<endl;
907 return kFALSE;
908 }
909
fabf4d8e 910 cout<<"fit done"<<endl;
0566386c 911 //reset value of fMass and fSigmaSgn to those found from fit
912 fMass=funcmass->GetParameter(nFitPars-2);
913 fSigmaSgn=funcmass->GetParameter(nFitPars-1);
fabf4d8e 914
915 for(Int_t i=0;i<nFitPars;i++){
916 fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
917 fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
56cbeefb 918 //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
fabf4d8e 919 }
920 /*
921 //check: cout parameters
922 for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
923 cout<<i<<"\t"<<fFitPars[i]<<endl;
924 }
925 */
74179930 926
a60f573d 927// if(draw){
928// TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
929// TH1F *fhistocopy=new TH1F(*fhistoInvMass);
930// canvas->cd();
931// fhistocopy->DrawClone();
932// if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
933// else funcbkg->DrawClone("sames");
934// funcmass->DrawClone("sames");
935// cout<<"Drawn"<<endl;
936// }
937
938 if(funcmass->GetParameter(nFitPars-1) <0 || funcmass->GetParameter(nFitPars-2) <0 || funcmass->GetParameter(nFitPars-3) <0 ) {
939 cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
b7d4bc49 940 return kFALSE;
941 }
0566386c 942
56cbeefb 943 //increase counter of number of fits done
944 fcounter++;
945
ad73ca8f 946 //contour plots
947 if(draw){
948
949 for (Int_t kpar=1; kpar<nFitPars;kpar++){
950
951 for(Int_t jpar=kpar+1;jpar<nFitPars;jpar++){
952 cout<<"Par "<<kpar<<" and "<<jpar<<endl;
953
954 // produce 2 contours per couple of parameters
955 TGraph* cont[2] = {0x0, 0x0};
956 const Double_t errDef[2] = {1., 4.};
957 for (Int_t i=0; i<2; i++) {
958 gMinuit->SetErrorDef(errDef[i]);
959 cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
6321ee46 960 cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
ad73ca8f 961 }
962
963 if(!cont[0] || !cont[1]){
964 cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
965 continue;
966 }
967
968 // set graph titles and add them to the list
969 TString title = "Contour plot";
970 TString titleX = funcmass->GetParName(kpar);
971 TString titleY = funcmass->GetParName(jpar);
972 for (Int_t i=0; i<2; i++) {
973 cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
974 cont[i]->SetTitle(title);
975 cont[i]->GetXaxis()->SetTitle(titleX);
976 cont[i]->GetYaxis()->SetTitle(titleY);
977 cont[i]->GetYaxis()->SetLabelSize(0.033);
978 cont[i]->GetYaxis()->SetTitleSize(0.033);
979 cont[i]->GetYaxis()->SetTitleOffset(1.67);
980
981 fContourGraph->Add(cont[i]);
982 }
983
984 // plot them
985 TString cvname = Form("c%d%d", kpar, jpar);
986 TCanvas *c4=new TCanvas(cvname,cvname,600,600);
987 c4->cd();
988 cont[1]->SetFillColor(38);
989 cont[1]->Draw("alf");
990 cont[0]->SetFillColor(9);
991 cont[0]->Draw("lf");
992
993 }
994
995 }
996
997 }
998
a60f573d 999 if (ftypeOfFit4Sgn == 1 && funcbkg1) {
1000 delete funcbkg1;
1001 funcbkg1=NULL;
1002 }
1003 if (funcbkg) {
1004 delete funcbkg;
1005 funcbkg=NULL;
1006 }
1007 if (funcmass) {
1008 delete funcmass;
1009 funcmass=NULL;
1010 }
16856d6e 1011
1012 AddFunctionsToHisto();
a60f573d 1013 if (draw) DrawFit();
16856d6e 1014
1015
34c79b83 1016 return kTRUE;
1017}
1018//_________________________________________________________________________
1019Double_t AliHFMassFitter::GetChiSquare() const{
1020 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1021 return funcmass->GetChisquare();
fabf4d8e 1022}
56cbeefb 1023
b7d4bc49 1024//_________________________________________________________________________
1025Double_t AliHFMassFitter::GetReducedChiSquare() const{
1026 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1027 return funcmass->GetChisquare()/funcmass->GetNDF();
1028}
1029
56cbeefb 1030//*********output
1031
fabf4d8e 1032//_________________________________________________________________________
2f328d65 1033void AliHFMassFitter::GetFitPars(Float_t *vector) const {
fabf4d8e 1034 // Return fit parameters
1035
2f328d65 1036 for(Int_t i=0;i<fParsSize;i++){
1037 vector[i]=fFitPars[i];
fabf4d8e 1038 }
1039}
fabf4d8e 1040
fabf4d8e 1041
6321ee46 1042//_________________________________________________________________________
dc262918 1043void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1044
1045 //gives the integral of signal obtained from fit parameters
1046
6321ee46 1047 Int_t index=fParsSize/2 - 3;
1048 valuewitherror[0]=fFitPars[index];
1049 index=fParsSize - 3;
1050 valuewitherror[1]=fFitPars[index];
1051 }
1052
1053
16856d6e 1054//_________________________________________________________________________
1055void AliHFMassFitter::AddFunctionsToHisto(){
fabf4d8e 1056
dc262918 1057 //Add the background function in the complete range to the list of functions attached to the histogram
1058
16856d6e 1059 cout<<"AddFunctionsToHisto called"<<endl;
56cbeefb 1060 TString bkgname = "funcbkg";
ab0fa301 1061 Int_t np=-99;
56cbeefb 1062 switch (ftypeOfFit4Bkg){
16856d6e 1063 case 0: //expo
56cbeefb 1064 np=2;
1065 break;
16856d6e 1066 case 1: //linear
56cbeefb 1067 np=2;
1068 break;
16856d6e 1069 case 2: //pol2
56cbeefb 1070 np=3;
1071 break;
16856d6e 1072 case 3: //no bkg
56cbeefb 1073 np=1;
1074 break;
1075 }
1076 if (ftypeOfFit4Sgn == 1) {
1077 bkgname += 1;
1078 np+=3;
1079 }
16856d6e 1080
1081 TString bkgnamesave=bkgname;
1082 TString testname=bkgname;
1083 testname += "FullRange";
1084 TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1085 if(testfunc){
1086 cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1087 return;
1088 }
1089
1090 TList *hlist=fhistoInvMass->GetListOfFunctions();
34c79b83 1091 hlist->ls();
16856d6e 1092
34c79b83 1093 TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
a60f573d 1094 if(!b){
1095 cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1096 return;
1097 }
1098
56cbeefb 1099 bkgname += "FullRange";
16856d6e 1100 TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
1101 //cout<<bfullrange->GetName()<<endl;
56cbeefb 1102 for(Int_t i=0;i<np;i++){
16856d6e 1103 //cout<<i<<" di "<<np<<endl;
1104 bfullrange->SetParName(i,b->GetParName(i));
1105 bfullrange->SetParameter(i,b->GetParameter(i));
1106 bfullrange->SetParError(i,b->GetParError(i));
1107 }
a60f573d 1108 bfullrange->SetLineStyle(4);
1109 bfullrange->SetLineColor(14);
16856d6e 1110
1111 bkgnamesave += "Recalc";
1112
1113 TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
1114
1115 TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1116
1117 if (!mass){
1118 cout<<"funcmass doesn't exist "<<endl;
1119 return;
1120 }
1121
1122 blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(np));
1123 blastpar->SetParError(0,mass->GetParError(np));
1124 if (np>=2) {
1125 blastpar->SetParameter(1,mass->GetParameter(1));
1126 blastpar->SetParError(1,mass->GetParError(1));
56cbeefb 1127 }
16856d6e 1128 if (np==3) {
1129 blastpar->SetParameter(2,mass->GetParameter(2));
1130 blastpar->SetParError(2,mass->GetParError(2));
1131 }
1132
a60f573d 1133 blastpar->SetLineStyle(1);
16856d6e 1134 blastpar->SetLineColor(2);
1135
1136 hlist->Add(bfullrange->Clone());
1137 hlist->Add(blastpar->Clone());
1138 hlist->ls();
34c79b83 1139
a60f573d 1140 if(bfullrange) {
1141 delete bfullrange;
1142 bfullrange=NULL;
1143 }
1144 if(blastpar) {
1145 delete blastpar;
1146 blastpar=NULL;
1147 }
16856d6e 1148}
1149
1150//_________________________________________________________________________
1151
1152TH1F* AliHFMassFitter::GetHistoClone() const{
1153
1154 TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1155 return hout;
1156}
1157//_________________________________________________________________________
1158
dc262918 1159void AliHFMassFitter::WriteHisto(TString path) const {
1160
1161 //Write the histogram in the default file HFMassFitterOutput.root
1162
16856d6e 1163 if (fcounter == 0) {
1164 cout<<"Use MassFitter method before WriteHisto"<<endl;
1165 return;
1166 }
1167 TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1168
56cbeefb 1169 path += "HFMassFitterOutput.root";
1170 TFile *output;
34c79b83 1171
56cbeefb 1172 if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1173 else output = new TFile(path.Data(),"update");
1174 output->cd();
1175 hget->Write();
ad73ca8f 1176 fContourGraph->Write();
56cbeefb 1177 output->Close();
fabf4d8e 1178
56cbeefb 1179 cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
fabf4d8e 1180
a60f573d 1181 if(output) {
1182 delete output;
1183 output=NULL;
1184 }
fabf4d8e 1185
fabf4d8e 1186}
1187
1188//_________________________________________________________________________
1189
56cbeefb 1190void AliHFMassFitter::WriteNtuple(TString path) const{
a60f573d 1191 //TNtuple* nget=(TNtuple*)fntuParam->Clone();
56cbeefb 1192 path += "HFMassFitterOutput.root";
1193 TFile *output = new TFile(path.Data(),"update");
1194 output->cd();
a60f573d 1195 fntuParam->Write();
1196 //nget->Write();
56cbeefb 1197 output->Close();
a60f573d 1198 //cout<<nget->GetName()<<" written in "<<path<<endl;
1199 cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1200 /*
1201 if(nget) {
1202 //delete nget;
1203 nget=NULL;
1204 }
1205 */
1206 if(output) {
1207 delete output;
1208 output=NULL;
1209 }
1210}
1211
1212//_________________________________________________________________________
1213
1214void AliHFMassFitter::DrawFit() const{
1215
dc262918 1216 //draws histogram together with fit functions with default nice colors
1217
a60f573d 1218 TString cvtitle="fit of ";
1219 cvtitle+=fhistoInvMass->GetName();
1220 TString cvname="c";
1221 cvname+=fcounter;
1222
1223 TCanvas *c=new TCanvas(cvname,cvtitle);
1224 c->cd();
1225 GetHistoClone()->Draw("hist");
1226 GetHistoClone()->GetFunction("funcbkgFullRange")->Draw("same");
1227 GetHistoClone()->GetFunction("funcbkgRecalc")->Draw("same");
1228 GetHistoClone()->GetFunction("funcmass")->Draw("same");
1229
1230
fabf4d8e 1231}
fabf4d8e 1232
fabf4d8e 1233
ad73ca8f 1234//_________________________________________________________________________
1235
1236void AliHFMassFitter::PrintParTitles() const{
dc262918 1237
1238 //prints to screen the parameters names
1239
ad73ca8f 1240 TF1 *f=fhistoInvMass->GetFunction("funcmass");
1241 if(!f) {
1242 cout<<"Fit function not found"<<endl;
1243 return;
1244 }
1245
1246 Int_t np=0;
1247 switch (ftypeOfFit4Bkg){
1248 case 0: //expo
1249 np=2;
1250 break;
1251 case 1: //linear
1252 np=2;
1253 break;
1254 case 2: //pol2
1255 np=3;
1256 break;
1257 case 3: //no bkg
1258 np=1;
1259 break;
1260 }
1261
1262 np+=3; //3 parameter for signal
1263 if (ftypeOfFit4Sgn == 1) np+=3;
1264
1265 cout<<"Parameter Titles \n";
1266 for(Int_t i=0;i<np;i++){
1267 cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1268 }
1269 cout<<endl;
1270
1271}
1272
56cbeefb 1273//************ significance
fabf4d8e 1274
1275//_________________________________________________________________________
1276
1277void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1278 // Return signal integral in mean+- n sigma
1279
16856d6e 1280 if(fcounter==0) {
1281 cout<<"Use MassFitter method before Signal"<<endl;
1282 return;
1283 }
56cbeefb 1284
1285 //functions names
16856d6e 1286 TString bkgname= "funcbkgRecalc";
1287 TString bkg1name="funcbkg1Recalc";
56cbeefb 1288 TString massname="funcmass";
1289
b7d4bc49 1290
1291 TF1 *funcbkg=0;
56cbeefb 1292 TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
16856d6e 1293 if(!funcmass){
1294 cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1295 return;
1296 }
1297
b7d4bc49 1298 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1299 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
b7d4bc49 1300
16856d6e 1301 if(!funcbkg){
1302 cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1303 return;
1304 }
b7d4bc49 1305
16856d6e 1306 Int_t np=-99;
1307 switch (ftypeOfFit4Bkg){
1308 case 0: //expo
1309 np=2;
1310 break;
1311 case 1: //linear
1312 np=2;
1313 break;
1314 case 2: //pol2
1315 np=3;
1316 break;
1317 case 3: //no bkg
1318 np=1;
1319 break;
b7d4bc49 1320 }
1321
fabf4d8e 1322 Double_t intS,intSerr;
1323
b7d4bc49 1324 //relative error evaluation
16856d6e 1325 intS=funcmass->GetParameter(np);
1326 intSerr=funcmass->GetParError(np);
b7d4bc49 1327 cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
16856d6e 1328 Double_t background,errbackground;
1329 Background(nOfSigma,background,errbackground);
1330
b7d4bc49 1331 //signal +/- error in nsigma
1332 Double_t min=fMass-nOfSigma*fSigmaSgn;
1333 Double_t max=fMass+nOfSigma*fSigmaSgn;
16856d6e 1334
1335 Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1336 signal=mass - background;
1337 errsignal=TMath::Sqrt((intSerr/intS*mass)*(intSerr/intS*mass)/*assume relative error is the same as for total integral*/ + errbackground*errbackground);
fabf4d8e 1338 return;
1339}
16856d6e 1340
fabf4d8e 1341//_________________________________________________________________________
1342
16856d6e 1343void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1344
1345 // Return signal integral in a range
1346
1347 if(fcounter==0) {
1348 cout<<"Use MassFitter method before Signal"<<endl;
1349 return;
1350 }
1351
56cbeefb 1352 //functions names
16856d6e 1353 TString bkgname="funcbkgRecalc";
1354 TString bkg1name="funcbkg1Recalc";
56cbeefb 1355 TString massname="funcmass";
fabf4d8e 1356
b7d4bc49 1357
1358 TF1 *funcbkg=0;
56cbeefb 1359 TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
16856d6e 1360 if(!funcmass){
1361 cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1362 return;
1363 }
1364
b7d4bc49 1365 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1366 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1367
16856d6e 1368 if(!funcbkg){
1369 cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1370 return;
1371 }
1372
1373 Int_t np=-99;
1374 switch (ftypeOfFit4Bkg){
1375 case 0: //expo
1376 np=2;
1377 break;
1378 case 1: //linear
1379 np=2;
1380 break;
1381 case 2: //pol2
1382 np=3;
1383 break;
1384 case 3: //no bkg
1385 np=1;
1386 break;
1387 }
b7d4bc49 1388
16856d6e 1389 Double_t intS,intSerr;
b7d4bc49 1390
16856d6e 1391 //relative error evaluation
1392 intS=funcmass->GetParameter(np);
1393 intSerr=funcmass->GetParError(np);
fabf4d8e 1394
16856d6e 1395 cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1396 Double_t background,errbackground;
1397 Background(min,max,background,errbackground);
b7d4bc49 1398
16856d6e 1399 //signal +/- error in the range
1400
1401 Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1402 signal=mass - background;
a60f573d 1403 errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
16856d6e 1404
1405}
1406
1407//_________________________________________________________________________
1408
1409void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1410 // Return background integral in mean+- n sigma
1411
1412 if(fcounter==0) {
1413 cout<<"Use MassFitter method before Background"<<endl;
1414 return;
b7d4bc49 1415 }
1416
16856d6e 1417 //functions names
1418 TString bkgname="funcbkgRecalc";
1419 TString bkg1name="funcbkg1Recalc";
1420
1421 TF1 *funcbkg=0;
1422 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1423 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1424 if(!funcbkg){
1425 cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1426 return;
1427 }
1428
1429 Double_t intB,intBerr;
1430
1431 //relative error evaluation: from final parameters of the fit
1432 if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1433 else{
1434 intB=funcbkg->GetParameter(0);
1435 intBerr=funcbkg->GetParError(0);
1436 cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
fabf4d8e 1437 }
1438
b7d4bc49 1439 Double_t min=fMass-nOfSigma*fSigmaSgn;
1440 Double_t max=fMass+nOfSigma*fSigmaSgn;
1441
1442 //relative error evaluation: from histo
1443
1444 intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1445 Double_t sum2=0;
1446 for(Int_t i=1;i<=fSideBandl;i++){
1447 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1448 }
1449 for(Int_t i=fSideBandr;i<=fNbin;i++){
1450 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1451 }
56cbeefb 1452
b7d4bc49 1453 intBerr=TMath::Sqrt(sum2);
1454 cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1455
16856d6e 1456 cout<<"Last estimation of bkg error is used"<<endl;
b7d4bc49 1457
1458 //backround +/- error in nsigma
34c79b83 1459 if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1460 background = 0;
1461 errbackground = 0;
1462 }
1463 else{
1464 background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
b7d4bc49 1465 errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
16856d6e 1466 //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1467 }
1468 return;
1469
1470}
1471//___________________________________________________________________________
1472
1473void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1474 // Return background integral in a range
1475
1476 if(fcounter==0) {
1477 cout<<"Use MassFitter method before Background"<<endl;
1478 return;
1479 }
1480
1481 //functions names
1482 TString bkgname="funcbkgRecalc";
1483 TString bkg1name="funcbkg1Recalc";
1484
1485 TF1 *funcbkg=0;
1486 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1487 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1488 if(!funcbkg){
1489 cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1490 return;
1491 }
1492
1493
1494 Double_t intB,intBerr;
1495
1496 //relative error evaluation: from final parameters of the fit
1497 if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1498 else{
1499 intB=funcbkg->GetParameter(0);
1500 intBerr=funcbkg->GetParError(0);
1501 cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1502 }
1503
1504 //relative error evaluation: from histo
1505
1506 intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1507 Double_t sum2=0;
1508 for(Int_t i=1;i<=fSideBandl;i++){
1509 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1510 }
1511 for(Int_t i=fSideBandr;i<=fNbin;i++){
1512 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1513 }
1514
1515 intBerr=TMath::Sqrt(sum2);
1516 cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1517
1518 cout<<"Last estimation of bkg error is used"<<endl;
1519
1520 //backround +/- error in the range
1521 if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1522 background = 0;
1523 errbackground = 0;
1524 }
1525 else{
1526 background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1527 errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1528 //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
34c79b83 1529 }
fabf4d8e 1530 return;
1531
1532}
1533
16856d6e 1534
fabf4d8e 1535//__________________________________________________________________________
1536
b7d4bc49 1537void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
fabf4d8e 1538 // Return significance in mean+- n sigma
1539
1540 Double_t signal,errsignal,background,errbackground;
1541 Signal(nOfSigma,signal,errsignal);
1542 Background(nOfSigma,background,errbackground);
1543
1544 significance = signal/TMath::Sqrt(signal+background);
b7d4bc49 1545
fabf4d8e 1546 errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
b7d4bc49 1547
fabf4d8e 1548 return;
1549}
1550
1551//__________________________________________________________________________
74179930 1552
16856d6e 1553void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1554 // Return significance integral in a range
1555
1556 Double_t signal,errsignal,background,errbackground;
1557 Signal(min, max,signal,errsignal);
1558 Background(min, max,background,errbackground);
1559
1560 significance = signal/TMath::Sqrt(signal+background);
1561
1562 errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1563
1564 return;
1565}
ad73ca8f 1566
1567