1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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)
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
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)
32 // photon rates for QGP: Phys.Rep., 364, 175 (2002), section 2.1.1
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)
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
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 //-------------------------------------------------------------------------
62 <img src="picts/AliGeneratorClass.gif">
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>.
73 ///////////////////////////////////////////////////////////////////
81 #include "AliGenEventHeader.h"
82 #include "AliGenThermalPhotons.h"
86 ClassImp(AliGenThermalPhotons)
88 // -----------------------------------------------------------------------------------------------------
89 static Double_t rateQGP(Double_t *x, Double_t *par) {
90 //---------------------------------------------------
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
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];
105 Double_t pT=par[0],y=par[1],tau0=par[2],t0=par[3],tc=par[4];
106 Int_t iProcQGP=Int_t(par[5]), nFl=3;
108 Double_t e=pT*TMath::CosH(yprime-y),t=t0*TMath::Power(tau0/tau,1./3.);
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 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
121 rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t);
125 // 2-loop: bremsstrahlung
126 rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t;
130 // 2-loop: annihilation with scattering
131 rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e;
135 printf("NO iProcQGP=%i \n",iProcQGP);
141 // -----------------------------------------------------------------------------------------------------
142 static Double_t fromQGP(Double_t *x, Double_t *par) {
143 //---------------------------------------------------
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
155 // d^{2}N/(dp_t dy) (1/GeV)
156 //---------------------------------------------------
158 Double_t tau0=par[0],t0=par[1],tauCQGP=par[2],yNucl=par[3],tc=par[4],y=par[5];
159 Int_t iProcQGP=Int_t(par[6]);
161 TF2 frateQGP("frateQGP",&rateQGP,tau0,tauCQGP,-yNucl,yNucl,6);
162 frateQGP.SetParameters(pT,y,tau0,t0,tc,iProcQGP);
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);
167 // -----------------------------------------------------------------------------------------------------
168 static Double_t rateMixQ(Double_t *x, Double_t *par) {
169 //---------------------------------------------------
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
178 // EdR/d^3p = EdN/d^4xd^3p (fm fm^{-4}GeV^{-2})
179 //---------------------------------------------------
180 Double_t yprime=x[0];
181 Double_t pT=par[0],y=par[1],tc=par[2];
182 Int_t iProcQGP=Int_t(par[3]),nFl=3;
184 Double_t e=pT*TMath::CosH(yprime-y),t=tc;
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 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
197 rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t);
201 // 2-loop: bremsstrahlung
202 rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t;
206 // 2-loop: annihilation with scattering
207 rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e;
211 printf("NO iProcQGP=%i \n",iProcQGP);
217 // -----------------------------------------------------------------------------------------------------
218 static Double_t fromMixQ(Double_t *x, Double_t *par) {
219 //---------------------------------------------------
221 // x[0] - p_T (GeV), photon p_T
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
229 // d^{2}N/(dp_t dy) (1/GeV)
230 //---------------------------------------------------
232 Double_t lamQGP=par[0],yNucl=par[1],tc=par[2],y=par[3];
233 Int_t iProcQGP=Int_t(par[4]);
235 TF1 frateMixQ("frateMixQ",&rateMixQ,-yNucl,yNucl,4);
236 frateMixQ.SetParameters(pT,y,tc,iProcQGP);
237 frateMixQ.SetParNames("transverse momentum","rapidity","critical temperature","process number");
238 return TMath::TwoPi()*pT*lamQGP*frateMixQ.Integral(-yNucl,yNucl);
241 // -----------------------------------------------------------------------------------------------------
242 static Double_t rateMixH(Double_t *x, Double_t *par) {
243 //---------------------------------------------------
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
252 // EdR/d^3p = EdN/d^4xd^3p (fm^{-4}GeV^{-2})
253 //---------------------------------------------------
254 Double_t yprime=x[0];
255 Double_t pT=par[0],y=par[1],tc=par[2];
256 Int_t iProcHHG=Int_t(par[3]);
258 Double_t e=pT*TMath::CosH(yprime-y),t=tc;
259 const Double_t mPi=0.135;
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();
263 const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm
268 // pi rho --> pi gamma
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);
274 // pi pi --> rho gamma
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);
280 // rho --> pi pi gamma
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);
286 // omega --> pi gamma
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;
293 printf("NO iProcHHG=%i \n",iProcHHG);
299 // -----------------------------------------------------------------------------------------------------
300 static Double_t fromMixH(Double_t *x, Double_t *par) {
301 //---------------------------------------------------
303 // x[0] - p_T (GeV), photon p_T
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
311 // d^{2}N/(dp_t dy) (1/GeV)
312 //---------------------------------------------------
314 Double_t lamHHG=par[0],yNucl=par[1],tc=par[2],y=par[3];
315 Int_t iProcHHG=Int_t(par[4]);
317 TF1 frateMixH("frateMixH",&rateMixH,-yNucl,yNucl,4);
318 frateMixH.SetParameters(pT,y,tc,iProcHHG);
319 frateMixH.SetParNames("transverse momentum","rapidity","critical temperature","process number");
320 return TMath::TwoPi()*pT*lamHHG*frateMixH.Integral(-yNucl,yNucl);
323 // -----------------------------------------------------------------------------------------------------
324 static Double_t rateHHG(Double_t *x, Double_t *par) {
325 //---------------------------------------------------
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
336 // EdR/d^3p = EdN/d^4xd^3p (fm^{-4}GeV^{-2})
337 //---------------------------------------------------
338 Double_t tau=x[0],yprime=x[1];
339 Double_t pT=par[0],y=par[1],tauCHHG=par[2],tc=par[3];
340 Int_t iProcHHG=Int_t(par[4]);
342 Double_t e=pT*TMath::CosH(yprime-y),t=tc*TMath::Power(tauCHHG/tau,1./3.);
343 const Double_t mPi=0.135;
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();
347 const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm
352 // pi rho --> pi gamma
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);
358 // pi pi --> rho gamma
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);
364 // rho --> pi pi gamma
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);
370 // omega --> pi gamma
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;
377 printf("NO iProcHHG=%i \n",iProcHHG);
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 //---------------------------------------------------
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
398 // d^{2}N/(dp_t dy) (1/GeV)
399 //---------------------------------------------------
401 Double_t tauCHHG=par[0],tauF=par[1],yNucl=par[2],tc=par[3],y=par[4],iProcHHG=par[5];
403 TF2 frateHHG("frateHHG",&rateHHG,tauCHHG,tauF,-yNucl,yNucl,5);
404 frateHHG.SetParameters(pT,y,tauCHHG,tc,iProcHHG);
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);
409 // -----------------------------------------------------------------------------------------------------
410 static Double_t fOverlapAB(Double_t *x, Double_t *par)
412 //-------------------------------------------------------------------------
413 // overlap area at the impact parameter b
415 // x[0] - impact parameter b < RA+RB
416 // par[0] - radius of A
417 // par[1] - radius of B
418 //-------------------------------------------------------------------------
420 Double_t b=x[0], ra=par[0], rb=par[1];
421 if(rb>ra) {ra=par[1]; rb=par[0];} // ra > rb
428 return TMath::Pi()*rb*rb;
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;
440 //_____________________________________________________________________________
441 AliGenThermalPhotons::AliGenThermalPhotons()
453 // Default constructor
463 //_____________________________________________________________________________
464 AliGenThermalPhotons::AliGenThermalPhotons(Int_t npart)
465 :AliGenerator(npart),
476 // Standard constructor
479 fName="ThermalPhotons";
480 fTitle="Direct thermal photons in 1+1 Bjorken hydrodynamics";
490 //_____________________________________________________________________________
491 AliGenThermalPhotons::~AliGenThermalPhotons()
494 // Standard destructor
499 //_____________________________________________________________________________
500 void AliGenThermalPhotons::Init()
503 const Double_t step=0.1;
504 Int_t nPt=Int_t((fPtMax-fPtMin)/step);
506 fSumPt = new TH1F("fSumPt","thermal #gamma from QGP",nPt,fPtMin,fPtMax);
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;
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;
525 fSumPt->AddBinContent(i+1,w);
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;
537 fSumPt->AddBinContent(i+1,w);
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;
549 fSumPt->AddBinContent(i+1,w);
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;
561 fSumPt->AddBinContent(i+1,w);
567 //_____________________________________________________________________________
568 void AliGenThermalPhotons::Generate()
571 // Generate thermal photons of a event
574 Float_t polar[3]= {0,0,0};
580 for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
582 if(fVertexSmear==kPerEvent) {
584 for (j=0;j<3;j++) origin[j]=fVertex[j];
589 eventVertex[0] = origin[0];
590 eventVertex[1] = origin[1];
591 eventVertex[2] = origin[2];
594 Float_t impPar,area,pt,rapidity,phi,ww;
595 Float_t r0=1.3,ra=r0*TMath::Power(fAProjectile,1./3.),rb=r0*TMath::Power(fATarget,1./3.);
597 TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2);
598 funcOL->SetParameters(ra,rb);
599 funcOL->SetParNames("radiusA","radiusB");
601 impPar=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
602 area=funcOL->Eval(impPar);
604 ww=area*(fYMax-fYMin)*fSumPt->GetBinWidth(1)*fSumPt->GetSumOfWeights();
606 if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
609 for(Int_t i=0; i<nGam; i++) {
610 pt=fSumPt->GetRandom();
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);
618 if(fVertexSmear==kPerTrack) {
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]));
626 PushTrack(fTrackIt,-1,22,p,origin,polar,0,kPPrimary,nt,1.);
632 AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons");
634 header->SetPrimaryVertex(eventVertex);
635 header->SetNProduced(fNpart);
636 gAlice->SetGenEventHeader(header);
640 void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
641 AliGenerator::SetPtRange(ptmin, ptmax);
644 void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) {
645 AliGenerator::SetYRange(ymin, ymax);