]>
Commit | Line | Data |
---|---|---|
d36ecf07 | 1 | /************************************************************************** |
2 | * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | /////////////////////////////////////////////////////////////////////// | |
19 | // Class with the fits algorithms to be used in the identified // | |
20 | // spectra analysis using the ITS in stand-alone mode // | |
21 | // Author: E.Biolcati, biolcati@to.infn.it // | |
22 | // F.Prino, prino@to.infn.it // | |
23 | /////////////////////////////////////////////////////////////////////// | |
24 | ||
25 | #include <Riostream.h> | |
26 | #include <TLatex.h> | |
27 | #include <TH1F.h> | |
28 | #include <TF1.h> | |
29 | #include <TH1D.h> | |
30 | #include <TLine.h> | |
31 | #include <TH2F.h> | |
32 | #include <TMath.h> | |
33 | #include <TGraph.h> | |
34 | #include <TGraphErrors.h> | |
35 | #include <TLegend.h> | |
36 | #include <TLegendEntry.h> | |
37 | #include <TStyle.h> | |
38 | #include <Rtypes.h> | |
39 | #include "AliITSsadEdxFitter.h" | |
40 | ||
41 | ||
42 | ClassImp(AliITSsadEdxFitter) | |
43 | //______________________________________________________________________ | |
44 | AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(), | |
45 | fExpPionMean(0.), | |
46 | fExpKaonMean(0.), | |
47 | fExpProtonMean(0.), | |
48 | fExpPionMeanRange(0.), | |
49 | fExpKaonMeanRange(0.), | |
50 | fExpProtonMeanRange(0.), | |
51 | fExpPionSigma(0.), | |
52 | fExpKaonSigma(0.), | |
53 | fExpProtonSigma(0.), | |
54 | fExpPionAmpl(0.), | |
55 | fIsMC(kFALSE), | |
56 | fOptInit(0), | |
d7cc7766 | 57 | fFitOption("M0R+"), |
d36ecf07 | 58 | fITSpid(0) |
59 | { | |
60 | // standard constructor | |
61 | for(Int_t i=0; i<9; i++) fFitPars[i] = 0.; | |
62 | for(Int_t i=0; i<9; i++) fFitParErrors[i] = 0.; | |
63 | for(Int_t i=0; i<3; i++) fFitpar[i] = 0.; | |
64 | for(Int_t i=0; i<3; i++) fFitparErr[i] = 0.; | |
65 | SetRangeStep1(); | |
66 | SetRangeStep2(); | |
67 | SetRangeStep3(); | |
68 | SetRangeFinalStep(); | |
69 | SetLimitsOnSigmaPion(); | |
70 | SetLimitsOnSigmaKaon(); | |
71 | SetLimitsOnSigmaProton(); | |
72 | SetBinsUsedPion(); | |
73 | SetBinsUsedKaon(); | |
74 | SetBinsUsedProton(); | |
75 | fITSpid=new AliITSPIDResponse(kFALSE); | |
76 | }; | |
77 | //______________________________________________________________________ | |
78 | AliITSsadEdxFitter::AliITSsadEdxFitter(Bool_t isMC):TObject(), | |
79 | fExpPionMean(0.), | |
80 | fExpKaonMean(0.), | |
81 | fExpProtonMean(0.), | |
82 | fExpPionMeanRange(0.), | |
83 | fExpKaonMeanRange(0.), | |
84 | fExpProtonMeanRange(0.), | |
85 | fExpPionSigma(0.), | |
86 | fExpKaonSigma(0.), | |
87 | fExpProtonSigma(0.), | |
88 | fExpPionAmpl(0.), | |
89 | fIsMC(isMC), | |
90 | fOptInit(0), | |
d7cc7766 | 91 | fFitOption("M0R+"), |
d36ecf07 | 92 | fITSpid(0) |
93 | { | |
94 | // standard constructor | |
95 | for(Int_t i=0; i<9; i++) fFitPars[i] = 0.; | |
96 | for(Int_t i=0; i<9; i++) fFitParErrors[i] = 0.; | |
97 | for(Int_t i=0; i<3; i++) fFitpar[i] = 0.; | |
98 | for(Int_t i=0; i<3; i++) fFitparErr[i] = 0.; | |
99 | SetRangeStep1(); | |
100 | SetRangeStep2(); | |
101 | SetRangeStep3(); | |
102 | SetRangeFinalStep(); | |
103 | SetLimitsOnSigmaPion(); | |
104 | SetLimitsOnSigmaKaon(); | |
105 | SetLimitsOnSigmaProton(); | |
106 | SetBinsUsedPion(); | |
107 | SetBinsUsedKaon(); | |
108 | SetBinsUsedProton(); | |
109 | fITSpid=new AliITSPIDResponse(isMC); | |
110 | }; | |
111 | ||
112 | //________________________________________________________ | |
113 | Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x) const { | |
114 | // calculate the sigma 12-ott-2010 | |
115 | Double_t p[2]={0.}; | |
116 | Double_t xMinKaon=0.15; //minimum pt value to consider the kaon peak | |
117 | Double_t xMinProton=0.3;//minimum pt value to consider the proton peak | |
118 | if(fIsMC){ | |
119 | if(code==211){ | |
120 | p[0] = -1.20337e-04; | |
121 | p[1] = 1.13060e-01; | |
122 | } | |
123 | else if(code==321 && x>xMinKaon){ | |
124 | p[0] = -2.39631e-03; | |
125 | p[1] = 1.15723e-01; | |
126 | } | |
127 | else if(code==2212 && x>xMinProton){ | |
128 | p[0] = -8.34576e-03; | |
129 | p[1] = 1.34237e-01; | |
130 | } | |
131 | else return -1; | |
132 | } | |
133 | else{ | |
134 | if(code==211){ | |
135 | p[0] = -6.55200e-05; | |
136 | p[1] = 1.26657e-01; | |
137 | } | |
138 | else if(code==321 && x>xMinKaon){ | |
139 | p[0] = -6.22639e-04; | |
140 | p[1] = 1.43289e-01; | |
141 | } | |
142 | else if(code==2212 && x>xMinProton){ | |
143 | p[0] = -2.13243e-03; | |
144 | p[1] = 1.68614e-01; | |
145 | } | |
146 | else return -1; | |
147 | } | |
148 | return p[0]/(x*x)*TMath::Log(x)+p[1]; | |
149 | } | |
150 | ||
151 | //_______________________________________________________ | |
152 | Int_t AliITSsadEdxFitter::InitializeMeanPosParam(Int_t code, Float_t x){ | |
153 | // Set initial values for peak positions from ad-hoc parameterizations (Melinda, 5 april 2010) | |
154 | ||
155 | Double_t p1[5]={0.}; | |
156 | Double_t p2[5]={0.}; | |
157 | Double_t p3[5]={0.}; | |
158 | Double_t ep1[5]={0.}; | |
159 | Double_t ep2[5]={0.}; | |
160 | Double_t ep3[5]={0.}; | |
161 | if(!fIsMC){ // data parameters | |
162 | if(code==211){ | |
163 | ||
164 | p1[0]=-4.21001e-04; p1[1]=-3.41133e-02; p1[2]=0.; p1[3]=0.; p1[4]=0.; | |
165 | ep1[0]= 3.62105e-06; ep1[1]= 1.52077e-04; ep1[2]= 0.; ep1[3]= 0.; ep1[4]= 0.; | |
166 | ||
167 | p2[0] = 8.04038e+00; p2[1] = 1.51139e-01; p2[2] = 8.89984e-02; p2[3]= 0.; p2[4]= 0.; | |
168 | ep2[0]= 6.47330e-02; ep2[1]= 3.09501e-03; ep2[2]= 6.03642e-04; ep2[3]= 0.; ep2[4]= 0.; | |
169 | ||
170 | p3[0] = 7.41939e+00; p3[1] =-1.31031e+00; p3[2] = 9.44093e-01; p3[3]= 0.; p3[4]= 0.; | |
171 | ep3[0]= 1.29862e-01; ep3[1]= 1.37370e-02; ep3[2]= 4.74995e-03; ep3[3]= 0.; ep3[4]= 0.; | |
172 | ||
173 | fExpPionMean=p1[0]/(x*x)*TMath::Log(x)+p1[1]; | |
174 | fExpPionMeanRange=TMath::Sqrt(TMath::Power(1/(x*x)*TMath::Log(x),2.)*ep1[0]*ep1[0]+ep1[1]*ep1[1] ); | |
175 | fExpKaonMean=p2[0]*TMath::Landau(x,p2[1],p2[2],kFALSE); | |
176 | fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]); | |
177 | fExpProtonMean=-9999.; | |
178 | fExpProtonMeanRange=0.; | |
179 | if(x>0.2){ | |
180 | fExpProtonMean=p3[0]*TMath::Exp(-0.5*TMath::Power((x-p3[1])/p3[2],2.)) ; | |
181 | fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]); | |
182 | } | |
183 | ||
184 | }else if(code==321){ | |
185 | ||
186 | p1[0] = 8.35645e-01; p1[1] = 3.61638e+00; p1[2] =-2.51068e-01; p1[3]= 0.; p1[4]= 0.; | |
187 | ep1[0]= 6.63101e-04; ep1[1]= 2.33140e-02; ep1[2]= 1.93276e-03; ep1[3]= 0.; ep1[4]= 0.; | |
188 | ||
189 | p2[0] = -4.75573e-02; p2[1] = -7.05764e-03; p2[2] = 7.96719e+00; p2[3]= -7.48333e-02; p2[4]= 0.; | |
190 | ep2[0]= 1.70759e-02; ep2[1]= 2.95064e-04; ep2[2]= 2.70493e+00; ep2[3]= 1.14466e-02; ep2[4]= 0.; | |
191 | ||
192 | p3[0] = 4.20188e+00; p3[1] = 3.95855e-01; p3[2] = 2.01159e-01; p3[3]= 0.; p3[4]= 0.; | |
193 | ep3[0]= 2.18608e-02; ep3[1]= 4.15293e-03; ep3[2]= 8.91846e-03; ep3[3]= 0.; ep3[4]= 0.; | |
194 | ||
195 | fExpPionMean= p1[0]*TMath::Log(x)*TMath::Exp(-x*x*p1[1])+p1[2]; | |
196 | fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]); | |
197 | fExpKaonMean = p2[0]*TMath::Log(x)+p2[1]*TMath::Exp(-x*x*p2[2])/(x*x)+p2[3]; | |
198 | fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]); | |
199 | fExpProtonMean = p3[0]*TMath::Landau(x, p3[1], p3[2],kFALSE); | |
200 | fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]); | |
201 | ||
202 | }else if(code==2212){ | |
203 | ||
204 | p1[0] =1.48277e+00; p1[1] =1.14050e+00; p1[2] =1.01751e+02; p1[3]= -2.52768e-01; p1[4]= 0.; | |
205 | ep1[0]= 1.85491e-03; ep1[1]= 4.23200e-02; ep1[2]= 2.66855e+00; ep1[3]= 1.82240e-03; ep1[4]= 0.; | |
206 | ||
207 | p2[0]= 3.80541e-01 ; p2[1]=-3.79158e-02; p2[2]=-4.81653e-01; p2[3]= 0.; p2[4]= 0.; | |
208 | ep2[0]= 1.43643e-01; ep2[1]= 2.47156e-02; ep2[2]= 1.06836e-01; ep2[3]= 0.; ep2[4]= 0.; | |
209 | ||
210 | p3[0] = 1.81168e-01; p3[1] = -6.69364e-01; p3[2] = 2.10685e+01; p3[3]= 5.14868e-02; p3[4]= 0.; | |
211 | ep3[0]= 4.61718e-02; ep3[1]= 8.21949e-02; ep3[2]= 4.46134e+00; ep3[3]= 2.69751e-02; ep3[4]= 0.; | |
212 | ||
213 | fExpPionMean= p1[0]*TMath::Log(x)+p1[1]*TMath::Exp(-x*x*p1[2])+p1[3]; | |
214 | fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]); | |
215 | fExpKaonMean = p2[0]*TMath::Log(x)+p2[1]/(x)+p2[2]; | |
216 | fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]+(x*x*x*x*x*x)*ep2[3]*ep2[3]+(x*x*x*x*x*x*x*x)*ep2[4]*ep2[4]); | |
217 | fExpProtonMean= p3[0]*TMath::Log(x)+p3[1]*TMath::Exp(-x*x*p3[2])+p3[3]; | |
218 | fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]+(x*x*x*x*x*x)*ep3[3]*ep3[3]+(x*x*x*x*x*x*x*x)*ep3[4]*ep3[4]); | |
219 | }else{ | |
220 | printf("ERROR: Wrong particle code %d\n",code); | |
221 | return -1; | |
222 | } | |
223 | }else{ // MC parameters | |
224 | if(code==211){ | |
225 | ||
226 | p1[0]= -3.28744e-04; p1[1]= -1.78123e-02; p1[2]= 0.; p1[3]= 0.; p1[4]= 0.; | |
227 | ep1[0]= 2.03102e-06; ep1[1]= 7.65859e-05; ep1[2]= 0.; ep1[3]= 0.; ep1[4]= 0.; | |
228 | ||
229 | p2[0] = 8.84991e+00 ; p2[1] = 1.32096e-01 ; p2[2] = 9.42394e-02 ; p2[3]=0.; p2[4]=0.; | |
230 | ep2[0]= 4.78123e-02; ep2[1]= 1.99037e-03; ep2[2]= 2.85136e-04; ep2[3]= 0.; ep2[4]= 0.; | |
231 | ||
232 | p3[0] = 6.25025e+00; p3[1] = -1.08315e+00; p3[2] = 8.86474e-01; p3[3]= 0.; p3[4]= 0.; | |
233 | ep3[0]= 1.42379e+00; ep3[1]= 2.34807e-01; ep3[2]= 6.78015e-02; ep3[3]= 0.; ep3[4]= 0.; | |
234 | ||
235 | fExpPionMean= p1[0]/(x*x)*TMath::Log(x)+p1[1]; | |
236 | fExpPionMeanRange=TMath::Sqrt(TMath::Power(1/(x*x)*TMath::Log(x),2.)*ep1[0]*ep1[0]+ep1[1]*ep1[1] ); | |
237 | fExpKaonMean = p2[0]*TMath::Landau(x,p2[1],p2[2],kFALSE); | |
238 | fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]); | |
239 | fExpProtonMean=-9999.; | |
240 | fExpProtonMeanRange=0.; | |
241 | if(x>0.2){ | |
242 | fExpProtonMean = p3[0]*TMath::Exp(-0.5*TMath::Power((x-p3[1])/p3[2],2.)); | |
243 | fExpProtonMeanRange = TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]); | |
244 | } | |
245 | ||
246 | }else if(code==321){ | |
247 | ||
248 | p1[0] = 9.19679e-01; p1[1] = 5.07077e-03; p1[2] = 5.48262e-05; p1[3]= 8.65674e-02; p1[4]= 0.; | |
249 | ep1[0]= 4.03325e-04; ep1[1]= 8.85421e-06; ep1[2]= 1.77437e+00; ep1[3]= 9.00492e-03; ep1[4]= 0.; | |
250 | ||
251 | p2[0] = -2.05127e-01; p2[1] = -2.05303e-03; p2[2] = 3.10090e-06; p2[3]= -2.19936e-01; p2[4]= 0.; | |
252 | ep2[0]= 3.41607e-03; ep2[1]= 1.03661e-04; ep2[2]= 9.50776e-01; ep2[3]= 3.41148e-03; ep2[4]= 0.; | |
253 | ||
254 | p3[0] = 4.68022e+00; p3[1] = 3.02178e-01; p3[2] = 2.47185e-01; p3[3]= 0.; p3[4]= 0.; | |
255 | ep3[0]= 1.63926e-02; ep3[1]= 6.42072e-03; ep3[2]= 7.24079e-03; ep3[3]= 0.; ep3[4]= 0.; | |
256 | ||
257 | fExpPionMean=p1[0]*TMath::Log(x)+p1[1]*TMath::Exp(-x*x*p1[2])/(x*x)+p1[3]; | |
258 | fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]); | |
259 | fExpKaonMean=p2[0]*TMath::Log(x)+p2[1]*TMath::Exp(-x*x*p2[2])/(x*x+p2[3]); | |
260 | fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]); | |
261 | fExpProtonMean= p3[0]*TMath::Landau(x, p3[1], p3[2]); | |
262 | fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]); | |
263 | ||
264 | }else if(code==2212){ | |
265 | ||
266 | p1[0] = 1.05826e+00; p1[1]=8.59088e-01; p1[2]=9.39446e+01; p1[3]= -4.86736e-01; p1[4]= 0.; | |
267 | ep1[0]= 1.07813e-03; ep1[1]= 1.15536e-02; ep1[2]= 9.74843e-01; ep1[3]= 1.24205e-03; ep1[4]= 0.; | |
268 | ||
269 | p2[0]=4.91075e-01; p2[1]=6.27855e-01; p2[2] = 1.09965e+01; p2[3] =-4.24656e-01; p2[4] = 0.; | |
270 | ep2[0]=3.00310e-02; ep2[1]= 4.42650e-02; ep2[2]= 3.40324e-01; ep2[3]= 2.46448e-02; ep2[4]= 0. ; | |
271 | ||
272 | p3[0] = 6.75789e-01; p3[1] = 1.23129e+00; p3[2] = 3.46294e+00; p3[3]= -2.48336e-02; p3[4]= 0.; | |
273 | ep3[0]= 2.90342e-02; ep3[1]= 4.34524e-02; ep3[2]= 7.18743e-02; ep3[3]= 1.08694e-02; ep3[4]= 0.; | |
274 | ||
275 | fExpPionMean= p1[0]*TMath::Log(x)+p1[1]*TMath::Exp(-x*x*p1[2])+p1[3]; | |
276 | fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]); | |
277 | fExpKaonMean = p2[0]*TMath::Log(x)+p2[1]*TMath::Exp(-x*x*p2[2])+p2[3]; | |
278 | fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]+(x*x*x*x*x*x)*ep2[3]*ep2[3]+(x*x*x*x*x*x*x*x)*ep2[4]*ep2[4]); | |
279 | fExpProtonMean= p3[0]*TMath::Log(x)+p3[1]*TMath::Exp(-x*x*p3[2])+p3[3]; | |
280 | fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]+(x*x*x*x*x*x)*ep3[3]*ep3[3]+(x*x*x*x*x*x*x*x)*ep3[4]*ep3[4]); | |
281 | }else{ | |
282 | printf("ERROR: Wrong particle code %d\n",code); | |
283 | return -1; | |
284 | } | |
285 | } | |
286 | return 0; | |
287 | } | |
288 | ||
289 | //_______________________________________________________ | |
290 | Int_t AliITSsadEdxFitter::InitializeMeanPosBB(Int_t code, Float_t pt){ | |
291 | // computes peak positions from parametrized BB | |
292 | ||
293 | Double_t massp=AliPID::ParticleMass(AliPID::kProton); | |
294 | Double_t massk=AliPID::ParticleMass(AliPID::kKaon); | |
295 | Double_t masspi=AliPID::ParticleMass(AliPID::kPion); | |
296 | Bool_t isSA=kTRUE; | |
297 | Float_t sinthmed=0.8878; | |
298 | Float_t p=pt/sinthmed; | |
299 | Double_t bethep=fITSpid->Bethe(p,massp,isSA); | |
300 | Double_t bethek=fITSpid->Bethe(p,massk,isSA); | |
301 | Double_t bethepi=fITSpid->Bethe(p,masspi,isSA); | |
302 | Double_t betheref=bethepi; | |
303 | if(TMath::Abs(code)==321) betheref=bethek; | |
304 | else if(TMath::Abs(code)==2212) betheref=bethep; | |
305 | fExpPionMean=TMath::Log(bethepi)-TMath::Log(betheref); | |
306 | fExpKaonMean=TMath::Log(bethek)-TMath::Log(betheref); | |
307 | fExpProtonMean=TMath::Log(bethep)-TMath::Log(betheref); | |
308 | return 0; | |
309 | } | |
310 | //________________________________________________________ | |
311 | Bool_t AliITSsadEdxFitter::IsGoodBin(Int_t bin,Int_t code){ | |
312 | //method to select the bins used for the analysis only | |
313 | Bool_t retvalue=kTRUE; | |
314 | TLine *l[2]; //for the cross | |
315 | l[0]=new TLine(-2.1,0,2.,100.); | |
316 | l[1]=new TLine(-1.9,120,2.,0.); | |
317 | for(Int_t j=0;j<2;j++){ | |
318 | l[j]->SetLineColor(4); | |
319 | l[j]->SetLineWidth(5); | |
320 | } | |
321 | ||
322 | if(code==211 && (bin<fBinsUsedPion[0] || bin>fBinsUsedPion[1])){ | |
323 | for(Int_t j=0;j<2;j++) l[j]->Draw("same"); | |
324 | retvalue=kFALSE; | |
325 | } | |
326 | if(code==321 && (bin<fBinsUsedKaon[0] || bin>fBinsUsedKaon[1])){ | |
327 | for(Int_t j=0;j<2;j++) l[j]->Draw("same"); | |
328 | retvalue=kFALSE; | |
329 | } | |
330 | if(code==2212 && (bin<fBinsUsedProton[0] || bin>fBinsUsedProton[1])){ | |
331 | for(Int_t j=0;j<2;j++) l[j]->Draw("same"); | |
332 | retvalue=kFALSE; | |
333 | } | |
334 | return retvalue; | |
335 | } | |
336 | ||
337 | //________________________________________________________ | |
338 | Double_t SingleGausTail(const Double_t *x, const Double_t *par){ | |
339 | //single gaussian with exponential tail | |
340 | Double_t s2pi=TMath::Sqrt(2*TMath::Pi()); | |
341 | Double_t xx = x[0]; | |
342 | Double_t mean1 = par[1]; | |
343 | Double_t sigma1 = par[2]; | |
344 | Double_t xNormSquare1 = (xx - mean1) * (xx - mean1); | |
345 | Double_t sigmaSquare1 = sigma1 * sigma1; | |
346 | Double_t xdiv1 = mean1 + par[3] * sigma1; | |
347 | Double_t y1=0.0; | |
348 | if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1); | |
349 | else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3])); | |
350 | return y1; | |
351 | } | |
352 | ||
353 | //________________________________________________________ | |
354 | Double_t DoubleGausTail(const Double_t *x, const Double_t *par){ | |
355 | //sum of two gaussians with exponential tail | |
356 | Double_t s2pi=TMath::Sqrt(2*TMath::Pi()); | |
357 | Double_t xx = x[0]; | |
358 | Double_t mean1 = par[1]; | |
359 | Double_t sigma1 = par[2]; | |
360 | Double_t xNormSquare1 = (xx - mean1) * (xx - mean1); | |
361 | Double_t sigmaSquare1 = sigma1 * sigma1; | |
362 | Double_t xdiv1 = mean1 + par[3] * sigma1; | |
363 | Double_t y1=0.0; | |
364 | if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1); | |
365 | else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3])); | |
366 | ||
367 | Double_t mean2 = par[6]; | |
368 | Double_t sigma2 = par[7]; | |
369 | Double_t xNormSquare2 = (xx - mean2) * (xx - mean2); | |
370 | Double_t sigmaSquare2 = sigma2 * sigma2; | |
371 | Double_t xdiv2 = mean2 + par[8] * sigma2; | |
372 | Double_t y2=0.0; | |
373 | if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2); | |
374 | else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8])); | |
375 | return y1+y2; | |
376 | } | |
377 | ||
378 | //________________________________________________________ | |
379 | Double_t FinalGausTail(const Double_t *x, const Double_t *par){ | |
380 | //sum of three gaussians with exponential tail | |
381 | Double_t s2pi=TMath::Sqrt(2*TMath::Pi()); | |
382 | Double_t xx = x[0]; | |
383 | Double_t mean1 = par[1]; | |
384 | Double_t sigma1 = par[2]; | |
385 | Double_t xNormSquare1 = (xx - mean1) * (xx - mean1); | |
386 | Double_t sigmaSquare1 = sigma1 * sigma1; | |
387 | Double_t xdiv1 = mean1 + par[3] * sigma1; | |
388 | Double_t y1=0.0; | |
389 | if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1); | |
390 | else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3])); | |
391 | ||
392 | Double_t mean2 = par[6]; | |
393 | Double_t sigma2 = par[7]; | |
394 | Double_t xNormSquare2 = (xx - mean2) * (xx - mean2); | |
395 | Double_t sigmaSquare2 = sigma2 * sigma2; | |
396 | Double_t xdiv2 = mean2 + par[8] * sigma2; | |
397 | Double_t y2=0.0; | |
398 | if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2); | |
399 | else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8])); | |
400 | ||
401 | Double_t mean3 = par[11]; | |
402 | Double_t sigma3 = par[12]; | |
403 | Double_t xNormSquare3 = (xx - mean3) * (xx - mean3); | |
404 | Double_t sigmaSquare3 = sigma3 * sigma3; | |
405 | Double_t xdiv3 = mean3 + par[13] * sigma3; | |
406 | Double_t y3=0.0; | |
407 | if(xx < xdiv3) y3 = par[10]/(s2pi*par[12]) * TMath::Exp(-0.5 * xNormSquare3 / sigmaSquare3); | |
408 | else y3 = TMath::Exp(-par[14]*(xx-xdiv3)) * par[10] / (s2pi*par[12]) * TMath::Exp(-0.5*(par[13]*par[13])); | |
409 | return y1+y2+y3; | |
410 | } | |
411 | ||
412 | //______________________________________________________________________ | |
413 | void AliITSsadEdxFitter::CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const{ | |
414 | //code to calculate the difference fit function and histogram point (residual) | |
415 | //to use as goodness test for the fit | |
416 | Int_t ipt=0; | |
417 | Double_t x=0.,yhis=0.,yfun=0.; | |
418 | for(Int_t i=0;i<h->GetNbinsX();i++){ | |
419 | x=(h->GetBinLowEdge(i+2)+h->GetBinLowEdge(i+1))/2; | |
420 | yfun=fun->Eval(x); | |
421 | yhis=h->GetBinContent(i+1); | |
422 | if(yhis>0. && yfun>0.) { | |
423 | gres->SetPoint(ipt,x,(yhis-yfun)/yhis); | |
424 | ipt++; | |
425 | } | |
426 | } | |
427 | return; | |
428 | } | |
429 | ||
430 | ||
431 | //________________________________________________________ | |
432 | Double_t SingleGausStep(const Double_t *x, const Double_t *par){ | |
433 | //single normalized gaussian | |
434 | Double_t s2pi=TMath::Sqrt(2*TMath::Pi()); | |
435 | Double_t xx = x[0]; | |
436 | Double_t mean1 = par[1]; | |
437 | Double_t sigma1 = par[2]; | |
438 | Double_t xNorm1Square = (xx - mean1) * (xx - mean1); | |
439 | Double_t sigma1Square = sigma1 * sigma1; | |
440 | Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNorm1Square / sigma1Square); | |
441 | return step1; | |
442 | } | |
443 | ||
444 | //________________________________________________________ | |
445 | Double_t DoubleGausStep(const Double_t *x, const Double_t *par){ | |
446 | //sum of two normalized gaussians | |
447 | Double_t s2pi=TMath::Sqrt(2*TMath::Pi()); | |
448 | Double_t xx = x[0]; | |
449 | Double_t mean1 = par[1]; | |
450 | Double_t sigma1 = par[2]; | |
451 | Double_t xNorm1Square = (xx - mean1) * (xx - mean1); | |
452 | Double_t sigma1Square = sigma1 * sigma1; | |
453 | Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square ); | |
454 | Double_t mean2 = par[4]; | |
455 | Double_t sigma2 = par[5]; | |
456 | Double_t xNorm2Square = (xx - mean2) * (xx - mean2); | |
457 | Double_t sigma2Square = sigma2 * sigma2; | |
458 | Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square ); | |
459 | return step1+step2; | |
460 | } | |
461 | ||
462 | //________________________________________________________ | |
463 | Double_t FinalGausStep(const Double_t *x, const Double_t *par){ | |
464 | //sum of three normalized gaussians | |
465 | Double_t s2pi=TMath::Sqrt(2*TMath::Pi()); | |
466 | Double_t xx = x[0]; | |
467 | Double_t mean1 = par[1]; | |
468 | Double_t sigma1 = par[2]; | |
469 | Double_t xNorm1Square = (xx - mean1) * (xx - mean1); | |
470 | Double_t sigma1Square = sigma1 * sigma1; | |
471 | Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square ); | |
472 | Double_t mean2 = par[4]; | |
473 | Double_t sigma2 = par[5]; | |
474 | Double_t xNorm2Square = (xx - mean2) * (xx - mean2); | |
475 | Double_t sigma2Square = sigma2 * sigma2; | |
476 | Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square ); | |
477 | Double_t mean3 = par[7]; | |
478 | Double_t sigma3 = par[8]; | |
479 | Double_t xNorm3Square = (xx - mean3) * (xx - mean3); | |
480 | Double_t sigma3Square = sigma3 * sigma3; | |
481 | Double_t step3 = par[6]/(s2pi*par[8]) * TMath::Exp( - 0.5 * xNorm3Square / sigma3Square); | |
482 | return step1+step2+step3; | |
483 | } | |
484 | ||
485 | //______________________________________________________________________ | |
486 | 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{ | |
487 | //gaussian with an exponential tail on the right side | |
488 | Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi())); | |
489 | Double_t returnvalue=0.0; | |
490 | 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); | |
491 | if (x<mean+cut*rms) returnvalue=TMath::Exp(-1.0*(x-mean)*(x-mean)/(2.0*rms*rms))*factor/TMath::Abs(rms); | |
492 | else returnvalue=TMath::Exp(slope*(mean+cut*rms-x))*TMath::Exp(-cut*cut*0.5)*factor/TMath::Abs(rms); | |
493 | return c*returnvalue/n; | |
494 | } | |
495 | ||
496 | ||
497 | //______________________________________________________________________ | |
498 | Double_t AliITSsadEdxFitter::GausOnBackground(const Double_t* x, const Double_t *par) const { | |
499 | //gaussian with a background parametrisation | |
500 | //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<x[0]<< endl; | |
501 | Double_t returnvalue=0.0; | |
502 | Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi())); | |
503 | 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]); | |
504 | else returnvalue+=GausPlusTail(x[0], par[0], par[1],par[2], par[6], 1.2); | |
505 | returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]); | |
506 | return returnvalue; | |
507 | } | |
508 | ||
509 | //______________________________________________________________________ | |
510 | void AliITSsadEdxFitter::DrawFitFunction(TF1 *fun) const { | |
511 | //code to draw the final fit function and the single gaussians used | |
512 | //to extract the yields for the 3 species | |
513 | TF1 *fdraw1=new TF1("fdraw1",SingleGausStep,-3.5,3.5,3); | |
514 | TF1 *fdraw2=new TF1("fdraw2",SingleGausStep,-3.5,3.5,3); | |
515 | TF1 *fdraw3=new TF1("fdraw3",SingleGausStep,-3.5,3.5,3); | |
516 | fdraw1->SetParameter(0,fun->GetParameter(0)); | |
517 | fdraw1->SetParameter(1,fun->GetParameter(1)); | |
518 | fdraw1->SetParameter(2,fun->GetParameter(2)); | |
519 | fdraw2->SetParameter(0,fun->GetParameter(3)); | |
520 | fdraw2->SetParameter(1,fun->GetParameter(4)); | |
521 | fdraw2->SetParameter(2,fun->GetParameter(5)); | |
522 | fdraw3->SetParameter(0,fun->GetParameter(6)); | |
523 | fdraw3->SetParameter(1,fun->GetParameter(7)); | |
524 | fdraw3->SetParameter(2,fun->GetParameter(8)); | |
525 | ||
526 | fdraw1->SetLineColor(6);//color code | |
527 | fdraw2->SetLineColor(2); | |
528 | fdraw3->SetLineColor(4); | |
529 | ||
530 | fdraw1->SetLineStyle(2);//dot lines | |
531 | fdraw2->SetLineStyle(2); | |
532 | fdraw3->SetLineStyle(2); | |
533 | ||
534 | fun->SetLineWidth(3); | |
535 | fdraw1->DrawCopy("same"); | |
536 | fdraw2->DrawCopy("same"); | |
537 | fdraw3->DrawCopy("same"); | |
538 | fun->DrawCopy("same"); | |
539 | ||
540 | TLatex *ltx[3]; | |
541 | for(Int_t j=0;j<3;j++){ | |
542 | ltx[0]=new TLatex(0.13,0.25,"pions"); | |
543 | ltx[1]=new TLatex(0.13,0.20,"kaons"); | |
544 | ltx[2]=new TLatex(0.13,0.15,"protons"); | |
545 | ltx[0]->SetTextColor(6); | |
546 | ltx[1]->SetTextColor(2); | |
547 | ltx[2]->SetTextColor(4); | |
548 | ltx[j]->SetNDC(); | |
549 | ltx[j]->Draw(); | |
550 | } | |
551 | return; | |
552 | } | |
553 | ||
554 | //______________________________________________________________________ | |
555 | void AliITSsadEdxFitter::Initialize(TH1F* h, Int_t code,Int_t bin){ | |
556 | //code to get the expected values to use for fitting | |
557 | ||
558 | 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}; | |
559 | Double_t pt=(xbins[bin+1]+xbins[bin])/2; | |
560 | ||
561 | //draw and label | |
562 | h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],code)); | |
563 | h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}"); | |
564 | h->GetYaxis()->SetTitle("counts"); | |
565 | h->Draw("e"); | |
566 | h->SetFillColor(11); | |
567 | ||
568 | //expected values | |
569 | fExpPionAmpl = h->GetMaximum()/(h->GetRMS()*TMath::Sqrt(2*TMath::Pi())); | |
570 | if(fOptInit==0) InitializeMeanPosBB(code, pt); | |
571 | else if(fOptInit==1) InitializeMeanPosParam(code, pt); | |
572 | fExpPionSigma = CalcSigma(211,pt); //expected sigma values | |
573 | fExpKaonSigma = CalcSigma(321,pt); | |
574 | fExpProtonSigma = CalcSigma(2212,pt); | |
575 | ||
576 | } | |
577 | ||
578 | ||
579 | //________________________________________________________ | |
580 | void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){ | |
581 | // 3-gaussian fit to log(dedx)-log(dedxBB) histogram | |
582 | // pt bin from 0 to 20, code={211,321,2212} | |
583 | // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed | |
584 | // final step: refit all using the parameters and tollerance limits (+-20%) | |
585 | ||
586 | TF1 *fstep1, *fstep2, *fstep3, *fstepTot; | |
d36ecf07 | 587 | Int_t code=TMath::Abs(signedcode); |
588 | Initialize(h,code,bin); | |
589 | PrintInitialValues(); | |
590 | if(!IsGoodBin(bin,code)) return; | |
591 | ||
592 | printf("___________________________________________________________________ First Step: pions\n"); | |
593 | fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3); | |
594 | fstep1->SetParameter(0,fExpPionAmpl); //initial ampl pion | |
595 | fstep1->SetParameter(1,fExpPionMean); //initial mean pion | |
596 | fstep1->SetParameter(2,fExpPionSigma); //initial sigma pion | |
597 | fstep1->SetParLimits(0,0.,fExpPionAmpl*1.2); //limits ampl pion | |
598 | fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]); //limits mean pion (dummy) | |
599 | fstep1->SetParLimits(2,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]); //limits sigma pion | |
600 | ||
d7cc7766 | 601 | if(fExpPionSigma>0) h->Fit(fstep1,fFitOption.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//first fit |
d36ecf07 | 602 | else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001); |
603 | ||
604 | printf("___________________________________________________________________ Second Step: kaons\n"); | |
605 | fstep2 = new TF1("fstep2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6); | |
606 | fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion | |
607 | fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion | |
608 | fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion | |
609 | fstep2->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon | |
610 | fstep2->SetParameter(4,fExpKaonMean); //initial mean kaon | |
611 | fstep2->SetParameter(3,fExpKaonSigma); //initial sigma kaon | |
612 | fstep2->SetParLimits(3,0.,fstep1->GetParameter(0)); //limits ampl kaon | |
613 | fstep2->SetParLimits(4,fstep1->GetParameter(1),fRangeStep4[1]); //limits mean kaon | |
614 | fstep2->SetParLimits(5,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon | |
615 | ||
d7cc7766 | 616 | if(fExpKaonSigma>0) h->Fit(fstep2,fFitOption.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//second fit |
d36ecf07 | 617 | else for(Int_t npar=3;npar<6;npar++) fstep2->FixParameter(npar,0.00001); |
618 | ||
619 | ||
620 | printf("___________________________________________________________________ Third Step: protons\n"); | |
621 | fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9); | |
622 | fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion | |
623 | fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion | |
624 | fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion | |
625 | fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl kaon | |
626 | fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean kaon | |
627 | fstep3->FixParameter(5,fstep2->GetParameter(5)); //fidex sigma kaon | |
628 | fstep3->SetParameter(6,fstep2->GetParameter(0)/16.); //initial ampl proton | |
629 | fstep3->SetParameter(7,fExpProtonMean); //initial mean proton | |
630 | fstep3->SetParameter(8,fExpProtonSigma); //initial sigma proton | |
631 | fstep3->SetParLimits(6,0.,fstep2->GetParameter(0)); //limits ampl proton | |
632 | fstep3->SetParLimits(7,fstep2->GetParameter(4),fRangeStep4[1]); //limits mean proton | |
633 | fstep3->SetParLimits(8,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton | |
634 | ||
d7cc7766 | 635 | if(fExpProtonSigma>0) h->Fit(fstep3,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//third fit |
d36ecf07 | 636 | else for(Int_t npar=6;npar<9;npar++) fstep3->FixParameter(npar,0.00001); |
637 | ||
638 | printf("___________________________________________________________________ Final Step: refit all\n"); | |
639 | fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9); | |
640 | fstepTot->SetLineColor(1); | |
641 | ||
642 | Double_t initialParametersStepTot[9]; | |
643 | initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian | |
644 | initialParametersStepTot[1] = fstep1->GetParameter(1); | |
645 | initialParametersStepTot[2] = fstep1->GetParameter(2); | |
646 | ||
647 | initialParametersStepTot[3] = fstep2->GetParameter(3);//second gaussian | |
648 | initialParametersStepTot[4] = fstep2->GetParameter(4); | |
649 | initialParametersStepTot[5] = fstep2->GetParameter(5); | |
650 | ||
651 | initialParametersStepTot[6] = fstep3->GetParameter(6);//third gaussian | |
652 | initialParametersStepTot[7] = fstep3->GetParameter(7); | |
653 | initialParametersStepTot[8] = fstep3->GetParameter(8); | |
654 | ||
655 | fstepTot->SetParameters(initialParametersStepTot); //initial parameters | |
656 | fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion | |
657 | fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion | |
658 | fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion | |
659 | fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon | |
660 | fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon | |
661 | fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon | |
662 | fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton | |
663 | fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton | |
664 | fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton | |
665 | ||
d7cc7766 | 666 | h->Fit(fstepTot,fFitOption.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all |
d36ecf07 | 667 | |
668 | for(Int_t j=0;j<9;j++) { | |
669 | fFitPars[j] = fstepTot->GetParameter(j); | |
670 | fFitParErrors[j] = fstepTot->GetParError(j); | |
671 | } | |
672 | ||
673 | //************************************* storing parameter to calculate the yields ******* | |
674 | ||
675 | Int_t chpa=0; | |
676 | if(code==321) chpa=3; | |
677 | if(code==2212) chpa=6; | |
678 | for(Int_t j=0;j<3;j++) { | |
679 | fFitpar[j] = fstepTot->GetParameter(j+chpa); | |
680 | fFitparErr[j] = fstepTot->GetParError(j+chpa); | |
681 | } | |
682 | ||
683 | DrawFitFunction(fstepTot); | |
684 | CalcResidual(h,fstepTot,gres); | |
685 | return; | |
686 | } | |
687 | ||
688 | //________________________________________________________ | |
689 | void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){ | |
690 | // 3-gaussian fit to log(dedx)-log(dedxBB) histogram | |
691 | // pt bin from 0 to 20, code={211,321,2212} | |
692 | // first step: pion peak, second step: proton peak, third step: kaon peak | |
693 | // final step: refit all using the parameters | |
694 | TF1 *fstep1, *fstep2, *fstep3, *fstepTot; | |
d36ecf07 | 695 | Int_t code=TMath::Abs(signedcode); |
696 | Initialize(h,code,bin); | |
697 | PrintInitialValues(); | |
698 | if(!IsGoodBin(bin,code)) return; | |
699 | ||
700 | printf("___________________________________________________________________ First Step: pion\n"); | |
701 | fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3); | |
702 | fstep1->SetParameter(0,fExpPionAmpl); //initial ampl pion | |
703 | fstep1->SetParameter(1,fExpPionMean); //initial mean pion | |
704 | fstep1->SetParameter(2,fExpPionSigma); //initial sigma pion | |
705 | fstep1->SetParLimits(0,0,fExpPionAmpl*1.2); //limits ampl pion | |
706 | fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]); //limits mean pion (dummy) | |
707 | fstep1->SetParLimits(2,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]); //limits sigma pion | |
708 | ||
d7cc7766 | 709 | if(fExpPionSigma>0) h->Fit(fstep1,fFitOption.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//first fit |
d36ecf07 | 710 | else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001); |
711 | ||
712 | printf("___________________________________________________________________ Second Step: proton\n"); | |
713 | fstep2 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3); | |
714 | fstep2->SetParameter(0,fstep1->GetParameter(0)/16.);//initial ampl proton | |
715 | fstep2->SetParameter(1,fExpProtonMean); //initial mean proton | |
716 | fstep2->SetParameter(2,fExpProtonSigma); //initial sigma proton | |
717 | fstep2->SetParLimits(0,0.,fstep1->GetParameter(0)); //limits ampl proton | |
718 | fstep2->SetParLimits(1,fstep1->GetParameter(1),fRangeStep4[1]); //limits mean proton | |
719 | fstep2->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton | |
720 | ||
d7cc7766 | 721 | if(fExpProtonSigma>0) h->Fit(fstep2,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//second fit |
d36ecf07 | 722 | else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001); |
723 | ||
724 | printf("___________________________________________________________________ Third Step: kaon\n"); | |
725 | fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9); | |
726 | fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion | |
727 | fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion | |
728 | fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion | |
729 | fstep3->FixParameter(6,fstep2->GetParameter(0)); //fixed ampl proton | |
730 | fstep3->FixParameter(7,fstep2->GetParameter(1)); //fixed mean proton | |
731 | fstep3->FixParameter(8,fstep2->GetParameter(2)); //fixed sigma proton | |
732 | fstep3->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon | |
733 | fstep3->SetParameter(4,fExpKaonMean); //initial mean kaon | |
734 | fstep3->SetParameter(5,fExpKaonSigma); //initial sigma kaon | |
735 | fstep3->SetParLimits(3,fstep2->GetParameter(0),fstep1->GetParameter(0)); //limits ampl kaon | |
736 | fstep3->SetParLimits(4,fstep1->GetParameter(1),fstep2->GetParameter(1)); //limits mean kaon | |
737 | fstep3->SetParLimits(5,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon | |
d7cc7766 | 738 | if(fExpKaonSigma>0) h->Fit(fstep3,fFitOption.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//third fit |
d36ecf07 | 739 | else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001); |
740 | ||
741 | printf("___________________________________________________________________ Final Step: refit all\n"); | |
742 | fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9); | |
743 | fstepTot->SetLineColor(1); | |
744 | ||
745 | Double_t initialParametersStepTot[9]; | |
746 | initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian | |
747 | initialParametersStepTot[1] = fstep1->GetParameter(1); | |
748 | initialParametersStepTot[2] = fstep1->GetParameter(2); | |
749 | ||
750 | initialParametersStepTot[3] = fstep3->GetParameter(3);//second gaussian | |
751 | initialParametersStepTot[4] = fstep3->GetParameter(4); | |
752 | initialParametersStepTot[5] = fstep3->GetParameter(5); | |
753 | ||
754 | initialParametersStepTot[6] = fstep2->GetParameter(0);//third gaussian | |
755 | initialParametersStepTot[7] = fstep2->GetParameter(1); | |
756 | initialParametersStepTot[8] = fstep2->GetParameter(2); | |
757 | ||
758 | fstepTot->SetParameters(initialParametersStepTot); //initial parameters | |
759 | fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion | |
760 | fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion | |
761 | fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion | |
762 | fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon | |
763 | fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon | |
764 | fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon | |
765 | fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton | |
766 | fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton | |
767 | fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton | |
768 | ||
d7cc7766 | 769 | h->Fit(fstepTot,fFitOption.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all |
d36ecf07 | 770 | |
771 | for(Int_t j=0;j<9;j++) { | |
772 | fFitPars[j] = fstepTot->GetParameter(j); | |
773 | fFitParErrors[j] = fstepTot->GetParError(j); | |
774 | } | |
775 | //************************************* storing parameter to calculate the yields ******* | |
776 | Int_t chpa=0; | |
777 | if(code==321) chpa=3; | |
778 | if(code==2212) chpa=6; | |
779 | for(Int_t j=0;j<3;j++) { | |
780 | fFitpar[j] = fstepTot->GetParameter(j+chpa); | |
781 | fFitparErr[j] = fstepTot->GetParError(j+chpa); | |
782 | } | |
783 | ||
784 | DrawFitFunction(fstepTot); | |
785 | CalcResidual(h,fstepTot,gres); | |
786 | return; | |
787 | } | |
788 | ||
789 | //________________________________________________________ | |
790 | void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){ | |
791 | // 3-gaussian fit to log(dedx)-log(dedxBB) histogram | |
792 | // pt bin from 0 to 20, code={211,321,2212} | |
793 | // first step: proton peak, second step: pion peak, third step: kaon peak | |
794 | // final step: refit all using the parameters | |
795 | TF1 *fstep1, *fstep2, *fstep3, *fstepTot; | |
d36ecf07 | 796 | Int_t code=TMath::Abs(signedcode); |
797 | Initialize(h,code,bin); | |
798 | PrintInitialValues(); | |
799 | if(!IsGoodBin(bin,code)) return; | |
800 | ||
801 | printf("___________________________________________________________________ First Step: proton\n"); | |
802 | fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3); | |
803 | fstep1->SetParameter(0,fExpPionAmpl/16.); //initial ampl proton` | |
804 | fstep1->SetParameter(1,fExpProtonMean); //initial mean proton | |
805 | fstep1->SetParameter(2,fExpProtonSigma); //initial sigma proton | |
806 | fstep1->SetParLimits(0,0,fExpPionAmpl); //limits ampl proton | |
807 | fstep1->SetParLimits(1,fExpPionMean,fRangeStep4[1]); //limits mean proton (dummy) | |
808 | fstep1->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton | |
809 | ||
d7cc7766 | 810 | if(fExpProtonSigma>0) h->Fit(fstep1,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//first fit |
d36ecf07 | 811 | else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001); |
812 | ||
813 | printf("___________________________________________________________________ Second Step: pion\n"); | |
814 | fstep2 = new TF1("step2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6); | |
815 | fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton | |
816 | fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton | |
817 | fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton | |
818 | fstep2->SetParameter(3,fExpPionAmpl); //initial ampl pion | |
819 | fstep2->SetParameter(4,fExpPionMean); //initial mean pion | |
820 | fstep2->SetParameter(5,fExpPionSigma); //initial sigma pion | |
821 | fstep2->SetParLimits(3,0.,1.2*fExpPionAmpl); //limits ampl pion | |
822 | fstep2->SetParLimits(4,fRangeStep4[0],fstep1->GetParameter(1)); //limits mean pion | |
823 | fstep2->SetParLimits(5,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]); //limits sigma pion | |
824 | ||
d7cc7766 | 825 | if(fExpPionSigma>0) h->Fit(fstep2,fFitOption.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//second fit |
d36ecf07 | 826 | else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001); |
827 | ||
828 | printf("___________________________________________________________________ Third Step: kaon\n"); | |
829 | fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9); | |
830 | fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton | |
831 | fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton | |
832 | fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton | |
833 | fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl pion | |
834 | fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean pion | |
835 | fstep3->FixParameter(5,fstep2->GetParameter(5)); //fixed sigma pion | |
836 | fstep3->SetParameter(6,fstep2->GetParameter(0)/8.); //initial ampl kaon | |
837 | fstep3->SetParameter(7,fExpKaonMean); //initial mean kaon | |
838 | fstep3->SetParameter(8,fExpKaonSigma); //initial sigma kaon | |
839 | fstep3->SetParLimits(6,fstep1->GetParameter(0),fstep2->GetParameter(3)); //limits ampl kaon | |
840 | fstep3->SetParLimits(7,fstep2->GetParameter(4),fstep1->GetParameter(1)); //limits mean kaon | |
841 | fstep3->SetParLimits(8,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon | |
d7cc7766 | 842 | if(fExpKaonSigma>0) h->Fit(fstep3,fFitOption.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//third fit |
d36ecf07 | 843 | else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001); |
844 | ||
845 | printf("___________________________________________________________________ Final Step: refit all\n"); | |
846 | fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9); | |
847 | fstepTot->SetLineColor(1); | |
848 | ||
849 | Double_t initialParametersStepTot[9]; | |
850 | initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian | |
851 | initialParametersStepTot[1] = fstep1->GetParameter(1); | |
852 | initialParametersStepTot[2] = fstep1->GetParameter(2); | |
853 | ||
854 | initialParametersStepTot[3] = fstep3->GetParameter(6);//second gaussian | |
855 | initialParametersStepTot[4] = fstep3->GetParameter(7); | |
856 | initialParametersStepTot[5] = fstep3->GetParameter(8); | |
857 | ||
858 | initialParametersStepTot[6] = fstep2->GetParameter(3);//third gaussian | |
859 | initialParametersStepTot[7] = fstep2->GetParameter(4); | |
860 | initialParametersStepTot[8] = fstep2->GetParameter(5); | |
861 | ||
862 | fstepTot->SetParameters(initialParametersStepTot); //initial parameters | |
863 | fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl proton | |
864 | fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion | |
865 | fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion | |
866 | fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon | |
867 | fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon | |
868 | fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon | |
869 | fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl pion | |
870 | fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton | |
871 | fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton | |
872 | ||
d7cc7766 | 873 | h->Fit(fstepTot,fFitOption.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all |
d36ecf07 | 874 | |
875 | for(Int_t j=0;j<9;j++) { | |
876 | fFitPars[j] = fstepTot->GetParameter(j); | |
877 | fFitParErrors[j] = fstepTot->GetParError(j); | |
878 | } | |
879 | //************************************* storing parameter to calculate the yields ******* | |
880 | Int_t chpa=0; | |
881 | if(code==321) chpa=3; | |
882 | if(code==211) chpa=6; | |
883 | for(Int_t j=0;j<3;j++) { | |
884 | fFitpar[j] = fstepTot->GetParameter(j+chpa); | |
885 | fFitparErr[j] = fstepTot->GetParError(j+chpa); | |
886 | } | |
887 | ||
888 | DrawFitFunction(fstepTot); | |
889 | CalcResidual(h,fstepTot,gres); | |
890 | return; | |
891 | } | |
892 | ||
893 | ||
894 | //________________________________________________________ | |
895 | void AliITSsadEdxFitter::DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode){ | |
896 | // single-gaussian fit to log(dedx)-log(dedxBB) histogram | |
897 | TF1 *fstep1; | |
d36ecf07 | 898 | Int_t code=TMath::Abs(signedcode); |
899 | Initialize(h,code,bin); | |
900 | PrintInitialValues(); | |
901 | if(!IsGoodBin(bin,code)) return; | |
902 | ||
903 | printf("___________________________________________________________________ Single Step\n"); | |
904 | fstep1 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3); | |
905 | fstep1->SetParameter(0,fExpPionAmpl/16.); //initial ampl | |
906 | fstep1->SetParameter(1,fExpProtonMean); //initial mean | |
907 | fstep1->SetParameter(2,fExpProtonSigma); //initial sigma | |
908 | fstep1->SetParLimits(0,0.,fExpPionAmpl); //limits ampl proton | |
909 | fstep1->SetParLimits(1,fExpPionMean,fRangeStep4[1]); //limits mean proton | |
910 | //fstep1->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton | |
911 | ||
d7cc7766 | 912 | if(fExpProtonSigma>0) h->Fit(fstep1,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//fit |
d36ecf07 | 913 | else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001); |
914 | ||
915 | fstep1->SetLineColor(1); | |
916 | fstep1->Draw("same"); | |
917 | ||
918 | fFitpar[0] = fstep1->GetParameter(0); | |
919 | fFitparErr[0] = fstep1->GetParError(0); | |
920 | return; | |
921 | } | |
922 | ||
923 | //________________________________________________________ | |
924 | void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t signedcode){ | |
925 | // 3-gaussian fit to log(dedx)-log(dedxBB) histogram | |
926 | // pt bin from 0 to 20, code={211,321,2212} | |
927 | // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed | |
928 | // final step: refit all using the parameters and tollerance limits (+-20%) | |
929 | // WARNING: exponential tail added in the right of the Gaussian shape | |
930 | Int_t code=TMath::Abs(signedcode); | |
931 | if(!IsGoodBin(bin,code)) return; | |
932 | ||
933 | TF1 *fstep1, *fstep2, *fstep3, *fstepTot; | |
d36ecf07 | 934 | Initialize(h,code,bin); |
935 | PrintInitialValues(); | |
936 | ||
937 | printf("\n___________________________________________________________________\n First Step: pions\n\n"); | |
938 | fstep1 = new TF1("step1",SingleGausTail,-3.5,3.5,5); | |
939 | fstep1->SetParameter(0,fExpPionAmpl);//initial | |
940 | fstep1->SetParameter(1,fExpPionMean); | |
941 | fstep1->SetParameter(3,1.2); | |
942 | fstep1->SetParameter(4,10.); | |
943 | ||
944 | fstep1->SetParLimits(0,0,fExpPionAmpl*1.2); | |
945 | fstep1->SetParLimits(1,-3.5,3.5); | |
946 | fstep1->SetParLimits(2,0.1,0.25); | |
947 | fstep1->SetParLimits(4,5.,20.); | |
948 | if(bin<8) fstep1->SetParLimits(4,13.,25.); | |
949 | ||
d7cc7766 | 950 | h->Fit(fstep1,fFitOption.Data(),"",fExpPionMean-0.45,fExpPionMean+0.45);//first fit |
d36ecf07 | 951 | |
952 | printf("\n___________________________________________________________________\n Second Step: kaons\n\n"); | |
953 | fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10); | |
954 | fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed | |
955 | fstep2->FixParameter(1,fstep1->GetParameter(1)); | |
956 | fstep2->FixParameter(2,fstep1->GetParameter(2)); | |
957 | fstep2->FixParameter(3,fstep1->GetParameter(3)); | |
958 | fstep2->FixParameter(4,fstep1->GetParameter(4)); | |
959 | ||
960 | fstep2->SetParameter(5,fstep1->GetParameter(0)/8);//initial | |
961 | //fstep2->SetParameter(6,CalcP(code,322,bin)); | |
962 | fstep2->SetParameter(7,0.145); | |
963 | fstep2->FixParameter(8,1.2); | |
964 | fstep2->SetParameter(9,13.); | |
965 | ||
966 | fstep2->SetParLimits(5,fstep1->GetParameter(0)/100,fstep1->GetParameter(0));//limits | |
967 | fstep2->SetParLimits(6,-3.5,3.5); | |
968 | fstep2->SetParLimits(7,0.12,0.2); | |
969 | fstep2->SetParLimits(9,9.,20.); | |
970 | if(bin<9) fstep2->SetParLimits(9,13.,25.); | |
971 | ||
d36ecf07 | 972 | if(bin<6 || bin>12) for(Int_t npar=5;npar<10;npar++) fstep2->FixParameter(npar,-0.0000000001); |
973 | ||
974 | printf("\n____________________________________________________________________\n Third Step: protons \n\n"); | |
975 | fstep3= new TF1("fstep3",FinalGausTail,-3.5,3.5,15); | |
976 | fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed | |
977 | fstep3->FixParameter(1,fstep1->GetParameter(1)); | |
978 | fstep3->FixParameter(2,fstep1->GetParameter(2)); | |
979 | fstep3->FixParameter(3,fstep1->GetParameter(3)); | |
980 | fstep3->FixParameter(4,fstep1->GetParameter(4)); | |
981 | fstep3->FixParameter(5,fstep2->GetParameter(5)); | |
982 | fstep3->FixParameter(6,fstep2->GetParameter(6)); | |
983 | fstep3->FixParameter(7,fstep2->GetParameter(7)); | |
984 | fstep3->FixParameter(8,fstep2->GetParameter(8)); | |
985 | fstep3->FixParameter(9,fstep2->GetParameter(9)); | |
986 | ||
987 | fstep3->SetParameter(10,fstep2->GetParameter(0)/8);//initial | |
988 | //fstep3->SetParameter(11,CalcP(code,2212,bin)); | |
989 | fstep3->SetParameter(12,0.145); | |
990 | fstep3->FixParameter(13,1.2); | |
991 | fstep3->SetParameter(14,10.); | |
992 | ||
993 | fstep3->SetParLimits(10,fstep2->GetParameter(0)/100,fstep2->GetParameter(0));//limits | |
994 | fstep3->SetParLimits(11,-3.5,3.5); | |
995 | fstep3->SetParLimits(12,0.12,0.2); | |
996 | fstep3->SetParLimits(14,11.,25.); | |
997 | ||
998 | printf("\n_____________________________________________________________________\n Final Step: refit all \n\n"); | |
999 | fstepTot = new TF1("funztot",FinalGausTail,-3.5,3.5,15); | |
1000 | fstepTot->SetLineColor(1); | |
1001 | ||
1002 | Double_t initialParametersStepTot[15]; | |
1003 | initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian | |
1004 | initialParametersStepTot[1] = fstep1->GetParameter(1); | |
1005 | initialParametersStepTot[2] = fstep1->GetParameter(2); | |
1006 | initialParametersStepTot[3] = fstep1->GetParameter(3); | |
1007 | initialParametersStepTot[4] = fstep1->GetParameter(4); | |
1008 | ||
1009 | initialParametersStepTot[5] = fstep2->GetParameter(5);//second gaussian | |
1010 | initialParametersStepTot[6] = fstep2->GetParameter(6); | |
1011 | initialParametersStepTot[7] = fstep2->GetParameter(7); | |
1012 | initialParametersStepTot[8] = fstep2->GetParameter(8); | |
1013 | initialParametersStepTot[9] = fstep2->GetParameter(9); | |
1014 | ||
1015 | initialParametersStepTot[10] = fstep3->GetParameter(10);//third gaussian | |
1016 | initialParametersStepTot[11] = fstep3->GetParameter(11); | |
1017 | initialParametersStepTot[12] = fstep3->GetParameter(12); | |
1018 | initialParametersStepTot[13] = fstep3->GetParameter(13); | |
1019 | initialParametersStepTot[14] = fstep3->GetParameter(14); | |
1020 | ||
1021 | fstepTot->SetParameters(initialParametersStepTot);//initial parameter | |
1022 | ||
1023 | ||
1024 | fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tollerance limit | |
1025 | fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1); | |
1026 | fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1); | |
1027 | fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); | |
1028 | fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1); | |
1029 | fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1); | |
1030 | fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); | |
1031 | fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1); | |
1032 | fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1); | |
1033 | fstepTot->SetParLimits(9,initialParametersStepTot[9]*0.9,initialParametersStepTot[9]*1.1); | |
1034 | fstepTot->SetParLimits(10,initialParametersStepTot[10]*0.9,initialParametersStepTot[10]*1.1); | |
1035 | fstepTot->SetParLimits(11,initialParametersStepTot[11]*0.9,initialParametersStepTot[11]*1.1); | |
1036 | fstepTot->SetParLimits(12,initialParametersStepTot[12]*0.9,initialParametersStepTot[12]*1.1); | |
1037 | fstepTot->SetParLimits(13,initialParametersStepTot[13]*0.9,initialParametersStepTot[13]*1.1); | |
1038 | fstepTot->SetParLimits(14,initialParametersStepTot[14]*0.9,initialParametersStepTot[14]*1.1); | |
1039 | ||
1040 | if(bin<9) for(Int_t npar=10;npar<15;npar++) fstepTot->FixParameter(npar,-0.00000000001); | |
d7cc7766 | 1041 | h->Fit(fstepTot,fFitOption.Data(),"",-3.5,3.5); //refit all |
d36ecf07 | 1042 | |
1043 | for(Int_t j=0;j<9;j++) { | |
1044 | fFitPars[j] = fstepTot->GetParameter(j); | |
1045 | fFitParErrors[j] = fstepTot->GetParError(j); | |
1046 | } | |
1047 | ||
1048 | //************************************* storing parameter to calculate the yields ******* | |
1049 | Int_t chpa=0; | |
1050 | if(code==321) chpa=5; | |
1051 | if(code==2212) chpa=10; | |
1052 | for(Int_t j=0;j<3;j++) { | |
1053 | fFitpar[j] = fstepTot->GetParameter(j+chpa); | |
1054 | fFitparErr[j] = fstepTot->GetParError(j+chpa); | |
1055 | } | |
1056 | ||
1057 | DrawFitFunction(fstepTot); | |
1058 | return; | |
1059 | } | |
1060 | ||
1061 | //________________________________________________________ | |
1062 | void AliITSsadEdxFitter::FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code){ | |
1063 | // fill the spectra histo calculating the yield | |
1064 | // first bin has to be 1 | |
1065 | Double_t yield = 0; | |
1066 | Double_t err = 0; | |
1067 | Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin); | |
1068 | ||
1069 | if(IsGoodBin(bin-1,code)) { | |
1070 | yield = fFitpar[0] / ptbin / binsize; | |
1071 | err = fFitparErr[0] / ptbin / binsize; | |
1072 | } | |
1073 | ||
1074 | hsps->SetBinContent(bin,yield); | |
1075 | hsps->SetBinError(bin,err); | |
1076 | return; | |
1077 | } | |
1078 | ||
1079 | //________________________________________________________ | |
1080 | void AliITSsadEdxFitter::FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h){ | |
1081 | // fill the spectra histo calculating the yield (using the MC truth) | |
1082 | // first bin has to be 1 | |
1083 | Double_t yield = 0; | |
1084 | Double_t erryield=0; | |
1085 | Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin); | |
1086 | ||
1087 | if(IsGoodBin(bin-1,code)){ | |
1088 | yield = h->GetEntries() / ptbin; | |
1089 | erryield=TMath::Sqrt(h->GetEntries()) / ptbin; | |
1090 | } | |
1091 | hsps->SetBinContent(bin,yield); | |
1092 | hsps->SetBinError(bin,erryield); | |
1093 | return; | |
1094 | } | |
1095 | ||
1096 | //________________________________________________________ | |
1097 | void AliITSsadEdxFitter::GetFitPar(Double_t *fitpar, Double_t *fitparerr) const { | |
1098 | // getter of the fit parameters and the relative errors | |
1099 | for(Int_t i=0;i<3;i++) { | |
1100 | fitpar[i] = fFitpar[i]; | |
1101 | fitparerr[i] = fFitparErr[i]; | |
1102 | } | |
1103 | return; | |
1104 | } | |
1105 | ||
1106 | ||
1107 | //________________________________________________________ | |
1108 | void AliITSsadEdxFitter::PrintAll() const{ | |
1109 | // print parameters | |
1110 | ||
1111 | printf("Range 1 = %f %f\n",fRangeStep1[0],fRangeStep1[1]); | |
1112 | printf("Range 2 = %f %f\n",fRangeStep2[0],fRangeStep2[1]); | |
1113 | printf("Range 3 = %f %f\n",fRangeStep3[0],fRangeStep3[1]); | |
1114 | printf("Range F = %f %f\n",fRangeStep4[0],fRangeStep4[1]); | |
1115 | printf(" Sigma1 = %f %f\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1]); | |
1116 | printf(" Sigma2 = %f %f\n",fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1]); | |
1117 | printf(" Sigma3 = %f %f\n",fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]); | |
1118 | } | |
1119 | //______________________________________________________________________ | |
1120 | void AliITSsadEdxFitter::PrintInitialValues() const{ | |
1121 | // print initial values | |
1122 | ||
1123 | printf("mean values -> %f %f %f\n",fExpPionMean,fExpKaonMean,fExpProtonMean); | |
1124 | printf("integration ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n", | |
1125 | fRangeStep1[0],fRangeStep1[1],fRangeStep2[0],fRangeStep2[1],fRangeStep3[0],fRangeStep3[1]); | |
1126 | printf("sigma values -> %f %f %f\n",fExpPionSigma,fExpKaonSigma,fExpProtonSigma); | |
1127 | printf("sigma ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n", | |
1128 | fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1], | |
1129 | fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1], | |
1130 | fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]); | |
1131 | } |