]> 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
1- Storing in the output of the cuts object
[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=- 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 cout<<"fix1"<<endl;
1097 funcmass->FixParameter(0,totInt);
1098 }
1099 if(fFixPar[1]){
1100 cout<<"fix2"<<endl;
1101 funcmass->FixParameter(1,slope1);
1102 }
1103 if(fFixPar[2]){
1104 cout<<"fix3"<<endl;
1105 funcmass->FixParameter(2,sgnInt);
1106 }
1107 if(fFixPar[3]){
1108 cout<<"fix4"<<endl;
1109 funcmass->FixParameter(3,fMass);
1110 }
1111 if(fFixPar[4]){
1112 cout<<"fix5"<<endl;
1113 funcmass->FixParameter(4,fSigmaSgn);
1114 }
1115 }
1116 if (fNFinalPars==6){
1117 funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
1118 funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
1119
1120 //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1121 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1122 if(fFixPar[0])funcmass->FixParameter(0,totInt);
1123 if(fFixPar[1])funcmass->FixParameter(1,slope1);
1124 if(fFixPar[2])funcmass->FixParameter(2,conc1);
1125 if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
1126 if(fFixPar[4])funcmass->FixParameter(4,fMass);
1127 if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
1128 //
1129 //funcmass->FixParameter(2,sgnInt);
1130 }
1131 if(fNFinalPars==4){
1132 funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
1133 if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
1134 else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
1135 if(fFixPar[0]) funcmass->FixParameter(0,0.);
1136 //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1137 //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1138
1139 }
1140
1141 Int_t status;
1142
1143 status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
1144 if (status != 0){
1145 cout<<"Minuit returned "<<status<<endl;
1146 return kFALSE;
1147 }
1148
1149 cout<<"fit done"<<endl;
1150 //reset value of fMass and fSigmaSgn to those found from fit
1151 fMass=funcmass->GetParameter(fNFinalPars-2);
1152 fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
1153
1154 for(Int_t i=0;i<fNFinalPars;i++){
1155 fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1156 fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1157 //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1158 }
1159 /*
1160 //check: cout parameters
1161 for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
1162 cout<<i<<"\t"<<fFitPars[i]<<endl;
1163 }
1164 */
1165
1166 if(funcmass->GetParameter(fNFinalPars-1) <0 || funcmass->GetParameter(fNFinalPars-2) <0 || funcmass->GetParameter(fNFinalPars-3) <0 ) {
1167 cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1168 return kFALSE;
1169 }
1170
1171 //increase counter of number of fits done
1172 fcounter++;
1173
1174 //contour plots
1175 if(draw){
1176
1177 for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
1178
1179 for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
1180 cout<<"Par "<<kpar<<" and "<<jpar<<endl;
1181
1182 // produce 2 contours per couple of parameters
1183 TGraph* cont[2] = {0x0, 0x0};
1184 const Double_t errDef[2] = {1., 4.};
1185 for (Int_t i=0; i<2; i++) {
1186 gMinuit->SetErrorDef(errDef[i]);
1187 cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
1188 cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
1189 }
1190
1191 if(!cont[0] || !cont[1]){
1192 cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
1193 continue;
1194 }
1195
1196 // set graph titles and add them to the list
1197 TString title = "Contour plot";
1198 TString titleX = funcmass->GetParName(kpar);
1199 TString titleY = funcmass->GetParName(jpar);
1200 for (Int_t i=0; i<2; i++) {
1201 cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
1202 cont[i]->SetTitle(title);
1203 cont[i]->GetXaxis()->SetTitle(titleX);
1204 cont[i]->GetYaxis()->SetTitle(titleY);
1205 cont[i]->GetYaxis()->SetLabelSize(0.033);
1206 cont[i]->GetYaxis()->SetTitleSize(0.033);
1207 cont[i]->GetYaxis()->SetTitleOffset(1.67);
1208
1209 fContourGraph->Add(cont[i]);
1210 }
1211
1212 // plot them
1213 TString cvname = Form("c%d%d", kpar, jpar);
1214 TCanvas *c4=new TCanvas(cvname,cvname,600,600);
1215 c4->cd();
1216 cont[1]->SetFillColor(38);
1217 cont[1]->Draw("alf");
1218 cont[0]->SetFillColor(9);
1219 cont[0]->Draw("lf");
1220
1221 }
1222
1223 }
1224
1225 }
1226
1227 if (ftypeOfFit4Sgn == 1 && funcbkg1) {
1228 delete funcbkg1;
1229 funcbkg1=NULL;
1230 }
1231 if (funcbkg) {
1232 delete funcbkg;
1233 funcbkg=NULL;
1234 }
1235 if (funcmass) {
1236 delete funcmass;
1237 funcmass=NULL;
1238 }
1239
1240 AddFunctionsToHisto();
1241 if (draw) DrawFit();
1242
1243
1244 return kTRUE;
1245}
1246
1247//______________________________________________________________________________
1248
1249Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
1250
1251 //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
1252 //If you want to change the backgroud function or range use SetType or SetRangeFit before
1253
1254 TString bkgname="funcbkgonly";
1255 fSideBands = kFALSE;
1256
1257 TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1258
1259 funcbkg->SetLineColor(kBlue+3); //dark blue
1260
1261 Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
1262
1263 switch (ftypeOfFit4Bkg) {
1264 case 0: //gaus+expo
1265 funcbkg->SetParNames("BkgInt","Slope");
1266 funcbkg->SetParameters(integral,-2.);
1267 break;
1268 case 1:
1269 funcbkg->SetParNames("BkgInt","Slope");
1270 funcbkg->SetParameters(integral,-100.);
1271 break;
1272 case 2:
1273 funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1274 funcbkg->SetParameters(integral,-10.,5);
1275 break;
1276 case 3:
1277 cout<<"Warning! This choice does not have a lot of sense..."<<endl;
1278 if(ftypeOfFit4Sgn==0){
1279 funcbkg->SetParNames("Const");
1280 funcbkg->SetParameter(0,0.);
1281 funcbkg->FixParameter(0,0.);
1282 }
1283 break;
1284 default:
1285 cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1286 return kFALSE;
1287 break;
1288 }
1289
1290
1291 Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
1292 if (status != 0){
1293 cout<<"Minuit returned "<<status<<endl;
1294 return kFALSE;
1295 }
1296 AddFunctionsToHisto();
1297
1298 if(draw) DrawFit();
1299
1300 return kTRUE;
1301
1302}
1303//_________________________________________________________________________
1304Double_t AliHFMassFitter::GetChiSquare() const{
1305 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1306 if(!funcmass) {
1307 cout<<"funcmass not found"<<endl;
1308 return -1;
1309 }
1310 return funcmass->GetChisquare();
1311}
1312
1313//_________________________________________________________________________
1314Double_t AliHFMassFitter::GetReducedChiSquare() const{
1315 TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1316 if(!funcmass) {
1317 cout<<"funcmass not found"<<endl;
1318 return -1;
1319 }
1320 return funcmass->GetChisquare()/funcmass->GetNDF();
1321}
1322
1323//*********output
1324
1325//_________________________________________________________________________
1326void AliHFMassFitter::GetFitPars(Float_t *vector) const {
1327 // Return fit parameters
1328
1329 for(Int_t i=0;i<fParsSize;i++){
1330 vector[i]=fFitPars[i];
1331 }
1332}
1333
1334
1335//_________________________________________________________________________
1336void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1337
1338 //gives the integral of signal obtained from fit parameters
1339 if(!valuewitherror)valuewitherror=new Float_t[2];
1340
1341 Int_t index=fParsSize/2 - 3;
1342 valuewitherror[0]=fFitPars[index];
1343 index=fParsSize - 3;
1344 valuewitherror[1]=fFitPars[index];
1345 }
1346
1347
1348//_________________________________________________________________________
1349void AliHFMassFitter::AddFunctionsToHisto(){
1350
1351 //Add the background function in the complete range to the list of functions attached to the histogram
1352
1353 cout<<"AddFunctionsToHisto called"<<endl;
1354 TString bkgname = "funcbkg";
1355
1356 Bool_t done1=kFALSE,done2=kFALSE;
1357
1358 TString bkgnamesave=bkgname;
1359 TString testname=bkgname;
1360 testname += "FullRange";
1361 TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1362 if(testfunc){
1363 done1=kTRUE;
1364 testfunc=0x0;
1365 }
1366 testname="funcbkgonly";
1367 testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1368 if(testfunc){
1369 done2=kTRUE;
1370 testfunc=0x0;
1371 }
1372
1373 if(done1 && done2){
1374 cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1375 return;
1376 }
1377
1378 TList *hlist=fhistoInvMass->GetListOfFunctions();
1379 hlist->ls();
1380
1381 if(!done2){
1382 TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1383 if(!bonly){
1384 cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1385 }else{
1386 bonly->SetLineColor(kBlue+3);
1387 hlist->Add((TF1*)bonly->Clone());
1388 if(bonly) {
1389 delete bonly;
1390 bonly=NULL;
1391 }
1392 }
1393
1394 }
1395
1396 if(!done1){
1397 TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1398 if(!b){
1399 cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1400 return;
1401 }
1402
1403 bkgname += "FullRange";
1404 TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4Bkg");
1405 //cout<<bfullrange->GetName()<<endl;
1406 for(Int_t i=0;i<fNFinalPars;i++){
1407 //cout<<i<<" di "<<fNFinalPars<<endl;
1408 bfullrange->SetParName(i,b->GetParName(i));
1409 bfullrange->SetParameter(i,b->GetParameter(i));
1410 bfullrange->SetParError(i,b->GetParError(i));
1411 }
1412 bfullrange->SetLineStyle(4);
1413 bfullrange->SetLineColor(14);
1414
1415 bkgnamesave += "Recalc";
1416
1417 TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4Bkg");
1418
1419 TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1420
1421 if (!mass){
1422 cout<<"funcmass doesn't exist "<<endl;
1423 return;
1424 }
1425
1426 //intBkg=intTot-intS
1427 blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
1428 blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
1429 if (fNFinalPars>=2) {
1430 blastpar->SetParameter(1,mass->GetParameter(1));
1431 blastpar->SetParError(1,mass->GetParError(1));
1432 }
1433 if (fNFinalPars==3) {
1434 blastpar->SetParameter(2,mass->GetParameter(2));
1435 blastpar->SetParError(2,mass->GetParError(2));
1436 }
1437
1438 blastpar->SetLineStyle(1);
1439 blastpar->SetLineColor(2);
1440
1441 hlist->Add((TF1*)bfullrange->Clone());
1442 hlist->Add((TF1*)blastpar->Clone());
1443 hlist->ls();
1444
1445 if(bfullrange) {
1446 delete bfullrange;
1447 bfullrange=NULL;
1448 }
1449 if(blastpar) {
1450 delete blastpar;
1451 blastpar=NULL;
1452 }
1453 }
1454
1455
1456}
1457
1458//_________________________________________________________________________
1459
1460TH1F* AliHFMassFitter::GetHistoClone() const{
1461
1462 TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1463 return hout;
1464}
1465//_________________________________________________________________________
1466
1467void AliHFMassFitter::WriteHisto(TString path) const {
1468
1469 //Write the histogram in the default file HFMassFitterOutput.root
1470
1471 if (fcounter == 0) {
1472 cout<<"Use MassFitter method before WriteHisto"<<endl;
1473 return;
1474 }
1475 TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1476
1477 path += "HFMassFitterOutput.root";
1478 TFile *output;
1479
1480 if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1481 else output = new TFile(path.Data(),"update");
1482 output->cd();
1483 hget->Write();
1484 fContourGraph->Write();
1485
1486
1487 output->Close();
1488
1489 cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1490
1491 if(output) {
1492 delete output;
1493 output=NULL;
1494 }
1495
1496}
1497
1498//_________________________________________________________________________
1499
1500void AliHFMassFitter::WriteNtuple(TString path) const{
1501 //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1502 path += "HFMassFitterOutput.root";
1503 TFile *output = new TFile(path.Data(),"update");
1504 output->cd();
1505 fntuParam->Write();
1506 //nget->Write();
1507 output->Close();
1508 //cout<<nget->GetName()<<" written in "<<path<<endl;
1509 cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1510 /*
1511 if(nget) {
1512 //delete nget;
1513 nget=NULL;
1514 }
1515 */
1516 if(output) {
1517 delete output;
1518 output=NULL;
1519 }
1520}
1521
1522//_________________________________________________________________________
1523void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1524
1525 //write the canvas in a root file
1526
1527 gStyle->SetOptStat(0);
1528 gStyle->SetCanvasColor(0);
1529 gStyle->SetFrameFillColor(0);
1530
1531 TString type="";
1532
1533 switch (ftypeOfFit4Bkg){
1534 case 0:
1535 type="Exp"; //3+2
1536 break;
1537 case 1:
1538 type="Lin"; //3+2
1539 break;
1540 case 2:
1541 type="Pl2"; //3+3
1542 break;
1543 case 3:
1544 type="noB"; //3+1
1545 break;
1546 }
1547
1548 TString filename=Form("%sMassFit.root",type.Data());
1549 filename.Prepend(userIDstring);
1550 path.Append(filename);
1551
1552 TFile* outputcv=new TFile(path.Data(),"update");
1553
1554 TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1555 c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1556 if(draw)c->DrawClone();
1557 outputcv->cd();
1558 c->Write();
1559 outputcv->Close();
1560}
1561
1562//_________________________________________________________________________
1563
1564TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1565 //return a TVirtualPad with the fitted histograms and info
1566
1567 TString cvtitle="fit of ";
1568 cvtitle+=fhistoInvMass->GetName();
1569 TString cvname="c";
1570 cvname+=fcounter;
1571
1572 TCanvas *c=new TCanvas(cvname,cvtitle);
1573 PlotFit(c->cd(),nsigma,writeFitInfo);
1574 return c->cd();
1575}
1576//_________________________________________________________________________
1577
1578void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1579 //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1580 gStyle->SetOptStat(0);
1581 gStyle->SetCanvasColor(0);
1582 gStyle->SetFrameFillColor(0);
1583
1584 cout<<"nsigma = "<<nsigma<<endl;
1585 cout<<"Verbosity = "<<writeFitInfo<<endl;
1586
1587 TH1F* hdraw=GetHistoClone();
1588
1589 if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1590 cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1591 return;
1592 }
1593
1594 if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1595 cout<<"Drawing background fit only"<<endl;
1596 hdraw->SetMinimum(0);
1597 hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1598 pd->cd();
1599 hdraw->SetMarkerStyle(20);
1600 hdraw->DrawClone("PE");
1601 hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1602
1603 if(writeFitInfo > 0){
1604 TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1605 pinfo->SetBorderSize(0);
1606 pinfo->SetFillStyle(0);
1607 TF1* f=hdraw->GetFunction("funcbkgonly");
1608 for (Int_t i=0;i<fNFinalPars-3;i++){
1609 pinfo->SetTextColor(kBlue+3);
1610 TString str=Form("%s = %f #pm %f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1611 pinfo->AddText(str);
1612 }
1613
1614 pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1615 pd->cd();
1616 pinfo->DrawClone();
1617
1618
1619 }
1620
1621 return;
1622 }
1623
1624 hdraw->SetMinimum(0);
1625 hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1626 pd->cd();
1627 hdraw->SetMarkerStyle(20);
1628 hdraw->DrawClone("PE");
1629 if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1630 if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1631 if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1632
1633 if(writeFitInfo > 0){
1634 TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1635 TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1636 pinfob->SetBorderSize(0);
1637 pinfob->SetFillStyle(0);
1638 pinfom->SetBorderSize(0);
1639 pinfom->SetFillStyle(0);
1640 TF1* ff=fhistoInvMass->GetFunction("funcmass");
1641
1642 for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1643 pinfom->SetTextColor(kBlue);
1644 TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1645 if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1646 }
1647 pd->cd();
1648 pinfom->DrawClone();
1649
1650 TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1651 pinfo2->SetBorderSize(0);
1652 pinfo2->SetFillStyle(0);
1653
1654 Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1655
1656 Significance(nsigma,signif,errsignif);
1657 Signal(nsigma,signal,errsignal);
1658 Background(nsigma,bkg, errbkg);
1659 /*
1660 Significance(1.828,1.892,signif,errsignif);
1661 Signal(1.828,1.892,signal,errsignal);
1662 Background(1.828,1.892,bkg, errbkg);
1663 */
1664 TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1665 pinfo2->AddText(str);
1666 str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1667 pinfo2->AddText(str);
1668 str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1669 pinfo2->AddText(str);
1670
1671 pd->cd();
1672 pinfo2->Draw();
1673
1674 if(writeFitInfo > 1){
1675 for (Int_t i=0;i<fNFinalPars-3;i++){
1676 pinfob->SetTextColor(kRed);
1677 str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1678 pinfob->AddText(str);
1679 }
1680 pd->cd();
1681 pinfob->DrawClone();
1682 }
1683
1684
1685 }
1686 return;
1687}
1688
1689//_________________________________________________________________________
1690
1691void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1692 //draws histogram together with fit functions with default nice colors in user canvas
1693 PlotFit(pd,nsigma,writeFitInfo);
1694
1695 pd->Draw();
1696
1697}
1698//_________________________________________________________________________
1699void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1700
1701 //draws histogram together with fit functions with default nice colors
1702 gStyle->SetOptStat(0);
1703 gStyle->SetCanvasColor(0);
1704 gStyle->SetFrameFillColor(0);
1705
1706
1707 TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1708 c->Draw();
1709
1710}
1711
1712
1713//_________________________________________________________________________
1714
1715void AliHFMassFitter::PrintParTitles() const{
1716
1717 //prints to screen the parameters names
1718
1719 TF1 *f=fhistoInvMass->GetFunction("funcmass");
1720 if(!f) {
1721 cout<<"Fit function not found"<<endl;
1722 return;
1723 }
1724
1725 cout<<"Parameter Titles \n";
1726 for(Int_t i=0;i<fNFinalPars;i++){
1727 cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1728 }
1729 cout<<endl;
1730
1731}
1732
1733//************ significance
1734
1735//_________________________________________________________________________
1736
1737void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1738 // Return signal integral in mean+- n sigma
1739
1740 if(fcounter==0) {
1741 cout<<"Use MassFitter method before Signal"<<endl;
1742 return;
1743 }
1744
1745 Double_t min=fMass-nOfSigma*fSigmaSgn;
1746 Double_t max=fMass+nOfSigma*fSigmaSgn;
1747
1748 Signal(min,max,signal,errsignal);
1749
1750
1751 return;
1752
1753}
1754
1755//_________________________________________________________________________
1756
1757void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1758
1759 // Return signal integral in a range
1760
1761 if(fcounter==0) {
1762 cout<<"Use MassFitter method before Signal"<<endl;
1763 return;
1764 }
1765
1766 //functions names
1767 TString bkgname="funcbkgRecalc";
1768 TString bkg1name="funcbkg1Recalc";
1769 TString massname="funcmass";
1770
1771
1772 TF1 *funcbkg=0;
1773 TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1774 if(!funcmass){
1775 cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1776 return;
1777 }
1778
1779 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1780 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1781
1782 if(!funcbkg){
1783 cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1784 return;
1785 }
1786
1787 Int_t np=fNFinalPars-3;
1788
1789 Double_t intS,intSerr;
1790
1791 //relative error evaluation
1792 intS=funcmass->GetParameter(np);
1793 intSerr=funcmass->GetParError(np);
1794
1795 cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1796 Double_t background,errbackground;
1797 Background(min,max,background,errbackground);
1798
1799 //signal +/- error in the range
1800
1801 Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1802 signal=mass - background;
1803 errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1804
1805}
1806
1807//_________________________________________________________________________
1808
1809void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1810 // Return background integral in mean+- n sigma
1811
1812 if(fcounter==0) {
1813 cout<<"Use MassFitter method before Background"<<endl;
1814 return;
1815 }
1816 Double_t min=fMass-nOfSigma*fSigmaSgn;
1817 Double_t max=fMass+nOfSigma*fSigmaSgn;
1818
1819 Background(min,max,background,errbackground);
1820
1821 return;
1822
1823}
1824//___________________________________________________________________________
1825
1826void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1827 // Return background integral in a range
1828
1829 if(fcounter==0) {
1830 cout<<"Use MassFitter method before Background"<<endl;
1831 return;
1832 }
1833
1834 //functions names
1835 TString bkgname="funcbkgRecalc";
1836 TString bkg1name="funcbkg1Recalc";
1837
1838 TF1 *funcbkg=0;
1839 if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1840 else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1841 if(!funcbkg){
1842 cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1843 return;
1844 }
1845
1846
1847 Double_t intB,intBerr;
1848
1849 //relative error evaluation: from final parameters of the fit
1850 if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1851 else{
1852 intB=funcbkg->GetParameter(0);
1853 intBerr=funcbkg->GetParError(0);
1854 cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1855 }
1856
1857 //relative error evaluation: from histo
1858
1859 intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1860 Double_t sum2=0;
1861 for(Int_t i=1;i<=fSideBandl;i++){
1862 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1863 }
1864 for(Int_t i=fSideBandr;i<=fNbin;i++){
1865 sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1866 }
1867
1868 intBerr=TMath::Sqrt(sum2);
1869 cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1870
1871 cout<<"Last estimation of bkg error is used"<<endl;
1872
1873 //backround +/- error in the range
1874 if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1875 background = 0;
1876 errbackground = 0;
1877 }
1878 else{
1879 background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1880 errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1881 //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1882 }
1883 return;
1884
1885}
1886
1887
1888//__________________________________________________________________________
1889
1890void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
1891 // Return significance in mean+- n sigma
1892
1893 Double_t min=fMass-nOfSigma*fSigmaSgn;
1894 Double_t max=fMass+nOfSigma*fSigmaSgn;
1895 Significance(min, max, significance, errsignificance);
1896
1897 return;
1898}
1899
1900//__________________________________________________________________________
1901
1902void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1903 // Return significance integral in a range
1904
1905 Double_t signal,errsignal,background,errbackground;
1906 Signal(min, max,signal,errsignal);
1907 Background(min, max,background,errbackground);
1908
1909 if (signal+background <= 0.){
1910 cout<<"Cannot calculate significance because of div by 0!"<<endl;
1911 significance=-1;
1912 errsignificance=0;
1913 }
1914
1915 significance = signal/TMath::Sqrt(signal+background);
1916
1917 errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1918
1919 return;
1920}
1921
1922