// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // $Id: G4PolarizedComptonScattering.cc,v 1.18 2008/10/15 17:53:44 vnivanch Exp $ // GEANT4 tag $Name: geant4-09-03 $ // // //---------- G4PolarizedComptonScattering physics process ---------------------- // by Vicente Lara, March 1998 // // ----------------------------------------------------------------------------- // Corrections by Rui Curado da Silva (Nov. 2000) // - Sampling of Phi // - Depolarization probability // // 13-07-01, DoIt: suppression of production cut for the electron (mma) // 20-09-01, DoIt: fminimalEnergy = 1*eV (mma) // 04-05-05, Inheritance from ComptonScattering52 (V.Ivanchenko) // 30-01-06, DoIt : return G4ComptonScattering52::PostStepDoIt(aTrack,aStep) mma // // ----------------------------------------------------------------------------- #include "G4PolarizedComptonScattering.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... using namespace std; G4PolarizedComptonScattering::G4PolarizedComptonScattering(const G4String& pname) : G4ComptonScattering52 (pname) { } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4VParticleChange* G4PolarizedComptonScattering::PostStepDoIt( const G4Track& aTrack, const G4Step& aStep) // // The scattered gamma energy is sampled according to Klein - Nishina formula. // The random number techniques of Butcher & Messel are used // (Nuc Phys 20(1960),15). // GEANT4 internal units // // Note : Effects due to binding of atomic electrons are negliged. { aParticleChange.Initialize(aTrack); const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle(); G4ThreeVector GammaPolarization0 = aDynamicGamma->GetPolarization(); if (std::abs(GammaPolarization0.mag() - 1.e0) > 1.e-14) return G4ComptonScattering52::PostStepDoIt(aTrack,aStep); G4double GammaEnergy0 = aDynamicGamma->GetKineticEnergy(); G4double E0_m = GammaEnergy0 / electron_mass_c2; G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection(); // // sample the energy rate of the scattered gamma // G4double epsilon, epsilonsq, onecost, sint2, greject; G4double epsilon0 = 1./(1. + 2*E0_m) , epsilon0sq = epsilon0*epsilon0; G4double alpha1 = - log(epsilon0) , alpha2 = 0.5*(1.- epsilon0sq); do { if (alpha1/(alpha1+alpha2) > G4UniformRand()) { epsilon = exp(-alpha1*G4UniformRand()); // epsilon0**r epsilonsq = epsilon*epsilon; } else { epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand(); epsilon = sqrt(epsilonsq); }; onecost = (1.- epsilon)/(epsilon*E0_m); sint2 = onecost*(2.-onecost); greject = 1. - epsilon*sint2/(1.+ epsilonsq); } while (greject < G4UniformRand()); // // Phi determination // G4double minimum=0., maximum=twopi, middle=0., resolution=0.001; G4double Rand = G4UniformRand(); int j = 0; while ((j < 100) && (std::abs(SetPhi(epsilon,sint2,middle,Rand)) > resolution)) { middle = (maximum + minimum)/2; if (SetPhi(epsilon,sint2,middle,Rand)* SetPhi(epsilon,sint2,minimum,Rand)<0) maximum = middle; else minimum = middle; j++; } // // scattered gamma angles. ( Z - axis along the parent gamma) // G4double cosTeta = 1. - onecost , sinTeta = sqrt (sint2); G4double Phi = middle; G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta; // // update G4VParticleChange for the scattered gamma // G4double GammaEnergy1 = epsilon*GammaEnergy0; // New polarization // G4ThreeVector GammaPolarization1 = SetNewPolarization(epsilon,sint2,Phi, cosTeta, GammaPolarization0); // Set new direction G4ThreeVector GammaDirection1 ( dirx,diry,dirz ); // Change reference frame. SystemOfRefChange(GammaDirection0,GammaDirection1, GammaPolarization0,GammaPolarization1); G4double localEnergyDeposit = 0.; if (GammaEnergy1 > fminimalEnergy) { aParticleChange.ProposeEnergy(GammaEnergy1); } else { localEnergyDeposit += GammaEnergy1; aParticleChange.ProposeEnergy(0.) ; aParticleChange.ProposeTrackStatus(fStopAndKill); } // // kinematic of the scattered electron // G4double ElecKineEnergy = GammaEnergy0 - GammaEnergy1; if (ElecKineEnergy > fminimalEnergy) { G4double ElecMomentum = sqrt(ElecKineEnergy* (ElecKineEnergy+2.*electron_mass_c2)); G4ThreeVector ElecDirection ( (GammaEnergy0*GammaDirection0 - GammaEnergy1*GammaDirection1) *(1./ElecMomentum)); // create G4DynamicParticle object for the electron. G4DynamicParticle* aElectron= new G4DynamicParticle ( G4Electron::Electron(),ElecDirection,ElecKineEnergy); aParticleChange.SetNumberOfSecondaries(1); aParticleChange.AddSecondary( aElectron ); } else { aParticleChange.SetNumberOfSecondaries(0); localEnergyDeposit += ElecKineEnergy; } aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit); // Reset NbOfInteractionLengthLeft and return aParticleChange return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4double G4PolarizedComptonScattering::SetPhi(G4double EnergyRate, G4double sinsqrth, G4double phi, G4double rand) { G4double cosphi = cos(phi), sinphi = sin(phi); G4double PhiDetermination = ((twopi*rand - phi) *(EnergyRate + 1./EnergyRate - sinsqrth)) + (sinsqrth*sinphi*cosphi); return PhiDetermination; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4ThreeVector G4PolarizedComptonScattering::SetNewPolarization( G4double EnergyRate, G4double sinsqrth, G4double phi, G4double costheta, G4ThreeVector&) { G4double cosphi = cos(phi), sinphi = sin(phi); //// G4double ParallelIntensityPolar = EnergyRate + 1./EnergyRate //// + 2. - 4.*sinsqrth*cosphi*cosphi; G4double ParallelIntensityPolar = EnergyRate + 1./EnergyRate - 2.*sinsqrth*cosphi*cosphi; G4double PerpendiIntensityPolar = EnergyRate + 1./EnergyRate - 2.; G4double PolarizationDegree = sqrt(sinsqrth*sinphi*sinphi+costheta*costheta); G4double sintheta = sqrt(sinsqrth); G4ThreeVector GammaPolarization1; // depolarization probability (1-P) if ( G4UniformRand() > (PerpendiIntensityPolar/ParallelIntensityPolar) ) { // Parallel to initial polarization GammaPolarization1.setX(PolarizationDegree); GammaPolarization1.setY(-sinsqrth*sinphi*cosphi/PolarizationDegree); GammaPolarization1.setZ(-sintheta*costheta*cosphi/PolarizationDegree); } else { // Perpendicular to initial polarization GammaPolarization1.setX(0.); GammaPolarization1.setY(costheta/PolarizationDegree); GammaPolarization1.setZ(-sintheta*sinphi/PolarizationDegree); }; return GammaPolarization1; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void G4PolarizedComptonScattering::SystemOfRefChange(G4ThreeVector& Direction0, G4ThreeVector& Direction1, G4ThreeVector& Polarization0, G4ThreeVector& Polarization1) { // Angles for go back to the original RS G4double cosTeta0 = Direction0.cosTheta(), sinTeta0 = sin(Direction0.theta()); G4double cosPhi0 = cos(Direction0.phi()), sinPhi0 = sin(Direction0.phi()); G4double cosPsi, sinPsi; if (sinTeta0 != 0. ) { cosPsi = -Polarization0.z()/sinTeta0; if (cosPhi0 != 0.) sinPsi = (Polarization0.y() - cosTeta0*sinPhi0*cosPsi)/cosPhi0; else sinPsi = -Polarization0.x()/sinPhi0; } else { cosPsi = Polarization0.x()/cosTeta0; sinPsi = Polarization0.y(); } G4double Psi = atan(sinPsi/cosPsi); // Rotation along Z axe Direction1.rotateZ(Psi); // Direction1.rotateUz(Direction0); aParticleChange.ProposeMomentumDirection(Direction1); // 3 Euler angles rotation for scattered photon polarization Polarization1.rotateZ(Psi); Polarization1.rotateUz(Direction0); aParticleChange.ProposePolarization(Polarization1); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......