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

Last change on this file since 1285 was 1197, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 17.9 KB
RevLine 
[1197]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.