// // ******************************************************************** // * 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: G4OpMie.cc,v 1.14 2008/07/21 21:08:54 gunter Exp $ // GEANT4 tag $Name: geant4-09-01-patch-01 $ // // //////////////////////////////////////////////////////////////////////// // Optical Photon Mie Scattering Class Implementation //////////////////////////////////////////////////////////////////////// // // File: G4OpMie.cc // Description: Discrete Process -- Mie scattering of optical // photons // Phase Functions taken from Vladimir I. Haltrin, // http://dx.doi.org/10.1117/12.558313 // Version: 1.0 // Created: 2008-07-21 // Author: Apostolos Tsirigotis // Updated: // // mail: tsirigotis@eap.gr // //////////////////////////////////////////////////////////////////////// #include "G4ios.hh" #include "G4ExceptionHandler.hh" #include #include "G4OpMie.hh" #include "G4SystemOfUnits.hh" #include using namespace CLHEP; ///////////////////////// // Class Implementation ///////////////////////// ////////////// // Operators ////////////// // G4OpMie::operator=(const G4OpMie &right) // { // } ///////////////// // Constructors ///////////////// G4OpMie::G4OpMie(const G4String &processName, G4ProcessType type) : G4VDiscreteProcess(processName, type) { if (verboseLevel > 0) { G4cout << GetProcessName() << " is created " << G4endl; } BuildThePhysicsTable(); } // G4OpMie::G4OpMie(const G4OpMie &right) // { // } //////////////// // Destructors //////////////// G4OpMie::~G4OpMie() { delete PhaseRand; } //////////// // Methods //////////// // PostStepDoIt // ------------- // G4VParticleChange *G4OpMie::PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) { aParticleChange.Initialize(aTrack); // newmie if it has been emmitted as scattered by parametrization kill it // newmie const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle(); if (verboseLevel > 0) { G4cout << "Scattering Photon!" << G4endl; } // find polar angle of new direction w.r.t. old direction G4double CosTheta = std::cos(SampleAngle()); // G4double CosTheta=std::cos(PhaseRand->shoot()*pi); G4double SinTheta = std::sqrt(1. - CosTheta * CosTheta); // find azimuthal angle of new direction w.r.t. old direction G4double rand = G4UniformRand(); G4double Phi = twopi * rand; G4double SinPhi = std::sin(Phi); G4double CosPhi = std::cos(Phi); G4double unit_x = SinTheta * CosPhi; G4double unit_y = SinTheta * SinPhi; G4double unit_z = CosTheta; G4ThreeVector NewMomentumDirection(unit_x, unit_y, unit_z); G4ThreeVector OldMomentumDirection = aParticle->GetMomentumDirection(); NewMomentumDirection.rotateUz(OldMomentumDirection); // find new polarization G4ThreeVector OldPolarization = aParticle->GetPolarization(); G4ThreeVector NewPolarization = OldPolarization - OldPolarization.dot(NewMomentumDirection) * NewMomentumDirection; if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization; NewPolarization = NewPolarization.unit(); aParticleChange.ProposePolarization(NewPolarization); aParticleChange.ProposeMomentumDirection(NewMomentumDirection); if (verboseLevel > 0) { printf("OldNew All %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le " "%13.6le %13.6le %13.6le %13.6le %13.6le\n", OldMomentumDirection.x(), OldMomentumDirection.y(), OldMomentumDirection.z(), NewMomentumDirection.x(), NewMomentumDirection.y(), NewMomentumDirection.z(), OldPolarization.x(), OldPolarization.y(), OldPolarization.z(), NewPolarization.x(), NewPolarization.y(), NewPolarization.z()); } return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); } // BuildThePhysicsTable for the Mie Scattering process (reads only the phase // function model) // -------------------------------------------------------- // void G4OpMie::BuildThePhysicsTable() { // loop for materials const G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable(); G4int numOfMaterials = G4Material::GetNumberOfMaterials(); for (G4int i = 0; i < numOfMaterials; i++) { if ((*theMaterialTable)[i]->GetName() == "Water") { G4MaterialPropertiesTable *aMaterialPropertiesTable = (*theMaterialTable)[i]->GetMaterialPropertiesTable(); if (aMaterialPropertiesTable) { DParticleContribution = aMaterialPropertiesTable->GetConstProperty("MIEETA"); } } } // next build the array, step 0.1 degrees for (G4int i = 0; i <= 1800; i++) { G4double angle = G4double(i) / 10.0; PhaseArray[i] = PhaseFunction(angle * degree); } PhaseRand = new CLHEP::RandGeneral(PhaseArray, 1801); // Forward the seed to the used generator for the dist. sampling gen.seed(G4Random::getTheSeed()); // Sample distribution size_t n = 1000; std::vector angles(n); std::vector prob(n); for(size_t i=0;i(angles.begin(), angles.end(), prob.begin()); } G4double G4OpMie::SampleAngle(void) { G4double angle = (*phase_space_dist)(gen); return angle; } // phase function*sin(theta) with maximum 1 G4double G4OpMie::PhaseFunction(G4double angle) { G4double cosa = std::cos(angle); /* * Water scaattering model according to: * - Einstein-Smoluchowski model pure water scattering * - Henyey-Greenstein parametrisation for particle scattering * c1 - particle relative DParticleContribution * c3 - parameter in Einstein-Smoluchowski model b=(1-delta)/(1+delta) where delta is the H20 anisotropy (delta = 0.09) * c2 - normalisation parameter Integral(pure water scat over 4Pi) = 1 * c4 - HG parameter to fit Petzold data (taken correspondingly to the average cosine of the Petzolds’ angular scattering function) * More info here: https://simulation.pages.km3net.de/input_tables/Simulations_Description.pdf */ static const G4double c3 = 0.835; static const G4double c4 = 0.924; static const G4double c2 = 1.0 / (1.0 + c3/3.0) / (4*std::acos(-1.)); static const G4double c1 = DParticleContribution; G4double val = std::sin(angle) * (c1 * c2 * (1 + c3 * cosa * cosa) + (1 - c1) * (1 / (4 * pi)) * (1 - c4 * c4) / std::pow(1 + c4 * c4 - 2 * c4 * cosa, 1.5)); return val; } // GetMeanFreePath() // ----------------- // G4double G4OpMie::GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *) { const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle(); const G4Material *aMaterial = aTrack.GetMaterial(); G4double thePhotonMomentum = aParticle->GetTotalMomentum(); G4double AttenuationLength = DBL_MAX; G4MaterialPropertiesTable *aMaterialPropertyTable = aMaterial->GetMaterialPropertiesTable(); if (aMaterialPropertyTable) { G4MaterialPropertyVector *AttenuationLengthVector = aMaterialPropertyTable->GetProperty("MIELENGTH"); if (AttenuationLengthVector) { AttenuationLength = AttenuationLengthVector->Value(thePhotonMomentum); } else { // G4cout << "No Mie scattering length specified" << G4endl; } } else { // G4cout << "No Mie scattering length specified" << G4endl; } return AttenuationLength; }