Merge branch 'feature-doxygen'
[u/mrichter/AliRoot.git] / TEPEMGEN / AliGenQEDBg.cxx
CommitLineData
0ebd3869 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 for background from e+e- pair production in ultraperipheral PbPb collisions
28// at 5.5 TeV/nucleon, integrated over specific readout cycle
29// Derived from AliGenEpEmv1
30// Author: ruben.shahoyan@cern.ch
31//%
32// References:
33// [1] "Multiple electromagnetic electron positron pair production in
34// relativistic heavy ion collisions".
35// Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur,
36// Phys. Rev. A55 (1997) 396.
37// [2] K.Hencken, Yu.Kharlov, S.Sadovsky, Internal ALICE Note 2002-27.
38//%
39// Usage:
40// Initialization:
41// AliGenQEDBg *gener = new AliGenQEDBg();
42// gener->SetXXXRange(); // Set kinematics range
43// gener->SetLumiIntTime(double lumi, double sec); // luminosity and intergration time in seconds
44// gener->Init();
45// Event generation:
46// gener->Generate(); // Produce poissonian number of e+e- pair with average number
47// corresponding to requested integration time with given beam luminosity
48// Each pair has its own vertex, and time randomly distributed at originT
49// and originT+ integration time
50//
15a3d47f 51// For details of pair generation see AliGenEpEmv1.cxx by Yuri.Kharlov@cern.ch
52
53/*
54 The most useful way of using it is in the overlay with other generators (using the cocktail):
55 In the Config.C
56 // load the library:
4070f709 57 gSystem->Load("libTEPEMGEN");
15a3d47f 58 //
59 // add to AliGenCocktail as:
60 AliGenCocktail *cocktail = new AliGenCocktail();
61 // ... setup other stuff
62 //
63 // QED background
64 AliGenQEDBg* genBg = new AliGenQEDBg();
65 genBg->SetEnergyCMS(5500);
66 genBg->SetProjectile("A", 208, 82);
67 genBg->SetTarget ("A", 208, 82);
68 genBg->SetYRange(-6.,3);
69 genBg->SetPtRange(1.e-3,1.0); // Set pt limits (GeV) for e+-: 1MeV corresponds to max R=13.3mm at 5kGaus
70 genBg->SetLumiIntTime(6.e27,20e-6); // luminosity and integration time
71 //
72 cocktail->AddGenerator(genBg,"QEDep",1);
73 // We need to generate independent vertex for every event, this must be set after adding to cocktail
74 // Note that the IO origin and sigma's is transferred from AliGenCocktail to daughter generators
75 genBg->SetVertexSource(kInternal);
76*/
0ebd3869 77
78#include "AliLog.h"
79#include "AliGenQEDBg.h"
80#include <TParticle.h>
81#include <TParticlePDG.h>
82#include <TDatabasePDG.h>
83#include <TEpEmGen.h>
84
85ClassImp(AliGenQEDBg)
86
87//------------------------------------------------------------
88
89AliGenQEDBg::AliGenQEDBg()
90: fLumi(0)
91 ,fXSection(0)
92 ,fXSectionEps(1e-2)
93 ,fIntTime(0)
94 ,fPairsInt(-1)
95 ,fMinXSTest(1e3)
96 ,fMaxXSTest(1e7)
97{
98}
99
100//____________________________________________________________
101AliGenQEDBg::~AliGenQEDBg()
102{
103 // Destructor
104}
105
106//____________________________________________________________
107void AliGenQEDBg::Init()
108{
109 // Initialisation:
110 AliInfo(Form("Will estimate QED bg. for L=%e cm^-2*s^-1 and Integration Time of %e s.",fLumi,fIntTime));
111 if (fLumi<=0 || fIntTime<=0) {
112 AliWarning("One of parameters is not set properly, no pairs will be generated");
113 return;
114 }
115 //
116 // initialize the generator of e+e- pair production
117 AliGenEpEmv1::Init();
118 //
119 fPairsInt = 0;
120 int ngen = 0;
121 AliInfo(Form("Estimating x-section with min.relative precision of %f and min/max test: %d/%d",
122 fXSectionEps,int(fMinXSTest),int(fMaxXSTest)));
123 //
bdc5420c 124 double yElectron,yPositron,xElectron,xPositron,phi12,weight,err=0;
0ebd3869 125 fXSection = -1;
126 do {
127 fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
128 if (++ngen>fMinXSTest) { // ensure min number of tests
129 fXSection = fEpEmGen->GetXsection();
130 err = fEpEmGen->GetDsection();
131 }
132 } while(!((fXSection>0 && err/fXSection<fXSectionEps) || ngen>fMaxXSTest));
133 //
134 if (fXSection<=0) {
135 AliError(Form("X-section = %e after %d trials, cannot generate",fXSection,ngen));
136 return;
137 }
138 fPairsInt = fXSection*1e-21*fLumi*fIntTime; // xsestion is in kbarn!
139 AliInfo(Form("Estimated x-secion: %e+-%ekb in %d tests, <Npairs>=%e per %e time interval",
140 fXSection,err,ngen,fPairsInt,fIntTime));
141 //
142}
143
144//____________________________________________________________
145void AliGenQEDBg::Generate()
146{
147 //
148 // Generate poissian <fPairsInt> e+e- pairs, each one with its vertex
149 //
150 Float_t polar[3]= {0,0,0};
151 Float_t origin[3];
152 Float_t time = 0.;
153 Float_t p[3];
154 //
155 int npairs,nt,id;;
156 if (fPairsInt<=0 || (npairs=gRandom->Poisson(fPairsInt))<1) return;
157 //
158 Double_t ptElectron,ptPositron, phiElectron,phiPositron, mt, ms2 = fMass*fMass;
159 Double_t phi12=0,xElectron=0,xPositron=0,yElectron=0,yPositron=0,weight=0;
160 //
161 for (int i=0;i<npairs;i++) {
162 // each pair has its own vertex and time
163 Vertex();
164 for (int j=0;j<3;j++) origin[j] = fVertex[j];
165 time = fTimeOrigin+gRandom->Rndm()*fIntTime;
166 //
167 fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
168 ptElectron = TMath::Power(10,xElectron) * 1.e-03;;
169 ptPositron = TMath::Power(10,xPositron) * 1.e-03;;
170 phiElectron = fPhiMin + gRandom->Rndm() * (fPhiMax-fPhiMin);
171 phiPositron = phiElectron + phi12;
172 // Produce electron
173 mt = TMath::Sqrt(ptElectron*ptElectron + ms2);
174 p[0] = ptElectron*TMath::Cos(phiElectron);
175 p[1] = ptElectron*TMath::Sin(phiElectron);
176 p[2] = mt*TMath::SinH(yElectron);
177 id = 11;
178 PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,1);
179 //
180 // Produce positron
181 mt = TMath::Sqrt(ptPositron*ptPositron + ms2);
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 PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,1);
187 //
188 }
189 fEvent++;
190 //
191 fHeader.SetNProduced(2*npairs);
192 fHeader.SetEventWeight(1);
193 fHeader.SetInteractionTime(fTimeOrigin);
194 AddHeader(&fHeader);
195}
196
197//__________________________________________________________
198void AliGenQEDBg::SetLumiIntTime(double lumi, double intTime)
199{
200 // assign luminosity and integration time
201 if (lumi<=0) AliFatal(Form("Luminosity must be positive, in cm^-2*s^-1, %e asked",lumi));
202 if (intTime<=0) AliFatal(Form("Integration time must be positive, in seconnds, %e asked",intTime));
203 fLumi = lumi;
204 fIntTime = intTime;
205 //
206}
207
208//__________________________________________________________
209void AliGenQEDBg::SetMinMaxXSTest(double mn,double mx)
210{
211 // set min,max number of generator calls for xsection estimates
212 fMinXSTest = mn>100 ? mn : 100.;
213 fMaxXSTest = mx>mx ? mx : mx+100.;
214}