66 scatterFunctionData(0),
67 crossSectionHandler(0),fAtomDeexcitation(0)
78 G4cout <<
"Livermore Modified Compton model is constructed " <<
G4endl;
99 G4cout <<
"Calling G4LivermoreComptonModifiedModel::Initialise()" <<
G4endl;
111 G4String crossSectionFile =
"comp/ce-cs-";
115 G4String scatterFile =
"comp/ce-sf-";
127 G4cout <<
"Loaded cross section files for Livermore Modified Compton model" <<
G4endl;
138 G4cout <<
"Livermore modified Compton model is initialized " <<
G4endl
155 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" <<
G4endl;
189 G4cout <<
"G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
207 G4double epsilon0Local = 1. / (1. + 2. * e0m);
208 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
209 G4double alpha1 = -std::log(epsilon0Local);
210 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
231 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
232 epsilon = std::sqrt(epsilonSq);
235 oneCosT = (1. -
epsilon) / ( epsilon * e0m);
236 sinT2 = oneCosT * (2. - oneCosT);
237 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/
cm);
239 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
244 G4double sinTheta = std::sqrt (sinT2);
246 G4double dirx = sinTheta * std::cos(phi);
247 G4double diry = sinTheta * std::sin(phi);
256 G4int maxDopplerIterations = 1000;
258 G4double photonEoriginal = epsilon * photonEnergy0;
265 G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
266 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
291 }
while(Alpha >= (
pi/2.0));
293 ePAU = pSample / std::cos(Alpha);
297 G4double ePSI = ePAU * momentum_au_to_nat;
298 G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
302 systemE = eEIncident+photonEnergy0;
306 G4double pDoppler2 = pDoppler * pDoppler;
308 G4double var3 = var2*var2 - pDoppler2;
309 G4double var4 = var2 - pDoppler2 * cosTheta;
310 G4double var = var4*var4 - var3 + pDoppler2 * var3;
316 if (
G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
317 else { photonE = (var4 + varSqrt) * scale; }
323 }
while ( iteration <= maxDopplerIterations &&
324 (photonE < 0. || photonE > eMax ) );
335 if(eKineticEnergy < 0.0) {
336 G4cout <<
"Error, kinetic energy of electron less than zero" <<
G4endl;
343 G4double E_num = photonEnergy0 - photonE*cosTheta;
344 G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
346 G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
348 eDirX = sinThetaE * std::cos(phi);
349 eDirY = sinThetaE * std::sin(phi);
353 eDirection.
rotateUz(photonDirection0);
355 eDirection,eKineticEnergy) ;
356 fvect->push_back(dp);
362 if (iteration >= maxDopplerIterations)
364 photonE = photonEoriginal;
371 photonDirection1.
rotateUz(photonDirection0);
376 if (photonEnergy1 > 0.)
380 if (iteration < maxDopplerIterations)
383 eDirection.
rotateUz(photonDirection0);
385 eDirection,eKineticEnergy) ;
386 fvect->push_back(dp);
401 size_t nbefore = fvect->size();
405 size_t nafter = fvect->size();
406 if(nafter > nbefore) {
407 for (
size_t i=nbefore; i<nafter; ++i) {
408 bindingE -= ((*fvect)[i])->GetKineticEnergy();
413 if(bindingE < 0.0) { bindingE = 0.0; }