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