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

Last change on this file since 1347 was 1347, checked in by garnier, 14 years ago

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