/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
* Contributors are mentioned in the code where appropriate. *
* *
* Permission to use, copy, modify and distribute this software and its *
* documentation strictly for non-commercial purposes is hereby granted *
* without fee, provided that the above copyright notice appears in all *
* copies and that both the copyright notice and this permission notice *
* appear in the supporting documentation. The authors make no claims *
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
//-------------------------------------------------------------------------
// author: Sergey Kiselev, ITEP, Moscow
// e-mail: Sergey.Kiselev@cern.ch
// tel.: 007 495 129 95 45
//-------------------------------------------------------------------------
// Generator of prompt photons for the reaction A+B, sqrt(S)
//
// main assumptions:
// 1. flat rapidity distribution
// 2. all existing p+p(pbar) data at y_{c.m.} can be described by the function
// F(x_T) = (sqrt(s))^5 Ed^3sigma/d^3p, x_T = 2p_t/sqrt(s)
// all data points cover the region x_T: 0.01 - 0.6
// see Nucl.Phys.A783:577-582,2007, hep-ex/0609037
// 3. binary scaling: for A+B at the impact parameter b
// Ed^3N^{AB}(b)/d^3p = Ed^3sigma^{pp}/d^3p A B T_{AB}(b),
// T_{AB}(b) - nuclear overlapping fuction, calculated in the Glauber approach,
// nuclear density is parametrized by a Woods-Saxon with nuclear radius
// R_A = 1.19 A^{1/3} - 1.61 A^{-1/3} fm and surface thickness a=0.54 fm
// 4. nuclear effects (Cronin, shadowing, ...) are ignored
//
// input parameters:
// fAProjectile, fATarget - number of nucleons in a nucleus A and B
// fMinImpactParam - minimal impct parameter, fm
// fMaxImpactParam - maximal impct parameter, fm
// fEnergyCMS - sqrt(S) per nucleon pair, AGeV
//
// fYMin - minimal rapidity of photons
// fYMax - maximal rapidity of photons
// fPtMin - minimal p_t value of gamma, GeV/c
// fPtMax - maximal p_t value of gamma, GeV/c
//-------------------------------------------------------------------------
// comparison with SPS and RHIC data, prediction for LHC can be found in
// arXiv:0811.2634 [nucl-th]
//-------------------------------------------------------------------------
//Begin_Html
/*
The responsible person for this module is
Andreas Morsch.
*/ //End_Html // // /////////////////////////////////////////////////////////////////// #include#include #include "AliConst.h" #include "AliGenEventHeader.h" #include "AliGenPromptPhotons.h" #include "AliLog.h" #include "AliRun.h" ClassImp(AliGenPromptPhotons) TF1* AliGenPromptPhotons::fgDataPt = NULL; TF1* AliGenPromptPhotons::fgWSzA = NULL; TF1* AliGenPromptPhotons::fgWSzB = NULL; TF1* AliGenPromptPhotons::fgTA = NULL; TF1* AliGenPromptPhotons::fgTB = NULL; TF1* AliGenPromptPhotons::fgTAxTB = NULL; TF1* AliGenPromptPhotons::fgTAB = NULL; //_____________________________________________________________________________ AliGenPromptPhotons::AliGenPromptPhotons() :AliGenerator(-1), fAProjectile(0.), fATarget(0.), fEnergyCMS(0.), fMinImpactParam(0.), fMaxImpactParam(0.) { // // Default constructor // SetCutVertexZ(); SetPtRange(); SetYRange(); } //_____________________________________________________________________________ AliGenPromptPhotons::AliGenPromptPhotons(Int_t npart) :AliGenerator(npart), fAProjectile(208), fATarget(208), fEnergyCMS(5500.), fMinImpactParam(0.), fMaxImpactParam(0.) { // // Standard constructor // fName="PromptPhotons"; fTitle="Prompt photons from pp data fit + T_AB"; SetCutVertexZ(); SetPtRange(); SetYRange(); } //_____________________________________________________________________________ AliGenPromptPhotons::~AliGenPromptPhotons() { // // Standard destructor // delete fgDataPt; delete fgWSzA; delete fgWSzB; delete fgTA; delete fgTB; delete fgTAxTB; delete fgTAB; } //_____________________________________________________________________________ void AliGenPromptPhotons::Init() { // Initialisation fgDataPt = new TF1("fgDataPt",FitData,fPtMin,fPtMax,1); fgDataPt->SetParameter(0,fEnergyCMS); const Double_t d=0.54; // surface thickness (fm) const Double_t ra = 1.19*TMath::Power(fAProjectile,1./3.)-1.61/TMath::Power(fAProjectile,1./3.); const Double_t eps=0.01; // cut WS at ramax: WS(ramax)/WS(0)=eps const Double_t ramax=ra+d*TMath::Log((1.-eps+TMath::Exp(-ra/d))/eps); TF1 fWSforNormA("fWSforNormA",&WSforNorm,0,ramax,2); fWSforNormA.SetParameters(ra,d); fWSforNormA.SetParNames("RA","d"); const Double_t ca=1./fWSforNormA.Integral(0.,ramax); const Double_t rb=1.19*TMath::Power(fATarget,1./3.)-1.61/TMath::Power(fATarget,1./3.); const Double_t rbmax=rb+d*TMath::Log((1.-eps+TMath::Exp(-rb/d))/eps); TF1 fWSforNormB("fWSforNormB",&WSforNorm,0,rbmax,2); fWSforNormB.SetParameters(rb,d); fWSforNormB.SetParNames("RB","d"); const Double_t cb=1./fWSforNormB.Integral(0.,rbmax); fgWSzA = new TF1("fgWSzA",WSz,0.,ramax,4); fgWSzA->SetParameter(0,ra); fgWSzA->SetParameter(1,d); fgWSzA->SetParameter(2,ca); fgTA = new TF1("fgTA",TA,0.,ramax,1); fgTA->SetParameter(0,ramax); fgWSzB = new TF1("fgWSzB",WSz,0.,rbmax,4); fgWSzB->SetParameter(0,rb); fgWSzB->SetParameter(1,d); fgWSzB->SetParameter(2,cb); fgTB = new TF1("fgTB",TB,0.,TMath::Pi(),3); fgTB->SetParameter(0,rbmax); fgTAxTB = new TF1("fgTAxTB",TAxTB,0,ramax,2); fgTAxTB->SetParameter(0,rbmax); fgTAB = new TF1("fgTAB",TAB,0.,ramax+rbmax,2); fgTAB->SetParameter(0,ramax); fgTAB->SetParameter(1,rbmax); } //_____________________________________________________________________________ void AliGenPromptPhotons::Generate() { // // Generate thermal photons of a event // Float_t polar[3]= {0,0,0}; Float_t origin[3]; Float_t time; Float_t p[3]; Float_t random[6]; Int_t nt; for (Int_t j=0;j<3;j++) origin[j]=fOrigin[j]; time = fTimeOrigin; /* if(fVertexSmear==kPerEvent) { Vertex(); for (j=0;j<3;j++) origin[j]=fVertex[j]; time = fTime; } */ TArrayF eventVertex; eventVertex.Set(3); eventVertex[0] = origin[0]; eventVertex[1] = origin[1]; eventVertex[2] = origin[2]; Float_t eventTime = time; Int_t nGam; Float_t b,pt,rapidity,phi,ww; b=TMath::Sqrt(fMinImpactParam*fMinImpactParam+(fMaxImpactParam*fMaxImpactParam-fMinImpactParam*fMinImpactParam)*gRandom->Rndm()); Double_t dsdyPP=fgDataPt->Integral(fPtMin,fPtMax); // pb, fm^2 = 1e10 pb ww=dsdyPP*1.0e-10*(fYMax-fYMin)*fAProjectile*fATarget*fgTAB->Eval(b); nGam=Int_t(ww); if(gRandom->Rndm() < (ww-(Float_t)nGam)) nGam++; if(nGam) { for(Int_t i=0; i GetRandom(); Rndm(random,2); rapidity=(fYMax-fYMin)*random[0]+fYMin; phi=2.*TMath::Pi()*random[1]; p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*TMath::SinH(rapidity); if(fVertexSmear==kPerTrack) { Rndm(random,6); for (Int_t j=0;j<3;j++) { origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* TMath::Sqrt(-2*TMath::Log(random[2*j+1])); } Rndm(random,2); time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()* TMath::Cos(2*random[0]*TMath::Pi())* TMath::Sqrt(-2*TMath::Log(random[1])); } PushTrack(fTrackIt,-1,22,p,origin,polar,time,kPPrimary,nt,1.); } } // Header AliGenEventHeader* header = new AliGenEventHeader("PromptPhotons"); // Event Vertex header->SetPrimaryVertex(eventVertex); header->SetInteractionTime(eventTime); header->SetNProduced(fNpart); gAlice->SetGenEventHeader(header); } void AliGenPromptPhotons::SetPtRange(Float_t ptmin, Float_t ptmax) { AliGenerator::SetPtRange(ptmin, ptmax); } void AliGenPromptPhotons::SetYRange(Float_t ymin, Float_t ymax) { AliGenerator::SetYRange(ymin, ymax); } //********************************************************************************** Double_t AliGenPromptPhotons::FitData(const Double_t* x, const Double_t* par) { //--------------------------------------------------- // input: // x[0] - p_t (GeV). // par[0]=sqrt(s_NN) (GeV), // // output: // d^{2}#sigma/(dp_t dy) (pb/GeV) //--------------------------------------------------- // // d^{2}#sigma/(dp_t dy) = (2 pi p_t) Ed^{3}#sigma/d^{3}p // // data presentation: Nucl.Phys.A783:577-582,2007, hep-ex/0609037, fig.3 // F(x_t)=(#sqrt{s})^{5} Ed^{3}#sigma/d^{3}p //--------------------------------------------------- // approximate tabulation of F(x_t) const Int_t nMax=4; const Double_t log10x[nMax]={ -2., -1., -0.6, -0.3}; const Double_t log10F[nMax]={ 19., 13., 9.8, 7.}; const Double_t xT=2.*x[0]/par[0]; const Double_t log10xT=TMath::Log10(xT); Int_t i=0; while(log10xT>log10x[i] && i SetParameter(3,b); return 2.*fgWSzA->Integral(0.,TMath::Sqrt(ramax*ramax-b*b)); } //********************************************************************************** Double_t AliGenPromptPhotons::TB(const Double_t* x, const Double_t* par) { //--------------------------------------------------- // input: // x[0] - phi (rad) // par[0] - RBMAX (fm), max. value of target radius // par[1] - b (fm), impact parameter // par[2] - s (fm) // // output: // nuclear thickness function T_B(phi)=T_B(sqtr(s**2+b**2-2*s*b*cos(phi))) //--------------------------------------------------- const Double_t phi=x[0]; const Double_t rbmax=par[0],b=par[1],s=par[2]; const Double_t w=TMath::Sqrt(s*s+b*b-2.*s*b*TMath::Cos(phi)); fgWSzB->SetParameter(3,w); return 2.*fgWSzB->Integral(0.,TMath::Sqrt(rbmax*rbmax-w*w));; } //********************************************************************************** Double_t AliGenPromptPhotons::TAxTB(const Double_t* x, const Double_t* par) { //--------------------------------------------------- // input: // x[0] - s (fm) // par[0] - RBMAX (fm), max. value of target radius // par[1] - b (fm), impact parameter // // output: // s * TA(s) * 2 * Integral(0,phiMax) TB(phi(s,b)) //--------------------------------------------------- const Double_t s = x[0]; const Double_t rbmax=par[0],b=par[1]; if(s==0.) return 0.; fgTB->SetParameter(1,b); fgTB->SetParameter(2,s); Double_t phiMax; if(b Eval(s)*2.*fgTB->Integral(0.,phiMax); } // --------------------------------------------------------------------------------- Double_t AliGenPromptPhotons::TAB(const Double_t* x, const Double_t* par) { //--------------------------------------------------- // input: // x[0] - b (fm), impact parameter // par[0] - RAMAX (fm), max. value of projectile radius // par[1] - RAMAX (fm), max. value of target radius // // output: // overlap function TAB(b) (1/fm**2) //--------------------------------------------------- const Double_t b=x[0]; const Double_t ramax=par[0],rbmax=par[1]; Double_t sMin=0.,sMax=ramax; if(b>rbmax) sMin=b-rbmax; if(b<(ramax-rbmax)) sMax=b+rbmax; fgTAxTB->SetParameter(1,b); return fgTAxTB->Integral(sMin,sMax); }