G4OpRayleigh.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 //
00027 // $Id$
00028 //
00029 // 
00031 // Optical Photon Rayleigh Scattering Class Implementation
00033 //
00034 // File:        G4OpRayleigh.cc 
00035 // Description: Discrete Process -- Rayleigh scattering of optical 
00036 //              photons  
00037 // Version:     1.0
00038 // Created:     1996-05-31  
00039 // Author:      Juliet Armstrong
00040 // Updated:     2010-06-11 - Fix Bug 207; Thanks to Xin Qian
00041 //              (Kellogg Radiation Lab of Caltech)
00042 //              2005-07-28 - add G4ProcessType to constructor
00043 //              2001-10-18 by Peter Gumplinger
00044 //              eliminate unused variable warning on Linux (gcc-2.95.2)
00045 //              2001-09-18 by mma
00046 //              >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
00047 //              2001-01-30 by Peter Gumplinger
00048 //              > allow for positiv and negative CosTheta and force the
00049 //              > new momentum direction to be in the same plane as the
00050 //              > new and old polarization vectors
00051 //              2001-01-29 by Peter Gumplinger
00052 //              > fix calculation of SinTheta (from CosTheta)
00053 //              1997-04-09 by Peter Gumplinger
00054 //              > new physics/tracking scheme
00055 // mail:        gum@triumf.ca
00056 //
00058 
00059 #include "G4OpRayleigh.hh"
00060 
00061 #include "G4ios.hh"
00062 #include "G4PhysicalConstants.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "G4OpProcessSubType.hh"
00065 
00067 // Class Implementation
00069 
00071         // Operators
00073 
00074 // G4OpRayleigh::operator=(const G4OpRayleigh &right)
00075 // {
00076 // }
00077 
00079         // Constructors
00081 
00082 G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type)
00083            : G4VDiscreteProcess(processName, type)
00084 {
00085         SetProcessSubType(fOpRayleigh);
00086 
00087         thePhysicsTable = 0;
00088 
00089         DefaultWater = false;
00090 
00091         if (verboseLevel>0) {
00092            G4cout << GetProcessName() << " is created " << G4endl;
00093         }
00094 
00095         BuildThePhysicsTable();
00096 }
00097 
00098 // G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
00099 // {
00100 // }
00101 
00103         // Destructors
00105 
00106 G4OpRayleigh::~G4OpRayleigh()
00107 {
00108         if (thePhysicsTable!= 0) {
00109            thePhysicsTable->clearAndDestroy();
00110            delete thePhysicsTable;
00111         }
00112 }
00113 
00115         // Methods
00117 
00118 // PostStepDoIt
00119 // -------------
00120 //
00121 G4VParticleChange*
00122 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00123 {
00124         aParticleChange.Initialize(aTrack);
00125 
00126         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00127 
00128         if (verboseLevel>0) {
00129                 G4cout << "Scattering Photon!" << G4endl;
00130                 G4cout << "Old Momentum Direction: "
00131                        << aParticle->GetMomentumDirection() << G4endl;
00132                 G4cout << "Old Polarization: "
00133                        << aParticle->GetPolarization() << G4endl;
00134         }
00135 
00136         G4double cosTheta;
00137         G4ThreeVector OldMomentumDirection, NewMomentumDirection;
00138         G4ThreeVector OldPolarization, NewPolarization;
00139 
00140         do {
00141            // Try to simulate the scattered photon momentum direction
00142            // w.r.t. the initial photon momentum direction
00143 
00144            G4double CosTheta = G4UniformRand();
00145            G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
00146            // consider for the angle 90-180 degrees
00147            if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
00148 
00149            // simulate the phi angle
00150            G4double rand = twopi*G4UniformRand();
00151            G4double SinPhi = std::sin(rand);
00152            G4double CosPhi = std::cos(rand);
00153 
00154            // start constructing the new momentum direction
00155            G4double unit_x = SinTheta * CosPhi; 
00156            G4double unit_y = SinTheta * SinPhi;  
00157            G4double unit_z = CosTheta; 
00158            NewMomentumDirection.set (unit_x,unit_y,unit_z);
00159 
00160            // Rotate the new momentum direction into global reference system
00161            OldMomentumDirection = aParticle->GetMomentumDirection();
00162            OldMomentumDirection = OldMomentumDirection.unit();
00163            NewMomentumDirection.rotateUz(OldMomentumDirection);
00164            NewMomentumDirection = NewMomentumDirection.unit();
00165 
00166            // calculate the new polarization direction
00167            // The new polarization needs to be in the same plane as the new
00168            // momentum direction and the old polarization direction
00169            OldPolarization = aParticle->GetPolarization();
00170            G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
00171 
00172            NewPolarization = NewMomentumDirection + constant*OldPolarization;
00173            NewPolarization = NewPolarization.unit();
00174 
00175            // There is a corner case, where the Newmomentum direction
00176            // is the same as oldpolariztion direction:
00177            // random generate the azimuthal angle w.r.t. Newmomentum direction
00178            if (NewPolarization.mag() == 0.) {
00179               rand = G4UniformRand()*twopi;
00180               NewPolarization.set(std::cos(rand),std::sin(rand),0.);
00181               NewPolarization.rotateUz(NewMomentumDirection);
00182            } else {
00183               // There are two directions which are perpendicular
00184               // to the new momentum direction
00185               if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
00186            }
00187           
00188            // simulate according to the distribution cos^2(theta)
00189            cosTheta = NewPolarization.dot(OldPolarization);
00190         } while (std::pow(cosTheta,2) < G4UniformRand());
00191 
00192         aParticleChange.ProposePolarization(NewPolarization);
00193         aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
00194 
00195         if (verboseLevel>0) {
00196                 G4cout << "New Polarization: " 
00197                      << NewPolarization << G4endl;
00198                 G4cout << "Polarization Change: "
00199                      << *(aParticleChange.GetPolarization()) << G4endl;  
00200                 G4cout << "New Momentum Direction: " 
00201                      << NewMomentumDirection << G4endl;
00202                 G4cout << "Momentum Change: "
00203                      << *(aParticleChange.GetMomentumDirection()) << G4endl; 
00204         }
00205 
00206         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00207 }
00208 
00209 // BuildThePhysicsTable for the Rayleigh Scattering process
00210 // --------------------------------------------------------
00211 //
00212 void G4OpRayleigh::BuildThePhysicsTable()
00213 {
00214 //      Builds a table of scattering lengths for each material
00215 
00216         if (thePhysicsTable) return;
00217 
00218         const G4MaterialTable* theMaterialTable=
00219                                G4Material::GetMaterialTable();
00220         G4int numOfMaterials = G4Material::GetNumberOfMaterials();
00221 
00222         // create a new physics table
00223 
00224         thePhysicsTable = new G4PhysicsTable(numOfMaterials);
00225 
00226         // loop for materials
00227 
00228         for (G4int i=0 ; i < numOfMaterials; i++)
00229         {
00230             G4PhysicsOrderedFreeVector* ScatteringLengths = NULL;
00231 
00232             G4MaterialPropertiesTable *aMaterialPropertiesTable =
00233                          (*theMaterialTable)[i]->GetMaterialPropertiesTable();
00234                                                                                 
00235             if(aMaterialPropertiesTable){
00236 
00237               G4MaterialPropertyVector* AttenuationLengthVector =
00238                             aMaterialPropertiesTable->GetProperty("RAYLEIGH");
00239 
00240               if(!AttenuationLengthVector){
00241 
00242                 if ((*theMaterialTable)[i]->GetName() == "Water")
00243                 {
00244                    // Call utility routine to Generate
00245                    // Rayleigh Scattering Lengths
00246 
00247                    DefaultWater = true;
00248 
00249                    ScatteringLengths =
00250                    RayleighAttenuationLengthGenerator(aMaterialPropertiesTable);
00251                 }
00252               }
00253             }
00254 
00255             thePhysicsTable->insertAt(i,ScatteringLengths);
00256         } 
00257 }
00258 
00259 // GetMeanFreePath()
00260 // -----------------
00261 //
00262 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
00263                                      G4double ,
00264                                      G4ForceCondition* )
00265 {
00266         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00267         const G4Material* aMaterial = aTrack.GetMaterial();
00268 
00269         G4double thePhotonEnergy = aParticle->GetTotalEnergy();
00270 
00271         G4double AttenuationLength = DBL_MAX;
00272 
00273         if (aMaterial->GetName() == "Water" && DefaultWater){
00274 
00275            G4bool isOutRange;
00276 
00277            AttenuationLength =
00278                 (*thePhysicsTable)(aMaterial->GetIndex())->
00279                            GetValue(thePhotonEnergy, isOutRange);
00280         }
00281         else {
00282 
00283            G4MaterialPropertiesTable* aMaterialPropertyTable =
00284                            aMaterial->GetMaterialPropertiesTable();
00285 
00286            if(aMaterialPropertyTable){
00287              G4MaterialPropertyVector* AttenuationLengthVector =
00288                    aMaterialPropertyTable->GetProperty("RAYLEIGH");
00289              if(AttenuationLengthVector){
00290                AttenuationLength = AttenuationLengthVector ->
00291                                     Value(thePhotonEnergy);
00292              }
00293              else{
00294 //               G4cout << "No Rayleigh scattering length specified" << G4endl;
00295              }
00296            }
00297            else{
00298 //             G4cout << "No Rayleigh scattering length specified" << G4endl; 
00299            }
00300         }
00301 
00302         return AttenuationLength;
00303 }
00304 
00305 // RayleighAttenuationLengthGenerator()
00306 // ------------------------------------
00307 // Private method to compute Rayleigh Scattering Lengths (for water)
00308 //
00309 G4PhysicsOrderedFreeVector* 
00310 G4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT) 
00311 {
00312         // Physical Constants
00313 
00314         // isothermal compressibility of water
00315         G4double betat = 7.658e-23*m3/MeV;
00316 
00317         // K Boltzman
00318         G4double kboltz = 8.61739e-11*MeV/kelvin;
00319 
00320         // Temperature of water is 10 degrees celsius
00321         // conversion to kelvin:
00322         // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15
00323         G4double temp = 283.15*kelvin;
00324 
00325         // Retrieve vectors for refraction index
00326         // and photon energy from the material properties table
00327 
00328         G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");
00329 
00330         G4double refsq;
00331         G4double e;
00332         G4double xlambda;
00333         G4double c1, c2, c3, c4;
00334         G4double Dist;
00335         G4double refraction_index;
00336 
00337         G4PhysicsOrderedFreeVector *RayleighScatteringLengths = 
00338                                 new G4PhysicsOrderedFreeVector();
00339 
00340         if (Rindex ) {
00341 
00342            for (size_t i = 0; i < Rindex->GetVectorLength(); i++) {
00343 
00344                 e = Rindex->Energy(i);
00345 
00346                 refraction_index = (*Rindex)[i];
00347 
00348                 refsq = refraction_index*refraction_index;
00349                 xlambda = h_Planck*c_light/e;
00350 
00351                 if (verboseLevel>0) {
00352                         G4cout << Rindex->Energy(i) << " MeV\t";
00353                         G4cout << xlambda << " mm\t";
00354                 }
00355 
00356                 c1 = 1 / (6.0 * pi);
00357                 c2 = std::pow((2.0 * pi / xlambda), 4);
00358                 c3 = std::pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
00359                 c4 = betat * temp * kboltz;
00360 
00361                 Dist = 1.0 / (c1*c2*c3*c4);
00362 
00363                 if (verboseLevel>0) {
00364                         G4cout << Dist << " mm" << G4endl;
00365                 }
00366                 RayleighScatteringLengths->
00367                         InsertValues(Rindex->Energy(i), Dist);
00368            }
00369 
00370         }
00371 
00372         return RayleighScatteringLengths;
00373 }

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