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