| Geant4 Cross Reference |
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 // Author: Sebastien Incerti << 26 // $Id: G4LivermoreRayleighModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
27 // 31 March 2012 << 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $
28 // on base of G4LivermoreRayleighModel <<
29 // 28 //
30 29
31 #include "G4LivermoreRayleighModel.hh" 30 #include "G4LivermoreRayleighModel.hh"
32 31
33 #include "G4AutoLock.hh" <<
34 #include "G4EmParameters.hh" <<
35 #include "G4RayleighAngularGenerator.hh" <<
36 #include "G4SystemOfUnits.hh" <<
37 <<
38 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 33
40 namespace << 34 using namespace std;
41 { <<
42 G4Mutex LivermoreRayleighModelMutex = G4MUTEX_ <<
43 } <<
44 <<
45 G4PhysicsFreeVector* G4LivermoreRayleighModel: <<
46 G4String G4LivermoreRayleighModel::gDataDirect <<
47 35
48 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 37
50 G4LivermoreRayleighModel::G4LivermoreRayleighM << 38 G4LivermoreRayleighModel::G4LivermoreRayleighModel(const G4ParticleDefinition*,
51 { << 39 const G4String& nam)
52 fParticleChange = nullptr; << 40 :G4VEmModel(nam),isInitialised(false)
53 lowEnergyLimit = 10 * CLHEP::eV; << 41 {
54 << 42 lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
55 SetAngularDistribution(new G4RayleighAngular << 43 highEnergyLimit = 100 * GeV;
56 << 44
57 verboseLevel = 0; << 45 SetLowEnergyLimit(lowEnergyLimit);
58 // Verbosity scale for debugging purposes: << 46 SetHighEnergyLimit(highEnergyLimit);
59 // 0 = nothing << 47 //
60 // 1 = calculation of cross sections, file o << 48 verboseLevel= 0;
61 // 2 = entering in methods << 49 // Verbosity scale:
62 << 50 // 0 = nothing
63 if (verboseLevel > 0) { << 51 // 1 = warning for energy non-conservation
64 G4cout << "G4LivermoreRayleighModel is con << 52 // 2 = details of energy budget
65 } << 53 // 3 = calculation of cross sections, file openings, sampling of atoms
>> 54 // 4 = entering in methods
>> 55
>> 56 G4cout << "Livermore Rayleigh is constructed " << G4endl
>> 57 << "Energy range: "
>> 58 << lowEnergyLimit / keV << " keV - "
>> 59 << highEnergyLimit / GeV << " GeV"
>> 60 << G4endl;
66 } 61 }
67 62
68 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 64
70 G4LivermoreRayleighModel::~G4LivermoreRayleigh 65 G4LivermoreRayleighModel::~G4LivermoreRayleighModel()
71 { << 66 {
72 if (IsMaster()) { << 67 delete meanFreePathTable;
73 for (G4int i = 0; i <= maxZ; ++i) { << 68 delete crossSectionHandler;
74 if (nullptr != dataCS[i]) { << 69 delete formFactorData;
75 delete dataCS[i]; <<
76 dataCS[i] = nullptr; <<
77 } <<
78 } <<
79 } <<
80 } 70 }
81 71
82 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 73
84 void G4LivermoreRayleighModel::Initialise(cons 74 void G4LivermoreRayleighModel::Initialise(const G4ParticleDefinition* particle,
85 cons << 75 const G4DataVector& cuts)
86 { 76 {
87 if (verboseLevel > 1) { << 77 if (verboseLevel > 3)
88 G4cout << "Calling Initialise() of G4Liver << 78 G4cout << "Calling G4LivermoreRayleighModel::Initialise()" << G4endl;
89 << "Energy range: " << LowEnergyLim <<
90 << " GeV" << G4endl; <<
91 } <<
92 79
93 if (IsMaster()) { << 80 InitialiseElementSelectors(particle,cuts);
94 // Initialise element selector <<
95 InitialiseElementSelectors(particle, cuts) <<
96 <<
97 // Access to elements <<
98 const G4ElementTable* elemTable = G4Elemen <<
99 std::size_t numElems = (*elemTable).size() <<
100 for (std::size_t ie = 0; ie < numElems; ++ <<
101 const G4Element* elem = (*elemTable)[ie] <<
102 const G4int Z = std::min(maxZ, elem->Get <<
103 if (dataCS[Z] == nullptr) { <<
104 ReadData(Z); <<
105 } <<
106 } <<
107 } <<
108 if (isInitialised) { <<
109 return; <<
110 } <<
111 fParticleChange = GetParticleChangeForGamma( <<
112 isInitialised = true; <<
113 } <<
114 81
115 //....oooOO0OOooo........oooOO0OOooo........oo << 82 // Energy limits
>> 83
>> 84 if (LowEnergyLimit() < lowEnergyLimit)
>> 85 {
>> 86 G4cout << "G4LivermoreRayleighModel: low energy limit increased from " <<
>> 87 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
>> 88 SetLowEnergyLimit(lowEnergyLimit);
>> 89 }
>> 90
>> 91 if (HighEnergyLimit() > highEnergyLimit)
>> 92 {
>> 93 G4cout << "G4LivermoreRayleighModel: high energy limit decreased from " <<
>> 94 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
>> 95 SetHighEnergyLimit(highEnergyLimit);
>> 96 }
>> 97
>> 98 // Data are read for all materials
>> 99
>> 100 crossSectionHandler = new G4CrossSectionHandler;
>> 101 crossSectionHandler->Clear();
>> 102 G4String crossSectionFile = "rayl/re-cs-";
>> 103 crossSectionHandler->LoadData(crossSectionFile);
>> 104
>> 105 meanFreePathTable = 0;
>> 106 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
>> 107
>> 108 G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
>> 109 G4String formFactorFile = "rayl/re-ff-";
>> 110 formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
>> 111 formFactorData->LoadData(formFactorFile);
>> 112
>> 113 //
>> 114
>> 115 if (verboseLevel > 2)
>> 116 G4cout << "Loaded cross section files for Livermore Rayleigh model" << G4endl;
>> 117
>> 118 G4cout << "Livermore Rayleigh model is initialized " << G4endl
>> 119 << "Energy range: "
>> 120 << LowEnergyLimit() / keV << " keV - "
>> 121 << HighEnergyLimit() / GeV << " GeV"
>> 122 << G4endl;
>> 123
>> 124 if(isInitialised) return;
>> 125
>> 126 if(pParticleChange)
>> 127 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
>> 128 else
>> 129 fParticleChange = new G4ParticleChangeForGamma();
116 130
117 void G4LivermoreRayleighModel::InitialiseLocal << 131 isInitialised = true;
118 { <<
119 SetElementSelectors(masterModel->GetElementS <<
120 } <<
121 <<
122 //....oooOO0OOooo........oooOO0OOooo........oo <<
123 132
124 const G4String& G4LivermoreRayleighModel::Find <<
125 { <<
126 // no check in this method - environment var <<
127 if (gDataDirectory.empty()) { <<
128 auto param = G4EmParameters::Instance(); <<
129 std::ostringstream ost; <<
130 if (param->LivermoreDataDir() == "livermor <<
131 ost << param->GetDirLEDATA() << "/liverm <<
132 } <<
133 else { <<
134 ost << param->GetDirLEDATA() << "/epics2 <<
135 } <<
136 gDataDirectory = ost.str(); <<
137 } <<
138 return gDataDirectory; <<
139 } 133 }
140 134
141 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 136
143 void G4LivermoreRayleighModel::ReadData(const << 137 G4double G4LivermoreRayleighModel::ComputeCrossSectionPerAtom(
>> 138 const G4ParticleDefinition*,
>> 139 G4double GammaEnergy,
>> 140 G4double Z, G4double,
>> 141 G4double, G4double)
144 { 142 {
145 if (verboseLevel > 1) { << 143 if (verboseLevel > 3)
146 G4cout << "Calling ReadData() of G4Livermo << 144 G4cout << "Calling CrossSectionPerAtom() of G4LivermoreRayleighModel" << G4endl;
147 } <<
148 const G4int Z = std::min(ZZ, maxZ); <<
149 <<
150 if (nullptr != dataCS[Z]) { <<
151 return; <<
152 } <<
153 <<
154 dataCS[Z] = new G4PhysicsFreeVector(); <<
155 <<
156 std::ostringstream ostCS; <<
157 ostCS << FindDirectoryPath() << "re-cs-" << <<
158 145
159 std::ifstream finCS(ostCS.str().c_str()); << 146 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
160 << 147 return cs;
161 if (!finCS.is_open()) { <<
162 G4ExceptionDescription ed; <<
163 ed << "G4LivermoreRayleighModel data file <<
164 << G4endl; <<
165 G4Exception("G4LivermoreRayleighModel::Rea <<
166 "G4LEDATA version should be G4 <<
167 return; <<
168 } <<
169 else { <<
170 if (verboseLevel > 3) { <<
171 G4cout << "File " << ostCS.str() << " is <<
172 } <<
173 dataCS[Z]->Retrieve(finCS, true); <<
174 } <<
175 } 148 }
176 149
177 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 151
179 G4double G4LivermoreRayleighModel::ComputeCros << 152 void G4LivermoreRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
180 << 153 const G4MaterialCutsCouple* couple,
181 << 154 const G4DynamicParticle* aDynamicGamma,
>> 155 G4double,
>> 156 G4double)
182 { 157 {
183 if (verboseLevel > 1) { << 158 if (verboseLevel > 3)
184 G4cout << "G4LivermoreRayleighModel::Compu << 159 G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel" << G4endl;
185 } <<
186 160
187 if (GammaEnergy < lowEnergyLimit) { << 161 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
188 return 0.0; << 162
>> 163 if (photonEnergy0 <= lowEnergyLimit)
>> 164 {
>> 165 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 166 fParticleChange->SetProposedKineticEnergy(0.);
>> 167 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
>> 168 // SI - IS THE FOLLOWING RETURN NECESSARY ?
>> 169 return ;
189 } 170 }
190 171
191 G4double xs = 0.0; << 172 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
192 G4int intZ = G4lrint(Z); <<
193 if (intZ < 1 || intZ > maxZ) { <<
194 return xs; <<
195 } <<
196 173
197 G4PhysicsFreeVector* pv = dataCS[intZ]; << 174 // Select randomly one element in the current material
>> 175 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
198 176
199 // if element was not initialised << 177 // Sample the angle of the scattered photon
200 // do initialisation safely for MT mode <<
201 if (nullptr == pv) { <<
202 InitialiseForElement(nullptr, intZ); <<
203 pv = dataCS[intZ]; <<
204 if (nullptr == pv) { <<
205 return xs; <<
206 } <<
207 } <<
208 178
209 auto n = G4int(pv->GetVectorLength() - 1); << 179 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
210 G4double e = GammaEnergy / MeV; <<
211 if (e >= pv->Energy(n)) { <<
212 xs = (*pv)[n] / (e * e); <<
213 } <<
214 else if (e >= pv->Energy(0)) { <<
215 xs = pv->Value(e) / (e * e); <<
216 } <<
217 180
218 if (verboseLevel > 1) { << 181 G4double gReject,x,dataFormFactor;
219 G4cout << "****** DEBUG: tcs value for Z=" << 182 G4double randomFormFactor;
220 G4cout << " cs (Geant4 internal unit)=" < << 183 G4double cosTheta;
221 G4cout << " -> first E*E*cs value in CS << 184 G4double sinTheta;
222 G4cout << " -> last E*E*cs value in CS << 185 G4double fcostheta;
223 G4cout << "******************************* << 186
224 } << 187 do
225 return xs; << 188 {
>> 189 do
>> 190 {
>> 191 cosTheta = 2. * G4UniformRand() - 1.;
>> 192 fcostheta = ( 1. + cosTheta*cosTheta)/2.;
>> 193 } while (fcostheta < G4UniformRand());
>> 194
>> 195 G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
>> 196 x = sinThetaHalf / (wlPhoton/cm);
>> 197 if (x > 1.e+005)
>> 198 dataFormFactor = formFactorData->FindValue(x,Z-1);
>> 199 else
>> 200 dataFormFactor = formFactorData->FindValue(0.,Z-1);
>> 201 randomFormFactor = G4UniformRand() * Z * Z;
>> 202 sinTheta = std::sqrt(1. - cosTheta*cosTheta);
>> 203 gReject = dataFormFactor * dataFormFactor;
>> 204
>> 205 } while( gReject < randomFormFactor);
>> 206
>> 207 // Scattered photon angles. ( Z - axis along the parent photon)
>> 208 G4double phi = twopi * G4UniformRand() ;
>> 209 G4double dirX = sinTheta*std::cos(phi);
>> 210 G4double dirY = sinTheta*std::sin(phi);
>> 211 G4double dirZ = cosTheta;
>> 212
>> 213 // Update G4VParticleChange for the scattered photon
>> 214 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
>> 215 photonDirection1.rotateUz(photonDirection0);
>> 216 fParticleChange->ProposeMomentumDirection(photonDirection1);
>> 217
>> 218 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
226 } 219 }
227 220
228 //....oooOO0OOooo........oooOO0OOooo........oo 221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229 222
230 void G4LivermoreRayleighModel::SampleSecondari << 223 G4double G4LivermoreRayleighModel::GetMeanFreePath(const G4Track& track,
231 << 224 G4double, // previousStepSize
232 << 225 G4ForceCondition*)
233 << 226 {
234 { << 227 const G4DynamicParticle* photon = track.GetDynamicParticle();
235 if (verboseLevel > 1) { << 228 G4double energy = photon->GetKineticEnergy();
236 G4cout << "Calling SampleSecondaries() of << 229 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
237 } << 230 size_t materialIndex = couple->GetIndex();
238 G4double photonEnergy0 = aDynamicGamma->GetK << 231
239 << 232 G4double meanFreePath;
240 // Select randomly one element in the curren << 233 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
241 const G4ParticleDefinition* particle = aDyna << 234 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
242 const G4Element* elm = SelectRandomAtom(coup << 235 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
243 G4int Z = elm->GetZasInt(); << 236 return meanFreePath;
244 <<
245 // Sample the angle of the scattered photon <<
246 G4ThreeVector photonDirection = GetAngularDi <<
247 aDynamicGamma, photonEnergy0, Z, couple->G <<
248 fParticleChange->ProposeMomentumDirection(ph <<
249 } <<
250 <<
251 //....oooOO0OOooo........oooOO0OOooo........oo <<
252 <<
253 void G4LivermoreRayleighModel::InitialiseForEl <<
254 { <<
255 if (nullptr != dataCS[Z]) { <<
256 return; <<
257 } <<
258 G4AutoLock l(&LivermoreRayleighModelMutex); <<
259 if (nullptr == dataCS[Z]) { <<
260 ReadData(Z); <<
261 } <<
262 l.unlock(); <<
263 } 237 }
264 238
265 //....oooOO0OOooo........oooOO0OOooo........oo <<
266 239