]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenThermalPhotons.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / EVGEN / AliGenThermalPhotons.cxx
CommitLineData
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
86ClassImp(AliGenThermalPhotons)
87
88// -----------------------------------------------------------------------------------------------------
4a33c50d 89static Double_t rateQGP(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 142static Double_t fromQGP(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 168static Double_t rateMixQ(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 218static Double_t fromMixQ(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 242static Double_t rateMixH(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 300static Double_t fromMixH(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 324static Double_t rateHHG(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 383static Double_t fromHHG(const Double_t *x, const Double_t *par) {
e0e1d4c6 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// -----------------------------------------------------------------------------------------------------
4a33c50d 410static Double_t fOverlapAB(const Double_t *x, const Double_t *par)
e0e1d4c6 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//_____________________________________________________________________________
441AliGenThermalPhotons::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();
eb35583a 458 fAProjectile = 208;
459 fATarget = 208;
6274ddfc 460 fEnergyCMS = 5500.;
e0e1d4c6 461}
462
463//_____________________________________________________________________________
464AliGenThermalPhotons::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();
eb35583a 485 fAProjectile = 208;
486 fATarget = 208;
6274ddfc 487 fEnergyCMS = 5500.;
e0e1d4c6 488}
489
490//_____________________________________________________________________________
491AliGenThermalPhotons::~AliGenThermalPhotons()
492{
493 //
494 // Standard destructor
495 //
496 delete fSumPt;
497}
498
499//_____________________________________________________________________________
500void AliGenThermalPhotons::Init()
501{
4a33c50d 502 // Initialisation
e0e1d4c6 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//_____________________________________________________________________________
568void 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];
21391258 576 Float_t time;
e0e1d4c6 577 Float_t p[3];
578 Float_t random[6];
579 Int_t nt;
580
581 for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
21391258 582 time = fTimeOrigin;
e0e1d4c6 583/*
584 if(fVertexSmear==kPerEvent) {
585 Vertex();
586 for (j=0;j<3;j++) origin[j]=fVertex[j];
21391258 587 time = fTime;
e0e1d4c6 588 }
589*/
590 TArrayF eventVertex;
591 eventVertex.Set(3);
592 eventVertex[0] = origin[0];
593 eventVertex[1] = origin[1];
594 eventVertex[2] = origin[2];
21391258 595 Float_t eventTime = time;
e0e1d4c6 596
597 Int_t nGam;
598 Float_t impPar,area,pt,rapidity,phi,ww;
6274ddfc 599 Float_t r0=1.3,ra=r0*TMath::Power(fAProjectile,1./3.),rb=r0*TMath::Power(fATarget,1./3.);
e0e1d4c6 600
601 TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2);
6274ddfc 602 funcOL->SetParameters(ra,rb);
e0e1d4c6 603 funcOL->SetParNames("radiusA","radiusB");
604
605 impPar=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
606 area=funcOL->Eval(impPar);
607
608 ww=area*(fYMax-fYMin)*fSumPt->GetBinWidth(1)*fSumPt->GetSumOfWeights();
609 nGam=Int_t(ww);
610 if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
611
612 if(nGam) {
613 for(Int_t i=0; i<nGam; i++) {
614 pt=fSumPt->GetRandom();
615 Rndm(random,2);
616 rapidity=(fYMax-fYMin)*random[0]+fYMin;
617 phi=2.*TMath::Pi()*random[1];
618 p[0]=pt*TMath::Cos(phi);
619 p[1]=pt*TMath::Sin(phi);
620 p[2]=pt*TMath::SinH(rapidity);
621
622 if(fVertexSmear==kPerTrack) {
623 Rndm(random,6);
624 for (Int_t j=0;j<3;j++) {
625 origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
626 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
627 }
21391258 628 Rndm(random,2);
629 time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
630 TMath::Cos(2*random[0]*TMath::Pi())*
631 TMath::Sqrt(-2*TMath::Log(random[1]));
e0e1d4c6 632 }
633
21391258 634 PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.);
e0e1d4c6 635 }
636 }
637
638 delete funcOL;
639// Header
640 AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons");
641// Event Vertex
642 header->SetPrimaryVertex(eventVertex);
21391258 643 header->SetInteractionTime(eventTime);
e0e1d4c6 644 header->SetNProduced(fNpart);
645 gAlice->SetGenEventHeader(header);
646
647}
648
649void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
650 AliGenerator::SetPtRange(ptmin, ptmax);
651}
652
653void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) {
654 AliGenerator::SetYRange(ymin, ymax);
655}