ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermorePolarizedPhotoElectricGDModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermorePolarizedPhotoElectricGDModel.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 // Authors: G.Depaola & F.Longo
28 //
29 
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4LossTableManager.hh"
34 #include "G4VAtomDeexcitation.hh"
35 #include "G4AtomicShell.hh"
36 #include "G4Gamma.hh"
37 #include "G4Electron.hh"
38 #include "G4LPhysicsFreeVector.hh"
40 
41 #include <vector>
42 
45 std::vector<G4double>* G4LivermorePolarizedPhotoElectricGDModel::fParam[] = {0};
51 
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
55 using namespace std;
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 
60  const G4String& nam)
61  :G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
62  nShellLimit(100), fDeexcitationActive(false), isInitialised(false),
63  fAtomDeexcitation(nullptr)
64 {
65  verboseLevel= 0;
66  // Verbosity scale:
67  // 0 = nothing
68  // 1 = warning for energy non-conservation
69  // 2 = details of energy budget
70  // 3 = calculation of cross sections, file openings, sampling of atoms
71  // 4 = entering in methods
72 
75 
76  SetDeexcitationFlag(true);
77  fSandiaCof.resize(4,0.0);
78  fCurrSection = 0.0;
79 
80  if (verboseLevel > 0) {
81  G4cout << "Livermore Polarized PhotoElectric is constructed "
82  << " nShellLimit "
83  << nShellLimit << G4endl;
84  }
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
90 {
91  if(IsMaster()) {
92  delete fShellCrossSection;
93  for(G4int i=0; i<maxZ; ++i) {
94  delete fParam[i];
95  fParam[i] = 0;
96  delete fCrossSection[i];
97  fCrossSection[i] = 0;
98  delete fCrossSectionLE[i];
99  fCrossSectionLE[i] = 0;
100  }
101  }
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 
106 void
108  const G4ParticleDefinition*,
109  const G4DataVector&)
110 {
111  if (verboseLevel > 2) {
112  G4cout << "Calling G4LivermorePolarizedPhotoElectricGDModel::Initialise()" << G4endl;
113  }
114 
115  if(IsMaster()) {
116 
117  if(!fWater) {
118  fWater = G4Material::GetMaterial("G4_WATER", false);
119  if(fWater) { fWaterEnergyLimit = 13.6*eV; }
120  }
121 
123 
124  char* path = std::getenv("G4LEDATA");
125 
126  G4ProductionCutsTable* theCoupleTable =
128  G4int numOfCouples = theCoupleTable->GetTableSize();
129 
130  for(G4int i=0; i<numOfCouples; ++i) {
131  const G4MaterialCutsCouple* couple =
132  theCoupleTable->GetMaterialCutsCouple(i);
133  const G4Material* material = couple->GetMaterial();
134  const G4ElementVector* theElementVector = material->GetElementVector();
135  G4int nelm = material->GetNumberOfElements();
136 
137  for (G4int j=0; j<nelm; ++j) {
138  G4int Z = (G4int)(*theElementVector)[j]->GetZ();
139  if(Z < 1) { Z = 1; }
140  else if(Z > maxZ) { Z = maxZ; }
141  if(!fCrossSection[Z]) { ReadData(Z, path); }
142  }
143  }
144  }
145  //
146  if (verboseLevel > 2) {
147  G4cout << "Loaded cross section files for LivermorePhotoElectric model"
148  << G4endl;
149  }
150  if(!isInitialised) {
151  isInitialised = true;
153 
155  }
156  fDeexcitationActive = false;
157  if(fAtomDeexcitation) {
159  }
160 
161  if (verboseLevel > 0) {
162  G4cout << "LivermorePolarizedPhotoElectric model is initialized " << G4endl
163  << G4endl;
164  }
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168 
170  const G4ParticleDefinition*,
171  G4double GammaEnergy,
172  G4double ZZ, G4double,
174 {
175  if (verboseLevel > 3) {
176  G4cout << "G4LivermorePolarizedPhotoElectricGDModel::ComputeCrossSectionPerAtom():"
177  << " Z= " << ZZ << " R(keV)= " << GammaEnergy/keV << G4endl;
178  }
179  G4double cs = 0.0;
180  G4int Z = G4lrint(ZZ);
181  if(Z < 1 || Z >= maxZ) { return cs; }
182 
183  // if element was not initialised
184  // do initialisation safely for MT mode
185  if(!fCrossSection[Z]) {
186  InitialiseForElement(0, Z);
187  if(!fCrossSection[Z]) { return cs; }
188  }
189 
190  G4int idx = fNShells[Z]*6 - 4;
191  if (GammaEnergy < (*(fParam[Z]))[idx-1]) { GammaEnergy = (*(fParam[Z]))[idx-1]; }
192 
193  G4double x1 = 1.0/GammaEnergy;
194  G4double x2 = x1*x1;
195  G4double x3 = x2*x1;
196 
197  // parameterisation
198  if(GammaEnergy >= (*(fParam[Z]))[0]) {
199  G4double x4 = x2*x2;
200  cs = x1*((*(fParam[Z]))[idx] + x1*(*(fParam[Z]))[idx+1]
201  + x2*(*(fParam[Z]))[idx+2] + x3*(*(fParam[Z]))[idx+3]
202  + x4*(*(fParam[Z]))[idx+4]);
203  // high energy part
204  } else if (GammaEnergy >= (*(fParam[Z]))[1]) {
205  cs = x3*(fCrossSection[Z])->Value(GammaEnergy);
206 
207  // low energy part
208  } else {
209  cs = x3*(fCrossSectionLE[Z])->Value(GammaEnergy);
210  }
211  if (verboseLevel > 1) {
212  G4cout << "LivermorePolarizedPhotoElectricGDModel: E(keV)= " << GammaEnergy/keV
213  << " Z= " << Z << " cross(barn)= " << cs/barn << G4endl;
214  }
215  return cs;
216 }
217 
218 
219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220 
222  std::vector<G4DynamicParticle*>* fvect,
223  const G4MaterialCutsCouple* couple,
224  const G4DynamicParticle* aDynamicGamma,
225  G4double,
226  G4double)
227 {
228  if (verboseLevel > 3) {
229  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricGDModel"
230  << G4endl;
231  }
232 
233  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
234  if (verboseLevel > 3) {
235  G4cout << "G4LivermorePolarizedPhotoElectricGDModel::SampleSecondaries() Egamma(keV)= "
236  << photonEnergy/keV << G4endl;
237  }
238 
239  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
240  G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
241 
242  // kill incident photon
245 
246  // low-energy photo-effect in water - full absorption
247 
248  const G4Material* material = couple->GetMaterial();
249  if(fWater && (material == fWater ||
250  material->GetBaseMaterial() == fWater)) {
251  if(photonEnergy <= fWaterEnergyLimit) {
253  return;
254  }
255  }
256 
257  // Protection: a polarisation parallel to the
258  // direction causes problems;
259  // in that case find a random polarization
260 
261  // Make sure that the polarization vector is perpendicular to the
262  // gamma direction. If not
263 
264  if(!(gammaPolarization0.isOrthogonal(photonDirection, 1e-6))||(gammaPolarization0.mag()==0))
265  { // only for testing now
266  gammaPolarization0 = GetRandomPolarization(photonDirection);
267  }
268  else
269  {
270  if ( gammaPolarization0.howOrthogonal(photonDirection) != 0)
271  {
272  gammaPolarization0 = GetPerpendicularPolarization(photonDirection, gammaPolarization0);
273  }
274  }
275 
276  // End of Protection
277 
278  // G4double E0_m = photonEnergy / electron_mass_c2 ;
279 
280  // Shell
281 
282  // Select randomly one element in the current material
283  //G4cout << "Select random atom Egamma(keV)= " << photonEnergy/keV << G4endl;
284  const G4Element* elm = SelectRandomAtom(material, theGamma, photonEnergy);
285  G4int Z = G4lrint(elm->GetZ());
286 
287 
288  // Select the ionised shell in the current atom according to shell
289  // cross sections
290  // G4cout << "Select random shell Z= " << Z << G4endl;
291 
292  if(Z >= maxZ) { Z = maxZ-1; }
293 
294  // element was not initialised gamma should be absorbed
295  if(!fCrossSection[Z]) {
297  return;
298  }
299 
300  // shell index
301  size_t shellIdx = 0;
302  size_t nn = fNShellsUsed[Z];
303 
304  if(nn > 1) {
305  if(photonEnergy >= (*(fParam[Z]))[0]) {
306  G4double x1 = 1.0/photonEnergy;
307  G4double x2 = x1*x1;
308  G4double x3 = x2*x1;
309  G4double x4 = x3*x1;
310  G4int idx = nn*6 - 4;
311  // when do sampling common factors are not taken into account
312  // so cross section is not real
313  G4double cs0 = G4UniformRand()*((*(fParam[Z]))[idx]
314  + x1*(*(fParam[Z]))[idx+1]
315  + x2*(*(fParam[Z]))[idx+2]
316  + x3*(*(fParam[Z]))[idx+3]
317  + x4*(*(fParam[Z]))[idx+4]);
318  for(shellIdx=0; shellIdx<nn; ++shellIdx) {
319  idx = shellIdx*6 + 2;
320  if(photonEnergy > (*(fParam[Z]))[idx-1]) {
321  G4double cs = (*(fParam[Z]))[idx] + x1*(*(fParam[Z]))[idx+1]
322  + x2*(*(fParam[Z]))[idx+2] + x3*(*(fParam[Z]))[idx+3]
323  + x4*(*(fParam[Z]))[idx+4];
324  if(cs >= cs0) { break; }
325  }
326  }
327  if(shellIdx >= nn) { shellIdx = nn-1; }
328 
329  } else {
330 
331  // when do sampling common factors are not taken into account
332  // so cross section is not real
333  G4double cs = G4UniformRand();
334 
335  if(photonEnergy >= (*(fParam[Z]))[1]) {
336  cs *= (fCrossSection[Z])->Value(photonEnergy);
337  } else {
338  cs *= (fCrossSectionLE[Z])->Value(photonEnergy);
339  }
340 
341  for(size_t j=0; j<nn; ++j) {
342  shellIdx = (size_t)fShellCrossSection->GetComponentID(Z, j);
343  if(photonEnergy > (*(fParam[Z]))[6*shellIdx+1]) {
344  cs -= fShellCrossSection->GetValueForComponent(Z, j, photonEnergy);
345  }
346  if(cs <= 0.0 || j+1 == nn) { break; }
347  }
348  }
349  }
350 
351  G4double bindingEnergy = (*(fParam[Z]))[shellIdx*6 + 1];
352  //G4cout << "Z= " << Z << " shellIdx= " << shellIdx
353  // << " nShells= " << fNShells[Z]
354  // << " Ebind(keV)= " << bindingEnergy/keV
355  // << " Egamma(keV)= " << photonEnergy/keV << G4endl;
356 
357  const G4AtomicShell* shell = 0;
358 
359  // no de-excitation from the last shell
360  if(fDeexcitationActive && shellIdx + 1 < nn) {
362  shell = fAtomDeexcitation->GetAtomicShell(Z, as);
363  }
364 
365  // If binding energy of the selected shell is larger than photon energy
366  // do not generate secondaries
367  if(photonEnergy < bindingEnergy) {
369  return;
370  }
371 
372 
373  // Electron
374 
375  G4double eKineticEnergy = photonEnergy - bindingEnergy;
377 
378  G4double costheta = SetCosTheta(eKineticEnergy);
379  G4double sintheta = sqrt(1. - costheta*costheta);
380  G4double phi = SetPhi(photonEnergy,eKineticEnergy,costheta);
381  G4double dirX = sintheta*cos(phi);
382  G4double dirY = sintheta*sin(phi);
383  G4double dirZ = costheta;
384  G4ThreeVector electronDirection(dirX, dirY, dirZ);
385  SystemOfRefChange(photonDirection, electronDirection, gammaPolarization0);
387  electronDirection,
388  eKineticEnergy);
389  fvect->push_back(electron);
390 
391  // Deexcitation
392  // Sample deexcitation
393  if(shell) {
394  G4int index = couple->GetIndex();
396  G4int nbefore = fvect->size();
397 
398  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
399  G4int nafter = fvect->size();
400  if(nafter > nbefore) {
401  G4double esec = 0.0;
402  for (G4int j=nbefore; j<nafter; ++j) {
403 
404  G4double e = ((*fvect)[j])->GetKineticEnergy();
405  if(esec + e > edep) {
406  // correct energy in order to have energy balance
407  e = edep - esec;
408  ((*fvect)[j])->SetKineticEnergy(e);
409  esec += e;
410  // delete the rest of secondaries (should not happens)
411  for (G4int jj=nafter-1; jj>j; --jj) {
412  delete (*fvect)[jj];
413  fvect->pop_back();
414  }
415  break;
416  }
417  esec += e;
418  }
419  edep -= esec;
420  }
421  }
422  }
423  // energy balance - excitation energy left
424  if(edep > 0.0) {
426  }
427 }
428 
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
430 
432 {
433  if (verboseLevel > 1)
434  {
435  G4cout << "Calling ReadData() of G4LivermorePolarizedPhotoElectricGDModel"
436  << G4endl;
437  }
438 
439  if(fCrossSection[Z]) { return; }
440 
441  const char* datadir = path;
442 
443  if(!datadir)
444  {
445  datadir = std::getenv("G4LEDATA");
446  if(!datadir)
447  {
448  G4Exception("G4LivermorePhotoElectricModel::ReadData()",
449  "em0006",FatalException,
450  "Environment variable G4LEDATA not defined");
451  return;
452  }
453  }
454 
455  // spline for photoeffect total x-section above K-shell
457  fCrossSection[Z]->SetSpline(true);
458 
459  std::ostringstream ost;
460  ost << datadir << "/livermore/phot/pe-cs-" << Z <<".dat";
461  std::ifstream fin(ost.str().c_str());
462  if( !fin.is_open()) {
464  ed << "G4LivermorePolarizedPhotoElectricGDModel data file <" << ost.str().c_str()
465  << "> is not opened!" << G4endl;
466  G4Exception("G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
467  "em0003",FatalException,
468  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
469  return;
470  } else {
471  if(verboseLevel > 3) { G4cout << "File " << ost.str().c_str()
472  << " is opened by G4LivermorePolarizedPhotoElectricGDModel" << G4endl;}
473  fCrossSection[Z]->Retrieve(fin, true);
475  fin.close();
476  }
477 
478  fParam[Z] = new std::vector<G4double>;
479 
480  // read fit parameters
481  G4int n1 = 0;
482  G4int n2 = 0;
483  G4double x;
484  std::ostringstream ost1;
485  ost1 << datadir << "/livermore/phot/pe-" << Z <<".dat";
486  std::ifstream fin1(ost1.str().c_str());
487  if( !fin1.is_open()) {
489  ed << "G4LivermorePolarizedPhotoElectricGDModel data file <" << ost1.str().c_str()
490  << "> is not opened!" << G4endl;
491  G4Exception("G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
492  "em0003",FatalException,
493  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
494  return;
495  } else {
496  if(verboseLevel > 3) {
497  G4cout << "File " << ost1.str().c_str()
498  << " is opened by G4LivermorePolarizedPhotoElectricGDModel" << G4endl;
499  }
500  fin1 >> n1;
501  if(fin1.fail()) { return; }
502  if(0 > n1 || n1 >= INT_MAX) { n1 = 0; }
503 
504  fin1 >> n2;
505  if(fin1.fail()) { return; }
506  if(0 > n2 || n2 >= INT_MAX) { n2 = 0; }
507 
508  fin1 >> x;
509  if(fin1.fail()) { return; }
510 
511  fNShells[Z] = n1;
512  fParam[Z]->reserve(6*n1+1);
513  fParam[Z]->push_back(x*MeV);
514  for(G4int i=0; i<n1; ++i) {
515  for(G4int j=0; j<6; ++j) {
516  fin1 >> x;
517  if(0 == j) { x *= MeV; }
518  else { x *= barn; }
519  fParam[Z]->push_back(x);
520  }
521  }
522  fin1.close();
523  }
524  // there is a possibility to used only main shells
525  if(nShellLimit < n2) { n2 = nShellLimit; }
527  fNShellsUsed[Z] = n2;
528 
529  if(1 < n2) {
530  std::ostringstream ost2;
531  ost2 << datadir << "/livermore/phot/pe-ss-cs-" << Z <<".dat";
532  std::ifstream fin2(ost2.str().c_str());
533  if( !fin2.is_open()) {
535  ed << "G4LivermorePolarizedPhotoElectricGDModel data file <" << ost2.str().c_str()
536  << "> is not opened!" << G4endl;
537  G4Exception("G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
538  "em0003",FatalException,
539  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
540  return;
541  } else {
542  if(verboseLevel > 3) {
543  G4cout << "File " << ost2.str().c_str()
544  << " is opened by G4LivermorePolarizedPhotoElectricGDModel" << G4endl;
545  }
546 
547  G4int n3, n4;
548  G4double y;
549  for(G4int i=0; i<n2; ++i) {
550  fin2 >> x >> y >> n3 >> n4;
552  for(G4int j=0; j<n3; ++j) {
553  fin2 >> x >> y;
554  v->PutValues(j, x*MeV, y*barn);
555  }
556  fShellCrossSection->AddComponent(Z, n4, v);
557  }
558  fin2.close();
559  }
560  }
561 
562  // no spline for photoeffect total x-section below K-shell
563  if(1 < fNShells[Z]) {
565 
566  std::ostringstream ost3;
567  ost3 << datadir << "/livermore/phot/pe-le-cs-" << Z <<".dat";
568  std::ifstream fin3(ost3.str().c_str());
569  if( !fin3.is_open()) {
571  ed << "G4LivermorePolarizedPhotoElectricGDModel data file <" << ost3.str().c_str()
572  << "> is not opened!" << G4endl;
573  G4Exception("G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
574  "em0003",FatalException,
575  ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
576  return;
577  } else {
578  if(verboseLevel > 3) {
579  G4cout << "File " << ost3.str().c_str()
580  << " is opened by G4LivermorePolarizedPhotoElectricGDModel" << G4endl;
581  }
582  fCrossSectionLE[Z]->Retrieve(fin3, true);
584  fin3.close();
585  }
586  }
587 }
588 
589 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
590 
591 
593 {
594  G4double rand1,rand2,onemcost,greject;
595  G4double masarep = 510.99906*keV;
596 
597  G4double gamma = 1. + energyE/masarep;
598  G4double gamma2 = gamma*gamma;
599 
600  G4double beta = sqrt((gamma2 - 1.)/gamma2);
601 
602  G4double alfa = 1./beta - 1.;
603 
604  G4double g1 = 0.5*beta*gamma*(gamma-1.)*(gamma-2.);
605 
606  G4double alfap2 = alfa+2.;
607 
608  G4double grejectmax = 2.*(g1+1./alfa);
609 
610  do
611  {
612  rand1 = G4UniformRand();
613  onemcost = 2.*alfa*(2.*rand1 + alfap2 * sqrt(rand1))/
614  (alfap2*alfap2 - 4.*rand1);
615  greject = (2. - onemcost)*(g1+1./(alfa+onemcost));
616  rand2 = G4UniformRand();
617  }
618  while (rand2*grejectmax > greject);
619  G4double cosTheta = 1. - onemcost;
620  return cosTheta;
621 }
622 
623 
625  G4double E_energy,
626  G4double costheta)
627 {
628  G4double epsilon = E_energy/electron_mass_c2;
629  G4double k = Ph_energy/electron_mass_c2;
630  G4double gamma = 1. + epsilon;
631  G4double gamma2 = gamma*gamma;
632  G4double beta = sqrt((gamma2 - 1.)/gamma2);
633 
634  G4double d = (2./(k*gamma*(1-beta*costheta))-1)*(1/k);
635 
636  G4double norm_factor = 1 +2*d;
637 
638  G4double rnd1;
639  G4double rnd2;
640  G4double phi, phiprob;
641 
642  do
643  {
644  rnd1 =G4UniformRand();
645  rnd2 =G4UniformRand();
646  phi = rnd1*twopi;
647  phiprob = 1 +2*d*cos(phi)*cos(phi);
648  }
649  while (rnd2*norm_factor > phiprob);
650  return phi;
651 }
652 
653 
654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
655 
657 {
658  G4double dx = a.x();
659  G4double dy = a.y();
660  G4double dz = a.z();
661  G4double x = dx < 0.0 ? -dx : dx;
662  G4double y = dy < 0.0 ? -dy : dy;
663  G4double z = dz < 0.0 ? -dz : dz;
664  if (x < y) {
665  return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
666  }else{
667  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
668  }
669 }
670 
671 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
672 
674 {
675  G4ThreeVector d0 = direction0.unit();
676  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
677  G4ThreeVector a0 = a1.unit(); // unit vector
678 
680 
681  G4double angle = twopi*rand1; // random polar angle
682  G4ThreeVector b0 = d0.cross(a0); // cross product
683 
685 
686  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
687  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
688  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
689 
690  G4ThreeVector c0 = c.unit();
691 
692  return c0;
693 
694 }
695 
696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
697 
699 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
700 {
701 
702  //
703  // The polarization of a photon is always perpendicular to its momentum direction.
704  // Therefore this function removes those vector component of gammaPolarization, which
705  // points in direction of gammaDirection
706  //
707 
708  // Mathematically we search the projection of the vector a on the plane E, where n is the
709  // plains normal vector.
710  // The basic equation can be found in each geometry book (e.g. Bronstein):
711  // p = a - (a o n)/(n o n)*n
712 
713  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
714 
715 }
716 
717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
718 
719 
721 (G4ThreeVector& direction0,G4ThreeVector& direction1,
722  G4ThreeVector& polarization0)
723 {
724  // direction0 is the original photon direction ---> z
725  // polarization0 is the original photon polarization ---> x
726  // need to specify y axis in the real reference frame ---> y
727  G4ThreeVector Axis_Z0 = direction0.unit();
728  G4ThreeVector Axis_X0 = polarization0.unit();
729  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
730 
731  G4double direction_x = direction1.getX();
732  G4double direction_y = direction1.getY();
733  G4double direction_z = direction1.getZ();
734 
735  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
736 
737 }
738 
739 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
740 
741 #include "G4AutoLock.hh"
742 namespace { G4Mutex LivermorePolarizedPhotoElectricGDModelMutex = G4MUTEX_INITIALIZER; }
743 
745  const G4ParticleDefinition*, G4int Z)
746 {
747  G4AutoLock l(&LivermorePolarizedPhotoElectricGDModelMutex);
748  // G4cout << "G4LivermorePhotoElectricModel::InitialiseForElement Z= "
749  // << Z << G4endl;
750  if(!fCrossSection[Z]) { ReadData(Z); }
751  l.unlock();
752 }
753 
754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
755 
756