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.2)


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