Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEPPolarizedComptonModel.cc
Go to the documentation of this file.
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 // | G4LowEPPolarizedComptonModel-- Geant4 Monash University |
29 // | polarised low energy Compton scattering model. |
30 // | J. M. C. Brown, Monash University, Australia |
31 // | |
32 // | |
33 // *********************************************************************
34 // | |
35 // | The following is a Geant4 class to simulate the process of |
36 // | bound electron Compton scattering. General code structure is |
37 // | based on G4LowEnergyCompton.cc and |
38 // | G4LivermorePolarizedComptonModel.cc. |
39 // | Algorithms for photon energy, and ejected Compton electron |
40 // | direction taken from: |
41 // | |
42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin, |
43 // | "A low energy bound atomic electron Compton scattering model |
44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014. |
45 // | |
46 // | The author acknowledges the work of the Geant4 collaboration |
47 // | in developing the following algorithms that have been employed |
48 // | or adapeted for the present software: |
49 // | |
50 // | # sampling of photon scattering angle, |
51 // | # target element selection in composite materials, |
52 // | # target shell selection in element, |
53 // | # and sampling of bound electron momentum from Compton profiles. |
54 // | |
55 // *********************************************************************
56 // | |
57 // | History: |
58 // | -------- |
59 // | |
60 // | Jan. 2015 JMCB - 1st Version based on G4LowEPPComptonModel |
61 // | |
62 // *********************************************************************
63 
65 #include "G4PhysicalConstants.hh"
66 #include "G4SystemOfUnits.hh"
67 #include "G4Exp.hh"
68 
69 //****************************************************************************
70 
71 using namespace std;
72 
73 G4int G4LowEPPolarizedComptonModel::maxZ = 99;
74 G4LPhysicsFreeVector* G4LowEPPolarizedComptonModel::data[] = {0};
75 G4ShellData* G4LowEPPolarizedComptonModel::shellData = 0;
76 G4DopplerProfile* G4LowEPPolarizedComptonModel::profileData = 0;
77 
78 static const G4double ln10 = G4Log(10.);
79 
81  const G4String& nam)
82  : G4VEmModel(nam),isInitialised(false)
83 {
84  verboseLevel=1 ;
85  // Verbosity scale:
86  // 0 = nothing
87  // 1 = warning for energy non-conservation
88  // 2 = details of energy budget
89  // 3 = calculation of cross sections, file openings, sampling of atoms
90  // 4 = entering in methods
91 
92  if( verboseLevel>1 ) {
93  G4cout << "Low energy photon Compton model is constructed " << G4endl;
94  }
95 
96  //Mark this model as "applicable" for atomic deexcitation
97  SetDeexcitationFlag(true);
98 
99  fParticleChange = 0;
100  fAtomDeexcitation = 0;
101 }
102 
103 //****************************************************************************
104 
106 {
107  if(IsMaster()) {
108  delete shellData;
109  shellData = 0;
110  delete profileData;
111  profileData = 0;
112  }
113 }
114 
115 //****************************************************************************
116 
118  const G4DataVector& cuts)
119 {
120  if (verboseLevel > 1) {
121  G4cout << "Calling G4LowEPPolarizedComptonModel::Initialise()" << G4endl;
122  }
123 
124  // Initialise element selector
125 
126  if(IsMaster()) {
127 
128  // Access to elements
129 
130  char* path = getenv("G4LEDATA");
131 
132  G4ProductionCutsTable* theCoupleTable =
134  G4int numOfCouples = theCoupleTable->GetTableSize();
135 
136  for(G4int i=0; i<numOfCouples; ++i) {
137  const G4Material* material =
138  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
139  const G4ElementVector* theElementVector = material->GetElementVector();
140  G4int nelm = material->GetNumberOfElements();
141 
142  for (G4int j=0; j<nelm; ++j) {
143  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
144  if(Z < 1) { Z = 1; }
145  else if(Z > maxZ){ Z = maxZ; }
146 
147  if( (!data[Z]) ) { ReadData(Z, path); }
148  }
149  }
150 
151  // For Doppler broadening
152  if(!shellData) {
153  shellData = new G4ShellData();
154  shellData->SetOccupancyData();
155  G4String file = "/doppler/shell-doppler";
156  shellData->LoadData(file);
157  }
158  if(!profileData) { profileData = new G4DopplerProfile(); }
159 
160  InitialiseElementSelectors(particle, cuts);
161  }
162 
163  if (verboseLevel > 2) {
164  G4cout << "Loaded cross section files" << G4endl;
165  }
166 
167  if( verboseLevel>1 ) {
168  G4cout << "G4LowEPPolarizedComptonModel is initialized " << G4endl
169  << "Energy range: "
170  << LowEnergyLimit() / eV << " eV - "
171  << HighEnergyLimit() / GeV << " GeV"
172  << G4endl;
173  }
174 
175  if(isInitialised) { return; }
176 
177  fParticleChange = GetParticleChangeForGamma();
178  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
179  isInitialised = true;
180 }
181 
182 //****************************************************************************
183 
185  G4VEmModel* masterModel)
186 {
188 }
189 
190 //****************************************************************************
191 
192 void G4LowEPPolarizedComptonModel::ReadData(size_t Z, const char* path)
193 {
194  if (verboseLevel > 1)
195  {
196  G4cout << "G4LowEPPolarizedComptonModel::ReadData()"
197  << G4endl;
198  }
199  if(data[Z]) { return; }
200  const char* datadir = path;
201  if(!datadir)
202  {
203  datadir = getenv("G4LEDATA");
204  if(!datadir)
205  {
206  G4Exception("G4LowEPPolarizedComptonModel::ReadData()",
207  "em0006",FatalException,
208  "Environment variable G4LEDATA not defined");
209  return;
210  }
211  }
212 
213  data[Z] = new G4LPhysicsFreeVector();
214 
215  // Activation of spline interpolation
216  data[Z]->SetSpline(false);
217 
218  std::ostringstream ost;
219  ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat";
220  std::ifstream fin(ost.str().c_str());
221 
222  if( !fin.is_open())
223  {
225  ed << "G4LowEPPolarizedComptonModel data file <" << ost.str().c_str()
226  << "> is not opened!" << G4endl;
227  G4Exception("G4LowEPPolarizedComptonModel::ReadData()",
228  "em0003",FatalException,
229  ed,"G4LEDATA version should be G4EMLOW6.34 or later");
230  return;
231  } else {
232  if(verboseLevel > 3) {
233  G4cout << "File " << ost.str()
234  << " is opened by G4LowEPPolarizedComptonModel" << G4endl;
235  }
236  data[Z]->Retrieve(fin, true);
237  data[Z]->ScaleVector(MeV, MeV*barn);
238  }
239  fin.close();
240 }
241 
242 //****************************************************************************
243 
244 
245 G4double
247  G4double GammaEnergy,
248  G4double Z, G4double,
250 {
251  if (verboseLevel > 3) {
252  G4cout << "G4LowEPPolarizedComptonModel::ComputeCrossSectionPerAtom()"
253  << G4endl;
254  }
255  G4double cs = 0.0;
256 
257  if (GammaEnergy < LowEnergyLimit()) { return 0.0; }
258 
259  G4int intZ = G4lrint(Z);
260  if(intZ < 1 || intZ > maxZ) { return cs; }
261 
262  G4LPhysicsFreeVector* pv = data[intZ];
263 
264  // if element was not initialised
265  // do initialisation safely for MT mode
266  if(!pv)
267  {
268  InitialiseForElement(0, intZ);
269  pv = data[intZ];
270  if(!pv) { return cs; }
271  }
272 
273  G4int n = pv->GetVectorLength() - 1;
274  G4double e1 = pv->Energy(0);
275  G4double e2 = pv->Energy(n);
276 
277  if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
278  else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
279  else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
280 
281  return cs;
282 }
283 
284 //****************************************************************************
285 
286 void G4LowEPPolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
287  const G4MaterialCutsCouple* couple,
288  const G4DynamicParticle* aDynamicGamma,
290 {
291 
292  //Determine number of digits (in decimal base) that G4double can accurately represent
293  G4double g4d_order = G4double(numeric_limits<G4double>::digits10);
294  G4double g4d_limit = std::pow(10.,-g4d_order);
295 
296  // The scattered gamma energy is sampled according to Klein - Nishina formula.
297  // then accepted or rejected depending on the Scattering Function multiplied
298  // by factor from Klein - Nishina formula.
299  // Expression of the angular distribution as Klein Nishina
300  // angular and energy distribution and Scattering fuctions is taken from
301  // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
302  // Phys. Res. B 101 (1995). Method of sampling with form factors is different
303  // data are interpolated while in the article they are fitted.
304  // Reference to the article is from J. Stepanek New Photon, Positron
305  // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
306  // TeV (draft).
307  // The random number techniques of Butcher & Messel are used
308  // (Nucl Phys 20(1960),15).
309 
310 
311  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
312 
313  if (verboseLevel > 3) {
314  G4cout << "G4LowEPPolarizedComptonModel::SampleSecondaries() E(MeV)= "
315  << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
316  << G4endl;
317  }
318  // do nothing below the threshold
319  // should never get here because the XS is zero below the limit
320  if (photonEnergy0 < LowEnergyLimit())
321  return ;
322 
323  G4double e0m = photonEnergy0 / electron_mass_c2 ;
324  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
325 
326 
327  // Polarisation: check orientation of photon propagation direction and polarisation
328  // Fix if needed
329 
330  G4ThreeVector photonPolarization0 = aDynamicGamma->GetPolarization();
331 
332  // Check if polarisation vector is perpendicular and fix if not
333 
334  if (!(photonPolarization0.isOrthogonal(photonDirection0, 1e-6))||(photonPolarization0.mag()==0))
335  {
336  photonPolarization0 = GetRandomPolarization(photonDirection0);
337  }
338 
339  else
340  {
341  if ((photonPolarization0.howOrthogonal(photonDirection0) !=0) && (photonPolarization0.howOrthogonal(photonDirection0) > g4d_limit))
342  {
343  photonPolarization0 = GetPerpendicularPolarization(photonDirection0,photonPolarization0);
344  }
345  }
346 
347  // Select randomly one element in the current material
348 
349  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
350  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
351  G4int Z = (G4int)elm->GetZ();
352 
353  G4double LowEPPCepsilon0 = 1. / (1. + 2. * e0m);
354  G4double LowEPPCepsilon0Sq = LowEPPCepsilon0 * LowEPPCepsilon0;
355  G4double alpha1 = -std::log(LowEPPCepsilon0);
356  G4double alpha2 = 0.5 * (1. - LowEPPCepsilon0Sq);
357 
358  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
359 
360  // Sample the energy of the scattered photon
361  G4double LowEPPCepsilon;
362  G4double LowEPPCepsilonSq;
363  G4double oneCosT;
364  G4double sinT2;
365  G4double gReject;
366 
367  if (verboseLevel > 3) {
368  G4cout << "Started loop to sample gamma energy" << G4endl;
369  }
370 
371  do
372  {
373  if ( alpha1/(alpha1+alpha2) > G4UniformRand())
374  {
375  LowEPPCepsilon = G4Exp(-alpha1 * G4UniformRand());
376  LowEPPCepsilonSq = LowEPPCepsilon * LowEPPCepsilon;
377  }
378  else
379  {
380  LowEPPCepsilonSq = LowEPPCepsilon0Sq + (1. - LowEPPCepsilon0Sq) * G4UniformRand();
381  LowEPPCepsilon = std::sqrt(LowEPPCepsilonSq);
382  }
383 
384  oneCosT = (1. - LowEPPCepsilon) / ( LowEPPCepsilon * e0m);
385  sinT2 = oneCosT * (2. - oneCosT);
386  G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
387  G4double scatteringFunction = ComputeScatteringFunction(x, Z);
388  gReject = (1. - LowEPPCepsilon * sinT2 / (1. + LowEPPCepsilonSq)) * scatteringFunction;
389 
390  } while(gReject < G4UniformRand()*Z);
391 
392  G4double cosTheta = 1. - oneCosT;
393  G4double sinTheta = std::sqrt(sinT2);
394  G4double phi = SetPhi(LowEPPCepsilon,sinT2);
395  G4double dirx = sinTheta * std::cos(phi);
396  G4double diry = sinTheta * std::sin(phi);
397  G4double dirz = cosTheta ;
398 
399  // Set outgoing photon polarization
400 
401  G4ThreeVector photonPolarization1 = SetNewPolarization(LowEPPCepsilon,
402  sinT2,
403  phi,
404  cosTheta);
405 
406  // Scatter photon energy and Compton electron direction - Method based on:
407  // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin'
408  // "A low energy bound atomic electron Compton scattering model for Geant4"
409  // NIMB, Vol. 338, 77-88, 2014.
410 
411  // Set constants and initialize scattering parameters
412 
413  const G4double vel_c = c_light / (m/s);
414  const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s);
415  const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ;
416 
417  const G4int maxDopplerIterations = 1000;
418  G4double bindingE = 0.;
419  G4double pEIncident = photonEnergy0 ;
420  G4double pERecoil = -1.;
421  G4double eERecoil = -1.;
422  G4double e_alpha =0.;
423  G4double e_beta = 0.;
424 
425  G4double CE_emission_flag = 0.;
426  G4double ePAU = -1;
427  G4int shellIdx = 0;
428  G4double u_temp = 0;
429  G4double cosPhiE =0;
430  G4double sinThetaE =0;
431  G4double cosThetaE =0;
432  G4int iteration = 0;
433 
434  if (verboseLevel > 3) {
435  G4cout << "Started loop to sample photon energy and electron direction" << G4endl;
436  }
437 
438  do{
439 
440 
441  // ******************************************
442  // | Determine scatter photon energy |
443  // ******************************************
444 
445  do
446  {
447  iteration++;
448 
449 
450  // ********************************************
451  // | Sample bound electron information |
452  // ********************************************
453 
454  // Select shell based on shell occupancy
455 
456  shellIdx = shellData->SelectRandomShell(Z);
457  bindingE = shellData->BindingEnergy(Z,shellIdx)/MeV;
458 
459 
460  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
461  ePAU = profileData->RandomSelectMomentum(Z,shellIdx);
462 
463  // Convert to SI units
464  G4double ePSI = ePAU * momentum_au_to_nat;
465 
466  //Calculate bound electron velocity and normalise to natural units
467  u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
468 
469  // Sample incident electron direction, amorphous material, to scattering photon scattering plane
470 
471  e_alpha = pi*G4UniformRand();
472  e_beta = twopi*G4UniformRand();
473 
474  // Total energy of system
475 
476  G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
477  G4double systemE = eEIncident + pEIncident;
478 
479 
480  G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
481  G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
482  G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
483  G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
484  G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
485  pERecoil = (numerator/denominator);
486  eERecoil = systemE - pERecoil;
487  CE_emission_flag = pEIncident - pERecoil;
488  } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
489 
490 // End of recalculation of photon energy with Doppler broadening
491 
492 
493 
494  // *******************************************************
495  // | Determine ejected Compton electron direction |
496  // *******************************************************
497 
498  // Calculate velocity of ejected Compton electron
499 
500  G4double a_temp = eERecoil / electron_mass_c2;
501  G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
502 
503  // Coefficients and terms from simulatenous equations
504 
505  G4double sinAlpha = std::sin(e_alpha);
506  G4double cosAlpha = std::cos(e_alpha);
507  G4double sinBeta = std::sin(e_beta);
508  G4double cosBeta = std::cos(e_beta);
509 
510  G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
511  G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
512 
513  G4double var_A = pERecoil*u_p_temp*sinTheta;
514  G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
515  G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
516 
517  G4double var_D1 = gamma*electron_mass_c2*pERecoil;
518  G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
519  G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
520  G4double var_D = var_D1*var_D2 + var_D3;
521 
522  G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
523  G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
524  G4double var_E = var_E1 - var_E2;
525 
526  G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
527  G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
528  G4double var_F = var_F1 - var_F2;
529 
530  G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
531 
532  // Two equations form a quadratic form of Wx^2 + Yx + Z = 0
533  // Coefficents and solution to quadratic
534 
535  G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
536  G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
537  G4double var_W = var_W1 + var_W2;
538 
539  G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
540 
541  G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
542  G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
543  G4double var_Z = var_Z1 + var_Z2;
544  G4double diff1 = var_Y*var_Y;
545  G4double diff2 = 4*var_W*var_Z;
546  G4double diff = diff1 - diff2;
547 
548 
549  // Check if diff is less than zero, if so ensure it is due to FPE
550 
551  //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less
552  //than 10^(-g4d_order), then set diff to zero
553 
554  if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) )
555  {
556  diff = 0.0;
557  }
558 
559 
560  // Plus and minus of quadratic
561  G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
562  G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
563 
564  // Floating point precision protection
565  // Check if X_p and X_m are greater than or less than 1 or -1, if so clean up FPE
566  // Issue due to propagation of FPE and only impacts 8th sig fig onwards
567 
568  if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;}
569  if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;}
570 
571  // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron
572  G4double ThetaE = 0.;
573 
574 
575  G4double sol_select = G4UniformRand();
576 
577  if (sol_select < 0.5)
578  {
579  ThetaE = std::acos(X_p);
580  }
581  if (sol_select > 0.5)
582  {
583  ThetaE = std::acos(X_m);
584  }
585 
586  cosThetaE = std::cos(ThetaE);
587  sinThetaE = std::sin(ThetaE);
588  G4double Theta = std::acos(cosTheta);
589 
590  //Calculate electron Phi
591  G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
592  G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
593  G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
594  // Trigs
595  cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
596 
597  // End of calculation of ejection Compton electron direction
598 
599  //Fix for floating point errors
600 
601  } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
602 
603  // Revert to original if maximum number of iterations threshold has been reached
604  if (iteration >= maxDopplerIterations)
605  {
606  pERecoil = photonEnergy0 ;
607  bindingE = 0.;
608  dirx=0.0;
609  diry=0.0;
610  dirz=1.0;
611  }
612 
613  // Set "scattered" photon direction and energy
614 
615  G4ThreeVector photonDirection1(dirx,diry,dirz);
616  photonDirection1.rotateUz(photonDirection0);
617  photonPolarization1.rotateUz(photonDirection0);
618  if (pERecoil > 0.)
619  {
620  fParticleChange->SetProposedKineticEnergy(pERecoil) ;
621  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
622  fParticleChange->ProposePolarization(photonPolarization1);
623 
624  // Set ejected Compton electron direction and energy
625  G4double PhiE = std::acos(cosPhiE);
626  G4double eDirX = sinThetaE * std::cos(phi+PhiE);
627  G4double eDirY = sinThetaE * std::sin(phi+PhiE);
628  G4double eDirZ = cosThetaE;
629 
630  G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
631 
632  G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
633  eDirection.rotateUz(photonDirection0);
635  eDirection,eKineticEnergy) ;
636  fvect->push_back(dp);
637 
638  }
639  else
640  {
641  fParticleChange->SetProposedKineticEnergy(0.);
642  fParticleChange->ProposeTrackStatus(fStopAndKill);
643  }
644 
645  // sample deexcitation
646  //
647 
648  if (verboseLevel > 3) {
649  G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
650  }
651 
652  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
653  G4int index = couple->GetIndex();
654  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
655  size_t nbefore = fvect->size();
657  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
658  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
659  size_t nafter = fvect->size();
660  if(nafter > nbefore) {
661  for (size_t i=nbefore; i<nafter; ++i) {
662  //Check if there is enough residual energy
663  if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
664  {
665  //Ok, this is a valid secondary: keep it
666  bindingE -= ((*fvect)[i])->GetKineticEnergy();
667  }
668  else
669  {
670  //Invalid secondary: not enough energy to create it!
671  //Keep its energy in the local deposit
672  delete (*fvect)[i];
673  (*fvect)[i]=0;
674  }
675  }
676  }
677  }
678  }
679 
680  //This should never happen
681  if(bindingE < 0.0)
682  G4Exception("G4LowEPPolarizedComptonModel::SampleSecondaries()",
683  "em2051",FatalException,"Negative local energy deposit");
684 
685  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
686 
687 }
688 
689 //****************************************************************************
690 
691 G4double
692 G4LowEPPolarizedComptonModel::ComputeScatteringFunction(G4double x, G4int Z)
693 {
694  G4double value = Z;
695  if (x <= ScatFuncFitParam[Z][2]) {
696 
697  G4double lgq = G4Log(x)/ln10;
698 
699  if (lgq < ScatFuncFitParam[Z][1]) {
700  value = ScatFuncFitParam[Z][3] + lgq*ScatFuncFitParam[Z][4];
701  } else {
702  value = ScatFuncFitParam[Z][5] + lgq*ScatFuncFitParam[Z][6] +
703  lgq*lgq*ScatFuncFitParam[Z][7] + lgq*lgq*lgq*ScatFuncFitParam[Z][8];
704  }
705  value = G4Exp(value*ln10);
706  }
707  return value;
708 }
709 
710 
711 //****************************************************************************
712 
713 #include "G4AutoLock.hh"
714 namespace { G4Mutex LowEPPolarizedComptonModelMutex = G4MUTEX_INITIALIZER; }
715 
716 void
718  G4int Z)
719 {
720  G4AutoLock l(&LowEPPolarizedComptonModelMutex);
721  if(!data[Z]) { ReadData(Z); }
722  l.unlock();
723 }
724 
725 //****************************************************************************
726 
727 //Fitting data to compute scattering function
728 
729 const G4double G4LowEPPolarizedComptonModel::ScatFuncFitParam[101][9] = {
730 { 0, 0., 0., 0., 0., 0., 0., 0., 0.},
731 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143.374, 50.787, -5.951, 0.2304418},
732 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53.649, 13.892, -0.948, 0.006996759},
733 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62.090, 21.462, -2.453, 0.093416},
734 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127.906, 46.491, -5.614, 0.2262103},
735 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131.153, 47.132, -5.619, 0.2233819},
736 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128.143, 45.379, -5.325, 0.2083009},
737 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131.048, 46.314, -5.421, 0.2114925},
738 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128.225, 44.818, -5.183, 0.1997155},
739 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122.112, 42.103, -4.796, 0.1819099},
740 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110.143, 37.225, -4.143, 0.1532094},
741 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41.137, 12.313, -1.152, 0.03384553},
742 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53.549, 17.420, -1.840, 0.06431849},
743 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66.243, 22.297, -2.460, 0.09045854},
744 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78.271, 26.757, -3.008, 0.1128195},
745 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85.069, 29.164, -3.291, 0.1239113},
746 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93.640, 32.274, -3.665, 0.1388633},
747 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98.534, 33.958, -3.857, 0.1461557},
748 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100.077, 34.379, -3.891, 0.1468902},
749 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53.819, 17.528, -1.851, 0.0648722},
750 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52.221, 17.169, -1.832, 0.06502094},
751 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55.365, 18.276, -1.961, 0.07002778},
752 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57.412, 18.957, -2.036, 0.07278856},
753 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58.478, 19.270, -2.065, 0.07362722},
754 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62.192, 20.358, -2.167, 0.07666583},
755 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58.007, 18.924, -2.003, 0.0704305},
756 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61.176, 20.067, -2.141, 0.0760269},
757 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61.909, 20.271, -2.159, 0.07653559},
758 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62.402, 20.391, -2.167, 0.07664243},
759 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67.305, 21.954, -2.331, 0.0823267},
760 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62.064, 20.136, -2.122, 0.07437589},
761 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61.068, 19.808, -2.086, 0.07307488},
762 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63.126, 20.553, -2.175, 0.07660222},
763 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65.674, 21.445, -2.278, 0.08054694},
764 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69.457, 22.811, -2.442, 0.08709536},
765 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72.283, 23.808, -2.558, 0.09156808},
766 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74.696, 24.641, -2.653, 0.09516597},
767 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46.235, 14.519, -1.458, 0.04837659},
768 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41.784, 13.065, -1.300, 0.04267703},
769 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44.609, 14.114, -1.429, 0.0479348},
770 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47.142, 15.019, -1.536, 0.0521347},
771 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53.114, 17.052, -1.766, 0.06079426},
772 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54.590, 17.550, -1.822, 0.06290335},
773 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51.272, 16.423, -1.694, 0.05806108},
774 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58.314, 18.839, -1.969, 0.0684608},
775 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59.674, 19.295, -2.020, 0.07037294},
776 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72.228, 23.693, -2.532, 0.09017969},
777 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60.890, 19.647, -2.053, 0.07138694},
778 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56.152, 18.002, -1.863, 0.06410123},
779 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56.208, 18.052, -1.872, 0.06456884},
780 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58.478, 18.887, -1.973, 0.06860356},
781 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60.708, 19.676, -2.066, 0.07225841},
782 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63.341, 20.632, -2.180, 0.0767412},
783 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66.339, 21.716, -2.310, 0.08191981},
784 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67.649, 22.151, -2.357, 0.08357825},
785 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45.302, 14.219, -1.423, 0.04712317},
786 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39.825, 12.363, -1.214, 0.03931009},
787 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38.952, 11.982, -1.160, 0.03681554},
788 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41.959, 13.118, -1.302, 0.04271291},
789 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40.528, 12.555, -1.230, 0.03971407},
790 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39.986, 12.329, -1.200, 0.03843737},
791 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42.756, 13.362, -1.327, 0.0436124},
792 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42.682, 13.314, -1.319, 0.04322509},
793 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42.399, 13.185, -1.301, 0.04243861},
794 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43.226, 13.475, -1.335, 0.04377341},
795 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43.232, 13.456, -1.330, 0.04347536},
796 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42.047, 12.990, -1.270, 0.04095499},
797 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42.405, 13.106, -1.283, 0.04146024},
798 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41.974, 12.926, -1.259, 0.040435},
799 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42.132, 12.967, -1.262, 0.04048908},
800 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42.582, 13.122, -1.280, 0.04119599},
801 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42.586, 13.115, -1.278, 0.04107587},
802 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43.708, 13.509, -1.324, 0.04286491},
803 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44.838, 13.902, -1.369, 0.04457132},
804 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45.545, 14.137, -1.395, 0.04553459},
805 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46.426, 14.431, -1.428, 0.04678218},
806 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47.513, 14.806, -1.471, 0.04842566},
807 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48.225, 15.042, -1.497, 0.04938364},
808 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53.256, 16.739, -1.687, 0.05645173},
809 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53.900, 16.946, -1.709, 0.05723134},
810 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49.801, 15.536, -1.547, 0.05103522},
811 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49.651, 15.512, -1.548, 0.05123203},
812 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51.021, 16.018, -1.609, 0.05364831},
813 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52.693, 16.612, -1.679, 0.05638698},
814 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54.415, 17.238, -1.754, 0.05935566},
815 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56.083, 17.834, -1.824, 0.06206068},
816 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57.860, 18.463, -1.898, 0.0649633},
817 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39.973, 12.164, -1.162, 0.0364598},
818 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34.591, 10.338, -0.956, 0.0287409},
819 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36.002, 10.867, -1.021, 0.03136835},
820 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37.660, 11.475, -1.095, 0.03435334},
821 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37.268, 11.301, -1.071, 0.0330539},
822 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37.695, 11.438, -1.085, 0.03376669},
823 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39.010, 11.927, -1.146, 0.03630848},
824 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37.192, 11.229, -1.057, 0.0325621},
825 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37.589, 11.363, -1.072, 0.03312393},
826 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39.573, 12.095, -1.162, 0.03680527},
827 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40.007, 12.242, -1.178, 0.03737377},
828 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39.509, 12.041, -1.152, 0.03629023},
829 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39.939, 12.183, -1.168, 0.03690464},
830 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39.973, 12.180, -1.166, 0.036773}
831  };
832 
833 //****************************************************************************
834 
835 //Supporting functions for photon polarisation effects
836 
837 G4double G4LowEPPolarizedComptonModel::SetPhi(G4double energyRate,
838  G4double sinT2)
839 {
840  G4double rand1;
841  G4double rand2;
842  G4double phiProbability;
843  G4double phi;
844  G4double a, b;
845 
846  do
847  {
848  rand1 = G4UniformRand();
849  rand2 = G4UniformRand();
850  phiProbability=0.;
851  phi = twopi*rand1;
852 
853  a = 2*sinT2;
854  b = energyRate + 1/energyRate;
855 
856  phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
857 
858 
859 
860  }
861  while ( rand2 > phiProbability );
862  return phi;
863 }
864 
865 //****************************************************************************
866 
867 G4ThreeVector G4LowEPPolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
868 {
869  G4double dx = a.x();
870  G4double dy = a.y();
871  G4double dz = a.z();
872  G4double x = dx < 0.0 ? -dx : dx;
873  G4double y = dy < 0.0 ? -dy : dy;
874  G4double z = dz < 0.0 ? -dz : dz;
875  if (x < y) {
876  return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
877  }else{
878  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
879  }
880 }
881 
882 //****************************************************************************
883 
884 G4ThreeVector G4LowEPPolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
885 {
886  G4ThreeVector d0 = direction0.unit();
887  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
888  G4ThreeVector a0 = a1.unit(); // unit vector
889 
890  G4double rand1 = G4UniformRand();
891 
892  G4double angle = twopi*rand1; // random polar angle
893  G4ThreeVector b0 = d0.cross(a0); // cross product
894 
896 
897  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
898  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
899  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
900 
901  G4ThreeVector c0 = c.unit();
902 
903  return c0;
904 
905 }
906 
907 //****************************************************************************
908 
909 G4ThreeVector G4LowEPPolarizedComptonModel::GetPerpendicularPolarization
910 (const G4ThreeVector& photonDirection, const G4ThreeVector& photonPolarization) const
911 {
912 
913  //
914  // The polarization of a photon is always perpendicular to its momentum direction.
915  // Therefore this function removes those vector component of photonPolarization, which
916  // points in direction of photonDirection
917  //
918  // Mathematically we search the projection of the vector a on the plane E, where n is the
919  // plains normal vector.
920  // The basic equation can be found in each geometry book (e.g. Bronstein):
921  // p = a - (a o n)/(n o n)*n
922 
923  return photonPolarization - photonPolarization.dot(photonDirection)/photonDirection.dot(photonDirection) * photonDirection;
924 }
925 
926 //****************************************************************************
927 
928 G4ThreeVector G4LowEPPolarizedComptonModel::SetNewPolarization(G4double LowEPPCepsilon,
929  G4double sinT2,
930  G4double phi,
931  G4double costheta)
932 {
933  G4double rand1;
934  G4double rand2;
935  G4double cosPhi = std::cos(phi);
936  G4double sinPhi = std::sin(phi);
937  G4double sinTheta = std::sqrt(sinT2);
938  G4double cosP2 = cosPhi*cosPhi;
939  G4double normalisation = std::sqrt(1. - cosP2*sinT2);
940 
941 
942  // Method based on:
943  // D. Xu, Z. He and F. Zhang
944  // "Detection of Gamma Ray Polarization Using a 3-D Position Sensitive CdZnTe Detector"
945  // IEEE TNS, Vol. 52(4), 1160-1164, 2005.
946 
947  // Determination of Theta
948 
949  G4double theta;
950 
951  rand1 = G4UniformRand();
952  rand2 = G4UniformRand();
953 
954  if (rand1<(LowEPPCepsilon+1.0/LowEPPCepsilon-2)/(2.0*(LowEPPCepsilon+1.0/LowEPPCepsilon)-4.0*sinT2*cosP2))
955  {
956  if (rand2<0.5)
957  theta = pi/2.0;
958  else
959  theta = 3.0*pi/2.0;
960  }
961  else
962  {
963  if (rand2<0.5)
964  theta = 0;
965  else
966  theta = pi;
967  }
968  G4double cosBeta = std::cos(theta);
969  G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
970 
971  G4ThreeVector photonPolarization1;
972 
973  G4double xParallel = normalisation*cosBeta;
974  G4double yParallel = -(sinT2*cosPhi*sinPhi)*cosBeta/normalisation;
975  G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
976  G4double xPerpendicular = 0.;
977  G4double yPerpendicular = (costheta)*sinBeta/normalisation;
978  G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
979 
980  G4double xTotal = (xParallel + xPerpendicular);
981  G4double yTotal = (yParallel + yPerpendicular);
982  G4double zTotal = (zParallel + zPerpendicular);
983 
984  photonPolarization1.setX(xTotal);
985  photonPolarization1.setY(yTotal);
986  photonPolarization1.setZ(zTotal);
987 
988  return photonPolarization1;
989 
990 }
void SetOccupancyData()
Definition: G4ShellData.hh:70
const G4double a0
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
double x() const
double dot(const Hep3Vector &) const
float h_Planck
Definition: hepunit.py:263
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetZ() const
Definition: G4Element.hh:131
static G4double angle[DIM]
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
size_t GetVectorLength() const
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void setY(double)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
G4int SelectRandomShell(G4int Z) const
Definition: G4ShellData.cc:363
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
double z() const
void setZ(double)
void setX(double)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
void LoadData(const G4String &fileName)
Definition: G4ShellData.cc:234
static const G4double ln10
const XML_Char const XML_Char * data
Definition: expat.h:268
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
string material
Definition: eplot.py:19
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
const XML_Char * s
Definition: expat.h:262
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static constexpr double m
Definition: G4SIunits.hh:129
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ThreeVector & GetMomentumDirection() const
G4double BindingEnergy(G4int Z, G4int shellIndex) const
Definition: G4ShellData.cc:166
static constexpr double cm
Definition: G4SIunits.hh:119
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void ProposePolarization(const G4ThreeVector &dir)
static constexpr double eV
Definition: G4SIunits.hh:215
float electron_mass_c2
Definition: hepunit.py:274
const G4int n
static constexpr double kg
Definition: G4SIunits.hh:182
G4double Energy(size_t index) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4LowEPPolarizedComptonModel(const G4ParticleDefinition *p=0, const G4String &nam="LowEPComptonModel")
static G4ProductionCutsTable * GetProductionCutsTable()
G4int G4Mutex
Definition: G4Threading.hh:173
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:817
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
Hep3Vector unit() const
double y() const
const G4ThreeVector & GetPolarization() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static constexpr double GeV
Definition: G4SIunits.hh:217
tuple z
Definition: test.py:28
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
Hep3Vector cross(const Hep3Vector &) const
static constexpr double halfpi
Definition: G4SIunits.hh:77
G4VAtomDeexcitation * AtomDeexcitation()
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
tuple c
Definition: test.py:13
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788
static constexpr double barn
Definition: G4SIunits.hh:105
double mag() const
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
float c_light
Definition: hepunit.py:257
G4AtomicShellEnumerator
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
const G4Material * GetMaterial() const