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