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