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