| 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 76988 2013-11-20 09:54:40Z 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 "G4Exp.hh" <<
47 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
49 48
50 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 50
52 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstr << 51 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstrahlungFS(G4int verbosity) :
53 fReducedXSTable(nullptr),fEffectiveZSq(nullp << 52 theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
54 fPBcut(nullptr),fVerbosity(verbosity) << 53 thePBcut(0),fVerbosity(verbosity)
55 { 54 {
56 fCache.Put(0); 55 fCache.Put(0);
57 G4double tempvector[fNBinsX] = << 56 G4double tempvector[nBinsX] =
58 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15 57 {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 58 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 59 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 60 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
62 61
63 for (std::size_t ix=0;ix<fNBinsX;++ix) << 62 for (size_t ix=0;ix<nBinsX;ix++)
64 theXGrid[ix] = tempvector[ix]; 63 theXGrid[ix] = tempvector[ix];
65 64
66 for (std::size_t i=0;i<fNBinsE;++i) << 65 for (size_t i=0;i<nBinsE;i++)
67 theEGrid[i] = 0.; 66 theEGrid[i] = 0.;
68 67
69 fElementData = new std::map<G4int,G4DataVect << 68 theElementData = new std::map<G4int,G4DataVector*>;
70 } 69 }
71 70
72 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 72
74 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsst 73 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsstrahlungFS()
75 { 74 {
76 ClearTables(); << 75 ClearTables();
77 << 76
78 //The G4Physics*Vector pointers contained in << 77 //The G4Physics*Vector pointers contained in the fCache are automatically deleted by
79 //the G4AutoDelete so there is no need to ta << 78 //the G4Allocator, so there is no need to take care of them manually
80 << 79
81 //Clear manually fElementData << 80 //Clear manually theElementData
82 if (fElementData) << 81 if (theElementData)
83 { << 82 {
84 for (auto& item : (*fElementData)) << 83 std::map<G4int,G4DataVector*>::iterator i;
85 delete item.second; << 84 for (i=theElementData->begin(); i != theElementData->end(); i++)
86 delete fElementData; << 85 delete i->second;
87 fElementData = nullptr; << 86 delete theElementData;
>> 87 theElementData = 0;
88 } 88 }
>> 89
89 } 90 }
90 91
91 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
92 93
93 94
94 void G4PenelopeBremsstrahlungFS::ClearTables(G 95 void G4PenelopeBremsstrahlungFS::ClearTables(G4bool isMaster)
95 { << 96 {
96 //Just to check 97 //Just to check
97 if (!isMaster) 98 if (!isMaster)
98 G4Exception("G4PenelopeBremsstrahlungFS::C 99 G4Exception("G4PenelopeBremsstrahlungFS::ClearTables()",
99 "em0100",FatalException,"Worker thread in << 100 "em0100",FatalException,"Worker thread in this method");
>> 101
>> 102 std::map< std::pair<const G4Material*,G4double> ,G4PhysicsTable*>::iterator j;
100 103
101 if (fReducedXSTable) << 104 if (theReducedXSTable)
102 { 105 {
103 for (auto& item : (*fReducedXSTable)) << 106 for (j=theReducedXSTable->begin(); j != theReducedXSTable->end(); j++)
104 { 107 {
105 G4PhysicsTable* tab = item.second; << 108 G4PhysicsTable* tab = j->second;
106 tab->clearAndDestroy(); << 109 //tab->clearAndDestroy();
107 delete tab; << 110 delete tab;
108 } 111 }
109 fReducedXSTable->clear(); << 112 delete theReducedXSTable;
110 delete fReducedXSTable; << 113 theReducedXSTable = 0;
111 fReducedXSTable = nullptr; <<
112 } 114 }
113 115
114 if (fSamplingTable) << 116 if (theSamplingTable)
115 { 117 {
116 for (auto& item : (*fSamplingTable)) << 118 for (j=theSamplingTable->begin(); j != theSamplingTable->end(); j++)
117 { 119 {
118 G4PhysicsTable* tab = item.second; << 120 G4PhysicsTable* tab = j->second;
119 tab->clearAndDestroy(); << 121 // tab->clearAndDestroy();
120 delete tab; 122 delete tab;
121 } 123 }
122 fSamplingTable->clear(); << 124 delete theSamplingTable;
123 delete fSamplingTable; << 125 theSamplingTable = 0;
124 fSamplingTable = nullptr; << 126 }
125 } << 127 if (thePBcut)
126 if (fPBcut) <<
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 = 0;
135 } 136 }
136 137
137 if (fEffectiveZSq) << 138
>> 139 if (theEffectiveZSq)
138 { 140 {
139 delete fEffectiveZSq; << 141 delete theEffectiveZSq;
140 fEffectiveZSq = nullptr; << 142 theEffectiveZSq = 0;
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=std::exp((*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=std::exp((*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] = std::exp((*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);
>> 669 //The G4Physics*Vector pointers are automatically deleted by the G4Allocator,
>> 670 //so there is no need to take care of it manually
660 fCache.Put(theTempVec); 671 fCache.Put(theTempVec);
661 // The G4AutoDelete takes care here to c <<
662 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; <<
699 do 707 do
700 { << 708 {
701 G4double pt = pbcut + G4UniformRand()*(p 709 G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut);
702 nIterations++; <<
703 710
704 //find where it is 711 //find where it is
705 std::size_t ibin = 0; << 712 size_t ibin = 0;
706 if (pt < (*theTempVec)[0]) 713 if (pt < (*theTempVec)[0])
707 ibin = 0; 714 ibin = 0;
708 else if (pt > (*theTempVec)[fNBinsX-1]) << 715 else if (pt > (*theTempVec)[nBinsX-1])
709 { 716 {
710 //We observed problems due to numerical ro << 717 //We observed problems due to numerical rounding here (STT).
711 //delta here is a tiny positive number 718 //delta here is a tiny positive number
712 G4double delta = pt-(*theTempVec)[fNBinsX- << 719 G4double delta = pt-(*theTempVec)[nBinsX-1];
713 if (delta < pt*1e-10) // very small! Numer 720 if (delta < pt*1e-10) // very small! Numerical rounding only
714 { 721 {
715 ibin = fNBinsX-2; << 722 ibin = nBinsX-2;
716 G4ExceptionDescription ed; 723 G4ExceptionDescription ed;
717 ed << "Found that (pt > (*theTempVec)[ << 724 ed << "Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
718 " , (*theTempVec)[fNBinsX-1] = " << (*theT << 725 " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] << " and delta = " <<
719 (pt-(*theTempVec)[fNBinsX-1]) << G4endl; << 726 (pt-(*theTempVec)[nBinsX-1]) << G4endl;
720 ed << "Possible symptom of problem wit 727 ed << "Possible symptom of problem with numerical precision" << G4endl;
721 G4Exception("G4PenelopeBremsstrahlungF 728 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
722 "em2015",JustWarning,ed); 729 "em2015",JustWarning,ed);
723 } 730 }
724 else //real problem 731 else //real problem
725 { 732 {
726 G4ExceptionDescription ed; 733 G4ExceptionDescription ed;
727 ed << "Crash at (pt > (*theTempVec)[fN << 734 ed << "Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
728 " , (*theTempVec)[fNBinsX-1]=" << (*theTem << 735 " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] << " and nBinsX = " <<
729 fNBinsX << G4endl; << 736 nBinsX << G4endl;
730 ed << "Material: " << mat->GetName() < << 737 ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" <<
731 G4endl; << 738 G4endl;
732 G4Exception("G4PenelopeBremsstrahlungF 739 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
733 "em2015",FatalException,ed); << 740 "em2015",FatalException,ed);
734 } 741 }
735 } 742 }
736 else 743 else
737 { 744 {
738 std::size_t i=0; << 745 size_t i=0;
739 std::size_t j=fNBinsX-1; << 746 size_t j=nBinsX-1;
740 while ((j-i)>1) 747 while ((j-i)>1)
741 { 748 {
742 std::size_t k = (i+j)/2; << 749 size_t k = (i+j)/2;
743 if (pt > (*theTempVec)[k]) 750 if (pt > (*theTempVec)[k])
744 i = k; 751 i = k;
745 else 752 else
746 j = k; 753 j = k;
747 } 754 }
748 ibin = i; 755 ibin = i;
749 } 756 }
750 << 757
751 G4double w1 = theXGrid[ibin]; 758 G4double w1 = theXGrid[ibin];
752 G4double w2 = theXGrid[ibin+1]; << 759 G4double w2 = theXGrid[ibin+1];
753 760
754 const G4PhysicsFreeVector* v1 = (G4Physi 761 const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin];
755 const G4PhysicsFreeVector* v2 = (G4Physi << 762 const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1];
756 //Remember: the table fReducedXSTable ha << 763 //Remember: the table theReducedXSTable has a fake first point in energy
757 //so, it contains one more bin than fNBi << 764 //so, it contains one more bin than nBinsE.
758 G4double pdf1 = G4Exp((*v1)[eBin+1]); << 765 G4double pdf1 = std::exp((*v1)[eBin+1]);
759 G4double pdf2 = G4Exp((*v2)[eBin+1]); << 766 G4double pdf2 = std::exp((*v2)[eBin+1]);
760 G4double deltaW = w2-w1; 767 G4double deltaW = w2-w1;
761 G4double dpdfb = pdf2-pdf1; 768 G4double dpdfb = pdf2-pdf1;
762 G4double B = dpdfb/deltaW; 769 G4double B = dpdfb/deltaW;
763 G4double A = pdf1-B*w1; 770 G4double A = pdf1-B*w1;
764 //I already made an interpolation in ene << 771 //I already made an interpolation in energy, so I can use the actual value for the
765 //calculation of the wbcut, instead of t 772 //calculation of the wbcut, instead of the grid values (except for the last bin)
766 G4double wbcut = (cut < energy) ? cut/e 773 G4double wbcut = (cut < energy) ? cut/energy : 1.0;
767 if (firstOrLastBin) //this is an particu 774 if (firstOrLastBin) //this is an particular case: no interpolation available
768 wbcut = (cut < theEGrid[eBin]) ? cut/theEGr 775 wbcut = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
769 << 776
770 if (w1 < wbcut) 777 if (w1 < wbcut)
771 w1 = wbcut; 778 w1 = wbcut;
772 if (w2 < w1) 779 if (w2 < w1)
773 { 780 {
774 //This configuration can happen if initial << 781 //This configuration can happen if initially wbcut > w2 > w1. Due to the previous
775 //statement, (w1 = wbcut), it becomes wbcu << 782 //statement, (w1 = wbcut), it becomes wbcut = w1 > w2. In this case, it is not a
776 //real problem. It becomes a problem if w2 << 783 //real problem. It becomes a problem if w2 < w1 before the w1 = wbcut statement. Issue
777 //a warning only in this specific case. 784 //a warning only in this specific case.
778 if (w2 > wbcut) 785 if (w2 > wbcut)
779 { 786 {
780 G4ExceptionDescription ed; 787 G4ExceptionDescription ed;
781 ed << "Warning in G4PenelopeBremsstrah 788 ed << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl;
782 ed << "Conflicting end-point values: w 789 ed << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl;
783 ed << "wbcut = " << wbcut << " energy= 790 ed << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl;
784 ed << "cut = " << cut/keV << " keV" << 791 ed << "cut = " << cut/keV << " keV" << G4endl;
785 G4Exception("G4PenelopeBremsstrahlungF 792 G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()","em2015",
786 JustWarning,ed); 793 JustWarning,ed);
787 } 794 }
788 return w1*energy; 795 return w1*energy;
789 } 796 }
790 << 797
791 G4double pmax = std::max(A+B*w1,A+B*w2); 798 G4double pmax = std::max(A+B*w1,A+B*w2);
792 G4bool loopAgain = false; << 799 G4bool loopAgain = false;
793 do 800 do
794 { 801 {
795 loopAgain = false; 802 loopAgain = false;
796 eGamma = w1* std::pow((w2/w1),G4UniformRan 803 eGamma = w1* std::pow((w2/w1),G4UniformRand());
797 if (G4UniformRand()*pmax > (A+B*eGamma)) 804 if (G4UniformRand()*pmax > (A+B*eGamma))
798 loopAgain = true; << 805 loopAgain = true;
799 }while(loopAgain); << 806 }while(loopAgain);
800 eGamma *= energy; 807 eGamma *= energy;
801 if (nIterations > 100) //protection agai << 808 }while(eGamma < cut); //repeat if sampled sub-cut!
802 return eGamma; <<
803 }while(eGamma < cut); //repeat if sampled <<
804 809
805 return eGamma; 810 return eGamma;
806 } 811 }
>> 812
>> 813
>> 814
>> 815
807 816