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