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