| 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 // Author: Sebastien Incerti
27 // 31 March 2012 27 // 31 March 2012
28 // on base of G4LivermoreRayleighModel 28 // on base of G4LivermoreRayleighModel
29 // 29 //
30 30
31 #include "G4LivermoreRayleighModel.hh" 31 #include "G4LivermoreRayleighModel.hh"
32 <<
33 #include "G4AutoLock.hh" <<
34 #include "G4EmParameters.hh" <<
35 #include "G4RayleighAngularGenerator.hh" <<
36 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
>> 33 #include "G4RayleighAngularGenerator.hh"
37 34
38 //....oooOO0OOooo........oooOO0OOooo........oo << 35 using namespace std;
39 <<
40 namespace <<
41 { <<
42 G4Mutex LivermoreRayleighModelMutex = G4MUTEX_ <<
43 } <<
44 <<
45 G4PhysicsFreeVector* G4LivermoreRayleighModel: <<
46 G4String G4LivermoreRayleighModel::gDataDirect <<
47 36
48 //....oooOO0OOooo........oooOO0OOooo........oo 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 39
50 G4LivermoreRayleighModel::G4LivermoreRayleighM << 40 G4int G4LivermoreRayleighModel::maxZ = 100;
51 { << 41 G4LPhysicsFreeVector* G4LivermoreRayleighModel::dataCS[] = {0};
52 fParticleChange = nullptr; <<
53 lowEnergyLimit = 10 * CLHEP::eV; <<
54 42
>> 43 G4LivermoreRayleighModel::G4LivermoreRayleighModel()
>> 44 :G4VEmModel("LivermoreRayleigh"),isInitialised(false)
>> 45 {
>> 46 fParticleChange = 0;
>> 47 lowEnergyLimit = 10 * eV;
>> 48
55 SetAngularDistribution(new G4RayleighAngular 49 SetAngularDistribution(new G4RayleighAngularGenerator());
56 << 50
57 verboseLevel = 0; << 51 verboseLevel= 0;
58 // Verbosity scale for debugging purposes: 52 // Verbosity scale for debugging purposes:
59 // 0 = nothing << 53 // 0 = nothing
60 // 1 = calculation of cross sections, file o 54 // 1 = calculation of cross sections, file openings...
61 // 2 = entering in methods 55 // 2 = entering in methods
62 56
63 if (verboseLevel > 0) { << 57 if(verboseLevel > 0)
>> 58 {
64 G4cout << "G4LivermoreRayleighModel is con 59 G4cout << "G4LivermoreRayleighModel is constructed " << G4endl;
65 } 60 }
>> 61
66 } 62 }
67 63
68 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 65
70 G4LivermoreRayleighModel::~G4LivermoreRayleigh 66 G4LivermoreRayleighModel::~G4LivermoreRayleighModel()
71 { 67 {
72 if (IsMaster()) { << 68 if(IsMaster()) {
73 for (G4int i = 0; i <= maxZ; ++i) { << 69 for(G4int i=0; i<maxZ; ++i) {
74 if (nullptr != dataCS[i]) { << 70 if(dataCS[i]) {
75 delete dataCS[i]; << 71 delete dataCS[i];
76 dataCS[i] = nullptr; << 72 dataCS[i] = 0;
77 } 73 }
78 } 74 }
79 } 75 }
80 } 76 }
81 77
82 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 79
84 void G4LivermoreRayleighModel::Initialise(cons 80 void G4LivermoreRayleighModel::Initialise(const G4ParticleDefinition* particle,
85 cons << 81 const G4DataVector& cuts)
86 { 82 {
87 if (verboseLevel > 1) { << 83 if (verboseLevel > 1)
>> 84 {
88 G4cout << "Calling Initialise() of G4Liver 85 G4cout << "Calling Initialise() of G4LivermoreRayleighModel." << G4endl
89 << "Energy range: " << LowEnergyLim << 86 << "Energy range: "
90 << " GeV" << G4endl; << 87 << LowEnergyLimit() / eV << " eV - "
>> 88 << HighEnergyLimit() / GeV << " GeV"
>> 89 << G4endl;
91 } 90 }
92 91
93 if (IsMaster()) { << 92 if(IsMaster()) {
>> 93
94 // Initialise element selector 94 // Initialise element selector
95 InitialiseElementSelectors(particle, cuts) 95 InitialiseElementSelectors(particle, cuts);
96 96
97 // Access to elements 97 // Access to elements
98 const G4ElementTable* elemTable = G4Elemen << 98 char* path = std::getenv("G4LEDATA");
99 std::size_t numElems = (*elemTable).size() << 99 G4ProductionCutsTable* theCoupleTable =
100 for (std::size_t ie = 0; ie < numElems; ++ << 100 G4ProductionCutsTable::GetProductionCutsTable();
101 const G4Element* elem = (*elemTable)[ie] << 101 G4int numOfCouples = theCoupleTable->GetTableSize();
102 const G4int Z = std::min(maxZ, elem->Get << 102
103 if (dataCS[Z] == nullptr) { << 103 for(G4int i=0; i<numOfCouples; ++i)
104 ReadData(Z); << 104 {
>> 105 const G4MaterialCutsCouple* couple =
>> 106 theCoupleTable->GetMaterialCutsCouple(i);
>> 107 const G4Material* material = couple->GetMaterial();
>> 108 const G4ElementVector* theElementVector = material->GetElementVector();
>> 109 G4int nelm = material->GetNumberOfElements();
>> 110
>> 111 for (G4int j=0; j<nelm; ++j)
>> 112 {
>> 113 G4int Z = (*theElementVector)[j]->GetZasInt();
>> 114 if(Z < 1) { Z = 1; }
>> 115 else if(Z > maxZ) { Z = maxZ; }
>> 116 if( (!dataCS[Z]) ) { ReadData(Z, path); }
>> 117 }
105 } 118 }
106 } <<
107 } <<
108 if (isInitialised) { <<
109 return; <<
110 } 119 }
>> 120
>> 121 if(isInitialised) { return; }
111 fParticleChange = GetParticleChangeForGamma( 122 fParticleChange = GetParticleChangeForGamma();
112 isInitialised = true; 123 isInitialised = true;
>> 124
113 } 125 }
114 126
115 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 128
117 void G4LivermoreRayleighModel::InitialiseLocal << 129 void G4LivermoreRayleighModel::InitialiseLocal(const G4ParticleDefinition*,
>> 130 G4VEmModel* masterModel)
118 { 131 {
119 SetElementSelectors(masterModel->GetElementS 132 SetElementSelectors(masterModel->GetElementSelectors());
120 } 133 }
121 134
122 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 136
124 const G4String& G4LivermoreRayleighModel::Find << 137 void G4LivermoreRayleighModel::ReadData(size_t Z, const char* path)
125 { 138 {
126 // no check in this method - environment var << 139 if (verboseLevel > 1)
127 if (gDataDirectory.empty()) { << 140 {
128 auto param = G4EmParameters::Instance(); << 141 G4cout << "Calling ReadData() of G4LivermoreRayleighModel"
129 std::ostringstream ost; << 142 << G4endl;
130 if (param->LivermoreDataDir() == "livermor << 143 }
131 ost << param->GetDirLEDATA() << "/liverm << 144
>> 145 if(dataCS[Z]) { return; }
>> 146
>> 147 const char* datadir = path;
>> 148
>> 149 if(!datadir)
>> 150 {
>> 151 datadir = std::getenv("G4LEDATA");
>> 152 if(!datadir)
>> 153 {
>> 154 G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
>> 155 FatalException,
>> 156 "Environment variable G4LEDATA not defined");
>> 157 return;
132 } 158 }
133 else { <<
134 ost << param->GetDirLEDATA() << "/epics2 <<
135 } <<
136 gDataDirectory = ost.str(); <<
137 } 159 }
138 return gDataDirectory; <<
139 } <<
140 <<
141 //....oooOO0OOooo........oooOO0OOooo........oo <<
142 <<
143 void G4LivermoreRayleighModel::ReadData(const <<
144 { <<
145 if (verboseLevel > 1) { <<
146 G4cout << "Calling ReadData() of G4Livermo <<
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 160
>> 161 //
>> 162
>> 163 dataCS[Z] = new G4LPhysicsFreeVector();
>> 164
>> 165 // Activation of spline interpolation
>> 166 //dataCS[Z] ->SetSpline(true);
>> 167
156 std::ostringstream ostCS; 168 std::ostringstream ostCS;
157 ostCS << FindDirectoryPath() << "re-cs-" << << 169 ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
158 <<
159 std::ifstream finCS(ostCS.str().c_str()); 170 std::ifstream finCS(ostCS.str().c_str());
160 << 171
161 if (!finCS.is_open()) { << 172 if( !finCS .is_open() )
>> 173 {
162 G4ExceptionDescription ed; 174 G4ExceptionDescription ed;
163 ed << "G4LivermoreRayleighModel data file << 175 ed << "G4LivermoreRayleighModel data file <" << ostCS.str().c_str()
164 << G4endl; << 176 << "> is not opened!" << G4endl;
165 G4Exception("G4LivermoreRayleighModel::Rea << 177 G4Exception("G4LivermoreRayleighModel::ReadData()","em0003",FatalException,
166 "G4LEDATA version should be G4 << 178 ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
167 return; 179 return;
168 } << 180 }
169 else { << 181 else
170 if (verboseLevel > 3) { << 182 {
171 G4cout << "File " << ostCS.str() << " is << 183 if(verboseLevel > 3) {
>> 184 G4cout << "File " << ostCS.str()
>> 185 << " is opened by G4LivermoreRayleighModel" << G4endl;
172 } 186 }
173 dataCS[Z]->Retrieve(finCS, true); 187 dataCS[Z]->Retrieve(finCS, true);
174 } << 188 }
175 } 189 }
176 190
177 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 192
179 G4double G4LivermoreRayleighModel::ComputeCros << 193 G4double G4LivermoreRayleighModel::ComputeCrossSectionPerAtom(
180 << 194 const G4ParticleDefinition*,
181 << 195 G4double GammaEnergy,
>> 196 G4double Z, G4double,
>> 197 G4double, G4double)
182 { 198 {
183 if (verboseLevel > 1) { << 199 if (verboseLevel > 1)
184 G4cout << "G4LivermoreRayleighModel::Compu << 200 {
185 } << 201 G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()"
186 << 202 << G4endl;
187 if (GammaEnergy < lowEnergyLimit) { <<
188 return 0.0; <<
189 } 203 }
190 204
>> 205 if(GammaEnergy < lowEnergyLimit) { return 0.0; }
>> 206
191 G4double xs = 0.0; 207 G4double xs = 0.0;
>> 208
192 G4int intZ = G4lrint(Z); 209 G4int intZ = G4lrint(Z);
193 if (intZ < 1 || intZ > maxZ) { <<
194 return xs; <<
195 } <<
196 210
197 G4PhysicsFreeVector* pv = dataCS[intZ]; << 211 if(intZ < 1 || intZ > maxZ) { return xs; }
>> 212
>> 213 G4LPhysicsFreeVector* pv = dataCS[intZ];
198 214
199 // if element was not initialised 215 // if element was not initialised
200 // do initialisation safely for MT mode 216 // do initialisation safely for MT mode
201 if (nullptr == pv) { << 217 if(!pv) {
202 InitialiseForElement(nullptr, intZ); << 218 InitialiseForElement(0, intZ);
203 pv = dataCS[intZ]; 219 pv = dataCS[intZ];
204 if (nullptr == pv) { << 220 if(!pv) { return xs; }
205 return xs; <<
206 } <<
207 } 221 }
208 222
209 auto n = G4int(pv->GetVectorLength() - 1); << 223 G4int n = pv->GetVectorLength() - 1;
210 G4double e = GammaEnergy / MeV; << 224 G4double e = GammaEnergy/MeV;
211 if (e >= pv->Energy(n)) { << 225 if(e >= pv->Energy(n)) {
212 xs = (*pv)[n] / (e * e); << 226 xs = (*pv)[n]/(e*e);
213 } << 227 } else if(e >= pv->Energy(0)) {
214 else if (e >= pv->Energy(0)) { << 228 xs = pv->Value(e)/(e*e);
215 xs = pv->Value(e) / (e * e); << 229 }
216 } << 230
217 << 231 if(verboseLevel > 0)
218 if (verboseLevel > 1) { << 232 {
219 G4cout << "****** DEBUG: tcs value for Z=" << 233 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
220 G4cout << " cs (Geant4 internal unit)=" < << 234 << e << G4endl;
221 G4cout << " -> first E*E*cs value in CS << 235 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
222 G4cout << " -> last E*E*cs value in CS << 236 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
223 G4cout << "******************************* << 237 << G4endl;
>> 238 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
>> 239 << G4endl;
>> 240 G4cout << "*********************************************************"
>> 241 << G4endl;
224 } 242 }
225 return xs; 243 return xs;
226 } 244 }
227 245
228 //....oooOO0OOooo........oooOO0OOooo........oo 246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229 247
230 void G4LivermoreRayleighModel::SampleSecondari << 248 void G4LivermoreRayleighModel::SampleSecondaries(
231 << 249 std::vector<G4DynamicParticle*>*,
232 << 250 const G4MaterialCutsCouple* couple,
233 << 251 const G4DynamicParticle* aDynamicGamma,
>> 252 G4double, G4double)
234 { 253 {
235 if (verboseLevel > 1) { 254 if (verboseLevel > 1) {
236 G4cout << "Calling SampleSecondaries() of << 255 G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel"
>> 256 << G4endl;
237 } 257 }
238 G4double photonEnergy0 = aDynamicGamma->GetK 258 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
239 259
>> 260 // absorption of low-energy gamma
>> 261 /*
>> 262 if (photonEnergy0 <= lowEnergyLimit)
>> 263 {
>> 264 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 265 fParticleChange->SetProposedKineticEnergy(0.);
>> 266 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
>> 267 return ;
>> 268 }
>> 269 */
240 // Select randomly one element in the curren 270 // Select randomly one element in the current material
241 const G4ParticleDefinition* particle = aDyna << 271 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
242 const G4Element* elm = SelectRandomAtom(coup << 272 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
243 G4int Z = elm->GetZasInt(); << 273 G4int Z = G4lrint(elm->GetZ());
244 274
245 // Sample the angle of the scattered photon 275 // Sample the angle of the scattered photon
246 G4ThreeVector photonDirection = GetAngularDi << 276
247 aDynamicGamma, photonEnergy0, Z, couple->G << 277 G4ThreeVector photonDirection =
>> 278 GetAngularDistribution()->SampleDirection(aDynamicGamma,
>> 279 photonEnergy0,
>> 280 Z, couple->GetMaterial());
248 fParticleChange->ProposeMomentumDirection(ph 281 fParticleChange->ProposeMomentumDirection(photonDirection);
249 } 282 }
250 283
251 //....oooOO0OOooo........oooOO0OOooo........oo 284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252 285
253 void G4LivermoreRayleighModel::InitialiseForEl << 286 #include "G4AutoLock.hh"
>> 287 namespace { G4Mutex LivermoreRayleighModelMutex = G4MUTEX_INITIALIZER; }
>> 288
>> 289 void
>> 290 G4LivermoreRayleighModel::InitialiseForElement(const G4ParticleDefinition*,
>> 291 G4int Z)
254 { 292 {
255 if (nullptr != dataCS[Z]) { <<
256 return; <<
257 } <<
258 G4AutoLock l(&LivermoreRayleighModelMutex); 293 G4AutoLock l(&LivermoreRayleighModelMutex);
259 if (nullptr == dataCS[Z]) { << 294 // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= "
260 ReadData(Z); << 295 // << Z << G4endl;
261 } << 296 if(!dataCS[Z]) { ReadData(Z); }
262 l.unlock(); 297 l.unlock();
263 } 298 }
264 299
265 //....oooOO0OOooo........oooOO0OOooo........oo 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
266 301