Warnings corrected.
[u/mrichter/AliRoot.git] / EVGEN / AliGenPromptPhotons.cxx
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::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     :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 fgDataPt;
130     delete fgWSzA;
131     delete fgWSzB;
132     delete fgTA;
133     delete fgTB;
134     delete fgTAxTB;
135     delete fgTAB;
136 }
137
138 //_____________________________________________________________________________
139 void AliGenPromptPhotons::Init()
140 {
141
142   fgDataPt = new TF1("fgDataPt",FitData,fPtMin,fPtMax,1);
143   fgDataPt->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   fgWSzA = new TF1("fgWSzA",WSz,0.,ramax,4);
164   fgWSzA->SetParameter(0,ra);
165   fgWSzA->SetParameter(1,d);
166   fgWSzA->SetParameter(2,ca);
167
168   fgTA = new TF1("fgTA",TA,0.,ramax,1);
169   fgTA->SetParameter(0,ramax);
170
171   fgWSzB = new TF1("fgWSzB",WSz,0.,rbmax,4);
172   fgWSzB->SetParameter(0,rb);
173   fgWSzB->SetParameter(1,d);
174   fgWSzB->SetParameter(2,cb);
175
176   fgTB = new TF1("fgTB",TB,0.,TMath::Pi(),3);
177   fgTB->SetParameter(0,rbmax);
178
179   fgTAxTB = new TF1("fgTAxTB",TAxTB,0,ramax,2);
180   fgTAxTB->SetParameter(0,rbmax);
181
182   fgTAB = new TF1("fgTAB",TAB,0.,ramax+rbmax,2);
183   fgTAB->SetParameter(0,ramax);
184   fgTAB->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=fgDataPt->Integral(fPtMin,fPtMax); // pb, fm^2 = 1e10 pb
220     ww=dsdyPP*1.0e-10*(fYMax-fYMin)*fAProjectile*fATarget*fgTAB->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=fgDataPt->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 bigR=par[0],d=par[1];
314
315   return 4.*TMath::Pi()*r*r/(1.+TMath::Exp((r-bigR)/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 bigR=par[0],d=par[1],C=par[2],b=par[3];
334
335   return C/(1.+TMath::Exp((TMath::Sqrt(b*b+z*z)-bigR)/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   fgWSzA->SetParameter(3,b);
352
353   return 2.*fgWSzA->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   fgWSzB->SetParameter(3,w);
374
375   return 2.*fgWSzB->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   fgTB->SetParameter(1,b);
395   fgTB->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*fgTA->Eval(s)*2.*fgTB->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   fgTAxTB->SetParameter(1,b);
426
427   return fgTAxTB->Integral(sMin,sMax);
428 }