]>
Commit | Line | Data |
---|---|---|
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 | |
42 | ClassImp(AliITSsadEdxFitter)\r | |
7d629884 | 43 | //______________________________________________________________________\r |
44 | AliITSsadEdxFitter::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 | |
61 | Double_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 | |
100 | Int_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 | |
133 | Bool_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 | |
160 | Double_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 | |
176 | Double_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 | |
201 | Double_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 | |
235 | void 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 | |
254 | Double_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 | |
267 | Double_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 | |
285 | Double_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 | |
308 | Double_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 | |
320 | Double_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 | |
332 | void 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 | |
377 | void AliITSsadEdxFitter::GetInitialParam(TH1F* h,Bool_t mc,Int_t code,Int_t bin, Float_t &pt, Float_t &l, 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 | |
409 | void 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 | |
521 | void 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 | |
628 | void 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 | |
739 | void 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 | |
769 | void 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 | |
906 | void 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 | |
924 | void 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 | |
941 | void 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 | |
952 | void 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 |