Line data Source code
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 6 : ClassImp(AliGenThermalPhotons)
87 :
88 : // -----------------------------------------------------------------------------------------------------
89 : static Double_t rateQGP(const Double_t *x, const 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 0 : Double_t tau=x[0],yprime=x[1];
105 0 : Double_t pT=par[0],y=par[1],tau0=par[2],t0=par[3],tc=par[4];
106 0 : Int_t iProcQGP=Int_t(par[5]), nFl=3;
107 :
108 0 : Double_t e=pT*TMath::CosH(yprime-y),t=t0*TMath::Power(tau0/tau,1./3.);
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;
113 0 : Double_t alphaS=3.*TMath::TwoPi()/((33.-2.*nFl)*TMath::Log(8.*t/tc));
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 0 : switch (iProcQGP) {
118 :
119 : case 0:
120 : // 1-loop
121 0 : rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t);
122 0 : break ;
123 :
124 : case 1:
125 : // 2-loop: bremsstrahlung
126 0 : rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t;
127 0 : break ;
128 :
129 : case 2:
130 : // 2-loop: annihilation with scattering
131 0 : rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e;
132 0 : break ;
133 :
134 : default:
135 0 : printf("NO iProcQGP=%i \n",iProcQGP);
136 0 : }
137 :
138 0 : return tau*rate;
139 : }
140 :
141 : // -----------------------------------------------------------------------------------------------------
142 : static Double_t fromQGP(const Double_t *x, const 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 0 : Double_t pT=x[0];
158 0 : Double_t tau0=par[0],t0=par[1],tauCQGP=par[2],yNucl=par[3],tc=par[4],y=par[5];
159 0 : Int_t iProcQGP=Int_t(par[6]);
160 :
161 0 : TF2 frateQGP("frateQGP",&rateQGP,tau0,tauCQGP,-yNucl,yNucl,6);
162 0 : frateQGP.SetParameters(pT,y,tau0,t0,tc,iProcQGP);
163 0 : frateQGP.SetParNames("transverse momentum","rapidity","initial time","initial temperature","critical temperature","process number");
164 0 : return TMath::TwoPi()*pT*frateQGP.Integral(tau0,tauCQGP,-yNucl,yNucl,1e-6);
165 0 : }
166 :
167 : // -----------------------------------------------------------------------------------------------------
168 : static Double_t rateMixQ(const Double_t *x, const 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 0 : Double_t yprime=x[0];
181 0 : Double_t pT=par[0],y=par[1],tc=par[2];
182 0 : Int_t iProcQGP=Int_t(par[3]),nFl=3;
183 :
184 0 : Double_t e=pT*TMath::CosH(yprime-y),t=tc;
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;
189 0 : Double_t alphaS=3.*TMath::TwoPi()/((33.-2.*nFl)*TMath::Log(8.*t/tc));
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 0 : switch (iProcQGP) {
194 :
195 : case 0:
196 : // 1-loop
197 0 : rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t);
198 0 : break ;
199 :
200 : case 1:
201 : // 2-loop: bremsstrahlung
202 0 : rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t;
203 0 : break ;
204 :
205 : case 2:
206 : // 2-loop: annihilation with scattering
207 0 : rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e;
208 0 : break ;
209 :
210 : default:
211 0 : printf("NO iProcQGP=%i \n",iProcQGP);
212 0 : }
213 :
214 0 : return rate;
215 : }
216 :
217 : // -----------------------------------------------------------------------------------------------------
218 : static Double_t fromMixQ(const Double_t *x, const 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 0 : Double_t pT=x[0];
232 0 : Double_t lamQGP=par[0],yNucl=par[1],tc=par[2],y=par[3];
233 0 : Int_t iProcQGP=Int_t(par[4]);
234 :
235 0 : TF1 frateMixQ("frateMixQ",&rateMixQ,-yNucl,yNucl,4);
236 0 : frateMixQ.SetParameters(pT,y,tc,iProcQGP);
237 0 : frateMixQ.SetParNames("transverse momentum","rapidity","critical temperature","process number");
238 0 : return TMath::TwoPi()*pT*lamQGP*frateMixQ.Integral(-yNucl,yNucl);
239 0 : }
240 :
241 : // -----------------------------------------------------------------------------------------------------
242 : static Double_t rateMixH(const Double_t *x, const 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 0 : Double_t yprime=x[0];
255 0 : Double_t pT=par[0],y=par[1],tc=par[2];
256 0 : Int_t iProcHHG=Int_t(par[3]);
257 :
258 0 : Double_t e=pT*TMath::CosH(yprime-y),t=tc;
259 : const Double_t mPi=0.135;
260 0 : Double_t xx=t/mPi,yy=e/mPi;
261 : Double_t f,rate=1.,emin,factor;
262 0 : const Double_t mOm=0.783,width=0.00844,br=0.085,e0=0.379,pi=TMath::Pi();
263 : const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm
264 :
265 0 : switch (iProcHHG) {
266 :
267 : case 0:
268 : // pi rho --> pi gamma
269 0 : 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 0 : rate=t*t*TMath::Exp(-e/t+f);
271 0 : break ;
272 :
273 : case 1:
274 : // pi pi --> rho gamma
275 0 : 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 0 : rate=t*t*TMath::Exp(-e/t+f);
277 0 : break ;
278 :
279 : case 2:
280 : // rho --> pi pi gamma
281 0 : f=-6.295+1.6459*xx+(-0.4015+0.089*xx)*yy+(-0.954+2.05777*xx)/yy;
282 0 : rate=t*t*TMath::Exp(-e/t+f);
283 0 : break ;
284 :
285 : case 3:
286 : // omega --> pi gamma
287 0 : emin=mOm*(e*e+e0*e0)/(2.*e*e0);
288 0 : factor=(3.*mOm*width*br)/(16.*pi*pi*pi*e0);
289 0 : rate=factor*t*(emin+t)*TMath::Exp(-emin/t)/e/hc4;
290 0 : break ;
291 :
292 : default:
293 0 : printf("NO iProcHHG=%i \n",iProcHHG);
294 0 : }
295 :
296 0 : return rate;
297 : }
298 :
299 : // -----------------------------------------------------------------------------------------------------
300 : static Double_t fromMixH(const Double_t *x, const 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 0 : Double_t pT=x[0];
314 0 : Double_t lamHHG=par[0],yNucl=par[1],tc=par[2],y=par[3];
315 0 : Int_t iProcHHG=Int_t(par[4]);
316 :
317 0 : TF1 frateMixH("frateMixH",&rateMixH,-yNucl,yNucl,4);
318 0 : frateMixH.SetParameters(pT,y,tc,iProcHHG);
319 0 : frateMixH.SetParNames("transverse momentum","rapidity","critical temperature","process number");
320 0 : return TMath::TwoPi()*pT*lamHHG*frateMixH.Integral(-yNucl,yNucl);
321 0 : }
322 :
323 : // -----------------------------------------------------------------------------------------------------
324 : static Double_t rateHHG(const Double_t *x, const 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 0 : Double_t tau=x[0],yprime=x[1];
339 0 : Double_t pT=par[0],y=par[1],tauCHHG=par[2],tc=par[3];
340 0 : Int_t iProcHHG=Int_t(par[4]);
341 :
342 0 : Double_t e=pT*TMath::CosH(yprime-y),t=tc*TMath::Power(tauCHHG/tau,1./3.);
343 : const Double_t mPi=0.135;
344 0 : Double_t xx=t/mPi,yy=e/mPi;
345 : Double_t f,rate=1.,emin,factor;
346 0 : const Double_t mOm=0.783,width=0.00844,br=0.085,e0=0.379,pi=TMath::Pi();
347 : const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm
348 :
349 0 : switch (iProcHHG) {
350 :
351 : case 0:
352 : // pi rho --> pi gamma
353 0 : 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 0 : rate=t*t*TMath::Exp(-e/t+f);
355 0 : break ;
356 :
357 : case 1:
358 : // pi pi --> rho gamma
359 0 : 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 0 : rate=t*t*TMath::Exp(-e/t+f);
361 0 : break ;
362 :
363 : case 2:
364 : // rho --> pi pi gamma
365 0 : f=-6.295+1.6459*xx+(-0.4015+0.089*xx)*yy+(-0.954+2.05777*xx)/yy;
366 0 : rate=t*t*TMath::Exp(-e/t+f);
367 0 : break ;
368 :
369 : case 3:
370 : // omega --> pi gamma
371 0 : emin=mOm*(e*e+e0*e0)/(2.*e*e0);
372 0 : factor=(3.*mOm*width*br)/(16.*pi*pi*pi*e0);
373 0 : rate=factor*t*(emin+t)*TMath::Exp(-emin/t)/e/hc4;
374 0 : break ;
375 :
376 : default:
377 0 : printf("NO iProcHHG=%i \n",iProcHHG);
378 0 : }
379 0 : return tau*rate;
380 : }
381 :
382 : // -----------------------------------------------------------------------------------------------------
383 : static Double_t fromHHG(const Double_t *x, const 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 0 : Double_t pT=x[0];
401 0 : Double_t tauCHHG=par[0],tauF=par[1],yNucl=par[2],tc=par[3],y=par[4],iProcHHG=par[5];
402 :
403 0 : TF2 frateHHG("frateHHG",&rateHHG,tauCHHG,tauF,-yNucl,yNucl,5);
404 0 : frateHHG.SetParameters(pT,y,tauCHHG,tc,iProcHHG);
405 0 : frateHHG.SetParNames("transverse momentum","rapidity","start of HHG","criti temperature","process number");
406 0 : return TMath::TwoPi()*pT*frateHHG.Integral(tauCHHG,tauF,-yNucl,yNucl,1e-6);
407 0 : }
408 :
409 : // -----------------------------------------------------------------------------------------------------
410 : static Double_t fOverlapAB(const Double_t *x, const 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 :
420 0 : Double_t b=x[0], ra=par[0], rb=par[1];
421 0 : if(rb>ra) {ra=par[1]; rb=par[0];} // ra > rb
422 :
423 0 : if(b>(ra+rb)) {
424 0 : return 0.;
425 : }
426 :
427 0 : if(b<=(ra-rb)) {
428 0 : return TMath::Pi()*rb*rb;
429 : }
430 :
431 0 : Double_t p=0.5*(b+ra+rb), S=TMath::Sqrt(p*(p-b)*(p-ra)*(p-rb)), h=2.*S/b;
432 0 : Double_t sA=ra*ra*TMath::ASin(h/ra)-h*TMath::Sqrt(ra*ra-h*h);
433 0 : Double_t sB=rb*rb*TMath::ASin(h/rb)-h*TMath::Sqrt(rb*rb-h*h);
434 0 : if(ra>rb && b*b<ra*ra-rb*rb) sB=TMath::Pi()*rb*rb-sB;
435 :
436 0 : return sA+sB;
437 :
438 0 : }
439 :
440 : //_____________________________________________________________________________
441 : AliGenThermalPhotons::AliGenThermalPhotons()
442 0 : :AliGenerator(-1),
443 0 : fMinImpactParam(0.),
444 0 : fMaxImpactParam(0.),
445 0 : fTau0(0.),
446 0 : fT0(0.),
447 0 : fTc(0.),
448 0 : fTf(0.),
449 0 : fGhhg(0),
450 0 : fSumPt()
451 0 : {
452 : //
453 : // Default constructor
454 : //
455 0 : SetCutVertexZ();
456 0 : SetPtRange();
457 0 : SetYRange();
458 0 : fAProjectile = 208;
459 0 : fATarget = 208;
460 0 : fEnergyCMS = 5500.;
461 0 : }
462 :
463 : //_____________________________________________________________________________
464 : AliGenThermalPhotons::AliGenThermalPhotons(Int_t npart)
465 0 : :AliGenerator(npart),
466 0 : fMinImpactParam(0.),
467 0 : fMaxImpactParam(0.),
468 0 : fTau0(0.1),
469 0 : fT0(0.650),
470 0 : fTc(0.170),
471 0 : fTf(0.100),
472 0 : fGhhg(8),
473 0 : fSumPt()
474 0 : {
475 : //
476 : // Standard constructor
477 : //
478 :
479 0 : fName="ThermalPhotons";
480 0 : fTitle="Direct thermal photons in 1+1 Bjorken hydrodynamics";
481 :
482 0 : SetCutVertexZ();
483 0 : SetPtRange();
484 0 : SetYRange();
485 0 : fAProjectile = 208;
486 0 : fATarget = 208;
487 0 : fEnergyCMS = 5500.;
488 0 : }
489 :
490 : //_____________________________________________________________________________
491 0 : AliGenThermalPhotons::~AliGenThermalPhotons()
492 0 : {
493 : //
494 : // Standard destructor
495 : //
496 0 : delete fSumPt;
497 0 : }
498 :
499 : //_____________________________________________________________________________
500 : void AliGenThermalPhotons::Init()
501 : {
502 : // Initialisation
503 : const Double_t step=0.1;
504 0 : Int_t nPt=Int_t((fPtMax-fPtMin)/step);
505 :
506 0 : 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 0 : Double_t yNucl=TMath::ACosH(fEnergyCMS/2.);
512 0 : Double_t tauCQGP=TMath::Power(fT0/fTc,3.)*fTau0,tauCHHG=gQGP*tauCQGP/fGhhg,tauF=tauCHHG*TMath::Power(fTc/fTf,3.);
513 0 : Double_t lambda1=tauCQGP*(gQGP/(gQGP-fGhhg)),lambda2=-fGhhg/(gQGP-fGhhg);
514 0 : 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 0 : for(Int_t j=0; j<3; j++) {
519 0 : TF1 func("func",&fromQGP,fPtMin,fPtMax,7);
520 0 : func.SetParameters(fTau0,fT0,tauCQGP,yNucl,fTc,yRap,j);
521 0 : func.SetParNames("nuclear radius","initial time","initial temperature","end of pure QGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
522 0 : for(Int_t i=0; i<nPt; i++) {
523 0 : pt=fPtMin+(i+0.5)*step;
524 0 : w=func.Eval(pt);
525 0 : fSumPt->AddBinContent(i+1,w);
526 : }
527 0 : }
528 :
529 : // photons from mixed QGP phase
530 0 : for(Int_t j=0; j<3; j++) {
531 0 : TF1 func("func",&fromMixQ,fPtMin,fPtMax,5);
532 0 : func.SetParameters(lamQGP,yNucl,fTc,yRap,j);
533 0 : func.SetParNames("lamQGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
534 0 : for(Int_t i=0; i<nPt; i++) {
535 0 : pt=fPtMin+(i+0.5)*step;
536 0 : w=func.Eval(pt);
537 0 : fSumPt->AddBinContent(i+1,w);
538 : }
539 0 : }
540 :
541 : // photons from mixed HHG phase
542 0 : for(Int_t j=0; j<4; j++) {
543 0 : TF1 func("func",&fromMixH,fPtMin,fPtMax,5);
544 0 : func.SetParameters(lamHHG,yNucl,fTc,yRap,j);
545 0 : func.SetParNames("lamQGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
546 0 : for(Int_t i=0; i<nPt; i++) {
547 0 : pt=fPtMin+(i+0.5)*step;
548 0 : w=func.Eval(pt);
549 0 : fSumPt->AddBinContent(i+1,w);
550 : }
551 0 : }
552 :
553 : // photons from pure HHG phase
554 0 : for(Int_t j=0; j<4; j++) {
555 0 : TF1 func("func",&fromHHG,fPtMin,fPtMax,6);
556 0 : func.SetParameters(tauCHHG,tauF,yNucl,fTc,yRap,j);
557 0 : func.SetParNames("nuclear radius","start HHG","freeze-out proper time","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
558 0 : for(Int_t i=0; i<nPt; i++) {
559 0 : pt=fPtMin+(i+0.5)*step;
560 0 : w=func.Eval(pt);
561 0 : fSumPt->AddBinContent(i+1,w);
562 : }
563 0 : }
564 :
565 0 : }
566 :
567 : //_____________________________________________________________________________
568 : void AliGenThermalPhotons::Generate()
569 : {
570 : //
571 : // Generate thermal photons of a event
572 : //
573 :
574 0 : Float_t polar[3]= {0,0,0};
575 0 : Float_t origin[3];
576 : Float_t time;
577 0 : Float_t p[3];
578 0 : Float_t random[6];
579 0 : Int_t nt;
580 :
581 0 : for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
582 0 : time = fTimeOrigin;
583 : /*
584 : if(fVertexSmear==kPerEvent) {
585 : Vertex();
586 : for (j=0;j<3;j++) origin[j]=fVertex[j];
587 : time = fTime;
588 : }
589 : */
590 0 : TArrayF eventVertex;
591 0 : eventVertex.Set(3);
592 0 : eventVertex[0] = origin[0];
593 0 : eventVertex[1] = origin[1];
594 0 : eventVertex[2] = origin[2];
595 : Float_t eventTime = time;
596 :
597 : Int_t nGam;
598 : Float_t impPar,area,pt,rapidity,phi,ww;
599 0 : Float_t r0=1.3,ra=r0*TMath::Power(fAProjectile,1./3.),rb=r0*TMath::Power(fATarget,1./3.);
600 :
601 0 : TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2);
602 0 : funcOL->SetParameters(ra,rb);
603 0 : funcOL->SetParNames("radiusA","radiusB");
604 :
605 0 : impPar=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
606 0 : area=funcOL->Eval(impPar);
607 :
608 0 : ww=area*(fYMax-fYMin)*fSumPt->GetBinWidth(1)*fSumPt->GetSumOfWeights();
609 0 : nGam=Int_t(ww);
610 0 : if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
611 :
612 0 : if(nGam) {
613 0 : for(Int_t i=0; i<nGam; i++) {
614 0 : pt=fSumPt->GetRandom();
615 0 : Rndm(random,2);
616 0 : rapidity=(fYMax-fYMin)*random[0]+fYMin;
617 0 : phi=2.*TMath::Pi()*random[1];
618 0 : p[0]=pt*TMath::Cos(phi);
619 0 : p[1]=pt*TMath::Sin(phi);
620 0 : p[2]=pt*TMath::SinH(rapidity);
621 :
622 0 : if(fVertexSmear==kPerTrack) {
623 0 : Rndm(random,6);
624 0 : for (Int_t j=0;j<3;j++) {
625 0 : origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
626 0 : TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
627 : }
628 0 : Rndm(random,2);
629 0 : time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
630 0 : TMath::Cos(2*random[0]*TMath::Pi())*
631 0 : TMath::Sqrt(-2*TMath::Log(random[1]));
632 0 : }
633 :
634 0 : PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.);
635 : }
636 0 : }
637 :
638 0 : delete funcOL;
639 : // Header
640 0 : AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons");
641 : // Event Vertex
642 0 : header->SetPrimaryVertex(eventVertex);
643 0 : header->SetInteractionTime(eventTime);
644 0 : header->SetNProduced(fNpart);
645 0 : gAlice->SetGenEventHeader(header);
646 :
647 0 : }
648 :
649 : void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
650 0 : AliGenerator::SetPtRange(ptmin, ptmax);
651 0 : }
652 :
653 : void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) {
654 0 : AliGenerator::SetYRange(ymin, ymax);
655 0 : }
|
