| 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 // 15 Feb 2010 L Pandola Implementation
32 // 18 Mar 2010 L Pandola Removed GetAtomsPe
33 // to G4PenelopeOsc
34 // 01 Feb 2011 L Pandola Suppress fake ener
35 // active.
36 // Make sure that flu
37 // if above thresho
38 // 24 May 2011 L Pandola Renamed (make v200
39 // 10 Jun 2011 L Pandola Migrate atomic dee
40 // 09 Oct 2013 L Pandola Migration to MT
41 // 25 Jul 2023 D Iuso Fix for possible i
42 //
43 #include "G4PenelopeComptonModel.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4ParticleDefinition.hh"
47 #include "G4MaterialCutsCouple.hh"
48 #include "G4DynamicParticle.hh"
49 #include "G4VEMDataSet.hh"
50 #include "G4PhysicsTable.hh"
51 #include "G4PhysicsLogVector.hh"
52 #include "G4AtomicTransitionManager.hh"
53 #include "G4AtomicShell.hh"
54 #include "G4Gamma.hh"
55 #include "G4Electron.hh"
56 #include "G4PenelopeOscillatorManager.hh"
57 #include "G4PenelopeOscillator.hh"
58 #include "G4LossTableManager.hh"
59 #include "G4Exp.hh"
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 G4PenelopeComptonModel::G4PenelopeComptonModel
64 const G4String& nam)
65 :G4VEmModel(nam),fParticleChange(nullptr),fP
66 fAtomDeexcitation(nullptr),
67 fOscManager(nullptr),fIsInitialised(false)
68 {
69 fIntrinsicLowEnergyLimit = 100.0*eV;
70 fIntrinsicHighEnergyLimit = 100.0*GeV;
71 SetHighEnergyLimit(fIntrinsicHighEnergyLimit
72 //
73 fOscManager = G4PenelopeOscillatorManager::G
74
75 if (part)
76 SetParticle(part);
77
78 fVerboseLevel= 0;
79 // Verbosity scale:
80 // 0 = nothing
81 // 1 = warning for energy non-conservation
82 // 2 = details of energy budget
83 // 3 = calculation of cross sections, file o
84 // 4 = entering in methods
85
86 //Mark this model as "applicable" for atomic
87 SetDeexcitationFlag(true);
88
89 fTransitionManager = G4AtomicTransitionManag
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 G4PenelopeComptonModel::~G4PenelopeComptonMode
95 {;}
96
97 //....oooOO0OOooo........oooOO0OOooo........oo
98
99 void G4PenelopeComptonModel::Initialise(const
100 const G4DataVector&)
101 {
102 if (fVerboseLevel > 3)
103 G4cout << "Calling G4PenelopeComptonModel:
104
105 fAtomDeexcitation = G4LossTableManager::Inst
106 //Issue warning if the AtomicDeexcitation ha
107 if (!fAtomDeexcitation)
108 {
109 G4cout << G4endl;
110 G4cout << "WARNING from G4PenelopeCompto
111 G4cout << "Atomic de-excitation module i
112 G4cout << "any fluorescence/Auger emissi
113 G4cout << "Please make sure this is inte
114 }
115
116 SetParticle(part);
117
118 if (IsMaster() && part == fParticle)
119 {
120
121 if (fVerboseLevel > 0)
122 {
123 G4cout << "Penelope Compton model v2008 is
124 << "Energy range: "
125 << LowEnergyLimit() / keV << " keV - "
126 << HighEnergyLimit() / GeV << " GeV";
127 }
128 //Issue a warning, if the model is going
129 //energy which is outside the validity o
130 if (LowEnergyLimit() < fIntrinsicLowEner
131 {
132 G4ExceptionDescription ed;
133 ed << "Using the Penelope Compton model ou
134 << G4endl;
135 ed << "-> LowEnergyLimit() in process = "
136 ed << "-> Instrinsic low-energy limit = "
137 << G4endl;
138 ed << "Result of the simulation have to be
139 G4Exception("G4PenelopeComptonModel::Initi
140 "em2100",JustWarning,ed);
141 }
142 }
143
144 if(fIsInitialised) return;
145 fParticleChange = GetParticleChangeForGamma(
146 fIsInitialised = true;
147
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oo
151
152 void G4PenelopeComptonModel::InitialiseLocal(c
153
154 {
155 if (fVerboseLevel > 3)
156 G4cout << "Calling G4PenelopeComptonModel
157 //
158 //Check that particle matches: one might hav
159 //for e+ and e-).
160 //
161 if (part == fParticle)
162 {
163 //Get the const table pointers from the
164 const G4PenelopeComptonModel* theModel =
165 static_cast<G4PenelopeComptonModel*> (
166
167 //Same verbosity for all workers, as the
168 fVerboseLevel = theModel->fVerboseLevel;
169 }
170 return;
171 }
172
173
174 //....oooOO0OOooo........oooOO0OOooo........oo
175
176 G4double G4PenelopeComptonModel::CrossSectionP
177
178
179
180
181 {
182 // Penelope model v2008 to calculate the Com
183 // D. Brusa et al., Nucl. Instrum. Meth. A 3
184 //
185 // The cross section for Compton scattering
186 // formula for energy > 5 MeV.
187 // For E < 5 MeV it is used a parametrizatio
188 // which is integrated numerically in cos(th
189 // The parametrization includes the J(p)
190 // distribution profiles for the atomic shel
191 // from F. Biggs et al., At. Data Nucl. Data
192 //
193 if (fVerboseLevel > 3)
194 G4cout << "Calling CrossSectionPerVolume()
195 SetupForMaterial(p, material, energy);
196
197 G4double cs = 0;
198 //Force null cross-section if below the low-
199 if (energy < LowEnergyLimit())
200 return cs;
201
202 //Retrieve the oscillator table for this mat
203 G4PenelopeOscillatorTable* theTable = fOscMa
204
205 if (energy < 5*MeV) //explicit calculation f
206 {
207 size_t numberOfOscillators = theTable->s
208 for (size_t i=0;i<numberOfOscillators;i+
209 {
210 G4PenelopeOscillator* theOsc = (*theTable)
211 //sum contributions coming from each oscil
212 cs += OscillatorTotalCrossSection(energy,t
213 }
214 }
215 else //use Klein-Nishina for E>5 MeV
216 cs = KleinNishinaCrossSection(energy,mater
217
218 //cross sections are in units of pi*classic_
219 cs *= pi*classic_electr_radius*classic_elect
220
221 //Now, cs is the cross section *per molecule
222 //cross section per volume
223 G4double atomDensity = material->GetTotNbOfA
224 G4double atPerMol = fOscManager->GetAtomsPe
225
226 if (fVerboseLevel > 3)
227 G4cout << "Material " << material->GetName
228 "atoms per molecule" << G4endl;
229
230 G4double moleculeDensity = 0.;
231
232 if (atPerMol)
233 moleculeDensity = atomDensity/atPerMol;
234
235 G4double csvolume = cs*moleculeDensity;
236
237 if (fVerboseLevel > 2)
238 G4cout << "Compton mean free path at " <<
239 material->GetName() << " = " << (1
240 return csvolume;
241 }
242
243 //....oooOO0OOooo........oooOO0OOooo........oo
244
245 //This is a dummy method. Never inkoved by the
246 //a warning if one tries to get Cross Sections
247 //G4EmCalculator.
248 G4double G4PenelopeComptonModel::ComputeCrossS
249
250
251
252
253
254 {
255 G4cout << "*** G4PenelopeComptonModel -- WAR
256 G4cout << "Penelope Compton model v2008 does
257 G4cout << "so the result is always zero. For
258 G4cout << "GetCrossSectionPerVolume() or Get
259 return 0;
260 }
261
262 //....oooOO0OOooo........oooOO0OOooo........oo
263
264 void G4PenelopeComptonModel::SampleSecondaries
265 const G4MaterialCutsCouple* c
266 const G4DynamicParticle* aDyna
267 G4double,
268 G4double)
269 {
270 // Penelope model v2008 to sample the Compto
271 // D. Brusa et al., Nucl. Instrum. Meth. A 3
272 // The model determines also the original sh
273 // in order to produce fluorescence de-excit
274 //
275 // The final state for Compton scattering is
276 // formula for energy > 5 MeV. In this case,
277 // one can assume that the target electron i
278 // For E < 5 MeV it is used the parametrizat
279 // to sample the scattering angle and the en
280 // G4PenelopeComptonModel::DifferentialCross
281 // used to sample cos(theta). The efficiency
282 // nearly independent on the Z; typical valu
283 // respectively.
284 // The parametrization includes the J(p) dis
285 // tabulated
286 // from Hartree-Fock calculations from F. Bi
287 // Doppler broadening is included.
288 //
289
290 if (fVerboseLevel > 3)
291 G4cout << "Calling SampleSecondaries() of
292
293 G4double photonEnergy0 = aDynamicGamma->GetK
294
295 // do nothing below the threshold
296 // should never get here because the XS is z
297 if(photonEnergy0 < LowEnergyLimit())
298 return;
299
300 G4ParticleMomentum photonDirection0 = aDynam
301 const G4Material* material = couple->GetMate
302
303 G4PenelopeOscillatorTable* theTable = fOscMa
304
305 const G4int nmax = 64;
306 G4double rn[nmax]={0.0};
307 G4double pac[nmax]={0.0};
308
309 G4double S=0.0;
310 G4double epsilon = 0.0;
311 G4double cosTheta = 1.0;
312 G4double hartreeFunc = 0.0;
313 G4double oscStren = 0.0;
314 size_t numberOfOscillators = theTable->size(
315 size_t targetOscillator = 0;
316 G4double ionEnergy = 0.0*eV;
317
318 G4double ek = photonEnergy0/electron_mass_c2
319 G4double ek2 = 2.*ek+1.0;
320 G4double eks = ek*ek;
321 G4double ek1 = eks-ek2-1.0;
322
323 G4double taumin = 1.0/ek2;
324 //This is meant to fix a possible (rare) flo
325 //causing an infinite loop. The maximum of t
326 //be represented (i.e. ~ 1. - 1e-16). Fix by
327 static G4double taumax = std::nexttoward(1.0
328 if (fVerboseLevel > 3)
329 G4cout << "G4PenelopeComptonModel: maximum
330 //To here.
331 G4double a1 = G4Log(ek2);
332 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
333
334 G4double TST = 0;
335 G4double tau = 0.;
336
337 //If the incoming photon is above 5 MeV, the
338 //pure Klein-Nishina formula is used
339 if (photonEnergy0 > 5*MeV)
340 {
341 do{
342 do{
343 if ((a2*G4UniformRand()) < a1)
344 tau = std::pow(taumin,G4UniformRand());
345 else
346 tau = std::sqrt(1.0+G4UniformRand()*(tau
347 //rejection function
348 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(e
349 }while (G4UniformRand()> TST);
350 if (tau > taumax) tau = taumax; //prevent FP
351 epsilon=tau;
352 cosTheta = 1.0 - (1.0-tau)/(ek*tau);
353
354 //Target shell electrons
355 TST = fOscManager->GetTotalZ(material)*G4Uni
356 targetOscillator = numberOfOscillators-1; //
357 S=0.0;
358 G4bool levelFound = false;
359 for (size_t j=0;j<numberOfOscillators && !le
360 {
361 S += (*theTable)[j]->GetOscillatorStreng
362 if (S > TST)
363 {
364 targetOscillator = j;
365 levelFound = true;
366 }
367 }
368 //check whether the level is valid
369 ionEnergy = (*theTable)[targetOscillator]->G
370 }while((epsilon*photonEnergy0-photonEner
371 }
372 else //photonEnergy0 < 5 MeV
373 {
374 //Incoherent scattering function for the
375 G4double s0=0.0;
376 G4double pzomc=0.0;
377 G4double rni=0.0;
378 G4double aux=0.0;
379 for (size_t i=0;i<numberOfOscillators;i+
380 {
381 ionEnergy = (*theTable)[i]->GetIonisationE
382 if (photonEnergy0 > ionEnergy)
383 {
384 G4double aux2 = photonEnergy0*(photonE
385 hartreeFunc = (*theTable)[i]->GetHartr
386 oscStren = (*theTable)[i]->GetOscillat
387 pzomc = hartreeFunc*(aux2-electron_mas
388 (electron_mass_c2*std::sqrt(2.0*aux2+ionEn
389 if (pzomc > 0)
390 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+st
391 (std::sqrt(0.5)+std::sqrt(2.0)*
392 else
393 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::s
394 (std::sqrt(0.5)-std::sqrt(2.0)*pzom
395 s0 += oscStren*rni;
396 }
397 }
398 //Sampling tau
399 G4double cdt1 = 0.;
400 do
401 {
402 if ((G4UniformRand()*a2) < a1)
403 tau = std::pow(taumin,G4UniformRand());
404 else
405 tau = std::sqrt(1.0+G4UniformRand()*(tau
406 if (tau > taumax) tau = taumax; //prevent
407 cdt1 = (1.0-tau)/(ek*tau);
408 //Incoherent scattering function
409 S = 0.;
410 for (size_t i=0;i<numberOfOscillators;i++)
411 {
412 ionEnergy = (*theTable)[i]->GetIonisat
413 if (photonEnergy0 > ionEnergy) //sum o
414 {
415 aux = photonEnergy0*(photonEnergy0-ionEn
416 hartreeFunc = (*theTable)[i]->GetHartree
417 oscStren = (*theTable)[i]->GetOscillator
418 pzomc = hartreeFunc*(aux-electron_mass_c
419 (electron_mass_c2*std::sqrt(2.0*aux+io
420 if (pzomc > 0)
421 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0
422 (std::sqrt(0.5)+std::sqrt(2.0)*
423 else
424 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)-
425 (std::sqrt(0.5)-std::sqrt(2.0)*pzom
426 S += oscStren*rn[i];
427 pac[i] = S;
428 }
429 else
430 pac[i] = S-1e-6;
431 }
432 //Rejection function
433 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/
434 }while ((G4UniformRand()*s0) > TST);
435
436 cosTheta = 1.0 - cdt1;
437 G4double fpzmax=0.0,fpz=0.0;
438 G4double A=0.0;
439 //Target electron shell
440 do
441 {
442 do
443 {
444 TST = S*G4UniformRand();
445 targetOscillator = numberOfOscillators
446 G4bool levelFound = false;
447 for (size_t i=0;i<numberOfOscillators
448 {
449 if (pac[i]>TST)
450 {
451 targetOscillator = i;
452 levelFound = true;
453 }
454 }
455 A = G4UniformRand()*rn[targetOscillato
456 hartreeFunc = (*theTable)[targetOscill
457 oscStren = (*theTable)[targetOscillato
458 if (A < 0.5)
459 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-G4Lo
460 (std::sqrt(2.0)*hartreeFunc);
461 else
462 pzomc = (std::sqrt(0.5-G4Log(2.0-2.0*A))-s
463 (std::sqrt(2.0)*hartreeFunc);
464 } while (pzomc < -1);
465
466 // F(EP) rejection
467 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
468 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau
469 if (AF > 0)
470 fpzmax = 1.0+AF*0.2;
471 else
472 fpzmax = 1.0-AF*0.2;
473 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),
474 }while ((fpzmax*G4UniformRand())>fpz);
475
476 //Energy of the scattered photon
477 G4double T = pzomc*pzomc;
478 G4double b1 = 1.0-T*tau*tau;
479 G4double b2 = 1.0-T*tau*cosTheta;
480 if (pzomc > 0.0)
481 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2
482 else
483 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2
484 } //energy < 5 MeV
485
486 //Ok, the kinematics has been calculated.
487 G4double sinTheta = std::sqrt(1-cosTheta*cos
488 G4double phi = twopi * G4UniformRand() ;
489 G4double dirx = sinTheta * std::cos(phi);
490 G4double diry = sinTheta * std::sin(phi);
491 G4double dirz = cosTheta ;
492
493 // Update G4VParticleChange for the scattere
494 G4ThreeVector photonDirection1(dirx,diry,dir
495 photonDirection1.rotateUz(photonDirection0);
496 fParticleChange->ProposeMomentumDirection(ph
497
498 G4double photonEnergy1 = epsilon * photonEne
499
500 if (photonEnergy1 > 0.)
501 fParticleChange->SetProposedKineticEnergy(
502 else
503 {
504 fParticleChange->SetProposedKineticEnergy(
505 fParticleChange->ProposeTrackStatus(fStopA
506 }
507
508 //Create scattered electron
509 G4double diffEnergy = photonEnergy0*(1-epsil
510 ionEnergy = (*theTable)[targetOscillator]->G
511
512 G4double Q2 =
513 photonEnergy0*photonEnergy0+photonEnergy1*
514 G4double cosThetaE = 0.; //scattering angle
515
516 if (Q2 > 1.0e-12)
517 cosThetaE = (photonEnergy0-photonEnergy1*c
518 else
519 cosThetaE = 1.0;
520 G4double sinThetaE = std::sqrt(1-cosThetaE*c
521
522 //Now, try to handle fluorescence
523 //Notice: merged levels are indicated with Z
524 G4int shFlag = (*theTable)[targetOscillator]
525 G4int Z = (G4int) (*theTable)[targetOscillat
526
527 //initialize here, then check photons create
528 G4double bindingEnergy = 0.*eV;
529 const G4AtomicShell* shell = 0;
530
531 //Real level
532 if (Z > 0 && shFlag<30)
533 {
534 shell = fTransitionManager->Shell(Z,shFl
535 bindingEnergy = shell->BindingEnergy();
536 }
537
538 G4double ionEnergyInPenelopeDatabase = ionEn
539 //protection against energy non-conservation
540 ionEnergy = std::max(bindingEnergy,ionEnergy
541
542 //subtract the excitation energy. If not emi
543 //the ionization energy is deposited as loca
544 G4double eKineticEnergy = diffEnergy - ionEn
545 G4double localEnergyDeposit = ionEnergy;
546 G4double energyInFluorescence = 0.; //testin
547 G4double energyInAuger = 0; //testing purpos
548
549 if (eKineticEnergy < 0)
550 {
551 //It means that there was some problem/m
552 //Try to make it work
553 //In this case available Energy (diffEne
554 //Full residual energy is deposited loca
555 localEnergyDeposit = diffEnergy;
556 eKineticEnergy = 0.0;
557 }
558
559 //the local energy deposit is what remains:
560 //Notice: shell might be nullptr (invalid!)
561 //Now, take care of fluorescence, if require
562 if (fAtomDeexcitation && shell)
563 {
564 G4int index = couple->GetIndex();
565 if (fAtomDeexcitation->CheckDeexcitation
566 {
567 size_t nBefore = fvect->size();
568 fAtomDeexcitation->GenerateParticles(fvect
569 size_t nAfter = fvect->size();
570
571 if (nAfter > nBefore) //actual production
572 {
573 for (size_t j=nBefore;j<nAfter;j++) //
574 {
575 G4double itsEnergy = ((*fvect)[j])->GetK
576 if (itsEnergy < localEnergyDeposit) // v
577 {
578 localEnergyDeposit -= itsEnergy;
579 if (((*fvect)[j])->GetParticleDefini
580 energyInFluorescence += itsEnergy;
581 else if (((*fvect)[j])->GetParticleD
582 G4Electron::Definition())
583 energyInAuger += itsEnergy;
584 }
585 else //invalid secondary: takes more tha
586 {
587 delete (*fvect)[j];
588 (*fvect)[j] = nullptr;
589 }
590 }
591 }
592
593 }
594 }
595
596 //Always produce explicitly the electron
597 G4DynamicParticle* electron = 0;
598
599 G4double xEl = sinThetaE * std::cos(phi+pi);
600 G4double yEl = sinThetaE * std::sin(phi+pi);
601 G4double zEl = cosThetaE;
602 G4ThreeVector eDirection(xEl,yEl,zEl); //ele
603 eDirection.rotateUz(photonDirection0);
604 electron = new G4DynamicParticle (G4Electron
605 eDirection,eKineticEnergy) ;
606 fvect->push_back(electron);
607
608 if (localEnergyDeposit < 0) //Should not be:
609 {
610 G4Exception("G4PenelopeComptonModel::Sam
611 "em2099",JustWarning,"WARNING: Negative
612 localEnergyDeposit=0.;
613 }
614 fParticleChange->ProposeLocalEnergyDeposit(l
615
616 G4double electronEnergy = 0.;
617 if (electron)
618 electronEnergy = eKineticEnergy;
619 if (fVerboseLevel > 1)
620 {
621 G4cout << "-----------------------------
622 G4cout << "Energy balance from G4Penelop
623 G4cout << "Incoming photon energy: " <<
624 G4cout << "-----------------------------
625 G4cout << "Scattered photon: " << photon
626 G4cout << "Scattered electron " << elect
627 if (energyInFluorescence)
628 G4cout << "Fluorescence x-rays: " << energyI
629 if (energyInAuger)
630 G4cout << "Auger electrons: " << energyInAug
631 G4cout << "Local energy deposit " << loc
632 G4cout << "Total final state: " << (phot
633 localEnergyDeposit+energyInAuger)/
634 " keV" << G4endl;
635 G4cout << "-----------------------------
636 }
637 if (fVerboseLevel > 0)
638 {
639 G4double energyDiff = std::fabs(photonEn
640 electronEnergy+energyInFluoresce
641 localEnergyDeposit+energyInAuger
642 if (energyDiff > 0.05*keV)
643 G4cout << "Warning from G4PenelopeCompton: p
644 (photonEnergy1+electronEnergy+energyInFluo
645 " keV (final) vs. " <<
646 photonEnergy0/keV << " keV (initial)" << G
647 }
648 }
649
650 //....oooOO0OOooo........oooOO0OOooo........oo
651
652 G4double G4PenelopeComptonModel::DifferentialC
653 G4PenelopeOscillator* osc)
654 {
655 //
656 // Penelope model v2008. Single differential
657 // for photon Compton scattering by
658 // electrons in the given atomic oscillator,
659 // scattering photon. This is in units of pi
660 //
661 // D. Brusa et al., Nucl. Instrum. Meth. A 3
662 // The parametrization includes the J(p) dis
663 // that are tabulated from Hartree-Fock calc
664 // from F. Biggs et al., At. Data Nucl. Data
665 //
666 G4double ionEnergy = osc->GetIonisationEnerg
667 G4double harFunc = osc->GetHartreeFactor();
668
669 static const G4double k2 = std::sqrt(2.);
670 static const G4double k1 = 1./k2;
671
672 if (energy < ionEnergy)
673 return 0;
674
675 //energy of the Compton line
676 G4double cdt1 = 1.0-cosTheta;
677 G4double EOEC = 1.0+(energy/electron_mass_c2
678 G4double ECOE = 1.0/EOEC;
679
680 //Incoherent scattering function (analytical
681 G4double aux = energy*(energy-ionEnergy)*cdt
682 G4double Pzimax =
683 (aux - electron_mass_c2*ionEnergy)/(electr
684 G4double sia = 0.0;
685 G4double x = harFunc*Pzimax;
686 if (x > 0)
687 sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x
688 else
689 sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x));
690
691 //1st order correction, integral of Pz times
692 //Calculated approximately using a free-elec
693 G4double pf = 3.0/(4.0*harFunc);
694 if (std::fabs(Pzimax) < pf)
695 {
696 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*
697 G4double p2 = Pzimax*Pzimax;
698 G4double dspz = std::sqrt(QCOE2)*
699 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
700 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
701 sia += std::max(dspz,-1.0*sia);
702 }
703
704 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosThe
705
706 //Differential cross section (per electron,
707 G4double diffCS = ECOE*ECOE*XKN*sia;
708
709 return diffCS;
710 }
711
712 //....oooOO0OOooo........oooOO0OOooo........oo
713
714 G4double G4PenelopeComptonModel::OscillatorTot
715 {
716 //Total cross section (integrated) for the g
717 //pi*classic_electr_radius^2
718
719 //Integrate differential cross section for e
720 G4double stre = osc->GetOscillatorStrength()
721
722 // here one uses the using the 20-point
723 // Gauss quadrature method with an adaptive
724 const G4int npoints=10;
725 const G4int ncallsmax=20000;
726 const G4int nst=256;
727 static const G4double Abscissas[10] = {7.652
728 5.1086700195082710e-01,6.36053680726
729 8.3911697182221882e-01,9.12234428251
730 9.9312859918509492e-01};
731 static const G4double Weights[10] = {1.52753
732 1.3168863844917663e-01,1.1819453196151
733 8.3276741576704749e-02,6.2672048334109
734 1.7614007139152118e-02};
735
736 G4double MaxError = 1e-5;
737 //Error control
738 G4double Ctol = std::min(std::max(MaxError,1
739 G4double Ptol = 0.01*Ctol;
740 G4double Err=1e35;
741
742 //Gauss integration from -1 to 1
743 G4double LowPoint = -1.0;
744 G4double HighPoint = 1.0;
745
746 G4double h=HighPoint-LowPoint;
747 G4double sumga=0.0;
748 G4double a=0.5*(HighPoint-LowPoint);
749 G4double b=0.5*(HighPoint+LowPoint);
750 G4double c=a*Abscissas[0];
751 G4double d= Weights[0]*
752 (DifferentialCrossSection(b+c,energy,osc)+
753 for (G4int i=2;i<=npoints;i++)
754 {
755 c=a*Abscissas[i-1];
756 d += Weights[i-1]*
757 (DifferentialCrossSection(b+c,energy,osc)+Di
758 }
759 G4int icall = 2*npoints;
760 G4int LH=1;
761 G4double S[nst],x[nst],sn[nst],xrn[nst];
762 S[0]=d*a;
763 x[0]=LowPoint;
764
765 G4bool loopAgain = true;
766
767 //Adaptive bipartition scheme
768 do{
769 G4double h0=h;
770 h=0.5*h; //bipartition
771 G4double sumr=0;
772 G4int LHN=0;
773 G4double si,xa,xb,xc;
774 for (G4int i=1;i<=LH;i++){
775 si=S[i-1];
776 xa=x[i-1];
777 xb=xa+h;
778 xc=xa+h0;
779 a=0.5*(xb-xa);
780 b=0.5*(xb+xa);
781 c=a*Abscissas[0];
782 G4double dLocal = Weights[0]*
783 (DifferentialCrossSection(b+c,energy,osc)+Di
784
785 for (G4int j=1;j<npoints;j++)
786 {
787 c=a*Abscissas[j];
788 dLocal += Weights[j]*
789 (DifferentialCrossSection(b+c,energy,osc
790 }
791 G4double s1=dLocal*a;
792 a=0.5*(xc-xb);
793 b=0.5*(xc+xb);
794 c=a*Abscissas[0];
795 dLocal=Weights[0]*
796 (DifferentialCrossSection(b+c,energy,osc)+Di
797
798 for (G4int j=1;j<npoints;j++)
799 {
800 c=a*Abscissas[j];
801 dLocal += Weights[j]*
802 (DifferentialCrossSection(b+c,energy,osc
803 }
804 G4double s2=dLocal*a;
805 icall=icall+4*npoints;
806 G4double s12=s1+s2;
807 if (std::abs(s12-si)<std::max(Ptol*std::
808 sumga += s12;
809 else
810 {
811 sumr += s12;
812 LHN += 2;
813 sn[LHN-1]=s2;
814 xrn[LHN-1]=xb;
815 sn[LHN-2]=s1;
816 xrn[LHN-2]=xa;
817 }
818
819 if (icall>ncallsmax || LHN>nst)
820 {
821 G4cout << "G4PenelopeComptonModel: " << G4
822 G4cout << "LowPoint: " << LowPoint << ", H
823 G4cout << "Tolerance: " << MaxError << G4e
824 G4cout << "Calls: " << icall << ", Integra
825 G4cout << "Number of open subintervals: "
826 G4cout << "WARNING: the required accuracy
827 loopAgain = false;
828 }
829 }
830 Err=std::abs(sumr)/std::max(std::abs(sumr+
831 if (Err < Ctol || LHN == 0)
832 loopAgain = false; //end of cycle
833 LH=LHN;
834 for (G4int i=0;i<LH;i++)
835 {
836 S[i]=sn[i];
837 x[i]=xrn[i];
838 }
839 }while(Ctol < 1.0 && loopAgain);
840
841 G4double xs = stre*sumga;
842
843 return xs;
844 }
845
846 //....oooOO0OOooo........oooOO0OOooo........oo
847
848 G4double G4PenelopeComptonModel::KleinNishinaC
849 const G4Material* material)
850 {
851 // use Klein-Nishina formula
852 // total cross section in units of pi*classi
853 G4double cs = 0;
854
855 G4double ek =energy/electron_mass_c2;
856 G4double eks = ek*ek;
857 G4double ek2 = 1.0+ek+ek;
858 G4double ek1 = eks-ek2-1.0;
859
860 G4double t0 = 1.0/ek2;
861 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*G4Lo
862
863 G4PenelopeOscillatorTable* theTable = fOscMa
864
865 for (size_t i=0;i<theTable->size();i++)
866 {
867 G4PenelopeOscillator* theOsc = (*theTabl
868 G4double ionEnergy = theOsc->GetIonisati
869 G4double tau=(energy-ionEnergy)/energy;
870 if (tau > t0)
871 {
872 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1
873 G4double stre = theOsc->GetOscillatorStren
874
875 cs += stre*(csu-csl);
876 }
877 }
878 cs /= (ek*eks);
879
880 return cs;
881
882 }
883
884 //....oooOO0OOooo........oooOO0OOooo........oo
885
886 void G4PenelopeComptonModel::SetParticle(const
887 {
888 if(!fParticle) {
889 fParticle = p;
890 }
891 }
892