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