| 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: G4PenelopeBremsstrahlungFS.cc 101820 2016-12-01 08:34:05Z gcosmo $
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 23 Nov 2010 L Pandola First complete i 32 // 23 Nov 2010 L Pandola First complete implementation
32 // 02 May 2011 L.Pandola Remove dependenc 33 // 02 May 2011 L.Pandola Remove dependency on CLHEP::HepMatrix
33 // 24 May 2011 L.Pandola Renamed (make v2 34 // 24 May 2011 L.Pandola Renamed (make v2008 as default Penelope)
34 // 03 Oct 2013 L.Pandola Migration to MT 35 // 03 Oct 2013 L.Pandola Migration to MT
35 // 30 Oct 2013 L.Pandola Use G4Cache to a 36 // 30 Oct 2013 L.Pandola Use G4Cache to avoid new/delete of the
36 // data vector on 37 // data vector on the fly in SampleGammaEnergy()
37 // 38 //
38 39
39 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 #include "G4PenelopeBremsstrahlungFS.hh" 41 #include "G4PenelopeBremsstrahlungFS.hh"
41 #include "G4PhysicsFreeVector.hh" 42 #include "G4PhysicsFreeVector.hh"
42 #include "G4PhysicsTable.hh" 43 #include "G4PhysicsTable.hh"
43 #include "G4Material.hh" 44 #include "G4Material.hh"
44 #include "Randomize.hh" 45 #include "Randomize.hh"
45 #include "G4AutoDelete.hh" 46 #include "G4AutoDelete.hh"
46 #include "G4Exp.hh" 47 #include "G4Exp.hh"
47 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
49 50
50 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 52
52 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstr 53 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstrahlungFS(G4int verbosity) :
53 fReducedXSTable(nullptr),fEffectiveZSq(nullp << 54 theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
54 fPBcut(nullptr),fVerbosity(verbosity) << 55 thePBcut(0),fVerbosity(verbosity)
55 { 56 {
56 fCache.Put(0); 57 fCache.Put(0);
57 G4double tempvector[fNBinsX] = << 58 G4double tempvector[nBinsX] =
58 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15 59 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
59 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6 60 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
60 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0 61 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
61 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e 62 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
62 63
63 for (std::size_t ix=0;ix<fNBinsX;++ix) << 64 for (size_t ix=0;ix<nBinsX;ix++)
64 theXGrid[ix] = tempvector[ix]; 65 theXGrid[ix] = tempvector[ix];
65 66
66 for (std::size_t i=0;i<fNBinsE;++i) << 67 for (size_t i=0;i<nBinsE;i++)
67 theEGrid[i] = 0.; 68 theEGrid[i] = 0.;
68 69
69 fElementData = new std::map<G4int,G4DataVect << 70 theElementData = new std::map<G4int,G4DataVector*>;
70 } 71 }
71 72
72 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 74
74 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsst 75 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsstrahlungFS()
75 { 76 {
76 ClearTables(); 77 ClearTables();
77 78
78 //The G4Physics*Vector pointers contained in 79 //The G4Physics*Vector pointers contained in the fCache are automatically deleted by
79 //the G4AutoDelete so there is no need to ta 80 //the G4AutoDelete so there is no need to take care of them manually
80 81
81 //Clear manually fElementData << 82 //Clear manually theElementData
82 if (fElementData) << 83 if (theElementData)
83 { 84 {
84 for (auto& item : (*fElementData)) << 85 for (auto& item : (*theElementData))
85 delete item.second; 86 delete item.second;
86 delete fElementData; << 87 delete theElementData;
87 fElementData = nullptr; << 88 theElementData = nullptr;
88 } 89 }
>> 90
89 } 91 }
90 92
91 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
92 94
93 95
94 void G4PenelopeBremsstrahlungFS::ClearTables(G 96 void G4PenelopeBremsstrahlungFS::ClearTables(G4bool isMaster)
95 { 97 {
96 //Just to check 98 //Just to check
97 if (!isMaster) 99 if (!isMaster)
98 G4Exception("G4PenelopeBremsstrahlungFS::C 100 G4Exception("G4PenelopeBremsstrahlungFS::ClearTables()",
99 "em0100",FatalException,"Worker thread in 101 "em0100",FatalException,"Worker thread in this method");
100 102
101 if (fReducedXSTable) << 103
>> 104 if (theReducedXSTable)
102 { 105 {
103 for (auto& item : (*fReducedXSTable)) << 106 for (auto& item : (*theReducedXSTable))
104 { 107 {
105 G4PhysicsTable* tab = item.second; << 108 //G4PhysicsTable* tab = item.second;
106 tab->clearAndDestroy(); << 109 //tab->clearAndDestroy();
107 delete tab; << 110 delete item.second;
108 } 111 }
109 fReducedXSTable->clear(); << 112 delete theReducedXSTable;
110 delete fReducedXSTable; << 113 theReducedXSTable = nullptr;
111 fReducedXSTable = nullptr; <<
112 } 114 }
113 115
114 if (fSamplingTable) << 116 if (theSamplingTable)
115 { 117 {
116 for (auto& item : (*fSamplingTable)) << 118 for (auto& item : (*theSamplingTable))
117 { 119 {
118 G4PhysicsTable* tab = item.second; << 120 //G4PhysicsTable* tab = item.second;
119 tab->clearAndDestroy(); << 121 // tab->clearAndDestroy();
120 delete tab; << 122 delete item.second;
121 } 123 }
122 fSamplingTable->clear(); << 124 delete theSamplingTable;
123 delete fSamplingTable; << 125 theSamplingTable = nullptr;
124 fSamplingTable = nullptr; <<
125 } 126 }
126 if (fPBcut) << 127 if (thePBcut)
127 { 128 {
128 /* 129 /*
129 std::map< std::pair<const G4Material*,G4doub 130 std::map< std::pair<const G4Material*,G4double> ,G4PhysicsFreeVector*>::iterator kk;
130 for (kk=fPBcut->begin(); kk != fPBcut->end() << 131 for (kk=thePBcut->begin(); kk != thePBcut->end(); kk++)
131 delete kk->second; 132 delete kk->second;
132 */ 133 */
133 delete fPBcut; << 134 delete thePBcut;
134 fPBcut = nullptr; << 135 thePBcut = nullptr;
135 } 136 }
136 137
137 if (fEffectiveZSq) << 138
>> 139 if (theEffectiveZSq)
138 { 140 {
139 delete fEffectiveZSq; << 141 delete theEffectiveZSq;
140 fEffectiveZSq = nullptr; << 142 theEffectiveZSq = nullptr;
141 } 143 }
142 144
143 return; 145 return;
144 } 146 }
145 147
146 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147 149
148 G4double G4PenelopeBremsstrahlungFS::GetEffect 150 G4double G4PenelopeBremsstrahlungFS::GetEffectiveZSquared(const G4Material* material) const
149 { 151 {
150 if (!fEffectiveZSq) << 152 if (!theEffectiveZSq)
151 { 153 {
152 G4ExceptionDescription ed; 154 G4ExceptionDescription ed;
153 ed << "The container for the <Z^2> value 155 ed << "The container for the <Z^2> values is not initialized" << G4endl;
154 G4Exception("G4PenelopeBremsstrahlungFS: 156 G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
155 "em2007",FatalException,ed); 157 "em2007",FatalException,ed);
156 return 0; 158 return 0;
157 } 159 }
158 //found in the table: return it 160 //found in the table: return it
159 if (fEffectiveZSq->count(material)) << 161 if (theEffectiveZSq->count(material))
160 return fEffectiveZSq->find(material)->seco << 162 return theEffectiveZSq->find(material)->second;
161 else 163 else
162 { 164 {
163 G4ExceptionDescription ed; 165 G4ExceptionDescription ed;
164 ed << "The value of <Z^2> is not proper 166 ed << "The value of <Z^2> is not properly set for material " <<
165 material->GetName() << G4endl; 167 material->GetName() << G4endl;
166 //requires running of BuildScaledXSTable 168 //requires running of BuildScaledXSTable()
167 G4Exception("G4PenelopeBremsstrahlungFS: 169 G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
168 "em2008",FatalException,ed); 170 "em2008",FatalException,ed);
169 } 171 }
170 return 0; 172 return 0;
171 } 173 }
172 174
173 //....oooOO0OOooo........oooOO0OOooo........oo 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174 176
175 void G4PenelopeBremsstrahlungFS::BuildScaledXS 177 void G4PenelopeBremsstrahlungFS::BuildScaledXSTable(const G4Material* material,
176 G4double cut,G4bool isMaster) 178 G4double cut,G4bool isMaster)
177 { 179 {
178 //Corresponds to subroutines EBRaW and EBRaR 180 //Corresponds to subroutines EBRaW and EBRaR of PENELOPE
179 /* 181 /*
180 This method generates the table of the sca 182 This method generates the table of the scaled energy-loss cross section from
181 bremsstrahlung emission for the given mate 183 bremsstrahlung emission for the given material. Original data are read from
182 file. The table is normalized according to 184 file. The table is normalized according to the Berger-Seltzer cross section.
183 */ 185 */
184 186
185 //Just to check 187 //Just to check
186 if (!isMaster) 188 if (!isMaster)
187 G4Exception("G4PenelopeBremsstrahlungFS::B 189 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
188 "em0100",FatalException,"Worker thread in 190 "em0100",FatalException,"Worker thread in this method");
189 191
190 if (fVerbosity > 2) 192 if (fVerbosity > 2)
191 { 193 {
192 G4cout << "Entering in G4PenelopeBremsst 194 G4cout << "Entering in G4PenelopeBremsstrahlungFS::BuildScaledXSTable for " <<
193 material->GetName() << G4endl; 195 material->GetName() << G4endl;
194 G4cout << "Threshold = " << cut/keV << " 196 G4cout << "Threshold = " << cut/keV << " keV, isMaster= " << isMaster <<
195 G4endl; 197 G4endl;
196 } 198 }
197 199
198 //This method should be accessed by the mast 200 //This method should be accessed by the master only
199 if (!fSamplingTable) << 201 if (!theSamplingTable)
200 fSamplingTable = << 202 theSamplingTable =
201 new std::map< std::pair<const G4Material 203 new std::map< std::pair<const G4Material*,G4double> , G4PhysicsTable*>;
202 if (!fPBcut) << 204 if (!thePBcut)
203 fPBcut = << 205 thePBcut =
204 new std::map< std::pair<const G4Material 206 new std::map< std::pair<const G4Material*,G4double> , G4PhysicsFreeVector* >;
205 207
206 //check if the container exists (if not, cre 208 //check if the container exists (if not, create it)
207 if (!fReducedXSTable) << 209 if (!theReducedXSTable)
208 fReducedXSTable = new std::map< std::pair< << 210 theReducedXSTable = new std::map< std::pair<const G4Material*,G4double> ,
209 G4PhysicsTable*>; 211 G4PhysicsTable*>;
210 if (!fEffectiveZSq) << 212 if (!theEffectiveZSq)
211 fEffectiveZSq = new std::map<const G4Mater << 213 theEffectiveZSq = new std::map<const G4Material*,G4double>;
212 214
213 //****************************************** << 215
>> 216
>> 217 //*********************************************************************
214 //Determine the equivalent atomic number <Z^ 218 //Determine the equivalent atomic number <Z^2>
215 //****************************************** 219 //*********************************************************************
216 std::vector<G4double> *StechiometricFactors 220 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
217 std::size_t nElements = material->GetNumberO << 221 G4int nElements = material->GetNumberOfElements();
218 const G4ElementVector* elementVector = mater 222 const G4ElementVector* elementVector = material->GetElementVector();
219 const G4double* fractionVector = material->G 223 const G4double* fractionVector = material->GetFractionVector();
220 for (std::size_t i=0;i<nElements;i++) << 224 for (G4int i=0;i<nElements;i++)
221 { 225 {
222 G4double fraction = fractionVector[i]; 226 G4double fraction = fractionVector[i];
223 G4double atomicWeigth = (*elementVector) 227 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
224 StechiometricFactors->push_back(fraction 228 StechiometricFactors->push_back(fraction/atomicWeigth);
225 } 229 }
226 //Find max 230 //Find max
227 G4double MaxStechiometricFactor = 0.; 231 G4double MaxStechiometricFactor = 0.;
228 for (std::size_t i=0;i<nElements;i++) << 232 for (G4int i=0;i<nElements;i++)
229 { 233 {
230 if ((*StechiometricFactors)[i] > MaxStec 234 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
231 MaxStechiometricFactor = (*Stechiometr 235 MaxStechiometricFactor = (*StechiometricFactors)[i];
232 } 236 }
233 //Normalize 237 //Normalize
234 for (std::size_t i=0;i<nElements;i++) << 238 for (G4int i=0;i<nElements;i++)
235 (*StechiometricFactors)[i] /= MaxStechiom 239 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
236 240
237 G4double sumz2 = 0; 241 G4double sumz2 = 0;
238 G4double sums = 0; 242 G4double sums = 0;
239 for (std::size_t i=0;i<nElements;i++) << 243 for (G4int i=0;i<nElements;i++)
240 { 244 {
241 G4double Z = (*elementVector)[i]->GetZ() 245 G4double Z = (*elementVector)[i]->GetZ();
242 sumz2 += (*StechiometricFactors)[i]*Z*Z; 246 sumz2 += (*StechiometricFactors)[i]*Z*Z;
243 sums += (*StechiometricFactors)[i]; 247 sums += (*StechiometricFactors)[i];
244 } 248 }
245 G4double ZBR2 = sumz2/sums; 249 G4double ZBR2 = sumz2/sums;
246 250
247 fEffectiveZSq->insert(std::make_pair(materia << 251 theEffectiveZSq->insert(std::make_pair(material,ZBR2));
>> 252
248 253
249 //****************************************** 254 //*********************************************************************
250 // loop on elements and read data files 255 // loop on elements and read data files
251 //****************************************** 256 //*********************************************************************
252 G4DataVector* tempData = new G4DataVector(fN << 257 G4DataVector* tempData = new G4DataVector(nBinsE);
253 G4DataVector* tempMatrix = new G4DataVector( << 258 G4DataVector* tempMatrix = new G4DataVector(nBinsE*nBinsX,0.);
254 259
255 for (std::size_t iel=0;iel<nElements;iel++) << 260 for (G4int iel=0;iel<nElements;iel++)
256 { 261 {
257 G4double Z = (*elementVector)[iel]->GetZ 262 G4double Z = (*elementVector)[iel]->GetZ();
258 G4int iZ = (G4int) Z; 263 G4int iZ = (G4int) Z;
259 G4double wgt = (*StechiometricFactors)[i 264 G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
260 << 265 //
>> 266
261 //the element is not already loaded 267 //the element is not already loaded
262 if (!fElementData->count(iZ)) << 268 if (!theElementData->count(iZ))
263 { 269 {
264 ReadDataFile(iZ); 270 ReadDataFile(iZ);
265 if (!fElementData->count(iZ)) << 271 if (!theElementData->count(iZ))
266 { 272 {
267 G4ExceptionDescription ed; 273 G4ExceptionDescription ed;
268 ed << "Error in G4PenelopeBremsstrahlu 274 ed << "Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" << G4endl;
269 ed << "Unable to retrieve data for ele 275 ed << "Unable to retrieve data for element " << iZ << G4endl;
270 G4Exception("G4PenelopeBremsstrahlungF 276 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
271 "em2009",FatalException,ed); 277 "em2009",FatalException,ed);
272 } 278 }
273 } 279 }
274 280
275 G4DataVector* atomData = fElementData->f << 281 G4DataVector* atomData = theElementData->find(iZ)->second;
276 282
277 for (std::size_t ie=0;ie<fNBinsE;++ie) << 283 for (size_t ie=0;ie<nBinsE;ie++)
278 { 284 {
279 (*tempData)[ie] += wgt*(*atomData)[ie*(fNB << 285 (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX]; //last column contains total XS
280 for (std::size_t ix=0;ix<fNBinsX;++ix) << 286 for (size_t ix=0;ix<nBinsX;ix++)
281 (*tempMatrix)[ie*fNBinsX+ix] += wgt*(*at << 287 (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
282 } 288 }
283 } 289 }
284 290
285 //****************************************** 291 //*********************************************************************
286 // the total energy loss spectrum is re-norm 292 // the total energy loss spectrum is re-normalized to reproduce the total
287 // scaled cross section of Berger and Seltze 293 // scaled cross section of Berger and Seltzer
288 //****************************************** 294 //*********************************************************************
289 for (std::size_t ie=0;ie<fNBinsE;++ie) << 295 for (size_t ie=0;ie<nBinsE;ie++)
290 { 296 {
291 //for each energy, calculate integral of 297 //for each energy, calculate integral of dSigma/dx over dx
292 G4double* tempData2 = new G4double[fNBin << 298 G4double* tempData2 = new G4double[nBinsX];
293 for (std::size_t ix=0;ix<fNBinsX;++ix) << 299 for (size_t ix=0;ix<nBinsX;ix++)
294 tempData2[ix] = (*tempMatrix)[ie*fNBinsX+ix] << 300 tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
295 G4double rsum = GetMomentumIntegral(temp 301 G4double rsum = GetMomentumIntegral(tempData2,1.0,0);
296 delete[] tempData2; 302 delete[] tempData2;
297 G4double fact = millibarn*(theEGrid[ie]+ 303 G4double fact = millibarn*(theEGrid[ie]+electron_mass_c2)*(1./fine_structure_const)/
298 (classic_electr_radius*classic_electr_radius 304 (classic_electr_radius*classic_electr_radius*(theEGrid[ie]+2.0*electron_mass_c2));
299 G4double fnorm = (*tempData)[ie]/(rsum*f 305 G4double fnorm = (*tempData)[ie]/(rsum*fact);
300 G4double TST = 100.*std::fabs(fnorm-1.0) 306 G4double TST = 100.*std::fabs(fnorm-1.0);
301 if (TST > 1.0) 307 if (TST > 1.0)
302 { 308 {
303 G4ExceptionDescription ed; 309 G4ExceptionDescription ed;
304 ed << "G4PenelopeBremsstrahlungFS. Corrupt 310 ed << "G4PenelopeBremsstrahlungFS. Corrupted data files?" << G4endl;
305 G4cout << "TST= " << TST << "; fnorm = " < 311 G4cout << "TST= " << TST << "; fnorm = " << fnorm << G4endl;
306 G4cout << "rsum = " << rsum << G4endl; 312 G4cout << "rsum = " << rsum << G4endl;
307 G4cout << "fact = " << fact << G4endl; 313 G4cout << "fact = " << fact << G4endl;
308 G4cout << ie << " " << theEGrid[ie]/keV << 314 G4cout << ie << " " << theEGrid[ie]/keV << " " << (*tempData)[ie]/barn << G4endl;
309 G4Exception("G4PenelopeBremsstrahlungFS::B 315 G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
310 "em2010",FatalException,ed); 316 "em2010",FatalException,ed);
311 } 317 }
312 for (std::size_t ix=0;ix<fNBinsX;++ix) << 318 for (size_t ix=0;ix<nBinsX;ix++)
313 (*tempMatrix)[ie*fNBinsX+ix] *= fnorm; << 319 (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
314 } 320 }
315 321
316 //****************************************** 322 //*********************************************************************
317 // create and fill the tables 323 // create and fill the tables
318 //****************************************** 324 //*********************************************************************
319 G4PhysicsTable* thePhysicsTable = new G4Phys 325 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
320 // the table will contain 32 G4PhysicsFreeVe 326 // the table will contain 32 G4PhysicsFreeVectors with different
321 // values of x. Each of the G4PhysicsFreeVec 327 // values of x. Each of the G4PhysicsFreeVectors has a profile of
322 // log(XS) vs. log(E) 328 // log(XS) vs. log(E)
323 329
324 //reserve space of the vectors. Everything i 330 //reserve space of the vectors. Everything is log-log
325 //I add one extra "fake" point at low energy 331 //I add one extra "fake" point at low energy, since the Penelope
326 //table starts at 1 keV 332 //table starts at 1 keV
327 for (std::size_t i=0;i<fNBinsX;++i) << 333 for (size_t i=0;i<nBinsX;i++)
328 thePhysicsTable->push_back(new G4PhysicsFr << 334 thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsE+1));
329 335
330 for (std::size_t ix=0;ix<fNBinsX;++ix) << 336 for (size_t ix=0;ix<nBinsX;ix++)
331 { 337 {
332 G4PhysicsFreeVector* theVec = 338 G4PhysicsFreeVector* theVec =
333 (G4PhysicsFreeVector*) ((*thePhysicsTable)[i 339 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]);
334 for (std::size_t ie=0;ie<fNBinsE;++ie) << 340 for (size_t ie=0;ie<nBinsE;ie++)
335 { 341 {
336 G4double logene = G4Log(theEGrid[ie]); << 342 G4double logene = std::log(theEGrid[ie]);
337 G4double aValue = (*tempMatrix)[ie*fNBinsX << 343 G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
338 if (aValue < 1e-20*millibarn) //protection 344 if (aValue < 1e-20*millibarn) //protection against log(0)
339 aValue = 1e-20*millibarn; 345 aValue = 1e-20*millibarn;
340 theVec->PutValues(ie+1,logene,G4Log(aValue << 346 theVec->PutValue(ie+1,logene,std::log(aValue));
341 } 347 }
342 //Add fake point at 1 eV using an extrap 348 //Add fake point at 1 eV using an extrapolation with the derivative
343 //at the first valid point (Penelope app 349 //at the first valid point (Penelope approach)
344 G4double derivative = ((*theVec)[2]-(*th 350 G4double derivative = ((*theVec)[2]-(*theVec)[1])/(theVec->Energy(2) - theVec->Energy(1));
345 G4double log1eV = G4Log(1*eV); << 351 G4double log1eV = std::log(1*eV);
346 G4double val1eV = (*theVec)[1]+derivativ 352 G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->Energy(1));
347 //fake point at very low energy 353 //fake point at very low energy
348 theVec->PutValues(0,log1eV,val1eV); << 354 theVec->PutValue(0,log1eV,val1eV);
349 } 355 }
350 std::pair<const G4Material*,G4double> theKey 356 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
351 fReducedXSTable->insert(std::make_pair(theKe << 357 theReducedXSTable->insert(std::make_pair(theKey,thePhysicsTable));
352 358
353 delete StechiometricFactors; 359 delete StechiometricFactors;
354 delete tempData; 360 delete tempData;
355 delete tempMatrix; 361 delete tempMatrix;
356 362
357 //Do here also the initialization of the ene 363 //Do here also the initialization of the energy sampling
358 if (!(fSamplingTable->count(theKey))) << 364 if (!(theSamplingTable->count(theKey)))
359 InitializeEnergySampling(material,cut); 365 InitializeEnergySampling(material,cut);
360 366
361 return; 367 return;
362 } 368 }
363 369
364 //....oooOO0OOooo........oooOO0OOooo........oo 370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365 371
366 void G4PenelopeBremsstrahlungFS::ReadDataFile( 372 void G4PenelopeBremsstrahlungFS::ReadDataFile(G4int Z)
367 { 373 {
368 const char* path = G4FindDataDir("G4LEDATA") << 374
>> 375 char* path = getenv("G4LEDATA");
369 if (!path) 376 if (!path)
370 { 377 {
371 G4String excep = "G4PenelopeBremsstrahlu 378 G4String excep = "G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
372 G4Exception("G4PenelopeBremsstrahlungFS: 379 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
373 "em0006",FatalException,excep); 380 "em0006",FatalException,excep);
374 return; 381 return;
375 } 382 }
376 /* 383 /*
377 Read the cross section file 384 Read the cross section file
378 */ 385 */
379 std::ostringstream ost; 386 std::ostringstream ost;
380 if (Z>9) 387 if (Z>9)
381 ost << path << "/penelope/bremsstrahlung/p 388 ost << path << "/penelope/bremsstrahlung/pdebr" << Z << ".p08";
382 else 389 else
383 ost << path << "/penelope/bremsstrahlung/p 390 ost << path << "/penelope/bremsstrahlung/pdebr0" << Z << ".p08";
384 std::ifstream file(ost.str().c_str()); 391 std::ifstream file(ost.str().c_str());
385 if (!file.is_open()) 392 if (!file.is_open())
386 { 393 {
387 G4String excep = "G4PenelopeBremsstrahlu 394 G4String excep = "G4PenelopeBremsstrahlungFS - data file " +
388 G4String(ost.str()) + " not found!"; 395 G4String(ost.str()) + " not found!";
389 G4Exception("G4PenelopeBremsstrahlungFS: 396 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
390 "em0003",FatalException,excep); 397 "em0003",FatalException,excep);
391 return; 398 return;
392 } 399 }
393 400
394 G4int readZ =0; 401 G4int readZ =0;
395 file >> readZ; 402 file >> readZ;
396 403
397 //check the right file is opened. 404 //check the right file is opened.
398 if (readZ != Z) 405 if (readZ != Z)
399 { 406 {
400 G4ExceptionDescription ed; 407 G4ExceptionDescription ed;
401 ed << "Corrupted data file for Z=" << Z 408 ed << "Corrupted data file for Z=" << Z << G4endl;
402 G4Exception("G4PenelopeBremsstrahlungFS: 409 G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
403 "em0005",FatalException,ed); 410 "em0005",FatalException,ed);
404 return; 411 return;
405 } 412 }
406 413
407 G4DataVector* theMatrix = new G4DataVector(f << 414 G4DataVector* theMatrix = new G4DataVector(nBinsE*(nBinsX+1),0.); //initialized with zeros
408 415
409 for (std::size_t ie=0;ie<fNBinsE;++ie) << 416 for (size_t ie=0;ie<nBinsE;ie++)
410 { 417 {
411 G4double myDouble = 0; 418 G4double myDouble = 0;
412 file >> myDouble; //energy (eV) 419 file >> myDouble; //energy (eV)
413 if (!theEGrid[ie]) //fill only the first 420 if (!theEGrid[ie]) //fill only the first time
414 theEGrid[ie] = myDouble*eV; 421 theEGrid[ie] = myDouble*eV;
415 // 422 //
416 for (std::size_t ix=0;ix<fNBinsX;++ix) << 423 for (size_t ix=0;ix<nBinsX;ix++)
417 { 424 {
418 file >> myDouble; 425 file >> myDouble;
419 (*theMatrix)[ie*(fNBinsX+1)+ix] = myDouble << 426 (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*millibarn;
420 } 427 }
421 file >> myDouble; //total cross section 428 file >> myDouble; //total cross section
422 (*theMatrix)[ie*(fNBinsX+1)+fNBinsX] = m << 429 (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*millibarn;
423 } 430 }
424 431
425 if (fElementData) << 432 if (theElementData)
426 fElementData->insert(std::make_pair(Z,theM << 433 theElementData->insert(std::make_pair(Z,theMatrix));
427 else 434 else
428 delete theMatrix; 435 delete theMatrix;
429 file.close(); 436 file.close();
430 return; 437 return;
431 } 438 }
432 <<
433 //....oooOO0OOooo........oooOO0OOooo........oo 439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434 440
435 G4double G4PenelopeBremsstrahlungFS::GetMoment 441 G4double G4PenelopeBremsstrahlungFS::GetMomentumIntegral(G4double* y,
436 G4double xup,G4int momOrder) 442 G4double xup,G4int momOrder) const
437 //x is always the gridX 443 //x is always the gridX
438 { 444 {
439 //Corresponds to the function RLMOM of Penel 445 //Corresponds to the function RLMOM of Penelope
440 //This method performs the calculation of th 446 //This method performs the calculation of the integral of (x^momOrder)*y over the interval
441 //from x[0] to xup, obtained by linear inter 447 //from x[0] to xup, obtained by linear interpolation on a table of y.
442 //The independent variable is assumed to tak 448 //The independent variable is assumed to take positive values only.
443 // 449 //
444 std::size_t size = fNBinsX; << 450 size_t size = nBinsX;
445 const G4double eps = 1e-35; 451 const G4double eps = 1e-35;
446 452
447 //Check that the call is valid 453 //Check that the call is valid
448 if (momOrder<-1 || size<2 || theXGrid[0]<0) 454 if (momOrder<-1 || size<2 || theXGrid[0]<0)
449 { 455 {
450 G4Exception("G4PenelopeBremsstrahlungFS: 456 G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
451 "em2011",FatalException,"Invalid call"); 457 "em2011",FatalException,"Invalid call");
452 } 458 }
453 459
454 for (std::size_t i=1;i<size;++i) << 460 for (size_t i=1;i<size;i++)
455 { 461 {
456 if (theXGrid[i]<0 || theXGrid[i]<theXGri 462 if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
457 { 463 {
458 G4ExceptionDescription ed; 464 G4ExceptionDescription ed;
459 ed << "Invalid call for bin " << i << G4en 465 ed << "Invalid call for bin " << i << G4endl;
460 G4Exception("G4PenelopeBremsstrahlungFS::G 466 G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
461 "em2012",FatalException,ed); 467 "em2012",FatalException,ed);
462 } 468 }
463 } 469 }
464 470
465 //Compute the integral 471 //Compute the integral
466 G4double result = 0; 472 G4double result = 0;
467 if (xup < theXGrid[0]) 473 if (xup < theXGrid[0])
468 return result; 474 return result;
469 G4bool loopAgain = true; << 475 bool loopAgain = true;
470 G4double xt = std::min(xup,theXGrid[size-1]) 476 G4double xt = std::min(xup,theXGrid[size-1]);
471 G4double xtc = 0; 477 G4double xtc = 0;
472 for (std::size_t i=0;i<size-1;++i) << 478 for (size_t i=0;i<size-1;i++)
473 { 479 {
474 G4double x1 = std::max(theXGrid[i],eps); 480 G4double x1 = std::max(theXGrid[i],eps);
475 G4double y1 = y[i]; 481 G4double y1 = y[i];
476 G4double x2 = std::max(theXGrid[i+1],eps 482 G4double x2 = std::max(theXGrid[i+1],eps);
477 G4double y2 = y[i+1]; 483 G4double y2 = y[i+1];
478 if (xt < x2) 484 if (xt < x2)
479 { 485 {
480 xtc = xt; 486 xtc = xt;
481 loopAgain = false; 487 loopAgain = false;
482 } 488 }
483 else 489 else
484 xtc = x2; 490 xtc = x2;
485 G4double dx = x2-x1; 491 G4double dx = x2-x1;
486 G4double dy = y2-y1; 492 G4double dy = y2-y1;
487 G4double ds = 0; 493 G4double ds = 0;
488 if (std::fabs(dx)>1e-14*std::fabs(dy)) 494 if (std::fabs(dx)>1e-14*std::fabs(dy))
489 { 495 {
490 G4double b=dy/dx; 496 G4double b=dy/dx;
491 G4double a=y1-b*x1; 497 G4double a=y1-b*x1;
492 if (momOrder == -1) 498 if (momOrder == -1)
493 ds = a*G4Log(xtc/x1)+b*(xtc-x1); << 499 ds = a*std::log(xtc/x1)+b*(xtc-x1);
494 else if (momOrder == 0) //speed it up, not 500 else if (momOrder == 0) //speed it up, not using pow()
495 ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1); 501 ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
496 else 502 else
497 ds = a*(std::pow(xtc,momOrder+1)-std::po 503 ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((G4double) (momOrder + 1))
498 + b*(std::pow(xtc,momOrder+2)-std::pow 504 + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((G4double) (momOrder + 2));
499 } 505 }
500 else 506 else
501 ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOr 507 ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder);
502 result += ds; 508 result += ds;
503 if (!loopAgain) 509 if (!loopAgain)
504 return result; 510 return result;
505 } 511 }
506 return result; 512 return result;
507 } 513 }
508 514
509 //....oooOO0OOooo........oooOO0OOooo........oo 515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
510 516
511 const G4PhysicsTable* G4PenelopeBremsstrahlung 517 const G4PhysicsTable* G4PenelopeBremsstrahlungFS::GetScaledXSTable(const G4Material* mat,
512 const G4double cut) const 518 const G4double cut) const
513 { 519 {
514 //check if it already contains the entry 520 //check if it already contains the entry
515 std::pair<const G4Material*,G4double> theKey 521 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
516 522
517 if (!(fReducedXSTable->count(theKey))) << 523 if (!(theReducedXSTable->count(theKey)))
518 { 524 {
519 G4Exception("G4PenelopeBremsstrahlungFS: 525 G4Exception("G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
520 "em2013",FatalException,"Unable to retri 526 "em2013",FatalException,"Unable to retrieve the cross section table");
521 } 527 }
522 528
523 return fReducedXSTable->find(theKey)->second << 529 return theReducedXSTable->find(theKey)->second;
524 } 530 }
525 531
526 //....oooOO0OOooo........oooOO0OOooo........oo 532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
527 533
528 void G4PenelopeBremsstrahlungFS::InitializeEne 534 void G4PenelopeBremsstrahlungFS::InitializeEnergySampling(const G4Material* material,
529 G4double cut) 535 G4double cut)
530 { 536 {
531 if (fVerbosity > 2) 537 if (fVerbosity > 2)
532 G4cout << "Entering in G4PenelopeBremsstra 538 G4cout << "Entering in G4PenelopeBremsstrahlungFS::InitializeEnergySampling() for " <<
533 material->GetName() << G4endl; 539 material->GetName() << G4endl;
534 540
535 //This method should be accessed by the mast 541 //This method should be accessed by the master only
536 std::pair<const G4Material*,G4double> theKey 542 std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
537 543
538 G4PhysicsTable* thePhysicsTable = new G4Phys 544 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
539 // the table will contain 57 G4PhysicsFreeVe 545 // the table will contain 57 G4PhysicsFreeVectors with different
540 // values of E. 546 // values of E.
541 G4PhysicsFreeVector* thePBvec = new G4Physic << 547 G4PhysicsFreeVector* thePBvec = new G4PhysicsFreeVector(nBinsE);
542 548
543 //I reserve space of the vectors. 549 //I reserve space of the vectors.
544 for (std::size_t i=0;i<fNBinsE;++i) << 550 for (size_t i=0;i<nBinsE;i++)
545 thePhysicsTable->push_back(new G4PhysicsFr << 551 thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsX));
546 552
547 //Retrieve the table. Must already exist at 553 //Retrieve the table. Must already exist at this point, because this
548 //method is invoked by GetScaledXSTable() 554 //method is invoked by GetScaledXSTable()
549 if (!(fReducedXSTable->count(theKey))) << 555 if (!(theReducedXSTable->count(theKey)))
550 G4Exception("G4PenelopeBremsstrahlungFS::I 556 G4Exception("G4PenelopeBremsstrahlungFS::InitializeEnergySampling()",
551 "em2013",FatalException,"Unable to retriev 557 "em2013",FatalException,"Unable to retrieve the cross section table");
552 G4PhysicsTable* theTableReduced = fReducedXS << 558 G4PhysicsTable* theTableReduced = theReducedXSTable->find(theKey)->second;
553 559
554 for (std::size_t ie=0;ie<fNBinsE;++ie) << 560 for (size_t ie=0;ie<nBinsE;ie++)
555 { 561 {
556 G4PhysicsFreeVector* theVec = 562 G4PhysicsFreeVector* theVec =
557 (G4PhysicsFreeVector*) ((*thePhysicsTable)[i 563 (G4PhysicsFreeVector*) ((*thePhysicsTable)[ie]);
558 //Fill the table 564 //Fill the table
559 G4double value = 0; //first value 565 G4double value = 0; //first value
560 theVec->PutValues(0,theXGrid[0],value); << 566 theVec->PutValue(0,theXGrid[0],value);
561 for (std::size_t ix=1;ix<fNBinsX;++ix) << 567 for (size_t ix=1;ix<nBinsX;ix++)
562 { 568 {
563 //Here calculate the cumulative distributi 569 //Here calculate the cumulative distribution
564 // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx' 570 // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx'
565 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVe 571 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableReduced)[ix-1];
566 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVe 572 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
567 573
568 G4double x1=std::max(theXGrid[ix-1],1.0e-3 574 G4double x1=std::max(theXGrid[ix-1],1.0e-35);
569 //Remember: the table fReducedXSTable has << 575 //Remember: the table theReducedXSTable has a fake first point in energy
570 //so, it contains one more bin than fNBins << 576 //so, it contains one more bin than nBinsE.
571 G4double y1=G4Exp((*v1)[ie+1]); 577 G4double y1=G4Exp((*v1)[ie+1]);
572 G4double x2=std::max(theXGrid[ix],1.0e-35) 578 G4double x2=std::max(theXGrid[ix],1.0e-35);
573 G4double y2=G4Exp((*v2)[ie+1]); 579 G4double y2=G4Exp((*v2)[ie+1]);
574 G4double B = (y2-y1)/(x2-x1); 580 G4double B = (y2-y1)/(x2-x1);
575 G4double A = y1-B*x1; 581 G4double A = y1-B*x1;
576 G4double dS = A*G4Log(x2/x1)+B*(x2-x1); << 582 G4double dS = A*std::log(x2/x1)+B*(x2-x1);
577 value += dS; 583 value += dS;
578 theVec->PutValues(ix,theXGrid[ix],value); << 584 theVec->PutValue(ix,theXGrid[ix],value);
579 } 585 }
>> 586
580 //fill the PB vector 587 //fill the PB vector
581 G4double xc = cut/theEGrid[ie]; 588 G4double xc = cut/theEGrid[ie];
582 //Fill a temp data vector 589 //Fill a temp data vector
583 G4double* tempData = new G4double[fNBins << 590 G4double* tempData = new G4double[nBinsX];
584 for (std::size_t ix=0;ix<fNBinsX;++ix) << 591 for (size_t ix=0;ix<nBinsX;ix++)
585 { 592 {
586 G4PhysicsFreeVector* vv = (G4PhysicsFreeVe 593 G4PhysicsFreeVector* vv = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
587 tempData[ix] = G4Exp((*vv)[ie+1]); 594 tempData[ix] = G4Exp((*vv)[ie+1]);
588 } 595 }
589 G4double pbval = (xc<=1) ? 596 G4double pbval = (xc<=1) ?
590 GetMomentumIntegral(tempData,xc,-1) : 597 GetMomentumIntegral(tempData,xc,-1) :
591 GetMomentumIntegral(tempData,1.0,-1); 598 GetMomentumIntegral(tempData,1.0,-1);
592 thePBvec->PutValues(ie,theEGrid[ie],pbva << 599 thePBvec->PutValue(ie,theEGrid[ie],pbval);
593 delete[] tempData; 600 delete[] tempData;
594 } 601 }
595 602
596 fSamplingTable->insert(std::make_pair(theKey << 603 theSamplingTable->insert(std::make_pair(theKey,thePhysicsTable));
597 fPBcut->insert(std::make_pair(theKey,thePBve << 604 thePBcut->insert(std::make_pair(theKey,thePBvec));
598 return; 605 return;
599 } 606 }
600 607
601 //....oooOO0OOooo........oooOO0OOooo........oo 608 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
602 609
603 G4double G4PenelopeBremsstrahlungFS::SampleGam 610 G4double G4PenelopeBremsstrahlungFS::SampleGammaEnergy(G4double energy,const G4Material* mat,
604 const G4double cut) const 611 const G4double cut) const
605 { 612 {
606 std::pair<const G4Material*,G4double> theKey 613 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
607 if (!(fSamplingTable->count(theKey)) || !(fP << 614 if (!(theSamplingTable->count(theKey)) || !(thePBcut->count(theKey)))
608 { 615 {
609 G4ExceptionDescription ed; 616 G4ExceptionDescription ed;
610 ed << "Unable to retrieve the SamplingTa 617 ed << "Unable to retrieve the SamplingTable: " <<
611 fSamplingTable->count(theKey) << " " << << 618 theSamplingTable->count(theKey) << " " <<
612 fPBcut->count(theKey) << G4endl; << 619 thePBcut->count(theKey) << G4endl;
613 G4Exception("G4PenelopeBremsstrahlungFS: 620 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
614 "em2014",FatalException,ed); 621 "em2014",FatalException,ed);
615 } 622 }
616 const G4PhysicsTable* theTableInte = fSampli << 623
617 const G4PhysicsTable* theTableRed = fReduced << 624
>> 625 const G4PhysicsTable* theTableInte = theSamplingTable->find(theKey)->second;
>> 626 const G4PhysicsTable* theTableRed = theReducedXSTable->find(theKey)->second;
618 627
619 //Find the energy bin using bi-partition 628 //Find the energy bin using bi-partition
620 std::size_t eBin = 0; << 629 size_t eBin = 0;
621 G4bool firstOrLastBin = false; 630 G4bool firstOrLastBin = false;
622 631
623 if (energy < theEGrid[0]) //below first bin 632 if (energy < theEGrid[0]) //below first bin
624 { 633 {
625 eBin = 0; 634 eBin = 0;
626 firstOrLastBin = true; 635 firstOrLastBin = true;
627 } 636 }
628 else if (energy > theEGrid[fNBinsE-1]) //aft << 637 else if (energy > theEGrid[nBinsE-1]) //after last bin
629 { 638 {
630 eBin = fNBinsE-1; << 639 eBin = nBinsE-1;
631 firstOrLastBin = true; 640 firstOrLastBin = true;
632 } 641 }
633 else 642 else
634 { 643 {
635 std::size_t i=0; << 644 size_t i=0;
636 std::size_t j=fNBinsE-1; << 645 size_t j=nBinsE-1;
637 while ((j-i)>1) 646 while ((j-i)>1)
638 { 647 {
639 std::size_t k = (i+j)/2; << 648 size_t k = (i+j)/2;
640 if (energy > theEGrid[k]) 649 if (energy > theEGrid[k])
641 i = k; 650 i = k;
642 else 651 else
643 j = k; 652 j = k;
644 } 653 }
645 eBin = i; 654 eBin = i;
646 } 655 }
647 656
648 //Get the appropriate physics vector 657 //Get the appropriate physics vector
649 const G4PhysicsFreeVector* theVec1 = (G4Phys 658 const G4PhysicsFreeVector* theVec1 = (G4PhysicsFreeVector*) (*theTableInte)[eBin];
650 659
651 //Use a "temporary" vector which contains th 660 //Use a "temporary" vector which contains the linear interpolation of the x spectra
652 //in energy. The temporary vector is thread- 661 //in energy. The temporary vector is thread-local, so that there is no conflict.
653 //This is achieved via G4Cache. The theTempV 662 //This is achieved via G4Cache. The theTempVect is allocated only once per thread
654 //(member variable), but it is overwritten a 663 //(member variable), but it is overwritten at every call of this method
655 //(because the interpolation factors change! 664 //(because the interpolation factors change!)
656 G4PhysicsFreeVector* theTempVec = fCache.Get 665 G4PhysicsFreeVector* theTempVec = fCache.Get();
657 if (!theTempVec) //First time this thread ge 666 if (!theTempVec) //First time this thread gets the cache
658 { 667 {
659 theTempVec = new G4PhysicsFreeVector(fNB << 668 theTempVec = new G4PhysicsFreeVector(nBinsX);
660 fCache.Put(theTempVec); 669 fCache.Put(theTempVec);
661 // The G4AutoDelete takes care here to c 670 // The G4AutoDelete takes care here to clean up the vectors
662 G4AutoDelete::Register(theTempVec); 671 G4AutoDelete::Register(theTempVec);
663 if (fVerbosity > 4) 672 if (fVerbosity > 4)
664 G4cout << "Creating new instance of G4Physic 673 G4cout << "Creating new instance of G4PhysicsFreeVector() on the worker" << G4endl;
665 } 674 }
666 675
667 //theTempVect is allocated only once (member 676 //theTempVect is allocated only once (member variable), but it is overwritten at
668 //every call of this method (because the int 677 //every call of this method (because the interpolation factors change!)
669 if (!firstOrLastBin) 678 if (!firstOrLastBin)
670 { 679 {
671 const G4PhysicsFreeVector* theVec2 = (G4 680 const G4PhysicsFreeVector* theVec2 = (G4PhysicsFreeVector*) (*theTableInte)[eBin+1];
672 for (std::size_t iloop=0;iloop<fNBinsX;+ << 681 for (size_t iloop=0;iloop<nBinsX;iloop++)
673 { 682 {
674 G4double val = (*theVec1)[iloop]+(((*theVe 683 G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
675 (energy-theEGrid[eBin])/(theEGrid[eBin+1 684 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
676 theTempVec->PutValues(iloop,theXGrid[iloop << 685 theTempVec->PutValue(iloop,theXGrid[iloop],val);
677 } 686 }
678 } 687 }
679 else //first or last bin, no interpolation 688 else //first or last bin, no interpolation
680 { 689 {
681 for (std::size_t iloop=0;iloop<fNBinsX;+ << 690 for (size_t iloop=0;iloop<nBinsX;iloop++)
682 theTempVec->PutValues(iloop,theXGrid[iloop], << 691 theTempVec->PutValue(iloop,theXGrid[iloop],(*theVec1)[iloop]);
683 } 692 }
684 693
685 //Start the game 694 //Start the game
686 G4double pbcut = (*(fPBcut->find(theKey)->se << 695 G4double pbcut = (*(thePBcut->find(theKey)->second))[eBin];
687 696
688 if (!firstOrLastBin) //linear interpolation 697 if (!firstOrLastBin) //linear interpolation on pbcut as well
689 { 698 {
690 pbcut = (*(fPBcut->find(theKey)->second) << 699 pbcut = (*(thePBcut->find(theKey)->second))[eBin] +
691 ((*(fPBcut->find(theKey)->second))[eBin+1]-( << 700 ((*(thePBcut->find(theKey)->second))[eBin+1]-(*(thePBcut->find(theKey)->second))[eBin])*
692 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-th 701 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
693 } 702 }
694 703
695 G4double pCumulative = (*theTempVec)[fNBinsX << 704 G4double pCumulative = (*theTempVec)[nBinsX-1]; //last value
696 705
697 G4double eGamma = 0; 706 G4double eGamma = 0;
698 G4int nIterations = 0; 707 G4int nIterations = 0;
699 do 708 do
700 { 709 {
701 G4double pt = pbcut + G4UniformRand()*(p 710 G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut);
702 nIterations++; 711 nIterations++;
703 712
704 //find where it is 713 //find where it is
705 std::size_t ibin = 0; << 714 size_t ibin = 0;
706 if (pt < (*theTempVec)[0]) 715 if (pt < (*theTempVec)[0])
707 ibin = 0; 716 ibin = 0;
708 else if (pt > (*theTempVec)[fNBinsX-1]) << 717 else if (pt > (*theTempVec)[nBinsX-1])
709 { 718 {
710 //We observed problems due to numerical ro 719 //We observed problems due to numerical rounding here (STT).
711 //delta here is a tiny positive number 720 //delta here is a tiny positive number
712 G4double delta = pt-(*theTempVec)[fNBinsX- << 721 G4double delta = pt-(*theTempVec)[nBinsX-1];
713 if (delta < pt*1e-10) // very small! Numer 722 if (delta < pt*1e-10) // very small! Numerical rounding only
714 { 723 {
715 ibin = fNBinsX-2; << 724 ibin = nBinsX-2;
716 G4ExceptionDescription ed; 725 G4ExceptionDescription ed;
717 ed << "Found that (pt > (*theTempVec)[ << 726 ed << "Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
718 " , (*theTempVec)[fNBinsX-1] = " << (*theT << 727 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] << " and delta = " <<
719 (pt-(*theTempVec)[fNBinsX-1]) << G4endl; << 728 (pt-(*theTempVec)[nBinsX-1]) << G4endl;
720 ed << "Possible symptom of problem wit 729 ed << "Possible symptom of problem with numerical precision" << G4endl;
721 G4Exception("G4PenelopeBremsstrahlungF 730 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
722 "em2015",JustWarning,ed); 731 "em2015",JustWarning,ed);
723 } 732 }
724 else //real problem 733 else //real problem
725 { 734 {
726 G4ExceptionDescription ed; 735 G4ExceptionDescription ed;
727 ed << "Crash at (pt > (*theTempVec)[fN << 736 ed << "Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
728 " , (*theTempVec)[fNBinsX-1]=" << (*theTem << 737 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] << " and nBinsX = " <<
729 fNBinsX << G4endl; << 738 nBinsX << G4endl;
730 ed << "Material: " << mat->GetName() < 739 ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" <<
731 G4endl; 740 G4endl;
732 G4Exception("G4PenelopeBremsstrahlungF 741 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
733 "em2015",FatalException,ed); 742 "em2015",FatalException,ed);
734 } 743 }
735 } 744 }
736 else 745 else
737 { 746 {
738 std::size_t i=0; << 747 size_t i=0;
739 std::size_t j=fNBinsX-1; << 748 size_t j=nBinsX-1;
740 while ((j-i)>1) 749 while ((j-i)>1)
741 { 750 {
742 std::size_t k = (i+j)/2; << 751 size_t k = (i+j)/2;
743 if (pt > (*theTempVec)[k]) 752 if (pt > (*theTempVec)[k])
744 i = k; 753 i = k;
745 else 754 else
746 j = k; 755 j = k;
747 } 756 }
748 ibin = i; 757 ibin = i;
749 } 758 }
750 759
751 G4double w1 = theXGrid[ibin]; 760 G4double w1 = theXGrid[ibin];
752 G4double w2 = theXGrid[ibin+1]; 761 G4double w2 = theXGrid[ibin+1];
753 762
754 const G4PhysicsFreeVector* v1 = (G4Physi 763 const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin];
755 const G4PhysicsFreeVector* v2 = (G4Physi 764 const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1];
756 //Remember: the table fReducedXSTable ha << 765 //Remember: the table theReducedXSTable has a fake first point in energy
757 //so, it contains one more bin than fNBi << 766 //so, it contains one more bin than nBinsE.
758 G4double pdf1 = G4Exp((*v1)[eBin+1]); 767 G4double pdf1 = G4Exp((*v1)[eBin+1]);
759 G4double pdf2 = G4Exp((*v2)[eBin+1]); 768 G4double pdf2 = G4Exp((*v2)[eBin+1]);
760 G4double deltaW = w2-w1; 769 G4double deltaW = w2-w1;
761 G4double dpdfb = pdf2-pdf1; 770 G4double dpdfb = pdf2-pdf1;
762 G4double B = dpdfb/deltaW; 771 G4double B = dpdfb/deltaW;
763 G4double A = pdf1-B*w1; 772 G4double A = pdf1-B*w1;
764 //I already made an interpolation in ene 773 //I already made an interpolation in energy, so I can use the actual value for the
765 //calculation of the wbcut, instead of t 774 //calculation of the wbcut, instead of the grid values (except for the last bin)
766 G4double wbcut = (cut < energy) ? cut/e 775 G4double wbcut = (cut < energy) ? cut/energy : 1.0;
767 if (firstOrLastBin) //this is an particu 776 if (firstOrLastBin) //this is an particular case: no interpolation available
768 wbcut = (cut < theEGrid[eBin]) ? cut/theEGr 777 wbcut = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
769 778
770 if (w1 < wbcut) 779 if (w1 < wbcut)
771 w1 = wbcut; 780 w1 = wbcut;
772 if (w2 < w1) 781 if (w2 < w1)
773 { 782 {
774 //This configuration can happen if initial 783 //This configuration can happen if initially wbcut > w2 > w1. Due to the previous
775 //statement, (w1 = wbcut), it becomes wbcu 784 //statement, (w1 = wbcut), it becomes wbcut = w1 > w2. In this case, it is not a
776 //real problem. It becomes a problem if w2 785 //real problem. It becomes a problem if w2 < w1 before the w1 = wbcut statement. Issue
777 //a warning only in this specific case. 786 //a warning only in this specific case.
778 if (w2 > wbcut) 787 if (w2 > wbcut)
779 { 788 {
780 G4ExceptionDescription ed; 789 G4ExceptionDescription ed;
781 ed << "Warning in G4PenelopeBremsstrah 790 ed << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl;
782 ed << "Conflicting end-point values: w 791 ed << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl;
783 ed << "wbcut = " << wbcut << " energy= 792 ed << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl;
784 ed << "cut = " << cut/keV << " keV" << 793 ed << "cut = " << cut/keV << " keV" << G4endl;
785 G4Exception("G4PenelopeBremsstrahlungF 794 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()","em2015",
786 JustWarning,ed); 795 JustWarning,ed);
787 } 796 }
788 return w1*energy; 797 return w1*energy;
789 } 798 }
790 799
791 G4double pmax = std::max(A+B*w1,A+B*w2); 800 G4double pmax = std::max(A+B*w1,A+B*w2);
792 G4bool loopAgain = false; 801 G4bool loopAgain = false;
793 do 802 do
794 { 803 {
795 loopAgain = false; 804 loopAgain = false;
796 eGamma = w1* std::pow((w2/w1),G4UniformRan 805 eGamma = w1* std::pow((w2/w1),G4UniformRand());
797 if (G4UniformRand()*pmax > (A+B*eGamma)) 806 if (G4UniformRand()*pmax > (A+B*eGamma))
798 loopAgain = true; 807 loopAgain = true;
799 }while(loopAgain); 808 }while(loopAgain);
800 eGamma *= energy; 809 eGamma *= energy;
801 if (nIterations > 100) //protection agai 810 if (nIterations > 100) //protection against infinite loops
802 return eGamma; 811 return eGamma;
803 }while(eGamma < cut); //repeat if sampled 812 }while(eGamma < cut); //repeat if sampled sub-cut!
804 813
805 return eGamma; 814 return eGamma;
806 } 815 }
807 816