Fix for event mixing, when it was selecting events out of range of multiplicity cut
[u/mrichter/AliRoot.git] / 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(Double_t *x, 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(Double_t *x, 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(Double_t *x, 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(Double_t *x, 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(Double_t *x, 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(Double_t *x, 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(Double_t *x, 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(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 //---------------------------------------------------
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(Double_t *x, 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
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 p[3];
577     Float_t random[6];
578     Int_t nt;
579
580     for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
581 /*
582     if(fVertexSmear==kPerEvent) {
583       Vertex();
584       for (j=0;j<3;j++) origin[j]=fVertex[j];
585     }
586 */
587     TArrayF eventVertex;
588     eventVertex.Set(3);
589     eventVertex[0] = origin[0];
590     eventVertex[1] = origin[1];
591     eventVertex[2] = origin[2];
592
593     Int_t nGam;
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.);
596
597     TF1 *funcOL = new TF1("funcOL",fOverlapAB,fMinImpactParam,fMaxImpactParam,2); 
598     funcOL->SetParameters(ra,rb);
599     funcOL->SetParNames("radiusA","radiusB");
600
601     impPar=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
602     area=funcOL->Eval(impPar);
603
604       ww=area*(fYMax-fYMin)*fSumPt->GetBinWidth(1)*fSumPt->GetSumOfWeights();
605       nGam=Int_t(ww);
606       if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
607
608       if(nGam) {
609         for(Int_t i=0; i<nGam; i++) {
610           pt=fSumPt->GetRandom();
611           Rndm(random,2);
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);
617
618           if(fVertexSmear==kPerTrack) {
619             Rndm(random,6);
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]));
623             }
624           }
625
626           PushTrack(fTrackIt,-1,22,p,origin,polar,0,kPPrimary,nt,1.);
627         }
628       }
629
630     delete funcOL;
631 // Header
632     AliGenEventHeader* header = new AliGenEventHeader("ThermalPhotons");
633 // Event Vertex
634     header->SetPrimaryVertex(eventVertex);
635     header->SetNProduced(fNpart);
636     gAlice->SetGenEventHeader(header);
637
638 }
639
640 void AliGenThermalPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
641     AliGenerator::SetPtRange(ptmin, ptmax);
642 }
643
644 void AliGenThermalPhotons::SetYRange(Float_t ymin, Float_t ymax) {
645     AliGenerator::SetYRange(ymin, ymax);
646 }