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