]>
Commit | Line | Data |
---|---|---|
3bc7f4d6 | 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 prompt photons for the reaction A+B, sqrt(S) | |
22 | // | |
23 | // main assumptions: | |
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 | |
35 | // | |
36 | // input parameters: | |
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 | |
41 | // | |
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 | //------------------------------------------------------------------------- | |
50 | ||
51 | //Begin_Html | |
52 | /* | |
53 | <img src="picts/AliGeneratorClass.gif"> | |
54 | </pre> | |
55 | <br clear=left> | |
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>. | |
59 | </font> | |
60 | <pre> | |
61 | */ | |
62 | //End_Html | |
63 | // // | |
64 | /////////////////////////////////////////////////////////////////// | |
65 | ||
66 | #include <TArrayF.h> | |
67 | #include <TF1.h> | |
68 | ||
69 | #include "AliConst.h" | |
70 | #include "AliGenEventHeader.h" | |
71 | #include "AliGenPromptPhotons.h" | |
72 | #include "AliLog.h" | |
73 | #include "AliRun.h" | |
74 | ||
75 | ClassImp(AliGenPromptPhotons) | |
76 | ||
00d93180 | 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; | |
3bc7f4d6 | 84 | |
85 | //_____________________________________________________________________________ | |
86 | AliGenPromptPhotons::AliGenPromptPhotons() | |
87 | :AliGenerator(-1), | |
88 | fAProjectile(0.), | |
89 | fATarget(0.), | |
90 | fEnergyCMS(0.), | |
91 | fMinImpactParam(0.), | |
92 | fMaxImpactParam(0.) | |
93 | { | |
94 | // | |
95 | // Default constructor | |
96 | // | |
97 | SetCutVertexZ(); | |
98 | SetPtRange(); | |
99 | SetYRange(); | |
100 | } | |
101 | ||
102 | //_____________________________________________________________________________ | |
103 | AliGenPromptPhotons::AliGenPromptPhotons(Int_t npart) | |
104 | :AliGenerator(npart), | |
105 | fAProjectile(208), | |
106 | fATarget(208), | |
107 | fEnergyCMS(5500.), | |
108 | fMinImpactParam(0.), | |
109 | fMaxImpactParam(0.) | |
110 | { | |
111 | // | |
112 | // Standard constructor | |
113 | // | |
114 | ||
115 | fName="PromptPhotons"; | |
116 | fTitle="Prompt photons from pp data fit + T_AB"; | |
117 | ||
118 | SetCutVertexZ(); | |
119 | SetPtRange(); | |
120 | SetYRange(); | |
121 | } | |
122 | ||
123 | //_____________________________________________________________________________ | |
124 | AliGenPromptPhotons::~AliGenPromptPhotons() | |
125 | { | |
126 | // | |
127 | // Standard destructor | |
128 | // | |
00d93180 | 129 | delete fgDataPt; |
130 | delete fgWSzA; | |
131 | delete fgWSzB; | |
132 | delete fgTA; | |
133 | delete fgTB; | |
134 | delete fgTAxTB; | |
135 | delete fgTAB; | |
3bc7f4d6 | 136 | } |
137 | ||
138 | //_____________________________________________________________________________ | |
139 | void AliGenPromptPhotons::Init() | |
140 | { | |
141 | ||
00d93180 | 142 | fgDataPt = new TF1("fgDataPt",FitData,fPtMin,fPtMax,1); |
143 | fgDataPt->SetParameter(0,fEnergyCMS); | |
3bc7f4d6 | 144 | |
145 | const Double_t d=0.54; // surface thickness (fm) | |
00d93180 | 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); | |
3bc7f4d6 | 149 | |
00d93180 | 150 | TF1 fWSforNormA("fWSforNormA",&WSforNorm,0,ramax,2); |
151 | fWSforNormA.SetParameters(ra,d); | |
3bc7f4d6 | 152 | fWSforNormA.SetParNames("RA","d"); |
00d93180 | 153 | const Double_t ca=1./fWSforNormA.Integral(0.,ramax); |
3bc7f4d6 | 154 | |
00d93180 | 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); | |
3bc7f4d6 | 157 | |
00d93180 | 158 | TF1 fWSforNormB("fWSforNormB",&WSforNorm,0,rbmax,2); |
159 | fWSforNormB.SetParameters(rb,d); | |
3bc7f4d6 | 160 | fWSforNormB.SetParNames("RB","d"); |
00d93180 | 161 | const Double_t cb=1./fWSforNormB.Integral(0.,rbmax); |
3bc7f4d6 | 162 | |
00d93180 | 163 | fgWSzA = new TF1("fgWSzA",WSz,0.,ramax,4); |
164 | fgWSzA->SetParameter(0,ra); | |
165 | fgWSzA->SetParameter(1,d); | |
166 | fgWSzA->SetParameter(2,ca); | |
3bc7f4d6 | 167 | |
00d93180 | 168 | fgTA = new TF1("fgTA",TA,0.,ramax,1); |
169 | fgTA->SetParameter(0,ramax); | |
3bc7f4d6 | 170 | |
00d93180 | 171 | fgWSzB = new TF1("fgWSzB",WSz,0.,rbmax,4); |
172 | fgWSzB->SetParameter(0,rb); | |
173 | fgWSzB->SetParameter(1,d); | |
174 | fgWSzB->SetParameter(2,cb); | |
3bc7f4d6 | 175 | |
00d93180 | 176 | fgTB = new TF1("fgTB",TB,0.,TMath::Pi(),3); |
177 | fgTB->SetParameter(0,rbmax); | |
3bc7f4d6 | 178 | |
00d93180 | 179 | fgTAxTB = new TF1("fgTAxTB",TAxTB,0,ramax,2); |
180 | fgTAxTB->SetParameter(0,rbmax); | |
3bc7f4d6 | 181 | |
00d93180 | 182 | fgTAB = new TF1("fgTAB",TAB,0.,ramax+rbmax,2); |
183 | fgTAB->SetParameter(0,ramax); | |
184 | fgTAB->SetParameter(1,rbmax); | |
3bc7f4d6 | 185 | |
186 | } | |
187 | ||
188 | //_____________________________________________________________________________ | |
189 | void AliGenPromptPhotons::Generate() | |
190 | { | |
191 | // | |
192 | // Generate thermal photons of a event | |
193 | // | |
194 | ||
195 | Float_t polar[3]= {0,0,0}; | |
196 | Float_t origin[3]; | |
197 | Float_t p[3]; | |
198 | Float_t random[6]; | |
199 | Int_t nt; | |
200 | ||
201 | for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j]; | |
202 | /* | |
203 | if(fVertexSmear==kPerEvent) { | |
204 | Vertex(); | |
205 | for (j=0;j<3;j++) origin[j]=fVertex[j]; | |
206 | } | |
207 | */ | |
208 | TArrayF eventVertex; | |
209 | eventVertex.Set(3); | |
210 | eventVertex[0] = origin[0]; | |
211 | eventVertex[1] = origin[1]; | |
212 | eventVertex[2] = origin[2]; | |
213 | ||
214 | Int_t nGam; | |
215 | Float_t b,pt,rapidity,phi,ww; | |
216 | ||
217 | b=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm()); | |
218 | ||
00d93180 | 219 | Double_t dsdyPP=fgDataPt->Integral(fPtMin,fPtMax); // pb, fm^2 = 1e10 pb |
220 | ww=dsdyPP*1.0e-10*(fYMax-fYMin)*fAProjectile*fATarget*fgTAB->Eval(b); | |
3bc7f4d6 | 221 | nGam=Int_t(ww); |
222 | if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++; | |
223 | ||
224 | if(nGam) { | |
225 | for(Int_t i=0; i<nGam; i++) { | |
00d93180 | 226 | pt=fgDataPt->GetRandom(); |
3bc7f4d6 | 227 | Rndm(random,2); |
228 | rapidity=(fYMax-fYMin)*random[0]+fYMin; | |
229 | phi=2.*TMath::Pi()*random[1]; | |
230 | p[0]=pt*TMath::Cos(phi); | |
231 | p[1]=pt*TMath::Sin(phi); | |
232 | p[2]=pt*TMath::SinH(rapidity); | |
233 | ||
234 | if(fVertexSmear==kPerTrack) { | |
235 | Rndm(random,6); | |
236 | for (Int_t j=0;j<3;j++) { | |
237 | origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* | |
238 | TMath::Sqrt(-2*TMath::Log(random[2*j+1])); | |
239 | } | |
240 | } | |
241 | ||
242 | PushTrack(fTrackIt,-1,22,p,origin,polar,0,kPPrimary,nt,1.); | |
243 | } | |
244 | } | |
245 | ||
246 | // Header | |
247 | AliGenEventHeader* header = new AliGenEventHeader("PromptPhotons"); | |
248 | // Event Vertex | |
249 | header->SetPrimaryVertex(eventVertex); | |
250 | header->SetNProduced(fNpart); | |
251 | gAlice->SetGenEventHeader(header); | |
252 | ||
253 | } | |
254 | ||
255 | void AliGenPromptPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) { | |
256 | AliGenerator::SetPtRange(ptmin, ptmax); | |
257 | } | |
258 | ||
259 | void AliGenPromptPhotons::SetYRange(Float_t ymin, Float_t ymax) { | |
260 | AliGenerator::SetYRange(ymin, ymax); | |
261 | } | |
262 | ||
263 | //********************************************************************************** | |
00d93180 | 264 | Double_t AliGenPromptPhotons::FitData(Double_t* x, Double_t* par) { |
3bc7f4d6 | 265 | //--------------------------------------------------- |
266 | // input: | |
267 | // x[0] - p_t (GeV). | |
268 | // par[0]=sqrt(s_NN) (GeV), | |
269 | // | |
270 | // output: | |
271 | // d^{2}#sigma/(dp_t dy) (pb/GeV) | |
272 | //--------------------------------------------------- | |
273 | // | |
274 | // d^{2}#sigma/(dp_t dy) = (2 pi p_t) Ed^{3}#sigma/d^{3}p | |
275 | // | |
276 | // data presentation: Nucl.Phys.A783:577-582,2007, hep-ex/0609037, fig.3 | |
277 | // F(x_t)=(#sqrt{s})^{5} Ed^{3}#sigma/d^{3}p | |
278 | //--------------------------------------------------- | |
279 | // approximate tabulation of F(x_t) | |
280 | const Int_t nMax=4; | |
281 | const Double_t log10x[nMax]={ -2., -1., -0.6, -0.3}; | |
282 | const Double_t log10F[nMax]={ 19., 13., 9.8, 7.}; | |
283 | ||
284 | const Double_t xT=2.*x[0]/par[0]; | |
285 | const Double_t log10xT=TMath::Log10(xT); | |
286 | Int_t i=0; | |
287 | while(log10xT>log10x[i] && i<nMax) i++; | |
288 | if(i==0) i=1; | |
289 | if(i==nMax) i=nMax-1; | |
290 | const Double_t alpha=(log10F[i]-log10F[i-1])/(log10x[i]-log10x[i-1]); | |
291 | const Double_t beta=log10F[i-1]-alpha*log10x[i-1]; | |
292 | ||
293 | return (TMath::TwoPi()*x[0])*TMath::Power(10.,alpha*log10xT+beta)/TMath::Power(par[0],5.); | |
294 | } | |
295 | ||
296 | //********************************************************************************** | |
297 | Double_t AliGenPromptPhotons::WSforNorm(Double_t* x, Double_t* par) { | |
298 | //--------------------------------------------------- | |
299 | // input: | |
300 | // x[0] - r (fm) | |
301 | // par[0] - R (fm), radius | |
302 | // par[1] - d (fm), surface thickness | |
303 | // | |
304 | // output: | |
305 | // 4 pi r**2 /(1+exp((r-R)/d)) | |
306 | // | |
307 | // Wood Saxon (WS) C/(1+exp((r-RA)/d)) (nuclons/fm^3) | |
308 | // To get the normalization A = (Integral 4 pi r**2 dr WS): | |
309 | // C = A / (Integral 4 pi r**2 dr 1/(1+exp((r-RA)/d)) ) | |
310 | // Thus me need 4 pi r**2 /(1+exp((r-RA)/d)) (1/fm) | |
311 | //--------------------------------------------------- | |
312 | const Double_t r=x[0]; | |
00d93180 | 313 | const Double_t bigR=par[0],d=par[1]; |
3bc7f4d6 | 314 | |
00d93180 | 315 | return 4.*TMath::Pi()*r*r/(1.+TMath::Exp((r-bigR)/d)); |
3bc7f4d6 | 316 | } |
317 | ||
318 | //********************************************************************************** | |
319 | Double_t AliGenPromptPhotons::WSz(Double_t* x, Double_t* par) { | |
320 | //--------------------------------------------------- | |
321 | // input: | |
322 | // x[0] - z (fm) | |
323 | // par[0] - R (fm), radius | |
324 | // par[1] - d (fm), surface thickness | |
325 | // par[2] - C (nucleons/fm**2), normalization factor | |
326 | // par[3] - b (fm), impact parameter | |
327 | // | |
328 | // output: | |
329 | // Wood Saxon Parameterisation | |
330 | // as a function of z for fixed b (1/fm^3) | |
331 | //--------------------------------------------------- | |
332 | const Double_t z=x[0]; | |
00d93180 | 333 | const Double_t bigR=par[0],d=par[1],C=par[2],b=par[3]; |
3bc7f4d6 | 334 | |
00d93180 | 335 | return C/(1.+TMath::Exp((TMath::Sqrt(b*b+z*z)-bigR)/d)); |
3bc7f4d6 | 336 | } |
337 | ||
338 | //********************************************************************************** | |
339 | Double_t AliGenPromptPhotons::TA(Double_t* x, Double_t* par) { | |
340 | //--------------------------------------------------- | |
341 | // input: | |
342 | // x[0] - b (fm), impact parameter | |
343 | // par[0] - RAMAX (fm), max. value of projectile radius | |
344 | // | |
345 | // output: | |
346 | // nuclear thickness function T_A(b) (1/fm^2) | |
347 | //--------------------------------------------------- | |
348 | const Double_t b=x[0]; | |
00d93180 | 349 | const Double_t ramax=par[0]; |
3bc7f4d6 | 350 | |
00d93180 | 351 | fgWSzA->SetParameter(3,b); |
3bc7f4d6 | 352 | |
00d93180 | 353 | return 2.*fgWSzA->Integral(0.,TMath::Sqrt(ramax*ramax-b*b)); |
3bc7f4d6 | 354 | } |
355 | ||
356 | //********************************************************************************** | |
357 | Double_t AliGenPromptPhotons::TB(Double_t* x, Double_t* par) { | |
358 | //--------------------------------------------------- | |
359 | // input: | |
360 | // x[0] - phi (rad) | |
361 | // par[0] - RBMAX (fm), max. value of target radius | |
362 | // par[1] - b (fm), impact parameter | |
363 | // par[2] - s (fm) | |
364 | // | |
365 | // output: | |
366 | // nuclear thickness function T_B(phi)=T_B(sqtr(s**2+b**2-2*s*b*cos(phi))) | |
367 | //--------------------------------------------------- | |
368 | const Double_t phi=x[0]; | |
00d93180 | 369 | const Double_t rbmax=par[0],b=par[1],s=par[2]; |
3bc7f4d6 | 370 | |
371 | const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi)); | |
372 | ||
00d93180 | 373 | fgWSzB->SetParameter(3,w); |
3bc7f4d6 | 374 | |
00d93180 | 375 | return 2.*fgWSzB->Integral(0.,TMath::Sqrt(rbmax*rbmax-w*w));; |
3bc7f4d6 | 376 | } |
377 | ||
378 | //********************************************************************************** | |
379 | Double_t AliGenPromptPhotons::TAxTB(Double_t* x, Double_t* par) { | |
380 | //--------------------------------------------------- | |
381 | // input: | |
382 | // x[0] - s (fm) | |
383 | // par[0] - RBMAX (fm), max. value of target radius | |
384 | // par[1] - b (fm), impact parameter | |
385 | // | |
386 | // output: | |
387 | // s * TA(s) * 2 * Integral(0,phiMax) TB(phi(s,b)) | |
388 | //--------------------------------------------------- | |
389 | const Double_t s = x[0]; | |
00d93180 | 390 | const Double_t rbmax=par[0],b=par[1]; |
3bc7f4d6 | 391 | |
392 | if(s==0.) return 0.; | |
393 | ||
00d93180 | 394 | fgTB->SetParameter(1,b); |
395 | fgTB->SetParameter(2,s); | |
3bc7f4d6 | 396 | |
397 | Double_t phiMax; | |
00d93180 | 398 | if(b<rbmax && s<(rbmax-b)) { |
3bc7f4d6 | 399 | phiMax=TMath::Pi(); |
400 | } else { | |
00d93180 | 401 | phiMax=TMath::ACos((s*s+b*b-rbmax*rbmax)/(2.*s*b)); |
3bc7f4d6 | 402 | } |
403 | ||
00d93180 | 404 | return s*fgTA->Eval(s)*2.*fgTB->Integral(0.,phiMax); |
3bc7f4d6 | 405 | } |
406 | ||
407 | // --------------------------------------------------------------------------------- | |
408 | Double_t AliGenPromptPhotons::TAB(Double_t* x, Double_t* par) { | |
409 | //--------------------------------------------------- | |
410 | // input: | |
411 | // x[0] - b (fm), impact parameter | |
412 | // par[0] - RAMAX (fm), max. value of projectile radius | |
413 | // par[1] - RAMAX (fm), max. value of target radius | |
414 | // | |
415 | // output: | |
416 | // overlap function TAB(b) (1/fm**2) | |
417 | //--------------------------------------------------- | |
418 | const Double_t b=x[0]; | |
00d93180 | 419 | const Double_t ramax=par[0],rbmax=par[1]; |
3bc7f4d6 | 420 | |
00d93180 | 421 | Double_t sMin=0.,sMax=ramax; |
422 | if(b>rbmax) sMin=b-rbmax; | |
423 | if(b<(ramax-rbmax)) sMax=b+rbmax; | |
3bc7f4d6 | 424 | |
00d93180 | 425 | fgTAxTB->SetParameter(1,b); |
3bc7f4d6 | 426 | |
00d93180 | 427 | return fgTAxTB->Integral(sMin,sMax); |
3bc7f4d6 | 428 | } |