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