ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEnergyRayleigh.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LowEnergyRayleigh.cc
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 //
29 // Author: A. Forti
30 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31 //
32 // History:
33 // --------
34 // Added Livermore data table construction methods A. Forti
35 // Added BuildMeanFreePath A. Forti
36 // Added PostStepDoIt A. Forti
37 // Added SelectRandomAtom A. Forti
38 // Added map of the elements A.Forti
39 // 24.04.01 V.Ivanchenko remove RogueWave
40 // 11.08.2001 MGP - Major revision according to a design iteration
41 // 06.10.2001 MGP - Added strategy to test range for secondary generation
42 // 05.06.2002 F.Longo and G.Depaola - bug fixed in angular distribution
43 // 20.10.2002 G. Depaola - Change sampling method of theta
44 // 22.01.2003 V.Ivanchenko - Cut per region
45 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
46 //
47 // --------------------------------------------------------------------
48 
49 #include "G4LowEnergyRayleigh.hh"
50 #include "Randomize.hh"
51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4Track.hh"
55 #include "G4Step.hh"
56 #include "G4ForceCondition.hh"
57 #include "G4Gamma.hh"
58 #include "G4Electron.hh"
59 #include "G4DynamicParticle.hh"
60 #include "G4VParticleChange.hh"
61 #include "G4ThreeVector.hh"
64 #include "G4RDVEMDataSet.hh"
66 #include "G4RDVDataSetAlgorithm.hh"
68 
69 #include "G4MaterialCutsCouple.hh"
70 
72  : G4VDiscreteProcess(processName),
73  lowEnergyLimit(250*eV),
74  highEnergyLimit(100*GeV),
75  intrinsicLowEnergyLimit(10*eV),
76  intrinsicHighEnergyLimit(100*GeV)
77 {
80  {
81  G4Exception("G4LowEnergyRayleigh::G4LowEnergyRayleigh()",
82  "OutOfRange", FatalException,
83  "Energy limit outside intrinsic process validity range!");
84  }
85 
87 
88  G4RDVDataSetAlgorithm* ffInterpolation = new G4RDLogLogInterpolation;
89  G4String formFactorFile = "rayl/re-ff-";
90  formFactorData = new G4RDCompositeEMDataSet(ffInterpolation,1.,1.);
91  formFactorData->LoadData(formFactorFile);
92 
94 
95  if (verboseLevel > 0)
96  {
97  G4cout << GetProcessName() << " is created " << G4endl
98  << "Energy range: "
99  << lowEnergyLimit / keV << " keV - "
100  << highEnergyLimit / GeV << " GeV"
101  << G4endl;
102  }
103 }
104 
106 {
107  delete meanFreePathTable;
108  delete crossSectionHandler;
109  delete formFactorData;
110 }
111 
113 {
114 
116  G4String crossSectionFile = "rayl/re-cs-";
117  crossSectionHandler->LoadData(crossSectionFile);
118 
119  delete meanFreePathTable;
121 }
122 
124  const G4Step& aStep)
125 {
126 
127  aParticleChange.Initialize(aTrack);
128 
129  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
130  G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
131 
132  if (photonEnergy0 <= lowEnergyLimit)
133  {
137  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
138  }
139 
140  // G4double e0m = photonEnergy0 / electron_mass_c2 ;
141  G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
142 
143  // Select randomly one element in the current material
144  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
145  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
146 
147  // Sample the angle of the scattered photon
148 
149  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
150 
151  G4double gReject,x,dataFormFactor;
152  G4double randomFormFactor;
153  G4double cosTheta;
154  G4double sinTheta;
155  G4double fcostheta;
156 
157  do
158  {
159  do
160  {
161  cosTheta = 2. * G4UniformRand() - 1.;
162  fcostheta = ( 1. + cosTheta*cosTheta)/2.;
163  } while (fcostheta < G4UniformRand());
164 
165  G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
166  x = sinThetaHalf / (wlPhoton/cm);
167  if (x > 1.e+005)
168  dataFormFactor = formFactorData->FindValue(x,Z-1);
169  else
170  dataFormFactor = formFactorData->FindValue(0.,Z-1);
171  randomFormFactor = G4UniformRand() * Z * Z;
172  sinTheta = std::sqrt(1. - cosTheta*cosTheta);
173  gReject = dataFormFactor * dataFormFactor;
174 
175  } while( gReject < randomFormFactor);
176 
177  // Scattered photon angles. ( Z - axis along the parent photon)
179  G4double dirX = sinTheta*std::cos(phi);
180  G4double dirY = sinTheta*std::sin(phi);
181  G4double dirZ = cosTheta;
182 
183  // Update G4VParticleChange for the scattered photon
184  G4ThreeVector photonDirection1(dirX, dirY, dirZ);
185 
186  photonDirection1.rotateUz(photonDirection0);
187  aParticleChange.ProposeEnergy(photonEnergy0);
188  aParticleChange.ProposeMomentumDirection(photonDirection1);
189 
191 
192  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
193 }
194 
196 {
197  return ( &particle == G4Gamma::Gamma() );
198 }
199 
201  G4double, // previousStepSize
203 {
205  G4double energy = photon->GetKineticEnergy();
206  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
207  size_t materialIndex = couple->GetIndex();
208 
209  G4double meanFreePath;
210  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
211  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
212  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
213  return meanFreePath;
214 }