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.1.p1)


  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 84717 2014-10-20 07:39:47Z 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         do {
                                                   >> 143            // Try to simulate the scattered photon momentum direction
                                                   >> 144            // w.r.t. the initial photon momentum direction
                                                   >> 145 
                                                   >> 146            G4double CosTheta = G4UniformRand();
                                                   >> 147            G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
                                                   >> 148            // consider for the angle 90-180 degrees
                                                   >> 149            if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
                                                   >> 150 
                                                   >> 151            // simulate the phi angle
                                                   >> 152            G4double rand = twopi*G4UniformRand();
                                                   >> 153            G4double SinPhi = std::sin(rand);
                                                   >> 154            G4double CosPhi = std::cos(rand);
                                                   >> 155 
                                                   >> 156            // start constructing the new momentum direction
                                                   >> 157      G4double unit_x = SinTheta * CosPhi; 
                                                   >> 158      G4double unit_y = SinTheta * SinPhi;  
                                                   >> 159      G4double unit_z = CosTheta; 
                                                   >> 160      NewMomentumDirection.set (unit_x,unit_y,unit_z);
                                                   >> 161 
                                                   >> 162            // Rotate the new momentum direction into global reference system
                                                   >> 163            OldMomentumDirection = aParticle->GetMomentumDirection();
                                                   >> 164            OldMomentumDirection = OldMomentumDirection.unit();
                                                   >> 165            NewMomentumDirection.rotateUz(OldMomentumDirection);
                                                   >> 166            NewMomentumDirection = NewMomentumDirection.unit();
                                                   >> 167 
                                                   >> 168            // calculate the new polarization direction
                                                   >> 169            // The new polarization needs to be in the same plane as the new
                                                   >> 170            // momentum direction and the old polarization direction
                                                   >> 171            OldPolarization = aParticle->GetPolarization();
                                                   >> 172            G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
                                                   >> 173 
                                                   >> 174            NewPolarization = NewMomentumDirection + constant*OldPolarization;
                                                   >> 175            NewPolarization = NewPolarization.unit();
                                                   >> 176 
                                                   >> 177            // There is a corner case, where the Newmomentum direction
                                                   >> 178            // is the same as oldpolariztion direction:
                                                   >> 179            // random generate the azimuthal angle w.r.t. Newmomentum direction
                                                   >> 180            if (NewPolarization.mag() == 0.) {
                                                   >> 181               rand = G4UniformRand()*twopi;
                                                   >> 182               NewPolarization.set(std::cos(rand),std::sin(rand),0.);
                                                   >> 183               NewPolarization.rotateUz(NewMomentumDirection);
                                                   >> 184            } else {
                                                   >> 185               // There are two directions which are perpendicular
                                                   >> 186               // to the new momentum direction
                                                   >> 187               if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
                                                   >> 188            }
                                                   >> 189     
                                                   >> 190      // simulate according to the distribution cos^2(theta)
                                                   >> 191            cosTheta = NewPolarization.dot(OldPolarization);
                                                   >> 192         } while (std::pow(cosTheta,2) < G4UniformRand());
                                                   >> 193 
                                                   >> 194         aParticleChange.ProposePolarization(NewPolarization);
                                                   >> 195         aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
                                                   >> 196 
                                                   >> 197         if (verboseLevel>0) {
                                                   >> 198                 G4cout << "New Polarization: " 
                                                   >> 199                      << NewPolarization << G4endl;
                                                   >> 200                 G4cout << "Polarization Change: "
                                                   >> 201                      << *(aParticleChange.GetPolarization()) << G4endl;  
                                                   >> 202                 G4cout << "New Momentum Direction: " 
                                                   >> 203                      << NewMomentumDirection << G4endl;
                                                   >> 204                 G4cout << "Momentum Change: "
                                                   >> 205                      << *(aParticleChange.GetMomentumDirection()) << G4endl; 
                                                   >> 206         }
183                                                   207 
184   return G4VDiscreteProcess::PostStepDoIt(aTra << 208         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
185 }                                                 209 }
186                                                   210 
187 //....oooOO0OOooo........oooOO0OOooo........oo << 211 // BuildPhysicsTable for the Rayleigh Scattering process
                                                   >> 212 // --------------------------------------------------------
