ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eeToHadronsModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eeToHadronsModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 //
32 // File name: G4eeToHadronsModel
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 12.08.2003
37 //
38 // Modifications:
39 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
40 // 18-05-05 Use optimized interfaces (V.Ivantchenko)
41 //
42 //
43 // -------------------------------------------------------------------
44 //
45 
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
50 #include "G4eeToHadronsModel.hh"
51 #include "Randomize.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4Electron.hh"
55 #include "G4Gamma.hh"
56 #include "G4Positron.hh"
57 #include "G4PionPlus.hh"
58 #include "Randomize.hh"
59 #include "G4Vee2hadrons.hh"
60 #include "G4PhysicsVector.hh"
61 #include "G4PhysicsLogVector.hh"
62 #include "G4Log.hh"
63 #include "G4Exp.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
67 using namespace std;
68 
69 G4eeToHadronsModel::G4eeToHadronsModel(G4Vee2hadrons* mod, G4int ver,
70  const G4String& nam)
71  : G4VEmModel(nam),
72  model(mod),
73  crossPerElectron(0),
74  crossBornPerElectron(0),
75  isInitialised(false),
76  nbins(100),
77  verbose(ver)
78 {
79  theGamma = G4Gamma::Gamma();
80  highKinEnergy = HighEnergyLimit();
81  lowKinEnergy = LowEnergyLimit();
82  emin = lowKinEnergy;
83  emax = highKinEnergy;
84  peakKinEnergy = highKinEnergy;
85  epeak = emax;
86  //verbose = 1;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
91 G4eeToHadronsModel::~G4eeToHadronsModel()
92 {
93  delete model;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
98 void G4eeToHadronsModel::Initialise(const G4ParticleDefinition*,
99  const G4DataVector&)
100 {
101  if(isInitialised) { return; }
102  isInitialised = true;
103 
104  // CM system
105  emin = model->LowEnergy();
106  emax = model->HighEnergy();
107 
108  // peak energy
109  epeak = std::min(model->PeakEnergy(), emax);
110 
111  if(verbose>0) {
112  G4cout << "G4eeToHadronsModel::Initialise: " << G4endl;
113  G4cout << "CM System: emin(MeV)= " << emin/MeV
114  << " epeak(MeV)= " << epeak/MeV
115  << " emax(MeV)= " << emax/MeV
116  << G4endl;
117  }
118 
119  crossBornPerElectron = model->PhysicsVector();
120  crossPerElectron = model->PhysicsVector();
121  nbins = crossPerElectron->GetVectorLength();
122  for(G4int i=0; i<nbins; ++i) {
123  G4double e = crossPerElectron->Energy(i);
124  G4double cs = model->ComputeCrossSection(e);
125  crossBornPerElectron->PutValue(i, cs);
126  }
127  ComputeCMCrossSectionPerElectron();
128 
129  if(verbose>1) {
130  G4cout << "G4eeToHadronsModel: Cross sections per electron"
131  << " nbins= " << nbins
132  << " emin(MeV)= " << emin/MeV
133  << " emax(MeV)= " << emax/MeV
134  << G4endl;
135  for(G4int i=0; i<nbins; ++i) {
136  G4double e = crossPerElectron->Energy(i);
137  G4double s1 = crossPerElectron->Value(e);
138  G4double s2 = crossBornPerElectron->Value(e);
139  G4cout << "E(MeV)= " << e/MeV
140  << " cross(nb)= " << s1/nanobarn
141  << " crossBorn(nb)= " << s2/nanobarn
142  << G4endl;
143  }
144  }
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148 
149 G4double G4eeToHadronsModel::CrossSectionPerVolume(
150  const G4Material* mat,
151  const G4ParticleDefinition* p,
152  G4double kineticEnergy,
153  G4double, G4double)
154 {
155  return mat->GetElectronDensity()*
156  ComputeCrossSectionPerElectron(p, kineticEnergy);
157 }
158 
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160 
161 G4double G4eeToHadronsModel::ComputeCrossSectionPerAtom(
162  const G4ParticleDefinition* p,
163  G4double kineticEnergy,
164  G4double Z, G4double,
165  G4double, G4double)
166 {
167  return Z*ComputeCrossSectionPerElectron(p, kineticEnergy);
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 
172 G4double G4eeToHadronsModel::ComputeCrossSectionPerElectron(
173  const G4ParticleDefinition*,
174  G4double energy,
175  G4double, G4double)
176 {
177  return (crossPerElectron) ? crossPerElectron->Value(energy) : 0.0;
178 }
179 
180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
181 
182 void G4eeToHadronsModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
183  const G4MaterialCutsCouple*,
184  const G4DynamicParticle* dParticle,
185  G4double,
186  G4double)
187 {
188  if(crossPerElectron) {
189  G4double t = dParticle->GetKineticEnergy() + 2*electron_mass_c2;
190  G4LorentzVector inlv = dParticle->Get4Momentum() +
191  G4LorentzVector(0.0,0.0,0.0,electron_mass_c2);
192  G4double e = inlv.m();
193  G4ThreeVector inBoost = inlv.boostVector();
194  //G4cout << "G4eeToHadronsModel::SampleSecondaries e= " << e
195  // << " " << inlv << " " << inBoost <<G4endl;
196  if(e > emin) {
197  G4DynamicParticle* gamma = GenerateCMPhoton(e);
198  G4LorentzVector gLv = gamma->Get4Momentum();
199  G4LorentzVector lv(0.0,0.0,0.0,e);
200  lv -= gLv;
201  G4double mass = lv.m();
202  //G4cout << "mass= " << mass << " " << lv << G4endl;
203  G4ThreeVector boost = lv.boostVector();
204  //G4cout << "mass= " << mass << " " << boost << G4endl;
205  const G4ThreeVector dir = gamma->GetMomentumDirection();
206  model->SampleSecondaries(newp, mass, dir);
207  G4int np = newp->size();
208  for(G4int j=0; j<np; ++j) {
209  G4DynamicParticle* dp = (*newp)[j];
210  G4LorentzVector v = dp->Get4Momentum();
211  v.boost(boost);
212  //G4cout << j << ". " << v << G4endl;
213  v.boost(inBoost);
214  //G4cout << " " << v << G4endl;
215  dp->Set4Momentum(v);
216  t -= v.e();
217  }
218  //G4cout << "Gamma " << gLv << G4endl;
219  gLv.boost(inBoost);
220  //G4cout << " " << gLv << G4endl;
221  gamma->Set4Momentum(gLv);
222  t -= gLv.e();
223  newp->push_back(gamma);
224  if(std::abs(t) > MeV) {
225  G4cout << "G4eeToHadronsModel::SampleSecondaries: Ebalance(MeV)= "
226  << t/MeV << " primary 4-momentum: " << inlv << G4endl;
227  }
228  }
229  }
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
234 void G4eeToHadronsModel::ComputeCMCrossSectionPerElectron()
235 {
236  for(G4int i=0; i<nbins; i++) {
237  G4double e = crossPerElectron->Energy(i);
238  G4double cs = 0.0;
239  if(i > 0) {
240  G4double LL = 2.0*G4Log(e/electron_mass_c2);
241  G4double bt = 2.0*fine_structure_const*(LL - 1.0)/pi;
242  G4double btm1= bt - 1.0;
243  G4double del = 1. + fine_structure_const*(1.5*LL + pi*pi/3. -2.)/pi;
244  G4double s1 = crossBornPerElectron->Value(e);
245  G4double e1 = crossPerElectron->Energy(i-1);
246  G4double x1 = 1. - e1/e;
247  cs += s1*(del*G4Exp(G4Log(x1)*bt) - bt*(x1 - 0.25*x1*x1));
248  if(i > 1) {
249  G4double e2 = e1;
250  G4double x2 = x1;
251  G4double s2 = crossBornPerElectron->Value(e2);
252  G4double w2 = bt*(del*G4Exp(G4Log(x2)*btm1) - 1.0 + 0.5*x2);
253  G4double w1;
254 
255  for(G4int j=i-2; j>=0; --j) {
256  e1 = crossPerElectron->Energy(j);
257  x1 = 1. - e1/e;
258  s1 = crossBornPerElectron->Value(e1);
259  w1 = bt*(del*G4Exp(G4Log(x1)*btm1) - 1.0 + 0.5*x1);
260  cs += 0.5*(x1 - x2)*(w2*s2 + w1*s1);
261  e2 = e1;
262  x2 = x1;
263  s2 = s1;
264  w2 = w1;
265  }
266  }
267  }
268  crossPerElectron->PutValue(i, cs);
269  }
270 }
271 
272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
273 
274 G4DynamicParticle* G4eeToHadronsModel::GenerateCMPhoton(G4double e)
275 {
276  G4double x;
277  G4DynamicParticle* gamma = nullptr;
278  G4double LL = 2.0*G4Log(e/electron_mass_c2);
279  G4double bt = 2.0*fine_structure_const*(LL - 1.)/pi;
280  G4double btm1= bt - 1.0;
281  G4double del = 1. + fine_structure_const*(1.5*LL + pi*pi/3. -2.)/pi;
282 
283  G4double s0 = crossBornPerElectron->Value(e);
284  G4double de = (emax - emin)/(G4double)nbins;
285  G4double xmax = 0.5*(1.0 - (emin*emin)/(e*e));
286  G4double xmin = std::min(de/e, xmax);
287  G4double ds = s0*(del*G4Exp(G4Log(xmin)*bt) - bt*(xmin - 0.25*xmin*xmin));
288  G4double e1 = e*(1. - xmin);
289 
290  //G4cout << "e1= " << e1 << G4endl;
291  if(e1 < emax && s0*G4UniformRand()<ds) {
292  x = xmin*G4Exp(G4Log(G4UniformRand())/bt);
293  } else {
294 
295  x = xmin;
296  G4double s1 = crossBornPerElectron->Value(e1);
297  G4double w1 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
298  G4double grej = s1*w1;
299  G4double f;
300  /*
301  G4cout << "e(GeV)= " << e/GeV << " epeak(GeV)= " << epeak/GeV
302  << " s1= " << s1 << " w1= " << w1
303  << " grej= " << grej << G4endl;
304  */
305  // Above emax cross section is const
306  if(e1 > emax) {
307  x = 0.5*(1. - (emax*emax)/(e*e));
308  G4double s2 = crossBornPerElectron->Value(emax);
309  G4double w2 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
310  grej = s2*w2;
311  //G4cout << "emax= " << emax << " s2= " << s2 << " w2= " << w2
312  // << " grej= " << grej << G4endl;
313  }
314 
315  if(e1 > epeak) {
316  x = 0.5*(1.0 - (epeak*epeak)/(e*e));
317  G4double s2 = crossBornPerElectron->Value(epeak);
318  G4double w2 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
319  grej = std::max(grej,s2*w2);
320  //G4cout << "epeak= " << epeak << " s2= " << s2 << " w2= " << w2
321  // << " grej= " << grej << G4endl;
322  }
323  G4int ii = 0;
324  const G4int iimax = 1000;
325  do {
326  x = xmin + G4UniformRand()*(xmax - xmin);
327 
328  G4double s2 = crossBornPerElectron->Value(sqrt(1.0 - 2*x)*e);
329  G4double w2 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
330  /*
331  G4cout << "x= " << x << " xmin= " << xmin << " xmax= " << xmax
332  << " s2= " << s2 << " w2= " << w2 << G4endl;
333  */
334  f = s2*w2;
335  if(f > grej) {
336  G4cout << "G4DynamicParticle* G4eeToHadronsModel:WARNING "
337  << f << " > " << grej << " majorant is`small!"
338  << G4endl;
339  }
340  if(++ii >= iimax) { break; }
341  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
342  } while (f < grej*G4UniformRand());
343  }
344 
345  G4ThreeVector dir(0.0,0.0,1.0);
346  if(G4UniformRand() > 0.5) { dir.set(0.0,0.0,-1.0); }
347  //G4cout << "Egamma(MeV)= " << x*e << " " << dir << G4endl;
348  gamma = new G4DynamicParticle(theGamma,dir,x*e);
349  return gamma;
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353