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