188 void G4OpRayleigh::BuildPhysicsTable(const G4P    213 void G4OpRayleigh::BuildPhysicsTable(const G4ParticleDefinition&)
189 {                                                 214 {
190   if(thePhysicsTable)                          << 215   if (thePhysicsTable) {
191   {                                            << 216      thePhysicsTable->clearAndDestroy();
192     // thePhysicsTable->clearAndDestroy();     << 217      delete thePhysicsTable;
193     delete thePhysicsTable;                    << 218      thePhysicsTable = NULL;
194     thePhysicsTable = nullptr;                 << 
195   }                                               219   }
196                                                   220 
197   const G4MaterialTable* theMaterialTable = G4    221   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
198   const size_t numOfMaterials             = G4 << 222   const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
199   thePhysicsTable                         = ne << 
200                                                   223 
201   for(size_t i = 0; i < numOfMaterials; ++i)   << 224   thePhysicsTable = new G4PhysicsTable( numOfMaterials );
202   {                                            << 225   
203     G4Material* material               = (*the << 226   for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
204     G4MaterialPropertiesTable* matProp = mater << 227   {
205     G4PhysicsFreeVector* rayleigh = nullptr;   << 228       G4Material* material = (*theMaterialTable)[iMaterial];
206     if(matProp)                                << 229       G4MaterialPropertiesTable* materialProperties = 
207     {                                          << 230                                        material->GetMaterialPropertiesTable();
208       rayleigh = matProp->GetProperty(kRAYLEIG << 231       G4PhysicsOrderedFreeVector* rayleigh = NULL;
209       if(rayleigh == nullptr)                  << 232       if ( materialProperties != NULL ) {
210         rayleigh = CalculateRayleighMeanFreePa << 233          rayleigh = materialProperties->GetProperty( "RAYLEIGH" );
211     }                                          << 234          if ( rayleigh == NULL ) rayleigh = 
212     thePhysicsTable->insertAt(i, rayleigh);    << 235                                    CalculateRayleighMeanFreePaths( material );
                                                   >> 236       }
                                                   >> 237       thePhysicsTable->insertAt( iMaterial, rayleigh );
213   }                                               238   }
214 }                                                 239 }
215                                                   240 
216 //....oooOO0OOooo........oooOO0OOooo........oo << 241 // GetMeanFreePath()
217 G4double G4OpRayleigh::GetMeanFreePath(const G << 242 // -----------------
218                                        G4Force << 243 //
                                                   >> 244 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
                                                   >> 245                                        G4double ,
                                                   >> 246                                        G4ForceCondition* )
219 {                                                 247 {
220   auto rayleigh = static_cast<G4PhysicsFreeVec << 248   const G4DynamicParticle* particle = aTrack.GetDynamicParticle();
221       (*thePhysicsTable)(aTrack.GetMaterial()- << 249   const G4double photonMomentum = particle->GetTotalMomentum();
222                                                << 250   const G4Material* material = aTrack.GetMaterial();
                                                   >> 251 
                                                   >> 252   G4PhysicsOrderedFreeVector* rayleigh = 
                                                   >> 253                               static_cast<G4PhysicsOrderedFreeVector*>
                                                   >> 254                               ((*thePhysicsTable)(material->GetIndex()));
                                                   >> 255   
223   G4double rsLength = DBL_MAX;                    256   G4double rsLength = DBL_MAX;
224   if(rayleigh)                                 << 257   if( rayleigh != NULL ) rsLength = rayleigh->Value( photonMomentum );
225   {                                            << 
226     rsLength = rayleigh->Value(aTrack.GetDynam << 
227                                idx_rslength);  << 
228   }                                            << 
229   return rsLength;                                258   return rsLength;
230 }                                                 259 }
231                                                   260 
232 //....oooOO0OOooo........oooOO0OOooo........oo << 261 // CalculateRayleighMeanFreePaths()
233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa << 262 // --------------------------------
234   const G4Material* material) const            << 263 // Private method to compute Rayleigh Scattering Lengths
                                                   >> 264 G4PhysicsOrderedFreeVector* 
                                                   >> 265 G4OpRayleigh::CalculateRayleighMeanFreePaths( const G4Material* material ) const
235 {                                                 266 {
236   G4MaterialPropertiesTable* MPT = material->G << 267   G4MaterialPropertiesTable* materialProperties = 
                                                   >> 268                                        material->GetMaterialPropertiesTable();
237                                                   269 
238   // Retrieve the beta_T or isothermal compres    270   // Retrieve the beta_T or isothermal compressibility value. For backwards
239   // compatibility use a constant if the mater    271   // compatibility use a constant if the material is "Water". If the material
240   // doesn't have an ISOTHERMAL_COMPRESSIBILIT    272   // doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
241   G4double betat;                                 273   G4double betat;
242   if(material->GetName() == "Water")           << 274   if ( material->GetName() == "Water" )
243   {                                            << 275     betat = 7.658e-23*m3/MeV;
244     betat = 7.658e-23 * m3 / MeV;              << 276   else if(materialProperties->ConstPropertyExists("ISOTHERMAL_COMPRESSIBILITY"))
245   }                                            << 277     betat = materialProperties->GetConstProperty("ISOTHERMAL_COMPRESSIBILITY");
246   else if(MPT->ConstPropertyExists(kISOTHERMAL << 
247   {                                            << 
248     betat = MPT->GetConstProperty(kISOTHERMAL_ << 
249   }                                            << 
250   else                                            278   else
251   {                                            << 279     return NULL;
252     return nullptr;                            << 
253   }                                            << 
254                                                   280 
255   // If the material doesn't have a RINDEX pro    281   // If the material doesn't have a RINDEX property vector then return
256   G4MaterialPropertyVector* rIndex = MPT->GetP << 282   G4MaterialPropertyVector* rIndex = materialProperties->GetProperty("RINDEX");
257   if(rIndex == nullptr)                        << 283   if ( rIndex == NULL ) return NULL;
258     return nullptr;                            << 
259                                                   284 
260   // Retrieve the optional scale factor (scale << 285   // Retrieve the optional scale factor, (this just scales the scattering length
261   G4double scaleFactor = 1.0;                     286   G4double scaleFactor = 1.0;
262   if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR << 287   if( materialProperties->ConstPropertyExists( "RS_SCALE_FACTOR" ) )
263   {                                            << 288     scaleFactor= materialProperties->GetConstProperty("RS_SCALE_FACTOR" );
264     scaleFactor = MPT->GetConstProperty(kRS_SC << 
265   }                                            << 
266                                                   289 
267   // Retrieve the material temperature. For ba << 290   // Retrieve the material temperature. For backwards compatibility use a 
268   // constant if the material is "Water"          291   // constant if the material is "Water"
269   G4double temperature;                           292   G4double temperature;
270   if(material->GetName() == "Water")           << 293   if( material->GetName() == "Water" )
271   {                                            << 294     temperature = 283.15*kelvin; // Temperature of water is 10 degrees celsius
272     temperature =                              << 
273       283.15 * kelvin;  // Temperature of wate << 
274   }                                            << 
275   else                                            295   else
276   {                                            << 
277     temperature = material->GetTemperature();     296     temperature = material->GetTemperature();
278   }                                            << 
279                                                   297 
280   auto rayleighMFPs = new G4PhysicsFreeVector( << 298   G4PhysicsOrderedFreeVector* rayleighMeanFreePaths =
                                                   >> 299                                              new G4PhysicsOrderedFreeVector();
281   // This calculates the meanFreePath via the     300   // This calculates the meanFreePath via the Einstein-Smoluchowski formula
282   const G4double c1 =                          << 301   const G4double c1 = scaleFactor * betat * temperature * k_Boltzmann / 
283     scaleFactor * betat * temperature * k_Bolt << 302                       ( 6.0 * pi );
284                                                   303 
285   for(size_t uRIndex = 0; uRIndex < rIndex->Ge << 304   for( size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); uRIndex++ )
286   {                                               305   {
287     const G4double energy        = rIndex->Ene << 306      const G4double energy = rIndex->Energy( uRIndex );
288     const G4double rIndexSquared = (*rIndex)[u << 307      const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
289     const G4double xlambda       = h_Planck *  << 308      const G4double xlambda = h_Planck * c_light / energy;
290     const G4double c2            = std::pow(tw << 309      const G4double c2 = std::pow(twopi/xlambda,4);
291     const G4double c3 =                        << 310      const G4double c3 = 
292       std::pow(((rIndexSquared - 1.0) * (rInde << 311                     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                                                   312 
301     rayleighMFPs->InsertValues(energy, meanFre << 313      const G4double meanFreePath = 1.0 / ( c1 * c2 * c3 );
302   }                                            << 
303                                                   314 
304   return rayleighMFPs;                         << 315      if( verboseLevel>0 )
305 }                                              << 316        G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl;
306                                                   317 
307 //....oooOO0OOooo........oooOO0OOooo........oo << 318      rayleighMeanFreePaths->InsertValues( energy, meanFreePath );
308 void G4OpRayleigh::SetVerboseLevel(G4int verbo << 319   }
309 {                                              << 320 
310   verboseLevel = verbose;                      << 321   return rayleighMeanFreePaths;
311   G4OpticalParameters::Instance()->SetRayleigh << 
312 }                                                 322 }
313                                                   323