ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermorePolarizedPhotoElectricModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermorePolarizedPhotoElectricModel.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 //
27 //
28 // Author: Sebastien Incerti
29 // 30 October 2008
30 // on base of G4LowEnergyPhotoElectric developed by A.Forti and M.G.Pia
31 //
32 // 22 Oct 2012 A & V Ivanchenko Migration data structure to G4PhysicsVector
33 //
34 
36 #include "G4SystemOfUnits.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4LossTableManager.hh"
39 #include "G4Electron.hh"
40 #include "G4Gamma.hh"
42 #include "G4CrossSectionHandler.hh"
43 #include "G4LPhysicsFreeVector.hh"
44 #include "G4VAtomDeexcitation.hh"
47 #include "G4AtomicShell.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
53 std::vector<G4double>* G4LivermorePolarizedPhotoElectricModel::fParam[] = {nullptr};
59 
60 using namespace std;
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
65  : G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
66  nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
67  fAtomDeexcitation(nullptr)
68 {
69  verboseLevel= 0;
70  // Verbosity scale:
71  // 0 = nothing
72  // 1 = warning for energy non-conservation
73  // 2 = details of energy budget
74  // 3 = calculation of cross sections, file openings, sampling of atoms
75  // 4 = entering in methods
76 
79 
80  // default generator
81  // SetAngularDistribution(new G4SauterGavrilaAngularDistribution());
83 
84  if(verboseLevel>0) {
85  G4cout << "Livermore PhotoElectric is constructed "
86  << " nShellLimit= " << nShellLimit << G4endl;
87  }
88 
89  //Mark this model as "applicable" for atomic deexcitation
90  SetDeexcitationFlag(true);
91  fSandiaCof.resize(4,0.0);
92  fCurrSection = 0.0;
93 
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
99 {
100  if(IsMaster()) {
101  delete fShellCrossSection;
102  for(G4int i=0; i<maxZ; ++i) {
103  delete fParam[i];
104  fParam[i] = 0;
105  delete fCrossSection[i];
106  fCrossSection[i] = 0;
107  delete fCrossSectionLE[i];
108  fCrossSectionLE[i] = 0;
109  }
110  }
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114 
115 void
117  const G4DataVector&)
118 {
119  if (verboseLevel > 2) {
120  G4cout << "Calling G4LivermorePolarizedPhotoElectricModel::Initialise()" << G4endl;
121  }
122 
123  if(IsMaster()) {
124 
125  if(!fWater) {
126  fWater = G4Material::GetMaterial("G4_WATER", false);
127  if(fWater) { fWaterEnergyLimit = 13.6*eV; }
128  }
129 
131 
132  char* path = std::getenv("G4LEDATA");
133 
134  G4ProductionCutsTable* theCoupleTable =
136  G4int numOfCouples = theCoupleTable->GetTableSize();
137 
138  for(G4int i=0; i<numOfCouples; ++i) {
139  const G4MaterialCutsCouple* couple =
140  theCoupleTable->GetMaterialCutsCouple(i);
141  const G4Material* material = couple->GetMaterial();
142  const G4ElementVector* theElementVector = material->GetElementVector();
143  G4int nelm = material->GetNumberOfElements();
144 
145  for (G4int j=0; j<nelm; ++j) {
146  G4int Z = (G4int)(*theElementVector)[j]->GetZ();
147  if(Z < 1) { Z = 1; }
148  else if(Z > maxZ) { Z = maxZ; }
149  if(!fCrossSection[Z]) { ReadData(Z, path); }
150  }
151  }
152  }
153  //
154  if (verboseLevel > 2) {
155  G4cout << "Loaded cross section files for LivermorePolarizedPhotoElectric model"
156  << G4endl;
157  }
158  if(!isInitialised) {
159  isInitialised = true;
161 
163  }
164  fDeexcitationActive = false;
165  if(fAtomDeexcitation) {
167  }
168 
169  if (verboseLevel > 0) {
170  G4cout << "LivermorePolarizedPhotoElectric model is initialized " << G4endl
171  << G4endl;
172  }
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176 
178  const G4Material* material,
179  const G4ParticleDefinition* p,
182 {
183  fCurrSection = 0.0;
184  if(fWater && (material == fWater ||
185  material->GetBaseMaterial() == fWater)) {
186  if(energy <= fWaterEnergyLimit) {
188 
189  G4double energy2 = energy*energy;
190  G4double energy3 = energy*energy2;
191  G4double energy4 = energy2*energy2;
192 
193  fCurrSection = material->GetDensity()*
194  (fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
195  fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4);
196  }
197  }
198  if(0.0 == fCurrSection) {
199  fCurrSection = G4VEmModel::CrossSectionPerVolume(material, p, energy);
200  }
201  return fCurrSection;
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
205 
207  const G4ParticleDefinition*,
209  G4double ZZ, G4double,
211 {
212  if (verboseLevel > 3) {
213  G4cout << "G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom():"
214  << " Z= " << ZZ << " R(keV)= " << energy/keV << G4endl;
215  }
216  G4double cs = 0.0;
217  G4int Z = G4lrint(ZZ);
218  if(Z < 1 || Z >= maxZ) { return cs; }
219 
220  // if element was not initialised
221  // do initialisation safely for MT mode
222  if(!fCrossSection[Z]) {
223  InitialiseForElement(0, Z);
224  if(!fCrossSection[Z]) { return cs; }
225  }
226 
227  G4int idx = fNShells[Z]*6 - 4;
228  if (energy < (*(fParam[Z]))[idx-1]) { energy = (*(fParam[Z]))[idx-1]; }
229 
230  G4double x1 = 1.0/energy;
231  G4double x2 = x1*x1;
232  G4double x3 = x2*x1;
233 
234  // parameterisation
235  if(energy >= (*(fParam[Z]))[0]) {
236  G4double x4 = x2*x2;
237  cs = x1*((*(fParam[Z]))[idx] + x1*(*(fParam[Z]))[idx+1]
238  + x2*(*(fParam[Z]))[idx+2] + x3*(*(fParam[Z]))[idx+3]
239  + x4*(*(fParam[Z]))[idx+4]);
240  // high energy part
241  } else if(energy >= (*(fParam[Z]))[1]) {
242  cs = x3*(fCrossSection[Z])->Value(energy);
243 
244  // low energy part
245  } else {
246  cs = x3*(fCrossSectionLE[Z])->Value(energy);
247  }
248  if (verboseLevel > 1) {
249  G4cout << "LivermorePolarizedPhotoElectricModel: E(keV)= " << energy/keV
250  << " Z= " << Z << " cross(barn)= " << cs/barn << G4endl;
251  }
252  return cs;
253 }
254 
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
256 
257 void
259  std::vector<G4DynamicParticle*>* fvect,
260  const G4MaterialCutsCouple* couple,
261  const G4DynamicParticle* aDynamicGamma,
262  G4double,
263  G4double)
264 {
265  G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
266  if (verboseLevel > 3) {
267  G4cout << "G4LivermorePolarizedPhotoElectricModel::SampleSecondaries() Egamma(keV)= "
268  << gammaEnergy/keV << G4endl;
269  }
270 
271  // kill incident photon
274 
275  // low-energy photo-effect in water - full absorption
276  const G4Material* material = couple->GetMaterial();
277  if(fWater && (material == fWater ||
278  material->GetBaseMaterial() == fWater)) {
279  if(gammaEnergy <= fWaterEnergyLimit) {
281  return;
282  }
283  }
284 
285  // Returns the normalized direction of the momentum
286  G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
287 
288  // Select randomly one element in the current material
289  //G4cout << "Select random atom Egamma(keV)= " << gammaEnergy/keV << G4endl;
290  const G4Element* elm = SelectRandomAtom(material, theGamma, gammaEnergy);
291  G4int Z = G4lrint(elm->GetZ());
292 
293  // Select the ionised shell in the current atom according to shell
294  // cross sections
295  // G4cout << "Select random shell Z= " << Z << G4endl;
296 
297  if(Z >= maxZ) { Z = maxZ-1; }
298 
299  // element was not initialised gamma should be absorbed
300  if(!fCrossSection[Z]) {
302  return;
303  }
304 
305  // shell index
306  size_t shellIdx = 0;
307  size_t nn = fNShellsUsed[Z];
308 
309  if(nn > 1) {
310  if(gammaEnergy >= (*(fParam[Z]))[0]) {
311  G4double x1 = 1.0/gammaEnergy;
312  G4double x2 = x1*x1;
313  G4double x3 = x2*x1;
314  G4double x4 = x3*x1;
315  G4int idx = nn*6 - 4;
316  // when do sampling common factors are not taken into account
317  // so cross section is not real
318  G4double cs0 = G4UniformRand()*((*(fParam[Z]))[idx]
319  + x1*(*(fParam[Z]))[idx+1]
320  + x2*(*(fParam[Z]))[idx+2]
321  + x3*(*(fParam[Z]))[idx+3]
322  + x4*(*(fParam[Z]))[idx+4]);
323  for(shellIdx=0; shellIdx<nn; ++shellIdx) {
324  idx = shellIdx*6 + 2;
325  if(gammaEnergy > (*(fParam[Z]))[idx-1]) {
326  G4double cs = (*(fParam[Z]))[idx] + x1*(*(fParam[Z]))[idx+1]
327  + x2*(*(fParam[Z]))[idx+2] + x3*(*(fParam[Z]))[idx+3]
328  + x4*(*(fParam[Z]))[idx+4];
329  if(cs >= cs0) { break; }
330  }
331  }
332  if(shellIdx >= nn) { shellIdx = nn-1; }
333 
334  } else {
335 
336  // when do sampling common factors are not taken into account
337  // so cross section is not real
338  G4double cs = G4UniformRand();
339 
340  if(gammaEnergy >= (*(fParam[Z]))[1]) {
341  cs *= (fCrossSection[Z])->Value(gammaEnergy);
342  } else {
343  cs *= (fCrossSectionLE[Z])->Value(gammaEnergy);
344  }
345 
346  for(size_t j=0; j<nn; ++j) {
347  shellIdx = (size_t)fShellCrossSection->GetComponentID(Z, j);
348  if(gammaEnergy > (*(fParam[Z]))[6*shellIdx+1]) {
349  cs -= fShellCrossSection->GetValueForComponent(Z, j, gammaEnergy);
350  }
351  if(cs <= 0.0 || j+1 == nn) { break; }
352  }
353  }
354  }
355 
356  G4double bindingEnergy = (*(fParam[Z]))[shellIdx*6 + 1];
357  //G4cout << "Z= " << Z << " shellIdx= " << shellIdx
358  // << " nShells= " << fNShells[Z]
359  // << " Ebind(keV)= " << bindingEnergy/keV
360  // << " Egamma(keV)= " << gammaEnergy/keV << G4endl;
361 
362  const G4AtomicShell* shell = 0;
363 
364  // no de-excitation from the last shell
365  if(fDeexcitationActive && shellIdx + 1 < nn) {
367  shell = fAtomDeexcitation->GetAtomicShell(Z, as);
368  }
369 
370  // If binding energy of the selected shell is larger than photon energy
371  // do not generate secondaries
372  if(gammaEnergy < bindingEnergy) {
374  return;
375  }
376 
377  // Primary outcoming electron
378  G4double eKineticEnergy = gammaEnergy - bindingEnergy;
380 
381  // Calculate direction of the photoelectron
382  G4ThreeVector electronDirection =
383  GetAngularDistribution()->SampleDirection(aDynamicGamma,
384  eKineticEnergy,
385  shellIdx,
386  couple->GetMaterial());
387 
388  // The electron is created
390  electronDirection,
391  eKineticEnergy);
392  fvect->push_back(electron);
393 
394  // Sample deexcitation
395  if(shell) {
396  G4int index = couple->GetIndex();
398  G4int nbefore = fvect->size();
399 
400  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
401  G4int nafter = fvect->size();
402  if(nafter > nbefore) {
403  G4double esec = 0.0;
404  for (G4int j=nbefore; j<nafter; ++j) {
405 
406  G4double e = ((*fvect)[j])->GetKineticEnergy();
407  if(esec + e > edep) {
408  // correct energy in order to have energy balance
409  e = edep - esec;
410  ((*fvect)[j])->SetKineticEnergy(e);
411  esec += e;
412  // delete the rest of secondaries (should not happens)
413  for (G4int jj=nafter-1; jj>j; --jj) {
414  delete (*fvect)[jj];
415  fvect->pop_back();
416  }
417  break;
418  }
419  esec += e;
420  }
421  edep -= esec;
422  }
423  }
424  }
425  // energy balance - excitation energy left
426  if(edep > 0.0) {
428  }
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
432 
433 void
435 {
436  if (verboseLevel > 1)
437  {
438  G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
439  << G4endl;
440  }
441 
442  if(fCrossSection[Z]) { return; }
443 
444  const char* datadir = path;
445 
446  if(!datadir)
447  {
448  datadir = std::getenv("G4LEDATA");
449  if(!datadir)
450  {
451  G4Exception("G4LivermorePolarizedPhotoElectricModel::ReadData()",
452  "em0006",FatalException,
453  "Environment variable G4LEDATA not defined");
454  return;
455  }
456  }
457 
458  // spline for photoeffect total x-section above K-shell
460  fCrossSection[Z]->SetSpline(true);
461 
462  std::ostringstream ost;
463  ost << datadir << "/livermore/phot/pe-cs-" << Z <<".dat";
464  std::ifstream fin(ost.str().c_str());
465  if( !fin.is_open()) {
467  ed << "G4LivermorePolarizedPhotoElectricModel data file <" << ost.str().c_str()
468  << "> is not opened!" << G4endl;
469  G4Exception("G4LivermorePolarizedPhotoElectricModel::ReadData()",
470  "em0003",FatalException,
471  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
472  return;
473  } else {
474  if(verboseLevel > 3) { G4cout << "File " << ost.str().c_str()
475  << " is opened by G4LivermorePolarizedPhotoElectricModel" << G4endl;}
476  fCrossSection[Z]->Retrieve(fin, true);
478  fin.close();
479  }
480 
481  fParam[Z] = new std::vector<G4double>;
482 
483  // read fit parameters
484  G4int n1 = 0;
485  G4int n2 = 0;
486  G4double x;
487  std::ostringstream ost1;
488  ost1 << datadir << "/livermore/phot/pe-" << Z <<".dat";
489  std::ifstream fin1(ost1.str().c_str());
490  if( !fin1.is_open()) {
492  ed << "G4LivermorePolarizedPhotoElectricModel data file <" << ost1.str().c_str()
493  << "> is not opened!" << G4endl;
494  G4Exception("G4LivermorePolarizedPhotoElectricModel::ReadData()",
495  "em0003",FatalException,
496  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
497  return;
498  } else {
499  if(verboseLevel > 3) {
500  G4cout << "File " << ost1.str().c_str()
501  << " is opened by G4LivermorePolarizedPhotoElectricModel" << G4endl;
502  }
503  fin1 >> n1;
504  if(fin1.fail()) { return; }
505  if(0 > n1 || n1 >= INT_MAX) { n1 = 0; }
506 
507  fin1 >> n2;
508  if(fin1.fail()) { return; }
509  if(0 > n2 || n2 >= INT_MAX) { n2 = 0; }
510 
511  fin1 >> x;
512  if(fin1.fail()) { return; }
513 
514  fNShells[Z] = n1;
515  fParam[Z]->reserve(6*n1+1);
516  fParam[Z]->push_back(x*MeV);
517  for(G4int i=0; i<n1; ++i) {
518  for(G4int j=0; j<6; ++j) {
519  fin1 >> x;
520  if(0 == j) { x *= MeV; }
521  else { x *= barn; }
522  fParam[Z]->push_back(x);
523  }
524  }
525  fin1.close();
526  }
527  // there is a possibility to used only main shells
528  if(nShellLimit < n2) { n2 = nShellLimit; }
530  fNShellsUsed[Z] = n2;
531 
532  if(1 < n2) {
533  std::ostringstream ost2;
534  ost2 << datadir << "/livermore/phot/pe-ss-cs-" << Z <<".dat";
535  std::ifstream fin2(ost2.str().c_str());
536  if( !fin2.is_open()) {
538  ed << "G4LivermorePolarizedPhotoElectricModel data file <" << ost2.str().c_str()
539  << "> is not opened!" << G4endl;
540  G4Exception("G4LivermorePolarizedPhotoElectricModel::ReadData()",
541  "em0003",FatalException,
542  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
543  return;
544  } else {
545  if(verboseLevel > 3) {
546  G4cout << "File " << ost2.str().c_str()
547  << " is opened by G4LivermorePolarizedPhotoElectricModel" << G4endl;
548  }
549 
550  G4int n3, n4;
551  G4double y;
552  for(G4int i=0; i<n2; ++i) {
553  fin2 >> x >> y >> n3 >> n4;
555  for(G4int j=0; j<n3; ++j) {
556  fin2 >> x >> y;
557  v->PutValues(j, x*MeV, y*barn);
558  }
559  fShellCrossSection->AddComponent(Z, n4, v);
560  }
561  fin2.close();
562  }
563  }
564 
565  // no spline for photoeffect total x-section below K-shell
566  if(1 < fNShells[Z]) {
568 
569  std::ostringstream ost3;
570  ost3 << datadir << "/livermore/phot/pe-le-cs-" << Z <<".dat";
571  std::ifstream fin3(ost3.str().c_str());
572  if( !fin3.is_open()) {
574  ed << "G4LivermorePolarizedPhotoElectricModel data file <" << ost3.str().c_str()
575  << "> is not opened!" << G4endl;
576  G4Exception("G4LivermorePolarizedPhotoElectricModel::ReadData()",
577  "em0003",FatalException,
578  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
579  return;
580  } else {
581  if(verboseLevel > 3) {
582  G4cout << "File " << ost3.str().c_str()
583  << " is opened by G4LivermorePolarizedPhotoElectricModel" << G4endl;
584  }
585  fCrossSectionLE[Z]->Retrieve(fin3, true);
587  fin3.close();
588  }
589  }
590 }
591 
592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
593 
594 #include "G4AutoLock.hh"
595 namespace { G4Mutex LivermorePolarizedPhotoElectricModelMutex = G4MUTEX_INITIALIZER; }
596 
598  const G4ParticleDefinition*, G4int Z)
599 {
600  G4AutoLock l(&LivermorePolarizedPhotoElectricModelMutex);
601  // G4cout << "G4LivermorePolarizedPhotoElectricModel::InitialiseForElement Z= "
602  // << Z << G4endl;
603  if(!fCrossSection[Z]) { ReadData(Z); }
604  l.unlock();
605 }
606 
607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....