G4OpMieHG.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00028 //
00029 // File G4OpMieHG.hh
00030 // Description: Discrete Process -- Mie Scattering of Optical Photons
00031 // Created: 2010-07-03
00032 // Author: Xin Qian
00033 // Based on work from Vlasios Vasileiou
00034 //
00035 // This subroutine will mimic the Mie scattering based on 
00036 // Henyey-Greenstein phase function
00037 // Forward and backward angles are treated separately.
00038 //
00039 // mail:  gum@triumf.ca
00040 //
00042 
00043 #include "G4OpMieHG.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4OpProcessSubType.hh"
00046 
00047 G4OpMieHG::G4OpMieHG(const G4String& processName, G4ProcessType type)
00048            : G4VDiscreteProcess(processName, type)
00049 {
00050         if (verboseLevel>0) {
00051            G4cout << GetProcessName() << " is created " << G4endl;
00052         }
00053 
00054         SetProcessSubType(fOpMieHG);
00055 }
00056 
00057 G4OpMieHG::~G4OpMieHG(){}
00058 
00060         // Methods
00062 
00063 // PostStepDoIt
00064 // -------------
00065 //
00066 G4VParticleChange* 
00067 G4OpMieHG::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00068 {
00069         aParticleChange.Initialize(aTrack);
00070 
00071         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00072         const G4Material* aMaterial = aTrack.GetMaterial();
00073         G4MaterialPropertiesTable* aMaterialPropertyTable =
00074           aMaterial->GetMaterialPropertiesTable();
00075 
00076         G4double forward_g =
00077               aMaterialPropertyTable->GetConstProperty("MIEHG_FORWARD");
00078         G4double backward_g =
00079               aMaterialPropertyTable->GetConstProperty("MIEHG_BACKWARD");
00080         G4double ForwardRatio =
00081               aMaterialPropertyTable->GetConstProperty("MIEHG_FORWARD_RATIO");
00082 
00083         if (verboseLevel>0) {
00084                 G4cout << "MIE Scattering Photon!" << G4endl;
00085                 G4cout << "MIE Old Momentum Direction: "
00086                      << aParticle->GetMomentumDirection() << G4endl;
00087                 G4cout << "MIE Old Polarization: "
00088                      << aParticle->GetPolarization() << G4endl;
00089         }
00090 
00091         G4double gg;
00092         G4int direction;
00093         if (G4UniformRand()<=ForwardRatio){
00094            gg = forward_g;
00095            direction = 1;
00096         } else {
00097            gg = backward_g;
00098            direction = -1;
00099         }
00100 
00101         G4double r = G4UniformRand();
00102 
00103         G4double Theta;
00104         //sample the direction
00105         if (gg!=0) {
00106           Theta = std::acos(2*r*(1+gg)*(1+gg)*(1-gg+gg*r)/((1-gg+2*gg*r)*(1-gg+2*gg*r)) -1);
00107         } else {
00108           Theta = std::acos(2*r-1.);
00109         }
00110         G4double Phi = G4UniformRand()*2*pi;
00111 
00112         if (direction==-1) Theta = pi - Theta; //backward scattering
00113 
00114         G4ThreeVector NewMomentumDirection, OldMomentumDirection;
00115         G4ThreeVector OldPolarization, NewPolarization;
00116 
00117         NewMomentumDirection.set
00118                        (std::sin(Theta)*std::cos(Phi), std::sin(Theta)*std::sin(Phi), std::cos(Theta));
00119         OldMomentumDirection = aParticle->GetMomentumDirection();
00120         NewMomentumDirection.rotateUz(OldMomentumDirection);
00121         NewMomentumDirection = NewMomentumDirection.unit();
00122 
00123         OldPolarization = aParticle->GetPolarization();
00124         G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
00125 
00126         NewPolarization = NewMomentumDirection + constant*OldPolarization;
00127         NewPolarization = NewPolarization.unit();
00128 
00129         if (NewPolarization.mag()==0) {
00130            r = G4UniformRand()*twopi;
00131            NewPolarization.set(std::cos(r),std::sin(r),0.);
00132            NewPolarization.rotateUz(NewMomentumDirection);
00133         } else {
00134            // There are two directions which perpendicular
00135            // new momentum direction
00136            if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
00137         }
00138 
00139         aParticleChange.ProposePolarization(NewPolarization);
00140         aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
00141 
00142         if (verboseLevel>0) {
00143               G4cout << "MIE New Polarization: " 
00144                      << NewPolarization << G4endl;
00145               G4cout << "MIE Polarization Change: "
00146                      << *(aParticleChange.GetPolarization()) << G4endl;  
00147               G4cout << "MIE New Momentum Direction: " 
00148                      << NewMomentumDirection << G4endl;
00149               G4cout << "MIE Momentum Change: "
00150                      << *(aParticleChange.GetMomentumDirection()) << G4endl; 
00151         }
00152 
00153         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00154 }
00155 
00156 // GetMeanFreePath()
00157 // -----------------
00158 //
00159 G4double G4OpMieHG::GetMeanFreePath(const G4Track& aTrack,
00160                                     G4double ,
00161                                     G4ForceCondition* )
00162 {
00163         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00164         const G4Material* aMaterial = aTrack.GetMaterial();
00165 
00166         G4double thePhotonEnergy = aParticle->GetTotalEnergy();
00167 
00168         G4double AttenuationLength = DBL_MAX;
00169 
00170         G4MaterialPropertiesTable* aMaterialPropertyTable =
00171           aMaterial->GetMaterialPropertiesTable();
00172 
00173         if (aMaterialPropertyTable) {
00174            G4MaterialPropertyVector* AttenuationLengthVector =
00175                                  aMaterialPropertyTable->GetProperty("MIEHG");
00176            if (AttenuationLengthVector) {
00177               AttenuationLength = AttenuationLengthVector ->
00178                                     Value(thePhotonEnergy);
00179            } else {
00180 //              G4cout << "No Mie scattering length specified" << G4endl;
00181            }
00182         } else {
00183 //           G4cout << "No Mie scattering length specified" << G4endl; 
00184         }
00185 
00186 //      G4cout << thePhotonEnergy/GeV << " \t" << AttenuationLength/m << G4endl;
00187 
00188         return AttenuationLength;
00189 }

Generated on Mon May 27 17:49:13 2013 for Geant4 by  doxygen 1.4.7