| 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 // $Id$
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 08 Mar 2012 L Pandola First complete i 32 // 08 Mar 2012 L Pandola First complete implementation
32 // 33 //
33 34
34 #include "G4PenelopeIonisationXSHandler.hh" 35 #include "G4PenelopeIonisationXSHandler.hh"
35 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
37 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh"
38 #include "G4Electron.hh" 39 #include "G4Electron.hh"
39 #include "G4Positron.hh" 40 #include "G4Positron.hh"
40 #include "G4PenelopeOscillatorManager.hh" 41 #include "G4PenelopeOscillatorManager.hh"
41 #include "G4PenelopeOscillator.hh" 42 #include "G4PenelopeOscillator.hh"
42 #include "G4PenelopeCrossSection.hh" 43 #include "G4PenelopeCrossSection.hh"
43 #include "G4PhysicsFreeVector.hh" 44 #include "G4PhysicsFreeVector.hh"
44 #include "G4PhysicsLogVector.hh" 45 #include "G4PhysicsLogVector.hh"
45 46
46 G4PenelopeIonisationXSHandler::G4PenelopeIonis 47 G4PenelopeIonisationXSHandler::G4PenelopeIonisationXSHandler(size_t nb)
47 :fXSTableElectron(nullptr),fXSTablePositron( << 48 :XSTableElectron(0),XSTablePositron(0),
48 fDeltaTable(nullptr),fEnergyGrid(nullptr) << 49 theDeltaTable(0),energyGrid(0)
49 { 50 {
50 fNBins = nb; << 51 nBins = nb;
51 G4double LowEnergyLimit = 100.0*eV; 52 G4double LowEnergyLimit = 100.0*eV;
52 G4double HighEnergyLimit = 100.0*GeV; 53 G4double HighEnergyLimit = 100.0*GeV;
53 fOscManager = G4PenelopeOscillatorManager::G << 54 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
54 fXSTableElectron = new << 55 XSTableElectron = new
55 std::map< std::pair<const G4Material*,G4do 56 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
56 fXSTablePositron = new << 57 XSTablePositron = new
57 std::map< std::pair<const G4Material*,G4do 58 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
58 59
59 fDeltaTable = new std::map<const G4Material* << 60 theDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
60 fEnergyGrid = new G4PhysicsLogVector(LowEner << 61 energyGrid = new G4PhysicsLogVector(LowEnergyLimit,
61 HighEnergyLimit, 62 HighEnergyLimit,
62 fNBins-1); //one hidden bin is a << 63 nBins-1); //one hidden bin is added
63 fVerboseLevel = 0; << 64
>> 65 verboseLevel = 0;
64 } 66 }
65 67
66 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 69
68 G4PenelopeIonisationXSHandler::~G4PenelopeIoni 70 G4PenelopeIonisationXSHandler::~G4PenelopeIonisationXSHandler()
69 { 71 {
70 if (fXSTableElectron) << 72 if (XSTableElectron)
71 { 73 {
72 for (auto& item : (*fXSTableElectron)) << 74 for (auto& item : (*XSTableElectron))
73 { 75 {
74 //G4PenelopeCrossSection* tab = i->second; 76 //G4PenelopeCrossSection* tab = i->second;
75 delete item.second; 77 delete item.second;
76 } 78 }
77 delete fXSTableElectron; << 79 delete XSTableElectron;
78 fXSTableElectron = nullptr; << 80 XSTableElectron = nullptr;
79 } 81 }
80 82
81 if (fXSTablePositron) << 83 if (XSTablePositron)
82 { 84 {
83 for (auto& item : (*fXSTablePositron)) << 85 for (auto& item : (*XSTablePositron))
84 { 86 {
85 //G4PenelopeCrossSection* tab = i->second; 87 //G4PenelopeCrossSection* tab = i->second;
86 delete item.second; 88 delete item.second;
87 } 89 }
88 delete fXSTablePositron; << 90 delete XSTablePositron;
89 fXSTablePositron = nullptr; << 91 XSTablePositron = nullptr;
90 } 92 }
91 if (fDeltaTable) << 93 if (theDeltaTable)
92 { 94 {
93 for (auto& item : (*fDeltaTable)) << 95 for (auto& item : (*theDeltaTable))
94 delete item.second; 96 delete item.second;
95 delete fDeltaTable; << 97 delete theDeltaTable;
96 fDeltaTable = nullptr; << 98 theDeltaTable = nullptr;
97 } 99 }
98 if (fEnergyGrid) << 100 if (energyGrid)
99 delete fEnergyGrid; << 101 delete energyGrid;
100 102
101 if (fVerboseLevel > 2) << 103 if (verboseLevel > 2)
102 G4cout << "G4PenelopeIonisationXSHandler. 104 G4cout << "G4PenelopeIonisationXSHandler. Tables have been cleared"
103 << G4endl; 105 << G4endl;
104 } 106 }
105 107
106 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 109
108 const G4PenelopeCrossSection* 110 const G4PenelopeCrossSection*
109 G4PenelopeIonisationXSHandler::GetCrossSection 111 G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
110 const G4Material* mat, 112 const G4Material* mat,
111 const G4double cut) const 113 const G4double cut) const
112 { 114 {
113 if (part != G4Electron::Electron() && part ! 115 if (part != G4Electron::Electron() && part != G4Positron::Positron())
114 { 116 {
115 G4ExceptionDescription ed; 117 G4ExceptionDescription ed;
116 ed << "Invalid particle: " << part->GetP 118 ed << "Invalid particle: " << part->GetParticleName() << G4endl;
117 G4Exception("G4PenelopeIonisationXSHandl 119 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
118 "em0001",FatalException,ed); 120 "em0001",FatalException,ed);
119 return nullptr; 121 return nullptr;
120 } 122 }
121 123
122 if (part == G4Electron::Electron()) 124 if (part == G4Electron::Electron())
123 { 125 {
124 if (!fXSTableElectron) << 126 if (!XSTableElectron)
125 { 127 {
126 G4Exception("G4PenelopeIonisationXSHandler 128 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
127 "em0028",FatalException, 129 "em0028",FatalException,
128 "The Cross Section Table for e- was 130 "The Cross Section Table for e- was not initialized correctly!");
129 return nullptr; 131 return nullptr;
130 } 132 }
131 std::pair<const G4Material*,G4double> th 133 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
132 if (fXSTableElectron->count(theKey)) //t << 134 if (XSTableElectron->count(theKey)) //table already built
133 return fXSTableElectron->find(theKey)->secon << 135 return XSTableElectron->find(theKey)->second;
134 else 136 else
135 return nullptr; 137 return nullptr;
136 } 138 }
137 139
138 if (part == G4Positron::Positron()) 140 if (part == G4Positron::Positron())
139 { 141 {
140 if (!fXSTablePositron) << 142 if (!XSTablePositron)
141 { 143 {
142 G4Exception("G4PenelopeIonisationXSHandler 144 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
143 "em0028",FatalException, 145 "em0028",FatalException,
144 "The Cross Section Table for e+ was 146 "The Cross Section Table for e+ was not initialized correctly!");
145 return nullptr; 147 return nullptr;
146 } 148 }
147 std::pair<const G4Material*,G4double> th 149 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
148 if (fXSTablePositron->count(theKey)) //t << 150 if (XSTablePositron->count(theKey)) //table already built
149 return fXSTablePositron->find(theKey)->secon << 151 return XSTablePositron->find(theKey)->second;
150 else 152 else
151 return nullptr; 153 return nullptr;
152 } 154 }
153 return nullptr; 155 return nullptr;
154 } 156 }
155 157
156 //....oooOO0OOooo........oooOO0OOooo........oo 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 159
158 void G4PenelopeIonisationXSHandler::BuildXSTab 160 void G4PenelopeIonisationXSHandler::BuildXSTable(const G4Material* mat,G4double cut,
159 const G4ParticleDefinition* part, 161 const G4ParticleDefinition* part,
160 G4bool isMaster) 162 G4bool isMaster)
161 { 163 {
162 //Just to check 164 //Just to check
163 if (!isMaster) 165 if (!isMaster)
164 G4Exception("G4PenelopeIonisationXSHandler 166 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
165 "em0100",FatalException,"Worke 167 "em0100",FatalException,"Worker thread in this method");
166 168
167 // 169 //
168 //This method fills the G4PenelopeCrossSecti 170 //This method fills the G4PenelopeCrossSection containers for electrons or positrons
169 //and for the given material/cut couple. The 171 //and for the given material/cut couple. The calculation is done as sum over the
170 //individual shells. 172 //individual shells.
171 //Equivalent of subroutines EINaT and PINaT 173 //Equivalent of subroutines EINaT and PINaT of Penelope
172 // 174 //
173 if (fVerboseLevel > 2) << 175 if (verboseLevel > 2)
174 { 176 {
175 G4cout << "G4PenelopeIonisationXSHandler 177 G4cout << "G4PenelopeIonisationXSHandler: going to build cross section table " << G4endl;
176 G4cout << "for " << part->GetParticleNam 178 G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl;
177 G4cout << "Cut= " << cut/keV << " keV" < 179 G4cout << "Cut= " << cut/keV << " keV" << G4endl;
178 } 180 }
179 181
180 std::pair<const G4Material*,G4double> theKey 182 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
181 //Check if the table already exists 183 //Check if the table already exists
182 if (part == G4Electron::Electron()) 184 if (part == G4Electron::Electron())
183 { 185 {
184 if (fXSTableElectron->count(theKey)) //t << 186 if (XSTableElectron->count(theKey)) //table already built
185 return; 187 return;
186 } 188 }
187 if (part == G4Positron::Positron()) 189 if (part == G4Positron::Positron())
188 { 190 {
189 if (fXSTablePositron->count(theKey)) //t << 191 if (XSTablePositron->count(theKey)) //table already built
190 return; 192 return;
191 } 193 }
192 194
193 //check if the material has been built 195 //check if the material has been built
194 if (!(fDeltaTable->count(mat))) << 196 if (!(theDeltaTable->count(mat)))
195 BuildDeltaTable(mat); 197 BuildDeltaTable(mat);
196 198
197 199
198 //Tables have been already created (checked 200 //Tables have been already created (checked by GetCrossSectionTableForCouple)
199 G4PenelopeOscillatorTable* theTable = fOscMa << 201 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
200 size_t numberOfOscillators = theTable->size( 202 size_t numberOfOscillators = theTable->size();
201 203
202 if (fEnergyGrid->GetVectorLength() != fNBins << 204 if (energyGrid->GetVectorLength() != nBins)
203 { 205 {
204 G4ExceptionDescription ed; 206 G4ExceptionDescription ed;
205 ed << "Energy Grid looks not initialized 207 ed << "Energy Grid looks not initialized" << G4endl;
206 ed << fNBins << " " << fEnergyGrid->GetV << 208 ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
207 G4Exception("G4PenelopeIonisationXSHandl 209 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
208 "em2030",FatalException,ed); 210 "em2030",FatalException,ed);
209 } 211 }
210 212
211 G4PenelopeCrossSection* XSEntry = new G4Pene << 213 G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins,numberOfOscillators);
212 214
213 //loop on the energy grid 215 //loop on the energy grid
214 for (size_t bin=0;bin<fNBins;bin++) << 216 for (size_t bin=0;bin<nBins;bin++)
215 { 217 {
216 G4double energy = fEnergyGrid->GetLowEd << 218 G4double energy = energyGrid->GetLowEdgeEnergy(bin);
217 G4double XH0=0, XH1=0, XH2=0; 219 G4double XH0=0, XH1=0, XH2=0;
218 G4double XS0=0, XS1=0, XS2=0; 220 G4double XS0=0, XS1=0, XS2=0;
219 221
220 //oscillator loop 222 //oscillator loop
221 for (size_t iosc=0;iosc<numberOfOscilla 223 for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
222 { 224 {
223 G4DataVector* tempStorage = nullptr; << 225 G4DataVector* tempStorage = 0;
224 226
225 G4PenelopeOscillator* theOsc = (*theTable 227 G4PenelopeOscillator* theOsc = (*theTable)[iosc];
226 G4double delta = GetDensityCorrection(mat 228 G4double delta = GetDensityCorrection(mat,energy);
227 if (part == G4Electron::Electron()) 229 if (part == G4Electron::Electron())
228 tempStorage = ComputeShellCrossSections 230 tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta);
229 else if (part == G4Positron::Positron()) 231 else if (part == G4Positron::Positron())
230 tempStorage = ComputeShellCrossSections 232 tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta);
231 //check results are all right 233 //check results are all right
232 if (!tempStorage) 234 if (!tempStorage)
233 { 235 {
234 G4ExceptionDescription ed; 236 G4ExceptionDescription ed;
235 ed << "Problem in calculating the she 237 ed << "Problem in calculating the shell XS for shell # "
236 << iosc << G4endl; 238 << iosc << G4endl;
237 G4Exception("G4PenelopeIonisationXSHa 239 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
238 "em2031",FatalException,ed); 240 "em2031",FatalException,ed);
239 delete XSEntry; 241 delete XSEntry;
240 return; 242 return;
241 } 243 }
242 if (tempStorage->size() != 6) 244 if (tempStorage->size() != 6)
243 { 245 {
244 G4ExceptionDescription ed; 246 G4ExceptionDescription ed;
245 ed << "Problem in calculating the she 247 ed << "Problem in calculating the shell XS " << G4endl;
246 ed << "Result has dimension " << temp 248 ed << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl;
247 G4Exception("G4PenelopeIonisationXSHa 249 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
248 "em2031",FatalException,ed); 250 "em2031",FatalException,ed);
249 } 251 }
250 G4double stre = theOsc->GetOscillatorStre 252 G4double stre = theOsc->GetOscillatorStrength();
251 253
252 XH0 += stre*(*tempStorage)[0]; 254 XH0 += stre*(*tempStorage)[0];
253 XH1 += stre*(*tempStorage)[1]; 255 XH1 += stre*(*tempStorage)[1];
254 XH2 += stre*(*tempStorage)[2]; 256 XH2 += stre*(*tempStorage)[2];
255 XS0 += stre*(*tempStorage)[3]; 257 XS0 += stre*(*tempStorage)[3];
256 XS1 += stre*(*tempStorage)[4]; 258 XS1 += stre*(*tempStorage)[4];
257 XS2 += stre*(*tempStorage)[5]; 259 XS2 += stre*(*tempStorage)[5];
258 XSEntry->AddShellCrossSectionPoint(bin,io 260 XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]);
259 if (tempStorage) 261 if (tempStorage)
260 { 262 {
261 delete tempStorage; 263 delete tempStorage;
262 tempStorage = nullptr; << 264 tempStorage = 0;
263 } 265 }
264 } 266 }
265 XSEntry->AddCrossSectionPoint(bin,energ 267 XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
266 } 268 }
267 //Do (only once) the final normalization 269 //Do (only once) the final normalization
268 XSEntry->NormalizeShellCrossSections(); 270 XSEntry->NormalizeShellCrossSections();
269 271
270 //Insert in the appropriate table 272 //Insert in the appropriate table
271 if (part == G4Electron::Electron()) 273 if (part == G4Electron::Electron())
272 fXSTableElectron->insert(std::make_pair(th << 274 XSTableElectron->insert(std::make_pair(theKey,XSEntry));
273 else if (part == G4Positron::Positron()) 275 else if (part == G4Positron::Positron())
274 fXSTablePositron->insert(std::make_pair(th << 276 XSTablePositron->insert(std::make_pair(theKey,XSEntry));
275 else 277 else
276 delete XSEntry; 278 delete XSEntry;
277 279
278 return; 280 return;
279 } 281 }
280 282
281 283
282 //....oooOO0OOooo........oooOO0OOooo........oo 284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
283 285
284 G4double G4PenelopeIonisationXSHandler::GetDen 286 G4double G4PenelopeIonisationXSHandler::GetDensityCorrection(const G4Material* mat,
285 const G4double energy) cons 287 const G4double energy) const
286 { 288 {
287 G4double result = 0; 289 G4double result = 0;
288 if (!fDeltaTable) << 290 if (!theDeltaTable)
289 { 291 {
290 G4Exception("G4PenelopeIonisationXSHandl 292 G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
291 "em2032",FatalException, 293 "em2032",FatalException,
292 "Delta Table not initialized. Was Initia 294 "Delta Table not initialized. Was Initialise() run?");
293 return 0; 295 return 0;
294 } 296 }
295 if (energy <= 0*eV) 297 if (energy <= 0*eV)
296 { 298 {
297 G4cout << "G4PenelopeIonisationXSHandler 299 G4cout << "G4PenelopeIonisationXSHandler::GetDensityCorrection()" << G4endl;
298 G4cout << "Invalid energy " << energy/eV 300 G4cout << "Invalid energy " << energy/eV << " eV " << G4endl;
299 return 0; 301 return 0;
300 } 302 }
301 G4double logene = G4Log(energy); << 303 G4double logene = std::log(energy);
302 304
303 if (fDeltaTable->count(mat)) << 305 if (theDeltaTable->count(mat))
304 { 306 {
305 const G4PhysicsFreeVector* vec = fDeltaT << 307 const G4PhysicsFreeVector* vec = theDeltaTable->find(mat)->second;
306 result = vec->Value(logene); //the table 308 result = vec->Value(logene); //the table has delta vs. ln(E)
307 } 309 }
308 else 310 else
309 { 311 {
310 G4ExceptionDescription ed; 312 G4ExceptionDescription ed;
311 ed << "Unable to build table for " << ma 313 ed << "Unable to build table for " << mat->GetName() << G4endl;
312 G4Exception("G4PenelopeIonisationXSHandl 314 G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
313 "em2033",FatalException,ed); 315 "em2033",FatalException,ed);
314 } 316 }
315 317
316 return result; 318 return result;
317 } 319 }
318 320
319 //....oooOO0OOooo........oooOO0OOooo........oo 321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320 322
321 void G4PenelopeIonisationXSHandler::BuildDelta 323 void G4PenelopeIonisationXSHandler::BuildDeltaTable(const G4Material* mat)
322 { 324 {
323 G4PenelopeOscillatorTable* theTable = fOscMa << 325 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
324 G4double plasmaSq = fOscManager->GetPlasmaEn << 326 G4double plasmaSq = oscManager->GetPlasmaEnergySquared(mat);
325 G4double totalZ = fOscManager->GetTotalZ(mat << 327 G4double totalZ = oscManager->GetTotalZ(mat);
326 size_t numberOfOscillators = theTable->size( 328 size_t numberOfOscillators = theTable->size();
327 329
328 if (fEnergyGrid->GetVectorLength() != fNBins << 330 if (energyGrid->GetVectorLength() != nBins)
329 { 331 {
330 G4ExceptionDescription ed; 332 G4ExceptionDescription ed;
331 ed << "Energy Grid for Delta table looks 333 ed << "Energy Grid for Delta table looks not initialized" << G4endl;
332 ed << fNBins << " " << fEnergyGrid->GetV << 334 ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
333 G4Exception("G4PenelopeIonisationXSHandl 335 G4Exception("G4PenelopeIonisationXSHandler::BuildDeltaTable()",
334 "em2030",FatalException,ed); 336 "em2030",FatalException,ed);
335 } 337 }
336 338
337 G4PhysicsFreeVector* theVector = new G4Physi << 339 G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(nBins);
338 340
339 //loop on the energy grid 341 //loop on the energy grid
340 for (size_t bin=0;bin<fNBins;bin++) << 342 for (size_t bin=0;bin<nBins;bin++)
341 { 343 {
342 G4double delta = 0.; 344 G4double delta = 0.;
343 G4double energy = fEnergyGrid->GetLowEdg << 345 G4double energy = energyGrid->GetLowEdgeEnergy(bin);
344 346
345 //Here calculate delta 347 //Here calculate delta
346 G4double gam = 1.0+(energy/electron_mass 348 G4double gam = 1.0+(energy/electron_mass_c2);
347 G4double gamSq = gam*gam; 349 G4double gamSq = gam*gam;
348 350
349 G4double TST = totalZ/(gamSq*plasmaSq); 351 G4double TST = totalZ/(gamSq*plasmaSq);
350 G4double wl2 = 0; 352 G4double wl2 = 0;
351 G4double fdel = 0; 353 G4double fdel = 0;
352 354
353 //loop on oscillators 355 //loop on oscillators
354 for (size_t i=0;i<numberOfOscillators;i+ 356 for (size_t i=0;i<numberOfOscillators;i++)
355 { 357 {
356 G4PenelopeOscillator* theOsc = (*theTable) 358 G4PenelopeOscillator* theOsc = (*theTable)[i];
357 G4double wri = theOsc->GetResonanceEnergy( 359 G4double wri = theOsc->GetResonanceEnergy();
358 fdel += theOsc->GetOscillatorStrength()/(w 360 fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
359 } 361 }
360 if (fdel >= TST) //if fdel < TST, delta 362 if (fdel >= TST) //if fdel < TST, delta = 0
361 { 363 {
362 //get last oscillator 364 //get last oscillator
363 G4PenelopeOscillator* theOsc = (*theTable) 365 G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1];
364 wl2 = theOsc->GetResonanceEnergy()*theOsc- 366 wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy();
365 367
366 //First iteration 368 //First iteration
367 G4bool loopAgain = false; 369 G4bool loopAgain = false;
368 do 370 do
369 { 371 {
370 loopAgain = false; 372 loopAgain = false;
371 wl2 += wl2; 373 wl2 += wl2;
372 fdel = 0.; 374 fdel = 0.;
373 for (size_t i=0;i<numberOfOscillators; 375 for (size_t i=0;i<numberOfOscillators;i++)
374 { 376 {
375 G4PenelopeOscillator* theOscLocal1 = (*t 377 G4PenelopeOscillator* theOscLocal1 = (*theTable)[i];
376 G4double wri = theOscLocal1->GetResonanc 378 G4double wri = theOscLocal1->GetResonanceEnergy();
377 fdel += theOscLocal1->GetOscillatorStren 379 fdel += theOscLocal1->GetOscillatorStrength()/(wri*wri+wl2);
378 } 380 }
379 if (fdel > TST) 381 if (fdel > TST)
380 loopAgain = true; 382 loopAgain = true;
381 }while(loopAgain); 383 }while(loopAgain);
382 384
383 G4double wl2l = 0; 385 G4double wl2l = 0;
384 G4double wl2u = wl2; 386 G4double wl2u = wl2;
385 //second iteration 387 //second iteration
386 do 388 do
387 { 389 {
388 loopAgain = false; 390 loopAgain = false;
389 wl2 = 0.5*(wl2l+wl2u); 391 wl2 = 0.5*(wl2l+wl2u);
390 fdel = 0; 392 fdel = 0;
391 for (size_t i=0;i<numberOfOscillators; 393 for (size_t i=0;i<numberOfOscillators;i++)
392 { 394 {
393 G4PenelopeOscillator* theOscLocal2 = (*t 395 G4PenelopeOscillator* theOscLocal2 = (*theTable)[i];
394 G4double wri = theOscLocal2->GetResonanc 396 G4double wri = theOscLocal2->GetResonanceEnergy();
395 fdel += theOscLocal2->GetOscillatorStren 397 fdel += theOscLocal2->GetOscillatorStrength()/(wri*wri+wl2);
396 } 398 }
397 if (fdel > TST) 399 if (fdel > TST)
398 wl2l = wl2; 400 wl2l = wl2;
399 else 401 else
400 wl2u = wl2; 402 wl2u = wl2;
401 if ((wl2u-wl2l)>1e-12*wl2) 403 if ((wl2u-wl2l)>1e-12*wl2)
402 loopAgain = true; 404 loopAgain = true;
403 }while(loopAgain); 405 }while(loopAgain);
404 406
405 //Eventually get density correction 407 //Eventually get density correction
406 delta = 0.; 408 delta = 0.;
407 for (size_t i=0;i<numberOfOscillators;i++) 409 for (size_t i=0;i<numberOfOscillators;i++)
408 { 410 {
409 G4PenelopeOscillator* theOscLocal3 = ( 411 G4PenelopeOscillator* theOscLocal3 = (*theTable)[i];
410 G4double wri = theOscLocal3->GetResona 412 G4double wri = theOscLocal3->GetResonanceEnergy();
411 delta += theOscLocal3->GetOscillatorSt 413 delta += theOscLocal3->GetOscillatorStrength()*
412 G4Log(1.0+(wl2/(wri*wri))); << 414 std::log(1.0+(wl2/(wri*wri)));
413 } 415 }
414 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq 416 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq);
415 } 417 }
416 energy = std::max(1e-9*eV,energy); //pre 418 energy = std::max(1e-9*eV,energy); //prevents log(0)
417 theVector->PutValue(bin,G4Log(energy),de << 419 theVector->PutValue(bin,std::log(energy),delta);
418 } 420 }
419 fDeltaTable->insert(std::make_pair(mat,theVe << 421 theDeltaTable->insert(std::make_pair(mat,theVector));
420 return; 422 return;
421 } 423 }
422 424
423 //....oooOO0OOooo........oooOO0OOooo........oo 425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
424 G4DataVector* G4PenelopeIonisationXSHandler::C 426 G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsElectron(G4PenelopeOscillator* theOsc,
425 G4double energy, 427 G4double energy,
426 G4double cut, 428 G4double cut,
427 G4double delta) 429 G4double delta)
428 { 430 {
429 // 431 //
430 //This method calculates the hard and soft c 432 //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for
431 //the given oscillator/cut and at the given 433 //the given oscillator/cut and at the given energy.
432 //It returns a G4DataVector* with 6 entries 434 //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
433 //Equivalent of subroutines EINaT1 of Penelo 435 //Equivalent of subroutines EINaT1 of Penelope
434 // 436 //
435 // Results are _per target electron_ 437 // Results are _per target electron_
436 // 438 //
437 G4DataVector* result = new G4DataVector(); 439 G4DataVector* result = new G4DataVector();
438 for (size_t i=0;i<6;i++) 440 for (size_t i=0;i<6;i++)
439 result->push_back(0.); 441 result->push_back(0.);
440 G4double ionEnergy = theOsc->GetIonisationEn 442 G4double ionEnergy = theOsc->GetIonisationEnergy();
441 443
442 //return a set of zero's if the energy it to 444 //return a set of zero's if the energy it too low to excite the current oscillator
443 if (energy < ionEnergy) 445 if (energy < ionEnergy)
444 return result; 446 return result;
445 447
446 G4double H0=0.,H1=0.,H2=0.; 448 G4double H0=0.,H1=0.,H2=0.;
447 G4double S0=0.,S1=0.,S2=0.; 449 G4double S0=0.,S1=0.,S2=0.;
448 450
449 //Define useful constants to be used in the 451 //Define useful constants to be used in the calculation
450 G4double gamma = 1.0+energy/electron_mass_c2 452 G4double gamma = 1.0+energy/electron_mass_c2;
451 G4double gammaSq = gamma*gamma; 453 G4double gammaSq = gamma*gamma;
452 G4double beta = (gammaSq-1.0)/gammaSq; 454 G4double beta = (gammaSq-1.0)/gammaSq;
453 G4double pielr2 = pi*classic_electr_radius*c 455 G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
454 G4double constant = pielr2*2.0*electron_mass 456 G4double constant = pielr2*2.0*electron_mass_c2/beta;
455 G4double XHDT0 = G4Log(gammaSq)-beta; << 457 G4double XHDT0 = std::log(gammaSq)-beta;
456 458
457 G4double cpSq = energy*(energy+2.0*electron_ 459 G4double cpSq = energy*(energy+2.0*electron_mass_c2);
458 G4double cp = std::sqrt(cpSq); 460 G4double cp = std::sqrt(cpSq);
459 G4double amol = (energy/(energy+electron_mas 461 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
460 462
461 // 463 //
462 // Distant interactions 464 // Distant interactions
463 // 465 //
464 G4double resEne = theOsc->GetResonanceEnergy 466 G4double resEne = theOsc->GetResonanceEnergy();
465 G4double cutoffEne = theOsc->GetCutoffRecoil 467 G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
466 if (energy > resEne) 468 if (energy > resEne)
467 { 469 {
468 G4double cp1Sq = (energy-resEne)*(energy 470 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
469 G4double cp1 = std::sqrt(cp1Sq); 471 G4double cp1 = std::sqrt(cp1Sq);
470 472
471 //Distant longitudinal interactions 473 //Distant longitudinal interactions
472 G4double QM = 0; 474 G4double QM = 0;
473 if (resEne > 1e-6*energy) 475 if (resEne > 1e-6*energy)
474 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_ma 476 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
475 else 477 else
476 { 478 {
477 QM = resEne*resEne/(beta*2.0*electron_mass 479 QM = resEne*resEne/(beta*2.0*electron_mass_c2);
478 QM = QM*(1.0-0.5*QM/electron_mass_c2); 480 QM = QM*(1.0-0.5*QM/electron_mass_c2);
479 } 481 }
480 G4double SDL1 = 0; 482 G4double SDL1 = 0;
481 if (QM < cutoffEne) 483 if (QM < cutoffEne)
482 SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass << 484 SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
483 485
484 //Distant transverse interactions 486 //Distant transverse interactions
485 if (SDL1) 487 if (SDL1)
486 { 488 {
487 G4double SDT1 = std::max(XHDT0-delta,0.0); 489 G4double SDT1 = std::max(XHDT0-delta,0.0);
488 G4double SD1 = SDL1+SDT1; 490 G4double SD1 = SDL1+SDT1;
489 if (cut > resEne) 491 if (cut > resEne)
490 { 492 {
491 S1 = SD1; //XS1 493 S1 = SD1; //XS1
492 S0 = SD1/resEne; //XS0 494 S0 = SD1/resEne; //XS0
493 S2 = SD1*resEne; //XS2 495 S2 = SD1*resEne; //XS2
494 } 496 }
495 else 497 else
496 { 498 {
497 H1 = SD1; //XH1 499 H1 = SD1; //XH1
498 H0 = SD1/resEne; //XH0 500 H0 = SD1/resEne; //XH0
499 H2 = SD1*resEne; //XH2 501 H2 = SD1*resEne; //XH2
500 } 502 }
501 } 503 }
502 } 504 }
503 // 505 //
504 // Close collisions (Moller's cross section) 506 // Close collisions (Moller's cross section)
505 // 507 //
506 G4double wl = std::max(cut,cutoffEne); 508 G4double wl = std::max(cut,cutoffEne);
507 G4double ee = energy + ionEnergy; 509 G4double ee = energy + ionEnergy;
508 G4double wu = 0.5*ee; 510 G4double wu = 0.5*ee;
509 if (wl < wu-(1e-5*eV)) 511 if (wl < wu-(1e-5*eV))
510 { 512 {
511 H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1 513 H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) +
512 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/ << 514 (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
513 amol*(wu-wl)/(ee*ee); 515 amol*(wu-wl)/(ee*ee);
514 H1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee- << 516 H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
515 (2.0-amol)*G4Log((ee-wu)/(ee-wl)) + << 517 (2.0-amol)*std::log((ee-wu)/(ee-wl)) +
516 amol*(wu*wu-wl*wl)/(2.0*ee*ee); 518 amol*(wu*wu-wl*wl)/(2.0*ee*ee);
517 H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu) 519 H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) -
518 (wl*(2.0*ee-wl)/(ee-wl)) + 520 (wl*(2.0*ee-wl)/(ee-wl)) +
519 (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) + << 521 (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) +
520 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); 522 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
521 wu = wl; 523 wu = wl;
522 } 524 }
523 wl = cutoffEne; 525 wl = cutoffEne;
524 526
525 if (wl > wu-(1e-5*eV)) 527 if (wl > wu-(1e-5*eV))
526 { 528 {
527 (*result)[0] = constant*H0; 529 (*result)[0] = constant*H0;
528 (*result)[1] = constant*H1; 530 (*result)[1] = constant*H1;
529 (*result)[2] = constant*H2; 531 (*result)[2] = constant*H2;
530 (*result)[3] = constant*S0; 532 (*result)[3] = constant*S0;
531 (*result)[4] = constant*S1; 533 (*result)[4] = constant*S1;
532 (*result)[5] = constant*S2; 534 (*result)[5] = constant*S2;
533 return result; 535 return result;
534 } 536 }
535 537
536 S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) 538 S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) +
537 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu) << 539 (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
538 amol*(wu-wl)/(ee*ee); 540 amol*(wu-wl)/(ee*ee);
539 S1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) << 541 S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
540 (2.0-amol)*G4Log((ee-wu)/(ee-wl)) + << 542 (2.0-amol)*std::log((ee-wu)/(ee-wl)) +
541 amol*(wu*wu-wl*wl)/(2.0*ee*ee); 543 amol*(wu*wu-wl*wl)/(2.0*ee*ee);
542 S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee 544 S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) -
543 (wl*(2.0*ee-wl)/(ee-wl)) + 545 (wl*(2.0*ee-wl)/(ee-wl)) +
544 (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) + << 546 (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) +
545 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); 547 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
546 548
547 (*result)[0] = constant*H0; 549 (*result)[0] = constant*H0;
548 (*result)[1] = constant*H1; 550 (*result)[1] = constant*H1;
549 (*result)[2] = constant*H2; 551 (*result)[2] = constant*H2;
550 (*result)[3] = constant*S0; 552 (*result)[3] = constant*S0;
551 (*result)[4] = constant*S1; 553 (*result)[4] = constant*S1;
552 (*result)[5] = constant*S2; 554 (*result)[5] = constant*S2;
553 return result; 555 return result;
554 } 556 }
555 //....oooOO0OOooo........oooOO0OOooo........oo 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
556 G4DataVector* G4PenelopeIonisationXSHandler::C 558 G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsPositron(G4PenelopeOscillator* theOsc,
557 G4double energy, 559 G4double energy,
558 G4double cut, 560 G4double cut,
559 G4double delta) 561 G4double delta)
560 { 562 {
561 // 563 //
562 //This method calculates the hard and soft c 564 //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for
563 //the given oscillator/cut and at the given 565 //the given oscillator/cut and at the given energy.
564 //It returns a G4DataVector* with 6 entries 566 //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
565 //Equivalent of subroutines PINaT1 of Penelo 567 //Equivalent of subroutines PINaT1 of Penelope
566 // 568 //
567 // Results are _per target electron_ 569 // Results are _per target electron_
568 // 570 //
569 G4DataVector* result = new G4DataVector(); 571 G4DataVector* result = new G4DataVector();
570 for (size_t i=0;i<6;i++) 572 for (size_t i=0;i<6;i++)
571 result->push_back(0.); 573 result->push_back(0.);
572 G4double ionEnergy = theOsc->GetIonisationEn 574 G4double ionEnergy = theOsc->GetIonisationEnergy();
573 575
574 //return a set of zero's if the energy it to 576 //return a set of zero's if the energy it too low to excite the current oscillator
575 if (energy < ionEnergy) 577 if (energy < ionEnergy)
576 return result; 578 return result;
577 579
578 G4double H0=0.,H1=0.,H2=0.; 580 G4double H0=0.,H1=0.,H2=0.;
579 G4double S0=0.,S1=0.,S2=0.; 581 G4double S0=0.,S1=0.,S2=0.;
580 582
581 //Define useful constants to be used in the 583 //Define useful constants to be used in the calculation
582 G4double gamma = 1.0+energy/electron_mass_c2 584 G4double gamma = 1.0+energy/electron_mass_c2;
583 G4double gammaSq = gamma*gamma; 585 G4double gammaSq = gamma*gamma;
584 G4double beta = (gammaSq-1.0)/gammaSq; 586 G4double beta = (gammaSq-1.0)/gammaSq;
585 G4double pielr2 = pi*classic_electr_radius*c 587 G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2
586 G4double constant = pielr2*2.0*electron_mass 588 G4double constant = pielr2*2.0*electron_mass_c2/beta;
587 G4double XHDT0 = G4Log(gammaSq)-beta; << 589 G4double XHDT0 = std::log(gammaSq)-beta;
588 590
589 G4double cpSq = energy*(energy+2.0*electron_ 591 G4double cpSq = energy*(energy+2.0*electron_mass_c2);
590 G4double cp = std::sqrt(cpSq); 592 G4double cp = std::sqrt(cpSq);
591 G4double amol = (energy/(energy+electron_mas 593 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
592 G4double g12 = (gamma+1.0)*(gamma+1.0); 594 G4double g12 = (gamma+1.0)*(gamma+1.0);
593 //Bhabha coefficients 595 //Bhabha coefficients
594 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq- 596 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0);
595 G4double bha2 = amol*(3.0+1.0/g12); 597 G4double bha2 = amol*(3.0+1.0/g12);
596 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g 598 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12;
597 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0) 599 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12;
598 600
599 // 601 //
600 // Distant interactions 602 // Distant interactions
601 // 603 //
602 G4double resEne = theOsc->GetResonanceEnergy 604 G4double resEne = theOsc->GetResonanceEnergy();
603 G4double cutoffEne = theOsc->GetCutoffRecoil 605 G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
604 if (energy > resEne) 606 if (energy > resEne)
605 { 607 {
606 G4double cp1Sq = (energy-resEne)*(energy 608 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
607 G4double cp1 = std::sqrt(cp1Sq); 609 G4double cp1 = std::sqrt(cp1Sq);
608 610
609 //Distant longitudinal interactions 611 //Distant longitudinal interactions
610 G4double QM = 0; 612 G4double QM = 0;
611 if (resEne > 1e-6*energy) 613 if (resEne > 1e-6*energy)
612 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_ma 614 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
613 else 615 else
614 { 616 {
615 QM = resEne*resEne/(beta*2.0*electron_mass 617 QM = resEne*resEne/(beta*2.0*electron_mass_c2);
616 QM = QM*(1.0-0.5*QM/electron_mass_c2); 618 QM = QM*(1.0-0.5*QM/electron_mass_c2);
617 } 619 }
618 G4double SDL1 = 0; 620 G4double SDL1 = 0;
619 if (QM < cutoffEne) 621 if (QM < cutoffEne)
620 SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass << 622 SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
621 623
622 //Distant transverse interactions 624 //Distant transverse interactions
623 if (SDL1) 625 if (SDL1)
624 { 626 {
625 G4double SDT1 = std::max(XHDT0-delta,0.0); 627 G4double SDT1 = std::max(XHDT0-delta,0.0);
626 G4double SD1 = SDL1+SDT1; 628 G4double SD1 = SDL1+SDT1;
627 if (cut > resEne) 629 if (cut > resEne)
628 { 630 {
629 S1 = SD1; //XS1 631 S1 = SD1; //XS1
630 S0 = SD1/resEne; //XS0 632 S0 = SD1/resEne; //XS0
631 S2 = SD1*resEne; //XS2 633 S2 = SD1*resEne; //XS2
632 } 634 }
633 else 635 else
634 { 636 {
635 H1 = SD1; //XH1 637 H1 = SD1; //XH1
636 H0 = SD1/resEne; //XH0 638 H0 = SD1/resEne; //XH0
637 H2 = SD1*resEne; //XH2 639 H2 = SD1*resEne; //XH2
638 } 640 }
639 } 641 }
640 } 642 }
>> 643
641 // 644 //
642 // Close collisions (Bhabha's cross section) 645 // Close collisions (Bhabha's cross section)
643 // 646 //
644 G4double wl = std::max(cut,cutoffEne); 647 G4double wl = std::max(cut,cutoffEne);
645 G4double wu = energy; 648 G4double wu = energy;
646 G4double energySq = energy*energy; 649 G4double energySq = energy*energy;
647 if (wl < wu-(1e-5*eV)) 650 if (wl < wu-(1e-5*eV))
648 { 651 {
649 G4double wlSq = wl*wl; 652 G4double wlSq = wl*wl;
650 G4double wuSq = wu*wu; 653 G4double wuSq = wu*wu;
651 H0 += (1.0/wl) - (1.0/wu)- bha1*G4Log(wu << 654 H0 += (1.0/wl) - (1.0/wu)- bha1*std::log(wu/wl)/energy
652 + bha2*(wu-wl)/energySq 655 + bha2*(wu-wl)/energySq
653 - bha3*(wuSq-wlSq)/(2.0*energySq*energy) 656 - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
654 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energ 657 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
655 H1 += G4Log(wu/wl) - bha1*(wu-wl)/energy << 658 H1 += std::log(wu/wl) - bha1*(wu-wl)/energy
656 + bha2*(wuSq-wlSq)/(2.0*energySq) 659 + bha2*(wuSq-wlSq)/(2.0*energySq)
657 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energ 660 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
658 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*e 661 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
659 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*en 662 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
660 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) 663 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
661 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*e 664 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
662 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*ener 665 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
663 wu = wl; 666 wu = wl;
664 } 667 }
665 wl = cutoffEne; 668 wl = cutoffEne;
666 669
667 if (wl > wu-(1e-5*eV)) 670 if (wl > wu-(1e-5*eV))
668 { 671 {
669 (*result)[0] = constant*H0; 672 (*result)[0] = constant*H0;
670 (*result)[1] = constant*H1; 673 (*result)[1] = constant*H1;
671 (*result)[2] = constant*H2; 674 (*result)[2] = constant*H2;
672 (*result)[3] = constant*S0; 675 (*result)[3] = constant*S0;
673 (*result)[4] = constant*S1; 676 (*result)[4] = constant*S1;
674 (*result)[5] = constant*S2; 677 (*result)[5] = constant*S2;
675 return result; 678 return result;
676 } 679 }
677 680
678 G4double wlSq = wl*wl; 681 G4double wlSq = wl*wl;
679 G4double wuSq = wu*wu; 682 G4double wuSq = wu*wu;
680 683
681 S0 += (1.0/wl) - (1.0/wu) - bha1*G4Log(wu/wl << 684 S0 += (1.0/wl) - (1.0/wu) - bha1*std::log(wu/wl)/energy
682 + bha2*(wu-wl)/energySq 685 + bha2*(wu-wl)/energySq
683 - bha3*(wuSq-wlSq)/(2.0*energySq*energy) 686 - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
684 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*ene 687 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
685 688
686 S1 += G4Log(wu/wl) - bha1*(wu-wl)/energy << 689 S1 += std::log(wu/wl) - bha1*(wu-wl)/energy
687 + bha2*(wuSq-wlSq)/(2.0*energySq) 690 + bha2*(wuSq-wlSq)/(2.0*energySq)
688 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*ene 691 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy)
689 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq 692 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
690 693
691 S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy 694 S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
692 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) 695 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
693 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq 696 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
694 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*en 697 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
695 698
696 (*result)[0] = constant*H0; 699 (*result)[0] = constant*H0;
697 (*result)[1] = constant*H1; 700 (*result)[1] = constant*H1;
698 (*result)[2] = constant*H2; 701 (*result)[2] = constant*H2;
699 (*result)[3] = constant*S0; 702 (*result)[3] = constant*S0;
700 (*result)[4] = constant*S1; 703 (*result)[4] = constant*S1;
701 (*result)[5] = constant*S2; 704 (*result)[5] = constant*S2;
702 705
703 return result; 706 return result;
704 } 707 }
705 708