Skip to content

Instantly share code, notes, and snippets.

@pgjones

pgjones/G4OpRayleigh.cc

Last active August 29, 2015 14:05
Show Gist options
  • Save pgjones/738fa2a00c6e16d2ecbb to your computer and use it in GitHub Desktop.
Save pgjones/738fa2a00c6e16d2ecbb to your computer and use it in GitHub Desktop.
Download ZIP
Proposed new G4OpRayleigh process for Geant4
Raw
G4OpRayleigh.cc
//
// ********************************************************************
// * 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: G4OpRayleigh.cc 82117 2014-06-11 09:10:52Z gcosmo $
//
//
////////////////////////////////////////////////////////////////////////
// Optical Photon Rayleigh Scattering Class Implementation
////////////////////////////////////////////////////////////////////////
//
// File: G4OpRayleigh.cc
// Description: Discrete Process -- Rayleigh scattering of optical
// photons
// Version: 1.0
// Created: 1996-05-31
// Author: Juliet Armstrong
// Updated: 2010-06-11 - Fix Bug 207; Thanks to Xin Qian
// (Kellogg Radiation Lab of Caltech)
// 2005-07-28 - add G4ProcessType to constructor
// 2001-10-18 by Peter Gumplinger
// eliminate unused variable warning on Linux (gcc-2.95.2)
// 2001-09-18 by mma
// >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
// 2001-01-30 by Peter Gumplinger
// > allow for positiv and negative CosTheta and force the
// > new momentum direction to be in the same plane as the
// > new and old polarization vectors
// 2001-01-29 by Peter Gumplinger
// > fix calculation of SinTheta (from CosTheta)
// 1997-04-09 by Peter Gumplinger
// > new physics/tracking scheme
// mail: [email protected]
//
////////////////////////////////////////////////////////////////////////
#include "G4OpRayleigh.hh"
#include "G4ios.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4OpProcessSubType.hh"
/////////////////////////
// Class Implementation
/////////////////////////
//////////////
// Operators
//////////////
// G4OpRayleigh::operator=(const G4OpRayleigh &right)
// {
// }
/////////////////
// Constructors
/////////////////
G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type)
: G4VDiscreteProcess(processName, type)
{
SetProcessSubType(fOpRayleigh);
thePhysicsTable = NULL;
if (verboseLevel>0) {
G4cout << GetProcessName() << " is created " << G4endl;
}
}
// G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
// {
// }
////////////////
// Destructors
////////////////
G4OpRayleigh::~G4OpRayleigh()
{
if (thePhysicsTable) {
thePhysicsTable->clearAndDestroy();
delete thePhysicsTable;
}
}
////////////
// Methods
////////////
// PostStepDoIt
// -------------
//
G4VParticleChange*
G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
{
aParticleChange.Initialize(aTrack);
const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
if (verboseLevel>0) {
G4cout << "Scattering Photon!" << G4endl;
G4cout << "Old Momentum Direction: "
<< aParticle->GetMomentumDirection() << G4endl;
G4cout << "Old Polarization: "
<< aParticle->GetPolarization() << G4endl;
}
G4double cosTheta;
G4ThreeVector OldMomentumDirection, NewMomentumDirection;
G4ThreeVector OldPolarization, NewPolarization;
do {
// Try to simulate the scattered photon momentum direction
// w.r.t. the initial photon momentum direction
G4double CosTheta = G4UniformRand();
G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
// consider for the angle 90-180 degrees
if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
// simulate the phi angle
G4double rand = twopi*G4UniformRand();
G4double SinPhi = std::sin(rand);
G4double CosPhi = std::cos(rand);
// start constructing the new momentum direction
G4double unit_x = SinTheta * CosPhi;
G4double unit_y = SinTheta * SinPhi;
G4double unit_z = CosTheta;
NewMomentumDirection.set (unit_x,unit_y,unit_z);
// Rotate the new momentum direction into global reference system
OldMomentumDirection = aParticle->GetMomentumDirection();
OldMomentumDirection = OldMomentumDirection.unit();
NewMomentumDirection.rotateUz(OldMomentumDirection);
NewMomentumDirection = NewMomentumDirection.unit();
// calculate the new polarization direction
// The new polarization needs to be in the same plane as the new
// momentum direction and the old polarization direction
OldPolarization = aParticle->GetPolarization();
G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
NewPolarization = NewMomentumDirection + constant*OldPolarization;
NewPolarization = NewPolarization.unit();
// There is a corner case, where the Newmomentum direction
// is the same as oldpolariztion direction:
// random generate the azimuthal angle w.r.t. Newmomentum direction
if (NewPolarization.mag() == 0.) {
rand = G4UniformRand()*twopi;
NewPolarization.set(std::cos(rand),std::sin(rand),0.);
NewPolarization.rotateUz(NewMomentumDirection);
} else {
// There are two directions which are perpendicular
// to the new momentum direction
if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
}
// simulate according to the distribution cos^2(theta)
cosTheta = NewPolarization.dot(OldPolarization);
} while (std::pow(cosTheta,2) < G4UniformRand());
aParticleChange.ProposePolarization(NewPolarization);
aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
if (verboseLevel>0) {
G4cout << "New Polarization: "
<< NewPolarization << G4endl;
G4cout << "Polarization Change: "
<< *(aParticleChange.GetPolarization()) << G4endl;
G4cout << "New Momentum Direction: "
<< NewMomentumDirection << G4endl;
G4cout << "Momentum Change: "
<< *(aParticleChange.GetMomentumDirection()) << G4endl;
}
return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
// BuildPhysicsTable for the Rayleigh Scattering process
// --------------------------------------------------------
void G4OpRayleigh::BuildPhysicsTable(const G4ParticleDefinition&)
{
if (thePhysicsTable) {
thePhysicsTable->clearAndDestroy();
delete thePhysicsTable;
thePhysicsTable = NULL;
}
const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
fPhysicsTable = new G4PhysicsTable( numOfMaterials );
for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
{
G4Material* material = (*theMaterialTable)[iMaterial];
G4MaterialPropertiesTable* materialProperties = material->GetMaterialPropertiesTable();
G4PhysicsOrderedFreeVector* rayleigh = NULL;
if( materialProperties != NULL )
{
rayleigh = materialProperties->GetProperty( "RAYLEIGH" );
if( rayleigh == NULL )
rayleigh = CalculateRayleighMeanFreePaths( material );
}
fPhysicsTable->insertAt( iMaterial, rayleigh );
}
}
// GetMeanFreePath()
// -----------------
//
G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
G4double ,
G4ForceCondition* )
{
const G4DynamicParticle* particle = track.GetDynamicParticle();
const G4double photonMomentum = particle->GetTotalMomentum();
const G4Material* material = track.GetMaterial();
G4PhysicsOrderedFreeVector* rayleigh = static_cast<G4PhysicsOrderedFreeVector*>( (*fPhysicsTable)( material->GetIndex() ) );
G4double rsLength = DBL_MAX;
if( rayleigh != NULL )
rsLength = rayleigh->Value( photonMomentum );
return rsLength;
}
// CalculateRayleighMeanFreePaths()
// --------------------------------
// Private method to compute Rayleigh Scattering Lengths
G4PhysicsOrderedFreeVector*
G4OpRayleigh::CalculateRayleighMeanFreePaths( G4Material* const material ) const
{
G4MaterialPropertiesTable* materialProperties = material->GetMaterialPropertiesTable();
// Retrieve the beta_T or isothermal compressibility value. For backwards
// compatibility use a constant if the material is "Water". If the material
// doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
G4double betat;
if( material->GetName() == "Water" )
betat = 7.658e-23*m3/MeV;
else if( materialProperties->ConstPropertyExists( "ISOTHERMAL_COMPRESSIBILITY" ) )
betat = materialProperties->GetConstProperty( "ISOTHERMAL_COMPRESSIBILITY" );
else
return NULL;
// If the material doesn't have a RINDEX property vector then return
G4MaterialPropertyVector* rIndex = materialProperties->GetProperty("RINDEX");
if( rIndex == NULL )
return NULL;
// Retrieve the optional scale factor, (this just scales the scattering length
G4double scaleFactor = 1.0;
if( materialProperties->ConstPropertyExists( "RS_SCALE_FACTOR" ) )
scaleFactor= materialProperties->GetConstProperty("RS_SCALE_FACTOR" );
// Retrieve the material temperature. For backwards compatibility use a
// constant if the material is "Water"
G4double temperature;
if( material->GetName() == "Water" )
temperature = 283.15*kelvin; // Temperature of water is 10 degrees celsius
else
material->GetTemperature();
G4PhysicsOrderedFreeVector* rayleighMeanFreePaths = new G4PhysicsOrderedFreeVector();
// This calculates the meanFreePath via the Einstein-Smoluchowski formula
const G4double c1 = scaleFactor * betat * temperature * k_Boltzmann / ( 6.0 * pi );
for( size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); uRIndex++ )
{
const G4double energy = rIndex->Energy( uRIndex );
const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
const G4double xlambda = h_Planck * c_light / energy;
const G4double c2 = std::pow( twopi / xlambda, 4 );
const G4double c3 = std::pow( ( ( rIndexSquared - 1.0 ) * ( rIndexSquared + 2.0 ) / 3.0 ), 2 );
const G4double meanFreePath = 1.0 / ( c1 * c2 * c3 );
if( verboseLevel>0 )
G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl;
rayleighMeanFreePaths->InsertValues( energy, meanFreePath );
}
return rayleighMeanFreePaths;
}
Raw
G4OpRayleigh.hh
//
// ********************************************************************
// * 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: G4OpRayleigh.hh 82117 2014-06-11 09:10:52Z gcosmo $
//
//
////////////////////////////////////////////////////////////////////////
// Optical Photon Rayleigh Scattering Class Definition
////////////////////////////////////////////////////////////////////////
//
// File: G4OpRayleigh.hh
// Description: Discrete Process -- Rayleigh scattering of optical photons
// Version: 1.0
// Created: 1996-05-31
// Author: Juliet Armstrong
// Updated: 2014-08-20 allow for more material types
// 2005-07-28 add G4ProcessType to constructor
// 1999-10-29 add method and class descriptors
// 1997-04-09 by Peter Gumplinger
// > new physics/tracking scheme
// mail: [email protected]
//
////////////////////////////////////////////////////////////////////////
#ifndef G4OpRayleigh_h
#define G4OpRayleigh_h 1
/////////////
// Includes
/////////////
#include "globals.hh"
#include "templates.hh"
#include "Randomize.hh"
#include "G4ThreeVector.hh"
#include "G4ParticleMomentum.hh"
#include "G4Step.hh"
#include "G4VDiscreteProcess.hh"
#include "G4DynamicParticle.hh"
#include "G4Material.hh"
#include "G4OpticalPhoton.hh"
#include "G4PhysicsTable.hh"
#include "G4PhysicsOrderedFreeVector.hh"
// Class Description:
// Discrete Process -- Rayleigh scattering of optical photons.
// Class inherits publicly from G4VDiscreteProcess.
// Class Description - End:
/////////////////////
// Class Definition
/////////////////////
class G4OpRayleigh : public G4VDiscreteProcess
{
public:
////////////////////////////////
// Constructors and Destructor
////////////////////////////////
G4OpRayleigh(const G4String& processName = "OpRayleigh",
G4ProcessType type = fOptical);
~G4OpRayleigh();
private:
G4OpRayleigh(const G4OpRayleigh &right);
//////////////
// Operators
//////////////
G4OpRayleigh& operator=(const G4OpRayleigh &right);
public:
////////////
// Methods
////////////
G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
// Returns true -> 'is applicable' only for an optical photon.
void BuildPhysicsTable(const G4ParticleDefinition& aParticleType);
// Build thePhysicsTable at a right time
G4double GetMeanFreePath(const G4Track& aTrack,
G4double ,
G4ForceCondition* );
// Returns the mean free path for Rayleigh scattering in water.
// --- Not yet implemented for other materials! ---
G4VParticleChange* PostStepDoIt(const G4Track& aTrack,
const G4Step& aStep);
// This is the method implementing Rayleigh scattering.
G4PhysicsTable* GetPhysicsTable() const;
// Returns the address of the physics table.
void DumpPhysicsTable() const;
// Prints the physics table.
private:
/////////////////////
// Helper Functions
/////////////////////
/// Calculates the mean free paths for a material as a function of
/// photon energy
///
/// @param[in] material information
/// @return the mean free path vector
G4PhysicsOrderedFreeVector*
CalculateRayleighMeanFreePaths( G4Material* const material ) const;
///////////////////////
// Class Data Members
///////////////////////
protected:
G4PhysicsTable* thePhysicsTable;
// A Physics Table can be either a cross-sections table or
// an energy table (or can be used for other specific
// purposes).
private:
};
////////////////////
// Inline methods
////////////////////
inline
G4bool G4OpRayleigh::IsApplicable(const G4ParticleDefinition& aParticleType)
{
return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
}
inline
void G4OpRayleigh::DumpPhysicsTable() const
{
G4int PhysicsTableSize = thePhysicsTable->entries();
G4PhysicsOrderedFreeVector *v;
for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
{
v = (G4PhysicsOrderedFreeVector*)(*thePhysicsTable)[i];
v->DumpValues();
}
}
inline G4PhysicsTable* G4OpRayleigh::GetPhysicsTable() const
{
return thePhysicsTable;
}
#endif /* G4OpRayleigh_h */
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment