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