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 7.1)


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