Geant4  10.00.p02
G4LowEnergyRayleigh.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // --------------------------------------------------------------------
27 //
28 // $Id$
29 // GEANT4 tag $Name: geant4-09-01-ref-00 $
30 //
31 // Author: A. Forti
32 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
33 //
34 // History:
35 // --------
36 // Added Livermore data table construction methods A. Forti
37 // Added BuildMeanFreePath A. Forti
38 // Added PostStepDoIt A. Forti
39 // Added SelectRandomAtom A. Forti
40 // Added map of the elements A.Forti
41 // 24.04.01 V.Ivanchenko remove RogueWave
42 // 11.08.2001 MGP - Major revision according to a design iteration
43 // 06.10.2001 MGP - Added strategy to test range for secondary generation
44 // 05.06.2002 F.Longo and G.Depaola - bug fixed in angular distribution
45 // 20.10.2002 G. Depaola - Change sampling method of theta
46 // 22.01.2003 V.Ivanchenko - Cut per region
47 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
48 //
49 // --------------------------------------------------------------------
50 
51 #include "G4LowEnergyRayleigh.hh"
52 #include "Randomize.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4ParticleDefinition.hh"
56 #include "G4Track.hh"
57 #include "G4Step.hh"
58 #include "G4ForceCondition.hh"
59 #include "G4Gamma.hh"
60 #include "G4Electron.hh"
61 #include "G4DynamicParticle.hh"
62 #include "G4VParticleChange.hh"
63 #include "G4ThreeVector.hh"
66 #include "G4RDVEMDataSet.hh"
68 #include "G4RDVDataSetAlgorithm.hh"
70 
71 #include "G4MaterialCutsCouple.hh"
72 
74  : G4VDiscreteProcess(processName),
75  lowEnergyLimit(250*eV),
76  highEnergyLimit(100*GeV),
77  intrinsicLowEnergyLimit(10*eV),
78  intrinsicHighEnergyLimit(100*GeV)
79 {
82  {
83  G4Exception("G4LowEnergyRayleigh::G4LowEnergyRayleigh()",
84  "OutOfRange", FatalException,
85  "Energy limit outside intrinsic process validity range!");
86  }
87 
89 
90  G4RDVDataSetAlgorithm* ffInterpolation = new G4RDLogLogInterpolation;
91  G4String formFactorFile = "rayl/re-ff-";
92  formFactorData = new G4RDCompositeEMDataSet(ffInterpolation,1.,1.);
93  formFactorData->LoadData(formFactorFile);
94 
96 
97  if (verboseLevel > 0)
98  {
99  G4cout << GetProcessName() << " is created " << G4endl
100  << "Energy range: "
101  << lowEnergyLimit / keV << " keV - "
102  << highEnergyLimit / GeV << " GeV"
103  << G4endl;
104  }
105 }
106 
108 {
109  delete meanFreePathTable;
110  delete crossSectionHandler;
111  delete formFactorData;
112 }
113 
115 {
116 
118  G4String crossSectionFile = "rayl/re-cs-";
119  crossSectionHandler->LoadData(crossSectionFile);
120 
121  delete meanFreePathTable;
123 }
124 
126  const G4Step& aStep)
127 {
128 
129  aParticleChange.Initialize(aTrack);
130 
131  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
132  G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
133 
134  if (photonEnergy0 <= lowEnergyLimit)
135  {
139  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
140  }
141 
142  // G4double e0m = photonEnergy0 / electron_mass_c2 ;
143  G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
144 
145  // Select randomly one element in the current material
146  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
147  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
148 
149  // Sample the angle of the scattered photon
150 
151  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
152 
153  G4double gReject,x,dataFormFactor;
154  G4double randomFormFactor;
155  G4double cosTheta;
156  G4double sinTheta;
157  G4double fcostheta;
158 
159  do
160  {
161  do
162  {
163  cosTheta = 2. * G4UniformRand() - 1.;
164  fcostheta = ( 1. + cosTheta*cosTheta)/2.;
165  } while (fcostheta < G4UniformRand());
166 
167  G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
168  x = sinThetaHalf / (wlPhoton/cm);
169  if (x > 1.e+005)
170  dataFormFactor = formFactorData->FindValue(x,Z-1);
171  else
172  dataFormFactor = formFactorData->FindValue(0.,Z-1);
173  randomFormFactor = G4UniformRand() * Z * Z;
174  sinTheta = std::sqrt(1. - cosTheta*cosTheta);
175  gReject = dataFormFactor * dataFormFactor;
176 
177  } while( gReject < randomFormFactor);
178 
179  // Scattered photon angles. ( Z - axis along the parent photon)
180  G4double phi = twopi * G4UniformRand() ;
181  G4double dirX = sinTheta*std::cos(phi);
182  G4double dirY = sinTheta*std::sin(phi);
183  G4double dirZ = cosTheta;
184 
185  // Update G4VParticleChange for the scattered photon
186  G4ThreeVector photonDirection1(dirX, dirY, dirZ);
187 
188  photonDirection1.rotateUz(photonDirection0);
189  aParticleChange.ProposeEnergy(photonEnergy0);
190  aParticleChange.ProposeMomentumDirection(photonDirection1);
191 
193 
194  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
195 }
196 
198 {
199  return ( &particle == G4Gamma::Gamma() );
200 }
201 
203  G4double, // previousStepSize
205 {
207  G4double energy = photon->GetKineticEnergy();
208  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
209  size_t materialIndex = couple->GetIndex();
210 
211  G4double meanFreePath;
212  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
213  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
214  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
215  return meanFreePath;
216 }
static const double cm
Definition: G4SIunits.hh:106
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4double intrinsicLowEnergyLimit
G4RDVEMDataSet * formFactorData
const G4DynamicParticle * GetDynamicParticle() const
G4RDVEMDataSet * meanFreePathTable
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4RDVEMDataSet * BuildMeanFreePathForMaterials(const G4DataVector *energyCuts=0)
const G4double intrinsicHighEnergyLimit
int G4int
Definition: G4Types.hh:78
G4LowEnergyRayleigh(const G4String &processName="LowEnRayleigh")
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4bool IsApplicable(const G4ParticleDefinition &)
static const double GeV
Definition: G4SIunits.hh:196
Definition: G4Step.hh:76
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
void BuildPhysicsTable(const G4ParticleDefinition &photon)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void Initialize(const G4Track &)
static const double eV
Definition: G4SIunits.hh:194
G4double energy(const ThreeVector &p, const G4double m)
void SetNumberOfSecondaries(G4int totSecondaries)
G4int SelectRandomAtom(const G4MaterialCutsCouple *couple, G4double e) const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4RDVCrossSectionHandler * crossSectionHandler
#define G4endl
Definition: G4ios.hh:61
void LoadData(const G4String &dataFile)
static const double keV
Definition: G4SIunits.hh:195
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
G4ThreeVector G4ParticleMomentum
#define DBL_MAX
Definition: templates.hh:83
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
virtual G4bool LoadData(const G4String &fileName)=0