Line data Source code
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 6 : ClassImp(AliGenPromptPhotons)
76 :
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;
84 :
85 : //_____________________________________________________________________________
86 : AliGenPromptPhotons::AliGenPromptPhotons()
87 0 : :AliGenerator(-1),
88 0 : fAProjectile(0.),
89 0 : fATarget(0.),
90 0 : fEnergyCMS(0.),
91 0 : fMinImpactParam(0.),
92 0 : fMaxImpactParam(0.)
93 0 : {
94 : //
95 : // Default constructor
96 : //
97 0 : SetCutVertexZ();
98 0 : SetPtRange();
99 0 : SetYRange();
100 0 : }
101 :
102 : //_____________________________________________________________________________
103 : AliGenPromptPhotons::AliGenPromptPhotons(Int_t npart)
104 0 : :AliGenerator(npart),
105 0 : fAProjectile(208),
106 0 : fATarget(208),
107 0 : fEnergyCMS(5500.),
108 0 : fMinImpactParam(0.),
109 0 : fMaxImpactParam(0.)
110 0 : {
111 : //
112 : // Standard constructor
113 : //
114 :
115 0 : fName="PromptPhotons";
116 0 : fTitle="Prompt photons from pp data fit + T_AB";
117 :
118 0 : SetCutVertexZ();
119 0 : SetPtRange();
120 0 : SetYRange();
121 0 : }
122 :
123 : //_____________________________________________________________________________
124 0 : AliGenPromptPhotons::~AliGenPromptPhotons()
125 0 : {
126 : //
127 : // Standard destructor
128 : //
129 0 : delete fgDataPt;
130 0 : delete fgWSzA;
131 0 : delete fgWSzB;
132 0 : delete fgTA;
133 0 : delete fgTB;
134 0 : delete fgTAxTB;
135 0 : delete fgTAB;
136 0 : }
137 :
138 : //_____________________________________________________________________________
139 : void AliGenPromptPhotons::Init()
140 : {
141 : // Initialisation
142 0 : fgDataPt = new TF1("fgDataPt",FitData,fPtMin,fPtMax,1);
143 0 : fgDataPt->SetParameter(0,fEnergyCMS);
144 :
145 : const Double_t d=0.54; // surface thickness (fm)
146 0 : 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 0 : const Double_t ramax=ra+d*TMath::Log((1.-eps+TMath::Exp(-ra/d))/eps);
149 :
150 0 : TF1 fWSforNormA("fWSforNormA",&WSforNorm,0,ramax,2);
151 0 : fWSforNormA.SetParameters(ra,d);
152 0 : fWSforNormA.SetParNames("RA","d");
153 0 : const Double_t ca=1./fWSforNormA.Integral(0.,ramax);
154 :
155 0 : const Double_t rb=1.19*TMath::Power(fATarget,1./3.)-1.61/TMath::Power(fATarget,1./3.);
156 0 : const Double_t rbmax=rb+d*TMath::Log((1.-eps+TMath::Exp(-rb/d))/eps);
157 :
158 0 : TF1 fWSforNormB("fWSforNormB",&WSforNorm,0,rbmax,2);
159 0 : fWSforNormB.SetParameters(rb,d);
160 0 : fWSforNormB.SetParNames("RB","d");
161 0 : const Double_t cb=1./fWSforNormB.Integral(0.,rbmax);
162 :
163 0 : fgWSzA = new TF1("fgWSzA",WSz,0.,ramax,4);
164 0 : fgWSzA->SetParameter(0,ra);
165 0 : fgWSzA->SetParameter(1,d);
166 0 : fgWSzA->SetParameter(2,ca);
167 :
168 0 : fgTA = new TF1("fgTA",TA,0.,ramax,1);
169 0 : fgTA->SetParameter(0,ramax);
170 :
171 0 : fgWSzB = new TF1("fgWSzB",WSz,0.,rbmax,4);
172 0 : fgWSzB->SetParameter(0,rb);
173 0 : fgWSzB->SetParameter(1,d);
174 0 : fgWSzB->SetParameter(2,cb);
175 :
176 0 : fgTB = new TF1("fgTB",TB,0.,TMath::Pi(),3);
177 0 : fgTB->SetParameter(0,rbmax);
178 :
179 0 : fgTAxTB = new TF1("fgTAxTB",TAxTB,0,ramax,2);
180 0 : fgTAxTB->SetParameter(0,rbmax);
181 :
182 0 : fgTAB = new TF1("fgTAB",TAB,0.,ramax+rbmax,2);
183 0 : fgTAB->SetParameter(0,ramax);
184 0 : fgTAB->SetParameter(1,rbmax);
185 :
186 0 : }
187 :
188 : //_____________________________________________________________________________
189 : void AliGenPromptPhotons::Generate()
190 : {
191 : //
192 : // Generate thermal photons of a event
193 : //
194 :
195 0 : Float_t polar[3]= {0,0,0};
196 0 : Float_t origin[3];
197 : Float_t time;
198 0 : Float_t p[3];
199 0 : Float_t random[6];
200 0 : Int_t nt;
201 :
202 0 : for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j];
203 0 : time = fTimeOrigin;
204 : /*
205 : if(fVertexSmear==kPerEvent) {
206 : Vertex();
207 : for (j=0;j<3;j++) origin[j]=fVertex[j];
208 : time = fTime;
209 : }
210 : */
211 0 : TArrayF eventVertex;
212 0 : eventVertex.Set(3);
213 0 : eventVertex[0] = origin[0];
214 0 : eventVertex[1] = origin[1];
215 0 : eventVertex[2] = origin[2];
216 : Float_t eventTime = time;
217 :
218 : Int_t nGam;
219 : Float_t b,pt,rapidity,phi,ww;
220 :
221 0 : b=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm());
222 :
223 0 : Double_t dsdyPP=fgDataPt->Integral(fPtMin,fPtMax); // pb, fm^2 = 1e10 pb
224 0 : ww=dsdyPP*1.0e-10*(fYMax-fYMin)*fAProjectile*fATarget*fgTAB->Eval(b);
225 0 : nGam=Int_t(ww);
226 0 : if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++;
227 :
228 0 : if(nGam) {
229 0 : for(Int_t i=0; i<nGam; i++) {
230 0 : pt=fgDataPt->GetRandom();
231 0 : Rndm(random,2);
232 0 : rapidity=(fYMax-fYMin)*random[0]+fYMin;
233 0 : phi=2.*TMath::Pi()*random[1];
234 0 : p[0]=pt*TMath::Cos(phi);
235 0 : p[1]=pt*TMath::Sin(phi);
236 0 : p[2]=pt*TMath::SinH(rapidity);
237 :
238 0 : if(fVertexSmear==kPerTrack) {
239 0 : Rndm(random,6);
240 0 : for (Int_t j=0;j<3;j++) {
241 0 : origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
242 0 : TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
243 : }
244 0 : Rndm(random,2);
245 0 : time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
246 0 : TMath::Cos(2*random[0]*TMath::Pi())*
247 0 : TMath::Sqrt(-2*TMath::Log(random[1]));
248 0 : }
249 :
250 0 : PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.);
251 : }
252 0 : }
253 :
254 : // Header
255 0 : AliGenEventHeader* header = new AliGenEventHeader("PromptPhotons");
256 : // Event Vertex
257 0 : header->SetPrimaryVertex(eventVertex);
258 0 : header->SetInteractionTime(eventTime);
259 0 : header->SetNProduced(fNpart);
260 0 : gAlice->SetGenEventHeader(header);
261 :
262 0 : }
263 :
264 : void AliGenPromptPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
265 0 : AliGenerator::SetPtRange(ptmin, ptmax);
266 0 : }
267 :
268 : void AliGenPromptPhotons::SetYRange(Float_t ymin, Float_t ymax) {
269 0 : AliGenerator::SetYRange(ymin, ymax);
270 0 : }
271 :
272 : //**********************************************************************************
273 : Double_t AliGenPromptPhotons::FitData(const Double_t* x, const Double_t* par) {
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 0 : const Double_t xT=2.*x[0]/par[0];
294 0 : const Double_t log10xT=TMath::Log10(xT);
295 : Int_t i=0;
296 0 : while(log10xT>log10x[i] && i<nMax) i++;
297 0 : if(i==0) i=1;
298 0 : if(i==nMax) i=nMax-1;
299 0 : const Double_t alpha=(log10F[i]-log10F[i-1])/(log10x[i]-log10x[i-1]);
300 0 : const Double_t beta=log10F[i-1]-alpha*log10x[i-1];
301 :
302 0 : return (TMath::TwoPi()*x[0])*TMath::Power(10.,alpha*log10xT+beta)/TMath::Power(par[0],5.);
303 : }
304 :
305 : //**********************************************************************************
306 : Double_t AliGenPromptPhotons::WSforNorm(const Double_t* x, const Double_t* par) {
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 0 : const Double_t r=x[0];
322 0 : const Double_t bigR=par[0],d=par[1];
323 :
324 0 : return 4.*TMath::Pi()*r*r/(1.+TMath::Exp((r-bigR)/d));
325 : }
326 :
327 : //**********************************************************************************
328 : Double_t AliGenPromptPhotons::WSz(const Double_t* x, const Double_t* par) {
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 0 : const Double_t z=x[0];
342 0 : const Double_t bigR=par[0],d=par[1],C=par[2],b=par[3];
343 :
344 0 : return C/(1.+TMath::Exp((TMath::Sqrt(b*b+z*z)-bigR)/d));
345 : }
346 :
347 : //**********************************************************************************
348 : Double_t AliGenPromptPhotons::TA(const Double_t* x, const Double_t* par) {
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 0 : const Double_t b=x[0];
358 0 : const Double_t ramax=par[0];
359 :
360 0 : fgWSzA->SetParameter(3,b);
361 :
362 0 : return 2.*fgWSzA->Integral(0.,TMath::Sqrt(ramax*ramax-b*b));
363 : }
364 :
365 : //**********************************************************************************
366 : Double_t AliGenPromptPhotons::TB(const Double_t* x, const Double_t* par) {
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 0 : const Double_t phi=x[0];
378 0 : const Double_t rbmax=par[0],b=par[1],s=par[2];
379 :
380 0 : const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi));
381 :
382 0 : fgWSzB->SetParameter(3,w);
383 :
384 0 : return 2.*fgWSzB->Integral(0.,TMath::Sqrt(rbmax*rbmax-w*w));;
385 : }
386 :
387 : //**********************************************************************************
388 : Double_t AliGenPromptPhotons::TAxTB(const Double_t* x, const Double_t* par) {
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 0 : const Double_t s = x[0];
399 0 : const Double_t rbmax=par[0],b=par[1];
400 :
401 0 : if(s==0.) return 0.;
402 :
403 0 : fgTB->SetParameter(1,b);
404 0 : fgTB->SetParameter(2,s);
405 :
406 : Double_t phiMax;
407 0 : if(b<rbmax && s<(rbmax-b)) {
408 0 : phiMax=TMath::Pi();
409 0 : } else {
410 0 : phiMax=TMath::ACos((s*s+b*b-rbmax*rbmax)/(2.*s*b));
411 : }
412 :
413 0 : return s*fgTA->Eval(s)*2.*fgTB->Integral(0.,phiMax);
414 0 : }
415 :
416 : // ---------------------------------------------------------------------------------
417 : Double_t AliGenPromptPhotons::TAB(const Double_t* x, const Double_t* par) {
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 0 : const Double_t b=x[0];
428 0 : const Double_t ramax=par[0],rbmax=par[1];
429 :
430 : Double_t sMin=0.,sMax=ramax;
431 0 : if(b>rbmax) sMin=b-rbmax;
432 0 : if(b<(ramax-rbmax)) sMax=b+rbmax;
433 :
434 0 : fgTAxTB->SetParameter(1,b);
435 :
436 0 : return fgTAxTB->Integral(sMin,sMax);
437 : }
|
