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 prompt photons for the reaction A+B, sqrt(S)
24 // 1. flat rapidity distribution
25 // 2. all existing p+p(pbar) data at y_{c.m.} can be described by the function
26 // F(x_T) = (sqrt(s))^5 Ed^3sigma/d^3p, x_T = 2p_t/sqrt(s)
27 // all data points cover the region x_T: 0.01 - 0.6
28 // see Nucl.Phys.A783:577-582,2007, hep-ex/0609037
29 // 3. binary scaling: for A+B at the impact parameter b
30 // Ed^3N^{AB}(b)/d^3p = Ed^3sigma^{pp}/d^3p A B T_{AB}(b),
31 // T_{AB}(b) - nuclear overlapping fuction, calculated in the Glauber approach,
32 // nuclear density is parametrized by a Woods-Saxon with nuclear radius
33 // R_A = 1.19 A^{1/3} - 1.61 A^{-1/3} fm and surface thickness a=0.54 fm
34 // 4. nuclear effects (Cronin, shadowing, ...) are ignored
37 // fAProjectile, fATarget - number of nucleons in a nucleus A and B
38 // fMinImpactParam - minimal impct parameter, fm
39 // fMaxImpactParam - maximal impct parameter, fm
40 // fEnergyCMS - sqrt(S) per nucleon pair, AGeV
42 // fYMin - minimal rapidity of photons
43 // fYMax - maximal rapidity of photons
44 // fPtMin - minimal p_t value of gamma, GeV/c
45 // fPtMax - maximal p_t value of gamma, GeV/c
46 //-------------------------------------------------------------------------
47 // comparison with SPS and RHIC data, prediction for LHC can be found in
48 // arXiv:0811.2634 [nucl-th]
49 //-------------------------------------------------------------------------
53 <img src="picts/AliGeneratorClass.gif">
56 <font size=+2 color=red>
57 <p>The responsible person for this module is
58 <a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
64 ///////////////////////////////////////////////////////////////////
70 #include "AliGenEventHeader.h"
71 #include "AliGenPromptPhotons.h"
75 ClassImp(AliGenPromptPhotons)
77 TF1* AliGenPromptPhotons::fgDataPt = NULL;
78 TF1* AliGenPromptPhotons::fgWSzA = NULL;
79 TF1* AliGenPromptPhotons::fgWSzB = NULL;
80 TF1* AliGenPromptPhotons::fgTA = NULL;
81 TF1* AliGenPromptPhotons::fgTB = NULL;
82 TF1* AliGenPromptPhotons::fgTAxTB = NULL;
83 TF1* AliGenPromptPhotons::fgTAB = NULL;
85 //_____________________________________________________________________________
86 AliGenPromptPhotons::AliGenPromptPhotons()
95 // Default constructor
102 //_____________________________________________________________________________
103 AliGenPromptPhotons::AliGenPromptPhotons(Int_t npart)
104 :AliGenerator(npart),
112 // Standard constructor
115 fName="PromptPhotons";
116 fTitle="Prompt photons from pp data fit + T_AB";
123 //_____________________________________________________________________________
124 AliGenPromptPhotons::~AliGenPromptPhotons()
127 // Standard destructor
138 //_____________________________________________________________________________
139 void AliGenPromptPhotons::Init()
142 fgDataPt = new TF1("fgDataPt",FitData,fPtMin,fPtMax,1);
143 fgDataPt->SetParameter(0,fEnergyCMS);
145 const Double_t d=0.54; // surface thickness (fm)
146 const Double_t ra = 1.19*TMath::Power(fAProjectile,1./3.)-1.61/TMath::Power(fAProjectile,1./3.);
147 const Double_t eps=0.01; // cut WS at ramax: WS(ramax)/WS(0)=eps
148 const Double_t ramax=ra+d*TMath::Log((1.-eps+TMath::Exp(-ra/d))/eps);
150 TF1 fWSforNormA("fWSforNormA",&WSforNorm,0,ramax,2);
151 fWSforNormA.SetParameters(ra,d);
152 fWSforNormA.SetParNames("RA","d");
153 const Double_t ca=1./fWSforNormA.Integral(0.,ramax);
155 const Double_t rb=1.19*TMath::Power(fATarget,1./3.)-1.61/TMath::Power(fATarget,1./3.);
156 const Double_t rbmax=rb+d*TMath::Log((1.-eps+TMath::Exp(-rb/d))/eps);
158 TF1 fWSforNormB("fWSforNormB",&WSforNorm,0,rbmax,2);
159 fWSforNormB.SetParameters(rb,d);
160 fWSforNormB.SetParNames("RB","d");
161 const Double_t cb=1./fWSforNormB.Integral(0.,rbmax);
163 fgWSzA = new TF1("fgWSzA",WSz,0.,ramax,4);
164 fgWSzA->SetParameter(0,ra);
165 fgWSzA->SetParameter(1,d);
166 fgWSzA->SetParameter(2,ca);
168 fgTA = new TF1("fgTA",TA,0.,ramax,1);
169 fgTA->SetParameter(0,ramax);
171 fgWSzB = new TF1("fgWSzB",WSz,0.,rbmax,4);
172 fgWSzB->SetParameter(0,rb);
173 fgWSzB->SetParameter(1,d);
174 fgWSzB->SetParameter(2,cb);
176 fgTB = new TF1("fgTB",TB,0.,TMath::Pi(),3);
177 fgTB->SetParameter(0,rbmax);
179 fgTAxTB = new TF1("fgTAxTB",TAxTB,0,ramax,2);
180 fgTAxTB->SetParameter(0,rbmax);
182 fgTAB = new TF1("fgTAB",TAB,0.,ramax+rbmax,2);
183 fgTAB->SetParameter(0,ramax);
184 fgTAB->SetParameter(1,rbmax);
188 //_____________________________________________________________________________
189 void AliGenPromptPhotons::Generate()
192 // Generate thermal photons of a event
195 Float_t polar[3]= {0,0,0};
202 for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
205 if(fVertexSmear==kPerEvent) {
207 for (j=0;j<3;j++) origin[j]=fVertex[j];
213 eventVertex[0] = origin[0];
214 eventVertex[1] = origin[1];
215 eventVertex[2] = origin[2];
216 Float_t eventTime = time;
219 Float_t b,pt,rapidity,phi,ww;
221 b=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
223 Double_t dsdyPP=fgDataPt->Integral(fPtMin,fPtMax); // pb, fm^2 = 1e10 pb
224 ww=dsdyPP*1.0e-10*(fYMax-fYMin)*fAProjectile*fATarget*fgTAB->Eval(b);
226 if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
229 for(Int_t i=0; i<nGam; i++) {
230 pt=fgDataPt->GetRandom();
232 rapidity=(fYMax-fYMin)*random[0]+fYMin;
233 phi=2.*TMath::Pi()*random[1];
234 p[0]=pt*TMath::Cos(phi);
235 p[1]=pt*TMath::Sin(phi);
236 p[2]=pt*TMath::SinH(rapidity);
238 if(fVertexSmear==kPerTrack) {
240 for (Int_t j=0;j<3;j++) {
241 origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
242 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
245 time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
246 TMath::Cos(2*random[0]*TMath::Pi())*
247 TMath::Sqrt(-2*TMath::Log(random[1]));
250 PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.);
255 AliGenEventHeader* header = new AliGenEventHeader("PromptPhotons");
257 header->SetPrimaryVertex(eventVertex);
258 header->SetInteractionTime(eventTime);
259 header->SetNProduced(fNpart);
260 gAlice->SetGenEventHeader(header);
264 void AliGenPromptPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
265 AliGenerator::SetPtRange(ptmin, ptmax);
268 void AliGenPromptPhotons::SetYRange(Float_t ymin, Float_t ymax) {
269 AliGenerator::SetYRange(ymin, ymax);
272 //**********************************************************************************
273 Double_t AliGenPromptPhotons::FitData(const Double_t* x, const Double_t* par) {
274 //---------------------------------------------------
277 // par[0]=sqrt(s_NN) (GeV),
280 // d^{2}#sigma/(dp_t dy) (pb/GeV)
281 //---------------------------------------------------
283 // d^{2}#sigma/(dp_t dy) = (2 pi p_t) Ed^{3}#sigma/d^{3}p
285 // data presentation: Nucl.Phys.A783:577-582,2007, hep-ex/0609037, fig.3
286 // F(x_t)=(#sqrt{s})^{5} Ed^{3}#sigma/d^{3}p
287 //---------------------------------------------------
288 // approximate tabulation of F(x_t)
290 const Double_t log10x[nMax]={ -2., -1., -0.6, -0.3};
291 const Double_t log10F[nMax]={ 19., 13., 9.8, 7.};
293 const Double_t xT=2.*x[0]/par[0];
294 const Double_t log10xT=TMath::Log10(xT);
296 while(log10xT>log10x[i] && i<nMax) i++;
298 if(i==nMax) i=nMax-1;
299 const Double_t alpha=(log10F[i]-log10F[i-1])/(log10x[i]-log10x[i-1]);
300 const Double_t beta=log10F[i-1]-alpha*log10x[i-1];
302 return (TMath::TwoPi()*x[0])*TMath::Power(10.,alpha*log10xT+beta)/TMath::Power(par[0],5.);
305 //**********************************************************************************
306 Double_t AliGenPromptPhotons::WSforNorm(const Double_t* x, const Double_t* par) {
307 //---------------------------------------------------
310 // par[0] - R (fm), radius
311 // par[1] - d (fm), surface thickness
314 // 4 pi r**2 /(1+exp((r-R)/d))
316 // Wood Saxon (WS) C/(1+exp((r-RA)/d)) (nuclons/fm^3)
317 // To get the normalization A = (Integral 4 pi r**2 dr WS):
318 // C = A / (Integral 4 pi r**2 dr 1/(1+exp((r-RA)/d)) )
319 // Thus me need 4 pi r**2 /(1+exp((r-RA)/d)) (1/fm)
320 //---------------------------------------------------
321 const Double_t r=x[0];
322 const Double_t bigR=par[0],d=par[1];
324 return 4.*TMath::Pi()*r*r/(1.+TMath::Exp((r-bigR)/d));
327 //**********************************************************************************
328 Double_t AliGenPromptPhotons::WSz(const Double_t* x, const Double_t* par) {
329 //---------------------------------------------------
332 // par[0] - R (fm), radius
333 // par[1] - d (fm), surface thickness
334 // par[2] - C (nucleons/fm**2), normalization factor
335 // par[3] - b (fm), impact parameter
338 // Wood Saxon Parameterisation
339 // as a function of z for fixed b (1/fm^3)
340 //---------------------------------------------------
341 const Double_t z=x[0];
342 const Double_t bigR=par[0],d=par[1],C=par[2],b=par[3];
344 return C/(1.+TMath::Exp((TMath::Sqrt(b*b+z*z)-bigR)/d));
347 //**********************************************************************************
348 Double_t AliGenPromptPhotons::TA(const Double_t* x, const Double_t* par) {
349 //---------------------------------------------------
351 // x[0] - b (fm), impact parameter
352 // par[0] - RAMAX (fm), max. value of projectile radius
355 // nuclear thickness function T_A(b) (1/fm^2)
356 //---------------------------------------------------
357 const Double_t b=x[0];
358 const Double_t ramax=par[0];
360 fgWSzA->SetParameter(3,b);
362 return 2.*fgWSzA->Integral(0.,TMath::Sqrt(ramax*ramax-b*b));
365 //**********************************************************************************
366 Double_t AliGenPromptPhotons::TB(const Double_t* x, const Double_t* par) {
367 //---------------------------------------------------
370 // par[0] - RBMAX (fm), max. value of target radius
371 // par[1] - b (fm), impact parameter
375 // nuclear thickness function T_B(phi)=T_B(sqtr(s**2+b**2-2*s*b*cos(phi)))
376 //---------------------------------------------------
377 const Double_t phi=x[0];
378 const Double_t rbmax=par[0],b=par[1],s=par[2];
380 const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi));
382 fgWSzB->SetParameter(3,w);
384 return 2.*fgWSzB->Integral(0.,TMath::Sqrt(rbmax*rbmax-w*w));;
387 //**********************************************************************************
388 Double_t AliGenPromptPhotons::TAxTB(const Double_t* x, const Double_t* par) {
389 //---------------------------------------------------
392 // par[0] - RBMAX (fm), max. value of target radius
393 // par[1] - b (fm), impact parameter
396 // s * TA(s) * 2 * Integral(0,phiMax) TB(phi(s,b))
397 //---------------------------------------------------
398 const Double_t s = x[0];
399 const Double_t rbmax=par[0],b=par[1];
403 fgTB->SetParameter(1,b);
404 fgTB->SetParameter(2,s);
407 if(b<rbmax && s<(rbmax-b)) {
410 phiMax=TMath::ACos((s*s+b*b-rbmax*rbmax)/(2.*s*b));
413 return s*fgTA->Eval(s)*2.*fgTB->Integral(0.,phiMax);
416 // ---------------------------------------------------------------------------------
417 Double_t AliGenPromptPhotons::TAB(const Double_t* x, const Double_t* par) {
418 //---------------------------------------------------
420 // x[0] - b (fm), impact parameter
421 // par[0] - RAMAX (fm), max. value of projectile radius
422 // par[1] - RAMAX (fm), max. value of target radius
425 // overlap function TAB(b) (1/fm**2)
426 //---------------------------------------------------
427 const Double_t b=x[0];
428 const Double_t ramax=par[0],rbmax=par[1];
430 Double_t sMin=0.,sMax=ramax;
431 if(b>rbmax) sMin=b-rbmax;
432 if(b<(ramax-rbmax)) sMax=b+rbmax;
434 fgTAxTB->SetParameter(1,b);
436 return fgTAxTB->Integral(sMin,sMax);