Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/optical/src/G4OpRayleigh.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/optical/src/G4OpRayleigh.cc (Version 11.3.0) and /processes/optical/src/G4OpRayleigh.cc (Version 10.5)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 //                                                 27 //
 28 //                                             <<  28 // 
 29 //////////////////////////////////////////////     29 ////////////////////////////////////////////////////////////////////////
 30 // Optical Photon Rayleigh Scattering Class Im     30 // Optical Photon Rayleigh Scattering Class Implementation
 31 //////////////////////////////////////////////     31 ////////////////////////////////////////////////////////////////////////
 32 //                                                 32 //
 33 // File:        G4OpRayleigh.cc                    33 // File:        G4OpRayleigh.cc
 34 // Description: Discrete Process -- Rayleigh s     34 // Description: Discrete Process -- Rayleigh scattering of optical
 35 //    photons                                      35 //    photons
 36 // Version:     1.0                                36 // Version:     1.0
 37 // Created:     1996-05-31                         37 // Created:     1996-05-31
 38 // Author:      Juliet Armstrong                   38 // Author:      Juliet Armstrong
 39 // Updated:     2014-10-10 -  This version cal <<  39 // Updated:     2014-10-10 -  This version calculates the Rayleigh scattering   
 40 //              length for more materials than     40 //              length for more materials than just Water (although the Water
 41 //              default is kept). To do this t     41 //              default is kept). To do this the user would need to specify the
 42 //              ISOTHERMAL_COMPRESSIBILITY as      42 //              ISOTHERMAL_COMPRESSIBILITY as a material property and
 43 //              optionally an RS_SCALE_LENGTH      43 //              optionally an RS_SCALE_LENGTH (useful for testing). Code comes
 44 //              from Philip Graham (Queen Mary     44 //              from Philip Graham (Queen Mary University of London).
 45 //              2010-06-11 - Fix Bug 207; Than     45 //              2010-06-11 - Fix Bug 207; Thanks to Xin Qian
 46 //              (Kellogg Radiation Lab of Calt     46 //              (Kellogg Radiation Lab of Caltech)
 47 //              2005-07-28 - add G4ProcessType     47 //              2005-07-28 - add G4ProcessType to constructor
 48 //              2001-10-18 by Peter Gumplinger     48 //              2001-10-18 by Peter Gumplinger
 49 //              eliminate unused variable warn     49 //              eliminate unused variable warning on Linux (gcc-2.95.2)
 50 //              2001-09-18 by mma                  50 //              2001-09-18 by mma
 51 //              >numOfMaterials=G4Material::Ge <<  51 //    >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
 52 //              2001-01-30 by Peter Gumplinger     52 //              2001-01-30 by Peter Gumplinger
 53 //              > allow for positiv and negati     53 //              > allow for positiv and negative CosTheta and force the
 54 //              > new momentum direction to be     54 //              > new momentum direction to be in the same plane as the
 55 //              > new and old polarization vec     55 //              > new and old polarization vectors
 56 //              2001-01-29 by Peter Gumplinger     56 //              2001-01-29 by Peter Gumplinger
 57 //              > fix calculation of SinTheta      57 //              > fix calculation of SinTheta (from CosTheta)
 58 //              1997-04-09 by Peter Gumplinger     58 //              1997-04-09 by Peter Gumplinger
 59 //              > new physics/tracking scheme      59 //              > new physics/tracking scheme
                                                   >>  60 // mail:        gum@triumf.ca
 60 //                                                 61 //
 61 //////////////////////////////////////////////     62 ////////////////////////////////////////////////////////////////////////
 62                                                    63 
 63 #include "G4OpRayleigh.hh"                         64 #include "G4OpRayleigh.hh"
                                                   >>  65 
 64 #include "G4ios.hh"                                66 #include "G4ios.hh"
 65 #include "G4PhysicalConstants.hh"                  67 #include "G4PhysicalConstants.hh"
 66 #include "G4SystemOfUnits.hh"                      68 #include "G4SystemOfUnits.hh"
 67 #include "G4OpticalParameters.hh"              << 
 68 #include "G4OpProcessSubType.hh"                   69 #include "G4OpProcessSubType.hh"
 69                                                    70 
 70 //....oooOO0OOooo........oooOO0OOooo........oo <<  71 /////////////////////////
                                                   >>  72 // Class Implementation
                                                   >>  73 /////////////////////////
                                                   >>  74 
                                                   >>  75         //////////////
                                                   >>  76         // Operators
                                                   >>  77         //////////////
                                                   >>  78 
                                                   >>  79 // G4OpRayleigh::operator=(const G4OpRayleigh &right)
                                                   >>  80 // {
                                                   >>  81 // }
                                                   >>  82 
                                                   >>  83         /////////////////
                                                   >>  84         // Constructors
                                                   >>  85         /////////////////
                                                   >>  86 
 71 G4OpRayleigh::G4OpRayleigh(const G4String& pro     87 G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type)
 72   : G4VDiscreteProcess(processName, type)      <<  88            : G4VDiscreteProcess(processName, type)
 73 {                                                  89 {
 74   Initialise();                                <<  90         SetProcessSubType(fOpRayleigh);
 75   SetProcessSubType(fOpRayleigh);              << 
 76   thePhysicsTable = nullptr;                   << 
 77                                                    91 
 78   if(verboseLevel > 0)                         <<  92         thePhysicsTable = NULL;
 79   {                                            << 
 80     G4cout << GetProcessName() << " is created << 
 81   }                                            << 
 82 }                                              << 
 83                                                    93 
 84 //....oooOO0OOooo........oooOO0OOooo........oo <<  94         if (verboseLevel>0) {
 85 G4OpRayleigh::~G4OpRayleigh()                  <<  95            G4cout << GetProcessName() << " is created " << G4endl;
 86 {                                              <<  96         }
 87   // VI: inside this PhysicsTable all properti << 
 88   //     it is not possible to destroy         << 
 89   delete thePhysicsTable;                      << 
 90 }                                                  97 }
 91                                                    98 
 92 //....oooOO0OOooo........oooOO0OOooo........oo <<  99 // G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
 93 void G4OpRayleigh::PreparePhysicsTable(const G << 100 // {
 94 {                                              << 101 // }
 95   Initialise();                                << 102 
 96 }                                              << 103         ////////////////
                                                   >> 104         // Destructors
                                                   >> 105         ////////////////
 97                                                   106 
 98 //....oooOO0OOooo........oooOO0OOooo........oo << 107 G4OpRayleigh::~G4OpRayleigh()
 99 void G4OpRayleigh::Initialise()                << 
100 {                                                 108 {
101   SetVerboseLevel(G4OpticalParameters::Instanc << 109         if (thePhysicsTable) {
                                                   >> 110            thePhysicsTable->clearAndDestroy();
                                                   >> 111            delete thePhysicsTable;
                                                   >> 112         }
102 }                                                 113 }
103                                                   114 
104 //....oooOO0OOooo........oooOO0OOooo........oo << 115         ////////////
105 G4VParticleChange* G4OpRayleigh::PostStepDoIt( << 116         // Methods
106                                                << 117         ////////////
107 {                                              << 
108   aParticleChange.Initialize(aTrack);          << 
109   const G4DynamicParticle* aParticle = aTrack. << 
110                                                << 
111   if(verboseLevel > 1)                         << 
112   {                                            << 
113     G4cout << "OpRayleigh: Scattering Photon!" << 
114            << "Old Momentum Direction: " << aP << 
115            << G4endl << "Old Polarization: " < << 
116            << G4endl;                          << 
117   }                                            << 
118                                                << 
119   G4double cosTheta;                           << 
120   G4ThreeVector oldMomDir, newMomDir;          << 
121   G4ThreeVector oldPol, newPol;                << 
122   G4double rand;                               << 
123   G4double cost, sint, sinphi, cosphi;         << 
124                                                   118 
125   do                                           << 119 // PostStepDoIt
126   {                                            << 120 // -------------
127     // Try to simulate the scattered photon mo << 121 //
128     // w.r.t. the initial photon momentum dire << 122 G4VParticleChange*
129     cost = G4UniformRand();                    << 123 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
130     sint = std::sqrt(1. - cost * cost);        << 124 {
131     // consider for the angle 90-180 degrees   << 125         aParticleChange.Initialize(aTrack);
132     if(G4UniformRand() < 0.5)                  << 
133       cost = -cost;                            << 
134                                                << 
135     // simulate the phi angle                  << 
136     rand   = twopi * G4UniformRand();          << 
137     sinphi = std::sin(rand);                   << 
138     cosphi = std::cos(rand);                   << 
139                                                << 
140     // construct the new momentum direction    << 
141     newMomDir.set(sint * cosphi, sint * sinphi << 
142     oldMomDir = aParticle->GetMomentumDirectio << 
143     newMomDir.rotateUz(oldMomDir);             << 
144                                                << 
145     // calculate the new polarization directio << 
146     // The new polarization needs to be in the << 
147     // momentum direction and the old polariza << 
148     oldPol = aParticle->GetPolarization();     << 
149     newPol = (oldPol - newMomDir.dot(oldPol) * << 
150                                                << 
151     // There is a corner case, where the new m << 
152     // is the same as old polarization directi << 
153     // random generate the azimuthal angle w.r << 
154     if(newPol.mag() == 0.)                     << 
155     {                                          << 
156       rand = G4UniformRand() * twopi;          << 
157       newPol.set(std::cos(rand), std::sin(rand << 
158       newPol.rotateUz(newMomDir);              << 
159     }                                          << 
160     else                                       << 
161     {                                          << 
162       // There are two directions perpendicula << 
163       if(G4UniformRand() < 0.5)                << 
164         newPol = -newPol;                      << 
165     }                                          << 
166                                                << 
167     // simulate according to the distribution  << 
168     cosTheta = newPol.dot(oldPol);             << 
169     // Loop checking, 13-Aug-2015, Peter Gumpl << 
170   } while(std::pow(cosTheta, 2) < G4UniformRan << 
171                                                   126 
172   aParticleChange.ProposePolarization(newPol); << 127         const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
173   aParticleChange.ProposeMomentumDirection(new << 
174                                                   128 
175   if(verboseLevel > 1)                         << 129         if (verboseLevel>0) {
176   {                                            << 130                 G4cout << "Scattering Photon!" << G4endl;
177     G4cout << "New Polarization: " << newPol < << 131                 G4cout << "Old Momentum Direction: "
178            << "Polarization Change: " << *(aPa << 132                        << aParticle->GetMomentumDirection() << G4endl;
179            << G4endl << "New Momentum Directio << 133                 G4cout << "Old Polarization: "
180            << "Momentum Change: " << *(aPartic << 134                        << aParticle->GetPolarization() << G4endl;
181            << G4endl;                          << 135         }
182   }                                            << 136 
                                                   >> 137         G4double cosTheta;
                                                   >> 138         G4ThreeVector OldMomentumDirection, NewMomentumDirection;
                                                   >> 139         G4ThreeVector OldPolarization, NewPolarization;
                                                   >> 140 
                                                   >> 141         G4double rand, constant;
                                                   >> 142         G4double CosTheta, SinTheta, SinPhi, CosPhi, unit_x, unit_y, unit_z;
                                                   >> 143 
                                                   >> 144         do {
                                                   >> 145            // Try to simulate the scattered photon momentum direction
                                                   >> 146            // w.r.t. the initial photon momentum direction
                                                   >> 147 
                                                   >> 148            CosTheta = G4UniformRand();
                                                   >> 149            SinTheta = std::sqrt(1.-CosTheta*CosTheta);
                                                   >> 150            // consider for the angle 90-180 degrees
                                                   >> 151            if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
                                                   >> 152 
                                                   >> 153            // simulate the phi angle
                                                   >> 154            rand = twopi*G4UniformRand();
                                                   >> 155            SinPhi = std::sin(rand);
                                                   >> 156            CosPhi = std::cos(rand);
                                                   >> 157 
                                                   >> 158            // start constructing the new momentum direction
                                                   >> 159      unit_x = SinTheta * CosPhi; 
                                                   >> 160      unit_y = SinTheta * SinPhi;  
                                                   >> 161      unit_z = CosTheta; 
                                                   >> 162      NewMomentumDirection.set (unit_x,unit_y,unit_z);
                                                   >> 163 
                                                   >> 164            // Rotate the new momentum direction into global reference system
                                                   >> 165            OldMomentumDirection = aParticle->GetMomentumDirection();
                                                   >> 166            OldMomentumDirection = OldMomentumDirection.unit();
                                                   >> 167            NewMomentumDirection.rotateUz(OldMomentumDirection);
                                                   >> 168            NewMomentumDirection = NewMomentumDirection.unit();
                                                   >> 169 
                                                   >> 170            // calculate the new polarization direction
                                                   >> 171            // The new polarization needs to be in the same plane as the new
                                                   >> 172            // momentum direction and the old polarization direction
                                                   >> 173            OldPolarization = aParticle->GetPolarization();
                                                   >> 174            constant = -NewMomentumDirection.dot(OldPolarization);
                                                   >> 175 
                                                   >> 176            NewPolarization = OldPolarization + constant*NewMomentumDirection;
                                                   >> 177            NewPolarization = NewPolarization.unit();
                                                   >> 178 
                                                   >> 179            // There is a corner case, where the Newmomentum direction
                                                   >> 180            // is the same as oldpolariztion direction:
                                                   >> 181            // random generate the azimuthal angle w.r.t. Newmomentum direction
                                                   >> 182            if (NewPolarization.mag() == 0.) {
                                                   >> 183               rand = G4UniformRand()*twopi;
                                                   >> 184               NewPolarization.set(std::cos(rand),std::sin(rand),0.);
                                                   >> 185               NewPolarization.rotateUz(NewMomentumDirection);
                                                   >> 186            } else {
                                                   >> 187               // There are two directions which are perpendicular
                                                   >> 188               // to the new momentum direction
                                                   >> 189               if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
                                                   >> 190            }
                                                   >> 191     
                                                   >> 192      // simulate according to the distribution cos^2(theta)
                                                   >> 193            cosTheta = NewPolarization.dot(OldPolarization);
                                                   >> 194           // Loop checking, 13-Aug-2015, Peter Gumplinger
                                                   >> 195         } while (std::pow(cosTheta,2) < G4UniformRand());
                                                   >> 196 
                                                   >> 197         aParticleChange.ProposePolarization(NewPolarization);
                                                   >> 198         aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
                                                   >> 199 
                                                   >> 200         if (verboseLevel>0) {
                                                   >> 201                 G4cout << "New Polarization: " 
                                                   >> 202                      << NewPolarization << G4endl;
                                                   >> 203                 G4cout << "Polarization Change: "
                                                   >> 204                      << *(aParticleChange.GetPolarization()) << G4endl;  
                                                   >> 205                 G4cout << "New Momentum Direction: " 
                                                   >> 206                      << NewMomentumDirection << G4endl;
                                                   >> 207                 G4cout << "Momentum Change: "
                                                   >> 208                      << *(aParticleChange.GetMomentumDirection()) << G4endl; 
                                                   >> 209         }
183                                                   210 
184   return G4VDiscreteProcess::PostStepDoIt(aTra << 211         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
185 }                                                 212 }
186                                                   213 
187 //....oooOO0OOooo........oooOO0OOooo........oo << 214 // BuildPhysicsTable for the Rayleigh Scattering process
                                                   >> 215 // --------------------------------------------------------
188 void G4OpRayleigh::BuildPhysicsTable(const G4P    216 void G4OpRayleigh::BuildPhysicsTable(const G4ParticleDefinition&)
189 {                                                 217 {
190   if(thePhysicsTable)                          << 218   if (thePhysicsTable) {
191   {                                            << 219      thePhysicsTable->clearAndDestroy();
192     // thePhysicsTable->clearAndDestroy();     << 220      delete thePhysicsTable;
193     delete thePhysicsTable;                    << 221      thePhysicsTable = NULL;
194     thePhysicsTable = nullptr;                 << 
195   }                                               222   }
196                                                   223 
197   const G4MaterialTable* theMaterialTable = G4    224   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
198   const size_t numOfMaterials             = G4 << 225   const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
199   thePhysicsTable                         = ne << 
200                                                   226 
201   for(size_t i = 0; i < numOfMaterials; ++i)   << 227   thePhysicsTable = new G4PhysicsTable( numOfMaterials );
202   {                                            << 228   
203     G4Material* material               = (*the << 229   for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
204     G4MaterialPropertiesTable* matProp = mater << 230   {
205     G4PhysicsFreeVector* rayleigh = nullptr;   << 231       G4Material* material = (*theMaterialTable)[iMaterial];
206     if(matProp)                                << 232       G4MaterialPropertiesTable* materialProperties = 
207     {                                          << 233                                        material->GetMaterialPropertiesTable();
208       rayleigh = matProp->GetProperty(kRAYLEIG << 234       G4PhysicsOrderedFreeVector* rayleigh = NULL;
209       if(rayleigh == nullptr)                  << 235       if ( materialProperties != NULL ) {
210         rayleigh = CalculateRayleighMeanFreePa << 236          rayleigh = materialProperties->GetProperty( kRAYLEIGH );
211     }                                          << 237          if ( rayleigh == NULL ) rayleigh = 
212     thePhysicsTable->insertAt(i, rayleigh);    << 238                                    CalculateRayleighMeanFreePaths( material );
                                                   >> 239       }
                                                   >> 240       thePhysicsTable->insertAt( iMaterial, rayleigh );
213   }                                               241   }
214 }                                                 242 }
215                                                   243 
216 //....oooOO0OOooo........oooOO0OOooo........oo << 244 // GetMeanFreePath()
217 G4double G4OpRayleigh::GetMeanFreePath(const G << 245 // -----------------
218                                        G4Force << 246 //
                                                   >> 247 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
                                                   >> 248                                        G4double ,
                                                   >> 249                                        G4ForceCondition* )
219 {                                                 250 {
220   auto rayleigh = static_cast<G4PhysicsFreeVec << 251   const G4DynamicParticle* particle = aTrack.GetDynamicParticle();
221       (*thePhysicsTable)(aTrack.GetMaterial()- << 252   const G4double photonMomentum = particle->GetTotalMomentum();
222                                                << 253   const G4Material* material = aTrack.GetMaterial();
                                                   >> 254 
                                                   >> 255   G4PhysicsOrderedFreeVector* rayleigh = 
                                                   >> 256                               static_cast<G4PhysicsOrderedFreeVector*>
                                                   >> 257                               ((*thePhysicsTable)(material->GetIndex()));
                                                   >> 258   
223   G4double rsLength = DBL_MAX;                    259   G4double rsLength = DBL_MAX;
224   if(rayleigh)                                 << 260   if( rayleigh != NULL ) rsLength = rayleigh->Value( photonMomentum );
225   {                                            << 
226     rsLength = rayleigh->Value(aTrack.GetDynam << 
227                                idx_rslength);  << 
228   }                                            << 
229   return rsLength;                                261   return rsLength;
230 }                                                 262 }
231                                                   263 
232 //....oooOO0OOooo........oooOO0OOooo........oo << 264 // CalculateRayleighMeanFreePaths()
233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa << 265 // --------------------------------
234   const G4Material* material) const            << 266 // Private method to compute Rayleigh Scattering Lengths
                                                   >> 267 G4PhysicsOrderedFreeVector* 
                                                   >> 268 G4OpRayleigh::CalculateRayleighMeanFreePaths( const G4Material* material ) const
235 {                                                 269 {
236   G4MaterialPropertiesTable* MPT = material->G << 270   G4MaterialPropertiesTable* materialProperties = 
                                                   >> 271                                        material->GetMaterialPropertiesTable();
237                                                   272 
238   // Retrieve the beta_T or isothermal compres    273   // Retrieve the beta_T or isothermal compressibility value. For backwards
239   // compatibility use a constant if the mater    274   // compatibility use a constant if the material is "Water". If the material
240   // doesn't have an ISOTHERMAL_COMPRESSIBILIT    275   // doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
241   G4double betat;                                 276   G4double betat;
242   if(material->GetName() == "Water")           << 277   if ( material->GetName() == "Water" )
243   {                                            << 278     betat = 7.658e-23*m3/MeV;
244     betat = 7.658e-23 * m3 / MeV;              << 279   else if(materialProperties->ConstPropertyExists("ISOTHERMAL_COMPRESSIBILITY"))
245   }                                            << 280     betat = materialProperties->GetConstProperty(kISOTHERMAL_COMPRESSIBILITY);
246   else if(MPT->ConstPropertyExists(kISOTHERMAL << 
247   {                                            << 
248     betat = MPT->GetConstProperty(kISOTHERMAL_ << 
249   }                                            << 
250   else                                            281   else
251   {                                            << 282     return NULL;
252     return nullptr;                            << 
253   }                                            << 
254                                                   283 
255   // If the material doesn't have a RINDEX pro    284   // If the material doesn't have a RINDEX property vector then return
256   G4MaterialPropertyVector* rIndex = MPT->GetP << 285   G4MaterialPropertyVector* rIndex = materialProperties->GetProperty(kRINDEX);
257   if(rIndex == nullptr)                        << 286   if ( rIndex == NULL ) return NULL;
258     return nullptr;                            << 
259                                                   287 
260   // Retrieve the optional scale factor (scale << 288   // Retrieve the optional scale factor, (this just scales the scattering length
261   G4double scaleFactor = 1.0;                     289   G4double scaleFactor = 1.0;
262   if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR << 290   if( materialProperties->ConstPropertyExists( "RS_SCALE_FACTOR" ) )
263   {                                            << 291     scaleFactor= materialProperties->GetConstProperty(kRS_SCALE_FACTOR );
264     scaleFactor = MPT->GetConstProperty(kRS_SC << 
265   }                                            << 
266                                                   292 
267   // Retrieve the material temperature. For ba << 293   // Retrieve the material temperature. For backwards compatibility use a 
268   // constant if the material is "Water"          294   // constant if the material is "Water"
269   G4double temperature;                           295   G4double temperature;
270   if(material->GetName() == "Water")           << 296   if( material->GetName() == "Water" )
271   {                                            << 297     temperature = 283.15*kelvin; // Temperature of water is 10 degrees celsius
272     temperature =                              << 
273       283.15 * kelvin;  // Temperature of wate << 
274   }                                            << 
275   else                                            298   else
276   {                                            << 
277     temperature = material->GetTemperature();     299     temperature = material->GetTemperature();
278   }                                            << 
279                                                   300 
280   auto rayleighMFPs = new G4PhysicsFreeVector( << 301   G4PhysicsOrderedFreeVector* rayleighMeanFreePaths =
                                                   >> 302                                              new G4PhysicsOrderedFreeVector();
281   // This calculates the meanFreePath via the     303   // This calculates the meanFreePath via the Einstein-Smoluchowski formula
282   const G4double c1 =                          << 304   const G4double c1 = scaleFactor * betat * temperature * k_Boltzmann / 
283     scaleFactor * betat * temperature * k_Bolt << 305                       ( 6.0 * pi );
284                                                   306 
285   for(size_t uRIndex = 0; uRIndex < rIndex->Ge << 307   for( size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); uRIndex++ )
286   {                                               308   {
287     const G4double energy        = rIndex->Ene << 309      const G4double energy = rIndex->Energy( uRIndex );
288     const G4double rIndexSquared = (*rIndex)[u << 310      const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
289     const G4double xlambda       = h_Planck *  << 311      const G4double xlambda = h_Planck * c_light / energy;
290     const G4double c2            = std::pow(tw << 312      const G4double c2 = std::pow(twopi/xlambda,4);
291     const G4double c3 =                        << 313      const G4double c3 = 
292       std::pow(((rIndexSquared - 1.0) * (rInde << 314                     std::pow(((rIndexSquared-1.0)*(rIndexSquared+2.0 )/3.0),2);
293                                                << 
294     const G4double meanFreePath = 1.0 / (c1 *  << 
295                                                << 
296     if(verboseLevel > 0)                       << 
297     {                                          << 
298       G4cout << energy << "MeV\t" << meanFreeP << 
299     }                                          << 
300                                                   315 
301     rayleighMFPs->InsertValues(energy, meanFre << 316      const G4double meanFreePath = 1.0 / ( c1 * c2 * c3 );
302   }                                            << 
303                                                   317 
304   return rayleighMFPs;                         << 318      if( verboseLevel>0 )
305 }                                              << 319        G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl;
306                                                   320 
307 //....oooOO0OOooo........oooOO0OOooo........oo << 321      rayleighMeanFreePaths->InsertValues( energy, meanFreePath );
308 void G4OpRayleigh::SetVerboseLevel(G4int verbo << 322   }
309 {                                              << 323 
310   verboseLevel = verbose;                      << 324   return rayleighMeanFreePaths;
311   G4OpticalParameters::Instance()->SetRayleigh << 
312 }                                                 325 }
313                                                   326