| 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 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 13 Jan 2010 L Pandola First implementa
32 // 24 May 2011 L Pandola Renamed (make v2
33 // 18 Sep 2013 L Pandola Migration to MT
34 // data and create
35 //
36
37 #include "G4PenelopeGammaConversionModel.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4MaterialCutsCouple.hh"
42 #include "G4ProductionCutsTable.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4Element.hh"
45 #include "G4Gamma.hh"
46 #include "G4Electron.hh"
47 #include "G4Positron.hh"
48 #include "G4PhysicsFreeVector.hh"
49 #include "G4MaterialCutsCouple.hh"
50 #include "G4AutoLock.hh"
51 #include "G4Exp.hh"
52
53 //....oooOO0OOooo........oooOO0OOooo........oo
54 const G4int G4PenelopeGammaConversionModel::fM
55 G4PhysicsFreeVector* G4PenelopeGammaConversion
56 G4double G4PenelopeGammaConversionModel::fAtom
57 1.2281e+02,7.3167e+01,6.9
58 6.4696e+01,6.1228e+01,5.7
59 5.0787e+01,4.7851e+01,4.6
60 4.4503e+01,4.3815e+01,4.3
61 4.1586e+01,4.0953e+01,4.0
62 3.9756e+01,3.9144e+01,3.8
63 3.7174e+01,3.6663e+01,3.5
64 3.4688e+01,3.4197e+01,3.3
65 3.3068e+01,3.2740e+01,3.2
66 3.1884e+01,3.1622e+01,3.1
67 3.0950e+01,3.0758e+01,3.0
68 3.0097e+01,2.9832e+01,2.9
69 2.9247e+01,2.9085e+01,2.8
70 2.8580e+01,2.8442e+01,2.8
71 2.7973e+01,2.7819e+01,2.7
72 2.7285e+01,2.7093e+01,2.6
73 2.6516e+01,2.6304e+01,2.6
74 2.5730e+01,2.5577e+01,2.5
75 2.5100e+01,2.4941e+01,2.4
76 2.4506e+01,2.4391e+01,2.4
77 2.4039e+01,2.3922e+01,2.3
78 2.3621e+01,2.3523e+01,2.3
79 2.3238e+01,2.3139e+01,2.3
80 2.2833e+01,2.2694e+01,2.2
81 2.2446e+01,2.2358e+01,2.2
82
83 //....oooOO0OOooo........oooOO0OOooo........oo
84
85 G4PenelopeGammaConversionModel::G4PenelopeGamm
86 const G4String& nam)
87 :G4VEmModel(nam),fParticleChange(nullptr),fP
88 fEffectiveCharge(nullptr),fMaterialInvScree
89 fScreeningFunction(nullptr),fIsInitialised(
90 {
91 fIntrinsicLowEnergyLimit = 2.0*electron_mass
92 fIntrinsicHighEnergyLimit = 100.0*GeV;
93 fSmallEnergy = 1.1*MeV;
94
95 if (part)
96 SetParticle(part);
97
98 // SetLowEnergyLimit(fIntrinsicLowEnergyLim
99 SetHighEnergyLimit(fIntrinsicHighEnergyLimit
100 //
101 fVerboseLevel= 0;
102 // Verbosity scale:
103 // 0 = nothing
104 // 1 = warning for energy non-conservation
105 // 2 = details of energy budget
106 // 3 = calculation of cross sections, file o
107 // 4 = entering in methods
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oo
111
112 G4PenelopeGammaConversionModel::~G4PenelopeGam
113 {
114 //Delete shared tables, they exist only in t
115 if (IsMaster() || fLocalTable)
116 {
117 for(G4int i=0; i<=fMaxZ; ++i)
118 {
119 if(fLogAtomicCrossSection[i]) {
120 delete fLogAtomicCrossSection[i];
121 fLogAtomicCrossSection[i] = nullptr;
122 }
123 }
124 if (fEffectiveCharge)
125 delete fEffectiveCharge;
126 if (fMaterialInvScreeningRadius)
127 delete fMaterialInvScreeningRadius;
128 if (fScreeningFunction)
129 delete fScreeningFunction;
130 }
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oo
134
135 void G4PenelopeGammaConversionModel::Initialis
136 const G4DataVector&)
137 {
138 if (fVerboseLevel > 3)
139 G4cout << "Calling G4PenelopeGammaConvers
140
141 SetParticle(part);
142
143 //Only the master model creates/fills/destro
144 if (IsMaster() && part == fParticle)
145 {
146 //delete old material data...
147 if (fEffectiveCharge)
148 {
149 delete fEffectiveCharge;
150 fEffectiveCharge = nullptr;
151 }
152 if (fMaterialInvScreeningRadius)
153 {
154 delete fMaterialInvScreeningRadius;
155 fMaterialInvScreeningRadius = nullptr;
156 }
157 if (fScreeningFunction)
158 {
159 delete fScreeningFunction;
160 fScreeningFunction = nullptr;
161 }
162 //and create new ones
163 fEffectiveCharge = new std::map<const G4
164 fMaterialInvScreeningRadius = new std::m
165 fScreeningFunction = new std::map<const
166
167 G4ProductionCutsTable* theCoupleTable =
168 G4ProductionCutsTable::GetProductionCutsTabl
169
170 for (G4int i=0;i<(G4int)theCoupleTable->
171 {
172 const G4Material* material =
173 theCoupleTable->GetMaterialCutsCouple(i)
174 const G4ElementVector* theElementVector =
175
176 for (std::size_t j=0;j<material->GetNumber
177 {
178 G4int iZ = theElementVector->at(j)->Ge
179 //read data files only in the master
180 if (iZ <= fMaxZ && !fLogAtomicCrossSe
181 ReadDataFile(iZ);
182 }
183
184 //check if material data are available
185 if (!fEffectiveCharge->count(material))
186 InitializeScreeningFunctions(material);
187 }
188 if (fVerboseLevel > 0) {
189 G4cout << "Penelope Gamma Conversion model v
190 << "Energy range: "
191 << LowEnergyLimit() / MeV << " MeV -
192 << HighEnergyLimit() / GeV << " GeV"
193 << G4endl;
194 }
195 }
196 if(fIsInitialised) return;
197 fParticleChange = GetParticleChangeForGamma(
198 fIsInitialised = true;
199 }
200
201 //....oooOO0OOooo........oooOO0OOooo........oo
202
203 void G4PenelopeGammaConversionModel::Initialis
204 G4VEmModel *masterModel)
205 {
206 if (fVerboseLevel > 3)
207 G4cout << "Calling G4PenelopeGammaConvers
208 //
209 //Check that particle matches: one might hav
210 //for e+ and e-).
211 //
212 if (part == fParticle)
213 {
214 //Get the const table pointers from the
215 const G4PenelopeGammaConversionModel* th
216 static_cast<G4PenelopeGammaConversionModel*>
217
218 //Copy pointers to the data tables
219 fEffectiveCharge = theModel->fEffectiveC
220 fMaterialInvScreeningRadius = theModel->
221 fScreeningFunction = theModel->fScreenin
222 for(G4int i=0; i<=fMaxZ; ++i)
223 fLogAtomicCrossSection[i] = theModel->fLogAt
224
225 //Same verbosity for all workers, as the
226 fVerboseLevel = theModel->fVerboseLevel;
227 }
228
229 return;
230 }
231
232 //....oooOO0OOooo........oooOO0OOooo........oo
233 namespace { G4Mutex PenelopeGammaConversionMo
234
235 G4double G4PenelopeGammaConversionModel::Compu
236 const G4ParticleDefinition
237 G4double energy,
238 G4double Z, G4double,
239 G4double, G4double)
240 {
241 //
242 // Penelope model v2008.
243 // Cross section (including triplet producti
244 // through the G4CrossSectionHandler utility
245 // M.J. Berger and J.H. Hubbel (XCOM), Repor
246 //
247
248 if (energy < fIntrinsicLowEnergyLimit)
249 return 0;
250
251 G4int iZ = G4int(Z);
252
253 if (!fLogAtomicCrossSection[iZ])
254 {
255 //If we are here, it means that Initial
256 //not filled up. This can happen in a U
257 if (fVerboseLevel > 0)
258 {
259 //Issue a G4Exception (warning) only in v
260 G4ExceptionDescription ed;
261 ed << "Unable to retrieve the cross secti
262 ed << "This can happen only in Unit Tests
263 G4Exception("G4PenelopeGammaConversionMod
264 "em2018",JustWarning,ed);
265 }
266 //protect file reading via autolock
267 G4AutoLock lock(&PenelopeGammaConversio
268 ReadDataFile(iZ);
269 lock.unlock();
270 fLocalTable = true;
271 }
272 G4double cs = 0;
273 G4double logene = G4Log(energy);
274 G4PhysicsFreeVector* theVec = fLogAtomicCros
275 G4double logXS = theVec->Value(logene);
276 cs = G4Exp(logXS);
277
278 if (fVerboseLevel > 2)
279 G4cout << "Gamma conversion cross section
280 " = " << cs/barn << " barn" << G4endl;
281 return cs;
282 }
283
284 //....oooOO0OOooo........oooOO0OOooo........oo
285
286 void
287 G4PenelopeGammaConversionModel::SampleSecondar
288 const G4MaterialCutsCouple* coup
289 const G4DynamicParticle* aDynami
290 G4double,
291 G4double)
292 {
293 //
294 // Penelope model v2008.
295 // Final state is sampled according to the B
296 // corrections, according to the semi-empiri
297 // J. Baro' et al., Radiat. Phys. Chem. 44
298 //
299 // The model uses the high energy Coulomb co
300 // H. Davies et al., Phys. Rev. 93 (1954) 7
301 // and atomic screening radii tabulated from
302 // J.H. Hubbel et al., J. Phys. Chem. Ref.
303 // for Z= 1 to 92.
304 //
305 if (fVerboseLevel > 3)
306 G4cout << "Calling SamplingSecondaries() o
307
308 G4double photonEnergy = aDynamicGamma->GetKi
309
310 // Always kill primary
311 fParticleChange->ProposeTrackStatus(fStopAnd
312 fParticleChange->SetProposedKineticEnergy(0.
313
314 if (photonEnergy <= fIntrinsicLowEnergyLimit
315 {
316 fParticleChange->ProposeLocalEnergyDepos
317 return ;
318 }
319
320 G4ParticleMomentum photonDirection = aDynami
321 const G4Material* mat = couple->GetMaterial(
322
323 //Either Initialize() was not called, or we
324 //not invoked
325 if (!fEffectiveCharge)
326 {
327 //create a **thread-local** version of t
328 //Unit Tests
329 fLocalTable = true;
330 fEffectiveCharge = new std::map<const G4
331 fMaterialInvScreeningRadius = new std::m
332 fScreeningFunction = new std::map<const
333 }
334
335 if (!fEffectiveCharge->count(mat))
336 {
337 //If we are here, it means that Initiali
338 //not filled up. This can happen in a Un
339 if (fVerboseLevel > 0)
340 {
341 //Issue a G4Exception (warning) only in ve
342 G4ExceptionDescription ed;
343 ed << "Unable to allocate the EffectiveCha
344 mat->GetName() << G4endl;
345 ed << "This can happen only in Unit Tests"
346 G4Exception("G4PenelopeGammaConversionMode
347 "em2019",JustWarning,ed);
348 }
349 //protect file reading via autolock
350 G4AutoLock lock(&PenelopeGammaConversion
351 InitializeScreeningFunctions(mat);
352 lock.unlock();
353 }
354
355 // eps is the fraction of the photon energy
356 G4double eps = 0;
357 G4double eki = electron_mass_c2/photonEnergy
358
359 //Do it fast for photon energy < 1.1 MeV (cl
360 if (photonEnergy < fSmallEnergy)
361 eps = eki + (1.0-2.0*eki)*G4UniformRand();
362 else
363 {
364 //Complete calculation
365 G4double effC = fEffectiveCharge->find(m
366 G4double alz = effC*fine_structure_const
367 G4double T = std::sqrt(2.0*eki);
368 G4double F00=(-1.774-1.210e1*alz+1.118e1
369 +(8.523+7.326e1*alz-4.441e1*alz*alz)*
370 -(1.352e1+1.211e2*alz-9.641e1*alz*alz
371 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T
372
373 G4double F0b = fScreeningFunction->find(
374 G4double g0 = F0b + F00;
375 G4double invRad = fMaterialInvScreeningR
376 G4double bmin = 4.0*eki/invRad;
377 std::pair<G4double,G4double> scree = Ge
378 G4double g1 = scree.first;
379 G4double g2 = scree.second;
380 G4double g1min = g1+g0;
381 G4double g2min = g2+g0;
382 G4double xr = 0.5-eki;
383 G4double a1 = 2.*g1min*xr*xr/3.;
384 G4double p1 = a1/(a1+g2min);
385
386 G4bool loopAgain = false;
387 //Random sampling of eps
388 do{
389 loopAgain = false;
390 if (G4UniformRand() <= p1)
391 {
392 G4double ru2m1 = 2.0*G4UniformRand()-1.
393 if (ru2m1 < 0)
394 eps = 0.5-xr*std::pow(std::abs(ru2m1),
395 else
396 eps = 0.5+xr*std::pow(ru2m1,1./3.);
397 G4double B = eki/(invRad*eps*(1.0-eps));
398 scree = GetScreeningFunctions(B);
399 g1 = scree.first;
400 g1 = std::max(g1+g0,0.);
401 if (G4UniformRand()*g1min > g1)
402 loopAgain = true;
403 }
404 else
405 {
406 eps = eki+2.0*xr*G4UniformRand();
407 G4double B = eki/(invRad*eps*(1.0-eps));
408 scree = GetScreeningFunctions(B);
409 g2 = scree.second;
410 g2 = std::max(g2+g0,0.);
411 if (G4UniformRand()*g2min > g2)
412 loopAgain = true;
413 }
414 }while(loopAgain);
415 }
416 if (fVerboseLevel > 4)
417 G4cout << "Sampled eps = " << eps << G4end
418
419 G4double electronTotEnergy = eps*photonEnerg
420 G4double positronTotEnergy = (1.0-eps)*photo
421
422 // Scattered electron (positron) angles. ( Z
423
424 //electron kinematics
425 G4double electronKineEnergy = std::max(0.,el
426 G4double costheta_el = G4UniformRand()*2.0-1
427 G4double kk = std::sqrt(electronKineEnergy*(
428 costheta_el = (costheta_el*electronTotEnergy
429 G4double phi_el = twopi * G4UniformRand() ;
430 G4double dirX_el = std::sqrt(1.-costheta_el*
431 G4double dirY_el = std::sqrt(1.-costheta_el*
432 G4double dirZ_el = costheta_el;
433
434 //positron kinematics
435 G4double positronKineEnergy = std::max(0.,po
436 G4double costheta_po = G4UniformRand()*2.0-1
437 kk = std::sqrt(positronKineEnergy*(positronK
438 costheta_po = (costheta_po*positronTotEnergy
439 G4double phi_po = twopi * G4UniformRand() ;
440 G4double dirX_po = std::sqrt(1.-costheta_po*
441 G4double dirY_po = std::sqrt(1.-costheta_po*
442 G4double dirZ_po = costheta_po;
443
444 // Kinematics of the created pair:
445 // the electron and positron are assumed to
446 // distribution with respect to the Z axis a
447 G4double localEnergyDeposit = 0. ;
448
449 if (electronKineEnergy > 0.0)
450 {
451 G4ThreeVector electronDirection ( dirX_e
452 electronDirection.rotateUz(photonDirecti
453 G4DynamicParticle* electron = new G4Dyna
454 electronDirection,
455 electronKineEnergy);
456 fvect->push_back(electron);
457 }
458 else
459 {
460 localEnergyDeposit += electronKineEnergy
461 electronKineEnergy = 0;
462 }
463
464 //Generate the positron. Real particle in an
465 //threshold, produce it at rest
466 if (positronKineEnergy < 0.0)
467 {
468 localEnergyDeposit += positronKineEnergy
469 positronKineEnergy = 0; //produce it at
470 }
471 G4ThreeVector positronDirection(dirX_po,dirY
472 positronDirection.rotateUz(photonDirection);
473 G4DynamicParticle* positron = new G4DynamicP
474 positronDirection, positronK
475 fvect->push_back(positron);
476
477 //Add rest of energy to the local energy dep
478 fParticleChange->ProposeLocalEnergyDeposit(l
479
480 if (fVerboseLevel > 1)
481 {
482 G4cout << "-----------------------------
483 G4cout << "Energy balance from G4Penelop
484 G4cout << "Incoming photon energy: " <<
485 G4cout << "-----------------------------
486 if (electronKineEnergy)
487 G4cout << "Electron (explicitly produced) "
488 << G4endl;
489 if (positronKineEnergy)
490 G4cout << "Positron (not at rest) " << posit
491 G4cout << "Rest masses of e+/- " << 2.0*
492 if (localEnergyDeposit)
493 G4cout << "Local energy deposit " << localEn
494 G4cout << "Total final state: " << (elec
495 localEnergyDeposit+2.0*electron_ma
496 " keV" << G4endl;
497 G4cout << "-----------------------------
498 }
499 if (fVerboseLevel > 0)
500 {
501 G4double energyDiff = std::fabs(electron
502 localEnergyDeposit+2.0*electron_
503 if (energyDiff > 0.05*keV)
504 G4cout << "Warning from G4PenelopeGammaConve
505 << (electronKineEnergy+positronKineEn
506 localEnergyDeposit+2.0*electron_mass_c2
507 << " keV (final) vs. " << photonEnerg
508 }
509 }
510
511 //....oooOO0OOooo........oooOO0OOooo........oo
512
513 void G4PenelopeGammaConversionModel::ReadDataF
514 {
515 if (!IsMaster())
516 //Should not be here!
517 G4Exception("G4PenelopeGammaConversionMode
518 "em0100",FatalException,"Worker thread in
519
520 if (fVerboseLevel > 2)
521 {
522 G4cout << "G4PenelopeGammaConversionMode
523 G4cout << "Going to read Gamma Conversio
524 }
525
526 const char* path = G4FindDataDir("G4LEDATA
527 if(!path)
528 {
529 G4String excep =
530 "G4PenelopeGammaConversionModel - G4LEDATA e
531 G4Exception("G4PenelopeGammaConversionMo
532 "em0006",FatalException,excep);
533 return;
534 }
535
536 /*
537 Read the cross section file
538 */
539 std::ostringstream ost;
540 if (Z>9)
541 ost << path << "/penelope/pairproduction/p
542 else
543 ost << path << "/penelope/pairproduction/p
544 std::ifstream file(ost.str().c_str());
545 if (!file.is_open())
546 {
547 G4String excep = "G4PenelopeGammaConvers
548 G4String(ost.str()) + " not found!";
549 G4Exception("G4PenelopeGammaConversionMo
550 "em0003",FatalException,excep);
551 }
552
553 //I have to know in advance how many points
554 //to initialize the G4PhysicsFreeVector()
555 std::size_t ndata=0;
556 G4String line;
557 while( getline(file, line) )
558 ndata++;
559 ndata -= 1; //remove one header line
560
561 file.clear();
562 file.close();
563 file.open(ost.str().c_str());
564 G4int readZ =0;
565 file >> readZ;
566
567 if (fVerboseLevel > 3)
568 G4cout << "Element Z=" << Z << G4endl;
569
570 //check the right file is opened.
571 if (readZ != Z)
572 {
573 G4ExceptionDescription ed;
574 ed << "Corrupted data file for Z=" << Z
575 G4Exception("G4PenelopeGammaConversionMo
576 "em0005",FatalException,ed);
577 }
578
579 fLogAtomicCrossSection[Z] = new G4PhysicsFre
580 G4double ene=0,xs=0;
581 for (std::size_t i=0;i<ndata;++i)
582 {
583 file >> ene >> xs;
584 //dimensional quantities
585 ene *= eV;
586 xs *= barn;
587 if (xs < 1e-40*cm2) //protection against
588 xs = 1e-40*cm2;
589 fLogAtomicCrossSection[Z]->PutValue(i,G4
590 }
591 file.close();
592
593 return;
594 }
595
596 //....oooOO0OOooo........oooOO0OOooo........oo
597
598 void G4PenelopeGammaConversionModel::Initializ
599 {
600 // This is subroutine GPPa0 of Penelope
601 //
602 // 1) calculate the effective Z for the purp
603 //
604 G4double zeff = 0;
605 G4int intZ = 0;
606 G4int nElements = (G4int)material->GetNumber
607 const G4ElementVector* elementVector = mater
608
609 //avoid calculations if only one building el
610 if (nElements == 1)
611 {
612 zeff = (*elementVector)[0]->GetZ();
613 intZ = (G4int) zeff;
614 }
615 else // many elements...let's do the calcula
616 {
617 const G4double* fractionVector = materia
618
619 G4double atot = 0;
620 for (G4int i=0;i<nElements;i++)
621 {
622 G4double Zelement = (*elementVector)[i]->G
623 G4double Aelement = (*elementVector)[i]->G
624 atot += Aelement*fractionVector[i];
625 zeff += Zelement*Aelement*fractionVector[i
626 }
627 atot /= material->GetTotNbOfAtomsPerVolu
628 zeff /= (material->GetTotNbOfAtomsPerVol
629
630 intZ = (G4int) (zeff+0.25);
631 if (intZ <= 0)
632 intZ = 1;
633 if (intZ > fMaxZ)
634 intZ = fMaxZ;
635 }
636
637 if (fEffectiveCharge)
638 fEffectiveCharge->insert(std::make_pair(ma
639
640 //
641 // 2) Calculate Coulomb Correction
642 //
643 G4double alz = fine_structure_const*zeff;
644 G4double alzSquared = alz*alz;
645 G4double fc = alzSquared*(0.202059-alzSquar
646 (0.03693-alzSquared*
647 (0.00835-alzSquared*(0.00201-alzSq
648 (0.00049-alzSquared*
649 (0.00012-alzSquared*0.00003)))
650 +1.0/(alzSquared+1.0));
651 //
652 // 3) Screening functions and low-energy cor
653 //
654 G4double matRadius = 2.0/ fAtomicScreeningRa
655 if (fMaterialInvScreeningRadius)
656 fMaterialInvScreeningRadius->insert(std::m
657
658 std::pair<G4double,G4double> myPair(0,0);
659 G4double f0a = 4.0*G4Log(fAtomicScreeningRad
660 G4double f0b = f0a - 4.0*fc;
661 myPair.first = f0a;
662 myPair.second = f0b;
663
664 if (fScreeningFunction)
665 fScreeningFunction->insert(std::make_pair(
666
667 if (fVerboseLevel > 2)
668 {
669 G4cout << "Average Z for material " << m
670 zeff << G4endl;
671 G4cout << "Effective radius for material
672 fAtomicScreeningRadius[intZ] << " m_e*c/hbar
673 matRadius << G4endl;
674 G4cout << "Screening parameters F0 for m
675 f0a << "," << f0b << G4endl;
676 }
677 return;
678 }
679
680 //....oooOO0OOooo........oooOO0OOooo........oo
681
682 std::pair<G4double,G4double>
683 G4PenelopeGammaConversionModel::GetScreeningFu
684 {
685 // This is subroutine SCHIFF of Penelope
686 //
687 // Screening Functions F1(B) and F2(B) in th
688 // section for pair production
689 //
690 std::pair<G4double,G4double> result(0.,0.);
691 G4double BSquared = B*B;
692 G4double f1 = 2.0-2.0*G4Log(1.0+BSquared);
693 G4double f2 = f1 - 6.66666666e-1; // (-2/3)
694 if (B < 1.0e-10)
695 f1 = f1-twopi*B;
696 else
697 {
698 G4double a0 = 4.0*B*std::atan(1./B);
699 f1 = f1 - a0;
700 f2 += 2.0*BSquared*(4.0-a0-3.0*G4Log((1.
701 }
702 G4double g1 = 0.5*(3.0*f1-f2);
703 G4double g2 = 0.25*(3.0*f1+f2);
704
705 result.first = g1;
706 result.second = g2;
707
708 return result;
709 }
710
711 //....oooOO0OOooo........oooOO0OOooo........oo
712
713 void G4PenelopeGammaConversionModel::SetPartic
714 {
715 if(!fParticle) {
716 fParticle = p;
717 }
718 }
719