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