]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/SPECTRA/AliITSsadEdxFitter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity fixes for BUFFER_SIZE
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliITSsadEdxFitter.cxx
CommitLineData
fdfb7ea9 1/**************************************************************************\r
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15\r
16/* $Id$ */\r
17\r
18///////////////////////////////////////////////////////////////////////\r
19// Class with the fits algorithms to be used in the identified //\r
20// spectra analysis using the ITS in stand-alone mode //\r
21// Author: E.Biolcati, biolcati@to.infn.it //\r
22// F.Prino, prino@to.infn.it //\r
23///////////////////////////////////////////////////////////////////////\r
24\r
25#include <Riostream.h>\r
26#include <TLatex.h>\r
27#include <TH1F.h>\r
28#include <TF1.h>\r
29#include <TH1D.h>\r
30#include <TLine.h>\r
31#include <TH2F.h>\r
32#include <TMath.h>\r
33#include <TGraph.h>\r
34#include <TGraphErrors.h>\r
35#include <TLegend.h>\r
36#include <TLegendEntry.h>\r
37#include <TStyle.h>\r
38#include <Rtypes.h>\r
39#include "AliITSsadEdxFitter.h"\r
40\r
41\r
42ClassImp(AliITSsadEdxFitter)\r
7d629884 43//______________________________________________________________________\r
44AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(){\r
45 // standard constructor\r
46 for(Int_t i=0; i<3; i++) fFitpar[i] = 0.;\r
47 for(Int_t i=0; i<3; i++) fFitparErr[i] = 0.;\r
48 SetRangeStep1();\r
49 SetRangeStep2();\r
50 SetRangeStep3();\r
51 SetRangeFinalStep();\r
52 SetLimitsOnSigmaPion();\r
53 SetLimitsOnSigmaKaon();\r
54 SetLimitsOnSigmaProton();\r
55 SetBinsUsedPion();\r
56 SetBinsUsedKaon();\r
57 SetBinsUsedProton();\r
58};\r
fdfb7ea9 59\r
60//________________________________________________________\r
61Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x,Bool_t mc){\r
7d629884 62 // calculate the sigma 12-ott-2010 \r
63 Double_t p[2]={0.};\r
64 Double_t xMinKaon=0.15; //minimum pt value to consider the kaon peak\r
65 Double_t xMinProton=0.3;//minimum pt value to consider the proton peak\r
66 if(mc){\r
67 if(code==211){\r
68 p[0] = -1.20337e-04;\r
69 p[1] = 1.13060e-01;\r
70 } \r
71 else if(code==321 && x>xMinKaon){\r
72 p[0] = -2.39631e-03;\r
73 p[1] = 1.15723e-01;\r
74 } \r
75 else if(code==2212 && x>xMinProton){\r
76 p[0] = -8.34576e-03;\r
77 p[1] = 1.34237e-01;\r
78 } \r
79 else return -1;\r
80 }\r
81 else{\r
82 if(code==211){\r
83 p[0] = -6.55200e-05;\r
84 p[1] = 1.26657e-01;\r
85 } \r
86 else if(code==321 && x>xMinKaon){\r
87 p[0] = -6.22639e-04;\r
88 p[1] = 1.43289e-01;\r
89 } \r
90 else if(code==2212 && x>xMinProton){\r
91 p[0] = -2.13243e-03;\r
92 p[1] = 1.68614e-01;\r
93 } \r
94 else return -1;\r
95 }\r
96 return p[0]/(x*x)*TMath::Log(x)+p[1];\r
fdfb7ea9 97}\r
98\r
99//_______________________________________________________\r
100Int_t AliITSsadEdxFitter::CalcMean(Int_t code, Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2){\r
7d629884 101 // calculate the mean 12-ott-2010 \r
102 Double_t p1[4]={0.};\r
103 Double_t p2[4]={0.};\r
104 if(code==211){\r
105 p1[0] = 1.77049;\r
106 p1[1] = -2.65469;\r
107 p2[0] = 0.890856;\r
108 p2[1] = -0.276719;\r
109 mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;\r
110 mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;\r
111 }\r
112 else if(code==321){\r
113 p2[0] = 1.57664;\r
114 p2[1] = -6.88635;\r
115 p2[2] = 18.702;\r
116 p2[3] = -16.3385;\r
117 mean1 = 0.;\r
118 mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;\r
119 }\r
120 else if(code==2212){\r
121 p1[0] = 4.24861; \r
122 p1[1] = -19.178;\r
123 p1[2] = 31.5947;\r
124 p1[3] = -18.178;\r
125 mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;\r
126 mean2 = 0.; \r
127 }\r
128 else return -1;\r
129 return 0;\r
fdfb7ea9 130}\r
131\r
132//________________________________________________________\r
133Bool_t AliITSsadEdxFitter::IsGoodBin(Int_t bin,Int_t code){\r
7d629884 134 //method to select the bins used for the analysis only\r
135 Bool_t retvalue=kTRUE;\r
136 TLine *l[2]; //for the cross \r
137 l[0]=new TLine(-2.1,0,2.,100.);\r
138 l[1]=new TLine(-1.9,120,2.,0.);\r
139 for(Int_t j=0;j<2;j++){\r
140 l[j]->SetLineColor(4);\r
141 l[j]->SetLineWidth(5);\r
142 }\r
143\r
144 if(code==211 && (bin<fBinsUsedPion[0] || bin>fBinsUsedPion[1])){\r
145 for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
146 retvalue=kFALSE;\r
147 }\r
148 if(code==321 && (bin<fBinsUsedKaon[0] || bin>fBinsUsedKaon[1])){\r
149 for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
150 retvalue=kFALSE;\r
151 }\r
152 if(code==2212 && (bin<fBinsUsedProton[0] || bin>fBinsUsedProton[1])){\r
153 for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
154 retvalue=kFALSE;\r
155 }\r
156 return retvalue;\r
fdfb7ea9 157}\r
158\r
159//________________________________________________________\r
160Double_t SingleGausTail(const Double_t *x, const Double_t *par){\r
7d629884 161 //single gaussian with exponential tail\r
162 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
163 Double_t xx = x[0];\r
164 Double_t mean1 = par[1];\r
165 Double_t sigma1 = par[2];\r
166 Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
167 Double_t sigmaSquare1 = sigma1 * sigma1;\r
168 Double_t xdiv1 = mean1 + par[3] * sigma1;\r
169 Double_t y1=0.0;\r
170 if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
171 else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
172 return y1;\r
fdfb7ea9 173}\r
174\r
175//________________________________________________________\r
176Double_t DoubleGausTail(const Double_t *x, const Double_t *par){\r
7d629884 177 //sum of two gaussians with exponential tail\r
178 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
179 Double_t xx = x[0];\r
180 Double_t mean1 = par[1];\r
181 Double_t sigma1 = par[2];\r
182 Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
183 Double_t sigmaSquare1 = sigma1 * sigma1;\r
184 Double_t xdiv1 = mean1 + par[3] * sigma1;\r
185 Double_t y1=0.0;\r
186 if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
187 else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
188\r
189 Double_t mean2 = par[6];\r
190 Double_t sigma2 = par[7];\r
191 Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);\r
192 Double_t sigmaSquare2 = sigma2 * sigma2;\r
193 Double_t xdiv2 = mean2 + par[8] * sigma2;\r
194 Double_t y2=0.0;\r
195 if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);\r
196 else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));\r
197 return y1+y2;\r
fdfb7ea9 198}\r
199\r
200//________________________________________________________\r
201Double_t FinalGausTail(const Double_t *x, const Double_t *par){\r
7d629884 202 //sum of three gaussians with exponential tail\r
203 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
204 Double_t xx = x[0];\r
205 Double_t mean1 = par[1];\r
206 Double_t sigma1 = par[2];\r
207 Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
208 Double_t sigmaSquare1 = sigma1 * sigma1;\r
209 Double_t xdiv1 = mean1 + par[3] * sigma1;\r
210 Double_t y1=0.0;\r
211 if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
212 else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
213\r
214 Double_t mean2 = par[6];\r
215 Double_t sigma2 = par[7];\r
216 Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);\r
217 Double_t sigmaSquare2 = sigma2 * sigma2;\r
218 Double_t xdiv2 = mean2 + par[8] * sigma2;\r
219 Double_t y2=0.0;\r
220 if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);\r
221 else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));\r
222\r
223 Double_t mean3 = par[11];\r
224 Double_t sigma3 = par[12];\r
225 Double_t xNormSquare3 = (xx - mean3) * (xx - mean3);\r
226 Double_t sigmaSquare3 = sigma3 * sigma3;\r
227 Double_t xdiv3 = mean3 + par[13] * sigma3;\r
228 Double_t y3=0.0;\r
229 if(xx < xdiv3) y3 = par[10]/(s2pi*par[12]) * TMath::Exp(-0.5 * xNormSquare3 / sigmaSquare3);\r
230 else y3 = TMath::Exp(-par[14]*(xx-xdiv3)) * par[10] / (s2pi*par[12]) * TMath::Exp(-0.5*(par[13]*par[13]));\r
231 return y1+y2+y3;\r
fdfb7ea9 232}\r
233\r
234//______________________________________________________________________\r
235void AliITSsadEdxFitter::CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const{\r
7d629884 236 //code to calculate the difference fit function and histogram point (residual)\r
237 //to use as goodness test for the fit\r
238 Int_t ipt=0;\r
239 Double_t x=0.,yhis=0.,yfun=0.;\r
240 for(Int_t i=0;i<h->GetNbinsX();i++){\r
241 x=(h->GetBinLowEdge(i+2)+h->GetBinLowEdge(i+1))/2;\r
242 yfun=fun->Eval(x);\r
243 yhis=h->GetBinContent(i+1);\r
244 if(yhis>0. && yfun>0.) {\r
245 gres->SetPoint(ipt,x,(yhis-yfun)/yhis);\r
246 ipt++;\r
247 }\r
248 }\r
249 return;\r
fdfb7ea9 250}\r
251\r
252\r
253//________________________________________________________\r
254Double_t SingleGausStep(const Double_t *x, const Double_t *par){\r
7d629884 255 //single normalized gaussian\r
256 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
257 Double_t xx = x[0];\r
258 Double_t mean1 = par[1];\r
259 Double_t sigma1 = par[2];\r
260 Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
261 Double_t sigma1Square = sigma1 * sigma1;\r
262 Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNorm1Square / sigma1Square);\r
263 return step1;\r
fdfb7ea9 264}\r
265\r
266//________________________________________________________\r
267Double_t DoubleGausStep(const Double_t *x, const Double_t *par){\r
7d629884 268 //sum of two normalized gaussians\r
269 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
270 Double_t xx = x[0];\r
271 Double_t mean1 = par[1];\r
272 Double_t sigma1 = par[2];\r
273 Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
274 Double_t sigma1Square = sigma1 * sigma1;\r
275 Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );\r
276 Double_t mean2 = par[4];\r
277 Double_t sigma2 = par[5];\r
278 Double_t xNorm2Square = (xx - mean2) * (xx - mean2);\r
279 Double_t sigma2Square = sigma2 * sigma2;\r
280 Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );\r
281 return step1+step2;\r
fdfb7ea9 282}\r
283\r
284//________________________________________________________\r
285Double_t FinalGausStep(const Double_t *x, const Double_t *par){\r
7d629884 286 //sum of three normalized gaussians\r
287 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
288 Double_t xx = x[0];\r
289 Double_t mean1 = par[1];\r
290 Double_t sigma1 = par[2];\r
291 Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
292 Double_t sigma1Square = sigma1 * sigma1;\r
293 Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );\r
294 Double_t mean2 = par[4];\r
295 Double_t sigma2 = par[5];\r
296 Double_t xNorm2Square = (xx - mean2) * (xx - mean2);\r
297 Double_t sigma2Square = sigma2 * sigma2;\r
298 Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );\r
299 Double_t mean3 = par[7];\r
300 Double_t sigma3 = par[8];\r
301 Double_t xNorm3Square = (xx - mean3) * (xx - mean3);\r
302 Double_t sigma3Square = sigma3 * sigma3;\r
303 Double_t step3 = par[6]/(s2pi*par[8]) * TMath::Exp( - 0.5 * xNorm3Square / sigma3Square);\r
304 return step1+step2+step3;\r
fdfb7ea9 305}\r
306\r
307//______________________________________________________________________\r
308Double_t AliITSsadEdxFitter::GausPlusTail(const Double_t x, const Double_t mean, Double_t rms, Double_t c, Double_t slope, Double_t cut ) const{\r
7d629884 309 //gaussian with an exponential tail on the right side\r
310 Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));\r
311 Double_t returnvalue=0.0;\r
312 Double_t n=0.5*(1.0+TMath::Erf(cut/TMath::Sqrt(2.0)))+TMath::Exp(-cut*cut*0.5)*factor/(TMath::Abs(rms)*slope);\r
313 if (x<mean+cut*rms) returnvalue=TMath::Exp(-1.0*(x-mean)*(x-mean)/(2.0*rms*rms))*factor/TMath::Abs(rms);\r
314 else returnvalue=TMath::Exp(slope*(mean+cut*rms-x))*TMath::Exp(-cut*cut*0.5)*factor/TMath::Abs(rms);\r
315 return c*returnvalue/n;\r
fdfb7ea9 316}\r
317\r
318\r
319//______________________________________________________________________\r
320Double_t AliITSsadEdxFitter::GausOnBackground(const Double_t* x, const Double_t *par) const {\r
7d629884 321 //gaussian with a background parametrisation \r
322 //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<x[0]<< endl;\r
323 Double_t returnvalue=0.0;\r
324 Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));\r
325 if(par[6]<0.0) returnvalue+=TMath::Exp(-1.0*(x[0]-par[0])*(x[0]-par[0])/(2.0*par[1]*par[1]))*par[2]* factor/TMath::Abs(par[1]);\r
326 else returnvalue+=GausPlusTail(x[0], par[0], par[1],par[2], par[6], 1.2);\r
327 returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);\r
328 return returnvalue;\r
fdfb7ea9 329}\r
330\r
331//______________________________________________________________________\r
332void AliITSsadEdxFitter::DrawFitFunction(TF1 *fun) const {\r
7d629884 333 //code to draw the final fit function and the single gaussians used\r
334 //to extract the yields for the 3 species\r
335 TF1 *fdraw1=new TF1("fdraw1",SingleGausStep,-3.5,3.5,3);\r
336 TF1 *fdraw2=new TF1("fdraw2",SingleGausStep,-3.5,3.5,3);\r
337 TF1 *fdraw3=new TF1("fdraw3",SingleGausStep,-3.5,3.5,3);\r
338 fdraw1->SetParameter(0,fun->GetParameter(0));\r
339 fdraw1->SetParameter(1,fun->GetParameter(1));\r
340 fdraw1->SetParameter(2,fun->GetParameter(2));\r
341 fdraw2->SetParameter(0,fun->GetParameter(3));\r
342 fdraw2->SetParameter(1,fun->GetParameter(4));\r
343 fdraw2->SetParameter(2,fun->GetParameter(5));\r
344 fdraw3->SetParameter(0,fun->GetParameter(6));\r
345 fdraw3->SetParameter(1,fun->GetParameter(7));\r
346 fdraw3->SetParameter(2,fun->GetParameter(8));\r
347\r
348 fdraw1->SetLineColor(6);//color code\r
349 fdraw2->SetLineColor(2);\r
350 fdraw3->SetLineColor(4); \r
351\r
352 fdraw1->SetLineStyle(2);//dot lines\r
353 fdraw2->SetLineStyle(2);\r
354 fdraw3->SetLineStyle(2);\r
355\r
356 fun->SetLineWidth(3);\r
357 fdraw1->DrawCopy("same");\r
358 fdraw2->DrawCopy("same");\r
359 fdraw3->DrawCopy("same");\r
360 fun->DrawCopy("same");\r
361\r
362 TLatex *ltx[3];\r
363 for(Int_t j=0;j<3;j++){\r
364 ltx[0]=new TLatex(0.13,0.25,"pions");\r
365 ltx[1]=new TLatex(0.13,0.20,"kaons");\r
366 ltx[2]=new TLatex(0.13,0.15,"protons");\r
367 ltx[0]->SetTextColor(6);\r
368 ltx[1]->SetTextColor(2);\r
369 ltx[2]->SetTextColor(4);\r
370 ltx[j]->SetNDC();\r
371 ltx[j]->Draw();\r
372 }\r
373 return;\r
fdfb7ea9 374}\r
375\r
376//______________________________________________________________________\r
377void AliITSsadEdxFitter::GetInitialParam(TH1F* h,Bool_t mc,Int_t code,Int_t bin, Float_t &pt, Float_t &ampl, Float_t &mean1, Float_t &mean2, Float_t &mean3, Float_t &sigma1, Float_t &sigma2, Float_t &sigma3){\r
7d629884 378 //code to get the expected values to use for fitting\r
379 Double_t xbins[23]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};\r
380 pt=(xbins[bin+1]+xbins[bin])/2;\r
381\r
382 //draw and label\r
383 h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],code));\r
384 h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");\r
385 h->GetYaxis()->SetTitle("counts");\r
386 h->Draw("e");\r
387 h->SetFillColor(11);\r
fdfb7ea9 388 \r
7d629884 389 //expected values\r
390 Int_t xmax=-1,ymax=-1,zmax=-1;\r
391 h->GetMaximumBin(xmax,ymax,zmax);\r
392 printf("\n---------------------------------- BIN %d - hypothesis %d ----------------------------------\n",bin,code);\r
393 Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
394 ampl = h->GetMaximum()/(h->GetRMS()*s2pi);\r
395 mean1 = h->GetBinLowEdge(xmax); //expected mean values\r
396 Int_t calcmean=CalcMean(code,pt,mean1,mean2,mean3);\r
397 if(calcmean<0) cout<<"Error during mean calculation"<<endl;\r
398 printf("mean values -> %f %f %f\n",mean1,mean2,mean3);\r
399 printf("integration ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n",fRangeStep1[0],fRangeStep1[1],fRangeStep2[0],fRangeStep2[1],fRangeStep3[0],fRangeStep3[1]);\r
400 sigma1 = CalcSigma(211,pt,mc); //expected sigma values\r
401 sigma2 = CalcSigma(321,pt,mc);\r
402 sigma3 = CalcSigma(2212,pt,mc);\r
1f75ac7d 403 printf("sigma values -> %f %f %f\n",sigma1,sigma2,sigma3);\r
7d629884 404 printf("sigma ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1],fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1],fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);\r
fdfb7ea9 405 return;\r
406}\r
407\r
408//________________________________________________________\r
409void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){\r
7d629884 410 // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
411 // pt bin from 0 to 20, code={211,321,2212} \r
412 // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed\r
413 // final step: refit all using the parameters and tollerance limits (+-20%)\r
414 TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
415 TString modfit = "M0R+";\r
416 Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
417 Int_t code=TMath::Abs(signedcode);\r
418 GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
419 if(!IsGoodBin(bin,code)) return;\r
420\r
421 printf("___________________________________________________________________ First Step: pions\n");\r
422 fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
423 fstep1->SetParameter(0,ampl); //initial ampl pion\r
424 fstep1->SetParameter(1,mean); //initial mean pion\r
425 fstep1->SetParameter(2,expPionSig); //initial sigma pion\r
426 fstep1->SetParLimits(0,0.,ampl*1.2); //limits ampl pion\r
427 fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]); //limits mean pion (dummy)\r
428 fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]); //limits sigma pion\r
429\r
430 if(expPionSig>0) h->Fit(fstep1,modfit.Data(),"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//first fit\r
431 else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
432\r
433 printf("___________________________________________________________________ Second Step: kaons\n");\r
434 fstep2 = new TF1("fstep2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);\r
435 fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
436 fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
437 fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
438 fstep2->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon\r
439 fstep2->SetParameter(4,expKaonMean); //initial mean kaon\r
440 fstep2->SetParameter(3,expKaonSig); //initial sigma kaon\r
441 fstep2->SetParLimits(3,0.,fstep1->GetParameter(0)); //limits ampl kaon\r
442 fstep2->SetParLimits(4,fstep1->GetParameter(1),fRangeStep4[1]); //limits mean kaon \r
443 fstep2->SetParLimits(5,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
444\r
445 if(expKaonSig>0) h->Fit(fstep2,modfit.Data(),"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);//second fit\r
446 else for(Int_t npar=3;npar<6;npar++) fstep2->FixParameter(npar,0.00001);\r
447\r
448 /*TLine *l[3];\r
449 l[0] = new TLine(expKaonMean,0,expKaonMean,10000);\r
450 l[1] = new TLine(expKaonMean+fRangeStep2[0],0,expKaonMean+fRangeStep2[0],10000);\r
451 l[2] = new TLine(expKaonMean+fRangeStep2[1],0,expKaonMean+fRangeStep2[1],10000);\r
452 for(Int_t dp=0;dp<3;dp++) {\r
453 l[dp]->Draw("same");\r
454 l[dp]->SetLineColor(4);\r
455 l[dp]->SetLineWidth(3);\r
456 }*/\r
457\r
458 printf("___________________________________________________________________ Third Step: protons\n");\r
459 fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
460 fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
461 fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
462 fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
463 fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl kaon\r
464 fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean kaon\r
465 fstep3->FixParameter(5,fstep2->GetParameter(5)); //fidex sigma kaon\r
466 fstep3->SetParameter(6,fstep2->GetParameter(0)/16.); //initial ampl proton\r
467 fstep3->SetParameter(7,expProtonMean); //initial mean proton\r
468 fstep3->SetParameter(8,expProtonSig); //initial sigma proton\r
469 fstep3->SetParLimits(6,0.,fstep2->GetParameter(0)); //limits ampl proton\r
470 fstep3->SetParLimits(7,fstep2->GetParameter(4),fRangeStep4[1]); //limits mean proton\r
471 fstep3->SetParLimits(8,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
472\r
473 if(expProtonSig>0) h->Fit(fstep3,modfit.Data(),"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//third fit\r
474 else for(Int_t npar=6;npar<9;npar++) fstep3->FixParameter(npar,0.00001);\r
475\r
476 printf("___________________________________________________________________ Final Step: refit all\n");\r
477 fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
478 fstepTot->SetLineColor(1);\r
479\r
480 Double_t initialParametersStepTot[9];\r
481 initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
482 initialParametersStepTot[1] = fstep1->GetParameter(1);\r
483 initialParametersStepTot[2] = fstep1->GetParameter(2);\r
484\r
485 initialParametersStepTot[3] = fstep2->GetParameter(3);//second gaussian\r
486 initialParametersStepTot[4] = fstep2->GetParameter(4);\r
487 initialParametersStepTot[5] = fstep2->GetParameter(5);\r
488\r
489 initialParametersStepTot[6] = fstep3->GetParameter(6);//third gaussian\r
490 initialParametersStepTot[7] = fstep3->GetParameter(7);\r
491 initialParametersStepTot[8] = fstep3->GetParameter(8); \r
492\r
493 fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
494 fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion\r
495 fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
496 fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
497 fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
498 fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
499 fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
500 fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton\r
501 fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
502 fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
503\r
504 h->Fit(fstepTot,modfit.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
505\r
506 //************************************* storing parameter to calculate the yields *******\r
507 Int_t chpa=0;\r
508 if(code==321) chpa=3;\r
509 if(code==2212) chpa=6;\r
510 for(Int_t j=0;j<3;j++) {\r
511 fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
512 fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
513 }\r
514\r
515 DrawFitFunction(fstepTot);\r
516 CalcResidual(h,fstepTot,gres);\r
517 return;\r
fdfb7ea9 518}\r
519\r
520//________________________________________________________\r
521void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){\r
7d629884 522 // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
523 // pt bin from 0 to 20, code={211,321,2212} \r
524 // first step: pion peak, second step: proton peak, third step: kaon peak\r
525 // final step: refit all using the parameters\r
526 TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
527 TString modfit = "M0R+";\r
528 Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
529 Int_t code=TMath::Abs(signedcode);\r
530 GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
531 if(!IsGoodBin(bin,code)) return;\r
532\r
533 printf("___________________________________________________________________ First Step: pion\n");\r
534 fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
535 fstep1->SetParameter(0,ampl); //initial ampl pion\r
536 fstep1->SetParameter(1,mean); //initial mean pion\r
537 fstep1->SetParameter(2,expPionSig); //initial sigma pion\r
538 fstep1->SetParLimits(0,0,ampl*1.2); //limits ampl pion\r
539 fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]); //limits mean pion (dummy)\r
540 fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]); //limits sigma pion\r
541\r
542 if(expPionSig>0) h->Fit(fstep1,modfit,"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//first fit\r
543 else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
544\r
545 printf("___________________________________________________________________ Second Step: proton\n");\r
546 fstep2 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
547 fstep2->SetParameter(0,fstep1->GetParameter(0)/16.);//initial ampl proton\r
548 fstep2->SetParameter(1,expProtonMean); //initial mean proton\r
549 fstep2->SetParameter(2,expProtonSig); //initial sigma proton\r
550 fstep2->SetParLimits(0,0.,fstep1->GetParameter(0)); //limits ampl proton\r
551 fstep2->SetParLimits(1,fstep1->GetParameter(1),fRangeStep4[1]); //limits mean proton\r
552 fstep2->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
553\r
554 if(expProtonSig>0) h->Fit(fstep2,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//second fit\r
555 else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);\r
556\r
557 printf("___________________________________________________________________ Third Step: kaon\n");\r
558 fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
559 fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
560 fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
561 fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
562 fstep3->FixParameter(6,fstep2->GetParameter(0)); //fixed ampl proton\r
563 fstep3->FixParameter(7,fstep2->GetParameter(1)); //fixed mean proton\r
564 fstep3->FixParameter(8,fstep2->GetParameter(2)); //fixed sigma proton\r
565 fstep3->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon\r
566 fstep3->SetParameter(4,expKaonMean); //initial mean kaon\r
567 fstep3->SetParameter(5,expKaonSig); //initial sigma kaon\r
568 fstep3->SetParLimits(3,fstep2->GetParameter(0),fstep1->GetParameter(0)); //limits ampl kaon\r
569 fstep3->SetParLimits(4,fstep1->GetParameter(1),fstep2->GetParameter(1)); //limits mean kaon\r
570 fstep3->SetParLimits(5,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
571 /*TLine *l[3];\r
572 l[0] = new TLine(expProtonMean,0,expProtonMean,10000);\r
573 l[1] = new TLine(expProtonMean+fRangeStep3[0],0,expProtonMean+fRangeStep3[0],10000);\r
574 l[2] = new TLine(expProtonMean+fRangeStep3[1],0,expProtonMean+fRangeStep3[1],10000);\r
575 for(Int_t dp=0;dp<3;dp++) {\r
576 l[dp]->Draw("same");\r
577 l[dp]->SetLineColor(2);\r
578 l[dp]->SetLineWidth(4);\r
579 }*/\r
580 if(expKaonSig>0) h->Fit(fstep3,modfit,"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);//third fit\r
581 else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
582\r
583 printf("___________________________________________________________________ Final Step: refit all\n");\r
584 fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
585 fstepTot->SetLineColor(1);\r
586\r
587 Double_t initialParametersStepTot[9];\r
588 initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
589 initialParametersStepTot[1] = fstep1->GetParameter(1);\r
590 initialParametersStepTot[2] = fstep1->GetParameter(2);\r
591\r
592 initialParametersStepTot[3] = fstep3->GetParameter(3);//second gaussian\r
593 initialParametersStepTot[4] = fstep3->GetParameter(4);\r
594 initialParametersStepTot[5] = fstep3->GetParameter(5);\r
595\r
596 initialParametersStepTot[6] = fstep2->GetParameter(0);//third gaussian\r
597 initialParametersStepTot[7] = fstep2->GetParameter(1);\r
598 initialParametersStepTot[8] = fstep2->GetParameter(2); \r
599\r
600 fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
601 fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion\r
602 fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
603 fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
604 fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
605 fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
606 fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
607 fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton\r
608 fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
609 fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
610\r
611 h->Fit(fstepTot,modfit,"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
612\r
613 //************************************* storing parameter to calculate the yields *******\r
614 Int_t chpa=0;\r
615 if(code==321) chpa=3;\r
616 if(code==2212) chpa=6;\r
617 for(Int_t j=0;j<3;j++) {\r
618 fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
619 fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
620 }\r
621\r
622 DrawFitFunction(fstepTot);\r
623 CalcResidual(h,fstepTot,gres);\r
624 return;\r
fdfb7ea9 625}\r
626\r
627//________________________________________________________\r
628void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){\r
7d629884 629 // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
630 // pt bin from 0 to 20, code={211,321,2212} \r
631 // first step: proton peak, second step: pion peak, third step: kaon peak\r
632 // final step: refit all using the parameters\r
633 TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
634 TString modfit = "M0R+";\r
635 Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
636 Int_t code=TMath::Abs(signedcode);\r
637 GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
638 if(!IsGoodBin(bin,code)) return;\r
639\r
640 printf("___________________________________________________________________ First Step: proton\n");\r
641 fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
642 fstep1->SetParameter(0,ampl/16.); //initial ampl proton`\r
643 fstep1->SetParameter(1,expProtonMean); //initial mean proton\r
644 fstep1->SetParameter(2,expProtonSig); //initial sigma proton\r
645 fstep1->SetParLimits(0,0,ampl); //limits ampl proton\r
646 fstep1->SetParLimits(1,mean,fRangeStep4[1]); //limits mean proton (dummy)\r
647 fstep1->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
648\r
649 if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//first fit\r
650 else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
651\r
652 printf("___________________________________________________________________ Second Step: pion\n");\r
653 fstep2 = new TF1("step2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);\r
654 fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton \r
655 fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton\r
656 fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton\r
657 fstep2->SetParameter(3,ampl); //initial ampl pion\r
658 fstep2->SetParameter(4,mean); //initial mean pion\r
659 fstep2->SetParameter(5,expPionSig); //initial sigma pion\r
660 fstep2->SetParLimits(3,0.,ampl); //limits ampl pion\r
661 fstep2->SetParLimits(4,fRangeStep4[0],fstep1->GetParameter(1)); //limits mean pion\r
662 fstep2->SetParLimits(5,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]); //limits sigma pion\r
663\r
664 if(expPionSig>0) h->Fit(fstep2,modfit,"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//second fit\r
665 else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);\r
666\r
667 printf("___________________________________________________________________ Third Step: kaon\n");\r
668 fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
669 fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton\r
670 fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton\r
671 fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton\r
672 fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl pion\r
673 fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean pion\r
674 fstep3->FixParameter(5,fstep2->GetParameter(5)); //fixed sigma pion\r
675 fstep3->SetParameter(6,fstep2->GetParameter(0)/8.); //initial ampl kaon\r
676 fstep3->SetParameter(7,expKaonMean); //initial mean kaon\r
677 fstep3->SetParameter(8,expKaonSig); //initial sigma kaon\r
678 fstep3->SetParLimits(6,fstep1->GetParameter(0),fstep2->GetParameter(3)); //limits ampl kaon\r
679 fstep3->SetParLimits(7,fstep2->GetParameter(4),fstep1->GetParameter(1)); //limits mean kaon\r
680 fstep3->SetParLimits(8,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
681 /*TLine *l[3];\r
682 l[0] = new TLine(expProtonMean,0,expProtonMean,10000);\r
683 l[1] = new TLine(expProtonMean+fRangeStep3[0],0,expProtonMean+fRangeStep3[0],10000);\r
684 l[2] = new TLine(expProtonMean+fRangeStep3[1],0,expProtonMean+fRangeStep3[1],10000);\r
685 for(Int_t dp=0;dp<3;dp++) {\r
686 l[dp]->Draw("same");\r
687 l[dp]->SetLineColor(2);\r
688 l[dp]->SetLineWidth(4);\r
689 }*/\r
690 if(expKaonSig>0) h->Fit(fstep3,modfit,"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);//third fit\r
691 else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
692\r
693 printf("___________________________________________________________________ Final Step: refit all\n");\r
694 fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
695 fstepTot->SetLineColor(1);\r
696\r
697 Double_t initialParametersStepTot[9];\r
698 initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
699 initialParametersStepTot[1] = fstep1->GetParameter(1);\r
700 initialParametersStepTot[2] = fstep1->GetParameter(2);\r
701\r
702 initialParametersStepTot[3] = fstep3->GetParameter(6);//second gaussian\r
703 initialParametersStepTot[4] = fstep3->GetParameter(7);\r
704 initialParametersStepTot[5] = fstep3->GetParameter(8);\r
705\r
706 initialParametersStepTot[6] = fstep2->GetParameter(3);//third gaussian\r
707 initialParametersStepTot[7] = fstep2->GetParameter(4);\r
708 initialParametersStepTot[8] = fstep2->GetParameter(5); \r
709\r
710 fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
711 fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl proton\r
712 fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
713 fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
714 fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
715 fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
716 fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
717 fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl pion\r
718 fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
719 fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
720\r
721 h->Fit(fstepTot,modfit,"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
722\r
723 //************************************* storing parameter to calculate the yields *******\r
724 Int_t chpa=0;\r
725 if(code==321) chpa=3;\r
726 if(code==211) chpa=6;\r
727 for(Int_t j=0;j<3;j++) {\r
728 fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
729 fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
730 }\r
731\r
732 DrawFitFunction(fstepTot);\r
733 CalcResidual(h,fstepTot,gres);\r
734 return;\r
fdfb7ea9 735}\r
736\r
737\r
738//________________________________________________________\r
739void AliITSsadEdxFitter::DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc){\r
7d629884 740 // single-gaussian fit to log(dedx)-log(dedxBB) histogram\r
741 TF1 *fstep1;\r
742 TString modfit = "M0R+";\r
743 Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
744 Int_t code=TMath::Abs(signedcode);\r
745 GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
746 if(!IsGoodBin(bin,code)) return;\r
747\r
748 printf("___________________________________________________________________ Single Step\n");\r
749 fstep1 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
750 fstep1->SetParameter(0,ampl/16.); //initial ampl \r
751 fstep1->SetParameter(1,expProtonMean); //initial mean \r
752 fstep1->SetParameter(2,expProtonSig); //initial sigma \r
753 fstep1->SetParLimits(0,0.,ampl); //limits ampl proton\r
754 fstep1->SetParLimits(1,mean,fRangeStep4[1]); //limits mean proton\r
755 //fstep1->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
756\r
757 if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//fit\r
758 else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
759\r
760 fstep1->SetLineColor(1);\r
761 fstep1->Draw("same");\r
762\r
763 fFitpar[0] = fstep1->GetParameter(0);\r
764 fFitparErr[0] = fstep1->GetParError(0);\r
765 return;\r
fdfb7ea9 766}\r
767\r
768//________________________________________________________\r
769void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t signedcode){\r
7d629884 770 // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
771 // pt bin from 0 to 20, code={211,321,2212} \r
772 // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed\r
773 // final step: refit all using the parameters and tollerance limits (+-20%)\r
774 // WARNING: exponential tail added in the right of the Gaussian shape\r
775 Bool_t mc=kFALSE;\r
776 Int_t code=TMath::Abs(signedcode);\r
777 if(!IsGoodBin(bin,code)) return;\r
778\r
779 TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
780 TString modfit = "M0R+";\r
781 Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
782 GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
783\r
784 printf("\n___________________________________________________________________\n First Step: pions\n\n");\r
785 fstep1 = new TF1("step1",SingleGausTail,-3.5,3.5,5);\r
786 fstep1->SetParameter(0,ampl);//initial \r
787 fstep1->SetParameter(1,mean);\r
788 fstep1->SetParameter(3,1.2);\r
789 fstep1->SetParameter(4,10.);\r
790\r
791 fstep1->SetParLimits(0,0,ampl*1.2);\r
792 fstep1->SetParLimits(1,-3.5,3.5);\r
793 fstep1->SetParLimits(2,0.1,0.25);\r
794 fstep1->SetParLimits(4,5.,20.);\r
795 if(bin<8) fstep1->SetParLimits(4,13.,25.);\r
796\r
797 h->Fit(fstep1,modfit,"",mean-0.45,mean+0.45);//first fit\r
798\r
799 printf("\n___________________________________________________________________\n Second Step: kaons\n\n"); \r
800 fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10);\r
801 fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed\r
802 fstep2->FixParameter(1,fstep1->GetParameter(1));\r
803 fstep2->FixParameter(2,fstep1->GetParameter(2));\r
804 fstep2->FixParameter(3,fstep1->GetParameter(3));\r
805 fstep2->FixParameter(4,fstep1->GetParameter(4));\r
806\r
807 fstep2->SetParameter(5,fstep1->GetParameter(0)/8);//initial\r
808 //fstep2->SetParameter(6,CalcP(code,322,bin));\r
809 fstep2->SetParameter(7,0.145);\r
810 fstep2->FixParameter(8,1.2);\r
811 fstep2->SetParameter(9,13.);\r
812\r
813 fstep2->SetParLimits(5,fstep1->GetParameter(0)/100,fstep1->GetParameter(0));//limits\r
814 fstep2->SetParLimits(6,-3.5,3.5);\r
815 fstep2->SetParLimits(7,0.12,0.2);\r
816 fstep2->SetParLimits(9,9.,20.);\r
817 if(bin<9) fstep2->SetParLimits(9,13.,25.);\r
818\r
819 //h->Fit(fstep2,"M0R+","",CalcP(code,321,bin)-0.3,CalcP(code,321,bin)+0.3);//second fit\r
820 if(bin<6 || bin>12) for(Int_t npar=5;npar<10;npar++) fstep2->FixParameter(npar,-0.0000000001);\r
821\r
822 printf("\n____________________________________________________________________\n Third Step: protons \n\n");\r
823 fstep3= new TF1("fstep3",FinalGausTail,-3.5,3.5,15);\r
824 fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed\r
825 fstep3->FixParameter(1,fstep1->GetParameter(1));\r
826 fstep3->FixParameter(2,fstep1->GetParameter(2));\r
827 fstep3->FixParameter(3,fstep1->GetParameter(3));\r
828 fstep3->FixParameter(4,fstep1->GetParameter(4));\r
829 fstep3->FixParameter(5,fstep2->GetParameter(5));\r
830 fstep3->FixParameter(6,fstep2->GetParameter(6));\r
831 fstep3->FixParameter(7,fstep2->GetParameter(7));\r
832 fstep3->FixParameter(8,fstep2->GetParameter(8));\r
833 fstep3->FixParameter(9,fstep2->GetParameter(9));\r
834\r
835 fstep3->SetParameter(10,fstep2->GetParameter(0)/8);//initial\r
836 //fstep3->SetParameter(11,CalcP(code,2212,bin));\r
837 fstep3->SetParameter(12,0.145);\r
838 fstep3->FixParameter(13,1.2);\r
839 fstep3->SetParameter(14,10.);\r
840\r
841 fstep3->SetParLimits(10,fstep2->GetParameter(0)/100,fstep2->GetParameter(0));//limits\r
842 fstep3->SetParLimits(11,-3.5,3.5);\r
843 fstep3->SetParLimits(12,0.12,0.2);\r
844 fstep3->SetParLimits(14,11.,25.);\r
845\r
846 printf("\n_____________________________________________________________________\n Final Step: refit all \n\n"); \r
847 fstepTot = new TF1("funztot",FinalGausTail,-3.5,3.5,15);\r
848 fstepTot->SetLineColor(1);\r
849\r
850 Double_t initialParametersStepTot[15];\r
851 initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
852 initialParametersStepTot[1] = fstep1->GetParameter(1);\r
853 initialParametersStepTot[2] = fstep1->GetParameter(2);\r
854 initialParametersStepTot[3] = fstep1->GetParameter(3);\r
855 initialParametersStepTot[4] = fstep1->GetParameter(4);\r
856\r
857 initialParametersStepTot[5] = fstep2->GetParameter(5);//second gaussian\r
858 initialParametersStepTot[6] = fstep2->GetParameter(6);\r
859 initialParametersStepTot[7] = fstep2->GetParameter(7);\r
860 initialParametersStepTot[8] = fstep2->GetParameter(8);\r
861 initialParametersStepTot[9] = fstep2->GetParameter(9);\r
862\r
863 initialParametersStepTot[10] = fstep3->GetParameter(10);//third gaussian\r
864 initialParametersStepTot[11] = fstep3->GetParameter(11);\r
865 initialParametersStepTot[12] = fstep3->GetParameter(12); \r
866 initialParametersStepTot[13] = fstep3->GetParameter(13); \r
867 initialParametersStepTot[14] = fstep3->GetParameter(14); \r
868\r
869 fstepTot->SetParameters(initialParametersStepTot);//initial parameter\r
870\r
871\r
872 fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tollerance limit\r
873 fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);\r
874 fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);\r
875 fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);\r
876 fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);\r
877 fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);\r
878 fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);\r
879 fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);\r
880 fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);\r
881 fstepTot->SetParLimits(9,initialParametersStepTot[9]*0.9,initialParametersStepTot[9]*1.1);\r
882 fstepTot->SetParLimits(10,initialParametersStepTot[10]*0.9,initialParametersStepTot[10]*1.1);\r
883 fstepTot->SetParLimits(11,initialParametersStepTot[11]*0.9,initialParametersStepTot[11]*1.1);\r
884 fstepTot->SetParLimits(12,initialParametersStepTot[12]*0.9,initialParametersStepTot[12]*1.1);\r
885 fstepTot->SetParLimits(13,initialParametersStepTot[13]*0.9,initialParametersStepTot[13]*1.1);\r
886 fstepTot->SetParLimits(14,initialParametersStepTot[14]*0.9,initialParametersStepTot[14]*1.1);\r
887\r
888 if(bin<9) for(Int_t npar=10;npar<15;npar++) fstepTot->FixParameter(npar,-0.00000000001);\r
889 h->Fit(fstepTot,modfit,"",-3.5,3.5); //refit all\r
890\r
891\r
892 //************************************* storing parameter to calculate the yields *******\r
893 Int_t chpa=0;\r
894 if(code==321) chpa=5;\r
895 if(code==2212) chpa=10;\r
896 for(Int_t j=0;j<3;j++) {\r
897 fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
898 fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
899 }\r
900\r
901 DrawFitFunction(fstepTot);\r
902 return;\r
fdfb7ea9 903}\r
904\r
905//________________________________________________________\r
906void AliITSsadEdxFitter::FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code){\r
7d629884 907 // fill the spectra histo calculating the yield\r
908 // first bin has to be 1\r
909 Double_t yield = 0;\r
910 Double_t err = 0;\r
911 Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin); \r
912\r
913 if(IsGoodBin(bin-1,code)) {\r
914 yield = fFitpar[0] / ptbin / binsize;\r
915 err = fFitparErr[0] / ptbin / binsize;\r
916 }\r
917\r
918 hsps->SetBinContent(bin,yield);\r
919 hsps->SetBinError(bin,err);\r
920 return;\r
fdfb7ea9 921}\r
922\r
923//________________________________________________________\r
924void AliITSsadEdxFitter::FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h){\r
7d629884 925 // fill the spectra histo calculating the yield (using the MC truth)\r
926 // first bin has to be 1\r
927 Double_t yield = 0;\r
928 Double_t erryield=0;\r
929 Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);\r
930\r
931 if(IsGoodBin(bin-1,code)){\r
932 yield = h->GetEntries() / ptbin;\r
933 erryield=TMath::Sqrt(h->GetEntries()) / ptbin;\r
934 }\r
935 hsps->SetBinContent(bin,yield);\r
936 hsps->SetBinError(bin,erryield);\r
937 return;\r
fdfb7ea9 938}\r
939\r
940//________________________________________________________\r
941void AliITSsadEdxFitter::GetFitPar(Double_t *fitpar, Double_t *fitparerr) const {\r
7d629884 942 // getter of the fit parameters and the relative errors\r
943 for(Int_t i=0;i<3;i++) {\r
944 fitpar[i] = fFitpar[i];\r
945 fitparerr[i] = fFitparErr[i];\r
946 }\r
947 return;\r
fdfb7ea9 948}\r
949\r
950\r
951//________________________________________________________\r
952void AliITSsadEdxFitter::PrintAll() const{\r
953 //\r
7d629884 954 printf("Range 1 = %f %f\n",fRangeStep1[0],fRangeStep1[1]);\r
955 printf("Range 2 = %f %f\n",fRangeStep2[0],fRangeStep2[1]);\r
956 printf("Range 3 = %f %f\n",fRangeStep3[0],fRangeStep3[1]);\r
957 printf("Range F = %f %f\n",fRangeStep4[0],fRangeStep4[1]);\r
958 printf(" Sigma1 = %f %f\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1]);\r
959 printf(" Sigma2 = %f %f\n",fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1]);\r
960 printf(" Sigma3 = %f %f\n",fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);\r
fdfb7ea9 961}\r