1 /**************************************************************************
2 * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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. *
16 * Copyright(c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken *
17 * See $ALICE_ROOT/EpEmGen/diffcross.f for full Copyright notice *
20 * Copyright(c) 2002 Kai Hencken, Yuri Kharlov, Serguei Sadovsky *
21 * See $ALICE_ROOT/EpEmGen/epemgen.f for full Copyright notice *
23 **************************************************************************/
27 // Event generator of single e+e- pair production in ultraperipheral PbPb collisions
28 // at 5.5 TeV/nucleon.
29 // The generator is based on 5-dimentional differential cross section of the process.
32 // [1] "Multiple electromagnetic electron positron pair production in
33 // relativistic heavy ion collisions".
34 // Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur,
35 // Phys. Rev. A55 (1997) 396.
36 // [2] K.Hencken, Yu.Kharlov, S.Sadovsky, Internal ALICE Note 2002-27.
40 // AliGenEpEmv1 *gener = new AliGenEpEmv1();
41 // gener->SetXXXRange(); // Set kinematics range
44 // gener->Generate(); // Produce one e+e- pair with the event weight assigned
45 // // to each track. The sum of event weights, divided by
46 // // the total number of generated events, gives the
47 // // integral cross section of the process of e+e- pair
48 // // production in the above mentioned kinematics range.
49 // // Sum of the selected event weights, divided by the total
50 // // number of generated events, gives the integral cross
51 // // section corresponded to the set of selected events
53 // The generator consists of several modules:
54 // 1) $ALICE_ROOT/EpEmGen/diffcross.f:
55 // Exact calculation of the total differential e+ e- -pair production
56 // in Relativistic Heavy Ion Collisions for a point particle in an
57 // external field approach. See full comments in the mentioned file.
58 // 2) $ALICE_ROOT/EpEmGen/epemgen.f:
59 // Generator of e+e- pairs produced in PbPb collisions at LHC
60 // it generates events according to the parametrization of the
61 // differential cross section. Produces events have weights calculated
62 // by the exact differential cross section calculation (diffcross.f).
63 // See full comments in the mentioned file.
65 // Interface from the fortran event generator to ALIROOT
66 // 4) Class AliGenEpEmv1:
67 // The event generator to call within ALIROOT
69 // Author of this module: Yuri.Kharlov@cern.ch
72 #include "AliGenEpEmv1.h"
73 #include <TParticle.h>
74 #include <TParticlePDG.h>
75 #include <TDatabasePDG.h>
78 ClassImp(AliGenEpEmv1)
80 //------------------------------------------------------------
82 AliGenEpEmv1::AliGenEpEmv1()
84 // Default constructor
86 if (fPtMin == 0) fPtMin = 1.E-04;
89 //____________________________________________________________
90 AliGenEpEmv1::AliGenEpEmv1(const AliGenEpEmv1 & gen)
96 //____________________________________________________________
97 AliGenEpEmv1::~AliGenEpEmv1()
102 //____________________________________________________________
103 void AliGenEpEmv1::Init()
106 // 1) define a generator
107 // 2) initialize the generator of e+e- pair production
109 fMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
111 SetMC(new TEpEmGen());
112 fEpEmGen = (TEpEmGen*) fMCEvGen;
113 fEpEmGen ->Initialize(fYMin,fYMax,fPtMin,fPtMax);
117 //____________________________________________________________
118 void AliGenEpEmv1::Generate()
121 // Generate one e+e- pair
122 // Gaussian smearing on the vertex is done if selected.
124 // Each produced e+e- pair is defined by the following variables:
125 // rapidities of e-, e+ (yElectron,yPositron)
126 // log10(pt in MeV/c) of e-, e+ (xElectron,xPositron)
127 // azymuth angles between e- and e+ (phi12)
129 // On output an event weight is given (weight) which is assigned to each track.
130 // The sum of event weights, divided by the total number of generated events,
131 // gives the integral cross section of the e+e- pair production in the
132 // selected kinematics range.
135 Float_t polar[3]= {0,0,0};
139 Double_t ptElectron,ptPositron, phiElectron,phiPositron, mt;
140 Double_t phi12=0,xElectron=0,xPositron=0,yElectron=0,yPositron=0,weight=0;
144 fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,
145 yElectron,yPositron,xElectron,xPositron,phi12,weight);
147 printf("AliGenEpEmv1::Generate(): y=(%f,%f), x=(%f,%f), phi=%f\n",
148 yElectron,yPositron,xElectron,xPositron,phi12);
150 for (j=0;j<3;j++) origin[j]=fOrigin[j];
151 if(fVertexSmear==kPerEvent) {
154 origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
155 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
160 ptElectron = TMath::Power(10,xElectron) * 1.e-03;;
161 ptPositron = TMath::Power(10,xPositron) * 1.e-03;;
162 phiElectron = fPhiMin + random[0] * (fPhiMax-fPhiMin);
163 phiPositron = phiElectron + phi12;
166 mt = TMath::Sqrt(ptElectron*ptElectron + fMass*fMass);
167 p[0] = ptElectron*TMath::Cos(phiElectron);
168 p[1] = ptElectron*TMath::Sin(phiElectron);
169 p[2] = mt*TMath::SinH(yElectron);
172 printf("id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",id,p[0],p[1],p[2]);
173 PushTrack(fTrackIt,-1, id,p,origin,polar,0,kPPrimary,nt,weight);
176 mt = TMath::Sqrt(ptPositron*ptPositron + fMass*fMass);
177 p[0] = ptPositron*TMath::Cos(phiPositron);
178 p[1] = ptPositron*TMath::Sin(phiPositron);
179 p[2] = mt*TMath::SinH(yPositron);
182 printf("id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",id,p[0],p[1],p[2]);
183 PushTrack(fTrackIt,-1, id,p,origin,polar,0,kPPrimary,nt,weight);
186 if (fEvent%1000 == 0) {
187 printf("=====> AliGenEpEmv1::Generate(): \n Event %d, sigma=%f +- %f kb\n",
188 fEvent,fEpEmGen->GetXsection(),fEpEmGen->GetDsection());