/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

DsG4OpRayleigh.cc
Go to the documentation of this file.
00001 //
00002 //
00014 // ********************************************************************
00015 // * License and Disclaimer                                           *
00016 // *                                                                  *
00017 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00018 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00019 // * conditions of the Geant4 Software License,  included in the file *
00020 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00021 // * include a list of copyright holders.                             *
00022 // *                                                                  *
00023 // * Neither the authors of this software system, nor their employing *
00024 // * institutes,nor the agencies providing financial support for this *
00025 // * work  make  any representation or  warranty, express or implied, *
00026 // * regarding  this  software system or assume any liability for its *
00027 // * use.  Please see the license in the file  LICENSE  and URL above *
00028 // * for the full disclaimer and the limitation of liability.         *
00029 // *                                                                  *
00030 // * This  code  implementation is the result of  the  scientific and *
00031 // * technical work of the GEANT4 collaboration.                      *
00032 // * By using,  copying,  modifying or  distributing the software (or *
00033 // * any work based  on the software)  you  agree  to acknowledge its *
00034 // * use  in  resulting  scientific  publications,  and indicate your *
00035 // * acceptance of all terms of the Geant4 Software license.          *
00036 // ********************************************************************
00037 //
00038 //
00039 // $Id: G4OpRayleigh.cc,v 1.17 2008/10/24 19:51:12 gum Exp $
00040 // GEANT4 tag $Name: geant4-09-02 $
00041 //
00042 // 
00044 // Optical Photon Rayleigh Scattering Class Implementation
00046 //
00047 // File:        DsG4OpRayleigh.cc 
00048 // Description: Discrete Process -- Rayleigh scattering of optical 
00049 //              photons  
00050 // Version:     1.0
00051 // Created:     1996-05-31  
00052 // Author:      Juliet Armstrong
00053 // Updated:     2005-07-28 - add G4ProcessType to constructor
00054 //              2001-10-18 by Peter Gumplinger
00055 //              eliminate unused variable warning on Linux (gcc-2.95.2)
00056 //              2001-09-18 by mma
00057 //              >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
00058 //              2001-01-30 by Peter Gumplinger
00059 //              > allow for positiv and negative CosTheta and force the
00060 //              > new momentum direction to be in the same plane as the
00061 //              > new and old polarization vectors
00062 //              2001-01-29 by Peter Gumplinger
00063 //              > fix calculation of SinTheta (from CosTheta)
00064 //              1997-04-09 by Peter Gumplinger
00065 //              > new physics/tracking scheme
00066 // mail:        gum@triumf.ca
00067 //
00069 
00070 #include "G4ios.hh"
00071 #include "G4OpProcessSubType.hh"
00072 
00073 #include "DsG4OpRayleigh.h"
00074 
00076 // Class Implementation
00078 
00080         // Operators
00082 
00083 // DsG4OpRayleigh::operator=(const DsG4OpRayleigh &right)
00084 // {
00085 // }
00086 
00088         // Constructors
00090 
00091 DsG4OpRayleigh::DsG4OpRayleigh(const G4String& processName, G4ProcessType type)
00092            : G4VDiscreteProcess(processName, type)
00093 {
00094         SetProcessSubType(fOpRayleigh);
00095 
00096         thePhysicsTable = 0;
00097 
00098         DefaultWater = false;
00099 
00100         if (verboseLevel>0) {
00101            G4cout << GetProcessName() << " is created " << G4endl;
00102         }
00103 
00104         BuildThePhysicsTable();
00105 }
00106 
00107 // DsG4OpRayleigh::DsG4OpRayleigh(const DsG4OpRayleigh &right)
00108 // {
00109 // }
00110 
00112         // Destructors
00114 
00115 DsG4OpRayleigh::~DsG4OpRayleigh()
00116 {
00117         if (thePhysicsTable!= 0) {
00118            thePhysicsTable->clearAndDestroy();
00119            delete thePhysicsTable;
00120         }
00121 }
00122 
00124         // Methods
00126 
00127 // PostStepDoIt
00128 // -------------
00129 //
00130 G4VParticleChange* 
00131 DsG4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00132 {
00133         aParticleChange.Initialize(aTrack);
00134 
00135         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00136 
00137         if (verboseLevel>0) {
00138                 G4cout << "Scattering Photon!" << G4endl;
00139                 G4cout << "Old Momentum Direction: "
00140                      << aParticle->GetMomentumDirection() << G4endl;
00141                 G4cout << "Old Polarization: "
00142                      << aParticle->GetPolarization() << G4endl;
00143         }
00144 
00145         // find polar angle w.r.t. old polarization vector
00146 
00147         G4double rand = G4UniformRand();
00148 
00149         G4double CosTheta = std::pow(rand, 1./3.);
00150         G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
00151 
00152         if(G4UniformRand() < 0.5)CosTheta = -CosTheta;
00153 
00154         // find azimuthal angle w.r.t old polarization vector 
00155 
00156         rand = G4UniformRand();
00157 
00158         G4double Phi = twopi*rand;
00159         G4double SinPhi = std::sin(Phi); 
00160         G4double CosPhi = std::cos(Phi); 
00161         
00162         G4double unit_x = SinTheta * CosPhi; 
00163         G4double unit_y = SinTheta * SinPhi;  
00164         G4double unit_z = CosTheta; 
00165         
00166         G4ThreeVector NewPolarization (unit_x,unit_y,unit_z);
00167 
00168         // Rotate new polarization direction into global reference system 
00169 
00170         G4ThreeVector OldPolarization = aParticle->GetPolarization();
00171         OldPolarization = OldPolarization.unit();
00172 
00173         NewPolarization.rotateUz(OldPolarization);
00174         NewPolarization = NewPolarization.unit();
00175         
00176         // -- new momentum direction is normal to the new
00177         // polarization vector and in the same plane as the
00178         // old and new polarization vectors --
00179 
00180         G4ThreeVector NewMomentumDirection = 
00181                               OldPolarization - NewPolarization * CosTheta;
00182 
00183         if(G4UniformRand() < 0.5)NewMomentumDirection = -NewMomentumDirection;
00184         NewMomentumDirection = NewMomentumDirection.unit();
00185 
00186         aParticleChange.ProposePolarization(NewPolarization);
00187 
00188         aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
00189 
00190         if (verboseLevel>0) {
00191                 G4cout << "New Polarization: " 
00192                      << NewPolarization << G4endl;
00193                 G4cout << "Polarization Change: "
00194                      << *(aParticleChange.GetPolarization()) << G4endl;  
00195                 G4cout << "New Momentum Direction: " 
00196                      << NewMomentumDirection << G4endl;
00197                 G4cout << "Momentum Change: "
00198                      << *(aParticleChange.GetMomentumDirection()) << G4endl; 
00199         }
00200 
00201         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00202 }
00203 
00204 // BuildThePhysicsTable for the Rayleigh Scattering process
00205 // --------------------------------------------------------
00206 //
00207 void DsG4OpRayleigh::BuildThePhysicsTable()
00208 {
00209 //      Builds a table of scattering lengths for each material
00210 
00211         if (thePhysicsTable) return;
00212 
00213         const G4MaterialTable* theMaterialTable=
00214                                G4Material::GetMaterialTable();
00215         G4int numOfMaterials = G4Material::GetNumberOfMaterials();
00216 
00217         // create a new physics table
00218 
00219         thePhysicsTable = new G4PhysicsTable(numOfMaterials);
00220 
00221         // loop for materials
00222 
00223         for (G4int i=0 ; i < numOfMaterials; i++)
00224         {
00225             G4PhysicsOrderedFreeVector* ScatteringLengths =
00226                                 new G4PhysicsOrderedFreeVector();
00227 
00228             G4MaterialPropertiesTable *aMaterialPropertiesTable =
00229                          (*theMaterialTable)[i]->GetMaterialPropertiesTable();
00230                                                                                 
00231             if(aMaterialPropertiesTable){
00232 
00233               G4MaterialPropertyVector* AttenuationLengthVector =
00234                             aMaterialPropertiesTable->GetProperty("RAYLEIGH");
00235 
00236               if(!AttenuationLengthVector){
00237 
00238                 if ((*theMaterialTable)[i]->GetName() == "/dd/Materials/Water" || 
00239                     (*theMaterialTable)[i]->GetName() == "/dd/Materials/OwsWater" ||
00240                     (*theMaterialTable)[i]->GetName() == "/dd/Materials/IwsWater"
00241                     )
00242                 {
00243                    // Call utility routine to Generate
00244                    // Rayleigh Scattering Lengths
00245 
00246                    DefaultWater = true;
00247 
00248                    ScatteringLengths =
00249                    RayleighAttenuationLengthGenerator(aMaterialPropertiesTable);
00250                 }
00251               }
00252             }
00253 
00254             thePhysicsTable->insertAt(i,ScatteringLengths);
00255         } 
00256 }
00257 
00258 // GetMeanFreePath()
00259 // -----------------
00260 //
00261 G4double DsG4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
00262                                      G4double ,
00263                                      G4ForceCondition* )
00264 {
00265         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00266         const G4Material* aMaterial = aTrack.GetMaterial();
00267 
00268         G4double thePhotonEnergy = aParticle->GetTotalEnergy();
00269 
00270         G4double AttenuationLength = DBL_MAX;
00271 
00272         if ((strcmp(aMaterial->GetName(), "/dd/Materials/Water") == 0 ||
00273              strcmp(aMaterial->GetName(), "/dd/Materials/OwsWater") == 0 ||
00274              strcmp(aMaterial->GetName(), "/dd/Materials/IwsWater") == 0 )
00275             && DefaultWater){
00276 
00277            G4bool isOutRange;
00278 
00279            AttenuationLength =
00280                 (*thePhysicsTable)(aMaterial->GetIndex())->
00281                            GetValue(thePhotonEnergy, isOutRange);
00282         }
00283         else {
00284 
00285            G4MaterialPropertiesTable* aMaterialPropertyTable =
00286                            aMaterial->GetMaterialPropertiesTable();
00287 
00288            if(aMaterialPropertyTable){
00289              G4MaterialPropertyVector* AttenuationLengthVector =
00290                    aMaterialPropertyTable->GetProperty("RAYLEIGH");
00291              if(AttenuationLengthVector){
00292                AttenuationLength = AttenuationLengthVector ->
00293                                     GetProperty(thePhotonEnergy);
00294              }
00295              else{
00296 //               G4cout << "No Rayleigh scattering length specified" << G4endl;
00297              }
00298            }
00299            else{
00300 //             G4cout << "No Rayleigh scattering length specified" << G4endl; 
00301            }
00302         }
00303 
00304         return AttenuationLength;
00305 }
00306 
00307 // RayleighAttenuationLengthGenerator()
00308 // ------------------------------------
00309 // Private method to compute Rayleigh Scattering Lengths (for water)
00310 //
00311 G4PhysicsOrderedFreeVector* 
00312 DsG4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT) 
00313 {
00314         // Physical Constants
00315 
00316         // isothermal compressibility of water
00317         G4double betat = 7.658e-23*m3/MeV;
00318 
00319         // K Boltzman
00320         G4double kboltz = 8.61739e-11*MeV/kelvin;
00321 
00322         // Temperature of water is 10 degrees celsius
00323         // conversion to kelvin:
00324         // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15
00325         G4double temp = 283.15*kelvin;
00326 
00327         // Retrieve vectors for refraction index
00328         // and photon energy from the material properties table
00329 
00330         G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");
00331 
00332         G4double refsq;
00333         G4double e;
00334         G4double xlambda;
00335         G4double c1, c2, c3, c4;
00336         G4double Dist;
00337         G4double refraction_index;
00338 
00339         G4PhysicsOrderedFreeVector *RayleighScatteringLengths = 
00340                                 new G4PhysicsOrderedFreeVector();
00341 
00342         if (Rindex ) {
00343 
00344            Rindex->ResetIterator();
00345 
00346            while (++(*Rindex)) {
00347 
00348                 e = (Rindex->GetPhotonEnergy());
00349 
00350                 refraction_index = Rindex->GetProperty();
00351                 refsq = refraction_index*refraction_index;
00352                 xlambda = h_Planck*c_light/e;
00353 
00354                 if (verboseLevel>0) {
00355                         G4cout << Rindex->GetPhotonEnergy() << " MeV\t";
00356                         G4cout << xlambda << " mm\t";
00357                 }
00358 
00359                 c1 = 1 / (6.0 * pi);
00360                 c2 = std::pow((2.0 * pi / xlambda), 4);
00361                 c3 = std::pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
00362                 c4 = betat * temp * kboltz;
00363 
00364                 Dist = 1.0 / (c1*c2*c3*c4);
00365 
00366                 if (verboseLevel>0) {
00367                         G4cout << Dist << " mm" << G4endl;
00368                 }
00369                 RayleighScatteringLengths->
00370                         InsertValues(Rindex->GetPhotonEnergy(), Dist);
00371            }
00372 
00373         }
00374 
00375         return RayleighScatteringLengths;
00376 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:17:58 for DetSim by doxygen 1.7.4