Thermal Photon fast Generator (Serguei Kiselev)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 27 Nov 2008 10:35:16 +0000 (10:35 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 27 Nov 2008 10:35:16 +0000 (10:35 +0000)
EVGEN/AliGenThermalPhotons.cxx [new file with mode: 0644]
EVGEN/AliGenThermalPhotons.h [new file with mode: 0644]
EVGEN/EVGENLinkDef.h
EVGEN/libEVGEN.pkg

diff --git a/EVGEN/AliGenThermalPhotons.cxx b/EVGEN/AliGenThermalPhotons.cxx
new file mode 100644 (file)
index 0000000..090422a
--- /dev/null
@@ -0,0 +1,646 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+// author: Sergey Kiselev, ITEP, Moscow
+// e-mail: Sergey.Kiselev@cern.ch
+// tel.: 007 495 129 95 45
+//-------------------------------------------------------------------------
+// Generator of direct thermal photons for the reaction A+B, sqrt(S)
+// main assumptions:
+// 1+1 Bjorken scaling hydrodinamics.
+// 1st order phase transition
+// QGP + Mixed (QGP+HHG) + HHG (Hot Hadron Gas) phases, 
+// an ideal massless parton gas and ideal massless HHG 
+// see 
+// F.D.Steffen, nucl-th/9909035
+// F.D.Steffen and M.H.Thoma, Phys.Lett. B510, 98 (2001)
+// T.Peitzmann and M.H.Thoma, Phys.Rep., 364, 175 (2002) 
+//
+// photon rates for QGP: Phys.Rep., 364, 175 (2002), section 2.1.1
+//
+// photon rates for HHG
+// prates for i rho --> pi gamma, pi pi --> rho gamma and rho --> pi pi gamma:
+// Song and Fai, Phys.Rev. C58, 1689 (1998)
+// rates for omega --> pi gamma: Phys.Rev. D44, 2774 (1991)
+//
+// input parameters:
+//       fAProjectile, fATarget - number of nucleons in a nucleus A and B
+//       fMinImpactParam - minimal impct parameter, fm
+//       fMaxImpactParam - maximal impct parameter, fm
+//       fEnergyCMS - sqrt(S) per nucleon pair, AGeV
+//       fTau0 - initial proper time, fm/c
+//       fT0 - initial temperature, GeV
+//       fTc - critical temperature, GeV
+//       fTf - freeze-out temperature, GeV
+//       fGhhg - effective number of degrees of freedom in HHG
+//
+//       fYMin - minimal rapidity of photons 
+//       fYMax - maximal rapidity of photons
+//              in [fYMin,fYMax] uniform distribution of gamma is assumed
+//       fPtMin - minimal p_t value of gamma, GeV/c
+//       fPtMax - maximal p_t value of gamma, GeV/c
+//-------------------------------------------------------------------------
+// comparison with SPS and RHIC data, prediction for LHC can be found in
+// arXiv:0811.2634 [nucl-th]
+//-------------------------------------------------------------------------
+
+//Begin_Html
+/*
+<img src="picts/AliGeneratorClass.gif">
+</pre>
+<br clear=left>
+<font size=+2 color=red>
+<p>The responsible person for this module is
+<a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
+</font>
+<pre>
+*/
+//End_Html
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+#include <TArrayF.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TH1F.h>
+
+#include "AliConst.h"
+#include "AliGenEventHeader.h"
+#include "AliGenThermalPhotons.h"
+#include "AliLog.h"
+#include "AliRun.h"
+
+ClassImp(AliGenThermalPhotons)
+
+// -----------------------------------------------------------------------------------------------------
+static Double_t rateQGP(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - tau (fm), proper time
+// x[1] - yprime, space rapidity
+// par[0] - p_T (GeV), photon transverse momentum 
+// par[1] - y, photon rapidity in the c.m.s. A+A
+// par[2] - tau0 (fm), initial proper time 
+// par[3] - T_0 (GeV), initial temperature
+// par[4] - T_c (GeV), critical temperature
+// par[5] - iProcQGP, process number, 0, 1, 2
+//
+// output:
+// tau EdR/d^3p = tau EdN/d^4xd^3p (fm fm^{-4}GeV^{-2}) 
+//---------------------------------------------------
+  Double_t tau=x[0],yprime=x[1];
+  Double_t pT=par[0],y=par[1],tau0=par[2],T0=par[3],TC=par[4];
+  Int_t iProcQGP=Int_t(par[5]), nFl=3;
+
+  Double_t E=pT*TMath::CosH(yprime-y),T=T0*TMath::Power(tau0/tau,1./3.);
+
+  const Double_t alpha=1./137.;
+// factor to convert from GeV^2 to fm^{-4}GeV^{-2}: (1/hc)**4=(1/0.197)**4=659.921
+  const Double_t factor=659.921;
+  Double_t alphaS=3.*TMath::TwoPi()/((33.-2.*nFl)*TMath::Log(8.*T/TC));
+  const Double_t abc[3]={0.0338, 0.0281, 0.0135} ; // a, b, c for nFf=3
+  Double_t rate=1.;
+
+  switch (iProcQGP) {
+
+    case 0:
+// 1-loop
+      rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-E/T)*T*T*TMath::Log(0.2317*E/alphaS/T);
+    break ;
+
+    case 1:
+// 2-loop: bremsstrahlung
+      rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-E/T)*T*T;
+    break ;
+
+    case 2:
+// 2-loop: annihilation with scattering
+      rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-E/T)*T*E;
+    break ;
+  
+    default:
+      printf("NO iProcQGP=%i \n",iProcQGP);
+  }
+
+  return tau*rate;
+}   
+
+// -----------------------------------------------------------------------------------------------------
+static Double_t fromQGP(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - p_T (GeV), photon p_T
+// par[0] - tau0 (fm), initial proper time 
+// par[1] - T_0 (GeV), initial temperature
+// par[2] - tauCQGP (fm), end of QGP
+// par[3] - yNucl, rapidity of projectile nucleus
+// par[4] - T_c (GeV), critical temperature
+// par[5] - y, photon rapidity
+// par[6] - iProcQGP, process number
+//
+// output:
+// d^{2}N/(dp_t dy) (1/GeV)
+//---------------------------------------------------
+  Double_t pT=x[0];
+  Double_t tau0=par[0],T0=par[1],tauCQGP=par[2],yNucl=par[3],TC=par[4],y=par[5];
+  Int_t iProcQGP=Int_t(par[6]);
+
+  TF2 frateQGP("frateQGP",&rateQGP,tau0,tauCQGP,-yNucl,yNucl,6);
+  frateQGP.SetParameters(pT,y,tau0,T0,TC,iProcQGP);
+  frateQGP.SetParNames("transverse momentum","rapidity","initial time","initial temperature","critical temperature","process number");
+  return TMath::TwoPi()*pT*frateQGP.Integral(tau0,tauCQGP,-yNucl,yNucl,1e-6);
+}   
+         
+// -----------------------------------------------------------------------------------------------------
+static Double_t rateMixQ(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - yprime, space rapidity
+// par[0] - p_T (GeV), photon transverse momentum 
+// par[1] - y, photon rapidity in the c.m.s. A+A
+// par[2] - T_c (GeV), critical temperature
+// par[3] - iProcQGP, process number, 0, 1, 2
+//
+// output:
+// EdR/d^3p = EdN/d^4xd^3p (fm fm^{-4}GeV^{-2}) 
+//---------------------------------------------------
+  Double_t yprime=x[0];
+  Double_t pT=par[0],y=par[1],TC=par[2];
+  Int_t iProcQGP=Int_t(par[3]),nFl=3;
+
+  Double_t E=pT*TMath::CosH(yprime-y),T=TC;
+
+  const Double_t alpha=1./137.;
+// factor to convert from GeV^2 to fm^{-4}GeV^{-2}: (1/hc)**4=(1/0.197)**4=659.921
+  const Double_t factor=659.921;
+  Double_t alphaS=3.*TMath::TwoPi()/((33.-2.*nFl)*TMath::Log(8.*T/TC));
+  const Double_t abc[3]={0.0338, 0.0281, 0.0135}; // a, b, c for nF=3
+  Double_t rate=1.;
+
+  switch (iProcQGP) {
+
+    case 0:
+// 1-loop
+      rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-E/T)*T*T*TMath::Log(0.2317*E/alphaS/T);
+    break ;
+
+    case 1:
+// 2-loop: bremsstrahlung
+      rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-E/T)*T*T;
+    break ;
+
+    case 2:
+// 2-loop: annihilation with scattering
+      rate=factor*abc[iProcQGP]*alpha*alphaS*TMath::Exp(-E/T)*T*E;
+    break ;
+  
+    default:
+      printf("NO iProcQGP=%i \n",iProcQGP);
+  }
+
+  return rate;
+}   
+
+// -----------------------------------------------------------------------------------------------------
+static Double_t fromMixQ(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - p_T (GeV), photon p_T
+// par[0] - lamQGP
+// par[1] - yNucl, rapidity of projectile nucleus
+// par[2] - T_c (GeV), critical temperature
+// par[3] - y, photon rapidity
+// par[4] - iProcQGP, process number
+//
+// output:
+// d^{2}N/(dp_t dy) (1/GeV)
+//---------------------------------------------------
+  Double_t pT=x[0];
+  Double_t lamQGP=par[0],yNucl=par[1],TC=par[2],y=par[3];
+  Int_t iProcQGP=Int_t(par[4]);
+
+  TF1 frateMixQ("frateMixQ",&rateMixQ,-yNucl,yNucl,4);
+  frateMixQ.SetParameters(pT,y,TC,iProcQGP);
+  frateMixQ.SetParNames("transverse momentum","rapidity","critical temperature","process number");
+  return TMath::TwoPi()*pT*lamQGP*frateMixQ.Integral(-yNucl,yNucl);
+}   
+         
+// -----------------------------------------------------------------------------------------------------
+static Double_t rateMixH(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - yprime, space rapidity
+// par[0] - p_T (GeV), photon transverse momentum 
+// par[1] - y, photon rapidity in the c.m.s. A+A
+// par[2] - T_c (GeV), critical temperature
+// par[3] - iProcHHG, process number
+//
+// output:
+// EdR/d^3p = EdN/d^4xd^3p (fm^{-4}GeV^{-2}) 
+//---------------------------------------------------
+  Double_t yprime=x[0];
+  Double_t pT=par[0],y=par[1],TC=par[2];
+  Int_t iProcHHG=Int_t(par[3]);
+
+  Double_t E=pT*TMath::CosH(yprime-y),T=TC;
+  const Double_t mPi=0.135;
+  Double_t xx=T/mPi,yy=E/mPi;
+  Double_t F,rate=1.,Emin,factor;
+  const Double_t mOm=0.783,width=0.00844,Br=0.085,E0=0.379,pi=TMath::Pi();
+  const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm
+
+  switch (iProcHHG) {
+
+    case 0:
+// pi rho --> pi gamma
+      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);
+      rate=T*T*TMath::Exp(-E/T+F);
+    break ;
+
+    case 1:
+// pi pi --> rho gamma
+      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);
+      rate=T*T*TMath::Exp(-E/T+F);
+    break ;
+
+    case 2:
+// rho --> pi pi gamma
+      F=-6.295+1.6459*xx+(-0.4015+0.089*xx)*yy+(-0.954+2.05777*xx)/yy;
+      rate=T*T*TMath::Exp(-E/T+F);
+    break ;
+  
+    case 3:
+// omega --> pi gamma
+      Emin=mOm*(E*E+E0*E0)/(2.*E*E0);
+      factor=(3.*mOm*width*Br)/(16.*pi*pi*pi*E0);
+      rate=factor*T*(Emin+T)*TMath::Exp(-Emin/T)/E/hc4;
+    break ;
+  
+    default:
+      printf("NO iProcHHG=%i \n",iProcHHG);
+  }
+
+  return rate;
+}   
+
+// -----------------------------------------------------------------------------------------------------
+static Double_t fromMixH(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - p_T (GeV), photon p_T
+// par[0] - lamHHG
+// par[1] - yNucl, rapidity of projectile nucleus
+// par[2] - T_c (GeV), critical temperature
+// par[3] - y, photon rapidity
+// par[4] - iProcHHG, process number
+//
+// output:
+// d^{2}N/(dp_t dy) (1/GeV)
+//---------------------------------------------------
+  Double_t pT=x[0];
+  Double_t lamHHG=par[0],yNucl=par[1],TC=par[2],y=par[3];
+  Int_t iProcHHG=Int_t(par[4]);
+
+  TF1 frateMixH("frateMixH",&rateMixH,-yNucl,yNucl,4);
+  frateMixH.SetParameters(pT,y,TC,iProcHHG);
+  frateMixH.SetParNames("transverse momentum","rapidity","critical temperature","process number");
+  return TMath::TwoPi()*pT*lamHHG*frateMixH.Integral(-yNucl,yNucl);
+}   
+         
+// -----------------------------------------------------------------------------------------------------
+static Double_t rateHHG(Double_t *x, Double_t *par) {
+//---------------------------------------------------
+// input:
+// x[0] - tau (fm), proper time
+// x[1] - yprime, space rapidity
+// par[0] - p_T (GeV), photon transverse momentum 
+// par[1] - y, photon rapidity in the c.m.s. A+A
+// par[2] - tauCHHG (fm), start of HHG 
+// par[3] - T_c (GeV), critical temperature
+// par[4] - iProcHHG, process number
+//
+// output:
+// EdR/d^3p = EdN/d^4xd^3p (fm^{-4}GeV^{-2}) 
+//---------------------------------------------------
+  Double_t tau=x[0],yprime=x[1];
+  Double_t pT=par[0],y=par[1],tauCHHG=par[2],TC=par[3];
+  Int_t iProcHHG=Int_t(par[4]);
+
+  Double_t E=pT*TMath::CosH(yprime-y),T=TC*TMath::Power(tauCHHG/tau,1./3.);
+  const Double_t mPi=0.135;
+  Double_t xx=T/mPi,yy=E/mPi;
+  Double_t F,rate=1.,Emin,factor;
+  const Double_t mOm=0.783,width=0.00844,Br=0.085,E0=0.379,pi=TMath::Pi();
+  const Double_t hc=0.197,hc4=hc*hc*hc*hc; // GeV*fm
+
+  switch (iProcHHG) {
+
+    case 0:
+// pi rho --> pi gamma
+      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);
+      rate=T*T*TMath::Exp(-E/T+F);
+    break ;
+
+    case 1:
+// pi pi --> rho gamma
+      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);
+      rate=T*T*TMath::Exp(-E/T+F);
+    break ;
+
+    case 2:
+// rho --> pi pi gamma
+      F=-6.295+1.6459*xx+(-0.4015+0.089*xx)*yy+(-0.954+2.05777*xx)/yy;
+      rate=T*T*TMath::Exp(-E/T+F);
+    break ;
+  
+    case 3:
+// omega --> pi gamma
+      Emin=mOm*(E*E+E0*E0)/(2.*E*E0);
+      factor=(3.*mOm*width*Br)/(16.*pi*pi*pi*E0);
+      rate=factor*T*(Emin+T)*TMath::Exp(-Emin/T)/E/hc4;
+    break ;
+
+    default:
+      printf("NO iProcHHG=%i \n",iProcHHG);
+  }
+  return tau*rate;
+}   
+
+// -----------------------------------------------------------------------------------------------------
+static Double_t fromHHG(Double_t *x, Double_t *par) {
+// Thermal photon spectrum from Hot Hadron Gas (HHG)
+//  F.D.Steffen, nucl-th/9909035
+//  T.Peitzmann and M.H.Thoma, Phys.Rep., 364, 175 (2002), section 2.2.2 
+//---------------------------------------------------
+// input:
+// x[0] - p_T (GeV), photon p_T
+// par[0] - tauCHHG (fm), start of HHG 
+// par[1] - tauF (fm), freeze-out proper time 
+// par[2] - yNucl, rapidity of projectile nucleus
+// par[3] - T_c (GeV), critical temperature
+// par[4] - y, photon rapidity
+// par[5] - iProcHHG, process number
+//
+// output:
+// d^{2}N/(dp_t dy) (1/GeV)
+//---------------------------------------------------
+  Double_t pT=x[0];
+  Double_t tauCHHG=par[0],tauF=par[1],yNucl=par[2],TC=par[3],y=par[4],iProcHHG=par[5];
+
+  TF2 frateHHG("frateHHG",&rateHHG,tauCHHG,tauF,-yNucl,yNucl,5);
+  frateHHG.SetParameters(pT,y,tauCHHG,TC,iProcHHG);
+  frateHHG.SetParNames("transverse momentum","rapidity","start of HHG","criti temperature","process number");
+  return TMath::TwoPi()*pT*frateHHG.Integral(tauCHHG,tauF,-yNucl,yNucl,1e-6);
+}   
+
+// -----------------------------------------------------------------------------------------------------
+static Double_t fOverlapAB(Double_t *x, Double_t *par)
+{
+//-------------------------------------------------------------------------
+// overlap area at the impact parameter b
+// input:
+// x[0] - impact parameter b < RA+RB
+// par[0] - radius of A
+// par[1] - radius of B
+//-------------------------------------------------------------------------
+
+  Double_t b=x[0], RA=par[0], RB=par[1];
+  if(RB>RA) {RA=par[1]; RB=par[0];} // RA > RB
+
+  if(b>(RA+RB)) {
+    return 0.;
+  }
+
+  if(b<=(RA-RB)) {
+    return TMath::Pi()*RB*RB;
+  }
+
+  Double_t p=0.5*(b+RA+RB), S=TMath::Sqrt(p*(p-b)*(p-RA)*(p-RB)), h=2.*S/b;
+  Double_t sA=RA*RA*TMath::ASin(h/RA)-h*TMath::Sqrt(RA*RA-h*h);
+  Double_t sB=RB*RB*TMath::ASin(h/RB)-h*TMath::Sqrt(RB*RB-h*h);
+  if(RA>RB && b*b<RA*RA-RB*RB) sB=TMath::Pi()*RB*RB-sB;
+
+  return sA+sB;
+
+}
+
+//_____________________________________________________________________________
+AliGenThermalPhotons::AliGenThermalPhotons()
+    :AliGenerator(-1),
+        fAProjectile(0.),
+        fATarget(0.),
+        fEnergyCMS(0.),
+        fMinImpactParam(0.),
+        fMaxImpactParam(0.),
+        fTau0(0.),
+        fT0(0.),
+        fTc(0.),
+        fTf(0.),
+        fGhhg(0),
+        fSumPt()
+{
+    //
+    // Default constructor
+    //
+    SetCutVertexZ();
+    SetPtRange();
+    SetYRange();
+}
+
+//_____________________________________________________________________________
+AliGenThermalPhotons::AliGenThermalPhotons(Int_t npart)
+    :AliGenerator(npart),
+        fAProjectile(208),
+        fATarget(208),
+        fEnergyCMS(5500.),
+        fMinImpactParam(0.),
+        fMaxImpactParam(0.),
+        fTau0(0.1),
+        fT0(0.650),
+        fTc(0.170),
+        fTf(0.100),
+        fGhhg(8),
+        fSumPt()
+{
+  // 
+  // Standard constructor
+  //
+
+    fName="ThermalPhotons";
+    fTitle="Direct thermal photons in 1+1 Bjorken hydrodynamics";
+
+    SetCutVertexZ();
+    SetPtRange();
+    SetYRange();
+}
+
+//_____________________________________________________________________________
+AliGenThermalPhotons::~AliGenThermalPhotons()
+{
+  //
+  // Standard destructor
+  //
+    delete fSumPt;
+}
+
+//_____________________________________________________________________________
+void AliGenThermalPhotons::Init()
+{
+
+  const Double_t step=0.1; 
+  Int_t nPt=Int_t((fPtMax-fPtMin)/step);
+
+  fSumPt = new TH1F("fSumPt","thermal #gamma from QGP",nPt,fPtMin,fPtMax);
+
+  Double_t yRap=0.;
+  const Int_t nCo=3,nFl=3; //  number of colors for QGP
+  Double_t gQGP=2.*(nCo*nCo-1.)+(7./8.)*4.*nCo*nFl; //  number of degrees of freedom in QGP
+  Double_t yNucl=TMath::ACosH(fEnergyCMS/2.);
+  Double_t tauCQGP=TMath::Power(fT0/fTc,3.)*fTau0,tauCHHG=gQGP*tauCQGP/fGhhg,tauF=tauCHHG*TMath::Power(fTc/fTf,3.);
+  Double_t lambda1=tauCQGP*(gQGP/(gQGP-fGhhg)),lambda2=-fGhhg/(gQGP-fGhhg);
+  Double_t lamQGP=(tauCHHG-tauCQGP)*(lambda1+0.5*lambda2*(tauCHHG+tauCQGP)),lamHHG=0.5*(tauCHHG-tauCQGP)*(tauCHHG+tauCQGP)-lamQGP;
+
+  Double_t pt,w;
+// photons from pure QGP phase
+  for(Int_t j=0; j<3; j++) {
+    TF1 func("func",&fromQGP,fPtMin,fPtMax,7);
+    func.SetParameters(fTau0,fT0,tauCQGP,yNucl,fTc,yRap,j);
+    func.SetParNames("nuclear radius","initial time","initial temperature","end of pure QGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
+    for(Int_t i=0; i<nPt; i++) {
+      pt=fPtMin+(i+0.5)*step;
+      w=func.Eval(pt);
+      fSumPt->AddBinContent(i+1,w);
+    }
+  }
+
+// photons from mixed QGP phase
+  for(Int_t j=0; j<3; j++) {
+    TF1 func("func",&fromMixQ,fPtMin,fPtMax,5);
+    func.SetParameters(lamQGP,yNucl,fTc,yRap,j);
+    func.SetParNames("lamQGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
+    for(Int_t i=0; i<nPt; i++) {
+      pt=fPtMin+(i+0.5)*step;
+      w=func.Eval(pt);
+      fSumPt->AddBinContent(i+1,w);
+    }
+  }
+
+// photons from mixed HHG phase
+  for(Int_t j=0; j<4; j++) {
+    TF1 func("func",&fromMixH,fPtMin,fPtMax,5);
+    func.SetParameters(lamHHG,yNucl,fTc,yRap,j);
+    func.SetParNames("lamQGP","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
+    for(Int_t i=0; i<nPt; i++) {
+      pt=fPtMin+(i+0.5)*step;
+      w=func.Eval(pt);
+      fSumPt->AddBinContent(i+1,w);
+    }
+  }
+  
+// photons from pure HHG phase
+  for(Int_t j=0; j<4; j++) {
+    TF1 func("func",&fromHHG,fPtMin,fPtMax,6);
+    func.SetParameters(tauCHHG,tauF,yNucl,fTc,yRap,j);
+    func.SetParNames("nuclear radius","start HHG","freeze-out proper time","rapidity of projectile nucleus","critical temperature","photon rapidity","process number");
+    for(Int_t i=0; i<nPt; i++) {
+      pt=fPtMin+(i+0.5)*step;
+      w=func.Eval(pt);
+      fSumPt->AddBinContent(i+1,w);
+    }
+  }
+  
+}
+
+//_____________________________________________________________________________
+void AliGenThermalPhotons::Generate()
+{
+  //
+  // Generate thermal photons of a event 
+  //
+
+    Float_t polar[3]= {0,0,0};
+    Float_t origin[3];
+    Float_t p[3];
+    Float_t random[6];
+    Int_t nt;
+
+    for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
+/*
+    if(fVertexSmear==kPerEvent) {
+      Vertex();
+      for (j=0;j<3;j++) origin[j]=fVertex[j];
+    }
+*/
+    TArrayF eventVertex;
+    eventVertex.Set(3);
+    eventVertex[0] = origin[0];
+    eventVertex[1] = origin[1];
+    eventVertex[2] = origin[2];
+
+    Int_t nGam;
+    Float_t impPar,area,pt,rapidity,phi,ww;
+    Float_t r0=1.3,RA=r0*TMath::Power(fAProjectile,1./3.),RB=r0*TMath::Power(fATarget,1./3.);
+
+    TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2); 
+    funcOL->SetParameters(RA,RB);
+    funcOL->SetParNames("radiusA","radiusB");
+
+    impPar=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
+    area=funcOL->Eval(impPar);
+
+      ww=area*(fYMax-fYMin)*fSumPt->GetBinWidth(1)*fSumPt->GetSumOfWeights();
+      nGam=Int_t(ww);
+      if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
+
+      if(nGam) {
+        for(Int_t i=0; i<nGam; i++) {
+          pt=fSumPt->GetRandom();
+          Rndm(random,2);
+          rapidity=(fYMax-fYMin)*random[0]+fYMin;
+          phi=2.*TMath::Pi()*random[1];
+          p[0]=pt*TMath::Cos(phi);
+          p[1]=pt*TMath::Sin(phi);
+          p[2]=pt*TMath::SinH(rapidity);
+
+         if(fVertexSmear==kPerTrack) {
+            Rndm(random,6);
+           for (Int_t j=0;j<3;j++) {
+             origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
+             TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
+           }
+         }
+
+         PushTrack(fTrackIt,-1,22,p,origin,polar,0,kPPrimary,nt,1.);
+        }
+      }
+
+    delete funcOL;
+// Header
+    AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons");
+// Event Vertex
+    header->SetPrimaryVertex(eventVertex);
+    header->SetNProduced(fNpart);
+    gAlice->SetGenEventHeader(header);
+
+}
+
+void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
+    AliGenerator::SetPtRange(ptmin, ptmax);
+}
+
+void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) {
+    AliGenerator::SetYRange(ymin, ymax);
+}
diff --git a/EVGEN/AliGenThermalPhotons.h b/EVGEN/AliGenThermalPhotons.h
new file mode 100644 (file)
index 0000000..8013cab
--- /dev/null
@@ -0,0 +1,100 @@
+#ifndef ALIGENTHERMALPHOTONS_H
+#define ALIGENTHERMALPHOTONS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+// author: Sergey Kiselev, ITEP, Moscow
+// e-mail: Sergey.Kiselev@cern.ch
+// tel.: 007 495 129 95 45
+//-------------------------------------------------------------------------
+// Generator of direct thermal photons for the reaction A+B, sqrt(S)
+// main assumptions:
+// 1+1 Bjorken scaling hydrodinamics.
+// 1st order phase transition
+// QGP + Mixed (QGP+HHG) + HHG (Hot Hadron Gas) phases, 
+// an ideal massless parton gas and ideal massless HHG 
+// see 
+// F.D.Steffen, nucl-th/9909035
+// F.D.Steffen and M.H.Thoma, Phys.Lett. B510, 98 (2001)
+// T.Peitzmann and M.H.Thoma, Phys.Rep., 364, 175 (2002) 
+//
+// photon rates for QGP: Phys.Rep., 364, 175 (2002), section 2.1.1
+//
+// photon rates for HHG
+// prates for i rho --> pi gamma, pi pi --> rho gamma and rho --> pi pi gamma:
+// Song and Fai, Phys.Rev. C58, 1689 (1998)
+// rates for omega --> pi gamma: Phys.Rev. D44, 2774 (1991)
+//
+// input parameters:
+//       fAProjectile, fATarget - number of nucleons in a nucleus A and B
+//       fMinImpactParam - minimal impct parameter, fm
+//       fMaxImpactParam - maximal impct parameter, fm
+//       fEnergyCMS - sqrt(S) per nucleon pair, AGeV
+//       fTau0 - initial proper time, fm/c
+//       fT0 - initial temperature, GeV
+//       fTc - critical temperature, GeV
+//       fTf - freeze-out temperature, GeV
+//       fGhhg - effective number of degrees of freedom in HHG
+//
+//       fYMin - minimal rapidity of photons 
+//       fYMax - maximal rapidity of photons
+//              in [fYMin,fYMax] uniform distribution of gamma is assumed
+//       fPtMin - minimal p_t value of gamma, GeV/c
+//       fPtMax - maximal p_t value of gamma, GeV/c
+//-------------------------------------------------------------------------
+// comparison with SPS and RHIC data, prediction for LHC can be found in
+// arXiv:0811.2634 [nucl-th]
+//-------------------------------------------------------------------------
+
+class TH1F;
+
+#include "AliGenerator.h"
+
+class AliGenThermalPhotons : public AliGenerator
+{
+ public:
+
+  AliGenThermalPhotons();
+  AliGenThermalPhotons(Int_t npart);
+  virtual ~AliGenThermalPhotons();
+  virtual void Generate();
+  virtual void Init();
+  virtual void SetPtRange(Float_t ptmin = 0.1, Float_t ptmax=10.);
+  virtual void SetYRange(Float_t ymin = -1., Float_t ymax=1.);
+
+// Setters
+    virtual void SetAProjectile(Float_t a = 208) {fAProjectile = a;}
+    virtual void SetATarget(Float_t a = 208)     {fATarget     = a;}
+    virtual void SetEnergyCMS(Float_t energy = 5500.) {fEnergyCMS = energy;}
+    virtual void    SetImpactParameterRange(Float_t bmin = 0., Float_t bmax = 0.)
+       {fMinImpactParam=bmin; fMaxImpactParam=bmax;}
+    virtual void    SetTau0(Float_t tau0 = 0.1)             {fTau0   = tau0;}
+    virtual void    SetT0(Float_t   T0   = 0.650)           {fT0     = T0;}
+    virtual void    SetTc(Float_t   Tc   = 0.170)           {fTc     = Tc;}
+    virtual void    SetTf(Float_t   Tf   = 0.100)           {fTf     = Tf;}
+    virtual void    SetGhhg(Int_t   Ghhg = 8)               {fGhhg   = Ghhg;}
+
+ protected:
+  Float_t fAProjectile;     // Projectile nucleus mass number
+  Float_t fATarget;         // Target nucleus mass number
+  Float_t fEnergyCMS;       // Center of mass energy
+  Float_t fMinImpactParam;  // minimum impact parameter
+  Float_t fMaxImpactParam;  // maximum impact parameter        
+  Float_t fTau0;            // initial proper time, fm 
+  Float_t fT0;              // initial temperature, GeV        
+  Float_t fTc;              // critical temperature, GeV       
+  Float_t fTf;              // freeze-out temperature, GeV     
+  Int_t   fGhhg;            // number of degrees of freedom in HHG     
+
+  TH1F *fSumPt;             // histo with pt from all origins
+
+ private:
+
+  AliGenThermalPhotons(const AliGenThermalPhotons & ThermalPhotons);
+  AliGenThermalPhotons& operator = (const AliGenThermalPhotons & ThermalPhotons) ;
+
+
+  ClassDef(AliGenThermalPhotons, 1) // Direct thermal photon generator
+};
+#endif
index 07218a7..3e63b31 100644 (file)
@@ -57,4 +57,5 @@
 #pragma link C++ class  AliGenCorrHF+;
 #pragma link C++ class  AliGenCosmicsParam+;
 #pragma link C++ class  AliGenKrypton+;
+#pragma link C++ class  AliGenThermalPhotons+;
 #endif
index 1893620..d1bd41b 100644 (file)
@@ -20,7 +20,7 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliSlowNucleonModel.cxx AliSlowNucleonModelExp.cxx \
                AliGenMUONCocktail.cxx AliGenMUONCocktailpp.cxx AliGenHBTosl.cxx \
                AliGenReaderEMD.cxx AliDecayerPolarized.cxx AliGenCorrHF.cxx AliGenCosmicsParam.cxx \
-               AliGenKrypton.cxx
+               AliGenKrypton.cxx AliGenThermalPhotons.cxx
 
 # Headerfiles for this particular package (Path respect to own directory)
 HDRS= $(SRCS:.cxx=.h)