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