Geant4-11
G4PenelopeComptonModel.cc
Go to the documentation of this file.
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// Author: Luciano Pandola
28//
29// History:
30// --------
31// 15 Feb 2010 L Pandola Implementation
32// 18 Mar 2010 L Pandola Removed GetAtomsPerMolecule(), now demanded
33// to G4PenelopeOscillatorManager
34// 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is
35// active.
36// Make sure that fluorescence/Auger is generated only
37// if above threshold
38// 24 May 2011 L Pandola Renamed (make v2008 as default Penelope)
39// 10 Jun 2011 L Pandola Migrate atomic deexcitation interface
40// 09 Oct 2013 L Pandola Migration to MT
41//
44#include "G4SystemOfUnits.hh"
47#include "G4DynamicParticle.hh"
48#include "G4VEMDataSet.hh"
49#include "G4PhysicsTable.hh"
50#include "G4PhysicsLogVector.hh"
52#include "G4AtomicShell.hh"
53#include "G4Gamma.hh"
54#include "G4Electron.hh"
57#include "G4LossTableManager.hh"
58#include "G4Exp.hh"
59
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61
63 const G4String& nam)
64 :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
65 fAtomDeexcitation(nullptr),
66 fOscManager(nullptr),fIsInitialised(false)
67{
71 //
73
74 if (part)
75 SetParticle(part);
76
78 // Verbosity scale:
79 // 0 = nothing
80 // 1 = warning for energy non-conservation
81 // 2 = details of energy budget
82 // 3 = calculation of cross sections, file openings, sampling of atoms
83 // 4 = entering in methods
84
85 //Mark this model as "applicable" for atomic deexcitation
87
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
94{;}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97
99 const G4DataVector&)
100{
101 if (fVerboseLevel > 3)
102 G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
103
105 //Issue warning if the AtomicDeexcitation has not been declared
107 {
108 G4cout << G4endl;
109 G4cout << "WARNING from G4PenelopeComptonModel " << G4endl;
110 G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
111 G4cout << "any fluorescence/Auger emission." << G4endl;
112 G4cout << "Please make sure this is intended" << G4endl;
113 }
114
115 SetParticle(part);
116
117 if (IsMaster() && part == fParticle)
118 {
119
120 if (fVerboseLevel > 0)
121 {
122 G4cout << "Penelope Compton model v2008 is initialized " << G4endl
123 << "Energy range: "
124 << LowEnergyLimit() / keV << " keV - "
125 << HighEnergyLimit() / GeV << " GeV";
126 }
127 //Issue a warning, if the model is going to be used down to a
128 //energy which is outside the validity of the model itself
130 {
132 ed << "Using the Penelope Compton model outside its intrinsic validity range. "
133 << G4endl;
134 ed << "-> LowEnergyLimit() in process = " << LowEnergyLimit()/keV << "keV " << G4endl;
135 ed << "-> Instrinsic low-energy limit = " << fIntrinsicLowEnergyLimit/keV << "keV "
136 << G4endl;
137 ed << "Result of the simulation have to be taken with care" << G4endl;
138 G4Exception("G4PenelopeComptonModel::Initialise()",
139 "em2100",JustWarning,ed);
140 }
141 }
142
143 if(fIsInitialised) return;
145 fIsInitialised = true;
146
147}
148
149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150
152 G4VEmModel *masterModel)
153{
154 if (fVerboseLevel > 3)
155 G4cout << "Calling G4PenelopeComptonModel::InitialiseLocal()" << G4endl;
156 //
157 //Check that particle matches: one might have multiple master models (e.g.
158 //for e+ and e-).
159 //
160 if (part == fParticle)
161 {
162 //Get the const table pointers from the master to the workers
163 const G4PenelopeComptonModel* theModel =
164 static_cast<G4PenelopeComptonModel*> (masterModel);
165
166 //Same verbosity for all workers, as the master
167 fVerboseLevel = theModel->fVerboseLevel;
168 }
169 return;
170}
171
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174
176 const G4ParticleDefinition* p,
178 G4double,
179 G4double)
180{
181 // Penelope model v2008 to calculate the Compton scattering cross section:
182 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
183 //
184 // The cross section for Compton scattering is calculated according to the Klein-Nishina
185 // formula for energy > 5 MeV.
186 // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
187 // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
188 // The parametrization includes the J(p)
189 // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations
190 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
191 //
192 if (fVerboseLevel > 3)
193 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeComptonModel" << G4endl;
195
196 G4double cs = 0;
197 //Force null cross-section if below the low-energy edge of the table
198 if (energy < LowEnergyLimit())
199 return cs;
200
201 //Retrieve the oscillator table for this material
203
204 if (energy < 5*MeV) //explicit calculation for E < 5 MeV
205 {
206 size_t numberOfOscillators = theTable->size();
207 for (size_t i=0;i<numberOfOscillators;i++)
208 {
209 G4PenelopeOscillator* theOsc = (*theTable)[i];
210 //sum contributions coming from each oscillator
212 }
213 }
214 else //use Klein-Nishina for E>5 MeV
216
217 //cross sections are in units of pi*classic_electr_radius^2
219
220 //Now, cs is the cross section *per molecule*, let's calculate the
221 //cross section per volume
222 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
224
225 if (fVerboseLevel > 3)
226 G4cout << "Material " << material->GetName() << " has " << atPerMol <<
227 "atoms per molecule" << G4endl;
228
229 G4double moleculeDensity = 0.;
230
231 if (atPerMol)
232 moleculeDensity = atomDensity/atPerMol;
233
234 G4double csvolume = cs*moleculeDensity;
235
236 if (fVerboseLevel > 2)
237 G4cout << "Compton mean free path at " << energy/keV << " keV for material " <<
238 material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl;
239 return csvolume;
240}
241
242//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243
244//This is a dummy method. Never inkoved by the tracking, it just issues
245//a warning if one tries to get Cross Sections per Atom via the
246//G4EmCalculator.
248 G4double,
249 G4double,
250 G4double,
251 G4double,
252 G4double)
253{
254 G4cout << "*** G4PenelopeComptonModel -- WARNING ***" << G4endl;
255 G4cout << "Penelope Compton model v2008 does not calculate cross section _per atom_ " << G4endl;
256 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
257 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
258 return 0;
259}
260
261//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262
263void G4PenelopeComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
264 const G4MaterialCutsCouple* couple,
265 const G4DynamicParticle* aDynamicGamma,
266 G4double,
267 G4double)
268{
269 // Penelope model v2008 to sample the Compton scattering final state.
270 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
271 // The model determines also the original shell from which the electron is expelled,
272 // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
273 //
274 // The final state for Compton scattering is calculated according to the Klein-Nishina
275 // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and
276 // one can assume that the target electron is at rest.
277 // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
278 // to sample the scattering angle and the energy of the emerging electron, which is
279 // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is
280 // used to sample cos(theta). The efficiency increases monotonically with photon energy and is
281 // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV,
282 // respectively.
283 // The parametrization includes the J(p) distribution profiles for the atomic shells, that are
284 // tabulated
285 // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201.
286 // Doppler broadening is included.
287 //
288
289 if (fVerboseLevel > 3)
290 G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
291
292 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
293
294 // do nothing below the threshold
295 // should never get here because the XS is zero below the limit
296 if(photonEnergy0 < LowEnergyLimit())
297 return;
298
299 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
300 const G4Material* material = couple->GetMaterial();
301
303
304 const G4int nmax = 64;
305 G4double rn[nmax]={0.0};
306 G4double pac[nmax]={0.0};
307
308 G4double S=0.0;
309 G4double epsilon = 0.0;
310 G4double cosTheta = 1.0;
311 G4double hartreeFunc = 0.0;
312 G4double oscStren = 0.0;
313 size_t numberOfOscillators = theTable->size();
314 size_t targetOscillator = 0;
315 G4double ionEnergy = 0.0*eV;
316
317 G4double ek = photonEnergy0/electron_mass_c2;
318 G4double ek2 = 2.*ek+1.0;
319 G4double eks = ek*ek;
320 G4double ek1 = eks-ek2-1.0;
321
322 G4double taumin = 1.0/ek2;
323 G4double a1 = G4Log(ek2);
324 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
325
326 G4double TST = 0;
327 G4double tau = 0.;
328
329 //If the incoming photon is above 5 MeV, the quicker approach based on the
330 //pure Klein-Nishina formula is used
331 if (photonEnergy0 > 5*MeV)
332 {
333 do{
334 do{
335 if ((a2*G4UniformRand()) < a1)
336 tau = std::pow(taumin,G4UniformRand());
337 else
338 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
339 //rejection function
340 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
341 }while (G4UniformRand()> TST);
342 epsilon=tau;
343 cosTheta = 1.0 - (1.0-tau)/(ek*tau);
344
345 //Target shell electrons
347 targetOscillator = numberOfOscillators-1; //last level
348 S=0.0;
349 G4bool levelFound = false;
350 for (size_t j=0;j<numberOfOscillators && !levelFound; j++)
351 {
352 S += (*theTable)[j]->GetOscillatorStrength();
353 if (S > TST)
354 {
355 targetOscillator = j;
356 levelFound = true;
357 }
358 }
359 //check whether the level is valid
360 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
361 }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
362 }
363 else //photonEnergy0 < 5 MeV
364 {
365 //Incoherent scattering function for theta=PI
366 G4double s0=0.0;
367 G4double pzomc=0.0;
368 G4double rni=0.0;
369 G4double aux=0.0;
370 for (size_t i=0;i<numberOfOscillators;i++)
371 {
372 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
373 if (photonEnergy0 > ionEnergy)
374 {
375 G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
376 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
377 oscStren = (*theTable)[i]->GetOscillatorStrength();
378 pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
379 (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
380 if (pzomc > 0)
381 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
382 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
383 else
384 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
385 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
386 s0 += oscStren*rni;
387 }
388 }
389 //Sampling tau
390 G4double cdt1 = 0.;
391 do
392 {
393 if ((G4UniformRand()*a2) < a1)
394 tau = std::pow(taumin,G4UniformRand());
395 else
396 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
397 cdt1 = (1.0-tau)/(ek*tau);
398 //Incoherent scattering function
399 S = 0.;
400 for (size_t i=0;i<numberOfOscillators;i++)
401 {
402 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
403 if (photonEnergy0 > ionEnergy) //sum only on excitable levels
404 {
405 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
406 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
407 oscStren = (*theTable)[i]->GetOscillatorStrength();
408 pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
409 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
410 if (pzomc > 0)
411 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
412 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
413 else
414 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
415 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
416 S += oscStren*rn[i];
417 pac[i] = S;
418 }
419 else
420 pac[i] = S-1e-6;
421 }
422 //Rejection function
423 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
424 }while ((G4UniformRand()*s0) > TST);
425
426 cosTheta = 1.0 - cdt1;
427 G4double fpzmax=0.0,fpz=0.0;
428 G4double A=0.0;
429 //Target electron shell
430 do
431 {
432 do
433 {
434 TST = S*G4UniformRand();
435 targetOscillator = numberOfOscillators-1; //last level
436 G4bool levelFound = false;
437 for (size_t i=0;i<numberOfOscillators && !levelFound;i++)
438 {
439 if (pac[i]>TST)
440 {
441 targetOscillator = i;
442 levelFound = true;
443 }
444 }
445 A = G4UniformRand()*rn[targetOscillator];
446 hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
447 oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
448 if (A < 0.5)
449 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-G4Log(2.0*A)))/
450 (std::sqrt(2.0)*hartreeFunc);
451 else
452 pzomc = (std::sqrt(0.5-G4Log(2.0-2.0*A))-std::sqrt(0.5))/
453 (std::sqrt(2.0)*hartreeFunc);
454 } while (pzomc < -1);
455
456 // F(EP) rejection
457 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
458 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
459 if (AF > 0)
460 fpzmax = 1.0+AF*0.2;
461 else
462 fpzmax = 1.0-AF*0.2;
463 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
464 }while ((fpzmax*G4UniformRand())>fpz);
465
466 //Energy of the scattered photon
467 G4double T = pzomc*pzomc;
468 G4double b1 = 1.0-T*tau*tau;
469 G4double b2 = 1.0-T*tau*cosTheta;
470 if (pzomc > 0.0)
471 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
472 else
473 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
474 } //energy < 5 MeV
475
476 //Ok, the kinematics has been calculated.
477 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
478 G4double phi = twopi * G4UniformRand() ;
479 G4double dirx = sinTheta * std::cos(phi);
480 G4double diry = sinTheta * std::sin(phi);
481 G4double dirz = cosTheta ;
482
483 // Update G4VParticleChange for the scattered photon
484 G4ThreeVector photonDirection1(dirx,diry,dirz);
485 photonDirection1.rotateUz(photonDirection0);
486 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
487
488 G4double photonEnergy1 = epsilon * photonEnergy0;
489
490 if (photonEnergy1 > 0.)
492 else
493 {
496 }
497
498 //Create scattered electron
499 G4double diffEnergy = photonEnergy0*(1-epsilon);
500 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
501
502 G4double Q2 =
503 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
504 G4double cosThetaE = 0.; //scattering angle for the electron
505
506 if (Q2 > 1.0e-12)
507 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
508 else
509 cosThetaE = 1.0;
510 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
511
512 //Now, try to handle fluorescence
513 //Notice: merged levels are indicated with Z=0 and flag=30
514 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
515 G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
516
517 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
519 const G4AtomicShell* shell = 0;
520
521 //Real level
522 if (Z > 0 && shFlag<30)
523 {
524 shell = fTransitionManager->Shell(Z,shFlag-1);
525 bindingEnergy = shell->BindingEnergy();
526 }
527
528 G4double ionEnergyInPenelopeDatabase = ionEnergy;
529 //protection against energy non-conservation
530 ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
531
532 //subtract the excitation energy. If not emitted by fluorescence
533 //the ionization energy is deposited as local energy deposition
534 G4double eKineticEnergy = diffEnergy - ionEnergy;
535 G4double localEnergyDeposit = ionEnergy;
536 G4double energyInFluorescence = 0.; //testing purposes only
537 G4double energyInAuger = 0; //testing purposes
538
539 if (eKineticEnergy < 0)
540 {
541 //It means that there was some problem/mismatch between the two databases.
542 //Try to make it work
543 //In this case available Energy (diffEnergy) < ionEnergy
544 //Full residual energy is deposited locally
545 localEnergyDeposit = diffEnergy;
546 eKineticEnergy = 0.0;
547 }
548
549 //the local energy deposit is what remains: part of this may be spent for fluorescence.
550 //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected
551 //Now, take care of fluorescence, if required
552 if (fAtomDeexcitation && shell)
553 {
554 G4int index = couple->GetIndex();
556 {
557 size_t nBefore = fvect->size();
558 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
559 size_t nAfter = fvect->size();
560
561 if (nAfter > nBefore) //actual production of fluorescence
562 {
563 for (size_t j=nBefore;j<nAfter;j++) //loop on products
564 {
565 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
566 if (itsEnergy < localEnergyDeposit) // valid secondary, generate it
567 {
568 localEnergyDeposit -= itsEnergy;
569 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
570 energyInFluorescence += itsEnergy;
571 else if (((*fvect)[j])->GetParticleDefinition() ==
573 energyInAuger += itsEnergy;
574 }
575 else //invalid secondary: takes more than the available energy: delete it
576 {
577 delete (*fvect)[j];
578 (*fvect)[j] = nullptr;
579 }
580 }
581 }
582
583 }
584 }
585
586 //Always produce explicitly the electron
588
589 G4double xEl = sinThetaE * std::cos(phi+pi);
590 G4double yEl = sinThetaE * std::sin(phi+pi);
591 G4double zEl = cosThetaE;
592 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
593 eDirection.rotateUz(photonDirection0);
595 eDirection,eKineticEnergy) ;
596 fvect->push_back(electron);
597
598 if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning)
599 {
600 G4Exception("G4PenelopeComptonModel::SampleSecondaries()",
601 "em2099",JustWarning,"WARNING: Negative local energy deposit");
602 localEnergyDeposit=0.;
603 }
604 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
605
606 G4double electronEnergy = 0.;
607 if (electron)
608 electronEnergy = eKineticEnergy;
609 if (fVerboseLevel > 1)
610 {
611 G4cout << "-----------------------------------------------------------" << G4endl;
612 G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
613 G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
614 G4cout << "-----------------------------------------------------------" << G4endl;
615 G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
616 G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
617 if (energyInFluorescence)
618 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
619 if (energyInAuger)
620 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
621 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
622 G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
623 localEnergyDeposit+energyInAuger)/keV <<
624 " keV" << G4endl;
625 G4cout << "-----------------------------------------------------------" << G4endl;
626 }
627 if (fVerboseLevel > 0)
628 {
629 G4double energyDiff = std::fabs(photonEnergy1+
630 electronEnergy+energyInFluorescence+
631 localEnergyDeposit+energyInAuger-photonEnergy0);
632 if (energyDiff > 0.05*keV)
633 G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " <<
634 (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/keV <<
635 " keV (final) vs. " <<
636 photonEnergy0/keV << " keV (initial)" << G4endl;
637 }
638}
639
640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
641
644{
645 //
646 // Penelope model v2008. Single differential cross section *per electron*
647 // for photon Compton scattering by
648 // electrons in the given atomic oscillator, differential in the direction of the
649 // scattering photon. This is in units of pi*classic_electr_radius**2
650 //
651 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
652 // The parametrization includes the J(p) distribution profiles for the atomic shells,
653 // that are tabulated from Hartree-Fock calculations
654 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
655 //
656 G4double ionEnergy = osc->GetIonisationEnergy();
657 G4double harFunc = osc->GetHartreeFactor();
658
659 static const G4double k2 = std::sqrt(2.);
660 static const G4double k1 = 1./k2;
661
662 if (energy < ionEnergy)
663 return 0;
664
665 //energy of the Compton line
666 G4double cdt1 = 1.0-cosTheta;
667 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
668 G4double ECOE = 1.0/EOEC;
669
670 //Incoherent scattering function (analytical profile)
671 G4double aux = energy*(energy-ionEnergy)*cdt1;
672 G4double Pzimax =
673 (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
674 G4double sia = 0.0;
675 G4double x = harFunc*Pzimax;
676 if (x > 0)
677 sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x));
678 else
679 sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x));
680
681 //1st order correction, integral of Pz times the Compton profile.
682 //Calculated approximately using a free-electron gas profile
683 G4double pf = 3.0/(4.0*harFunc);
684 if (std::fabs(Pzimax) < pf)
685 {
686 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta;
687 G4double p2 = Pzimax*Pzimax;
688 G4double dspz = std::sqrt(QCOE2)*
689 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
690 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
691 sia += std::max(dspz,-1.0*sia);
692 }
693
694 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta;
695
696 //Differential cross section (per electron, in units of pi*classic_electr_radius**2)
697 G4double diffCS = ECOE*ECOE*XKN*sia;
698
699 return diffCS;
700}
701
702//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
703
705{
706 //Total cross section (integrated) for the given oscillator in units of
707 //pi*classic_electr_radius^2
708
709 //Integrate differential cross section for each oscillator
710 G4double stre = osc->GetOscillatorStrength();
711
712 // here one uses the using the 20-point
713 // Gauss quadrature method with an adaptive bipartition scheme
714 const G4int npoints=10;
715 const G4int ncallsmax=20000;
716 const G4int nst=256;
717 static const G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01,
718 5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01,
719 8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01,
720 9.9312859918509492e-01};
721 static const G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01,
722 1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01,
723 8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02,
724 1.7614007139152118e-02};
725
726 G4double MaxError = 1e-5;
727 //Error control
728 G4double Ctol = std::min(std::max(MaxError,1e-13),1e-02);
729 G4double Ptol = 0.01*Ctol;
730 G4double Err=1e35;
731
732 //Gauss integration from -1 to 1
733 G4double LowPoint = -1.0;
734 G4double HighPoint = 1.0;
735
736 G4double h=HighPoint-LowPoint;
737 G4double sumga=0.0;
738 G4double a=0.5*(HighPoint-LowPoint);
739 G4double b=0.5*(HighPoint+LowPoint);
740 G4double c=a*Abscissas[0];
741 G4double d= Weights[0]*
743 for (G4int i=2;i<=npoints;i++)
744 {
745 c=a*Abscissas[i-1];
746 d += Weights[i-1]*
748 }
749 G4int icall = 2*npoints;
750 G4int LH=1;
751 G4double S[nst],x[nst],sn[nst],xrn[nst];
752 S[0]=d*a;
753 x[0]=LowPoint;
754
755 G4bool loopAgain = true;
756
757 //Adaptive bipartition scheme
758 do{
759 G4double h0=h;
760 h=0.5*h; //bipartition
761 G4double sumr=0;
762 G4int LHN=0;
763 G4double si,xa,xb,xc;
764 for (G4int i=1;i<=LH;i++){
765 si=S[i-1];
766 xa=x[i-1];
767 xb=xa+h;
768 xc=xa+h0;
769 a=0.5*(xb-xa);
770 b=0.5*(xb+xa);
771 c=a*Abscissas[0];
772 G4double dLocal = Weights[0]*
774
775 for (G4int j=1;j<npoints;j++)
776 {
777 c=a*Abscissas[j];
778 dLocal += Weights[j]*
780 }
781 G4double s1=dLocal*a;
782 a=0.5*(xc-xb);
783 b=0.5*(xc+xb);
784 c=a*Abscissas[0];
785 dLocal=Weights[0]*
787
788 for (G4int j=1;j<npoints;j++)
789 {
790 c=a*Abscissas[j];
791 dLocal += Weights[j]*
793 }
794 G4double s2=dLocal*a;
795 icall=icall+4*npoints;
796 G4double s12=s1+s2;
797 if (std::abs(s12-si)<std::max(Ptol*std::abs(s12),1e-35))
798 sumga += s12;
799 else
800 {
801 sumr += s12;
802 LHN += 2;
803 sn[LHN-1]=s2;
804 xrn[LHN-1]=xb;
805 sn[LHN-2]=s1;
806 xrn[LHN-2]=xa;
807 }
808
809 if (icall>ncallsmax || LHN>nst)
810 {
811 G4cout << "G4PenelopeComptonModel: " << G4endl;
812 G4cout << "LowPoint: " << LowPoint << ", High Point: " << HighPoint << G4endl;
813 G4cout << "Tolerance: " << MaxError << G4endl;
814 G4cout << "Calls: " << icall << ", Integral: " << sumga << ", Error: " << Err << G4endl;
815 G4cout << "Number of open subintervals: " << LHN << G4endl;
816 G4cout << "WARNING: the required accuracy has not been attained" << G4endl;
817 loopAgain = false;
818 }
819 }
820 Err=std::abs(sumr)/std::max(std::abs(sumr+sumga),1e-35);
821 if (Err < Ctol || LHN == 0)
822 loopAgain = false; //end of cycle
823 LH=LHN;
824 for (G4int i=0;i<LH;i++)
825 {
826 S[i]=sn[i];
827 x[i]=xrn[i];
828 }
829 }while(Ctol < 1.0 && loopAgain);
830
831 G4double xs = stre*sumga;
832
833 return xs;
834}
835
836//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
837
839 const G4Material* material)
840{
841 // use Klein-Nishina formula
842 // total cross section in units of pi*classic_electr_radius^2
843 G4double cs = 0;
844
846 G4double eks = ek*ek;
847 G4double ek2 = 1.0+ek+ek;
848 G4double ek1 = eks-ek2-1.0;
849
850 G4double t0 = 1.0/ek2;
851 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*G4Log(t0)-(1.0/t0);
852
854
855 for (size_t i=0;i<theTable->size();i++)
856 {
857 G4PenelopeOscillator* theOsc = (*theTable)[i];
858 G4double ionEnergy = theOsc->GetIonisationEnergy();
859 G4double tau=(energy-ionEnergy)/energy;
860 if (tau > t0)
861 {
862 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*G4Log(tau)-(1.0/tau);
863 G4double stre = theOsc->GetOscillatorStrength();
864
865 cs += stre*(csu-csl);
866 }
867 }
868 cs /= (ek*eks);
869
870 return cs;
871
872}
873
874//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
875
877{
878 if(!fParticle) {
879 fParticle = p;
880 }
881}
G4double epsilon(G4double density, G4double temperature)
G4double S(G4double temp)
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4double G4Log(G4double x)
Definition: G4Log.hh:226
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
static constexpr double twopi
Definition: G4SIunits.hh:56
static constexpr double mm
Definition: G4SIunits.hh:95
static constexpr double keV
Definition: G4SIunits.hh:202
static constexpr double eV
Definition: G4SIunits.hh:201
static constexpr double GeV
Definition: G4SIunits.hh:203
static constexpr double MeV
Definition: G4SIunits.hh:200
static constexpr double pi
Definition: G4SIunits.hh:55
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
const G4double A[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
G4double BindingEnergy() const
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4AtomicTransitionManager * Instance()
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Definition()
Definition: G4Electron.cc:48
static G4Electron * Electron()
Definition: G4Electron.cc:93
static G4Gamma * Definition()
Definition: G4Gamma.cc:48
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void SetParticle(const G4ParticleDefinition *)
const G4AtomicTransitionManager * fTransitionManager
G4ParticleChangeForGamma * fParticleChange
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double KleinNishinaCrossSection(G4double energy, const G4Material *)
G4double OscillatorTotalCrossSection(G4double energy, G4PenelopeOscillator *osc)
const G4ParticleDefinition * fParticle
G4PenelopeComptonModel(const G4ParticleDefinition *p=nullptr, const G4String &processName="PenCompton")
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double) override
G4VAtomDeexcitation * fAtomDeexcitation
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) override
G4PenelopeOscillatorManager * fOscManager
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
G4double DifferentialCrossSection(G4double cdt, G4double energy, G4PenelopeOscillator *osc)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
static G4PenelopeOscillatorManager * GetOscillatorManager()
G4double GetTotalZ(const G4Material *)
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:767
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:123
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:662
G4bool IsMaster() const
Definition: G4VEmModel.hh:746
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:655
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:440
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:823
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double energy(const ThreeVector &p, const G4double m)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double bindingEnergy(G4int A, G4int Z)
string material
Definition: eplot.py:19
float electron_mass_c2
Definition: hepunit.py:273
int classic_electr_radius
Definition: hepunit.py:287