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