]>
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 | { | |
4a33c50d | 141 | // Initialisation |
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]; | |
21391258 | 197 | Float_t time; |
3bc7f4d6 | 198 | Float_t p[3]; |
199 | Float_t random[6]; | |
200 | Int_t nt; | |
201 | ||
202 | for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j]; | |
21391258 | 203 | time = fTimeOrigin; |
3bc7f4d6 | 204 | /* |
205 | if(fVertexSmear==kPerEvent) { | |
206 | Vertex(); | |
207 | for (j=0;j<3;j++) origin[j]=fVertex[j]; | |
21391258 | 208 | time = fTime; |
3bc7f4d6 | 209 | } |
210 | */ | |
211 | TArrayF eventVertex; | |
212 | eventVertex.Set(3); | |
213 | eventVertex[0] = origin[0]; | |
214 | eventVertex[1] = origin[1]; | |
215 | eventVertex[2] = origin[2]; | |
21391258 | 216 | Float_t eventTime = time; |
3bc7f4d6 | 217 | |
218 | Int_t nGam; | |
219 | Float_t b,pt,rapidity,phi,ww; | |
220 | ||
221 | b=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm()); | |
222 | ||
00d93180 | 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); | |
3bc7f4d6 | 225 | nGam=Int_t(ww); |
226 | if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++; | |
227 | ||
228 | if(nGam) { | |
229 | for(Int_t i=0; i<nGam; i++) { | |
00d93180 | 230 | pt=fgDataPt->GetRandom(); |
3bc7f4d6 | 231 | Rndm(random,2); |
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); | |
237 | ||
238 | if(fVertexSmear==kPerTrack) { | |
239 | Rndm(random,6); | |
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])); | |
243 | } | |
21391258 | 244 | Rndm(random,2); |
245 | time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()* | |
246 | TMath::Cos(2*random[0]*TMath::Pi())* | |
247 | TMath::Sqrt(-2*TMath::Log(random[1])); | |
3bc7f4d6 | 248 | } |
249 | ||
21391258 | 250 | PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.); |
3bc7f4d6 | 251 | } |
252 | } | |
253 | ||
254 | // Header | |
255 | AliGenEventHeader* header = new AliGenEventHeader("PromptPhotons"); | |
256 | // Event Vertex | |
257 | header->SetPrimaryVertex(eventVertex); | |
21391258 | 258 | header->SetInteractionTime(eventTime); |
3bc7f4d6 | 259 | header->SetNProduced(fNpart); |
260 | gAlice->SetGenEventHeader(header); | |
261 | ||
262 | } | |
263 | ||
264 | void AliGenPromptPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) { | |
265 | AliGenerator::SetPtRange(ptmin, ptmax); | |
266 | } | |
267 | ||
268 | void AliGenPromptPhotons::SetYRange(Float_t ymin, Float_t ymax) { | |
269 | AliGenerator::SetYRange(ymin, ymax); | |
270 | } | |
271 | ||
272 | //********************************************************************************** | |
4a33c50d | 273 | Double_t AliGenPromptPhotons::FitData(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 274 | //--------------------------------------------------- |
275 | // input: | |
276 | // x[0] - p_t (GeV). | |
277 | // par[0]=sqrt(s_NN) (GeV), | |
278 | // | |
279 | // output: | |
280 | // d^{2}#sigma/(dp_t dy) (pb/GeV) | |
281 | //--------------------------------------------------- | |
282 | // | |
283 | // d^{2}#sigma/(dp_t dy) = (2 pi p_t) Ed^{3}#sigma/d^{3}p | |
284 | // | |
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) | |
289 | const Int_t nMax=4; | |
290 | const Double_t log10x[nMax]={ -2., -1., -0.6, -0.3}; | |
291 | const Double_t log10F[nMax]={ 19., 13., 9.8, 7.}; | |
292 | ||
293 | const Double_t xT=2.*x[0]/par[0]; | |
294 | const Double_t log10xT=TMath::Log10(xT); | |
295 | Int_t i=0; | |
296 | while(log10xT>log10x[i] && i<nMax) i++; | |
297 | if(i==0) i=1; | |
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]; | |
301 | ||
302 | return (TMath::TwoPi()*x[0])*TMath::Power(10.,alpha*log10xT+beta)/TMath::Power(par[0],5.); | |
303 | } | |
304 | ||
305 | //********************************************************************************** | |
4a33c50d | 306 | Double_t AliGenPromptPhotons::WSforNorm(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 307 | //--------------------------------------------------- |
308 | // input: | |
309 | // x[0] - r (fm) | |
310 | // par[0] - R (fm), radius | |
311 | // par[1] - d (fm), surface thickness | |
312 | // | |
313 | // output: | |
314 | // 4 pi r**2 /(1+exp((r-R)/d)) | |
315 | // | |
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]; | |
00d93180 | 322 | const Double_t bigR=par[0],d=par[1]; |
3bc7f4d6 | 323 | |
00d93180 | 324 | return 4.*TMath::Pi()*r*r/(1.+TMath::Exp((r-bigR)/d)); |
3bc7f4d6 | 325 | } |
326 | ||
327 | //********************************************************************************** | |
4a33c50d | 328 | Double_t AliGenPromptPhotons::WSz(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 329 | //--------------------------------------------------- |
330 | // input: | |
331 | // x[0] - z (fm) | |
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 | |
336 | // | |
337 | // output: | |
338 | // Wood Saxon Parameterisation | |
339 | // as a function of z for fixed b (1/fm^3) | |
340 | //--------------------------------------------------- | |
341 | const Double_t z=x[0]; | |
00d93180 | 342 | const Double_t bigR=par[0],d=par[1],C=par[2],b=par[3]; |
3bc7f4d6 | 343 | |
00d93180 | 344 | return C/(1.+TMath::Exp((TMath::Sqrt(b*b+z*z)-bigR)/d)); |
3bc7f4d6 | 345 | } |
346 | ||
347 | //********************************************************************************** | |
4a33c50d | 348 | Double_t AliGenPromptPhotons::TA(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 349 | //--------------------------------------------------- |
350 | // input: | |
351 | // x[0] - b (fm), impact parameter | |
352 | // par[0] - RAMAX (fm), max. value of projectile radius | |
353 | // | |
354 | // output: | |
355 | // nuclear thickness function T_A(b) (1/fm^2) | |
356 | //--------------------------------------------------- | |
357 | const Double_t b=x[0]; | |
00d93180 | 358 | const Double_t ramax=par[0]; |
3bc7f4d6 | 359 | |
00d93180 | 360 | fgWSzA->SetParameter(3,b); |
3bc7f4d6 | 361 | |
00d93180 | 362 | return 2.*fgWSzA->Integral(0.,TMath::Sqrt(ramax*ramax-b*b)); |
3bc7f4d6 | 363 | } |
364 | ||
365 | //********************************************************************************** | |
4a33c50d | 366 | Double_t AliGenPromptPhotons::TB(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 367 | //--------------------------------------------------- |
368 | // input: | |
369 | // x[0] - phi (rad) | |
370 | // par[0] - RBMAX (fm), max. value of target radius | |
371 | // par[1] - b (fm), impact parameter | |
372 | // par[2] - s (fm) | |
373 | // | |
374 | // output: | |
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]; | |
00d93180 | 378 | const Double_t rbmax=par[0],b=par[1],s=par[2]; |
3bc7f4d6 | 379 | |
380 | const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi)); | |
381 | ||
00d93180 | 382 | fgWSzB->SetParameter(3,w); |
3bc7f4d6 | 383 | |
00d93180 | 384 | return 2.*fgWSzB->Integral(0.,TMath::Sqrt(rbmax*rbmax-w*w));; |
3bc7f4d6 | 385 | } |
386 | ||
387 | //********************************************************************************** | |
4a33c50d | 388 | Double_t AliGenPromptPhotons::TAxTB(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 389 | //--------------------------------------------------- |
390 | // input: | |
391 | // x[0] - s (fm) | |
392 | // par[0] - RBMAX (fm), max. value of target radius | |
393 | // par[1] - b (fm), impact parameter | |
394 | // | |
395 | // output: | |
396 | // s * TA(s) * 2 * Integral(0,phiMax) TB(phi(s,b)) | |
397 | //--------------------------------------------------- | |
398 | const Double_t s = x[0]; | |
00d93180 | 399 | const Double_t rbmax=par[0],b=par[1]; |
3bc7f4d6 | 400 | |
401 | if(s==0.) return 0.; | |
402 | ||
00d93180 | 403 | fgTB->SetParameter(1,b); |
404 | fgTB->SetParameter(2,s); | |
3bc7f4d6 | 405 | |
406 | Double_t phiMax; | |
00d93180 | 407 | if(b<rbmax && s<(rbmax-b)) { |
3bc7f4d6 | 408 | phiMax=TMath::Pi(); |
409 | } else { | |
00d93180 | 410 | phiMax=TMath::ACos((s*s+b*b-rbmax*rbmax)/(2.*s*b)); |
3bc7f4d6 | 411 | } |
412 | ||
00d93180 | 413 | return s*fgTA->Eval(s)*2.*fgTB->Integral(0.,phiMax); |
3bc7f4d6 | 414 | } |
415 | ||
416 | // --------------------------------------------------------------------------------- | |
4a33c50d | 417 | Double_t AliGenPromptPhotons::TAB(const Double_t* x, const Double_t* par) { |
3bc7f4d6 | 418 | //--------------------------------------------------- |
419 | // input: | |
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 | |
423 | // | |
424 | // output: | |
425 | // overlap function TAB(b) (1/fm**2) | |
426 | //--------------------------------------------------- | |
427 | const Double_t b=x[0]; | |
00d93180 | 428 | const Double_t ramax=par[0],rbmax=par[1]; |
3bc7f4d6 | 429 | |
00d93180 | 430 | Double_t sMin=0.,sMax=ramax; |
431 | if(b>rbmax) sMin=b-rbmax; | |
432 | if(b<(ramax-rbmax)) sMax=b+rbmax; | |
3bc7f4d6 | 433 | |
00d93180 | 434 | fgTAxTB->SetParameter(1,b); |
3bc7f4d6 | 435 | |
00d93180 | 436 | return fgTAxTB->Integral(sMin,sMax); |
3bc7f4d6 | 437 | } |