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


** Warning: Cannot open xref database.

  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 //                                                
 28 //                                                
 29 //////////////////////////////////////////////    
 30 // Optical Photon Rayleigh Scattering Class Im    
 31 //////////////////////////////////////////////    
 32 //                                                
 33 // File:        G4OpRayleigh.cc                   
 34 // Description: Discrete Process -- Rayleigh s    
 35 //    photons                                     
 36 // Version:     1.0                               
 37 // Created:     1996-05-31                        
 38 // Author:      Juliet Armstrong                  
 39 // Updated:     2014-10-10 -  This version cal    
 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    
 50 //              2001-09-18 by mma                 
 51 //              >numOfMaterials=G4Material::Ge    
 52 //              2001-01-30 by Peter Gumplinger    
 53 //              > allow for positiv and negati    
 54 //              > new momentum direction to be    
 55 //              > new and old polarization vec    
 56 //              2001-01-29 by Peter Gumplinger    
 57 //              > fix calculation of SinTheta     
 58 //              1997-04-09 by Peter Gumplinger    
 59 //              > new physics/tracking scheme     
 60 //                                                
 61 //////////////////////////////////////////////    
 62                                                   
 63 #include "G4OpRayleigh.hh"                        
 64 #include "G4ios.hh"                               
 65 #include "G4PhysicalConstants.hh"                 
 66 #include "G4SystemOfUnits.hh"                     
 67 #include "G4OpticalParameters.hh"                 
 68 #include "G4OpProcessSubType.hh"                  
 69                                                   
 70 //....oooOO0OOooo........oooOO0OOooo........oo    
 71 G4OpRayleigh::G4OpRayleigh(const G4String& pro    
 72   : G4VDiscreteProcess(processName, type)         
 73 {                                                 
 74   Initialise();                                   
 75   SetProcessSubType(fOpRayleigh);                 
 76   thePhysicsTable = nullptr;                      
 77                                                   
 78   if(verboseLevel > 0)                            
 79   {                                               
 80     G4cout << GetProcessName() << " is created    
 81   }                                               
 82 }                                                 
 83                                                   
 84 //....oooOO0OOooo........oooOO0OOooo........oo    
 85 G4OpRayleigh::~G4OpRayleigh()                     
 86 {                                                 
 87   // VI: inside this PhysicsTable all properti    
 88   //     it is not possible to destroy            
 89   delete thePhysicsTable;                         
 90 }                                                 
 91                                                   
 92 //....oooOO0OOooo........oooOO0OOooo........oo    
 93 void G4OpRayleigh::PreparePhysicsTable(const G    
 94 {                                                 
 95   Initialise();                                   
 96 }                                                 
 97                                                   
 98 //....oooOO0OOooo........oooOO0OOooo........oo    
 99 void G4OpRayleigh::Initialise()                   
100 {                                                 
101   SetVerboseLevel(G4OpticalParameters::Instanc    
102 }                                                 
103                                                   
104 //....oooOO0OOooo........oooOO0OOooo........oo    
105 G4VParticleChange* G4OpRayleigh::PostStepDoIt(    
106                                                   
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                                                   
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                                                   
184   return G4VDiscreteProcess::PostStepDoIt(aTra    
185 }                                                 
186                                                   
187 //....oooOO0OOooo........oooOO0OOooo........oo    
188 void G4OpRayleigh::BuildPhysicsTable(const G4P    
189 {                                                 
190   if(thePhysicsTable)                             
191   {                                               
192     // thePhysicsTable->clearAndDestroy();        
193     delete thePhysicsTable;                       
194     thePhysicsTable = nullptr;                    
195   }                                               
196                                                   
197   const G4MaterialTable* theMaterialTable = G4    
198   const size_t numOfMaterials             = G4    
199   thePhysicsTable                         = ne    
200                                                   
201   for(size_t i = 0; i < numOfMaterials; ++i)      
202   {                                               
203     G4Material* material               = (*the    
204     G4MaterialPropertiesTable* matProp = mater    
205     G4PhysicsFreeVector* rayleigh = nullptr;      
206     if(matProp)                                   
207     {                                             
208       rayleigh = matProp->GetProperty(kRAYLEIG    
209       if(rayleigh == nullptr)                     
210         rayleigh = CalculateRayleighMeanFreePa    
211     }                                             
212     thePhysicsTable->insertAt(i, rayleigh);       
213   }                                               
214 }                                                 
215                                                   
216 //....oooOO0OOooo........oooOO0OOooo........oo    
217 G4double G4OpRayleigh::GetMeanFreePath(const G    
218                                        G4Force    
219 {                                                 
220   auto rayleigh = static_cast<G4PhysicsFreeVec    
221       (*thePhysicsTable)(aTrack.GetMaterial()-    
222                                                   
223   G4double rsLength = DBL_MAX;                    
224   if(rayleigh)                                    
225   {                                               
226     rsLength = rayleigh->Value(aTrack.GetDynam    
227                                idx_rslength);     
228   }                                               
229   return rsLength;                                
230 }                                                 
231                                                   
232 //....oooOO0OOooo........oooOO0OOooo........oo    
233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa    
234   const G4Material* material) const               
235 {                                                 
236   G4MaterialPropertiesTable* MPT = material->G    
237                                                   
238   // Retrieve the beta_T or isothermal compres    
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                                                   
301     rayleighMFPs->InsertValues(energy, meanFre    
302   }                                               
303                                                   
304   return rayleighMFPs;                            
305 }                                                 
306                                                   
307 //....oooOO0OOooo........oooOO0OOooo........oo    
308 void G4OpRayleigh::SetVerboseLevel(G4int verbo    
309 {                                                 
310   verboseLevel = verbose;                         
311   G4OpticalParameters::Instance()->SetRayleigh    
312 }                                                 
313