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   // Initialisation
 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 time;
 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];
 203     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     TArrayF eventVertex;
 212     eventVertex.Set(3);
 213     eventVertex[0] = origin[0];
 214     eventVertex[1] = origin[1];
 215     eventVertex[2] = origin[2];
 216     Float_t eventTime = time;
 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
 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);
 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++) {
 230           pt=fgDataPt->GetRandom();
 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             }
 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]));
 248           }
 249
 250           PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.);
 251         }
 252       }
 253
 254 // Header
 255     AliGenEventHeader* header = new AliGenEventHeader("PromptPhotons");
 256 // Event Vertex
 257     header->SetPrimaryVertex(eventVertex);
 258     header->SetInteractionTime(eventTime);
 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 //**********************************************************************************
 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   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 //**********************************************************************************
 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   const Double_t r=x[0];
 322   const Double_t bigR=par[0],d=par[1];
 323
 324   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   const Double_t z=x[0];
 342   const Double_t bigR=par[0],d=par[1],C=par[2],b=par[3];
 343
 344   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   const Double_t b=x[0];
 358   const Double_t ramax=par[0];
 359
 360   fgWSzA->SetParameter(3,b);
 361
 362   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   const Double_t phi=x[0];
 378   const Double_t rbmax=par[0],b=par[1],s=par[2];
 379
 380   const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi));
 381
 382   fgWSzB->SetParameter(3,w);
 383
 384   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   const Double_t s  = x[0];
 399   const Double_t rbmax=par[0],b=par[1];
 400
 401   if(s==0.) return 0.;
 402
 403   fgTB->SetParameter(1,b);
 404   fgTB->SetParameter(2,s);
 405
 406   Double_t phiMax;
 407   if(b<rbmax && s<(rbmax-b)) {
 408     phiMax=TMath::Pi();
 409   } else {
 410     phiMax=TMath::ACos((s*s+b*b-rbmax*rbmax)/(2.*s*b));
 411   }
 412
 413   return s*fgTA->Eval(s)*2.*fgTB->Integral(0.,phiMax);
 414 }
 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   const Double_t b=x[0];
 428   const Double_t ramax=par[0],rbmax=par[1];
 429
 430   Double_t sMin=0.,sMax=ramax;
 431   if(b>rbmax) sMin=b-rbmax;
 432   if(b<(ramax-rbmax)) sMax=b+rbmax;
 433
 434   fgTAxTB->SetParameter(1,b);
 435
 436   return fgTAxTB->Integral(sMin,sMax);
 437 }