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