]>
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 | ||
77 | TF1* AliGenPromptPhotons::fDataPt = NULL; | |
78 | TF1* AliGenPromptPhotons::fWSzA = NULL; | |
79 | TF1* AliGenPromptPhotons::fWSzB = NULL; | |
80 | TF1* AliGenPromptPhotons::fTA = NULL; | |
81 | TF1* AliGenPromptPhotons::fTB = NULL; | |
82 | TF1* AliGenPromptPhotons::fTAxTB = NULL; | |
83 | TF1* AliGenPromptPhotons::fTAB = NULL; | |
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 | // | |
129 | delete fDataPt; | |
130 | delete fWSzA; | |
131 | delete fWSzB; | |
132 | delete fTA; | |
133 | delete fTB; | |
134 | delete fTAxTB; | |
135 | delete fTAB; | |
136 | } | |
137 | ||
138 | //_____________________________________________________________________________ | |
139 | void AliGenPromptPhotons::Init() | |
140 | { | |
141 | ||
142 | fDataPt = new TF1("fDataPt",fitData,fPtMin,fPtMax,1); | |
143 | fDataPt->SetParameter(0,fEnergyCMS); | |
144 | ||
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); | |
149 | ||
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); | |
154 | ||
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); | |
157 | ||
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); | |
162 | ||
163 | fWSzA = new TF1("fWSzA",WSz,0.,RAMAX,4); | |
164 | fWSzA->SetParameter(0,RA); | |
165 | fWSzA->SetParameter(1,d); | |
166 | fWSzA->SetParameter(2,CA); | |
167 | ||
168 | fTA = new TF1("fTA",TA,0.,RAMAX,1); | |
169 | fTA->SetParameter(0,RAMAX); | |
170 | ||
171 | fWSzB = new TF1("fWSzB",WSz,0.,RBMAX,4); | |
172 | fWSzB->SetParameter(0,RB); | |
173 | fWSzB->SetParameter(1,d); | |
174 | fWSzB->SetParameter(2,CB); | |
175 | ||
176 | fTB = new TF1("fTB",TB,0.,TMath::Pi(),3); | |
177 | fTB->SetParameter(0,RBMAX); | |
178 | ||
179 | fTAxTB = new TF1("fTAxTB",TAxTB,0,RAMAX,2); | |
180 | fTAxTB->SetParameter(0,RBMAX); | |
181 | ||
182 | fTAB = new TF1("fTAB",TAB,0.,RAMAX+RBMAX,2); | |
183 | fTAB->SetParameter(0,RAMAX); | |
184 | fTAB->SetParameter(1,RBMAX); | |
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 | ||
219 | Double_t dsdyPP=fDataPt->Integral(fPtMin,fPtMax); // pb, fm^2 = 1e10 pb | |
220 | ww=dsdyPP*1.0e-10*(fYMax-fYMin)*fAProjectile*fATarget*fTAB->Eval(b); | |
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++) { | |
226 | pt=fDataPt->GetRandom(); | |
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 | //********************************************************************************** | |
264 | Double_t AliGenPromptPhotons::fitData(Double_t* x, Double_t* par) { | |
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]; | |
313 | const Double_t R=par[0],d=par[1]; | |
314 | ||
315 | return 4.*TMath::Pi()*r*r/(1.+TMath::Exp((r-R)/d)); | |
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]; | |
333 | const Double_t R=par[0],d=par[1],C=par[2],b=par[3]; | |
334 | ||
335 | return C/(1.+TMath::Exp((TMath::Sqrt(b*b+z*z)-R)/d)); | |
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]; | |
349 | const Double_t RAMAX=par[0]; | |
350 | ||
351 | fWSzA->SetParameter(3,b); | |
352 | ||
353 | return 2.*fWSzA->Integral(0.,TMath::Sqrt(RAMAX*RAMAX-b*b)); | |
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]; | |
369 | const Double_t RBMAX=par[0],b=par[1],s=par[2]; | |
370 | ||
371 | const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi)); | |
372 | ||
373 | fWSzB->SetParameter(3,w); | |
374 | ||
375 | return 2.*fWSzB->Integral(0.,TMath::Sqrt(RBMAX*RBMAX-w*w));; | |
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]; | |
390 | const Double_t RBMAX=par[0],b=par[1]; | |
391 | ||
392 | if(s==0.) return 0.; | |
393 | ||
394 | fTB->SetParameter(1,b); | |
395 | fTB->SetParameter(2,s); | |
396 | ||
397 | Double_t phiMax; | |
398 | if(b<RBMAX && s<(RBMAX-b)) { | |
399 | phiMax=TMath::Pi(); | |
400 | } else { | |
401 | phiMax=TMath::ACos((s*s+b*b-RBMAX*RBMAX)/(2.*s*b)); | |
402 | } | |
403 | ||
404 | return s*fTA->Eval(s)*2.*fTB->Integral(0.,phiMax); | |
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]; | |
419 | const Double_t RAMAX=par[0],RBMAX=par[1]; | |
420 | ||
421 | Double_t sMin=0.,sMax=RAMAX; | |
422 | if(b>RBMAX) sMin=b-RBMAX; | |
423 | if(b<(RAMAX-RBMAX)) sMax=b+RBMAX; | |
424 | ||
425 | fTAxTB->SetParameter(1,b); | |
426 | ||
427 | return fTAxTB->Integral(sMin,sMax); | |
428 | } |