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 /* $Id$ */
  26
  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.
  30 //%
  31 // References:
  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.
  37 //%
  38 // Usage:
  39 // Initialization:
  40 //    AliGenEpEmv1 *gener = new AliGenEpEmv1();
  41 //    gener->SetXXXRange(); // Set kinematics range
  42 //    gener->Init();
  43 // Event generation:
  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
  52 //%
  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.
  64 // 3) Class TEpEmGen:
  65 //    Interface from the fortran event generator to ALIROOT
  66 // 4) Class AliGenEpEmv1:
  67 //    The event generator to call within ALIROOT
  68 //%
  69 // Author of this module: Yuri.Kharlov@cern.ch
  70 // 9 October 2002
  71
  72 #include "AliGenEpEmv1.h"
  73 #include <TParticle.h>
  74 #include <TParticlePDG.h>
  75 #include <TDatabasePDG.h>
  76 #include <TEpEmGen.h>
  77
  78 ClassImp(AliGenEpEmv1)
  79
  80 //------------------------------------------------------------
  81
  82 AliGenEpEmv1::AliGenEpEmv1():
  83   fMass(0),
  84   fEpEmGen(0),
  85   fDebug(0),
  86   fEvent(0),
  87   fHeader(AliGenEventHeader())
  88 {
  89   // Default constructor
  90   // Avoid zero pt
  91   if (fPtMin == 0) fPtMin = 1.E-04;
  92 }
  93
  94 //____________________________________________________________
  95 AliGenEpEmv1::~AliGenEpEmv1()
  96 {
  97   // Destructor
  98 }
  99
 100 //____________________________________________________________
 101 void AliGenEpEmv1::Init()
 102 {
 103   // Initialisation:
 104   // 1) define a generator
 105   // 2) initialize the generator of e+e- pair production
 106
 107   fMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
 108
 109   SetMC(new TEpEmGen());
 110   fEpEmGen = (TEpEmGen*) fMCEvGen;
 111   fEpEmGen ->Initialize(fYMin,fYMax,fPtMin,fPtMax);
 112   fEvent = 0;
 113 }
 114
 115 //____________________________________________________________
 116 void AliGenEpEmv1::Generate()
 117 {
 118   //
 119   // Generate one e+e- pair
 120   // Gaussian smearing on the vertex is done if selected.
 121   //%
 122   // Each produced e+e- pair is defined by the following variables:
 123   // rapidities of e-, e+ (yElectron,yPositron)
 124   // log10(pt in MeV/c) of e-, e+ (xElectron,xPositron)
 125   // azymuth angles between e- and e+ (phi12)
 126   //%
 127   // On output an event weight is given (weight) which is assigned to each track.
 128   // The sum of event weights, divided by the total number of generated events,
 129   // gives the integral cross section of the e+e- pair production in the
 130   // selected kinematics range.
 131   //
 132
 133   Float_t polar[3]= {0,0,0};
 134   Float_t origin[3];
 135   Float_t time = 0.;
 136   Float_t p[3];
 137
 138   Double_t ptElectron,ptPositron, phiElectron,phiPositron, mt;
 139   Double_t phi12=0,xElectron=0,xPositron=0,yElectron=0,yPositron=0,weight=0;
 140   Int_t   j, nt, id;
 141   Float_t random[6];
 142
 143   fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,
 144            yElectron,yPositron,xElectron,xPositron,phi12,weight);
 145   if (fDebug == 1)
 146     printf("AliGenEpEmv1::Generate(): y=(%f,%f), x=(%f,%f), phi=%f\n",
 147            yElectron,yPositron,xElectron,xPositron,phi12);
 148
 149   for (j=0;j<3;j++) origin[j]=fOrigin[j];
 150   time = fTimeOrigin;
 151   if(fVertexSmear==kPerEvent) {
 152     Rndm(random,6);
 153     for (j=0;j<3;j++) {
 154       origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 155         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 156     }
 157     Rndm(random,2);
 158     time += fOsigma[2]/TMath::Ccgs()*
 159       TMath::Cos(2*random[0]*TMath::Pi())*
 160       TMath::Sqrt(-2*TMath::Log(random[1]));
 161   }
 162
 163   Rndm(random,1);
 164   ptElectron  = TMath::Power(10,xElectron) * 1.e-03;;
 165   ptPositron  = TMath::Power(10,xPositron) * 1.e-03;;
 166   phiElectron = fPhiMin + random[0] * (fPhiMax-fPhiMin);
 167   phiPositron = phiElectron + phi12;
 168
 169   // Produce electron
 170   mt = TMath::Sqrt(ptElectron*ptElectron + fMass*fMass);
 171   p[0] = ptElectron*TMath::Cos(phiElectron);
 172   p[1] = ptElectron*TMath::Sin(phiElectron);
 173   p[2] = mt*TMath::SinH(yElectron);
 174   id =  11;
 175   if (fDebug == 2)
 176     printf("id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",id,p[0],p[1],p[2]);
 177   PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,weight);
 178
 179   // Produce positron
 180   mt = TMath::Sqrt(ptPositron*ptPositron + fMass*fMass);
 181   p[0] = ptPositron*TMath::Cos(phiPositron);
 182   p[1] = ptPositron*TMath::Sin(phiPositron);
 183   p[2] = mt*TMath::SinH(yPositron);
 184   id = -11;
 185   if (fDebug == 2)
 186     printf("id=%+3d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",id,p[0],p[1],p[2]);
 187   PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,weight);
 188
 189   fEvent++;
 190   if (fEvent%1000 == 0) {
 191     printf("=====> AliGenEpEmv1::Generate(): \n   Event %d, sigma=%f +- %f kb\n",
 192            fEvent,fEpEmGen->GetXsection(),fEpEmGen->GetDsection());
 193   }
 194   fHeader.SetEventWeight(weight);
 195   fHeader.SetInteractionTime(time);
 196   AddHeader(&fHeader);
 197 }
 198