fix flow calculation in scalar product
[u/mrichter/AliRoot.git] / EVGEN / AliGenPromptPhotons.cxx
CommitLineData
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
75ClassImp(AliGenPromptPhotons)
76
77TF1* AliGenPromptPhotons::fDataPt = NULL;
78TF1* AliGenPromptPhotons::fWSzA = NULL;
79TF1* AliGenPromptPhotons::fWSzB = NULL;
80TF1* AliGenPromptPhotons::fTA = NULL;
81TF1* AliGenPromptPhotons::fTB = NULL;
82TF1* AliGenPromptPhotons::fTAxTB = NULL;
83TF1* AliGenPromptPhotons::fTAB = NULL;
84
85//_____________________________________________________________________________
86AliGenPromptPhotons::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//_____________________________________________________________________________
103AliGenPromptPhotons::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//_____________________________________________________________________________
124AliGenPromptPhotons::~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//_____________________________________________________________________________
139void 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//_____________________________________________________________________________
189void 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
255void AliGenPromptPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) {
256 AliGenerator::SetPtRange(ptmin, ptmax);
257}
258
259void AliGenPromptPhotons::SetYRange(Float_t ymin, Float_t ymax) {
260 AliGenerator::SetYRange(ymin, ymax);
261}
262
263//**********************************************************************************
264Double_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//**********************************************************************************
297Double_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//**********************************************************************************
319Double_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//**********************************************************************************
339Double_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//**********************************************************************************
357Double_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//**********************************************************************************
379Double_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// ---------------------------------------------------------------------------------
408Double_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}