]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/totEt/macros/hadEt/CorrNeutralLevyFit.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Opening more the TDC check window
[u/mrichter/AliRoot.git] / PWG4 / totEt / macros / hadEt / CorrNeutralLevyFit.C
CommitLineData
bc92a3aa 1//Christine Nattrass, University of Tennessee at Knoxville
2//This macro is to calculate the contributions to Et from various particles based on Levy fits to ALICE data at 900 GeV
3//It uses d^2N/dpTdy from the papers, weighs it by Et, and integrates over all pT.
4//A=0 for mesons
5//A=1 for antinucleons since they would anihilate in the calorimeter if they were observed by a calorimeter
6//A=-1 for nucleons since they would not anihilate and their mass would not be measured
7//At the end this is used to calculate the neutral energy correction
8class AliAnalysisLevyPtModified{
9public:
10 virtual ~AliAnalysisLevyPtModified(){;};
11
12
13 Double_t Evaluate(Double_t *pt, Double_t *par)
14 {
15 Double_t ldNdy = par[0];
16 Double_t l2pi = 2*TMath::Pi();
17 Double_t lTemp = par[1];
18 Double_t lPower = par[2];
19 Double_t lMass = par[3];
20 Double_t A = par[4];
21 Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)+A*lMass;
22
23 Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
24 Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
25
26 return ldNdy *Et* pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
27 }
28 ClassDef(AliAnalysisLevyPtModified, 1);
29};
df8b9b1f 30class AliAnalysisLevyPtModifiedBaryonEnhanced{
31public:
32 virtual ~AliAnalysisLevyPtModifiedBaryonEnhanced(){;};
33
34
35 Double_t Evaluate(Double_t *pt, Double_t *par)
36 {
37 Double_t ldNdy = par[0];
38 Double_t l2pi = 2*TMath::Pi();
39 Double_t lTemp = par[1];
40 Double_t lPower = par[2];
41 Double_t lMass = par[3];
42 Double_t A = par[4];
43 Double_t lMassEt = par[3];//only used for calculating Et
44 Double_t lBary0 = par[6];
45 Double_t lBary1 = par[7];
46 Double_t lBary2 = par[8];
47 Double_t lBary3 = par[9];
48 Double_t lBary4 = par[10];
49 Double_t lBary5 = par[11];
50 //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
51 Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
52
53 Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
54 Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
55 //this is the baryon enhancement factor
56 Double_t baryonEnhancement = lBary0*pow(pt[0],lBary1)*exp(-pow(pt[0]/lBary2,lBary3))/(lBary4+lBary5*pt[0]);
57 if(baryonEnhancement<1.0) baryonEnhancement = 1.0;
58 //This is the density of particles times the Et weight
59 return ldNdy *Et* baryonEnhancement*pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower); }
60 ClassDef(AliAnalysisLevyPtModifiedBaryonEnhanced, 1);
61};
62
63class AliAnalysisLevyPtModifiedStrangeness{
64public:
65 virtual ~AliAnalysisLevyPtModifiedStrangeness(){;};
66
67
68 Double_t Evaluate(Double_t *pt, Double_t *par)
69 {
70 Double_t ldNdy = par[0];
71 Double_t l2pi = 2*TMath::Pi();
72 Double_t lTemp = par[1];
73 Double_t lPower = par[2];
74 Double_t lMass = par[3];
75 Double_t A = par[4];
76 Double_t lMassEt = par[3];//only used for calculating Et
77 //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
78 Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
79
80 Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
81 Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
82 Double_t strangeness = (0.4333333*pt[0]-0.0666666)/(0.2333333333*pt[0]-0.01666666);
83 if(strangeness<1.0 || pt[0]<0.1) strangeness = 1.0;
84 //This is the density of particles times the Et weight
85 return ldNdy *strangeness* Et* pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
86 }
87 ClassDef(AliAnalysisLevyPtModifiedStrangeness, 1);
88};
89class AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced{
90public:
91 virtual ~AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced(){;};
92
93
94 Double_t Evaluate(Double_t *pt, Double_t *par)
95 {
96 Double_t ldNdy = par[0];
97 Double_t l2pi = 2*TMath::Pi();
98 Double_t lTemp = par[1];
99 Double_t lPower = par[2];
100 Double_t lMass = par[3];
101 Double_t A = par[4];
102 Double_t lMassEt = par[3];//only used for calculating Et
103 Double_t lBary0 = par[6];
104 Double_t lBary1 = par[7];
105 Double_t lBary2 = par[8];
106 Double_t lBary3 = par[9];
107 Double_t lBary4 = par[10];
108 Double_t lBary5 = par[11];
109 //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
110 Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
111
112 Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
113 Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
114 //this is the baryon enhancement factor
115 Double_t baryonEnhancement = lBary0*pow(pt[0],lBary1)*exp(-pow(pt[0]/lBary2,lBary3))/(lBary4+lBary5*pt[0]);
116 if(baryonEnhancement<1.0) baryonEnhancement = 1.0;
117 Double_t strangeness = (0.4333333*pt[0]-0.0666666)/(0.2333333333*pt[0]-0.01666666);
118 if(strangeness<1.0 || pt[0]<0.1) strangeness = 1.0;
119 //This is the density of particles times the Et weight
120 return ldNdy *Et*strangeness* baryonEnhancement*pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower); }
121 ClassDef(AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced, 1);
122};
123
bc92a3aa 124
836ce14a 125void CorrNeutralLevyFit(bool hadronic = false){
126
127 float factor = 0.0;
128 float factorerr = 0.0;
129 if(hadronic){
130 factor = 0.548;
131 factorerr = 0.003;
132 }
bc92a3aa 133
134// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
135// $\pi^{+}+\pi^{-}$ & 2.977 $\pm$ 0.15 & 0.126 $\pm$ 0.001 & 7.82 $\pm$ 0.1 & & & \\
136 AliAnalysisLevyPtModified *function = new AliAnalysisLevyPtModified();
137 fPion = new TF1("fPion",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
138 fPion->SetParameter(0,2.977);//dN/dy
139 fPion->SetParameter(1,0.126);//T
140 fPion->SetParameter(2,7.82);//n
141 fPion->SetParameter(3,0.13957);//mass
142 fPion->SetParameter(4,0);//A=0 for pions
143 float integralPion = fPion->Integral(0,50);
144 float integralErrPion = integralPion*0.15/2.977;
145 float myerrorPionT = 0.0;
146 float myerrorPionn = 0.0;
147 float tmpPion;
148 float T = 0.126;
149 float Terr = 0.001;
150 float n = 7.82;
151 float nerr = 0.1;
152 fPion->SetParameter(1,T+Terr);//T
153 tmpPion = fPion->Integral(0,50);
154 myerrorPionT = TMath::Abs(integralPion-tmpPion);
155 fPion->SetParameter(1,T-Terr);//T
156 tmpPion = fPion->Integral(0,50);
157 if(TMath::Abs(integralPion-tmpPion)>myerrorPionT) myerrorPionT = TMath::Abs(integralPion-tmpPion);
158 fPion->SetParameter(1,T);//T
159 fPion->SetParameter(2,n+nerr);//n
160 tmpPion = fPion->Integral(0,50);
bc92a3aa 161 myerrorPionn = TMath::Abs(integralPion-tmpPion);
162 fPion->SetParameter(2,n-nerr);//n
163 tmpPion = fPion->Integral(0,50);
bc92a3aa 164 if(TMath::Abs(integralPion-tmpPion)>myerrorPionn) myerrorPionn = TMath::Abs(integralPion-tmpPion);
bc92a3aa 165 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
836ce14a 166 //To get the correct error one would have to fit the spectra data to get the covariance matrix...
bc92a3aa 167 integralErrPion = TMath::Sqrt(TMath::Power(integralErrPion,2)+TMath::Power(myerrorPionT,2)+TMath::Power(myerrorPionn,2));
168 cout<<"Pion Et = "<<integralPion<<"$\\pm$"<<integralErrPion<<endl;
169
170// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
171// $K^{+}+K^{-}$ & 0.366 $\pm$ 0.03 & 0.160 $\pm$ 0.006 & 6.087 $\pm$ 0.4 & & & \\
172 fKaon = new TF1("fKaon",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
173 fKaon->SetParameter(0,0.366);//dN/dy
174 fKaon->SetParameter(1,0.160);//T
175 fKaon->SetParameter(2,6.087);//n
176 fKaon->SetParameter(3,0.493677);//mass
177 fKaon->SetParameter(4,0);//A=0 for kaons
178 float integralKaon = fKaon->Integral(0,50);
179 float integralErrKaon = integralKaon*0.03/0.366;
df8b9b1f 180 fKaonStrange = new TF1("fKaonStrange",function, &AliAnalysisLevyPtModifiedStrangeness::Evaluate,0,50,5,"AliAnalysisLevyPtModifiedStrangeness","Evaluate");
181 for(int i=0; i<5;i++){fKaonStrange->SetParameter(i,fKaon->GetParameter(i));}
182 float integralKaonStrange = fKaonStrange->Integral(0,50);
bc92a3aa 183 float myerrorKaonT = 0.0;
184 float myerrorKaonn = 0.0;
185 float tmpKaon;
186 float T = 0.160;
187 float Terr = 0.006;
188 float n = 6.087;
189 float nerr = 0.4;
190 fKaon->SetParameter(1,T+Terr);//T
191 tmpKaon = fKaon->Integral(0,50);
192 myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
193 fKaon->SetParameter(1,T-Terr);//T
194 tmpKaon = fKaon->Integral(0,50);
195 if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonT) myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
196 fKaon->SetParameter(1,T);//T
197 fKaon->SetParameter(2,n+nerr);//n
198 tmpKaon = fKaon->Integral(0,50);
bc92a3aa 199 myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
200 fKaon->SetParameter(2,n-nerr);//n
201 tmpKaon = fKaon->Integral(0,50);
bc92a3aa 202 if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonn) myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
bc92a3aa 203 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
204 integralErrKaon = TMath::Sqrt(TMath::Power(integralErrKaon,2)+TMath::Power(myerrorKaonT,2)+TMath::Power(myerrorKaonn,2));
205 cout<<"Kaon Et = "<<integralKaon<<"$\\pm$"<<integralErrKaon<<endl;
df8b9b1f 206 cout<<"Kaon Et Strange = "<<integralKaonStrange<<endl;
bc92a3aa 207
208// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
209// $p + \bar{p}$ & 0.157 $\pm$ 0.012 & 0.196 $\pm$ 0.009 & 8.6 $\pm$ 1.1 & & & \\
210 fProton = new TF1("fProton",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
211 fProton->SetParameter(0,0.157/2.0);//dN/dy
212 fProton->SetParameter(1,0.196);//T
213 fProton->SetParameter(2,8.6);//n
214 fProton->SetParameter(3,0.938272);//mass
215 fProton->SetParameter(4,-1);//A=0 for protons
df8b9b1f 216 fProtonEnhanced = new TF1("fProtonEnhanced",function, &AliAnalysisLevyPtModifiedBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedBaryonEnhanced","Evaluate");
217 for(int i=0;i<6;i++){fProtonEnhanced->SetParameter(i,fProton->GetParameter(i));}//set all of the spectra parameters to their normal values
218 //and now set the baryon enhancement parameters
219 fProtonEnhanced->SetParameter(6,0.900878*1.2);
220 fProtonEnhanced->SetParameter(7,1.38882);
221 fProtonEnhanced->SetParameter(8,2.6361);
222 fProtonEnhanced->SetParameter(9,1.37751);
223 fProtonEnhanced->SetParameter(10,0.5);
224 fProtonEnhanced->SetParameter(11,-0.03);
225 float integralProtonEnhanced = fProtonEnhanced->Integral(0,50);
bc92a3aa 226 float integralProton = fProton->Integral(0,50);
227 float integralErrProton = integralProton*0.012/0.157;
228 float myerrorProtonT = 0.0;
229 float myerrorProtonn = 0.0;
230 float tmpProton;
231 float T = 0.196;
232 float Terr = 0.009;
233 float n = 8.6;
234 float nerr = 1.1;
235 fProton->SetParameter(1,T+Terr);//T
236 tmpProton = fProton->Integral(0,50);
237 myerrorProtonT = TMath::Abs(integralProton-tmpProton);
238 fProton->SetParameter(1,T-Terr);//T
239 tmpProton = fProton->Integral(0,50);
240 if(TMath::Abs(integralProton-tmpProton)>myerrorProtonT) myerrorProtonT = TMath::Abs(integralProton-tmpProton);
241 fProton->SetParameter(1,T);//T
242 fProton->SetParameter(2,n+nerr);//n
243 tmpProton = fProton->Integral(0,50);
bc92a3aa 244 myerrorProtonn = TMath::Abs(integralProton-tmpProton);
245 fProton->SetParameter(2,n-nerr);//n
246 tmpProton = fProton->Integral(0,50);
bc92a3aa 247 if(TMath::Abs(integralProton-tmpProton)>myerrorProtonn) myerrorProtonn = TMath::Abs(integralProton-tmpProton);
bc92a3aa 248 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
249 integralErrProton = TMath::Sqrt(TMath::Power(integralErrProton,2)+TMath::Power(myerrorProtonT,2)+TMath::Power(myerrorProtonn,2));
250 cout<<"Proton Et = "<<integralProton<<"$\\pm$"<<integralErrProton<<endl;
251
252
253
df8b9b1f 254
255
bc92a3aa 256 //Antiprotons...
df8b9b1f 257 fProton->SetParameter(0,0.157/2.0);//dN/dy
258 fProton->SetParameter(1,0.196);//T
259 fProton->SetParameter(2,8.6);//n
260 fProton->SetParameter(3,0.938272);//mass
bc92a3aa 261 fProton->SetParameter(2,n);//n
262 fProton->SetParameter(4,1);//A=0 for protons
263 float integralAntiProton = fProton->Integral(0,50);
264 float integralErrAntiProton = integralAntiProton*0.012/0.157;
df8b9b1f 265 fAntiProtonEnhanced = new TF1("fAntiProtonEnhanced",function, &AliAnalysisLevyPtModifiedBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedBaryonEnhanced","Evaluate");
266 for(int i=0;i<6;i++){fAntiProtonEnhanced->SetParameter(i,fProton->GetParameter(i));}//set all of the spectra parameters to their normal values
267 fAntiProtonEnhanced->SetParameter(2,n);//n
268 fAntiProtonEnhanced->SetParameter(4,1);//A=0 for protons
269 //and now set the baryon enhancement parameters
270 fAntiProtonEnhanced->SetParameter(6,0.900878*1.2);
271 fAntiProtonEnhanced->SetParameter(7,1.38882);
272 fAntiProtonEnhanced->SetParameter(8,2.6361);
273 fAntiProtonEnhanced->SetParameter(9,1.37751);
274 fAntiProtonEnhanced->SetParameter(10,0.5);
275 fAntiProtonEnhanced->SetParameter(11,-0.03);
276 float integralAntiProtonEnhanced = fAntiProtonEnhanced->Integral(0,50);
bc92a3aa 277 float myerrorAntiProtonT = 0.0;
278 float myerrorAntiProtonn = 0.0;
279 float tmpAntiProton;
280 float T = 0.196;
281 float Terr = 0.009;
282 float n = 8.6;
283 float nerr = 1.1;
284 fProton->SetParameter(1,T+Terr);//T
285 tmpAntiProton = fProton->Integral(0,50);
286 myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
287 fProton->SetParameter(1,T-Terr);//T
288 tmpAntiProton = fProton->Integral(0,50);
289 if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonT) myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
290 fProton->SetParameter(1,T);//T
291 fProton->SetParameter(2,n+nerr);//n
292 tmpAntiProton = fProton->Integral(0,50);
bc92a3aa 293 myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
294 fProton->SetParameter(2,n-nerr);//n
295 tmpAntiProton = fProton->Integral(0,50);
bc92a3aa 296 if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonn) myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
bc92a3aa 297 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
298 integralErrAntiProton = TMath::Sqrt(TMath::Power(integralErrAntiProton,2)+TMath::Power(myerrorAntiProtonT,2)+TMath::Power(myerrorAntiProtonn,2));
299 cout<<"AntiProton Et = "<<integralAntiProton<<"$\\pm$"<<integralErrAntiProton<<endl;
df8b9b1f 300 cout<<"AntiProton enhanced "<<integralAntiProtonEnhanced<<", "<<integralAntiProton/integralAntiProtonEnhanced*100.0<<endl;
301 cout<<"Proton enhanced "<<integralProtonEnhanced<<", "<<integralProton/integralProtonEnhanced*100.0<<endl;
302 cout<<"Proton and antiproton enhanced "<< (integralProtonEnhanced+integralAntiProtonEnhanced) <<", "
303 << (integralProtonEnhanced+integralAntiProtonEnhanced) / (integralProton+integralAntiProton)*100.0
304 <<endl;
bc92a3aa 305
306
307// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
308// \kzeroshort &0.184 $\pm$ 0.006 & 0.168 $\pm$ 0.005 & 6.6 $\pm$ 0.3 & & & \\
309 fK0S = new TF1("fK0S",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
310 fK0S->SetParameter(0,0.184);//dN/dy
311 fK0S->SetParameter(1,0.168);//T
312 fK0S->SetParameter(2,6.6);//n
313 fK0S->SetParameter(3,.497614);//mass
314 fK0S->SetParameter(4,0);//A=0 for kaons
315 float integralK0S = fK0S->Integral(0,50);
316 float integralErrK0S = integralK0S*0.006/0.184;
df8b9b1f 317 fK0Strange = new TF1("fK0Strange",function, &AliAnalysisLevyPtModifiedStrangeness::Evaluate,0,50,5,"AliAnalysisLevyPtModifiedStrangeness","Evaluate");
318 for(int i=0; i<5;i++){fK0Strange->SetParameter(i,fK0S->GetParameter(i));}
319 float integralK0SStrange = fK0Strange->Integral(0,50);
bc92a3aa 320 float myerrorK0ST = 0.0;
321 float myerrorK0Sn = 0.0;
322 float tmpK0S;
323 float T = 0.184;
324 float Terr = 0.006;
325 float n = 6.6;
326 float nerr = 0.3;
327 fK0S->SetParameter(1,T+Terr);//T
328 tmpK0S = fK0S->Integral(0,50);
329 myerrorK0ST = TMath::Abs(integralK0S-tmpK0S);
330 fK0S->SetParameter(1,T-Terr);//T
331 tmpK0S = fK0S->Integral(0,50);
332 if(TMath::Abs(integralK0S-tmpK0S)>myerrorK0ST) myerrorK0ST = TMath::Abs(integralK0S-tmpK0S);
333 fK0S->SetParameter(1,T);//T
334 fK0S->SetParameter(2,n+nerr);//n
335 tmpK0S = fK0S->Integral(0,50);
bc92a3aa 336 myerrorK0Sn = TMath::Abs(integralK0S-tmpK0S);
337 fK0S->SetParameter(2,n-nerr);//n
338 tmpK0S = fK0S->Integral(0,50);
bc92a3aa 339 if(TMath::Abs(integralK0S-tmpK0S)>myerrorK0Sn) myerrorK0Sn = TMath::Abs(integralK0S-tmpK0S);
bc92a3aa 340 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
341 integralErrK0S = TMath::Sqrt(TMath::Power(integralErrK0S,2)+TMath::Power(myerrorK0ST,2)+TMath::Power(myerrorK0Sn,2));
342 cout<<"K0S Et = "<<integralK0S<<"$\\pm$"<<integralErrK0S<<endl;
df8b9b1f 343 cout<<"K0S Et Strange = "<<integralK0SStrange<<endl;
bc92a3aa 344
345// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
346// \lam &0.048 $\pm$ 0.004 & 0.229 $\pm$ 0.015 & 10.8 $\pm$ 2.0 & & & \\
347 fLambda = new TF1("fLambda",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
348 fLambda->SetParameter(0,0.048);//dN/dy
349 fLambda->SetParameter(1,0.229);//T
350 fLambda->SetParameter(2,10.8);//n
351 fLambda->SetParameter(3,1.115683);//mass
352 fLambda->SetParameter(4,0);//A=0 for kaons
353 float integralLambda = fLambda->Integral(0,50);
354 float integralErrLambda = integralLambda*0.004/0.048;
df8b9b1f 355 fLambdaEnhanced = new TF1("fLambdaEnhanced",function, &AliAnalysisLevyPtModifiedBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedBaryonEnhanced","Evaluate");
356 for(int i=0;i<6;i++){fLambdaEnhanced->SetParameter(i,fLambda->GetParameter(i));}//set all of the spectra parameters to their normal values
357 //and now set the baryon enhancement parameters
358 fLambdaEnhanced->SetParameter(6,0.900878);
359 fLambdaEnhanced->SetParameter(7,1.38882);
360 fLambdaEnhanced->SetParameter(8,2.6361);
361 fLambdaEnhanced->SetParameter(9,1.37751);
362 fLambdaEnhanced->SetParameter(10,0.5);
363 fLambdaEnhanced->SetParameter(11,-0.03);
364 float integralLambdaEnhanced = fLambdaEnhanced->Integral(0,50);
365// fLambdaEnhanced->Draw();
366// return;
367 fLambdaStrange = new TF1("fLambdaStrange",function, &AliAnalysisLevyPtModifiedStrangeness::Evaluate,0,50,5,"AliAnalysisLevyPtModifiedStrangeness","Evaluate");
368 for(int i=0; i<5;i++){fLambdaStrange->SetParameter(i,fLambda->GetParameter(i));}
369 //for(int i=0;i<12;i++){cout<<"Lambda a"<<i<<" "<<fLambdaStrange->GetParameter(i)<<" b"<<i<<" "<<fLambda->GetParameter(i)<<endl;}
370 float integralLambdaStrange = fLambdaStrange->Integral(0,50);
371 //cout<<"test "<<integralLambdaStrange<<endl;
372 fLambdaStrangeEnhanced = new TF1("fLambdaStrangeEnhanced",function, &AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced","Evaluate");
373 for(int i=0; i<12;i++){fLambdaStrangeEnhanced->SetParameter(i,fLambdaEnhanced->GetParameter(i));}
374// //and now set the baryon enhancement parameters
375// fLambdaStrangeEnhanced->SetParameter(6,0.900878);
376// fLambdaStrangeEnhanced->SetParameter(7,1.38882);
377// fLambdaStrangeEnhanced->SetParameter(8,2.6361);
378// fLambdaStrangeEnhanced->SetParameter(9,1.37751);
379// fLambdaStrangeEnhanced->SetParameter(10,0.5);
380// fLambdaStrangeEnhanced->SetParameter(11,-0.03);
381 //for(int i=0;i<12;i++){cout<<"a"<<i<<" "<<fLambdaStrangeEnhanced->GetParameter(i)<<" b"<<i<<" "<<fLambdaEnhanced->GetParameter(i)<<" c"<<i<<" "<<fLambda->GetParameter(i)<<endl;}
382 float integralLambdaStrangeEnhanced = fLambdaStrangeEnhanced->Integral(0,50);
383 cout<<"Lambda enhanced "<<integralLambdaEnhanced<<", ";
384 if(integralLambdaEnhanced>0.0) cout<<integralLambda/integralLambdaEnhanced*100.0;
385 cout<<endl;
bc92a3aa 386 float myerrorLambdaT = 0.0;
387 float myerrorLambdan = 0.0;
388 float tmpLambda;
389 float T = 0.229;
390 float Terr = 0.015;
391 float n = 10.8;
392 float nerr = 2.0;
393 fLambda->SetParameter(1,T+Terr);//T
394 tmpLambda = fLambda->Integral(0,50);
395 myerrorLambdaT = TMath::Abs(integralLambda-tmpLambda);
396 fLambda->SetParameter(1,T-Terr);//T
397 tmpLambda = fLambda->Integral(0,50);
398 if(TMath::Abs(integralLambda-tmpLambda)>myerrorLambdaT) myerrorLambdaT = TMath::Abs(integralLambda-tmpLambda);
399 fLambda->SetParameter(1,T);//T
400 fLambda->SetParameter(2,n+nerr);//n
401 tmpLambda = fLambda->Integral(0,50);
bc92a3aa 402 myerrorLambdan = TMath::Abs(integralLambda-tmpLambda);
403 fLambda->SetParameter(2,n-nerr);//n
404 tmpLambda = fLambda->Integral(0,50);
bc92a3aa 405 if(TMath::Abs(integralLambda-tmpLambda)>myerrorLambdan) myerrorLambdan = TMath::Abs(integralLambda-tmpLambda);
bc92a3aa 406 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
407 integralErrLambda = TMath::Sqrt(TMath::Power(integralErrLambda,2)+TMath::Power(myerrorLambdaT,2)+TMath::Power(myerrorLambdan,2));
408 cout<<"Lambda Et = "<<integralLambda<<"$\\pm$"<<integralErrLambda<<endl;
df8b9b1f 409 cout<<"Lambda Et Strange = "<<integralLambdaStrange<<endl;
410 cout<<"Lambda Et Strange Enhanced = "<<integralLambdaStrangeEnhanced<<endl;
bc92a3aa 411
412
413// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
414// \alam & 0.047 $\pm$ 0.005 & 0.210 $\pm$ 0.015 & 9.2 $\pm$ 1.4 & & & \\ & & \\
415 fAntiLambda = new TF1("fAntiLambda",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
416 fAntiLambda->SetParameter(0,0.047);//dN/dy
417 fAntiLambda->SetParameter(1,0.210);//T
418 fAntiLambda->SetParameter(2,9.2);//n
419 fAntiLambda->SetParameter(3,1.115683);//mass
420 fAntiLambda->SetParameter(4,0);//A=0 for kaons
421 float integralAntiLambda = fAntiLambda->Integral(0,50);
422 float integralErrAntiLambda = integralAntiLambda*0.005/0.047;
df8b9b1f 423 fAntiLambdaEnhanced = new TF1("fAntiLambdaEnhanced",function, &AliAnalysisLevyPtModifiedBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedBaryonEnhanced","Evaluate");
424 for(int i=0;i<6;i++){fAntiLambdaEnhanced->SetParameter(i,fAntiLambda->GetParameter(i));}//set all of the spectra parameters to their normal values
425 //and now set the baryon enhancement parameters
426 fAntiLambdaEnhanced->SetParameter(6,0.900878*1.2);
427 fAntiLambdaEnhanced->SetParameter(7,1.38882);
428 fAntiLambdaEnhanced->SetParameter(8,2.6361);
429 fAntiLambdaEnhanced->SetParameter(9,1.37751);
430 fAntiLambdaEnhanced->SetParameter(10,0.5);
431 fAntiLambdaEnhanced->SetParameter(11,-0.03);
432 float integralAntiLambdaEnhanced = fAntiLambdaEnhanced->Integral(0,50);
433 cout<<"AntiLambda enhanced "<<integralAntiLambdaEnhanced<<", "<<integralAntiLambda/integralAntiLambdaEnhanced*100.0<<endl;
434 fAntiLambdaStrange = new TF1("fAntiLambdaStrange",function, &AliAnalysisLevyPtModifiedStrangeness::Evaluate,0,50,5,"AliAnalysisLevyPtModifiedStrangeness","Evaluate");
435 for(int i=0; i<5;i++){fAntiLambdaStrange->SetParameter(i,fAntiLambda->GetParameter(i));}
436 float integralAntiLambdaStrange = fAntiLambdaStrange->Integral(0,50);
437
438 fAntiLambdaStrangeEnhanced = new TF1("fAntiLambdaStrangeEnhanced",function, &AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced","Evaluate");
439 for(int i=0; i<12;i++){fAntiLambdaStrangeEnhanced->SetParameter(i,fAntiLambdaEnhanced->GetParameter(i));}
440 float integralAntiLambdaStrangeEnhanced = fAntiLambdaStrangeEnhanced->Integral(0,50);
441
bc92a3aa 442 float myerrorAntiLambdaT = 0.0;
443 float myerrorAntiLambdan = 0.0;
444 float tmpAntiLambda;
445 float T = 0.210;
446 float Terr = 0.015;
447 float n = 9.2;
448 float nerr = 1.4;
449 fAntiLambda->SetParameter(1,T+Terr);//T
450 tmpAntiLambda = fAntiLambda->Integral(0,50);
451 myerrorAntiLambdaT = TMath::Abs(integralAntiLambda-tmpAntiLambda);
452 fAntiLambda->SetParameter(1,T-Terr);//T
453 tmpAntiLambda = fAntiLambda->Integral(0,50);
454 if(TMath::Abs(integralAntiLambda-tmpAntiLambda)>myerrorAntiLambdaT) myerrorAntiLambdaT = TMath::Abs(integralAntiLambda-tmpAntiLambda);
455 fAntiLambda->SetParameter(1,T);//T
456 fAntiLambda->SetParameter(2,n+nerr);//n
457 tmpAntiLambda = fAntiLambda->Integral(0,50);
bc92a3aa 458 myerrorAntiLambdan = TMath::Abs(integralAntiLambda-tmpAntiLambda);
459 fAntiLambda->SetParameter(2,n-nerr);//n
460 tmpAntiLambda = fAntiLambda->Integral(0,50);
bc92a3aa 461 if(TMath::Abs(integralAntiLambda-tmpAntiLambda)>myerrorAntiLambdan) myerrorAntiLambdan = TMath::Abs(integralAntiLambda-tmpAntiLambda);
bc92a3aa 462 //This isn't strictly correct because the errors on the parameters should be correlated but it's close
463 integralErrAntiLambda = TMath::Sqrt(TMath::Power(integralErrAntiLambda,2)+TMath::Power(myerrorAntiLambdaT,2)+TMath::Power(myerrorAntiLambdan,2));
464 cout<<"AntiLambda Et = "<<integralAntiLambda<<"$\\pm$"<<integralErrAntiLambda<<endl;
df8b9b1f 465 cout<<"AntiLambda Et Strange = "<<integralAntiLambdaStrange<<endl;
466 cout<<"AntiLambda Et Strange Enhanced = "<<integralAntiLambdaStrangeEnhanced<<endl;
bc92a3aa 467
468 //now we apply various assumptions
469 float integralK0L = integralK0S;
470 float integralErrK0L = integralErrK0S;
471 float integralNeutron = integralProton;
472 float integralErrNeutron = integralErrProton;
473
836ce14a 474 float totalEt = (1.0+factor)*integralPion+integralKaon+2.0*integralProton+2.0*integralAntiProton+2.0*integralK0S+integralLambda+integralAntiLambda;
bc92a3aa 475 float measuredEt = integralPion+integralKaon+1.0*integralProton+1.0*integralAntiProton;
476 float fneutral = measuredEt/totalEt;
477 cout<<"fneutral = "<<fneutral<<endl;
df8b9b1f 478
479 float totalEtEnhanced = (1.0+factor)*integralPion+integralKaon+2.0*integralProtonEnhanced+2.0*integralAntiProtonEnhanced+2.0*integralK0S+integralLambdaEnhanced+integralAntiLambdaEnhanced;
480 float measuredEtEnhanced = integralPion+integralKaon+1.0*integralProtonEnhanced+1.0*integralAntiProtonEnhanced;
481 float fneutralEnhanced = measuredEtEnhanced/totalEtEnhanced;
482 cout<<"fneutralEnhanced = "<<fneutralEnhanced<<endl;
483
484
485 float totalEtStrange = (1.0+factor)*integralPion+integralKaonStrange+2.0*integralProton+2.0*integralAntiProton+2.0*integralK0SStrange+integralLambdaStrange+integralAntiLambdaStrange;
486 float measuredEtStrange = integralPion+integralKaonStrange+1.0*integralProton+1.0*integralAntiProton;
487 float fneutralStrange = measuredEtStrange/totalEtStrange;
488 cout<<"fneutralStrange = "<<fneutralStrange<<endl;
489
490
491
492 float totalEtStrangeEnhanced = (1.0+factor)*integralPion+integralKaonStrange+2.0*integralProtonEnhanced+2.0*integralAntiProtonEnhanced+2.0*integralK0SStrange+integralLambdaStrangeEnhanced+integralAntiLambdaStrangeEnhanced;
493 float measuredEtStrangeEnhanced = integralPion+integralKaonStrange+1.0*integralProtonEnhanced+1.0*integralAntiProtonEnhanced;
494 float fneutralStrangeEnhanced = measuredEtStrangeEnhanced/totalEtStrangeEnhanced;
495 cout<<"fneutralStrangeEnhanced = "<<fneutralStrangeEnhanced<<endl;
496
bc92a3aa 497 //this is from ugly derivative taking for error propagation
498 //form 1: pions, kaons, protons
499 //for f=(xa+A)/(ya+B)
500 //df/da=(xf-yf^2)/(xa+A)
501 //x=y=1; xa+A=measuredEt
502 float errPion = (fneutral-fneutral*fneutral)/measuredEt*integralErrPion;
836ce14a 503 if(hadronic){
504 //then we have x=1, y = 1+factor
505 //the error on the extra bit of et from pi0s
506 float extraerror = integralPion*factorerr;
507 integralErrPion = TMath::Sqrt(TMath::Power(integralErrPion,2)+TMath::Power(extraerror,2));
508 errPion = (fneutral-(1.0+factor)*fneutral*fneutral)/measuredEt*integralErrPion;
509 }
bc92a3aa 510 float errKaon = (fneutral-fneutral*fneutral)/measuredEt*integralErrKaon;
511 //x=1,y=2
512 float errProton = (fneutral-2.0*fneutral*fneutral)/measuredEt*integralErrProton;
513 float errAntiProton = (fneutral-2.0*fneutral*fneutral)/measuredEt*integralErrAntiProton;
514 //form 2: K0S, lambda, antilambda
515 //for f=A/(B+xa)
516 //df/da=-xf^2/A
517 //x=1,A=measuredEt
518 float errLambda = fneutral*fneutral/measuredEt*integralErrLambda;
519 float errAntiLambda = fneutral*fneutral/measuredEt*integralErrAntiLambda;
520 //x=2,A=measuredEt
521 float errK0S = 2.0*fneutral*fneutral/measuredEt*integralErrK0S;
522 cout<<"Errors Pion "<<errPion<<" Kaon "<<errKaon<<" Proton "<<errProton<<" Lambda "<<errLambda<<" AntiLambda "<<errAntiLambda<<" K0S "<<errK0S<<endl;
523 float totalErr = TMath::Sqrt( TMath::Power(errPion,2) + TMath::Power(errKaon,2) + TMath::Power(errProton,2)+ TMath::Power(errAntiProton,2) + TMath::Power(errLambda,2) + TMath::Power(errAntiLambda,2) + TMath::Power(errK0S,2) );
524
525 cout<<"fneutral = "<<fneutral<<"$\\pm$"<<totalErr<<endl;
526 cout<<"1/fneutral = "<<1.0/fneutral<<"$\\pm$"<<1.0/fneutral*totalErr/fneutral<<endl;
527 cout<<"Percentage error "<<totalErr/fneutral*100.0<<endl;
df8b9b1f 528
529 float fneutralPb = (fneutralStrangeEnhanced+fneutral)/2.0;
530 float fneutralPbErr = TMath::Sqrt( TMath::Power((fneutralStrangeEnhanced-fneutral)/2.0,2.0) + TMath::Power(totalErr,2.0));
531 cout<<"fneutral strange enhanced "<<fneutralPb<<"$\\pm$"<<fneutralPbErr<<endl;
532 cout<<"1/fneutral strange enhanced = "<<1.0/fneutralPb<<"$\\pm$"<<1.0/fneutralPb*fneutralPbErr/fneutralPb<<endl;
533 cout<<"Percentage error "<<fneutralPbErr/fneutralPb*100.0<<endl;
534
bc92a3aa 535}
536