source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedPhotoElectricModel.cc @ 1250

Last change on this file since 1250 was 1197, checked in by garnier, 15 years ago

nvx fichiers dans CVS
File size: 17.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26
27#include "G4LivermorePolarizedPhotoElectricModel.hh"
28
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
30
31using namespace std;
32
33//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34
35G4LivermorePolarizedPhotoElectricModel::G4LivermorePolarizedPhotoElectricModel(const G4ParticleDefinition*,
36                                             const G4String& nam)
37  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),crossSectionHandler(0), shellCrossSectionHandler(0)
38{
39  lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
40  highEnergyLimit = 100 * GeV;
41  SetLowEnergyLimit(lowEnergyLimit);
42  SetHighEnergyLimit(highEnergyLimit);
43 
44  verboseLevel= 0;
45  // Verbosity scale:
46  // 0 = nothing
47  // 1 = warning for energy non-conservation
48  // 2 = details of energy budget
49  // 3 = calculation of cross sections, file openings, sampling of atoms
50  // 4 = entering in methods
51
52  G4cout << "Livermore Polarized PhotoElectric is constructed " << G4endl
53         << "Energy range: "
54         << lowEnergyLimit / keV << " keV - "
55         << highEnergyLimit / GeV << " GeV"
56         << G4endl;
57
58}
59
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61
62G4LivermorePolarizedPhotoElectricModel::~G4LivermorePolarizedPhotoElectricModel()
63{ 
64  if (meanFreePathTable)   delete meanFreePathTable;
65  if (crossSectionHandler) delete crossSectionHandler;
66  if (shellCrossSectionHandler) delete shellCrossSectionHandler;
67}
68
69
70
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72
73void G4LivermorePolarizedPhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
74                                       const G4DataVector& cuts)
75{
76  if (verboseLevel > 3)
77    G4cout << "Calling G4LivermorePolarizedPhotoElectricModel::Initialise()" << G4endl;
78
79  if (crossSectionHandler)
80  {
81    crossSectionHandler->Clear();
82    delete crossSectionHandler;
83  }
84
85  if (shellCrossSectionHandler)
86    {
87      shellCrossSectionHandler->Clear();
88      delete shellCrossSectionHandler;
89    }
90
91
92  // Energy limits
93 
94  if (LowEnergyLimit() < lowEnergyLimit)
95    {
96      G4cout << "G4LivermorePolarizedPhotoElectricModel: low energy limit increased from " << 
97        LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
98      SetLowEnergyLimit(lowEnergyLimit);
99    }
100
101  if (HighEnergyLimit() > highEnergyLimit)
102    {
103      G4cout << "G4LivermorePolarizedPhotoElectricModel: high energy limit decreased from " << 
104        HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
105      SetHighEnergyLimit(highEnergyLimit);
106    }
107
108  // Reading of data files - all materials are read
109 
110  crossSectionHandler = new G4CrossSectionHandler;
111  crossSectionHandler->Clear();
112  G4String crossSectionFile = "phot/pe-cs-";
113  crossSectionHandler->LoadData(crossSectionFile);
114
115  meanFreePathTable = 0;
116  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
117
118  shellCrossSectionHandler = new G4CrossSectionHandler();
119  shellCrossSectionHandler->Clear();
120  G4String shellCrossSectionFile = "phot/pe-ss-cs-";
121  shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
122
123
124  //
125  if (verboseLevel > 2) 
126    G4cout << "Loaded cross section files for Livermore Polarized PhotoElectric model" << G4endl;
127
128  InitialiseElementSelectors(particle,cuts);
129
130  G4cout << "Livermore Polarized PhotoElectric model is initialized " << G4endl
131         << "Energy range: "
132         << LowEnergyLimit() / keV << " keV - "
133         << HighEnergyLimit() / GeV << " GeV"
134         << G4endl;
135
136  //
137   
138  if(isInitialised) return;
139
140  if(pParticleChange)
141    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
142  else
143    fParticleChange = new G4ParticleChangeForGamma();
144   
145  isInitialised = true;
146}
147
148//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149
150G4double G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom(
151                                       const G4ParticleDefinition*,
152                                             G4double GammaEnergy,
153                                             G4double Z, G4double,
154                                             G4double, G4double)
155{
156  if (verboseLevel > 3)
157    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedPhotoElectricModel" << G4endl;
158
159  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
160  return cs;
161}
162
163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164
165void G4LivermorePolarizedPhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
166                                                               const G4MaterialCutsCouple* couple,
167                                                               const G4DynamicParticle* aDynamicGamma,
168                                                               G4double,
169                                                               G4double)
170{
171
172  // Fluorescence generated according to:
173  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
174  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
175  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
176
177  if (verboseLevel > 3)
178    G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricModel" << G4endl;
179
180  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
181  // Within energy limit?
182
183  if(photonEnergy <= lowEnergyLimit)
184    {
185      fParticleChange->ProposeTrackStatus(fStopAndKill);
186      fParticleChange->SetProposedKineticEnergy(0.);
187      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
188      return;
189    }
190
191
192  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
193
194  // Protection: a polarisation parallel to the
195  // direction causes problems;
196  // in that case find a random polarization
197
198  G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
199
200  // Make sure that the polarization vector is perpendicular to the
201  // gamma direction. If not
202
203  if(!(gammaPolarization0.isOrthogonal(photonDirection, 1e-6))||(gammaPolarization0.mag()==0))
204    { // only for testing now
205      gammaPolarization0 = GetRandomPolarization(photonDirection);
206    }
207  else
208    {
209      if ( gammaPolarization0.howOrthogonal(photonDirection) != 0)
210        {
211          gammaPolarization0 = GetPerpendicularPolarization(photonDirection, gammaPolarization0);
212        }
213    }
214
215  // End of Protection
216
217  //  G4double E0_m = photonEnergy / electron_mass_c2 ;
218
219  // Select randomly one element in the current material
220
221  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
222
223  // Select the ionised shell in the current atom according to shell cross sections
224
225  size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
226
227  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
228  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
229  G4double bindingEnergy = shell->BindingEnergy();
230  G4int shellId = shell->ShellId();
231
232  // Create lists of pointers to DynamicParticles (photons and electrons)
233  // (Is the electron vector necessary? To be checked)
234  std::vector<G4DynamicParticle*>* photonVector = 0;
235  std::vector<G4DynamicParticle*> electronVector;
236
237  G4double energyDeposit = 0.0;
238
239  // Primary outgoing  electron
240
241  G4double eKineticEnergy = photonEnergy - bindingEnergy;
242
243
244  if (eKineticEnergy > 0.)
245    { 
246
247      G4double costheta = SetCosTheta(eKineticEnergy);
248      G4double sintheta = sqrt(1. - costheta*costheta);
249      G4double phi = SetPhi(photonEnergy,eKineticEnergy,costheta);
250      G4double dirX = sintheta*cos(phi);
251      G4double dirY = sintheta*sin(phi);
252      G4double dirZ = costheta;
253      G4ThreeVector electronDirection(dirX, dirY, dirZ);
254      SystemOfRefChange(photonDirection, electronDirection, gammaPolarization0);
255      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
256                                                           electronDirection, 
257                                                           eKineticEnergy);
258      electronVector.push_back(electron);
259    }
260  else
261    {
262      bindingEnergy = photonEnergy;
263    }
264 
265
266  G4int nElectrons = electronVector.size();
267  size_t nTotPhotons = 0;
268  G4int nPhotons=0;
269
270  const G4ProductionCutsTable* theCoupleTable=
271    G4ProductionCutsTable::GetProductionCutsTable();
272  size_t index = couple->GetIndex();
273  G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
274  cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
275
276  G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
277  cute = std::min(cutForLowEnergySecondaryPhotons,cute);
278
279  G4DynamicParticle* aPhoton;
280
281  // Generation of fluorescence
282  // Data in EADL are available only for Z > 5
283  // Protection to avoid generating photons in the unphysical case of
284  // shell binding energy > photon energy
285  if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
286    {
287      photonVector = deexcitationManager.GenerateParticles(Z,shellId);
288      nTotPhotons = photonVector->size();
289      for (size_t k=0; k<nTotPhotons; k++)
290        {
291          aPhoton = (*photonVector)[k];
292          if (aPhoton)
293            {
294              G4double itsCut = cutg;
295              if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
296
297              G4double itsEnergy = aPhoton->GetKineticEnergy();
298
299              if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
300                {
301                  nPhotons++;
302                  // Local energy deposit is given as the sum of the
303                  // energies of incident photons minus the energies
304                  // of the outcoming fluorescence photons
305                  bindingEnergy -= itsEnergy;
306
307                }
308              else
309                {
310                  delete aPhoton;
311                  (*photonVector)[k] = 0;
312                }
313            }
314        } 
315    }     
316  energyDeposit += bindingEnergy;
317  // Final state
318
319  for (G4int l = 0; l<nElectrons; l++ )
320    {
321      aPhoton = electronVector[l];
322      if(aPhoton) {
323        fvect->push_back(aPhoton);
324      }
325    }
326  for ( size_t ll = 0; ll < nTotPhotons; ll++)
327    {
328      aPhoton = (*photonVector)[ll];
329      if(aPhoton) {
330        fvect->push_back(aPhoton);
331      }
332    }
333
334  delete photonVector;
335
336  if (energyDeposit < 0)
337    {
338      G4cout << "WARNING - "
339             << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
340             << G4endl;
341      energyDeposit = 0;
342    }
343
344  // kill incident photon
345  fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
346  fParticleChange->SetProposedKineticEnergy(0.);
347  fParticleChange->ProposeTrackStatus(fStopAndKill);
348  fParticleChange->ProposeLocalEnergyDeposit(energyDeposit);
349
350}
351
352//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353
354void G4LivermorePolarizedPhotoElectricModel::SetCutForLowEnSecPhotons(G4double cut)
355{
356  cutForLowEnergySecondaryPhotons = cut;
357  deexcitationManager.SetCutForSecondaryPhotons(cut);
358}
359
360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
361
362void G4LivermorePolarizedPhotoElectricModel::SetCutForLowEnSecElectrons(G4double cut)
363{
364  cutForLowEnergySecondaryElectrons = cut;
365  deexcitationManager.SetCutForAugerElectrons(cut);
366}
367
368//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
369
370void G4LivermorePolarizedPhotoElectricModel::ActivateAuger(G4bool val)
371{
372  deexcitationManager.ActivateAugerElectronProduction(val);
373}
374
375//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
376
377
378G4double G4LivermorePolarizedPhotoElectricModel::SetCosTheta(G4double energyE)
379
380{
381  G4double rand1,rand2,onemcost,greject;
382  G4double masarep = 510.99906*keV;
383
384  G4double gamma = 1. + energyE/masarep;
385  G4double gamma2 = gamma*gamma;
386
387  G4double beta = sqrt((gamma2 - 1.)/gamma2);
388
389  G4double alfa = 1./beta - 1.;
390
391  G4double g1 = 0.5*beta*gamma*(gamma-1.)*(gamma-2.);
392
393  G4double alfap2 = alfa+2.;
394
395  G4double grejectmax = 2.*(g1+1./alfa);
396
397  do
398    {
399      rand1 = G4UniformRand();
400      onemcost = 2.*alfa*(2.*rand1 + alfap2 * sqrt(rand1))/
401        (alfap2*alfap2 - 4.*rand1);
402      greject = (2. - onemcost)*(g1+1./(alfa+onemcost));
403      rand2 = G4UniformRand();
404    }
405  while (rand2*grejectmax > greject);
406  G4double cosTheta = 1. - onemcost;
407  return cosTheta;
408} 
409
410
411
412G4double G4LivermorePolarizedPhotoElectricModel::SetPhi(G4double Ph_energy,
413                                                   G4double E_energy,
414                                                   G4double costheta)
415{
416  G4double epsilon = E_energy/electron_mass_c2;
417  G4double k = Ph_energy/electron_mass_c2;
418  G4double gamma = 1. + epsilon;
419  G4double gamma2 = gamma*gamma;
420  G4double beta = sqrt((gamma2 - 1.)/gamma2);
421
422  G4double d = (2./(k*gamma*(1-beta*costheta))-1)*(1/k);
423
424  G4double norm_factor = 1 +2*d;
425
426  G4double rnd1; 
427  G4double rnd2;
428  G4double phi, phiprob;
429
430  do
431    {
432      rnd1 =G4UniformRand();
433      rnd2 =G4UniformRand();
434      phi = rnd1*twopi;
435      phiprob = 1 +2*d*cos(phi)*cos(phi);
436    }
437  while (rnd2*norm_factor > phiprob);
438  return phi;
439}
440
441
442//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
443
444G4ThreeVector G4LivermorePolarizedPhotoElectricModel::SetPerpendicularVector(G4ThreeVector& a)
445{
446  G4double dx = a.x();
447  G4double dy = a.y();
448  G4double dz = a.z();
449  G4double x = dx < 0.0 ? -dx : dx;
450  G4double y = dy < 0.0 ? -dy : dy;
451  G4double z = dz < 0.0 ? -dz : dz;
452  if (x < y) {
453    return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
454  }else{
455    return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
456  }
457}
458
459//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
460
461G4ThreeVector G4LivermorePolarizedPhotoElectricModel::GetRandomPolarization(G4ThreeVector& direction0)
462{
463  G4ThreeVector d0 = direction0.unit();
464  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
465  G4ThreeVector a0 = a1.unit(); // unit vector
466
467  G4double rand1 = G4UniformRand();
468 
469  G4double angle = twopi*rand1; // random polar angle
470  G4ThreeVector b0 = d0.cross(a0); // cross product
471 
472  G4ThreeVector c;
473 
474  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
475  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
476  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
477 
478  G4ThreeVector c0 = c.unit();
479
480  return c0;
481 
482}
483
484//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
485
486G4ThreeVector G4LivermorePolarizedPhotoElectricModel::GetPerpendicularPolarization
487(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
488{
489
490  //
491  // The polarization of a photon is always perpendicular to its momentum direction.
492  // Therefore this function removes those vector component of gammaPolarization, which
493  // points in direction of gammaDirection
494  //
495  // Mathematically we search the projection of the vector a on the plane E, where n is the
496  // plains normal vector.
497  // The basic equation can be found in each geometry book (e.g. Bronstein):
498  // p = a - (a o n)/(n o n)*n
499 
500  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
501}
502
503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504
505
506void G4LivermorePolarizedPhotoElectricModel::SystemOfRefChange
507    (G4ThreeVector& direction0,G4ThreeVector& direction1,
508     G4ThreeVector& polarization0)
509{
510  // direction0 is the original photon direction ---> z
511  // polarization0 is the original photon polarization ---> x
512  // need to specify y axis in the real reference frame ---> y
513  G4ThreeVector Axis_Z0 = direction0.unit();
514  G4ThreeVector Axis_X0 = polarization0.unit();
515  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
516 
517  G4double direction_x = direction1.getX();
518  G4double direction_y = direction1.getY();
519  G4double direction_z = direction1.getZ();
520 
521  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 +  direction_z*Axis_Z0).unit();
522 
523}
524
525
526//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
527
528G4double G4LivermorePolarizedPhotoElectricModel::GetMeanFreePath(const G4Track& track,
529                                                      G4double,
530                                                      G4ForceCondition*)
531{
532
533  const G4DynamicParticle* photon = track.GetDynamicParticle();
534  G4double energy = photon->GetKineticEnergy();
535  G4Material* material = track.GetMaterial();
536  //  size_t materialIndex = material->GetIndex();
537
538  G4double meanFreePath = DBL_MAX;
539
540  //  if (energy > highEnergyLimit)
541  //    meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
542  //  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
543  //  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
544
545  G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
546  if(cross > 0.0) meanFreePath = 1.0/cross;
547
548  return meanFreePath;
549
550
551}
552
553
Note: See TracBrowser for help on using the repository browser.