]>
Commit | Line | Data |
---|---|---|
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 | ||
24 | #include <TCanvas.h> | |
ad73ca8f | 25 | #include <TMinuit.h> |
26 | #include <TGraph.h> | |
27 | #include <TVirtualFitter.h> | |
28 | #include <TClonesArray.h> | |
fabf4d8e | 29 | |
30 | #include "AliHFMassFitter.h" | |
31 | ||
32 | ||
33 | ClassImp(AliHFMassFitter) | |
34 | ||
56cbeefb | 35 | //************** constructors |
2f328d65 | 36 | AliHFMassFitter::AliHFMassFitter() : |
37 | TNamed(), | |
fabf4d8e | 38 | fhistoInvMass(0), |
39 | fminMass(0), | |
40 | fmaxMass(0), | |
34c79b83 | 41 | fNbin(0), |
2f328d65 | 42 | fParsSize(1), |
43 | fFitPars(0), | |
fabf4d8e | 44 | fWithBkg(0), |
45 | ftypeOfFit4Bkg(0), | |
46 | ftypeOfFit4Sgn(0), | |
47 | ffactor(1), | |
48 | fntuParam(0), | |
49 | fMass(1.85), | |
50 | fSigmaSgn(0.012), | |
56cbeefb | 51 | fSideBands(0), |
b7d4bc49 | 52 | fSideBandl(0), |
53 | fSideBandr(0), | |
ad73ca8f | 54 | fcounter(0), |
55 | fContourGraph(0) | |
fabf4d8e | 56 | { |
57 | // default constructor | |
58 | ||
59 | cout<<"Default constructor"<<endl; | |
ad73ca8f | 60 | |
fabf4d8e | 61 | } |
62 | ||
63 | //___________________________________________________________________________ | |
64 | ||
65 | AliHFMassFitter::AliHFMassFitter (TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes): | |
66 | TNamed(), | |
67 | fhistoInvMass(0), | |
68 | fminMass(0), | |
69 | fmaxMass(0), | |
34c79b83 | 70 | fNbin(0), |
2f328d65 | 71 | fParsSize(1), |
72 | fFitPars(0), | |
fabf4d8e | 73 | fWithBkg(0), |
74 | ftypeOfFit4Bkg(0), | |
75 | ftypeOfFit4Sgn(0), | |
76 | ffactor(1), | |
77 | fntuParam(0), | |
78 | fMass(1.85), | |
79 | fSigmaSgn(0.012), | |
56cbeefb | 80 | fSideBands(0), |
b7d4bc49 | 81 | fSideBandl(0), |
82 | fSideBandr(0), | |
ad73ca8f | 83 | fcounter(0), |
84 | fContourGraph(0) | |
fabf4d8e | 85 | { |
86 | // standard constructor | |
87 | ||
56cbeefb | 88 | fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass"); |
6321ee46 | 89 | fhistoInvMass->SetDirectory(0); |
fabf4d8e | 90 | fminMass=minvalue; |
91 | fmaxMass=maxvalue; | |
92 | if(rebin!=1) RebinMass(rebin); | |
93 | else fNbin=(Int_t)fhistoInvMass->GetNbinsX(); | |
94 | ftypeOfFit4Bkg=fittypeb; | |
95 | ftypeOfFit4Sgn=fittypes; | |
96 | if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE; | |
97 | else fWithBkg=kTRUE; | |
98 | if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl; | |
99 | else cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl; | |
2f328d65 | 100 | |
101 | ComputeParSize(); | |
102 | fFitPars=new Float_t[fParsSize]; | |
ad73ca8f | 103 | |
104 | fContourGraph=new TList(); | |
105 | fContourGraph->SetOwner(); | |
106 | ||
fabf4d8e | 107 | } |
108 | ||
74179930 | 109 | |
110 | AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit): | |
111 | TNamed(), | |
112 | fhistoInvMass(mfit.fhistoInvMass), | |
113 | fminMass(mfit.fminMass), | |
114 | fmaxMass(mfit.fmaxMass), | |
115 | fNbin(mfit.fNbin), | |
2f328d65 | 116 | fParsSize(mfit.fParsSize), |
117 | fFitPars(0), | |
74179930 | 118 | fWithBkg(mfit.fWithBkg), |
119 | ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg), | |
120 | ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn), | |
121 | ffactor(mfit.ffactor), | |
122 | fntuParam(mfit.fntuParam), | |
123 | fMass(mfit.fMass), | |
124 | fSigmaSgn(mfit.fSigmaSgn), | |
56cbeefb | 125 | fSideBands(mfit.fSideBands), |
b7d4bc49 | 126 | fSideBandl(mfit.fSideBandl), |
127 | fSideBandr(mfit.fSideBandr), | |
ad73ca8f | 128 | fcounter(mfit.fcounter), |
129 | fContourGraph(mfit.fContourGraph) | |
74179930 | 130 | { |
131 | //copy constructor | |
2f328d65 | 132 | |
133 | if(mfit.fParsSize > 0){ | |
134 | fFitPars=new Float_t[fParsSize]; | |
135 | memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t)); | |
136 | } | |
137 | //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i]; | |
138 | } | |
139 | ||
140 | //_________________________________________________________________________ | |
141 | ||
142 | AliHFMassFitter::~AliHFMassFitter() { | |
56cbeefb | 143 | cout<<"AliHFMassFitter destructor called"<<endl; |
a60f573d | 144 | if(fhistoInvMass) { |
145 | cout<<"deleting histogram..."<<endl; | |
146 | delete fhistoInvMass; | |
147 | fhistoInvMass=NULL; | |
148 | } | |
149 | if(fntuParam){ | |
150 | cout<<"deleting ntuple..."<<endl; | |
151 | delete fntuParam; | |
152 | ||
153 | fntuParam=NULL; | |
154 | } | |
155 | ||
156 | if(fFitPars) { | |
157 | delete[] fFitPars; | |
158 | cout<<"deleting parameter array..."<<endl; | |
159 | fFitPars=NULL; | |
160 | } | |
161 | ||
56cbeefb | 162 | fcounter = 0; |
74179930 | 163 | } |
164 | ||
165 | //_________________________________________________________________________ | |
166 | ||
167 | AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){ | |
168 | ||
169 | //assignment operator | |
170 | ||
171 | if(&mfit == this) return *this; | |
172 | fhistoInvMass= mfit.fhistoInvMass; | |
173 | fminMass= mfit.fminMass; | |
174 | fmaxMass= mfit.fmaxMass; | |
175 | fNbin= mfit.fNbin; | |
2f328d65 | 176 | fParsSize= mfit.fParsSize; |
74179930 | 177 | fWithBkg= mfit.fWithBkg; |
178 | ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg; | |
179 | ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn; | |
180 | ffactor= mfit.ffactor; | |
181 | fntuParam= mfit.fntuParam; | |
182 | fMass= mfit.fMass; | |
183 | fSigmaSgn= mfit.fSigmaSgn; | |
184 | fSideBands= mfit.fSideBands; | |
b7d4bc49 | 185 | fSideBandl= mfit.fSideBandl; |
186 | fSideBandr= mfit.fSideBandr; | |
56cbeefb | 187 | fcounter= mfit.fcounter; |
ad73ca8f | 188 | fContourGraph= mfit.fContourGraph; |
189 | ||
2f328d65 | 190 | if(mfit.fParsSize > 0){ |
a60f573d | 191 | if(fFitPars) { |
192 | delete[] fFitPars; | |
193 | fFitPars=NULL; | |
194 | } | |
2f328d65 | 195 | fFitPars=new Float_t[fParsSize]; |
196 | memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t)); | |
197 | } | |
198 | // fFitPars=new Float_t[fParsSize]; | |
199 | // for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i]; | |
74179930 | 200 | |
201 | return *this; | |
202 | } | |
56cbeefb | 203 | |
204 | //************ tools & settings | |
205 | ||
2f328d65 | 206 | //__________________________________________________________________________ |
207 | ||
208 | void AliHFMassFitter::ComputeParSize() { | |
209 | switch (ftypeOfFit4Bkg) {//npar backround func | |
210 | case 0: | |
211 | fParsSize = 2*3; | |
212 | break; | |
213 | case 1: | |
214 | fParsSize = 2*3; | |
215 | break; | |
216 | case 2: | |
217 | fParsSize = 3*3; | |
218 | break; | |
219 | case 3: | |
56cbeefb | 220 | fParsSize = 1*3; |
2f328d65 | 221 | break; |
222 | default: | |
56cbeefb | 223 | cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl; |
2f328d65 | 224 | break; |
225 | } | |
226 | ||
56cbeefb | 227 | fParsSize += 3; // npar refl |
228 | fParsSize += 3; // npar signal gaus | |
2f328d65 | 229 | |
56cbeefb | 230 | fParsSize*=2; // add errors |
2f328d65 | 231 | cout<<"Parameters array size "<<fParsSize<<endl; |
232 | } | |
233 | ||
56cbeefb | 234 | //___________________________________________________________________________ |
235 | void AliHFMassFitter::SetHisto(TH1F *histoToFit){ | |
a60f573d | 236 | //fhistoInvMass = (TH1F*)histoToFit->Clone(); |
237 | fhistoInvMass = new TH1F(*histoToFit); | |
6321ee46 | 238 | fhistoInvMass->SetDirectory(0); |
a60f573d | 239 | cout<<"SetHisto pointer "<<fhistoInvMass<<endl; |
56cbeefb | 240 | } |
241 | ||
2f328d65 | 242 | //___________________________________________________________________________ |
243 | ||
244 | void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) { | |
245 | ftypeOfFit4Bkg = fittypeb; | |
246 | ftypeOfFit4Sgn = fittypes; | |
247 | ComputeParSize(); | |
248 | fFitPars = new Float_t[fParsSize]; | |
249 | ||
250 | } | |
74179930 | 251 | |
fabf4d8e | 252 | //___________________________________________________________________________ |
253 | ||
56cbeefb | 254 | void AliHFMassFitter::Reset() { |
a60f573d | 255 | cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl; |
0566386c | 256 | fMass=1.85; |
257 | fSigmaSgn=0.012; | |
a60f573d | 258 | cout<<"Reset "<<fhistoInvMass<<endl; |
259 | if(fhistoInvMass) { | |
ad73ca8f | 260 | //cout<<"esiste"<<endl; |
6321ee46 | 261 | delete fhistoInvMass; |
a60f573d | 262 | fhistoInvMass=NULL; |
263 | cout<<fhistoInvMass<<endl; | |
264 | } | |
265 | else cout<<"histogram doesn't exist, do not delete"<<endl; | |
266 | ||
56cbeefb | 267 | } |
268 | ||
269 | //_________________________________________________________________________ | |
270 | ||
271 | void AliHFMassFitter::InitNtuParam(char *ntuname) { | |
272 | // Create ntuple to keep fit parameters | |
273 | fntuParam=0; | |
274 | 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"); | |
275 | ||
276 | } | |
277 | ||
278 | //_________________________________________________________________________ | |
279 | ||
280 | void AliHFMassFitter::FillNtuParam() { | |
281 | // Fill ntuple with fit parameters | |
282 | ||
283 | Float_t nothing=0.; | |
284 | ||
285 | if (ftypeOfFit4Bkg==2) { | |
286 | fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]); | |
287 | fntuParam->SetBranchAddress("slope1",&fFitPars[1]); | |
288 | fntuParam->SetBranchAddress("conc1",&fFitPars[2]); | |
289 | fntuParam->SetBranchAddress("intGB",&fFitPars[3]); | |
290 | fntuParam->SetBranchAddress("meanGB",&fFitPars[4]); | |
291 | fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]); | |
292 | fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]); | |
293 | fntuParam->SetBranchAddress("slope2",&fFitPars[7]); | |
294 | fntuParam->SetBranchAddress("conc2",&fFitPars[8]); | |
295 | fntuParam->SetBranchAddress("inttot",&fFitPars[9]); | |
296 | fntuParam->SetBranchAddress("slope3",&fFitPars[10]); | |
297 | fntuParam->SetBranchAddress("conc3",&fFitPars[11]); | |
298 | fntuParam->SetBranchAddress("intsgn",&fFitPars[12]); | |
299 | fntuParam->SetBranchAddress("meansgn",&fFitPars[13]); | |
300 | fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]); | |
301 | ||
302 | fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]); | |
303 | fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]); | |
304 | fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]); | |
305 | fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]); | |
306 | fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]); | |
307 | fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]); | |
308 | fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]); | |
309 | fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]); | |
310 | fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]); | |
311 | fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]); | |
312 | fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]); | |
313 | fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]); | |
314 | fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]); | |
315 | fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]); | |
316 | fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]); | |
317 | ||
318 | } else { | |
319 | ||
320 | if(ftypeOfFit4Bkg==3){ | |
321 | fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]); | |
322 | fntuParam->SetBranchAddress("slope1",¬hing); | |
323 | fntuParam->SetBranchAddress("conc1",¬hing); | |
324 | fntuParam->SetBranchAddress("intGB",&fFitPars[1]); | |
325 | fntuParam->SetBranchAddress("meanGB",&fFitPars[2]); | |
326 | fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]); | |
327 | fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]); | |
328 | fntuParam->SetBranchAddress("slope2",¬hing); | |
329 | fntuParam->SetBranchAddress("conc2",¬hing); | |
330 | fntuParam->SetBranchAddress("inttot",&fFitPars[6]); | |
331 | fntuParam->SetBranchAddress("slope3",¬hing); | |
332 | fntuParam->SetBranchAddress("conc3",¬hing); | |
333 | fntuParam->SetBranchAddress("intsgn",&fFitPars[6]); | |
334 | fntuParam->SetBranchAddress("meansgn",&fFitPars[7]); | |
335 | fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]); | |
336 | ||
337 | fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]); | |
338 | fntuParam->SetBranchAddress("slope1Err",¬hing); | |
339 | fntuParam->SetBranchAddress("conc1Err",¬hing); | |
340 | fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]); | |
341 | fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]); | |
342 | fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]); | |
343 | fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]); | |
344 | fntuParam->SetBranchAddress("slope2Err",¬hing); | |
345 | fntuParam->SetBranchAddress("conc2Err",¬hing); | |
346 | fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]); | |
347 | fntuParam->SetBranchAddress("slope3Err",¬hing); | |
348 | fntuParam->SetBranchAddress("conc3Err",¬hing); | |
349 | fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]); | |
350 | fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]); | |
351 | fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]); | |
352 | ||
353 | } | |
354 | else{ | |
355 | fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]); | |
356 | fntuParam->SetBranchAddress("slope1",&fFitPars[1]); | |
357 | fntuParam->SetBranchAddress("conc1",¬hing); | |
358 | fntuParam->SetBranchAddress("intGB",&fFitPars[2]); | |
359 | fntuParam->SetBranchAddress("meanGB",&fFitPars[3]); | |
360 | fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]); | |
361 | fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]); | |
362 | fntuParam->SetBranchAddress("slope2",&fFitPars[6]); | |
363 | fntuParam->SetBranchAddress("conc2",¬hing); | |
364 | fntuParam->SetBranchAddress("inttot",&fFitPars[7]); | |
365 | fntuParam->SetBranchAddress("slope3",&fFitPars[8]); | |
366 | fntuParam->SetBranchAddress("conc3",¬hing); | |
367 | fntuParam->SetBranchAddress("intsgn",&fFitPars[9]); | |
368 | fntuParam->SetBranchAddress("meansgn",&fFitPars[10]); | |
369 | fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]); | |
370 | ||
371 | fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]); | |
372 | fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]); | |
373 | fntuParam->SetBranchAddress("conc1Err",¬hing); | |
374 | fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]); | |
375 | fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]); | |
376 | fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]); | |
377 | fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]); | |
378 | fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]); | |
379 | fntuParam->SetBranchAddress("conc2Err",¬hing); | |
380 | fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]); | |
381 | fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]); | |
382 | fntuParam->SetBranchAddress("conc3Err",¬hing); | |
383 | fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]); | |
384 | fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]); | |
385 | fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]); | |
386 | } | |
387 | ||
388 | } | |
389 | fntuParam->TTree::Fill(); | |
390 | } | |
391 | ||
392 | //_________________________________________________________________________ | |
393 | ||
394 | TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){ | |
395 | // Create, fill and return ntuple with fit parameters | |
396 | ||
397 | InitNtuParam(ntuname); | |
398 | FillNtuParam(); | |
399 | return fntuParam; | |
400 | } | |
401 | //_________________________________________________________________________ | |
402 | ||
403 | void AliHFMassFitter::RebinMass(Int_t bingroup){ | |
404 | // Rebin invariant mass histogram | |
405 | ||
406 | if(bingroup<1){ | |
407 | cout<<"Error! Cannot group "<<bingroup<<" bins\n"; | |
408 | fNbin=fhistoInvMass->GetNbinsX(); | |
409 | cout<<"Kept original number of bins: "<<fNbin<<endl; | |
410 | } else{ | |
411 | fhistoInvMass->Rebin(bingroup); | |
412 | fNbin = fhistoInvMass->GetNbinsX(); | |
413 | cout<<"New number of bins: "<<fNbin<<endl; | |
414 | } | |
415 | ||
416 | ||
417 | } | |
418 | ||
419 | //************* fit | |
420 | ||
421 | //___________________________________________________________________________ | |
422 | ||
fabf4d8e | 423 | Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){ |
424 | // Fit function for signal+background | |
425 | ||
426 | ||
427 | //exponential or linear fit | |
428 | // | |
429 | // par[0] = tot integral | |
430 | // par[1] = slope | |
431 | // par[2] = gaussian integral | |
432 | // par[3] = gaussian mean | |
433 | // par[4] = gaussian sigma | |
434 | ||
435 | Double_t total,bkg=0,sgn=0; | |
436 | ||
437 | if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) { | |
438 | if(ftypeOfFit4Sgn == 0) { | |
439 | ||
440 | Double_t parbkg[2] = {par[0]-par[2], par[1]}; | |
441 | bkg = FitFunction4Bkg(x,parbkg); | |
442 | } | |
443 | if(ftypeOfFit4Sgn == 1) { | |
444 | Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]}; | |
445 | bkg = FitFunction4Bkg(x,parbkg); | |
446 | } | |
447 | ||
448 | sgn = FitFunction4Sgn(x,&par[2]); | |
449 | ||
450 | } | |
451 | ||
452 | //polynomial fit | |
453 | ||
454 | // par[0] = tot integral | |
455 | // par[1] = coef1 | |
456 | // par[2] = coef2 | |
457 | // par[3] = gaussian integral | |
458 | // par[4] = gaussian mean | |
459 | // par[5] = gaussian sigma | |
460 | ||
461 | if (ftypeOfFit4Bkg==2) { | |
462 | ||
463 | if(ftypeOfFit4Sgn == 0) { | |
56cbeefb | 464 | |
fabf4d8e | 465 | Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]}; |
466 | bkg = FitFunction4Bkg(x,parbkg); | |
467 | } | |
468 | if(ftypeOfFit4Sgn == 1) { | |
56cbeefb | 469 | |
fabf4d8e | 470 | Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]}; |
16856d6e | 471 | bkg = FitFunction4Bkg(x,parbkg); |
fabf4d8e | 472 | } |
473 | ||
474 | sgn = FitFunction4Sgn(x,&par[3]); | |
475 | } | |
476 | ||
477 | if (ftypeOfFit4Bkg==3) { | |
478 | ||
479 | if(ftypeOfFit4Sgn == 0) { | |
480 | bkg=FitFunction4Bkg(x,par); | |
481 | sgn=FitFunction4Sgn(x,&par[1]); | |
482 | } | |
483 | if(ftypeOfFit4Sgn == 1) { | |
484 | Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]}; | |
485 | bkg=FitFunction4Bkg(x,parbkg); | |
486 | sgn=FitFunction4Sgn(x,&par[1]); | |
487 | } | |
488 | } | |
489 | ||
490 | total = bkg + sgn; | |
491 | ||
492 | return total; | |
493 | } | |
494 | ||
495 | //_________________________________________________________________________ | |
496 | Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){ | |
497 | // Fit function for the signal | |
498 | ||
499 | //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2) | |
500 | //Par: | |
501 | // * [0] = integralSgn | |
502 | // * [1] = mean | |
503 | // * [2] = sigma | |
504 | //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])] | |
505 | ||
506 | 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]); | |
507 | ||
508 | } | |
509 | ||
510 | //__________________________________________________________________________ | |
511 | ||
512 | Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){ | |
513 | // Fit function for the background | |
514 | ||
515 | Double_t maxDeltaM = 4.*fSigmaSgn; | |
516 | if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) { | |
517 | TF1::RejectPoint(); | |
518 | return 0; | |
519 | } | |
520 | Int_t firstPar=0; | |
521 | Double_t gaus2=0,total=-1; | |
522 | if(ftypeOfFit4Sgn == 1){ | |
523 | firstPar=3; | |
524 | //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2) | |
525 | //Par: | |
526 | // * [0] = integralSgn | |
527 | // * [1] = mean | |
528 | // * [2] = sigma | |
529 | //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])] | |
530 | gaus2 = FitFunction4Sgn(x,par); | |
531 | } | |
532 | ||
533 | switch (ftypeOfFit4Bkg){ | |
534 | case 0: | |
535 | //exponential | |
536 | //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max) | |
537 | //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written | |
538 | //as integralTot- integralGaus (=par [2]) | |
539 | //Par: | |
540 | // * [0] = integralBkg; | |
541 | // * [1] = B; | |
542 | //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x) | |
543 | 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]); | |
544 | break; | |
545 | case 1: | |
546 | //linear | |
547 | //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) | |
548 | // * [0] = integralBkg; | |
549 | // * [1] = b; | |
550 | total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass)); | |
551 | break; | |
552 | case 2: | |
553 | //polynomial | |
554 | //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) + | |
555 | //+ 1/3*c*(max^3-min^3) -> | |
556 | //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min) | |
557 | // * [0] = integralBkg; | |
558 | // * [1] = b; | |
559 | // * [2] = c; | |
560 | 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)); | |
561 | break; | |
562 | case 3: | |
563 | total=par[0+firstPar]; | |
564 | break; | |
565 | // default: | |
566 | // Types of Fit Functions for Background: | |
567 | // * 0 = exponential; | |
568 | // * 1 = linear; | |
569 | // * 2 = polynomial 2nd order | |
570 | // * 3 = no background"<<endl; | |
571 | ||
572 | } | |
573 | return total+gaus2; | |
574 | } | |
575 | ||
b7d4bc49 | 576 | //__________________________________________________________________________ |
577 | Bool_t AliHFMassFitter::SideBandsBounds(){ | |
578 | ||
579 | Double_t width=fhistoInvMass->GetBinWidth(8); | |
580 | if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX(); | |
581 | Double_t minHisto=fhistoInvMass->GetBinLowEdge(1); | |
582 | Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width; | |
583 | ||
584 | if(fMass-fminMass < 0 || fmaxMass-fMass <0) { | |
585 | cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl; | |
586 | return kFALSE; | |
587 | } | |
588 | ||
589 | if((fminMass!=minHisto || fmaxMass!= maxHisto) && (fMass-4.*fSigmaSgn-fminMass) <= 0){ | |
590 | Double_t coeff = (fMass-fminMass)/fSigmaSgn; | |
591 | ||
592 | fSideBandl=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width); | |
593 | fSideBandr=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width); | |
594 | cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl; | |
595 | if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl; | |
596 | ||
597 | if (coeff<2) { | |
598 | cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl; | |
599 | ftypeOfFit4Bkg=3; | |
600 | //set binleft and right without considering SetRangeFit- anyway no bkg! | |
601 | fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width); | |
602 | fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width); | |
603 | } | |
604 | } | |
605 | else { | |
606 | fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width); | |
607 | fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width); | |
a60f573d | 608 | // cout<<"\tfMass = "<<fMass<<"\tfSigmaSgn = "<<fSigmaSgn<<"\tminHisto = "<<minHisto<<endl; |
609 | // cout<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl; | |
b7d4bc49 | 610 | } |
611 | ||
612 | if (fSideBandl==0) { | |
613 | cout<<"Error! Range too little"; | |
614 | return kFALSE; | |
615 | } | |
b7d4bc49 | 616 | return kTRUE; |
617 | } | |
618 | ||
16856d6e | 619 | //__________________________________________________________________________ |
620 | ||
b7d4bc49 | 621 | void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{ |
622 | if (fSideBandl==0 && fSideBandr==0){ | |
623 | cout<<"Use MassFitter method first"<<endl; | |
624 | return; | |
625 | } | |
626 | left=fSideBandl; | |
627 | right=fSideBandr; | |
628 | } | |
629 | ||
fabf4d8e | 630 | //__________________________________________________________________________ |
631 | ||
34c79b83 | 632 | Bool_t AliHFMassFitter::MassFitter(Bool_t draw){ |
fabf4d8e | 633 | // Main method of the class: performs the fit of the histogram |
ad73ca8f | 634 | |
635 | //Set default fitter Minuit in order to use gMinuit in the contour plots | |
636 | TVirtualFitter::SetDefaultFitter("Minuit"); | |
637 | ||
fabf4d8e | 638 | Int_t nFitPars=0; //total function's number of parameters |
639 | switch (ftypeOfFit4Bkg){ | |
640 | case 0: | |
641 | nFitPars=5; //3+2 | |
642 | break; | |
643 | case 1: | |
644 | nFitPars=5; //3+2 | |
645 | break; | |
646 | case 2: | |
647 | nFitPars=6; //3+3 | |
648 | break; | |
649 | case 3: | |
650 | nFitPars=4; //3+1 | |
651 | break; | |
652 | } | |
653 | ||
654 | Int_t bkgPar = nFitPars-3; //background function's number of parameters | |
655 | ||
656 | cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl; | |
657 | ||
ad73ca8f | 658 | |
659 | TString listname="contourplot"; | |
660 | listname+=fcounter; | |
661 | if(!fContourGraph){ | |
662 | fContourGraph=new TList(); | |
663 | fContourGraph->SetOwner(); | |
664 | } | |
665 | ||
666 | fContourGraph->SetName(listname); | |
667 | ||
668 | ||
56cbeefb | 669 | //function names |
670 | TString bkgname="funcbkg"; | |
671 | TString bkg1name="funcbkg1"; | |
672 | TString massname="funcmass"; | |
673 | ||
fabf4d8e | 674 | //Total integral |
675 | Double_t totInt = fhistoInvMass->Integral("width"); | |
ad73ca8f | 676 | cout<<"Integral"<<endl; |
fabf4d8e | 677 | |
678 | fSideBands = kTRUE; | |
34c79b83 | 679 | Double_t width=fhistoInvMass->GetBinWidth(8); |
ad73ca8f | 680 | cout<<"fNbin"<<fNbin<<endl; |
34c79b83 | 681 | if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX(); |
682 | Double_t minHisto=fhistoInvMass->GetBinLowEdge(1); | |
683 | Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width; | |
b7d4bc49 | 684 | Bool_t ok=SideBandsBounds(); |
685 | if(!ok) return kFALSE; | |
34c79b83 | 686 | |
fabf4d8e | 687 | //sidebands integral - first approx (from histo) |
b7d4bc49 | 688 | Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width"); |
689 | cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl; | |
fabf4d8e | 690 | cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl; |
b7d4bc49 | 691 | if (sideBandsInt<=0) { |
692 | cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl; | |
693 | return kFALSE; | |
694 | } | |
695 | ||
fabf4d8e | 696 | /*Fit Bkg*/ |
697 | ||
16856d6e | 698 | |
34c79b83 | 699 | TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg"); |
56cbeefb | 700 | cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl; |
fabf4d8e | 701 | |
702 | funcbkg->SetLineColor(2); //red | |
703 | ||
704 | //first fit for bkg: approx bkgint | |
705 | ||
706 | switch (ftypeOfFit4Bkg) { | |
707 | case 0: //gaus+expo | |
708 | funcbkg->SetParNames("BkgInt","Slope"); | |
709 | funcbkg->SetParameters(sideBandsInt,-2.); | |
710 | break; | |
711 | case 1: | |
712 | funcbkg->SetParNames("BkgInt","Slope"); | |
713 | funcbkg->SetParameters(sideBandsInt,-100.); | |
714 | break; | |
715 | case 2: | |
716 | funcbkg->SetParNames("BkgInt","Coef1","Coef2"); | |
717 | funcbkg->SetParameters(sideBandsInt,-10.,5); | |
718 | break; | |
719 | case 3: | |
720 | if(ftypeOfFit4Sgn==0){ | |
fabf4d8e | 721 | funcbkg->SetParNames("Const"); |
722 | funcbkg->SetParameter(0,0.); | |
723 | funcbkg->FixParameter(0,0.); | |
724 | } | |
725 | break; | |
726 | default: | |
727 | cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl; | |
34c79b83 | 728 | return kFALSE; |
fabf4d8e | 729 | break; |
730 | } | |
731 | cout<<"\nBACKGROUND FIT - only combinatorial"<<endl; | |
732 | Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn; | |
733 | ||
734 | Double_t intbkg1=0,slope1=0,conc1=0; | |
56cbeefb | 735 | //if only signal and reflection: skip |
fabf4d8e | 736 | if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) { |
737 | ftypeOfFit4Sgn=0; | |
56cbeefb | 738 | fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0"); |
fabf4d8e | 739 | |
740 | for(Int_t i=0;i<bkgPar;i++){ | |
741 | fFitPars[i]=funcbkg->GetParameter(i); | |
742 | //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t"; | |
743 | fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i); | |
744 | //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl; | |
745 | } | |
b7d4bc49 | 746 | fSideBands = kFALSE; |
34c79b83 | 747 | //intbkg1 = funcbkg->GetParameter(0); |
748 | funcbkg->SetRange(fminMass,fmaxMass); | |
749 | intbkg1 = funcbkg->Integral(fminMass,fmaxMass); | |
fabf4d8e | 750 | if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1); |
751 | if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2); | |
6321ee46 | 752 | cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl; |
56cbeefb | 753 | } |
754 | else cout<<"\t\t//"<<endl; | |
fabf4d8e | 755 | |
756 | ftypeOfFit4Sgn=ftypeOfFit4SgnBkp; | |
757 | TF1 *funcbkg1=0; | |
758 | if (ftypeOfFit4Sgn == 1) { | |
759 | cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl; | |
760 | bkgPar+=3; | |
761 | ||
56cbeefb | 762 | cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl; |
763 | ||
764 | funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg"); | |
765 | cout<<"Function name = "<<funcbkg1->GetName()<<endl; | |
766 | ||
fabf4d8e | 767 | funcbkg1->SetLineColor(2); //red |
768 | ||
769 | if(ftypeOfFit4Bkg==2){ | |
770 | cout<<"*** Polynomial Fit ***"<<endl; | |
771 | funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2"); | |
772 | funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1); | |
6321ee46 | 773 | |
774 | //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl; | |
775 | //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl; | |
fabf4d8e | 776 | } else{ |
777 | if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened | |
778 | { | |
779 | cout<<"*** No background Fit ***"<<endl; | |
780 | funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const"); | |
781 | funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.); | |
782 | funcbkg1->FixParameter(3,0.); | |
783 | } else{ //expo or linear | |
784 | if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl; | |
785 | if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl; | |
786 | funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope"); | |
787 | funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1); | |
788 | } | |
789 | } | |
6321ee46 | 790 | Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0"); |
791 | if (status != 0){ | |
792 | cout<<"Minuit returned "<<status<<endl; | |
793 | return kFALSE; | |
794 | } | |
795 | ||
fabf4d8e | 796 | for(Int_t i=0;i<bkgPar;i++){ |
797 | fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i); | |
798 | //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i); | |
799 | fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i); | |
56cbeefb | 800 | //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl; |
fabf4d8e | 801 | } |
802 | ||
803 | intbkg1=funcbkg1->GetParameter(3); | |
804 | if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4); | |
805 | if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5); | |
806 | ||
807 | } else { | |
808 | bkgPar+=3; | |
809 | ||
810 | for(Int_t i=0;i<3;i++){ | |
811 | fFitPars[bkgPar-3+i]=0.; | |
56cbeefb | 812 | cout<<bkgPar-3+i<<"\t"<<0.<<"\t"; |
fabf4d8e | 813 | fFitPars[nFitPars+3*bkgPar-6+i]= 0.; |
56cbeefb | 814 | cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl; |
fabf4d8e | 815 | } |
816 | ||
817 | for(Int_t i=0;i<bkgPar-3;i++){ | |
818 | fFitPars[bkgPar+i]=funcbkg->GetParameter(i); | |
56cbeefb | 819 | cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t"; |
fabf4d8e | 820 | fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i); |
56cbeefb | 821 | cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl; |
fabf4d8e | 822 | } |
823 | ||
824 | ||
825 | } | |
826 | ||
827 | //sidebands integral - second approx (from fit) | |
828 | fSideBands = kFALSE; | |
829 | Double_t bkgInt; | |
b7d4bc49 | 830 | cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = "; |
fabf4d8e | 831 | if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass); |
832 | else bkgInt=funcbkg->Integral(fminMass,fmaxMass); | |
b7d4bc49 | 833 | cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl; |
fabf4d8e | 834 | |
835 | //Signal integral - first approx | |
836 | Double_t sgnInt; | |
837 | sgnInt = totInt-bkgInt; | |
838 | cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl; | |
6321ee46 | 839 | if (sgnInt <= 0){ |
840 | cout<<"Setting sgnInt = - sgnInt"<<endl; | |
841 | sgnInt=- sgnInt; | |
842 | } | |
fabf4d8e | 843 | /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */ |
56cbeefb | 844 | TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr"); |
845 | cout<<"Function name = "<<funcmass->GetName()<<endl<<endl; | |
fabf4d8e | 846 | funcmass->SetLineColor(4); //blue |
847 | ||
848 | //Set parameters | |
849 | cout<<"\nTOTAL FIT"<<endl; | |
850 | ||
851 | if(nFitPars==5){ | |
852 | funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma"); | |
853 | funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn); | |
56cbeefb | 854 | //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl; |
6321ee46 | 855 | //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl; |
fabf4d8e | 856 | funcmass->FixParameter(0,totInt); |
857 | } | |
858 | if (nFitPars==6){ | |
859 | funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma"); | |
860 | funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn); | |
6321ee46 | 861 | //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl; |
862 | //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl; | |
fabf4d8e | 863 | funcmass->FixParameter(0,totInt); |
864 | } | |
865 | if(nFitPars==4){ | |
866 | funcmass->SetParNames("Const","SgnInt","Mean","Sigma"); | |
56cbeefb | 867 | if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn); |
868 | else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn); | |
fabf4d8e | 869 | funcmass->FixParameter(0,0.); |
870 | //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl; | |
871 | //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl; | |
872 | ||
873 | } | |
6321ee46 | 874 | |
875 | Int_t status; | |
876 | ||
877 | status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0"); | |
878 | if (status != 0){ | |
879 | cout<<"Minuit returned "<<status<<endl; | |
880 | return kFALSE; | |
881 | } | |
882 | ||
fabf4d8e | 883 | cout<<"fit done"<<endl; |
0566386c | 884 | //reset value of fMass and fSigmaSgn to those found from fit |
885 | fMass=funcmass->GetParameter(nFitPars-2); | |
886 | fSigmaSgn=funcmass->GetParameter(nFitPars-1); | |
fabf4d8e | 887 | |
888 | for(Int_t i=0;i<nFitPars;i++){ | |
889 | fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i); | |
890 | fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i); | |
56cbeefb | 891 | //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl; |
fabf4d8e | 892 | } |
893 | /* | |
894 | //check: cout parameters | |
895 | for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){ | |
896 | cout<<i<<"\t"<<fFitPars[i]<<endl; | |
897 | } | |
898 | */ | |
74179930 | 899 | |
a60f573d | 900 | // if(draw){ |
901 | // TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName()); | |
902 | // TH1F *fhistocopy=new TH1F(*fhistoInvMass); | |
903 | // canvas->cd(); | |
904 | // fhistocopy->DrawClone(); | |
905 | // if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames"); | |
906 | // else funcbkg->DrawClone("sames"); | |
907 | // funcmass->DrawClone("sames"); | |
908 | // cout<<"Drawn"<<endl; | |
909 | // } | |
910 | ||
911 | if(funcmass->GetParameter(nFitPars-1) <0 || funcmass->GetParameter(nFitPars-2) <0 || funcmass->GetParameter(nFitPars-3) <0 ) { | |
912 | cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl; | |
b7d4bc49 | 913 | return kFALSE; |
914 | } | |
0566386c | 915 | |
56cbeefb | 916 | //increase counter of number of fits done |
917 | fcounter++; | |
918 | ||
ad73ca8f | 919 | //contour plots |
920 | if(draw){ | |
921 | ||
922 | for (Int_t kpar=1; kpar<nFitPars;kpar++){ | |
923 | ||
924 | for(Int_t jpar=kpar+1;jpar<nFitPars;jpar++){ | |
925 | cout<<"Par "<<kpar<<" and "<<jpar<<endl; | |
926 | ||
927 | // produce 2 contours per couple of parameters | |
928 | TGraph* cont[2] = {0x0, 0x0}; | |
929 | const Double_t errDef[2] = {1., 4.}; | |
930 | for (Int_t i=0; i<2; i++) { | |
931 | gMinuit->SetErrorDef(errDef[i]); | |
932 | cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar); | |
6321ee46 | 933 | cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl; |
ad73ca8f | 934 | } |
935 | ||
936 | if(!cont[0] || !cont[1]){ | |
937 | cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl; | |
938 | continue; | |
939 | } | |
940 | ||
941 | // set graph titles and add them to the list | |
942 | TString title = "Contour plot"; | |
943 | TString titleX = funcmass->GetParName(kpar); | |
944 | TString titleY = funcmass->GetParName(jpar); | |
945 | for (Int_t i=0; i<2; i++) { | |
946 | cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) ); | |
947 | cont[i]->SetTitle(title); | |
948 | cont[i]->GetXaxis()->SetTitle(titleX); | |
949 | cont[i]->GetYaxis()->SetTitle(titleY); | |
950 | cont[i]->GetYaxis()->SetLabelSize(0.033); | |
951 | cont[i]->GetYaxis()->SetTitleSize(0.033); | |
952 | cont[i]->GetYaxis()->SetTitleOffset(1.67); | |
953 | ||
954 | fContourGraph->Add(cont[i]); | |
955 | } | |
956 | ||
957 | // plot them | |
958 | TString cvname = Form("c%d%d", kpar, jpar); | |
959 | TCanvas *c4=new TCanvas(cvname,cvname,600,600); | |
960 | c4->cd(); | |
961 | cont[1]->SetFillColor(38); | |
962 | cont[1]->Draw("alf"); | |
963 | cont[0]->SetFillColor(9); | |
964 | cont[0]->Draw("lf"); | |
965 | ||
966 | } | |
967 | ||
968 | } | |
969 | ||
970 | } | |
971 | ||
a60f573d | 972 | if (ftypeOfFit4Sgn == 1 && funcbkg1) { |
973 | delete funcbkg1; | |
974 | funcbkg1=NULL; | |
975 | } | |
976 | if (funcbkg) { | |
977 | delete funcbkg; | |
978 | funcbkg=NULL; | |
979 | } | |
980 | if (funcmass) { | |
981 | delete funcmass; | |
982 | funcmass=NULL; | |
983 | } | |
16856d6e | 984 | |
985 | AddFunctionsToHisto(); | |
a60f573d | 986 | if (draw) DrawFit(); |
16856d6e | 987 | |
988 | ||
34c79b83 | 989 | return kTRUE; |
990 | } | |
991 | //_________________________________________________________________________ | |
992 | Double_t AliHFMassFitter::GetChiSquare() const{ | |
993 | TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass"); | |
994 | return funcmass->GetChisquare(); | |
fabf4d8e | 995 | } |
56cbeefb | 996 | |
b7d4bc49 | 997 | //_________________________________________________________________________ |
998 | Double_t AliHFMassFitter::GetReducedChiSquare() const{ | |
999 | TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass"); | |
1000 | return funcmass->GetChisquare()/funcmass->GetNDF(); | |
1001 | } | |
1002 | ||
56cbeefb | 1003 | //*********output |
1004 | ||
fabf4d8e | 1005 | //_________________________________________________________________________ |
2f328d65 | 1006 | void AliHFMassFitter::GetFitPars(Float_t *vector) const { |
fabf4d8e | 1007 | // Return fit parameters |
1008 | ||
2f328d65 | 1009 | for(Int_t i=0;i<fParsSize;i++){ |
1010 | vector[i]=fFitPars[i]; | |
fabf4d8e | 1011 | } |
1012 | } | |
fabf4d8e | 1013 | |
fabf4d8e | 1014 | |
6321ee46 | 1015 | //_________________________________________________________________________ |
1016 | void AliHFMassFitter::IntS(Float_t *valuewitherror){ | |
1017 | Int_t index=fParsSize/2 - 3; | |
1018 | valuewitherror[0]=fFitPars[index]; | |
1019 | index=fParsSize - 3; | |
1020 | valuewitherror[1]=fFitPars[index]; | |
1021 | } | |
1022 | ||
1023 | ||
16856d6e | 1024 | //_________________________________________________________________________ |
1025 | void AliHFMassFitter::AddFunctionsToHisto(){ | |
fabf4d8e | 1026 | |
16856d6e | 1027 | cout<<"AddFunctionsToHisto called"<<endl; |
56cbeefb | 1028 | TString bkgname = "funcbkg"; |
ab0fa301 | 1029 | Int_t np=-99; |
56cbeefb | 1030 | switch (ftypeOfFit4Bkg){ |
16856d6e | 1031 | case 0: //expo |
56cbeefb | 1032 | np=2; |
1033 | break; | |
16856d6e | 1034 | case 1: //linear |
56cbeefb | 1035 | np=2; |
1036 | break; | |
16856d6e | 1037 | case 2: //pol2 |
56cbeefb | 1038 | np=3; |
1039 | break; | |
16856d6e | 1040 | case 3: //no bkg |
56cbeefb | 1041 | np=1; |
1042 | break; | |
1043 | } | |
1044 | if (ftypeOfFit4Sgn == 1) { | |
1045 | bkgname += 1; | |
1046 | np+=3; | |
1047 | } | |
16856d6e | 1048 | |
1049 | TString bkgnamesave=bkgname; | |
1050 | TString testname=bkgname; | |
1051 | testname += "FullRange"; | |
1052 | TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data()); | |
1053 | if(testfunc){ | |
1054 | cout<<"AddFunctionsToHisto already used: exiting...."<<endl; | |
1055 | return; | |
1056 | } | |
1057 | ||
1058 | TList *hlist=fhistoInvMass->GetListOfFunctions(); | |
34c79b83 | 1059 | hlist->ls(); |
16856d6e | 1060 | |
34c79b83 | 1061 | TF1 *b=(TF1*)hlist->FindObject(bkgname.Data()); |
a60f573d | 1062 | if(!b){ |
1063 | cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl; | |
1064 | return; | |
1065 | } | |
1066 | ||
56cbeefb | 1067 | bkgname += "FullRange"; |
16856d6e | 1068 | TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg"); |
1069 | //cout<<bfullrange->GetName()<<endl; | |
56cbeefb | 1070 | for(Int_t i=0;i<np;i++){ |
16856d6e | 1071 | //cout<<i<<" di "<<np<<endl; |
1072 | bfullrange->SetParName(i,b->GetParName(i)); | |
1073 | bfullrange->SetParameter(i,b->GetParameter(i)); | |
1074 | bfullrange->SetParError(i,b->GetParError(i)); | |
1075 | } | |
a60f573d | 1076 | bfullrange->SetLineStyle(4); |
1077 | bfullrange->SetLineColor(14); | |
16856d6e | 1078 | |
1079 | bkgnamesave += "Recalc"; | |
1080 | ||
1081 | TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg"); | |
1082 | ||
1083 | TF1 *mass=fhistoInvMass->GetFunction("funcmass"); | |
1084 | ||
1085 | if (!mass){ | |
1086 | cout<<"funcmass doesn't exist "<<endl; | |
1087 | return; | |
1088 | } | |
1089 | ||
1090 | blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(np)); | |
1091 | blastpar->SetParError(0,mass->GetParError(np)); | |
1092 | if (np>=2) { | |
1093 | blastpar->SetParameter(1,mass->GetParameter(1)); | |
1094 | blastpar->SetParError(1,mass->GetParError(1)); | |
56cbeefb | 1095 | } |
16856d6e | 1096 | if (np==3) { |
1097 | blastpar->SetParameter(2,mass->GetParameter(2)); | |
1098 | blastpar->SetParError(2,mass->GetParError(2)); | |
1099 | } | |
1100 | ||
a60f573d | 1101 | blastpar->SetLineStyle(1); |
16856d6e | 1102 | blastpar->SetLineColor(2); |
1103 | ||
1104 | hlist->Add(bfullrange->Clone()); | |
1105 | hlist->Add(blastpar->Clone()); | |
1106 | hlist->ls(); | |
34c79b83 | 1107 | |
a60f573d | 1108 | if(bfullrange) { |
1109 | delete bfullrange; | |
1110 | bfullrange=NULL; | |
1111 | } | |
1112 | if(blastpar) { | |
1113 | delete blastpar; | |
1114 | blastpar=NULL; | |
1115 | } | |
16856d6e | 1116 | } |
1117 | ||
1118 | //_________________________________________________________________________ | |
1119 | ||
1120 | TH1F* AliHFMassFitter::GetHistoClone() const{ | |
1121 | ||
1122 | TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName()); | |
1123 | return hout; | |
1124 | } | |
1125 | //_________________________________________________________________________ | |
1126 | ||
1127 | void AliHFMassFitter::WriteHisto(TString path) { | |
1128 | if (fcounter == 0) { | |
1129 | cout<<"Use MassFitter method before WriteHisto"<<endl; | |
1130 | return; | |
1131 | } | |
1132 | TH1F* hget=(TH1F*)fhistoInvMass->Clone(); | |
1133 | ||
56cbeefb | 1134 | path += "HFMassFitterOutput.root"; |
1135 | TFile *output; | |
34c79b83 | 1136 | |
56cbeefb | 1137 | if (fcounter == 1) output = new TFile(path.Data(),"recreate"); |
1138 | else output = new TFile(path.Data(),"update"); | |
1139 | output->cd(); | |
1140 | hget->Write(); | |
ad73ca8f | 1141 | fContourGraph->Write(); |
56cbeefb | 1142 | output->Close(); |
fabf4d8e | 1143 | |
56cbeefb | 1144 | cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl; |
fabf4d8e | 1145 | |
a60f573d | 1146 | if(output) { |
1147 | delete output; | |
1148 | output=NULL; | |
1149 | } | |
fabf4d8e | 1150 | |
fabf4d8e | 1151 | } |
1152 | ||
1153 | //_________________________________________________________________________ | |
1154 | ||
56cbeefb | 1155 | void AliHFMassFitter::WriteNtuple(TString path) const{ |
a60f573d | 1156 | //TNtuple* nget=(TNtuple*)fntuParam->Clone(); |
56cbeefb | 1157 | path += "HFMassFitterOutput.root"; |
1158 | TFile *output = new TFile(path.Data(),"update"); | |
1159 | output->cd(); | |
a60f573d | 1160 | fntuParam->Write(); |
1161 | //nget->Write(); | |
56cbeefb | 1162 | output->Close(); |
a60f573d | 1163 | //cout<<nget->GetName()<<" written in "<<path<<endl; |
1164 | cout<<fntuParam->GetName()<<" written in "<<path<<endl; | |
1165 | /* | |
1166 | if(nget) { | |
1167 | //delete nget; | |
1168 | nget=NULL; | |
1169 | } | |
1170 | */ | |
1171 | if(output) { | |
1172 | delete output; | |
1173 | output=NULL; | |
1174 | } | |
1175 | } | |
1176 | ||
1177 | //_________________________________________________________________________ | |
1178 | ||
1179 | void AliHFMassFitter::DrawFit() const{ | |
1180 | ||
1181 | TString cvtitle="fit of "; | |
1182 | cvtitle+=fhistoInvMass->GetName(); | |
1183 | TString cvname="c"; | |
1184 | cvname+=fcounter; | |
1185 | ||
1186 | TCanvas *c=new TCanvas(cvname,cvtitle); | |
1187 | c->cd(); | |
1188 | GetHistoClone()->Draw("hist"); | |
1189 | GetHistoClone()->GetFunction("funcbkgFullRange")->Draw("same"); | |
1190 | GetHistoClone()->GetFunction("funcbkgRecalc")->Draw("same"); | |
1191 | GetHistoClone()->GetFunction("funcmass")->Draw("same"); | |
1192 | ||
1193 | ||
fabf4d8e | 1194 | } |
fabf4d8e | 1195 | |
fabf4d8e | 1196 | |
ad73ca8f | 1197 | //_________________________________________________________________________ |
1198 | ||
1199 | void AliHFMassFitter::PrintParTitles() const{ | |
1200 | TF1 *f=fhistoInvMass->GetFunction("funcmass"); | |
1201 | if(!f) { | |
1202 | cout<<"Fit function not found"<<endl; | |
1203 | return; | |
1204 | } | |
1205 | ||
1206 | Int_t np=0; | |
1207 | switch (ftypeOfFit4Bkg){ | |
1208 | case 0: //expo | |
1209 | np=2; | |
1210 | break; | |
1211 | case 1: //linear | |
1212 | np=2; | |
1213 | break; | |
1214 | case 2: //pol2 | |
1215 | np=3; | |
1216 | break; | |
1217 | case 3: //no bkg | |
1218 | np=1; | |
1219 | break; | |
1220 | } | |
1221 | ||
1222 | np+=3; //3 parameter for signal | |
1223 | if (ftypeOfFit4Sgn == 1) np+=3; | |
1224 | ||
1225 | cout<<"Parameter Titles \n"; | |
1226 | for(Int_t i=0;i<np;i++){ | |
1227 | cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl; | |
1228 | } | |
1229 | cout<<endl; | |
1230 | ||
1231 | } | |
1232 | ||
56cbeefb | 1233 | //************ significance |
fabf4d8e | 1234 | |
1235 | //_________________________________________________________________________ | |
1236 | ||
1237 | void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const { | |
1238 | // Return signal integral in mean+- n sigma | |
1239 | ||
16856d6e | 1240 | if(fcounter==0) { |
1241 | cout<<"Use MassFitter method before Signal"<<endl; | |
1242 | return; | |
1243 | } | |
56cbeefb | 1244 | |
1245 | //functions names | |
16856d6e | 1246 | TString bkgname= "funcbkgRecalc"; |
1247 | TString bkg1name="funcbkg1Recalc"; | |
56cbeefb | 1248 | TString massname="funcmass"; |
1249 | ||
b7d4bc49 | 1250 | |
1251 | TF1 *funcbkg=0; | |
56cbeefb | 1252 | TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data()); |
16856d6e | 1253 | if(!funcmass){ |
1254 | cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl; | |
1255 | return; | |
1256 | } | |
1257 | ||
b7d4bc49 | 1258 | if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data()); |
1259 | else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data()); | |
b7d4bc49 | 1260 | |
16856d6e | 1261 | if(!funcbkg){ |
1262 | cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl; | |
1263 | return; | |
1264 | } | |
b7d4bc49 | 1265 | |
16856d6e | 1266 | Int_t np=-99; |
1267 | switch (ftypeOfFit4Bkg){ | |
1268 | case 0: //expo | |
1269 | np=2; | |
1270 | break; | |
1271 | case 1: //linear | |
1272 | np=2; | |
1273 | break; | |
1274 | case 2: //pol2 | |
1275 | np=3; | |
1276 | break; | |
1277 | case 3: //no bkg | |
1278 | np=1; | |
1279 | break; | |
b7d4bc49 | 1280 | } |
1281 | ||
fabf4d8e | 1282 | Double_t intS,intSerr; |
1283 | ||
b7d4bc49 | 1284 | //relative error evaluation |
16856d6e | 1285 | intS=funcmass->GetParameter(np); |
1286 | intSerr=funcmass->GetParError(np); | |
b7d4bc49 | 1287 | cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl; |
16856d6e | 1288 | Double_t background,errbackground; |
1289 | Background(nOfSigma,background,errbackground); | |
1290 | ||
b7d4bc49 | 1291 | //signal +/- error in nsigma |
1292 | Double_t min=fMass-nOfSigma*fSigmaSgn; | |
1293 | Double_t max=fMass+nOfSigma*fSigmaSgn; | |
16856d6e | 1294 | |
1295 | Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4); | |
1296 | signal=mass - background; | |
1297 | errsignal=TMath::Sqrt((intSerr/intS*mass)*(intSerr/intS*mass)/*assume relative error is the same as for total integral*/ + errbackground*errbackground); | |
fabf4d8e | 1298 | return; |
1299 | } | |
16856d6e | 1300 | |
fabf4d8e | 1301 | //_________________________________________________________________________ |
1302 | ||
16856d6e | 1303 | void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const { |
1304 | ||
1305 | // Return signal integral in a range | |
1306 | ||
1307 | if(fcounter==0) { | |
1308 | cout<<"Use MassFitter method before Signal"<<endl; | |
1309 | return; | |
1310 | } | |
1311 | ||
56cbeefb | 1312 | //functions names |
16856d6e | 1313 | TString bkgname="funcbkgRecalc"; |
1314 | TString bkg1name="funcbkg1Recalc"; | |
56cbeefb | 1315 | TString massname="funcmass"; |
fabf4d8e | 1316 | |
b7d4bc49 | 1317 | |
1318 | TF1 *funcbkg=0; | |
56cbeefb | 1319 | TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data()); |
16856d6e | 1320 | if(!funcmass){ |
1321 | cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl; | |
1322 | return; | |
1323 | } | |
1324 | ||
b7d4bc49 | 1325 | if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data()); |
1326 | else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data()); | |
1327 | ||
16856d6e | 1328 | if(!funcbkg){ |
1329 | cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl; | |
1330 | return; | |
1331 | } | |
1332 | ||
1333 | Int_t np=-99; | |
1334 | switch (ftypeOfFit4Bkg){ | |
1335 | case 0: //expo | |
1336 | np=2; | |
1337 | break; | |
1338 | case 1: //linear | |
1339 | np=2; | |
1340 | break; | |
1341 | case 2: //pol2 | |
1342 | np=3; | |
1343 | break; | |
1344 | case 3: //no bkg | |
1345 | np=1; | |
1346 | break; | |
1347 | } | |
b7d4bc49 | 1348 | |
16856d6e | 1349 | Double_t intS,intSerr; |
b7d4bc49 | 1350 | |
16856d6e | 1351 | //relative error evaluation |
1352 | intS=funcmass->GetParameter(np); | |
1353 | intSerr=funcmass->GetParError(np); | |
fabf4d8e | 1354 | |
16856d6e | 1355 | cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl; |
1356 | Double_t background,errbackground; | |
1357 | Background(min,max,background,errbackground); | |
b7d4bc49 | 1358 | |
16856d6e | 1359 | //signal +/- error in the range |
1360 | ||
1361 | Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4); | |
1362 | signal=mass - background; | |
a60f573d | 1363 | errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/ |
16856d6e | 1364 | |
1365 | } | |
1366 | ||
1367 | //_________________________________________________________________________ | |
1368 | ||
1369 | void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const { | |
1370 | // Return background integral in mean+- n sigma | |
1371 | ||
1372 | if(fcounter==0) { | |
1373 | cout<<"Use MassFitter method before Background"<<endl; | |
1374 | return; | |
b7d4bc49 | 1375 | } |
1376 | ||
16856d6e | 1377 | //functions names |
1378 | TString bkgname="funcbkgRecalc"; | |
1379 | TString bkg1name="funcbkg1Recalc"; | |
1380 | ||
1381 | TF1 *funcbkg=0; | |
1382 | if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data()); | |
1383 | else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data()); | |
1384 | if(!funcbkg){ | |
1385 | cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl; | |
1386 | return; | |
1387 | } | |
1388 | ||
1389 | Double_t intB,intBerr; | |
1390 | ||
1391 | //relative error evaluation: from final parameters of the fit | |
1392 | if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl; | |
1393 | else{ | |
1394 | intB=funcbkg->GetParameter(0); | |
1395 | intBerr=funcbkg->GetParError(0); | |
1396 | cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl; | |
fabf4d8e | 1397 | } |
1398 | ||
b7d4bc49 | 1399 | Double_t min=fMass-nOfSigma*fSigmaSgn; |
1400 | Double_t max=fMass+nOfSigma*fSigmaSgn; | |
1401 | ||
1402 | //relative error evaluation: from histo | |
1403 | ||
1404 | intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin); | |
1405 | Double_t sum2=0; | |
1406 | for(Int_t i=1;i<=fSideBandl;i++){ | |
1407 | sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i); | |
1408 | } | |
1409 | for(Int_t i=fSideBandr;i<=fNbin;i++){ | |
1410 | sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i); | |
1411 | } | |
56cbeefb | 1412 | |
b7d4bc49 | 1413 | intBerr=TMath::Sqrt(sum2); |
1414 | cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl; | |
1415 | ||
16856d6e | 1416 | cout<<"Last estimation of bkg error is used"<<endl; |
b7d4bc49 | 1417 | |
1418 | //backround +/- error in nsigma | |
34c79b83 | 1419 | if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) { |
1420 | background = 0; | |
1421 | errbackground = 0; | |
1422 | } | |
1423 | else{ | |
1424 | background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2); | |
b7d4bc49 | 1425 | errbackground=intBerr/intB*background; // assume relative error is the same as for total integral |
16856d6e | 1426 | //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl; |
1427 | } | |
1428 | return; | |
1429 | ||
1430 | } | |
1431 | //___________________________________________________________________________ | |
1432 | ||
1433 | void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const { | |
1434 | // Return background integral in a range | |
1435 | ||
1436 | if(fcounter==0) { | |
1437 | cout<<"Use MassFitter method before Background"<<endl; | |
1438 | return; | |
1439 | } | |
1440 | ||
1441 | //functions names | |
1442 | TString bkgname="funcbkgRecalc"; | |
1443 | TString bkg1name="funcbkg1Recalc"; | |
1444 | ||
1445 | TF1 *funcbkg=0; | |
1446 | if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data()); | |
1447 | else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data()); | |
1448 | if(!funcbkg){ | |
1449 | cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl; | |
1450 | return; | |
1451 | } | |
1452 | ||
1453 | ||
1454 | Double_t intB,intBerr; | |
1455 | ||
1456 | //relative error evaluation: from final parameters of the fit | |
1457 | if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl; | |
1458 | else{ | |
1459 | intB=funcbkg->GetParameter(0); | |
1460 | intBerr=funcbkg->GetParError(0); | |
1461 | cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl; | |
1462 | } | |
1463 | ||
1464 | //relative error evaluation: from histo | |
1465 | ||
1466 | intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin); | |
1467 | Double_t sum2=0; | |
1468 | for(Int_t i=1;i<=fSideBandl;i++){ | |
1469 | sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i); | |
1470 | } | |
1471 | for(Int_t i=fSideBandr;i<=fNbin;i++){ | |
1472 | sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i); | |
1473 | } | |
1474 | ||
1475 | intBerr=TMath::Sqrt(sum2); | |
1476 | cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl; | |
1477 | ||
1478 | cout<<"Last estimation of bkg error is used"<<endl; | |
1479 | ||
1480 | //backround +/- error in the range | |
1481 | if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) { | |
1482 | background = 0; | |
1483 | errbackground = 0; | |
1484 | } | |
1485 | else{ | |
1486 | background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2); | |
1487 | errbackground=intBerr/intB*background; // assume relative error is the same as for total integral | |
1488 | //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl; | |
34c79b83 | 1489 | } |
fabf4d8e | 1490 | return; |
1491 | ||
1492 | } | |
1493 | ||
16856d6e | 1494 | |
fabf4d8e | 1495 | //__________________________________________________________________________ |
1496 | ||
b7d4bc49 | 1497 | void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const { |
fabf4d8e | 1498 | // Return significance in mean+- n sigma |
1499 | ||
1500 | Double_t signal,errsignal,background,errbackground; | |
1501 | Signal(nOfSigma,signal,errsignal); | |
1502 | Background(nOfSigma,background,errbackground); | |
1503 | ||
1504 | significance = signal/TMath::Sqrt(signal+background); | |
b7d4bc49 | 1505 | |
fabf4d8e | 1506 | errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal); |
b7d4bc49 | 1507 | |
fabf4d8e | 1508 | return; |
1509 | } | |
1510 | ||
1511 | //__________________________________________________________________________ | |
74179930 | 1512 | |
16856d6e | 1513 | void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const { |
1514 | // Return significance integral in a range | |
1515 | ||
1516 | Double_t signal,errsignal,background,errbackground; | |
1517 | Signal(min, max,signal,errsignal); | |
1518 | Background(min, max,background,errbackground); | |
1519 | ||
1520 | significance = signal/TMath::Sqrt(signal+background); | |
1521 | ||
1522 | errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal); | |
1523 | ||
1524 | return; | |
1525 | } | |
ad73ca8f | 1526 | |
1527 |