132 G4cout <<
"Old Momentum Direction: "
134 G4cout <<
"Old Polarization: "
143 G4double CosTheta, SinTheta, SinPhi, CosPhi, unit_x, unit_y, unit_z;
150 SinTheta = std::sqrt(1.-CosTheta*CosTheta);
156 SinPhi = std::sin(rand);
157 CosPhi = std::cos(rand);
160 unit_x = SinTheta * CosPhi;
161 unit_y = SinTheta * SinPhi;
163 NewMomentumDirection.set (unit_x,unit_y,unit_z);
167 OldMomentumDirection = OldMomentumDirection.unit();
168 NewMomentumDirection.rotateUz(OldMomentumDirection);
169 NewMomentumDirection = NewMomentumDirection.unit();
175 constant = -NewMomentumDirection.dot(OldPolarization);
177 NewPolarization = OldPolarization + constant*NewMomentumDirection;
178 NewPolarization = NewPolarization.unit();
183 if (NewPolarization.mag() == 0.) {
185 NewPolarization.set(std::cos(rand),std::sin(rand),0.);
186 NewPolarization.rotateUz(NewMomentumDirection);
190 if (
G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
194 cosTheta = NewPolarization.dot(OldPolarization);
202 G4cout <<
"New Polarization: "
203 << NewPolarization <<
G4endl;
204 G4cout <<
"Polarization Change: "
206 G4cout <<
"New Momentum Direction: "
207 << NewMomentumDirection <<
G4endl;
208 G4cout <<
"Momentum Change: "
230 for(
G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
232 G4Material* material = (*theMaterialTable)[iMaterial];
236 if ( materialProperties != NULL ) {
237 rayleigh = materialProperties->
GetProperty(
"RAYLEIGH" );
238 if ( rayleigh == NULL ) rayleigh =
258 ((*thePhysicsTable)(material->
GetIndex()));
261 if( rayleigh != NULL ) rsLength = rayleigh->
Value( photonMomentum );
278 if ( material->
GetName() ==
"Water" )
279 betat = 7.658e-23*
m3/
MeV;
287 if ( rIndex == NULL )
return NULL;
297 if( material->
GetName() ==
"Water" )
298 temperature = 283.15*
kelvin;
305 const G4double c1 = scaleFactor * betat * temperature * k_Boltzmann /
308 for(
size_t uRIndex = 0; uRIndex < rIndex->
GetVectorLength(); uRIndex++ )
311 const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
315 std::pow(((rIndexSquared-1.0)*(rIndexSquared+2.0 )/3.0),2);
317 const G4double meanFreePath = 1.0 / ( c1 * c2 *
c3 );
320 G4cout << energy <<
"MeV\t" << meanFreePath <<
"mm" <<
G4endl;
322 rayleighMeanFreePaths->
InsertValues( energy, meanFreePath );
325 return rayleighMeanFreePaths;
static c2_factory< G4double > c2
G4MaterialPropertyVector * GetProperty(const char *key)
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
const G4String & GetName() const
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
void InsertValues(G4double energy, G4double value)
const G4ThreeVector * GetMomentumDirection() const
size_t GetVectorLength() const
G4PhysicsTable * thePhysicsTable
G4double GetTotalMomentum() const
G4GLOB_DLL std::ostream G4cout
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
void SetProcessSubType(G4int)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
static size_t GetNumberOfMaterials()
const G4String & GetProcessName() const
G4Material * GetMaterial() const
static const double kelvin
virtual void Initialize(const G4Track &)
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
G4double energy(const ThreeVector &p, const G4double m)
G4OpRayleigh(const G4String &processName="OpRayleigh", G4ProcessType type=fOptical)
const G4ThreeVector & GetPolarization() const
G4ParticleChange aParticleChange
G4bool ConstPropertyExists(const char *key)
void insertAt(size_t, G4PhysicsVector *)
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
G4double GetTemperature() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4PhysicsOrderedFreeVector * CalculateRayleighMeanFreePaths(const G4Material *material) const
Calculates the mean free paths for a material as a function of photon energy.
const G4ThreeVector * GetPolarization() const
G4double GetConstProperty(const char *key)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)