Transition to NewIO
[u/mrichter/AliRoot.git] / TEPEMGEN / AliGenEpEmv1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2002, 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  * Copyright(c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken          *
17  * See $ALICE_ROOT/EpEmGen/diffcross.f for full Copyright notice          *
18  *                                                                        *
19  *                                                                        *
20  * Copyright(c) 2002 Kai Hencken, Yuri Kharlov, Serguei Sadovsky          *
21  * See $ALICE_ROOT/EpEmGen/epemgen.f for full Copyright notice            *
22  *                                                                        *
23  **************************************************************************/
24
25 /*
26 $Log$
27 Revision 1.1  2002/11/06 10:26:45  hristov
28 Event generator for e+e- pair production (Yu.Kharlov)
29
30 */
31
32 // Event generator of single e+e- pair production in ultraperipheral PbPb collisions
33 // at 5.5 TeV/nucleon.
34 // The generator is based on 5-dimentional differential cross section of the process.
35 //%
36 // References:
37 // [1] "Multiple electromagnetic electron positron pair production in
38 //      relativistic heavy ion collisions".
39 //      Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur,
40 //      Phys. Rev. A55 (1997) 396.
41 // [2] K.Hencken, Yu.Kharlov, S.Sadovsky, Internal ALICE Note 2002-27.
42 //%
43 // Usage:
44 // Initialization:
45 //    AliGenEpEmv1 *gener = new AliGenEpEmv1();
46 //    gener->SetXXXRange(); // Set kinematics range
47 //    gener->Init();
48 // Event generation:
49 //    gener->Generate(); // Produce one e+e- pair with the event weight assigned 
50 //                       // to each track. The sum of event weights, divided by 
51 //                       // the total number of generated events, gives the 
52 //                       // integral cross section of the process of e+e- pair 
53 //                       // production in the above mentioned kinematics range.
54 //                       // Sum of the selected event weights, divided by the total 
55 //                       // number of generated events, gives the integral cross 
56 //                       // section corresponded to the set of selected events
57 //%
58 // The generator consists of several modules:
59 // 1) $ALICE_ROOT/EpEmGen/diffcross.f:
60 //    Exact calculation of the total differential e+ e- -pair production
61 //    in Relativistic Heavy Ion Collisions for a point particle in an
62 //    external field approach. See full comments in the mentioned file.
63 // 2) $ALICE_ROOT/EpEmGen/epemgen.f:
64 //    Generator of e+e- pairs produced in PbPb collisions at LHC
65 //    it generates events according to the parametrization of the
66 //    differential cross section. Produces events have weights calculated
67 //    by the exact differential cross section calculation (diffcross.f).
68 //    See full comments in the mentioned file.
69 // 3) Class TEpEmGen:
70 //    Interface from the fortran event generator to ALIROOT
71 // 4) Class AliGenEpEmv1:
72 //    The event generator to call within ALIROOT
73 //%
74 // Author of this module: Yuri.Kharlov@cern.ch
75 // 9 October 2002
76
77 #include "AliGenEpEmv1.h"
78 #include <TParticle.h>
79 #include <TParticlePDG.h>
80 #include <TDatabasePDG.h>
81 #include <TEpEmGen.h>
82
83 ClassImp(AliGenEpEmv1)
84
85 //------------------------------------------------------------
86
87 AliGenEpEmv1::AliGenEpEmv1()
88 {
89   // Default constructor
90   // Avoid zero pt
91   if (fPtMin == 0) fPtMin = 1.E-04;
92 }
93
94 //____________________________________________________________
95 AliGenEpEmv1::AliGenEpEmv1(const AliGenEpEmv1 & gen)
96 {
97   // copy constructor
98   gen.Copy(*this);
99 }
100
101 //____________________________________________________________
102 AliGenEpEmv1::~AliGenEpEmv1()
103 {
104   // Destructor
105 }
106
107 //____________________________________________________________
108 void AliGenEpEmv1::Init()
109 {
110   // Initialisation:
111   // 1) define a generator
112   // 2) initialize the generator of e+e- pair production
113
114   fMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
115
116   SetMC(new TEpEmGen());
117   fEpEmGen = (TEpEmGen*) fgMCEvGen;
118   fEpEmGen ->Initialize(fYMin,fYMax,fPtMin,fPtMax);
119   fEvent = 0;
120 }
121
122 //____________________________________________________________
123 void AliGenEpEmv1::Generate()
124 {
125   //
126   // Generate one e+e- pair
127   // Gaussian smearing on the vertex is done if selected. 
128   //%
129   // Each produced e+e- pair is defined by the following variables:
130   // rapidities of e-, e+ (yElectron,yPositron)
131   // log10(pt in MeV/c) of e-, e+ (xElectron,xPositron)
132   // azymuth angles between e- and e+ (phi12)
133   //%
134   // On output an event weight is given (weight) which is assigned to each track.
135   // The sum of event weights, divided by the total number of generated events, 
136   // gives the integral cross section of the e+e- pair production in the   
137   // selected kinematics range.   
138   //
139
140   Float_t polar[3]= {0,0,0};
141   Float_t origin[3];
142   Float_t p[3];
143
144   Double_t ptElectron,ptPositron, phiElectron,phiPositron, mt;
145   Double_t phi12=0,xElectron=0,xPositron=0,yElectron=0,yPositron=0,weight=0;
146   Int_t   j, nt, id;
147   Float_t random[6];
148
149   fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,
150            yElectron,yPositron,xElectron,xPositron,phi12,weight);
151   if (fDebug == 1)
152     printf("AliGenEpEmv1::Generate(): y=(%f,%f), x=(%f,%f), phi=%f\n",
153            yElectron,yPositron,xElectron,xPositron,phi12);
154
155   for (j=0;j<3;j++) origin[j]=fOrigin[j];
156   if(fVertexSmear==kPerEvent) {
157     Rndm(random,6);
158     for (j=0;j<3;j++) {
159       origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
160         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
161     }
162   }
163
164   Rndm(random,1);
165   ptElectron  = TMath::Power(10,xElectron) * 1.e-03;;
166   ptPositron  = TMath::Power(10,xPositron) * 1.e-03;;
167   phiElectron = fPhiMin + random[0] * (fPhiMax-fPhiMin);
168   phiPositron = phiElectron + phi12;
169
170   // Produce electron
171   mt = TMath::Sqrt(ptElectron*ptElectron + fMass*fMass);
172   p[0] = ptElectron*TMath::Cos(phiElectron);
173   p[1] = ptElectron*TMath::Sin(phiElectron);
174   p[2] = mt*TMath::SinH(yElectron);
175   id =  11;
176   if (fDebug == 2)
177     printf("id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",id,p[0],p[1],p[2]);
178   SetTrack(fTrackIt,-1, id,p,origin,polar,0,kPPrimary,nt,weight);
179
180   // Produce positron
181   mt = TMath::Sqrt(ptPositron*ptPositron + fMass*fMass);
182   p[0] = ptPositron*TMath::Cos(phiPositron);
183   p[1] = ptPositron*TMath::Sin(phiPositron);
184   p[2] = mt*TMath::SinH(yPositron);
185   id = -11;
186   if (fDebug == 2)
187     printf("id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",id,p[0],p[1],p[2]);
188   SetTrack(fTrackIt,-1, id,p,origin,polar,0,kPPrimary,nt,weight);
189   
190   fEvent++;
191   if (fEvent%1000 == 0) {
192     printf("=====> AliGenEpEmv1::Generate(): \n   Event %d, sigma=%f +- %f kb\n",
193            fEvent,fEpEmGen->GetXsection(),fEpEmGen->GetDsection());
194   }
195 }
196