]>
Commit | Line | Data |
---|---|---|
e0e1d4c6 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //------------------------------------------------------------------------- | |
17 | // author: Sergey Kiselev, ITEP, Moscow | |
18 | // e-mail: Sergey.Kiselev@cern.ch | |
19 | // tel.: 007 495 129 95 45 | |
20 | //------------------------------------------------------------------------- | |
21 | // Generator of direct thermal photons for the reaction A+B, sqrt(S) | |
22 | // main assumptions: | |
23 | // 1+1 Bjorken scaling hydrodinamics. | |
24 | // 1st order phase transition | |
25 | // QGP + Mixed (QGP+HHG) + HHG (Hot Hadron Gas) phases, | |
26 | // an ideal massless parton gas and ideal massless HHG | |
27 | // see | |
28 | // F.D.Steffen, nucl-th/9909035 | |
29 | // F.D.Steffen and M.H.Thoma, Phys.Lett. B510, 98 (2001) | |
30 | // T.Peitzmann and M.H.Thoma, Phys.Rep., 364, 175 (2002) | |
31 | // | |
32 | // photon rates for QGP: Phys.Rep., 364, 175 (2002), section 2.1.1 | |
33 | // | |
34 | // photon rates for HHG | |
35 | // prates for i rho --> pi gamma, pi pi --> rho gamma and rho --> pi pi gamma: | |
36 | // Song and Fai, Phys.Rev. C58, 1689 (1998) | |
37 | // rates for omega --> pi gamma: Phys.Rev. D44, 2774 (1991) | |
38 | // | |
39 | // input parameters: | |
40 | // fAProjectile, fATarget - number of nucleons in a nucleus A and B | |
41 | // fMinImpactParam - minimal impct parameter, fm | |
42 | // fMaxImpactParam - maximal impct parameter, fm | |
43 | // fEnergyCMS - sqrt(S) per nucleon pair, AGeV | |
44 | // fTau0 - initial proper time, fm/c | |
45 | // fT0 - initial temperature, GeV | |
46 | // fTc - critical temperature, GeV | |
47 | // fTf - freeze-out temperature, GeV | |
48 | // fGhhg - effective number of degrees of freedom in HHG | |
49 | // | |
50 | // fYMin - minimal rapidity of photons | |
51 | // fYMax - maximal rapidity of photons | |
52 | // in [fYMin,fYMax] uniform distribution of gamma is assumed | |
53 | // fPtMin - minimal p_t value of gamma, GeV/c | |
54 | // fPtMax - maximal p_t value of gamma, GeV/c | |
55 | //------------------------------------------------------------------------- | |
56 | // comparison with SPS and RHIC data, prediction for LHC can be found in | |
57 | // arXiv:0811.2634 [nucl-th] | |
58 | //------------------------------------------------------------------------- | |
59 | ||
60 | //Begin_Html | |
61 | /* | |
62 | <img src="picts/AliGeneratorClass.gif"> | |
63 | </pre> | |
64 | <br clear=left> | |
65 | <font size=+2 color=red> | |
66 | <p>The responsible person for this module is | |
67 | <a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>. | |
68 | </font> | |
69 | <pre> | |
70 | */ | |
71 | //End_Html | |
72 | // // | |
73 | /////////////////////////////////////////////////////////////////// | |
74 | ||
75 | #include <TArrayF.h> | |
76 | #include <TF1.h> | |
77 | #include <TF2.h> | |
78 | #include <TH1F.h> | |
79 | ||
80 | #include "AliConst.h" | |
81 | #include "AliGenEventHeader.h" | |
82 | #include "AliGenThermalPhotons.h" | |
83 | #include "AliLog.h" | |
84 | #include "AliRun.h" | |
85 | ||
86 | ClassImp(AliGenThermalPhotons) | |
87 | ||
88 | // ----------------------------------------------------------------------------------------------------- | |
89 | static Double_t rateQGP(Double_t *x, Double_t *par) { | |
90 | //--------------------------------------------------- | |
91 | // input: | |
92 | // x[0] - tau (fm), proper time | |
93 | // x[1] - yprime, space rapidity | |
94 | // par[0] - p_T (GeV), photon transverse momentum | |
95 | // par[1] - y, photon rapidity in the c.m.s. A+A | |
96 | // par[2] - tau0 (fm), initial proper time | |
97 | // par[3] - T_0 (GeV), initial temperature | |
98 | // par[4] - T_c (GeV), critical temperature | |
99 | // par[5] - iProcQGP, process number, 0, 1, 2 | |
100 | // | |
101 | // output: | |
102 | // tau EdR/d^3p = tau EdN/d^4xd^3p (fm fm^{-4}GeV^{-2}) | |
103 | //--------------------------------------------------- | |
104 | Double_t tau=x[0],yprime=x[1]; | |
6274ddfc | 105 | Double_t pT=par[0],y=par[1],tau0=par[2],t0=par[3],tc=par[4]; |
e0e1d4c6 | 106 | Int_t iProcQGP=Int_t(par[5]), nFl=3; |
107 | ||
6274ddfc | 108 | Double_t e=pT*TMath::CosH(yprime-y),t=t0*TMath::Power(tau0/tau,1./3.); |
e0e1d4c6 | 109 | |
110 | const Double_t alpha=1./137.; | |
111 | // factor to convert from GeV^2 to fm^{-4}GeV^{-2}: (1/hc)**4=(1/0.197)**4=659.921 | |
112 | const Double_t factor=659.921; | |
6274ddfc | 113 | Double_t alphaS=3.*TMath::TwoPi()/((33.-2.*nFl)*TMath::Log(8.*t/tc)); |
e0e1d4c6 | 114 | const Double_t abc[3]={0.0338, 0.0281, 0.0135} ; // a, b, c for nFf=3 |
115 | Double_t rate=1.; | |
116 | ||
117 | switch (iProcQGP) { | |
118 | ||
119 | case 0: | |
120 | // 1-loop | |
6274ddfc | 121 | rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t); |
e0e1d4c6 | 122 | break ; |
123 | ||
124 | case 1: | |
125 | // 2-loop: bremsstrahlung | |
6274ddfc | 126 | rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t; |
e0e1d4c6 | 127 | break ; |
128 | ||
129 | case 2: | |
130 | // 2-loop: annihilation with scattering | |
6274ddfc | 131 | rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e; |
e0e1d4c6 | 132 | break ; |
133 | ||
134 | default: | |
135 | printf("NO iProcQGP=%i \n",iProcQGP); | |
136 | } | |
137 | ||
138 | return tau*rate; | |
139 | } | |
140 | ||
141 | // ----------------------------------------------------------------------------------------------------- | |
142 | static Double_t fromQGP(Double_t *x, Double_t *par) { | |
143 | //--------------------------------------------------- | |
144 | // input: | |
145 | // x[0] - p_T (GeV), photon p_T | |
146 | // par[0] - tau0 (fm), initial proper time | |
147 | // par[1] - T_0 (GeV), initial temperature | |
148 | // par[2] - tauCQGP (fm), end of QGP | |
149 | // par[3] - yNucl, rapidity of projectile nucleus | |
150 | // par[4] - T_c (GeV), critical temperature | |
151 | // par[5] - y, photon rapidity | |
152 | // par[6] - iProcQGP, process number | |
153 | // | |
154 | // output: | |
155 | // d^{2}N/(dp_t dy) (1/GeV) | |
156 | //--------------------------------------------------- | |
157 | Double_t pT=x[0]; | |
6274ddfc | 158 | Double_t tau0=par[0],t0=par[1],tauCQGP=par[2],yNucl=par[3],tc=par[4],y=par[5]; |
e0e1d4c6 | 159 | Int_t iProcQGP=Int_t(par[6]); |
160 | ||
161 | TF2 frateQGP("frateQGP",&rateQGP,tau0,tauCQGP,-yNucl,yNucl,6); | |
6274ddfc | 162 | frateQGP.SetParameters(pT,y,tau0,t0,tc,iProcQGP); |
e0e1d4c6 | 163 | frateQGP.SetParNames("transverse momentum","rapidity","initial time","initial temperature","critical temperature","process number"); |
164 | return TMath::TwoPi()*pT*frateQGP.Integral(tau0,tauCQGP,-yNucl,yNucl,1e-6); | |
165 | } | |
166 | ||
167 | // ----------------------------------------------------------------------------------------------------- | |
168 | static Double_t rateMixQ(Double_t *x, Double_t *par) { | |
169 | //--------------------------------------------------- | |
170 | // input: | |
171 | // x[0] - yprime, space rapidity | |
172 | // par[0] - p_T (GeV), photon transverse momentum | |
173 | // par[1] - y, photon rapidity in the c.m.s. A+A | |
174 | // par[2] - T_c (GeV), critical temperature | |
175 | // par[3] - iProcQGP, process number, 0, 1, 2 | |
176 | // | |
177 | // output: | |
178 | // EdR/d^3p = EdN/d^4xd^3p (fm fm^{-4}GeV^{-2}) | |
179 | //--------------------------------------------------- | |
180 | Double_t yprime=x[0]; | |
6274ddfc | 181 | Double_t pT=par[0],y=par[1],tc=par[2]; |
e0e1d4c6 | 182 | Int_t iProcQGP=Int_t(par[3]),nFl=3; |
183 | ||
6274ddfc | 184 | Double_t e=pT*TMath::CosH(yprime-y),t=tc; |
e0e1d4c6 | 185 | |
186 | const Double_t alpha=1./137.; | |
187 | // factor to convert from GeV^2 to fm^{-4}GeV^{-2}: (1/hc)**4=(1/0.197)**4=659.921 | |
188 | const Double_t factor=659.921; | |
6274ddfc | 189 | Double_t alphaS=3.*TMath::TwoPi()/((33.-2.*nFl)*TMath::Log(8.*t/tc)); |
e0e1d4c6 | 190 | const Double_t abc[3]={0.0338, 0.0281, 0.0135}; // a, b, c for nF=3 |
191 | Double_t rate=1.; | |
192 | ||
193 | switch (iProcQGP) { | |
194 | ||
195 | case 0: | |
196 | // 1-loop | |
6274ddfc | 197 | rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t); |
e0e1d4c6 | 198 | break ; |
199 | ||
200 | case 1: | |
201 | // 2-loop: bremsstrahlung | |
6274ddfc | 202 | rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t; |
e0e1d4c6 | 203 | break ; |
204 | ||
205 | case 2: | |
206 | // 2-loop: annihilation with scattering | |
6274ddfc | 207 | rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e; |
e0e1d4c6 | 208 | break ; |
209 | ||
210 | default: | |
211 | printf("NO iProcQGP=%i \n",iProcQGP); | |
212 | } | |
213 | ||
214 | return rate; | |
215 | } | |
216 | ||
217 | // ----------------------------------------------------------------------------------------------------- | |
218 | static Double_t fromMixQ(Double_t *x, Double_t *par) { | |
219 | //--------------------------------------------------- | |
220 | // input: | |
221 | // x[0] - p_T (GeV), photon p_T | |
222 | // par[0] - lamQGP | |
223 | // par[1] - yNucl, rapidity of projectile nucleus | |
224 | // par[2] - T_c (GeV), critical temperature | |
225 | // par[3] - y, photon rapidity | |
226 | // par[4] - iProcQGP, process number | |
227 | // | |
228 | // output: | |
229 | // d^{2}N/(dp_t dy) (1/GeV) | |
230 | //--------------------------------------------------- | |
231 | Double_t pT=x[0]; | |
6274ddfc | 232 | Double_t lamQGP=par[0],yNucl=par[1],tc=par[2],y=par[3]; |
e0e1d4c6 | 233 | Int_t iProcQGP=Int_t(par[4]); |
234 | ||
235 | TF1 frateMixQ("frateMixQ",&rateMixQ,-yNucl,yNucl,4); | |
6274ddfc | 236 | frateMixQ.SetParameters(pT,y,tc,iProcQGP); |
e0e1d4c6 | 237 | frateMixQ.SetParNames("transverse momentum","rapidity","critical temperature","process number"); |
238 | return TMath::TwoPi()*pT*lamQGP*frateMixQ.Integral(-yNucl,yNucl); | |
239 | } | |
240 | ||
241 | // ----------------------------------------------------------------------------------------------------- | |
242 | static Double_t rateMixH(Double_t *x, Double_t *par) { | |
243 | //--------------------------------------------------- | |
244 | // input: | |
245 | // x[0] - yprime, space rapidity | |
246 | // par[0] - p_T (GeV), photon transverse momentum | |
247 | // par[1] - y, photon rapidity in the c.m.s. A+A | |
248 | // par[2] - T_c (GeV), critical temperature | |
249 | // par[3] - iProcHHG, process number | |
250 | // | |
251 | // output: | |
252 | // EdR/d^3p = EdN/d^4xd^3p (fm^{-4}GeV^{-2}) | |
253 | //--------------------------------------------------- | |
254 | Double_t yprime=x[0]; | |
6274ddfc | 255 | Double_t pT=par[0],y=par[1],tc=par[2]; |
e0e1d4c6 | 256 | Int_t iProcHHG=Int_t(par[3]); |
257 | ||
6274ddfc | 258 | Double_t e=pT*TMath::CosH(yprime-y),t=tc; |
e0e1d4c6 | 259 | const Double_t mPi=0.135; |
6274ddfc | 260 | Double_t xx=t/mPi,yy=e/mPi; |
261 | Double_t f,rate=1.,emin,factor; | |
262 | const Double_t mOm=0.783,width=0.00844,br=0.085,e0=0.379,pi=TMath::Pi(); | |
e0e1d4c6 | 263 | const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm |
264 | ||
265 | switch (iProcHHG) { | |
266 | ||
267 | case 0: | |
268 | // pi rho --> pi gamma | |
6274ddfc | 269 | f=-2.447+0.796*xx+(0.0338+0.0528*xx)*yy+(-21.447+8.2179*xx)/yy+(1.52436-0.38562*xx)/(yy*yy); |
270 | rate=t*t*TMath::Exp(-e/t+f); | |
e0e1d4c6 | 271 | break ; |
272 | ||
273 | case 1: | |
274 | // pi pi --> rho gamma | |
6274ddfc | 275 | f=-12.055+4.387*xx+(0.3755+0.00826*xx)*yy+(-0.00777+0.000279*xx)*yy*yy+(5.7869-1.0258*xx)/yy+(-1.979+0.58*xx)/(yy*yy); |
276 | rate=t*t*TMath::Exp(-e/t+f); | |
e0e1d4c6 | 277 | break ; |
278 | ||
279 | case 2: | |
280 | // rho --> pi pi gamma | |
6274ddfc | 281 | f=-6.295+1.6459*xx+(-0.4015+0.089*xx)*yy+(-0.954+2.05777*xx)/yy; |
282 | rate=t*t*TMath::Exp(-e/t+f); | |
e0e1d4c6 | 283 | break ; |
284 | ||
285 | case 3: | |
286 | // omega --> pi gamma | |
6274ddfc | 287 | emin=mOm*(e*e+e0*e0)/(2.*e*e0); |
288 | factor=(3.*mOm*width*br)/(16.*pi*pi*pi*e0); | |
289 | rate=factor*t*(emin+t)*TMath::Exp(-emin/t)/e/hc4; | |
e0e1d4c6 | 290 | break ; |
291 | ||
292 | default: | |
293 | printf("NO iProcHHG=%i \n",iProcHHG); | |
294 | } | |
295 | ||
296 | return rate; | |
297 | } | |
298 | ||
299 | // ----------------------------------------------------------------------------------------------------- | |
300 | static Double_t fromMixH(Double_t *x, Double_t *par) { | |
301 | //--------------------------------------------------- | |
302 | // input: | |
303 | // x[0] - p_T (GeV), photon p_T | |
304 | // par[0] - lamHHG | |
305 | // par[1] - yNucl, rapidity of projectile nucleus | |
306 | // par[2] - T_c (GeV), critical temperature | |
307 | // par[3] - y, photon rapidity | |
308 | // par[4] - iProcHHG, process number | |
309 | // | |
310 | // output: | |
311 | // d^{2}N/(dp_t dy) (1/GeV) | |
312 | //--------------------------------------------------- | |
313 | Double_t pT=x[0]; | |
6274ddfc | 314 | Double_t lamHHG=par[0],yNucl=par[1],tc=par[2],y=par[3]; |
e0e1d4c6 | 315 | Int_t iProcHHG=Int_t(par[4]); |
316 | ||
317 | TF1 frateMixH("frateMixH",&rateMixH,-yNucl,yNucl,4); | |
6274ddfc | 318 | frateMixH.SetParameters(pT,y,tc,iProcHHG); |
e0e1d4c6 | 319 | frateMixH.SetParNames("transverse momentum","rapidity","critical temperature","process number"); |
320 | return TMath::TwoPi()*pT*lamHHG*frateMixH.Integral(-yNucl,yNucl); | |
321 | } | |
322 | ||
323 | // ----------------------------------------------------------------------------------------------------- | |
324 | static Double_t rateHHG(Double_t *x, Double_t *par) { | |
325 | //--------------------------------------------------- | |
326 | // input: | |
327 | // x[0] - tau (fm), proper time | |
328 | // x[1] - yprime, space rapidity | |
329 | // par[0] - p_T (GeV), photon transverse momentum | |
330 | // par[1] - y, photon rapidity in the c.m.s. A+A | |
331 | // par[2] - tauCHHG (fm), start of HHG | |
332 | // par[3] - T_c (GeV), critical temperature | |
333 | // par[4] - iProcHHG, process number | |
334 | // | |
335 | // output: | |
336 | // EdR/d^3p = EdN/d^4xd^3p (fm^{-4}GeV^{-2}) | |
337 | //--------------------------------------------------- | |
338 | Double_t tau=x[0],yprime=x[1]; | |
6274ddfc | 339 | Double_t pT=par[0],y=par[1],tauCHHG=par[2],tc=par[3]; |
e0e1d4c6 | 340 | Int_t iProcHHG=Int_t(par[4]); |
341 | ||
6274ddfc | 342 | Double_t e=pT*TMath::CosH(yprime-y),t=tc*TMath::Power(tauCHHG/tau,1./3.); |
e0e1d4c6 | 343 | const Double_t mPi=0.135; |
6274ddfc | 344 | Double_t xx=t/mPi,yy=e/mPi; |
345 | Double_t f,rate=1.,emin,factor; | |
346 | const Double_t mOm=0.783,width=0.00844,br=0.085,e0=0.379,pi=TMath::Pi(); | |
e0e1d4c6 | 347 | const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm |
348 | ||
349 | switch (iProcHHG) { | |
350 | ||
351 | case 0: | |
352 | // pi rho --> pi gamma | |
6274ddfc | 353 | f=-2.447+0.796*xx+(0.0338+0.0528*xx)*yy+(-21.447+8.2179*xx)/yy+(1.52436-0.38562*xx)/(yy*yy); |
354 | rate=t*t*TMath::Exp(-e/t+f); | |
e0e1d4c6 | 355 | break ; |
356 | ||
357 | case 1: | |
358 | // pi pi --> rho gamma | |
6274ddfc | 359 | f=-12.055+4.387*xx+(0.3755+0.00826*xx)*yy+(-0.00777+0.000279*xx)*yy*yy+(5.7869-1.0258*xx)/yy+(-1.979+0.58*xx)/(yy*yy); |
360 | rate=t*t*TMath::Exp(-e/t+f); | |
e0e1d4c6 | 361 | break ; |
362 | ||
363 | case 2: | |
364 | // rho --> pi pi gamma | |
6274ddfc | 365 | f=-6.295+1.6459*xx+(-0.4015+0.089*xx)*yy+(-0.954+2.05777*xx)/yy; |
366 | rate=t*t*TMath::Exp(-e/t+f); | |
e0e1d4c6 | 367 | break ; |
368 | ||
369 | case 3: | |
370 | // omega --> pi gamma | |
6274ddfc | 371 | emin=mOm*(e*e+e0*e0)/(2.*e*e0); |
372 | factor=(3.*mOm*width*br)/(16.*pi*pi*pi*e0); | |
373 | rate=factor*t*(emin+t)*TMath::Exp(-emin/t)/e/hc4; | |
e0e1d4c6 | 374 | break ; |
375 | ||
376 | default: | |
377 | printf("NO iProcHHG=%i \n",iProcHHG); | |
378 | } | |
379 | return tau*rate; | |
380 | } | |
381 | ||
382 | // ----------------------------------------------------------------------------------------------------- | |
383 | static Double_t fromHHG(Double_t *x, Double_t *par) { | |
384 | // Thermal photon spectrum from Hot Hadron Gas (HHG) | |
385 | // F.D.Steffen, nucl-th/9909035 | |
386 | // T.Peitzmann and M.H.Thoma, Phys.Rep., 364, 175 (2002), section 2.2.2 | |
387 | //--------------------------------------------------- | |
388 | // input: | |
389 | // x[0] - p_T (GeV), photon p_T | |
390 | // par[0] - tauCHHG (fm), start of HHG | |
391 | // par[1] - tauF (fm), freeze-out proper time | |
392 | // par[2] - yNucl, rapidity of projectile nucleus | |
393 | // par[3] - T_c (GeV), critical temperature | |
394 | // par[4] - y, photon rapidity | |
395 | // par[5] - iProcHHG, process number | |
396 | // | |
397 | // output: | |
398 | // d^{2}N/(dp_t dy) (1/GeV) | |
399 | //--------------------------------------------------- | |
400 | Double_t pT=x[0]; | |
6274ddfc | 401 | Double_t tauCHHG=par[0],tauF=par[1],yNucl=par[2],tc=par[3],y=par[4],iProcHHG=par[5]; |
e0e1d4c6 | 402 | |
403 | TF2 frateHHG("frateHHG",&rateHHG,tauCHHG,tauF,-yNucl,yNucl,5); | |
6274ddfc | 404 | frateHHG.SetParameters(pT,y,tauCHHG,tc,iProcHHG); |
e0e1d4c6 | 405 | frateHHG.SetParNames("transverse momentum","rapidity","start of HHG","criti temperature","process number"); |
406 | return TMath::TwoPi()*pT*frateHHG.Integral(tauCHHG,tauF,-yNucl,yNucl,1e-6); | |
407 | } | |
408 | ||
409 | // ----------------------------------------------------------------------------------------------------- | |
410 | static Double_t fOverlapAB(Double_t *x, Double_t *par) | |
411 | { | |
412 | //------------------------------------------------------------------------- | |
413 | // overlap area at the impact parameter b | |
414 | // input: | |
415 | // x[0] - impact parameter b < RA+RB | |
416 | // par[0] - radius of A | |
417 | // par[1] - radius of B | |
418 | //------------------------------------------------------------------------- | |
419 | ||
6274ddfc | 420 | Double_t b=x[0], ra=par[0], rb=par[1]; |
421 | if(rb>ra) {ra=par[1]; rb=par[0];} // ra > rb | |
e0e1d4c6 | 422 | |
6274ddfc | 423 | if(b>(ra+rb)) { |
e0e1d4c6 | 424 | return 0.; |
425 | } | |
426 | ||
6274ddfc | 427 | if(b<=(ra-rb)) { |
428 | return TMath::Pi()*rb*rb; | |
e0e1d4c6 | 429 | } |
430 | ||
6274ddfc | 431 | Double_t p=0.5*(b+ra+rb), S=TMath::Sqrt(p*(p-b)*(p-ra)*(p-rb)), h=2.*S/b; |
432 | Double_t sA=ra*ra*TMath::ASin(h/ra)-h*TMath::Sqrt(ra*ra-h*h); | |
433 | Double_t sB=rb*rb*TMath::ASin(h/rb)-h*TMath::Sqrt(rb*rb-h*h); | |
434 | if(ra>rb && b*b<ra*ra-rb*rb) sB=TMath::Pi()*rb*rb-sB; | |
e0e1d4c6 | 435 | |
436 | return sA+sB; | |
437 | ||
438 | } | |
439 | ||
440 | //_____________________________________________________________________________ | |
441 | AliGenThermalPhotons::AliGenThermalPhotons() | |
442 | :AliGenerator(-1), | |
e0e1d4c6 | 443 | fMinImpactParam(0.), |
444 | fMaxImpactParam(0.), | |
445 | fTau0(0.), | |
446 | fT0(0.), | |
447 | fTc(0.), | |
448 | fTf(0.), | |
449 | fGhhg(0), | |
450 | fSumPt() | |
451 | { | |
452 | // | |
453 | // Default constructor | |
454 | // | |
455 | SetCutVertexZ(); | |
456 | SetPtRange(); | |
457 | SetYRange(); | |
6274ddfc | 458 | fProjectile = 208; |
459 | fTarget = 208; | |
460 | fEnergyCMS = 5500.; | |
e0e1d4c6 | 461 | } |
462 | ||
463 | //_____________________________________________________________________________ | |
464 | AliGenThermalPhotons::AliGenThermalPhotons(Int_t npart) | |
465 | :AliGenerator(npart), | |
e0e1d4c6 | 466 | fMinImpactParam(0.), |
467 | fMaxImpactParam(0.), | |
468 | fTau0(0.1), | |
469 | fT0(0.650), | |
470 | fTc(0.170), | |
471 | fTf(0.100), | |
472 | fGhhg(8), | |
473 | fSumPt() | |
474 | { | |
475 | // | |
476 | // Standard constructor | |
477 | // | |
478 | ||
479 | fName="ThermalPhotons"; | |
480 | fTitle="Direct thermal photons in 1+1 Bjorken hydrodynamics"; | |
481 | ||
482 | SetCutVertexZ(); | |
483 | SetPtRange(); | |
484 | SetYRange(); | |
6274ddfc | 485 | fProjectile = 208; |
486 | fTarget = 208; | |
487 | fEnergyCMS = 5500.; | |
e0e1d4c6 | 488 | } |
489 | ||
490 | //_____________________________________________________________________________ | |
491 | AliGenThermalPhotons::~AliGenThermalPhotons() | |
492 | { | |
493 | // | |
494 | // Standard destructor | |
495 | // | |
496 | delete fSumPt; | |
497 | } | |
498 | ||
499 | //_____________________________________________________________________________ | |
500 | void AliGenThermalPhotons::Init() | |
501 | { | |
502 | ||
503 | const Double_t step=0.1; | |
504 | Int_t nPt=Int_t((fPtMax-fPtMin)/step); | |
505 | ||
506 | fSumPt = new TH1F("fSumPt","thermal #gamma from QGP",nPt,fPtMin,fPtMax); | |
507 | ||
508 | Double_t yRap=0.; | |
509 | const Int_t nCo=3,nFl=3; // number of colors for QGP | |
510 | Double_t gQGP=2.*(nCo*nCo-1.)+(7./8.)*4.*nCo*nFl; // number of degrees of freedom in QGP | |
511 | Double_t yNucl=TMath::ACosH(fEnergyCMS/2.); | |
512 | Double_t tauCQGP=TMath::Power(fT0/fTc,3.)*fTau0,tauCHHG=gQGP*tauCQGP/fGhhg,tauF=tauCHHG*TMath::Power(fTc/fTf,3.); | |
513 | Double_t lambda1=tauCQGP*(gQGP/(gQGP-fGhhg)),lambda2=-fGhhg/(gQGP-fGhhg); | |
514 | Double_t lamQGP=(tauCHHG-tauCQGP)*(lambda1+0.5*lambda2*(tauCHHG+tauCQGP)),lamHHG=0.5*(tauCHHG-tauCQGP)*(tauCHHG+tauCQGP)-lamQGP; | |
515 | ||
516 | Double_t pt,w; | |
517 | // photons from pure QGP phase | |
518 | for(Int_t j=0; j<3; j++) { | |
519 | TF1 func("func",&fromQGP,fPtMin,fPtMax,7); | |
520 | func.SetParameters(fTau0,fT0,tauCQGP,yNucl,fTc,yRap,j); | |
521 | func.SetParNames("nuclear radius","initial time","initial temperature","end of pure QGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number"); | |
522 | for(Int_t i=0; i<nPt; i++) { | |
523 | pt=fPtMin+(i+0.5)*step; | |
524 | w=func.Eval(pt); | |
525 | fSumPt->AddBinContent(i+1,w); | |
526 | } | |
527 | } | |
528 | ||
529 | // photons from mixed QGP phase | |
530 | for(Int_t j=0; j<3; j++) { | |
531 | TF1 func("func",&fromMixQ,fPtMin,fPtMax,5); | |
532 | func.SetParameters(lamQGP,yNucl,fTc,yRap,j); | |
533 | func.SetParNames("lamQGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number"); | |
534 | for(Int_t i=0; i<nPt; i++) { | |
535 | pt=fPtMin+(i+0.5)*step; | |
536 | w=func.Eval(pt); | |
537 | fSumPt->AddBinContent(i+1,w); | |
538 | } | |
539 | } | |
540 | ||
541 | // photons from mixed HHG phase | |
542 | for(Int_t j=0; j<4; j++) { | |
543 | TF1 func("func",&fromMixH,fPtMin,fPtMax,5); | |
544 | func.SetParameters(lamHHG,yNucl,fTc,yRap,j); | |
545 | func.SetParNames("lamQGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number"); | |
546 | for(Int_t i=0; i<nPt; i++) { | |
547 | pt=fPtMin+(i+0.5)*step; | |
548 | w=func.Eval(pt); | |
549 | fSumPt->AddBinContent(i+1,w); | |
550 | } | |
551 | } | |
552 | ||
553 | // photons from pure HHG phase | |
554 | for(Int_t j=0; j<4; j++) { | |
555 | TF1 func("func",&fromHHG,fPtMin,fPtMax,6); | |
556 | func.SetParameters(tauCHHG,tauF,yNucl,fTc,yRap,j); | |
557 | func.SetParNames("nuclear radius","start HHG","freeze-out proper time","rapidity of projectile nucleus","critical temperature","photon rapidity","process number"); | |
558 | for(Int_t i=0; i<nPt; i++) { | |
559 | pt=fPtMin+(i+0.5)*step; | |
560 | w=func.Eval(pt); | |
561 | fSumPt->AddBinContent(i+1,w); | |
562 | } | |
563 | } | |
564 | ||
565 | } | |
566 | ||
567 | //_____________________________________________________________________________ | |
568 | void AliGenThermalPhotons::Generate() | |
569 | { | |
570 | // | |
571 | // Generate thermal photons of a event | |
572 | // | |
573 | ||
574 | Float_t polar[3]= {0,0,0}; | |
575 | Float_t origin[3]; | |
576 | Float_t p[3]; | |
577 | Float_t random[6]; | |
578 | Int_t nt; | |
579 | ||
580 | for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j]; | |
581 | /* | |
582 | if(fVertexSmear==kPerEvent) { | |
583 | Vertex(); | |
584 | for (j=0;j<3;j++) origin[j]=fVertex[j]; | |
585 | } | |
586 | */ | |
587 | TArrayF eventVertex; | |
588 | eventVertex.Set(3); | |
589 | eventVertex[0] = origin[0]; | |
590 | eventVertex[1] = origin[1]; | |
591 | eventVertex[2] = origin[2]; | |
592 | ||
593 | Int_t nGam; | |
594 | Float_t impPar,area,pt,rapidity,phi,ww; | |
6274ddfc | 595 | Float_t r0=1.3,ra=r0*TMath::Power(fAProjectile,1./3.),rb=r0*TMath::Power(fATarget,1./3.); |
e0e1d4c6 | 596 | |
597 | TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2); | |
6274ddfc | 598 | funcOL->SetParameters(ra,rb); |
e0e1d4c6 | 599 | funcOL->SetParNames("radiusA","radiusB"); |
600 | ||
601 | impPar=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm()); | |
602 | area=funcOL->Eval(impPar); | |
603 | ||
604 | ww=area*(fYMax-fYMin)*fSumPt->GetBinWidth(1)*fSumPt->GetSumOfWeights(); | |
605 | nGam=Int_t(ww); | |
606 | if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++; | |
607 | ||
608 | if(nGam) { | |
609 | for(Int_t i=0; i<nGam; i++) { | |
610 | pt=fSumPt->GetRandom(); | |
611 | Rndm(random,2); | |
612 | rapidity=(fYMax-fYMin)*random[0]+fYMin; | |
613 | phi=2.*TMath::Pi()*random[1]; | |
614 | p[0]=pt*TMath::Cos(phi); | |
615 | p[1]=pt*TMath::Sin(phi); | |
616 | p[2]=pt*TMath::SinH(rapidity); | |
617 | ||
618 | if(fVertexSmear==kPerTrack) { | |
619 | Rndm(random,6); | |
620 | for (Int_t j=0;j<3;j++) { | |
621 | origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* | |
622 | TMath::Sqrt(-2*TMath::Log(random[2*j+1])); | |
623 | } | |
624 | } | |
625 | ||
626 | PushTrack(fTrackIt,-1,22,p,origin,polar,0,kPPrimary,nt,1.); | |
627 | } | |
628 | } | |
629 | ||
630 | delete funcOL; | |
631 | // Header | |
632 | AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons"); | |
633 | // Event Vertex | |
634 | header->SetPrimaryVertex(eventVertex); | |
635 | header->SetNProduced(fNpart); | |
636 | gAlice->SetGenEventHeader(header); | |
637 | ||
638 | } | |
639 | ||
640 | void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) { | |
641 | AliGenerator::SetPtRange(ptmin, ptmax); | |
642 | } | |
643 | ||
644 | void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) { | |
645 | AliGenerator::SetYRange(ymin, ymax); | |
646 | } |