| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Author: Sebastien Incerti
27 // 22 January 2012
28 // on base of G4BoldyshevTripletModel
29 // and G4LivermoreRayleighModel (MT ve
30
31 #include "G4BoldyshevTripletModel.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Log.hh"
35 #include "G4Exp.hh"
36
37 //....oooOO0OOooo........oooOO0OOooo........oo
38
39 using namespace std;
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 const G4int G4BoldyshevTripletModel::maxZ;
44 G4PhysicsFreeVector* G4BoldyshevTripletModel::
45
46 G4BoldyshevTripletModel::G4BoldyshevTripletMod
47 :G4VEmModel(nam),smallEnergy(4.*MeV)
48 {
49 fParticleChange = nullptr;
50
51 lowEnergyLimit = 4.0*electron_mass_c2;
52 momentumThreshold_c = energyThreshold = xb =
53
54 verboseLevel= 0;
55 // Verbosity scale for debugging purposes:
56 // 0 = nothing
57 // 1 = calculation of cross sections, file o
58 // 2 = entering in methods
59
60 if(verboseLevel > 0)
61 {
62 G4cout << "G4BoldyshevTripletModel is cons
63 }
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67
68 G4BoldyshevTripletModel::~G4BoldyshevTripletMo
69 {
70 if(IsMaster()) {
71 for(G4int i=0; i<maxZ; ++i) {
72 if(data[i]) {
73 delete data[i];
74 data[i] = nullptr;
75 }
76 }
77 }
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81
82 void G4BoldyshevTripletModel::Initialise(const
83 const G4DataVector&)
84 {
85 if (verboseLevel > 1)
86 {
87 G4cout << "Calling Initialise() of G4Boldy
88 << G4endl
89 << "Energy range: "
90 << LowEnergyLimit() / MeV << " MeV - "
91 << HighEnergyLimit() / GeV << " GeV isMas
92 << G4endl;
93 }
94 // compute values only once
95 energyThreshold = 1.1*electron_mass_c2;
96 momentumThreshold_c = std::sqrt(energyThresh
97 - electron_mas
98 G4double momentumThreshold_N = momentumThres
99 G4double t = 0.5*G4Log(momentumThreshold_N +
100 std::sqrt(momentumThreshold_N*momentumT
101 //G4cout << 0.5*asinh(momentumThreshold_N) <
102 G4double sinht = std::sinh(t);
103 G4double cosht = std::cosh(t);
104 G4double logsinht = G4Log(2.*sinht);
105 G4double J1 = 0.5*(t*cosht/sinht - logsinht)
106 G4double J2 = (-2./3.)*logsinht + t*cosht/si
107 + (sinht - t*cosht*cosht*cosht)/(3.*sinht*
108
109 xb = 2.*(J1-J2)/J1;
110 xn = 1. - xb/6.;
111
112 if(IsMaster())
113 {
114 // Access to elements
115 const char* path = G4FindDataDir("G4LEDATA
116
117 G4ProductionCutsTable* theCoupleTable =
118 G4ProductionCutsTable::GetProductionCuts
119
120 G4int numOfCouples = (G4int)theCoupleTable
121
122 for(G4int i=0; i<numOfCouples; ++i)
123 {
124 const G4Material* material =
125 theCoupleTable->GetMaterialCutsCouple(
126 const G4ElementVector* theElementVector
127 std::size_t nelm = material->GetNumberOf
128
129 for (std::size_t j=0; j<nelm; ++j)
130 {
131 G4int Z = std::min((*theElementVector)
132 if(!data[Z]) { ReadData(Z, path); }
133 }
134 }
135 }
136 if(!fParticleChange) {
137 fParticleChange = GetParticleChangeForGamm
138 }
139 }
140
141 //....oooOO0OOooo........oooOO0OOooo........oo
142
143 G4double
144 G4BoldyshevTripletModel::MinPrimaryEnergy(cons
145 const G4ParticleDefinition*,
146 G4double)
147 {
148 return lowEnergyLimit;
149 }
150
151 //....oooOO0OOooo........oooOO0OOooo........oo
152
153 void G4BoldyshevTripletModel::ReadData(size_t
154 {
155 if (verboseLevel > 1)
156 {
157 G4cout << "Calling ReadData() of G4Boldysh
158 << G4endl;
159 }
160
161 if(data[Z]) { return; }
162
163 const char* datadir = path;
164
165 if(!datadir)
166 {
167 datadir = G4FindDataDir("G4LEDATA");
168 if(!datadir)
169 {
170 G4Exception("G4BoldyshevTripletModel::Re
171 "em0006",FatalException,
172 "Environment variable G4LEDATA not defin
173 return;
174 }
175 }
176
177 data[Z] = new G4PhysicsFreeVector(0,/*spline
178 std::ostringstream ost;
179 ost << datadir << "/livermore/tripdata/pp-tr
180 std::ifstream fin(ost.str().c_str());
181
182 if( !fin.is_open())
183 {
184 G4ExceptionDescription ed;
185 ed << "G4BoldyshevTripletModel data file <
186 << "> is not opened!" << G4endl;
187 G4Exception("G4BoldyshevTripletModel::Read
188 "em0003",FatalException,
189 ed,"G4LEDATA version should be G4EMLOW6.27
190 return;
191 }
192
193 else
194 {
195
196 if(verboseLevel > 3) { G4cout << "File " <
197 << " is opened by G4BoldyshevTripletMod
198
199 data[Z]->Retrieve(fin, true);
200 }
201
202 // Activation of spline interpolation
203 data[Z]->FillSecondDerivatives();
204 }
205
206 //....oooOO0OOooo........oooOO0OOooo........oo
207
208 G4double G4BoldyshevTripletModel::ComputeCross
209 const G4ParticleDefinition* part,
210 G4double GammaEnergy, G4double Z, G4double,
211 {
212 if (verboseLevel > 1)
213 {
214 G4cout << "Calling ComputeCrossSectionPerA
215 << G4endl;
216 }
217
218 if (GammaEnergy < lowEnergyLimit) { return 0
219
220 G4double xs = 0.0;
221 G4int intZ = std::max(1, std::min(G4lrint(Z)
222 G4PhysicsFreeVector* pv = data[intZ];
223
224 // if element was not initialised
225 // do initialisation safely for MT mode
226 if(!pv)
227 {
228 InitialiseForElement(part, intZ);
229 pv = data[intZ];
230 if(!pv) { return xs; }
231 }
232 // x-section is taken from the table
233 xs = pv->Value(GammaEnergy);
234
235 if(verboseLevel > 1)
236 {
237 G4cout << "*** Triplet conversion xs for
238 << GammaEnergy/MeV << " cs=" << xs/mil
239 }
240 return xs;
241 }
242
243 //....oooOO0OOooo........oooOO0OOooo........oo
244
245 void G4BoldyshevTripletModel::SampleSecondarie
246 std::vector<G
247 const G4MaterialCutsCouple* /*couple*
248 const G4DynamicParticle* aDynamicGamm
249 G4double, G4double)
250 {
251
252 // The energies of the secondary particles a
253 // a modified Wheeler-Lamb model (see PhysRe
254 if (verboseLevel > 1) {
255 G4cout << "Calling SampleSecondaries() of
256 << G4endl;
257 }
258
259 G4double photonEnergy = aDynamicGamma->GetKi
260 G4ParticleMomentum photonDirection = aDynami
261
262 G4double epsilon;
263
264 CLHEP::HepRandomEngine* rndmEngine = G4Rando
265
266 // recoil electron thould be 3d particle
267 G4DynamicParticle* particle3 = nullptr;
268 static const G4double costlim = std::cos(4.4
269
270 G4double loga, f1_re, greject, cost;
271 G4double cosThetaMax = (energyThreshold - el
272 + electron_mass_c2*(energyThreshold + ele
273 /momentumThreshold_c;
274 if (cosThetaMax > 1.) {
275 //G4cout << "G4BoldyshevTripletModel::Samp
276 // << cosThetaMax << G4endl;
277 cosThetaMax = 1.0;
278 }
279
280 G4double logcostm = G4Log(cosThetaMax);
281 do {
282 cost = G4Exp(logcostm*rndmEngine->flat());
283 G4double are = 1./(14.*cost*cost);
284 G4double bre = (1.-5.*cost*cost)/(2.*cost)
285 loga = G4Log((1.+ cost)/(1.- cost));
286 f1_re = 1. - bre*loga;
287 greject = (cost < costlim) ? are*f1_re : 1
288 } while(greject < rndmEngine->flat());
289
290 // Calculo de phi - elecron de recoil
291 G4double sint2 = (1. - cost)*(1. + cost);
292 G4double fp = 1. - sint2*loga/(2.*cost) ;
293 G4double rt, phi_re;
294 do {
295 phi_re = twopi*rndmEngine->flat();
296 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/t
297 } while(rt < rndmEngine->flat());
298
299 // Calculo de la energia - elecron de recoil
300 G4double S = electron_mass_c2*(2.* photonEn
301 G4double P2 = S - electron_mass_c2*electron_
302
303 G4double D2 = 4.*S * electron_mass_c2*electr
304 G4double ener_re = electron_mass_c2 * (S + e
305
306 if(ener_re >= energyThreshold)
307 {
308 G4double electronRKineEnergy = ener_re -
309 G4double sint = std::sqrt(sint2);
310 G4ThreeVector electronRDirection (sint*s
311 electronRDirection.rotateUz(photonDirect
312 particle3 = new G4DynamicParticle (G4Ele
313 electronRDirection,
314 electronRKineEnergy);
315 }
316 else
317 {
318 // deposito la energia ener_re - electr
319 // G4cout << "electron de retroceso " <<
320 fParticleChange->ProposeLocalEnergyDepos
321 ener_re = 0.0;
322 }
323
324 // Depaola (2004) suggested distribution for
325 // VI: very suspect that 1 random number is
326 // and sampling below is not correct - s
327 G4double re = rndmEngine->flat();
328
329 G4double a = std::sqrt(16./xb - 3. - 36.*re
330 G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn +
331 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0
332
333 G4double photonEnergy1 = photonEnergy - ener
334 // resto al foton la energia del electron de
335 G4double positronTotEnergy = std::max(epsilo
336 G4double electronTotEnergy = std::max(photon
337
338 static const G4double a1 = 1.6;
339 static const G4double a2 = 0.5333333333;
340 G4double uu = -G4Log(rndmEngine->flat()*rndm
341 G4double u = (0.25 > rndmEngine->flat()) ? u
342
343 G4double thetaEle = u*electron_mass_c2/elect
344 G4double sinte = std::sin(thetaEle);
345 G4double coste = std::cos(thetaEle);
346
347 G4double thetaPos = u*electron_mass_c2/posit
348 G4double sintp = std::sin(thetaPos);
349 G4double costp = std::cos(thetaPos);
350
351 G4double phi = twopi * rndmEngine->flat();
352 G4double sinp = std::sin(phi);
353 G4double cosp = std::cos(phi);
354
355 // Kinematics of the created pair:
356 // the electron and positron are assumed to
357 // distribution with respect to the Z axis a
358
359 G4double electronKineEnergy = electronTotEne
360
361 G4ThreeVector electronDirection (sinte*cosp,
362 electronDirection.rotateUz(photonDirection);
363
364 G4DynamicParticle* particle1 = new G4Dynamic
365
366 electronKineEnergy);
367
368 G4double positronKineEnergy = positronTotEne
369
370 G4ThreeVector positronDirection (-sintp*cosp
371 positronDirection.rotateUz(photonDirection);
372
373 // Create G4DynamicParticle object for the p
374 G4DynamicParticle* particle2 = new G4Dynamic
375
376 // Fill output vector
377
378 fvect->push_back(particle1);
379 fvect->push_back(particle2);
380
381 if(particle3) { fvect->push_back(particle3);
382
383 // kill incident photon
384 fParticleChange->SetProposedKineticEnergy(0.
385 fParticleChange->ProposeTrackStatus(fStopAnd
386 }
387
388 //....oooOO0OOooo........oooOO0OOooo........oo
389
390 #include "G4AutoLock.hh"
391 namespace { G4Mutex BoldyshevTripletModelMutex
392
393 void G4BoldyshevTripletModel::InitialiseForEle
394 const G4ParticleDefinition*, G4int Z)
395 {
396 G4AutoLock l(&BoldyshevTripletModelMutex);
397 // G4cout << "G4BoldyshevTripletModel::Initi
398 // << Z << G4endl;
399 if(!data[Z]) { ReadData(Z); }
400 l.unlock();
401 }
402
403 //....oooOO0OOooo........oooOO0OOooo........oo
404