EVGEN/AliGenThermalPhotons.cxx

   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(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   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;
 107
 108   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   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   switch (iProcQGP) {
 118
 119     case 0:
 120 // 1-loop
 121       rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t);
 122     break ;
 123
 124     case 1:
 125 // 2-loop: bremsstrahlung
 126       rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t;
 127     break ;
 128
 129     case 2:
 130 // 2-loop: annihilation with scattering
 131       rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e;
 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(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   Double_t pT=x[0];
 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]);
 160
 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);
 165 }
 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   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;
 183
 184   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   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   switch (iProcQGP) {
 194
 195     case 0:
 196 // 1-loop
 197       rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t*TMath::Log(0.2317*e/alphaS/t);
 198     break ;
 199
 200     case 1:
 201 // 2-loop: bremsstrahlung
 202       rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*t;
 203     break ;
 204
 205     case 2:
 206 // 2-loop: annihilation with scattering
 207       rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-e/t)*t*e;
 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(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   Double_t pT=x[0];
 232   Double_t lamQGP=par[0],yNucl=par[1],tc=par[2],y=par[3];
 233   Int_t iProcQGP=Int_t(par[4]);
 234
 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);
 239 }
 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   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]);
 257
 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
 264
 265   switch (iProcHHG) {
 266
 267     case 0:
 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);
 271     break ;
 272
 273     case 1:
 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);
 277     break ;
 278
 279     case 2:
 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);
 283     break ;
 284
 285     case 3:
 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;
 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(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   Double_t pT=x[0];
 314   Double_t lamHHG=par[0],yNucl=par[1],tc=par[2],y=par[3];
 315   Int_t iProcHHG=Int_t(par[4]);
 316
 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);
 321 }
 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   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]);
 341
 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
 348
 349   switch (iProcHHG) {
 350
 351     case 0:
 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);
 355     break ;
 356
 357     case 1:
 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);
 361     break ;
 362
 363     case 2:
 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);
 367     break ;
 368
 369     case 3:
 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;
 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(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   Double_t pT=x[0];
 401   Double_t tauCHHG=par[0],tauF=par[1],yNucl=par[2],tc=par[3],y=par[4],iProcHHG=par[5];
 402
 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);
 407 }
 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   Double_t b=x[0], ra=par[0], rb=par[1];
 421   if(rb>ra) {ra=par[1]; rb=par[0];} // ra > rb
 422
 423   if(b>(ra+rb)) {
 424     return 0.;
 425   }
 426
 427   if(b<=(ra-rb)) {
 428     return TMath::Pi()*rb*rb;
 429   }
 430
 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;
 435
 436   return sA+sB;
 437
 438 }
 439
 440 //_____________________________________________________________________________
 441 AliGenThermalPhotons::AliGenThermalPhotons()
 442     :AliGenerator(-1),
 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();
 458     fAProjectile = 208;
 459     fATarget     = 208;
 460     fEnergyCMS  = 5500.;
 461 }
 462
 463 //_____________________________________________________________________________
 464 AliGenThermalPhotons::AliGenThermalPhotons(Int_t npart)
 465     :AliGenerator(npart),
 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();
 485     fAProjectile = 208;
 486     fATarget     = 208;
 487     fEnergyCMS  = 5500.;
 488 }
 489
 490 //_____________________________________________________________________________
 491 AliGenThermalPhotons::~AliGenThermalPhotons()
 492 {
 493   //
 494   // Standard destructor
 495   //
 496     delete fSumPt;
 497 }
 498
 499 //_____________________________________________________________________________
 500 void AliGenThermalPhotons::Init()
 501 {
 502     // Initialisation
 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 time;
 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];
 582     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     TArrayF eventVertex;
 591     eventVertex.Set(3);
 592     eventVertex[0] = origin[0];
 593     eventVertex[1] = origin[1];
 594     eventVertex[2] = origin[2];
 595     Float_t eventTime = time;
 596
 597     Int_t nGam;
 598     Float_t impPar,area,pt,rapidity,phi,ww;
 599     Float_t r0=1.3,ra=r0*TMath::Power(fAProjectile,1./3.),rb=r0*TMath::Power(fATarget,1./3.);
 600
 601     TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2);
 602     funcOL->SetParameters(ra,rb);
 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             }
 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]));
 632           }
 633
 634           PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.);
 635         }
 636       }
 637
 638     delete funcOL;
 639 // Header
 640     AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons");
 641 // Event Vertex
 642     header->SetPrimaryVertex(eventVertex);
 643     header->SetInteractionTime(eventTime);
 644     header->SetNProduced(fNpart);
 645     gAlice->SetGenEventHeader(header);
 646
 647 }
 648
 649 void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
 650     AliGenerator::SetPtRange(ptmin, ptmax);
 651 }
 652
 653 void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) {
 654     AliGenerator::SetYRange(ymin, ymax);
 655 }