| 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: G4PenelopeRayleighModel.cc 83584 2014-09-02 08:45:37Z gcosmo $
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 03 Dec 2009 L Pandola First implementa 32 // 03 Dec 2009 L Pandola First implementation
32 // 25 May 2011 L.Pandola Renamed (make v2 33 // 25 May 2011 L.Pandola Renamed (make v2008 as default Penelope)
33 // 19 Sep 2013 L.Pandola Migration to MT 34 // 19 Sep 2013 L.Pandola Migration to MT
34 // 35 //
35 36
36 #include "G4PenelopeRayleighModel.hh" 37 #include "G4PenelopeRayleighModel.hh"
37 #include "G4PhysicalConstants.hh" 38 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh"
39 #include "G4PenelopeSamplingData.hh" 40 #include "G4PenelopeSamplingData.hh"
40 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh"
41 #include "G4MaterialCutsCouple.hh" 42 #include "G4MaterialCutsCouple.hh"
42 #include "G4ProductionCutsTable.hh" 43 #include "G4ProductionCutsTable.hh"
43 #include "G4DynamicParticle.hh" 44 #include "G4DynamicParticle.hh"
44 #include "G4PhysicsTable.hh" 45 #include "G4PhysicsTable.hh"
45 #include "G4ElementTable.hh" 46 #include "G4ElementTable.hh"
46 #include "G4Element.hh" 47 #include "G4Element.hh"
47 #include "G4PhysicsFreeVector.hh" 48 #include "G4PhysicsFreeVector.hh"
48 #include "G4AutoLock.hh" 49 #include "G4AutoLock.hh"
49 #include "G4Exp.hh" <<
50 50
51 //....oooOO0OOooo........oooOO0OOooo........oo <<
52 <<
53 const G4int G4PenelopeRayleighModel::fMaxZ; <<
54 G4PhysicsFreeVector* G4PenelopeRayleighModel:: <<
55 G4PhysicsFreeVector* G4PenelopeRayleighModel:: <<
56 51
57 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 53
59 G4PenelopeRayleighModel::G4PenelopeRayleighMod 54 G4PenelopeRayleighModel::G4PenelopeRayleighModel(const G4ParticleDefinition* part,
60 const G4String& nam) 55 const G4String& nam)
61 :G4VEmModel(nam),fParticleChange(nullptr),fP << 56 :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
62 fLogFormFactorTable(nullptr),fPMaxTable(nul << 57 logAtomicCrossSection(0),atomicFormFactor(0),logFormFactorTable(0),
63 fIsInitialised(false),fLocalTable(false) << 58 pMaxTable(0),samplingTable(0),fLocalTable(false)
64 { 59 {
65 fIntrinsicLowEnergyLimit = 100.0*eV; 60 fIntrinsicLowEnergyLimit = 100.0*eV;
66 fIntrinsicHighEnergyLimit = 100.0*GeV; 61 fIntrinsicHighEnergyLimit = 100.0*GeV;
67 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 62 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
68 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 63 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
69 64
70 if (part) 65 if (part)
71 SetParticle(part); 66 SetParticle(part);
72 67
73 // 68 //
74 fVerboseLevel= 0; << 69 verboseLevel= 0;
75 // Verbosity scale: 70 // Verbosity scale:
76 // 0 = nothing << 71 // 0 = nothing
77 // 1 = warning for energy non-conservation << 72 // 1 = warning for energy non-conservation
78 // 2 = details of energy budget 73 // 2 = details of energy budget
79 // 3 = calculation of cross sections, file o 74 // 3 = calculation of cross sections, file openings, sampling of atoms
80 // 4 = entering in methods << 75 // 4 = entering in methods
81 76
82 //build the energy grid. It is the same for 77 //build the energy grid. It is the same for all materials
83 G4double logenergy = G4Log(fIntrinsicLowEner << 78 G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.);
84 G4double logmaxenergy = G4Log(1.5*fIntrinsic << 79 G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit);
85 //finer grid below 160 keV 80 //finer grid below 160 keV
86 G4double logtransitionenergy = G4Log(160*keV << 81 G4double logtransitionenergy = std::log(160*keV);
87 G4double logfactor1 = G4Log(10.)/250.; << 82 G4double logfactor1 = std::log(10.)/250.;
88 G4double logfactor2 = logfactor1*10; 83 G4double logfactor2 = logfactor1*10;
89 fLogEnergyGridPMax.push_back(logenergy); << 84 logEnergyGridPMax.push_back(logenergy);
90 do{ 85 do{
91 if (logenergy < logtransitionenergy) 86 if (logenergy < logtransitionenergy)
92 logenergy += logfactor1; 87 logenergy += logfactor1;
93 else 88 else
94 logenergy += logfactor2; 89 logenergy += logfactor2;
95 fLogEnergyGridPMax.push_back(logenergy); << 90 logEnergyGridPMax.push_back(logenergy);
96 }while (logenergy < logmaxenergy); 91 }while (logenergy < logmaxenergy);
97 } 92 }
98 93
99 //....oooOO0OOooo........oooOO0OOooo........oo 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 95
101 G4PenelopeRayleighModel::~G4PenelopeRayleighMo 96 G4PenelopeRayleighModel::~G4PenelopeRayleighModel()
102 { 97 {
103 if (IsMaster() || fLocalTable) 98 if (IsMaster() || fLocalTable)
104 { 99 {
105 << 100 std::map <G4int,G4PhysicsFreeVector*>::iterator i;
106 for(G4int i=0; i<=fMaxZ; ++i) << 101 if (logAtomicCrossSection)
107 { 102 {
108 if(fLogAtomicCrossSection[i]) << 103 for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
109 { << 104 if (i->second) delete i->second;
110 delete fLogAtomicCrossSection[i]; << 105
111 fLogAtomicCrossSection[i] = nullptr; << 106 delete logAtomicCrossSection;
112 } << 107 logAtomicCrossSection = 0;
113 if(fAtomicFormFactor[i]) << 108 }
114 { << 109 if (atomicFormFactor)
115 delete fAtomicFormFactor[i]; << 110 {
116 fAtomicFormFactor[i] = nullptr; << 111 for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++)
117 } << 112 if (i->second) delete i->second;
>> 113 delete atomicFormFactor;
>> 114 atomicFormFactor = 0;
118 } 115 }
>> 116
119 ClearTables(); 117 ClearTables();
120 } 118 }
121 } 119 }
122 120
123 //....oooOO0OOooo........oooOO0OOooo........oo 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 void G4PenelopeRayleighModel::ClearTables() 122 void G4PenelopeRayleighModel::ClearTables()
125 { 123 {
126 if (fLogFormFactorTable) << 124 /*
127 { << 125 if (!IsMaster())
128 for (auto& item : (*fLogFormFactorTable << 126 //Should not be here!
129 if (item.second) delete item.second; << 127 G4Exception("G4PenelopeRayleighModel::ClearTables()",
130 delete fLogFormFactorTable; << 128 "em0100",FatalException,"Worker thread in this method");
131 fLogFormFactorTable = nullptr; //zero e << 129 */
>> 130
>> 131 std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i;
>> 132
>> 133 if (logFormFactorTable)
>> 134 {
>> 135 for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++)
>> 136 if (i->second) delete i->second;
>> 137 delete logFormFactorTable;
>> 138 logFormFactorTable = 0; //zero explicitely
132 } 139 }
133 if (fPMaxTable) << 140
134 { << 141 if (pMaxTable)
135 for (auto& item : (*fPMaxTable)) << 142 {
136 if (item.second) delete item.second; << 143 for (i=pMaxTable->begin(); i != pMaxTable->end(); i++)
137 delete fPMaxTable; << 144 if (i->second) delete i->second;
138 fPMaxTable = nullptr; //zero explicitly << 145 delete pMaxTable;
>> 146 pMaxTable = 0; //zero explicitely
139 } 147 }
140 if (fSamplingTable) << 148
>> 149 std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii;
>> 150 if (samplingTable)
141 { 151 {
142 for (auto& item : (*fSamplingTable)) << 152 for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++)
143 if (item.second) delete item.second; << 153 if (ii->second) delete ii->second;
144 delete fSamplingTable; << 154 delete samplingTable;
145 fSamplingTable = nullptr; //zero explic << 155 samplingTable = 0; //zero explicitely
146 } << 156 }
>> 157
147 return; 158 return;
148 } 159 }
149 160
150 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 162
152 void G4PenelopeRayleighModel::Initialise(const 163 void G4PenelopeRayleighModel::Initialise(const G4ParticleDefinition* part,
153 const G4DataVector& ) 164 const G4DataVector& )
154 { 165 {
155 if (fVerboseLevel > 3) << 166 if (verboseLevel > 3)
156 G4cout << "Calling G4PenelopeRayleighModel 167 G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
157 168
158 SetParticle(part); 169 SetParticle(part);
159 170
160 //Only the master model creates/fills/destro 171 //Only the master model creates/fills/destroys the tables
161 if (IsMaster() && part == fParticle) << 172 if (IsMaster() && part == fParticle)
162 { 173 {
163 //clear tables depending on materials, n 174 //clear tables depending on materials, not the atomic ones
164 ClearTables(); 175 ClearTables();
165 176
166 if (fVerboseLevel > 3) << 177 if (verboseLevel > 3)
167 G4cout << "Calling G4PenelopeRayleighModel:: 178 G4cout << "Calling G4PenelopeRayleighModel::Initialise() [master]" << G4endl;
>> 179
>> 180 //create new tables
>> 181 //
>> 182 // logAtomicCrossSection and atomicFormFactor are created only once,
>> 183 // since they are never cleared
>> 184 if (!logAtomicCrossSection)
>> 185 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
>> 186 if (!atomicFormFactor)
>> 187 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
>> 188
>> 189 if (!logFormFactorTable)
>> 190 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
>> 191 if (!pMaxTable)
>> 192 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
>> 193 if (!samplingTable)
>> 194 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
168 195
169 //create new tables <<
170 if (!fLogFormFactorTable) <<
171 fLogFormFactorTable = new std::map<const G4M <<
172 if (!fPMaxTable) <<
173 fPMaxTable = new std::map<const G4Material*, <<
174 if (!fSamplingTable) <<
175 fSamplingTable = new std::map<const G4Materi <<
176 196
177 G4ProductionCutsTable* theCoupleTable = << 197 G4ProductionCutsTable* theCoupleTable =
178 G4ProductionCutsTable::GetProductionCutsTabl 198 G4ProductionCutsTable::GetProductionCutsTable();
179 << 199
180 for (G4int i=0;i<(G4int)theCoupleTable-> << 200 for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
181 { 201 {
182 const G4Material* material = << 202 const G4Material* material =
183 theCoupleTable->GetMaterialCutsCouple(i) 203 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
184 const G4ElementVector* theElementVector = 204 const G4ElementVector* theElementVector = material->GetElementVector();
185 << 205
186 for (std::size_t j=0;j<material->GetNumber << 206 for (size_t j=0;j<material->GetNumberOfElements();j++)
187 { 207 {
188 G4int iZ = theElementVector->at(j)->Ge << 208 G4int iZ = (G4int) theElementVector->at(j)->GetZ();
189 //read data files only in the master 209 //read data files only in the master
190 if (!fLogAtomicCrossSection[iZ]) << 210 if (!logAtomicCrossSection->count(iZ))
191 ReadDataFile(iZ); << 211 ReadDataFile(iZ);
192 } 212 }
193 213
194 //1) If the table has not been built for t 214 //1) If the table has not been built for the material, do it!
195 if (!fLogFormFactorTable->count(material)) << 215 if (!logFormFactorTable->count(material))
196 BuildFormFactorTable(material); 216 BuildFormFactorTable(material);
197 217
198 //2) retrieve or build the sampling table 218 //2) retrieve or build the sampling table
199 if (!(fSamplingTable->count(material))) << 219 if (!(samplingTable->count(material)))
200 InitializeSamplingAlgorithm(material); 220 InitializeSamplingAlgorithm(material);
201 221
202 //3) retrieve or build the pMax data 222 //3) retrieve or build the pMax data
203 if (!fPMaxTable->count(material)) << 223 if (!pMaxTable->count(material))
204 GetPMaxTable(material); << 224 GetPMaxTable(material);
205 } <<
206 225
207 if (fVerboseLevel > 1) { << 226 }
>> 227
>> 228 if (verboseLevel > 1) {
208 G4cout << "Penelope Rayleigh model v2008 is 229 G4cout << "Penelope Rayleigh model v2008 is initialized " << G4endl
209 << "Energy range: " 230 << "Energy range: "
210 << LowEnergyLimit() / keV << " keV - 231 << LowEnergyLimit() / keV << " keV - "
211 << HighEnergyLimit() / GeV << " GeV" 232 << HighEnergyLimit() / GeV << " GeV"
212 << G4endl; 233 << G4endl;
213 } 234 }
214 } 235 }
215 236
216 if(fIsInitialised) return; << 237 if(isInitialised) return;
217 fParticleChange = GetParticleChangeForGamma( 238 fParticleChange = GetParticleChangeForGamma();
218 fIsInitialised = true; << 239 isInitialised = true;
219 } 240 }
220 241
221 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 243
223 void G4PenelopeRayleighModel::InitialiseLocal( 244 void G4PenelopeRayleighModel::InitialiseLocal(const G4ParticleDefinition* part,
224 G4VEmModel *masterModel) 245 G4VEmModel *masterModel)
225 { 246 {
226 if (fVerboseLevel > 3) << 247 if (verboseLevel > 3)
227 G4cout << "Calling G4PenelopeRayleighMode 248 G4cout << "Calling G4PenelopeRayleighModel::InitialiseLocal()" << G4endl;
>> 249
228 // 250 //
229 //Check that particle matches: one might hav << 251 //Check that particle matches: one might have multiple master models (e.g.
230 //for e+ and e-). 252 //for e+ and e-).
231 // 253 //
232 if (part == fParticle) 254 if (part == fParticle)
233 { 255 {
234 //Get the const table pointers from the 256 //Get the const table pointers from the master to the workers
235 const G4PenelopeRayleighModel* theModel << 257 const G4PenelopeRayleighModel* theModel =
236 static_cast<G4PenelopeRayleighModel*> (maste 258 static_cast<G4PenelopeRayleighModel*> (masterModel);
237 << 259
238 //Copy pointers to the data tables 260 //Copy pointers to the data tables
239 for(G4int i=0; i<=fMaxZ; ++i) << 261 logAtomicCrossSection = theModel->logAtomicCrossSection;
240 { << 262 atomicFormFactor = theModel->atomicFormFactor;
241 fLogAtomicCrossSection[i] = theModel->fLog << 263 logFormFactorTable = theModel->logFormFactorTable;
242 fAtomicFormFactor[i] = theModel->fAtomicFo << 264 pMaxTable = theModel->pMaxTable;
243 } << 265 samplingTable = theModel->samplingTable;
244 fLogFormFactorTable = theModel->fLogForm << 266
245 fPMaxTable = theModel->fPMaxTable; <<
246 fSamplingTable = theModel->fSamplingTabl <<
247 <<
248 //copy the G4DataVector with the grid 267 //copy the G4DataVector with the grid
249 fLogQSquareGrid = theModel->fLogQSquareG << 268 logQSquareGrid = theModel->logQSquareGrid;
250 << 269
251 //Same verbosity for all workers, as the 270 //Same verbosity for all workers, as the master
252 fVerboseLevel = theModel->fVerboseLevel; << 271 verboseLevel = theModel->verboseLevel;
253 } 272 }
254 273
255 return; 274 return;
256 } 275 }
257 276
258 277
259 //....oooOO0OOooo........oooOO0OOooo........oo 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260 namespace { G4Mutex PenelopeRayleighModelMute 279 namespace { G4Mutex PenelopeRayleighModelMutex = G4MUTEX_INITIALIZER; }
261 G4double G4PenelopeRayleighModel::ComputeCross 280 G4double G4PenelopeRayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
262 G4double energy, 281 G4double energy,
263 G4double Z, 282 G4double Z,
264 G4double, 283 G4double,
265 G4double, 284 G4double,
266 G4double) 285 G4double)
267 { 286 {
268 // Cross section of Rayleigh scattering in P << 287 // Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97
269 // tabulation, Cuellen et al. (1997), with n << 288 // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel
270 // et al. J. Phys. Chem. Ref. Data 4 (1975) 289 // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.
271 290
272 if (fVerboseLevel > 3) << 291 if (verboseLevel > 3)
273 G4cout << "Calling CrossSectionPerAtom() o 292 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl;
274 << 293
275 G4int iZ = G4int(Z); << 294 G4int iZ = (G4int) Z;
276 << 295
277 if (!fLogAtomicCrossSection[iZ]) << 296 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
>> 297 //not invoked
>> 298 if (!logAtomicCrossSection)
>> 299 {
>> 300 //create a **thread-local** version of the table. Used only for G4EmCalculator and
>> 301 //Unit Tests
>> 302 fLocalTable = true;
>> 303 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
>> 304 }
>> 305 //now it should be ok
>> 306 if (!logAtomicCrossSection->count(iZ))
278 { 307 {
279 //If we are here, it means that Initial << 308 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
280 //not filled up. This can happen in a U 309 //not filled up. This can happen in a UnitTest or via G4EmCalculator
281 if (fVerboseLevel > 0) << 310 if (verboseLevel > 0)
282 { 311 {
283 //Issue a G4Exception (warning) only in ve 312 //Issue a G4Exception (warning) only in verbose mode
284 G4ExceptionDescription ed; 313 G4ExceptionDescription ed;
285 ed << "Unable to retrieve the cross sectio 314 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
286 ed << "This can happen only in Unit Tests 315 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
287 G4Exception("G4PenelopeRayleighModel::Comp 316 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
288 "em2040",JustWarning,ed); 317 "em2040",JustWarning,ed);
289 } 318 }
290 //protect file reading via autolock 319 //protect file reading via autolock
291 G4AutoLock lock(&PenelopeRayleighModelM 320 G4AutoLock lock(&PenelopeRayleighModelMutex);
292 ReadDataFile(iZ); 321 ReadDataFile(iZ);
293 lock.unlock(); << 322 lock.unlock();
294 } 323 }
295 324
296 G4double cross = 0; 325 G4double cross = 0;
297 G4PhysicsFreeVector* atom = fLogAtomicCross << 326
>> 327 G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
298 if (!atom) 328 if (!atom)
299 { 329 {
300 G4ExceptionDescription ed; 330 G4ExceptionDescription ed;
301 ed << "Unable to find Z=" << iZ << " in 331 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
302 G4Exception("G4PenelopeRayleighModel::C 332 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
303 "em2041",FatalException,ed); 333 "em2041",FatalException,ed);
304 return 0; 334 return 0;
305 } 335 }
306 G4double logene = G4Log(energy); << 336 G4double logene = std::log(energy);
307 G4double logXS = atom->Value(logene); 337 G4double logXS = atom->Value(logene);
308 cross = G4Exp(logXS); << 338 cross = std::exp(logXS);
309 339
310 if (fVerboseLevel > 2) << 340 if (verboseLevel > 2)
311 { << 341 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z <<
312 G4cout << "Rayleigh cross section at " << 342 " = " << cross/barn << " barn" << G4endl;
313 " = " << cross/barn << " barn" << G4endl; << 343 return cross;
314 } <<
315 return cross; <<
316 } 344 }
317 345
318 346
319 //....oooOO0OOooo........oooOO0OOooo........oo 347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320 void G4PenelopeRayleighModel::BuildFormFactorT 348 void G4PenelopeRayleighModel::BuildFormFactorTable(const G4Material* material)
321 { 349 {
322 /* << 350
323 1) get composition and equivalent molecula << 351 /*
324 */ << 352 1) get composition and equivalent molecular density
325 std::size_t nElements = material->GetNumberO << 353 */
>> 354
>> 355 G4int nElements = material->GetNumberOfElements();
326 const G4ElementVector* elementVector = mater 356 const G4ElementVector* elementVector = material->GetElementVector();
327 const G4double* fractionVector = material->G 357 const G4double* fractionVector = material->GetFractionVector();
328 358
329 std::vector<G4double> *StechiometricFactors 359 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
330 for (std::size_t i=0;i<nElements;++i) << 360 for (G4int i=0;i<nElements;i++)
331 { 361 {
332 G4double fraction = fractionVector[i]; 362 G4double fraction = fractionVector[i];
333 G4double atomicWeigth = (*elementVector) 363 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
334 StechiometricFactors->push_back(fraction 364 StechiometricFactors->push_back(fraction/atomicWeigth);
335 } 365 }
336 //Find max 366 //Find max
337 G4double MaxStechiometricFactor = 0.; 367 G4double MaxStechiometricFactor = 0.;
338 for (std::size_t i=0;i<nElements;++i) << 368 for (G4int i=0;i<nElements;i++)
339 { 369 {
340 if ((*StechiometricFactors)[i] > MaxStec 370 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
341 MaxStechiometricFactor = (*Stechiometr 371 MaxStechiometricFactor = (*StechiometricFactors)[i];
342 } 372 }
343 if (MaxStechiometricFactor<1e-16) 373 if (MaxStechiometricFactor<1e-16)
344 { 374 {
345 G4ExceptionDescription ed; 375 G4ExceptionDescription ed;
346 ed << "Inconsistent data of atomic compo << 376 ed << "Inconsistent data of atomic composition for " <<
347 material->GetName() << G4endl; 377 material->GetName() << G4endl;
348 G4Exception("G4PenelopeRayleighModel::Bu 378 G4Exception("G4PenelopeRayleighModel::BuildFormFactorTable()",
349 "em2042",FatalException,ed); 379 "em2042",FatalException,ed);
350 } 380 }
351 //Normalize 381 //Normalize
352 for (std::size_t i=0;i<nElements;++i) << 382 for (G4int i=0;i<nElements;i++)
353 (*StechiometricFactors)[i] /= MaxStechiom 383 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
354 << 384
>> 385 // Equivalent atoms per molecule
>> 386 G4double atomsPerMolecule = 0;
>> 387 for (G4int i=0;i<nElements;i++)
>> 388 atomsPerMolecule += (*StechiometricFactors)[i];
>> 389
355 /* 390 /*
356 CREATE THE FORM FACTOR TABLE 391 CREATE THE FORM FACTOR TABLE
357 */ 392 */
358 G4PhysicsFreeVector* theFFVec = new G4Physic << 393 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size());
>> 394 theFFVec->SetSpline(true);
359 395
360 for (std::size_t k=0;k<fLogQSquareGrid.size( << 396 for (size_t k=0;k<logQSquareGrid.size();k++)
361 { 397 {
362 G4double ff2 = 0; //squared form factor 398 G4double ff2 = 0; //squared form factor
363 for (std::size_t i=0;i<nElements;++i) << 399 for (G4int i=0;i<nElements;i++)
364 { 400 {
365 G4int iZ = (*elementVector)[i]->GetZasInt( << 401 G4int iZ = (G4int) (*elementVector)[i]->GetZ();
366 G4PhysicsFreeVector* theAtomVec = fAtomicF << 402 G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second;
367 G4double f = (*theAtomVec)[k]; //the q-gri << 403 G4double f = (*theAtomVec)[k]; //the q-grid is always the same
368 ff2 += f*f*(*StechiometricFactors)[i]; 404 ff2 += f*f*(*StechiometricFactors)[i];
369 } 405 }
370 if (ff2) 406 if (ff2)
371 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo << 407 theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2)
372 } 408 }
373 theFFVec->FillSecondDerivatives(); //vector << 409 logFormFactorTable->insert(std::make_pair(material,theFFVec));
374 fLogFormFactorTable->insert(std::make_pair(m <<
375 410
376 delete StechiometricFactors; 411 delete StechiometricFactors;
>> 412
377 return; 413 return;
378 } 414 }
379 415
>> 416
380 //....oooOO0OOooo........oooOO0OOooo........oo 417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 418
382 void G4PenelopeRayleighModel::SampleSecondarie 419 void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
383 const G4MaterialCutsCouple* couple 420 const G4MaterialCutsCouple* couple,
384 const G4DynamicParticle* aDynamicG 421 const G4DynamicParticle* aDynamicGamma,
385 G4double, 422 G4double,
386 G4double) 423 G4double)
387 { 424 {
388 // Sampling of the Rayleigh final state (nam << 425 // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
389 // from the Penelope2008 model. The scatteri << 426 // from the Penelope2008 model. The scattering angle is sampled from the atomic
390 // cross section dOmega/d(cosTheta) from Bor << 427 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
391 // anomalous scattering effects. The Form Fa << 428 // anomalous scattering effects. The Form Factor F(Q) function which appears in the
392 // analytical cross section is retrieved via << 429 // analytical cross section is retrieved via the method GetFSquared(); atomic data
393 // are tabulated for F(Q). Form factor for c << 430 // are tabulated for F(Q). Form factor for compounds is calculated according to
394 // the additivity rule. The sampling from th << 431 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
395 // Transform with Aliasing (RITA) algorithm; << 432 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
396 // for each material and managed by G4Penelo 433 // for each material and managed by G4PenelopeSamplingData objects.
397 // The sampling algorithm (rejection method) << 434 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
398 // increases with energy. For E=100 keV the << 435 // increases with energy. For E=100 keV the efficiency is 100% and 86% for
399 // hydrogen and uranium, respectively. 436 // hydrogen and uranium, respectively.
400 437
401 if (fVerboseLevel > 3) << 438 if (verboseLevel > 3)
402 G4cout << "Calling SamplingSecondaries() o 439 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
403 440
404 G4double photonEnergy0 = aDynamicGamma->GetK 441 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
405 << 442
406 if (photonEnergy0 <= fIntrinsicLowEnergyLimi 443 if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
407 { 444 {
408 fParticleChange->ProposeTrackStatus(fSto 445 fParticleChange->ProposeTrackStatus(fStopAndKill);
409 fParticleChange->SetProposedKineticEnerg 446 fParticleChange->SetProposedKineticEnergy(0.);
410 fParticleChange->ProposeLocalEnergyDepos 447 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
411 return ; 448 return ;
412 } 449 }
413 450
414 G4ParticleMomentum photonDirection0 = aDynam 451 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
415 << 452
416 const G4Material* theMat = couple->GetMateri 453 const G4Material* theMat = couple->GetMaterial();
417 454
>> 455
418 //1) Verify if tables are ready 456 //1) Verify if tables are ready
419 //Either Initialize() was not called, or we << 457 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
420 //not invoked 458 //not invoked
421 if (!fPMaxTable || !fSamplingTable || !fLogF << 459 if (!pMaxTable || !samplingTable || !logAtomicCrossSection || !atomicFormFactor ||
>> 460 !logFormFactorTable)
422 { 461 {
423 //create a **thread-local** version of t << 462 //create a **thread-local** version of the table. Used only for G4EmCalculator and
424 //Unit Tests 463 //Unit Tests
425 fLocalTable = true; 464 fLocalTable = true;
426 if (!fLogFormFactorTable) << 465 if (!logAtomicCrossSection)
427 fLogFormFactorTable = new std::map<const G4M << 466 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
428 if (!fPMaxTable) << 467 if (!atomicFormFactor)
429 fPMaxTable = new std::map<const G4Material*, << 468 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
430 if (!fSamplingTable) << 469 if (!logFormFactorTable)
431 fSamplingTable = new std::map<const G4Materi << 470 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
432 } << 471 if (!pMaxTable)
433 << 472 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
434 if (!fSamplingTable->count(theMat)) << 473 if (!samplingTable)
435 { << 474 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
436 //If we are here, it means that Initiali << 475 }
437 //not filled up. This can happen in a Un << 476
438 if (fVerboseLevel > 0) << 477 if (!samplingTable->count(theMat))
>> 478 {
>> 479 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
>> 480 //not filled up. This can happen in a UnitTest
>> 481 if (verboseLevel > 0)
439 { 482 {
440 //Issue a G4Exception (warning) only in ve 483 //Issue a G4Exception (warning) only in verbose mode
441 G4ExceptionDescription ed; 484 G4ExceptionDescription ed;
442 ed << "Unable to find the fSamplingTable d << 485 ed << "Unable to find the samplingTable data for " <<
443 theMat->GetName() << G4endl; 486 theMat->GetName() << G4endl;
444 ed << "This can happen only in Unit Tests" << 487 ed << "This can happen only in Unit Tests" << G4endl;
445 G4Exception("G4PenelopeRayleighModel::Samp 488 G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
446 "em2019",JustWarning,ed); 489 "em2019",JustWarning,ed);
447 } 490 }
448 const G4ElementVector* theElementVector 491 const G4ElementVector* theElementVector = theMat->GetElementVector();
449 //protect file reading via autolock 492 //protect file reading via autolock
450 G4AutoLock lock(&PenelopeRayleighModelMu 493 G4AutoLock lock(&PenelopeRayleighModelMutex);
451 for (std::size_t j=0;j<theMat->GetNumber << 494 for (size_t j=0;j<theMat->GetNumberOfElements();j++)
452 { 495 {
453 G4int iZ = theElementVector->at(j)->GetZas << 496 G4int iZ = (G4int) theElementVector->at(j)->GetZ();
454 if (!fLogAtomicCrossSection[iZ]) << 497 if (!logAtomicCrossSection->count(iZ))
455 { 498 {
456 lock.lock(); 499 lock.lock();
457 ReadDataFile(iZ); 500 ReadDataFile(iZ);
458 lock.unlock(); << 501 lock.unlock();
459 } << 502 }
460 } 503 }
461 lock.lock(); 504 lock.lock();
462 //1) If the table has not been built for 505 //1) If the table has not been built for the material, do it!
463 if (!fLogFormFactorTable->count(theMat)) << 506 if (!logFormFactorTable->count(theMat))
464 BuildFormFactorTable(theMat); 507 BuildFormFactorTable(theMat);
465 508
466 //2) retrieve or build the sampling tabl 509 //2) retrieve or build the sampling table
467 if (!(fSamplingTable->count(theMat))) << 510 if (!(samplingTable->count(theMat)))
468 InitializeSamplingAlgorithm(theMat); 511 InitializeSamplingAlgorithm(theMat);
469 << 512
470 //3) retrieve or build the pMax data 513 //3) retrieve or build the pMax data
471 if (!fPMaxTable->count(theMat)) << 514 if (!pMaxTable->count(theMat))
472 GetPMaxTable(theMat); 515 GetPMaxTable(theMat);
473 lock.unlock(); 516 lock.unlock();
474 } 517 }
475 518
476 //Ok, restart the job 519 //Ok, restart the job
477 G4PenelopeSamplingData* theDataTable = fSamp << 520
478 G4PhysicsFreeVector* thePMax = fPMaxTable->f << 521 G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
>> 522 G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;
479 523
480 G4double cosTheta = 1.0; 524 G4double cosTheta = 1.0;
481 << 525
482 //OK, ready to go! 526 //OK, ready to go!
483 G4double qmax = 2.0*photonEnergy0/electron_m 527 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
484 528
485 if (qmax < 1e-10) //very low momentum transf 529 if (qmax < 1e-10) //very low momentum transfer
486 { 530 {
487 G4bool loopAgain=false; 531 G4bool loopAgain=false;
488 do 532 do
489 { 533 {
490 loopAgain = false; 534 loopAgain = false;
491 cosTheta = 1.0-2.0*G4UniformRand(); 535 cosTheta = 1.0-2.0*G4UniformRand();
492 G4double G = 0.5*(1+cosTheta*cosTheta); 536 G4double G = 0.5*(1+cosTheta*cosTheta);
493 if (G4UniformRand()>G) 537 if (G4UniformRand()>G)
494 loopAgain = true; 538 loopAgain = true;
495 }while(loopAgain); 539 }while(loopAgain);
496 } 540 }
497 else //larger momentum transfer 541 else //larger momentum transfer
498 { 542 {
499 std::size_t nData = theDataTable->GetNum << 543 size_t nData = theDataTable->GetNumberOfStoredPoints();
500 G4double LastQ2inTheTable = theDataTable 544 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
501 G4double q2max = std::min(qmax*qmax,Last 545 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
502 546
503 G4bool loopAgain = false; 547 G4bool loopAgain = false;
504 G4double MaxPValue = thePMax->Value(phot 548 G4double MaxPValue = thePMax->Value(photonEnergy0);
505 G4double xx=0; 549 G4double xx=0;
506 << 550
507 //Sampling by rejection method. The reje << 551 //Sampling by rejection method. The rejection function is
508 //G = 0.5*(1+cos^2(theta)) 552 //G = 0.5*(1+cos^2(theta))
509 // 553 //
510 do{ 554 do{
511 loopAgain = false; 555 loopAgain = false;
512 G4double RandomMax = G4UniformRand()*MaxPVal 556 G4double RandomMax = G4UniformRand()*MaxPValue;
513 xx = theDataTable->SampleValue(RandomMax); 557 xx = theDataTable->SampleValue(RandomMax);
514 //xx is a random value of q^2 in (0,q2max),s << 558 //xx is a random value of q^2 in (0,q2max),sampled according to
515 //F(Q^2) via the RITA algorithm 559 //F(Q^2) via the RITA algorithm
516 if (xx > q2max) 560 if (xx > q2max)
517 loopAgain = true; 561 loopAgain = true;
518 cosTheta = 1.0-2.0*xx/q2max; 562 cosTheta = 1.0-2.0*xx/q2max;
519 G4double G = 0.5*(1+cosTheta*cosTheta); 563 G4double G = 0.5*(1+cosTheta*cosTheta);
520 if (G4UniformRand()>G) 564 if (G4UniformRand()>G)
521 loopAgain = true; 565 loopAgain = true;
522 }while(loopAgain); 566 }while(loopAgain);
523 } 567 }
524 << 568
525 G4double sinTheta = std::sqrt(1-cosTheta*cos 569 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
526 << 570
527 // Scattered photon angles. ( Z - axis along 571 // Scattered photon angles. ( Z - axis along the parent photon)
528 G4double phi = twopi * G4UniformRand() ; 572 G4double phi = twopi * G4UniformRand() ;
529 G4double dirX = sinTheta*std::cos(phi); 573 G4double dirX = sinTheta*std::cos(phi);
530 G4double dirY = sinTheta*std::sin(phi); 574 G4double dirY = sinTheta*std::sin(phi);
531 G4double dirZ = cosTheta; 575 G4double dirZ = cosTheta;
532 << 576
533 // Update G4VParticleChange for the scattere << 577 // Update G4VParticleChange for the scattered photon
534 G4ThreeVector photonDirection1(dirX, dirY, d 578 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
535 579
536 photonDirection1.rotateUz(photonDirection0); 580 photonDirection1.rotateUz(photonDirection0);
537 fParticleChange->ProposeMomentumDirection(ph 581 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
538 fParticleChange->SetProposedKineticEnergy(ph 582 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
539 << 583
540 return; 584 return;
541 } 585 }
542 586
543 587
544 //....oooOO0OOooo........oooOO0OOooo........oo 588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
545 589
546 void G4PenelopeRayleighModel::ReadDataFile(con 590 void G4PenelopeRayleighModel::ReadDataFile(const G4int Z)
547 { 591 {
548 if (fVerboseLevel > 2) << 592
>> 593 if (verboseLevel > 2)
549 { 594 {
550 G4cout << "G4PenelopeRayleighModel::Read 595 G4cout << "G4PenelopeRayleighModel::ReadDataFile()" << G4endl;
551 G4cout << "Going to read Rayleigh data f 596 G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
552 } 597 }
553 const char* path = G4FindDataDir("G4LEDATA << 598
554 if(!path) << 599 char* path = getenv("G4LEDATA");
>> 600 if (!path)
555 { 601 {
556 G4String excep = "G4LEDATA environment v 602 G4String excep = "G4LEDATA environment variable not set!";
557 G4Exception("G4PenelopeRayleighModel::Re 603 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
558 "em0006",FatalException,excep); 604 "em0006",FatalException,excep);
559 return; 605 return;
560 } 606 }
561 607
562 /* 608 /*
563 Read first the cross section file 609 Read first the cross section file
564 */ 610 */
565 std::ostringstream ost; 611 std::ostringstream ost;
566 if (Z>9) 612 if (Z>9)
567 ost << path << "/penelope/rayleigh/pdgra" 613 ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
568 else 614 else
569 ost << path << "/penelope/rayleigh/pdgra0" 615 ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
570 std::ifstream file(ost.str().c_str()); 616 std::ifstream file(ost.str().c_str());
571 if (!file.is_open()) 617 if (!file.is_open())
572 { 618 {
573 G4String excep = "Data file " + G4String 619 G4String excep = "Data file " + G4String(ost.str()) + " not found!";
574 G4Exception("G4PenelopeRayleighModel::Re 620 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
575 "em0003",FatalException,excep); 621 "em0003",FatalException,excep);
576 } 622 }
577 G4int readZ =0; 623 G4int readZ =0;
578 std::size_t nPoints= 0; << 624 size_t nPoints= 0;
579 file >> readZ >> nPoints; 625 file >> readZ >> nPoints;
580 //check the right file is opened. 626 //check the right file is opened.
581 if (readZ != Z || nPoints <= 0 || nPoints >= 627 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
582 { 628 {
583 G4ExceptionDescription ed; 629 G4ExceptionDescription ed;
584 ed << "Corrupted data file for Z=" << Z 630 ed << "Corrupted data file for Z=" << Z << G4endl;
585 G4Exception("G4PenelopeRayleighModel::Re 631 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
586 "em0005",FatalException,ed); 632 "em0005",FatalException,ed);
587 return; 633 return;
588 } << 634 }
589 << 635 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints);
590 fLogAtomicCrossSection[Z] = new G4PhysicsFre <<
591 G4double ene=0,f1=0,f2=0,xs=0; 636 G4double ene=0,f1=0,f2=0,xs=0;
592 for (std::size_t i=0;i<nPoints;++i) << 637 for (size_t i=0;i<nPoints;i++)
593 { 638 {
594 file >> ene >> f1 >> f2 >> xs; 639 file >> ene >> f1 >> f2 >> xs;
595 //dimensional quantities 640 //dimensional quantities
596 ene *= eV; 641 ene *= eV;
597 xs *= cm2; 642 xs *= cm2;
598 fLogAtomicCrossSection[Z]->PutValue(i,G4 << 643 theVec->PutValue(i,std::log(ene),std::log(xs));
599 if (file.eof() && i != (nPoints-1)) //fi 644 if (file.eof() && i != (nPoints-1)) //file ended too early
600 { 645 {
601 G4ExceptionDescription ed ; << 646 G4ExceptionDescription ed ;
602 ed << "Corrupted data file for Z=" << Z << 647 ed << "Corrupted data file for Z=" << Z << G4endl;
603 ed << "Found less than " << nPoints << "en 648 ed << "Found less than " << nPoints << "entries " <<G4endl;
604 G4Exception("G4PenelopeRayleighModel::Read 649 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
605 "em0005",FatalException,ed); 650 "em0005",FatalException,ed);
606 } 651 }
607 } 652 }
>> 653 if (!logAtomicCrossSection)
>> 654 {
>> 655 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
>> 656 "em2044",FatalException,"Unable to allocate the atomic cross section table");
>> 657 delete theVec;
>> 658 return;
>> 659 }
>> 660 logAtomicCrossSection->insert(std::make_pair(Z,theVec));
608 file.close(); 661 file.close();
609 662
610 /* 663 /*
611 Then read the form factor file 664 Then read the form factor file
612 */ 665 */
613 std::ostringstream ost2; 666 std::ostringstream ost2;
614 if (Z>9) 667 if (Z>9)
615 ost2 << path << "/penelope/rayleigh/pdaff" 668 ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
616 else 669 else
617 ost2 << path << "/penelope/rayleigh/pdaff0 670 ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
618 file.open(ost2.str().c_str()); 671 file.open(ost2.str().c_str());
619 if (!file.is_open()) 672 if (!file.is_open())
620 { 673 {
621 G4String excep = "Data file " + G4String 674 G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
622 G4Exception("G4PenelopeRayleighModel::Re 675 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
623 "em0003",FatalException,excep); 676 "em0003",FatalException,excep);
624 } 677 }
625 file >> readZ >> nPoints; 678 file >> readZ >> nPoints;
626 //check the right file is opened. 679 //check the right file is opened.
627 if (readZ != Z || nPoints <= 0 || nPoints >= 680 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
628 { 681 {
629 G4ExceptionDescription ed; 682 G4ExceptionDescription ed;
630 ed << "Corrupted data file for Z=" << Z 683 ed << "Corrupted data file for Z=" << Z << G4endl;
631 G4Exception("G4PenelopeRayleighModel::Re 684 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
632 "em0005",FatalException,ed); 685 "em0005",FatalException,ed);
633 return; 686 return;
634 } << 687 }
635 fAtomicFormFactor[Z] = new G4PhysicsFreeVect << 688 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
636 G4double q=0,ff=0,incoh=0; 689 G4double q=0,ff=0,incoh=0;
637 G4bool fillQGrid = false; 690 G4bool fillQGrid = false;
638 //fill this vector only the first time. 691 //fill this vector only the first time.
639 if (!fLogQSquareGrid.size()) << 692 if (!logQSquareGrid.size())
640 fillQGrid = true; 693 fillQGrid = true;
641 for (std::size_t i=0;i<nPoints;++i) << 694 for (size_t i=0;i<nPoints;i++)
642 { 695 {
643 file >> q >> ff >> incoh; 696 file >> q >> ff >> incoh;
644 //q and ff are dimensionless (q is in un 697 //q and ff are dimensionless (q is in units of (m_e*c)
645 fAtomicFormFactor[Z]->PutValue(i,q,ff); << 698 theFFVec->PutValue(i,q,ff);
646 if (fillQGrid) 699 if (fillQGrid)
647 { 700 {
648 fLogQSquareGrid.push_back(2.0*G4Log(q)); << 701 logQSquareGrid.push_back(2.0*std::log(q));
649 } 702 }
650 if (file.eof() && i != (nPoints-1)) //fi 703 if (file.eof() && i != (nPoints-1)) //file ended too early
651 { 704 {
652 G4ExceptionDescription ed; << 705 G4ExceptionDescription ed;
653 ed << "Corrupted data file for Z=" << Z << 706 ed << "Corrupted data file for Z=" << Z << G4endl;
654 ed << "Found less than " << nPoints << "en 707 ed << "Found less than " << nPoints << "entries " <<G4endl;
655 G4Exception("G4PenelopeRayleighModel::Read 708 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
656 "em0005",FatalException,ed); 709 "em0005",FatalException,ed);
657 } 710 }
658 } 711 }
>> 712 if (!atomicFormFactor)
>> 713 {
>> 714 G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
>> 715 "em2045",FatalException,
>> 716 "Unable to allocate the atomicFormFactor data table");
>> 717 delete theFFVec;
>> 718 return;
>> 719 }
>> 720 atomicFormFactor->insert(std::make_pair(Z,theFFVec));
659 file.close(); 721 file.close();
660 return; 722 return;
661 } 723 }
662 724
663 //....oooOO0OOooo........oooOO0OOooo........oo 725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
664 726
665 G4double G4PenelopeRayleighModel::GetFSquared( 727 G4double G4PenelopeRayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared)
666 { 728 {
667 G4double f2 = 0; 729 G4double f2 = 0;
668 //Input value QSquared could be zero: protec << 730 //Input value QSquared could be zero: protect the log() below against
669 //the FPE exception 731 //the FPE exception
670 //If Q<1e-10, set Q to 1e-10 732 //If Q<1e-10, set Q to 1e-10
671 G4double logQSquared = (QSquared>1e-10) ? G4 << 733 G4double logQSquared = (QSquared>1e-10) ? std::log(QSquared) : -23.;
672 //last value of the table 734 //last value of the table
673 G4double maxlogQ2 = fLogQSquareGrid[fLogQSqu << 735 G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1];
674 << 736
>> 737
675 //now it should be all right 738 //now it should be all right
676 G4PhysicsFreeVector* theVec = fLogFormFactor << 739 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
677 740
678 if (!theVec) 741 if (!theVec)
679 { 742 {
680 G4ExceptionDescription ed; 743 G4ExceptionDescription ed;
681 ed << "Unable to retrieve F squared tabl 744 ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
682 G4Exception("G4PenelopeRayleighModel::Ge 745 G4Exception("G4PenelopeRayleighModel::GetFSquared()",
683 "em2046",FatalException,ed); 746 "em2046",FatalException,ed);
684 return 0; 747 return 0;
685 } 748 }
686 if (logQSquared < -20) // Q < 1e-9 749 if (logQSquared < -20) // Q < 1e-9
687 { 750 {
688 G4double logf2 = (*theVec)[0]; //first v 751 G4double logf2 = (*theVec)[0]; //first value of the table
689 f2 = G4Exp(logf2); << 752 f2 = std::exp(logf2);
690 } 753 }
691 else if (logQSquared > maxlogQ2) 754 else if (logQSquared > maxlogQ2)
692 f2 =0; 755 f2 =0;
693 else 756 else
694 { 757 {
695 //log(Q^2) vs. log(F^2) 758 //log(Q^2) vs. log(F^2)
696 G4double logf2 = theVec->Value(logQSquar 759 G4double logf2 = theVec->Value(logQSquared);
697 f2 = G4Exp(logf2); << 760 f2 = std::exp(logf2);
698 761
699 } 762 }
700 if (fVerboseLevel > 3) << 763 if (verboseLevel > 3)
701 { 764 {
702 G4cout << "G4PenelopeRayleighModel::GetF << 765 G4cout << "G4PenelopeRayleighModel::GetFSquared() in " << mat->GetName() << G4endl;
703 G4cout << "Q^2 = " << QSquared << " (un 766 G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl;
704 } 767 }
705 return f2; 768 return f2;
706 } 769 }
707 770
708 //....oooOO0OOooo........oooOO0OOooo........oo 771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
709 772
710 void G4PenelopeRayleighModel::InitializeSampli 773 void G4PenelopeRayleighModel::InitializeSamplingAlgorithm(const G4Material* mat)
711 { 774 {
>> 775
712 G4double q2min = 0; 776 G4double q2min = 0;
713 G4double q2max = 0; 777 G4double q2max = 0;
714 const std::size_t np = 150; //hard-coded in << 778 const size_t np = 150; //hard-coded in Penelope
715 //G4cout << "Init N= " << fLogQSquareGrid.si << 779 //G4cout << "Init N= " << logQSquareGrid.size() << G4endl;
716 for (std::size_t i=1;i<fLogQSquareGrid.size( << 780 for (size_t i=1;i<logQSquareGrid.size();i++)
717 { 781 {
718 G4double Q2 = G4Exp(fLogQSquareGrid[i]); << 782 G4double Q2 = std::exp(logQSquareGrid[i]);
719 if (GetFSquared(mat,Q2) > 1e-35) 783 if (GetFSquared(mat,Q2) > 1e-35)
720 { 784 {
721 q2max = G4Exp(fLogQSquareGrid[i-1]); << 785 q2max = std::exp(logQSquareGrid[i-1]);
722 } 786 }
723 //G4cout << "Q2= " << Q2 << " q2max= " < 787 //G4cout << "Q2= " << Q2 << " q2max= " << q2max << G4endl;
724 } 788 }
725 << 789
726 std::size_t nReducedPoints = np/4; << 790 size_t nReducedPoints = np/4;
727 791
728 //check for errors 792 //check for errors
729 if (np < 16) 793 if (np < 16)
730 { 794 {
731 G4Exception("G4PenelopeRayleighModel::In 795 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
732 "em2047",FatalException, 796 "em2047",FatalException,
733 "Too few points to initialize the sampli 797 "Too few points to initialize the sampling algorithm");
734 } 798 }
735 if (q2min > (q2max-1e-10)) 799 if (q2min > (q2max-1e-10))
736 { 800 {
737 G4cout << "q2min= " << q2min << " q2max= 801 G4cout << "q2min= " << q2min << " q2max= " << q2max << G4endl;
738 G4Exception("G4PenelopeRayleighModel::In 802 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
739 "em2048",FatalException, 803 "em2048",FatalException,
740 "Too narrow grid to initialize the sampl 804 "Too narrow grid to initialize the sampling algorithm");
741 } 805 }
742 806
743 //This is subroutine RITAI0 of Penelope 807 //This is subroutine RITAI0 of Penelope
744 //Create an object of type G4PenelopeRayleig 808 //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material*
745 809
746 //temporary vectors --> Then everything is s 810 //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
747 G4DataVector* x = new G4DataVector(); 811 G4DataVector* x = new G4DataVector();
748 << 812
749 /******************************************* 813 /*******************************************************************************
750 Start with a grid of NUNIF points uniforml 814 Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
751 ******************************************** 815 ********************************************************************************/
752 std::size_t NUNIF = std::min(std::max(((std: << 816 size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2);
753 const G4int nip = 51; //hard-coded in Penelo 817 const G4int nip = 51; //hard-coded in Penelope
754 818
755 G4double dx = (q2max-q2min)/((G4double) NUNI << 819 G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
756 x->push_back(q2min); << 820 x->push_back(q2min);
757 for (std::size_t i=1;i<NUNIF-1;++i) << 821 for (size_t i=1;i<NUNIF-1;i++)
758 { 822 {
759 G4double app = q2min + i*dx; 823 G4double app = q2min + i*dx;
760 x->push_back(app); //increase << 824 x->push_back(app); //increase
761 } 825 }
762 x->push_back(q2max); 826 x->push_back(q2max);
763 << 827
764 if (fVerboseLevel> 3) << 828 if (verboseLevel> 3)
765 G4cout << "Vector x has " << x->size() << 829 G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
766 830
767 //Allocate temporary storage vectors 831 //Allocate temporary storage vectors
768 G4DataVector* area = new G4DataVector(); 832 G4DataVector* area = new G4DataVector();
769 G4DataVector* a = new G4DataVector(); 833 G4DataVector* a = new G4DataVector();
770 G4DataVector* b = new G4DataVector(); 834 G4DataVector* b = new G4DataVector();
771 G4DataVector* c = new G4DataVector(); 835 G4DataVector* c = new G4DataVector();
772 G4DataVector* err = new G4DataVector(); 836 G4DataVector* err = new G4DataVector();
773 837
774 for (std::size_t i=0;i<NUNIF-1;++i) //build << 838 for (size_t i=0;i<NUNIF-1;i++) //build all points but the last
775 { << 839 {
776 //Temporary vectors for this loop 840 //Temporary vectors for this loop
777 G4DataVector* pdfi = new G4DataVector(); 841 G4DataVector* pdfi = new G4DataVector();
778 G4DataVector* pdfih = new G4DataVector() 842 G4DataVector* pdfih = new G4DataVector();
779 G4DataVector* sumi = new G4DataVector(); 843 G4DataVector* sumi = new G4DataVector();
780 844
781 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4do 845 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
782 G4double pdfmax = 0; 846 G4double pdfmax = 0;
783 for (G4int k=0;k<nip;k++) 847 for (G4int k=0;k<nip;k++)
784 { 848 {
785 G4double xik = (*x)[i]+k*dxi; << 849 G4double xik = (*x)[i]+k*dxi;
786 G4double pdfk = std::max(GetFSquared(mat,x 850 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
787 pdfi->push_back(pdfk); 851 pdfi->push_back(pdfk);
788 pdfmax = std::max(pdfmax,pdfk); << 852 pdfmax = std::max(pdfmax,pdfk);
789 if (k < (nip-1)) 853 if (k < (nip-1))
790 { 854 {
791 G4double xih = xik + 0.5*dxi; 855 G4double xih = xik + 0.5*dxi;
792 G4double pdfIK = std::max(GetFSquared( 856 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
793 pdfih->push_back(pdfIK); 857 pdfih->push_back(pdfIK);
794 pdfmax = std::max(pdfmax,pdfIK); 858 pdfmax = std::max(pdfmax,pdfIK);
795 } 859 }
796 } 860 }
797 << 861
798 //Simpson's integration 862 //Simpson's integration
799 G4double cons = dxi*0.5*(1./3.); 863 G4double cons = dxi*0.5*(1./3.);
800 sumi->push_back(0.); 864 sumi->push_back(0.);
801 for (G4int k=1;k<nip;k++) 865 for (G4int k=1;k<nip;k++)
802 { 866 {
803 G4double previous = (*sumi)[k-1]; 867 G4double previous = (*sumi)[k-1];
804 G4double next = previous + cons*((*pdfi)[k 868 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
805 sumi->push_back(next); 869 sumi->push_back(next);
806 } 870 }
807 << 871
808 G4double lastIntegral = (*sumi)[sumi->si 872 G4double lastIntegral = (*sumi)[sumi->size()-1];
809 area->push_back(lastIntegral); 873 area->push_back(lastIntegral);
810 //Normalize cumulative function 874 //Normalize cumulative function
811 G4double factor = 1.0/lastIntegral; 875 G4double factor = 1.0/lastIntegral;
812 for (std::size_t k=0;k<sumi->size();++k) << 876 for (size_t k=0;k<sumi->size();k++)
813 (*sumi)[k] *= factor; 877 (*sumi)[k] *= factor;
814 << 878
815 //When the PDF vanishes at one of the in 879 //When the PDF vanishes at one of the interval end points, its value is modified
816 if ((*pdfi)[0] < 1e-35) << 880 if ((*pdfi)[0] < 1e-35)
817 (*pdfi)[0] = 1e-5*pdfmax; 881 (*pdfi)[0] = 1e-5*pdfmax;
818 if ((*pdfi)[pdfi->size()-1] < 1e-35) 882 if ((*pdfi)[pdfi->size()-1] < 1e-35)
819 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 883 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
820 884
821 G4double pli = (*pdfi)[0]*factor; 885 G4double pli = (*pdfi)[0]*factor;
822 G4double pui = (*pdfi)[pdfi->size()-1]*f 886 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
823 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx 887 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
824 G4double A_temp = (1.0/(pli*dx))-1.0-B_t 888 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
825 G4double C_temp = 1.0+A_temp+B_temp; 889 G4double C_temp = 1.0+A_temp+B_temp;
826 if (C_temp < 1e-35) 890 if (C_temp < 1e-35)
827 { 891 {
828 a->push_back(0.); 892 a->push_back(0.);
829 b->push_back(0.); 893 b->push_back(0.);
830 c->push_back(1.); << 894 c->push_back(1.);
831 } 895 }
832 else 896 else
833 { 897 {
834 a->push_back(A_temp); 898 a->push_back(A_temp);
835 b->push_back(B_temp); 899 b->push_back(B_temp);
836 c->push_back(C_temp); 900 c->push_back(C_temp);
837 } 901 }
838 902
839 //OK, now get ERR(I), the integral of th << 903 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
840 //and the true pdf, extended over the in 904 //and the true pdf, extended over the interval (X(I),X(I+1))
841 G4int icase = 1; //loop code 905 G4int icase = 1; //loop code
842 G4bool reLoop = false; 906 G4bool reLoop = false;
843 err->push_back(0.); 907 err->push_back(0.);
844 do 908 do
845 { 909 {
846 reLoop = false; 910 reLoop = false;
847 (*err)[i] = 0.; //zero variable 911 (*err)[i] = 0.; //zero variable
848 for (G4int k=0;k<nip;k++) 912 for (G4int k=0;k<nip;k++)
849 { 913 {
850 G4double rr = (*sumi)[k]; 914 G4double rr = (*sumi)[k];
851 G4double pap = (*area)[i]*(1.0+((*a)[i 915 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
852 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(* 916 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
853 if (k == 0 || k == nip-1) 917 if (k == 0 || k == nip-1)
854 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]) 918 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
855 else 919 else
856 (*err)[i] += std::fabs(pap-(*pdfi)[k]); 920 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
857 } 921 }
858 (*err)[i] *= dxi; 922 (*err)[i] *= dxi;
859 << 923
860 //If err(I) is too large, the pdf is appro 924 //If err(I) is too large, the pdf is approximated by a uniform distribution
861 if ((*err)[i] > 0.1*(*area)[i] && icase == << 925 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
862 { 926 {
863 (*b)[i] = 0; 927 (*b)[i] = 0;
864 (*a)[i] = 0; 928 (*a)[i] = 0;
865 (*c)[i] = 1.; 929 (*c)[i] = 1.;
866 icase = 2; 930 icase = 2;
867 reLoop = true; 931 reLoop = true;
868 } 932 }
869 }while(reLoop); 933 }while(reLoop);
>> 934
870 delete pdfi; 935 delete pdfi;
871 delete pdfih; 936 delete pdfih;
872 delete sumi; 937 delete sumi;
873 } //end of first loop over i 938 } //end of first loop over i
874 939
875 //Now assign last point 940 //Now assign last point
876 (*x)[x->size()-1] = q2max; 941 (*x)[x->size()-1] = q2max;
877 a->push_back(0.); 942 a->push_back(0.);
878 b->push_back(0.); 943 b->push_back(0.);
879 c->push_back(0.); 944 c->push_back(0.);
880 err->push_back(0.); 945 err->push_back(0.);
881 area->push_back(0.); 946 area->push_back(0.);
882 947
883 if (x->size() != NUNIF || a->size() != NUNIF << 948 if (x->size() != NUNIF || a->size() != NUNIF ||
884 err->size() != NUNIF || area->size() != 949 err->size() != NUNIF || area->size() != NUNIF)
885 { 950 {
886 G4ExceptionDescription ed; 951 G4ExceptionDescription ed;
887 ed << "Problem in building the Table for 952 ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
888 G4Exception("G4PenelopeRayleighModel::In 953 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
889 "em2049",FatalException,ed); 954 "em2049",FatalException,ed);
890 } 955 }
891 << 956
892 /******************************************* 957 /*******************************************************************************
893 New grid points are added by halving the su 958 New grid points are added by halving the sub-intervals with the largest absolute error
894 This is done up to np=150 points in the grid 959 This is done up to np=150 points in the grid
895 ******************************************** 960 ********************************************************************************/
896 do 961 do
897 { 962 {
898 G4double maxError = 0.0; 963 G4double maxError = 0.0;
899 std::size_t iErrMax = 0; << 964 size_t iErrMax = 0;
900 for (std::size_t i=0;i<err->size()-2;++i << 965 for (size_t i=0;i<err->size()-2;i++)
901 { 966 {
902 //maxError is the lagest of the interval e 967 //maxError is the lagest of the interval errors err[i]
903 if ((*err)[i] > maxError) 968 if ((*err)[i] > maxError)
904 { 969 {
905 maxError = (*err)[i]; 970 maxError = (*err)[i];
906 iErrMax = i; 971 iErrMax = i;
907 } 972 }
908 } 973 }
909 << 974
910 //OK, now I have to insert one new point 975 //OK, now I have to insert one new point in the position iErrMax
911 G4double newx = 0.5*((*x)[iErrMax]+(*x)[ 976 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
912 << 977
913 x->insert(x->begin()+iErrMax+1,newx); 978 x->insert(x->begin()+iErrMax+1,newx);
914 //Add place-holders in the other vectors 979 //Add place-holders in the other vectors
915 area->insert(area->begin()+iErrMax+1,0.) 980 area->insert(area->begin()+iErrMax+1,0.);
916 a->insert(a->begin()+iErrMax+1,0.); 981 a->insert(a->begin()+iErrMax+1,0.);
917 b->insert(b->begin()+iErrMax+1,0.); 982 b->insert(b->begin()+iErrMax+1,0.);
918 c->insert(c->begin()+iErrMax+1,0.); 983 c->insert(c->begin()+iErrMax+1,0.);
919 err->insert(err->begin()+iErrMax+1,0.); 984 err->insert(err->begin()+iErrMax+1,0.);
920 << 985
921 //Now calculate the other parameters 986 //Now calculate the other parameters
922 for (std::size_t i=iErrMax;i<=iErrMax+1; << 987 for (size_t i=iErrMax;i<=iErrMax+1;i++)
923 { 988 {
924 //Temporary vectors for this loop 989 //Temporary vectors for this loop
925 G4DataVector* pdfi = new G4DataVector(); 990 G4DataVector* pdfi = new G4DataVector();
926 G4DataVector* pdfih = new G4DataVector(); 991 G4DataVector* pdfih = new G4DataVector();
927 G4DataVector* sumi = new G4DataVector(); 992 G4DataVector* sumi = new G4DataVector();
928 << 993
929 G4double dxLocal = (*x)[i+1]-(*x)[i]; 994 G4double dxLocal = (*x)[i+1]-(*x)[i];
930 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4doub 995 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
931 G4double pdfmax = 0; 996 G4double pdfmax = 0;
932 for (G4int k=0;k<nip;k++) 997 for (G4int k=0;k<nip;k++)
933 { 998 {
934 G4double xik = (*x)[i]+k*dxi; 999 G4double xik = (*x)[i]+k*dxi;
935 G4double pdfk = std::max(GetFSquared(m 1000 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
936 pdfi->push_back(pdfk); 1001 pdfi->push_back(pdfk);
937 pdfmax = std::max(pdfmax,pdfk); << 1002 pdfmax = std::max(pdfmax,pdfk);
938 if (k < (nip-1)) 1003 if (k < (nip-1))
939 { 1004 {
940 G4double xih = xik + 0.5*dxi; 1005 G4double xih = xik + 0.5*dxi;
941 G4double pdfIK = std::max(GetFSquared(ma 1006 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
942 pdfih->push_back(pdfIK); 1007 pdfih->push_back(pdfIK);
943 pdfmax = std::max(pdfmax,pdfIK); 1008 pdfmax = std::max(pdfmax,pdfIK);
944 } 1009 }
945 } 1010 }
946 << 1011
947 //Simpson's integration 1012 //Simpson's integration
948 G4double cons = dxi*0.5*(1./3.); 1013 G4double cons = dxi*0.5*(1./3.);
949 sumi->push_back(0.); 1014 sumi->push_back(0.);
950 for (G4int k=1;k<nip;k++) 1015 for (G4int k=1;k<nip;k++)
951 { 1016 {
952 G4double previous = (*sumi)[k-1]; 1017 G4double previous = (*sumi)[k-1];
953 G4double next = previous + cons*((*pdf 1018 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
954 sumi->push_back(next); 1019 sumi->push_back(next);
955 } 1020 }
956 G4double lastIntegral = (*sumi)[sumi->size 1021 G4double lastIntegral = (*sumi)[sumi->size()-1];
957 (*area)[i] = lastIntegral; 1022 (*area)[i] = lastIntegral;
958 << 1023
959 //Normalize cumulative function 1024 //Normalize cumulative function
960 G4double factor = 1.0/lastIntegral; 1025 G4double factor = 1.0/lastIntegral;
961 for (std::size_t k=0;k<sumi->size();++k) << 1026 for (size_t k=0;k<sumi->size();k++)
962 (*sumi)[k] *= factor; 1027 (*sumi)[k] *= factor;
963 << 1028
964 //When the PDF vanishes at one of the inte 1029 //When the PDF vanishes at one of the interval end points, its value is modified
965 if ((*pdfi)[0] < 1e-35) << 1030 if ((*pdfi)[0] < 1e-35)
966 (*pdfi)[0] = 1e-5*pdfmax; 1031 (*pdfi)[0] = 1e-5*pdfmax;
967 if ((*pdfi)[pdfi->size()-1] < 1e-35) 1032 if ((*pdfi)[pdfi->size()-1] < 1e-35)
968 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 1033 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
969 << 1034
970 G4double pli = (*pdfi)[0]*factor; 1035 G4double pli = (*pdfi)[0]*factor;
971 G4double pui = (*pdfi)[pdfi->size()-1]*fac 1036 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
972 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal 1037 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
973 G4double A_temp = (1.0/(pli*dxLocal))-1.0- 1038 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
974 G4double C_temp = 1.0+A_temp+B_temp; 1039 G4double C_temp = 1.0+A_temp+B_temp;
975 if (C_temp < 1e-35) 1040 if (C_temp < 1e-35)
976 { 1041 {
977 (*a)[i]= 0.; 1042 (*a)[i]= 0.;
978 (*b)[i] = 0.; 1043 (*b)[i] = 0.;
979 (*c)[i] = 1; 1044 (*c)[i] = 1;
980 } 1045 }
981 else 1046 else
982 { 1047 {
983 (*a)[i]= A_temp; 1048 (*a)[i]= A_temp;
984 (*b)[i] = B_temp; 1049 (*b)[i] = B_temp;
985 (*c)[i] = C_temp; 1050 (*c)[i] = C_temp;
986 } 1051 }
987 //OK, now get ERR(I), the integral of the << 1052 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
988 //and the true pdf, extended over the inte 1053 //and the true pdf, extended over the interval (X(I),X(I+1))
989 G4int icase = 1; //loop code 1054 G4int icase = 1; //loop code
990 G4bool reLoop = false; 1055 G4bool reLoop = false;
991 do 1056 do
992 { 1057 {
993 reLoop = false; 1058 reLoop = false;
994 (*err)[i] = 0.; //zero variable 1059 (*err)[i] = 0.; //zero variable
995 for (G4int k=0;k<nip;k++) 1060 for (G4int k=0;k<nip;k++)
996 { 1061 {
997 G4double rr = (*sumi)[k]; << 1062 G4double rr = (*sumi)[k];
998 G4double pap = (*area)[i]*(1.0+((*a)[i]+ 1063 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
999 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1 1064 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1000 if (k == 0 || k == nip-1) 1065 if (k == 0 || k == nip-1)
1001 (*err)[i] += 0.5*std::fabs(pap-(*pdfi 1066 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1002 else 1067 else
1003 (*err)[i] += std::fabs(pap-(*pdfi)[k] 1068 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1004 } 1069 }
1005 (*err)[i] *= dxi; 1070 (*err)[i] *= dxi;
1006 << 1071
1007 //If err(I) is too large, the pdf is 1072 //If err(I) is too large, the pdf is approximated by a uniform distribution
1008 if ((*err)[i] > 0.1*(*area)[i] && ica << 1073 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1009 { 1074 {
1010 (*b)[i] = 0; 1075 (*b)[i] = 0;
1011 (*a)[i] = 0; 1076 (*a)[i] = 0;
1012 (*c)[i] = 1.; 1077 (*c)[i] = 1.;
1013 icase = 2; 1078 icase = 2;
1014 reLoop = true; 1079 reLoop = true;
1015 } 1080 }
1016 }while(reLoop); 1081 }while(reLoop);
1017 delete pdfi; 1082 delete pdfi;
1018 delete pdfih; 1083 delete pdfih;
1019 delete sumi; 1084 delete sumi;
1020 } 1085 }
1021 }while(x->size() < np); 1086 }while(x->size() < np);
1022 1087
1023 if (x->size() != np || a->size() != np || << 1088 if (x->size() != np || a->size() != np ||
1024 err->size() != np || area->size() != np 1089 err->size() != np || area->size() != np)
1025 { 1090 {
1026 G4Exception("G4PenelopeRayleighModel::I 1091 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1027 "em2050",FatalException, 1092 "em2050",FatalException,
1028 "Problem in building the extended Table 1093 "Problem in building the extended Table for Sampling: array dimensions do not match ");
1029 } 1094 }
1030 1095
1031 /****************************************** 1096 /*******************************************************************************
1032 Renormalization 1097 Renormalization
1033 ******************************************* 1098 ********************************************************************************/
1034 G4double ws = 0; 1099 G4double ws = 0;
1035 for (std::size_t i=0;i<np-1;++i) << 1100 for (size_t i=0;i<np-1;i++)
1036 ws += (*area)[i]; 1101 ws += (*area)[i];
1037 ws = 1.0/ws; 1102 ws = 1.0/ws;
1038 G4double errMax = 0; 1103 G4double errMax = 0;
1039 for (std::size_t i=0;i<np-1;++i) << 1104 for (size_t i=0;i<np-1;i++)
1040 { 1105 {
1041 (*area)[i] *= ws; 1106 (*area)[i] *= ws;
1042 (*err)[i] *= ws; 1107 (*err)[i] *= ws;
1043 errMax = std::max(errMax,(*err)[i]); 1108 errMax = std::max(errMax,(*err)[i]);
1044 } 1109 }
1045 1110
1046 //Vector with the normalized cumulative dis 1111 //Vector with the normalized cumulative distribution
1047 G4DataVector* PAC = new G4DataVector(); 1112 G4DataVector* PAC = new G4DataVector();
1048 PAC->push_back(0.); 1113 PAC->push_back(0.);
1049 for (std::size_t i=0;i<np-1;++i) << 1114 for (size_t i=0;i<np-1;i++)
1050 { 1115 {
1051 G4double previous = (*PAC)[i]; 1116 G4double previous = (*PAC)[i];
1052 PAC->push_back(previous+(*area)[i]); 1117 PAC->push_back(previous+(*area)[i]);
1053 } 1118 }
1054 (*PAC)[PAC->size()-1] = 1.; 1119 (*PAC)[PAC->size()-1] = 1.;
1055 << 1120
1056 /****************************************** 1121 /*******************************************************************************
1057 Pre-calculated limits for the initial binar 1122 Pre-calculated limits for the initial binary search for subsequent sampling
1058 ******************************************* 1123 ********************************************************************************/
1059 std::vector<std::size_t> *ITTL = new std::v << 1124
1060 std::vector<std::size_t> *ITTU = new std::v << 1125 //G4DataVector* ITTL = new G4DataVector();
>> 1126 std::vector<size_t> *ITTL = new std::vector<size_t>;
>> 1127 std::vector<size_t> *ITTU = new std::vector<size_t>;
1061 1128
1062 //Just create place-holders 1129 //Just create place-holders
1063 for (std::size_t i=0;i<np;++i) << 1130 for (size_t i=0;i<np;i++)
1064 { 1131 {
1065 ITTL->push_back(0); 1132 ITTL->push_back(0);
1066 ITTU->push_back(0); 1133 ITTU->push_back(0);
1067 } 1134 }
1068 1135
1069 G4double bin = 1.0/(np-1); 1136 G4double bin = 1.0/(np-1);
1070 (*ITTL)[0]=0; 1137 (*ITTL)[0]=0;
1071 for (std::size_t i=1;i<(np-1);++i) << 1138 for (size_t i=1;i<(np-1);i++)
1072 { 1139 {
1073 G4double ptst = i*bin; << 1140 G4double ptst = i*bin;
1074 G4bool found = false; 1141 G4bool found = false;
1075 for (std::size_t j=(*ITTL)[i-1];j<np && << 1142 for (size_t j=(*ITTL)[i-1];j<np && !found;j++)
1076 { 1143 {
1077 if ((*PAC)[j] > ptst) 1144 if ((*PAC)[j] > ptst)
1078 { 1145 {
1079 (*ITTL)[i] = j-1; 1146 (*ITTL)[i] = j-1;
1080 (*ITTU)[i-1] = j; 1147 (*ITTU)[i-1] = j;
1081 found = true; 1148 found = true;
1082 } 1149 }
1083 } 1150 }
1084 } 1151 }
1085 (*ITTU)[ITTU->size()-2] = ITTU->size()-1; 1152 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1086 (*ITTU)[ITTU->size()-1] = ITTU->size()-1; 1153 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1087 (*ITTL)[ITTL->size()-1] = ITTU->size()-2; 1154 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1088 1155
1089 if (ITTU->size() != np || ITTU->size() != n 1156 if (ITTU->size() != np || ITTU->size() != np)
1090 { 1157 {
1091 G4Exception("G4PenelopeRayleighModel::I 1158 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1092 "em2051",FatalException, 1159 "em2051",FatalException,
1093 "Problem in building the Limit Tables f 1160 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1094 } 1161 }
1095 1162
>> 1163
1096 /****************************************** 1164 /********************************************************************************
1097 Copy tables 1165 Copy tables
1098 ******************************************* 1166 ********************************************************************************/
1099 G4PenelopeSamplingData* theTable = new G4Pe 1167 G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
1100 for (std::size_t i=0;i<np;++i) << 1168 for (size_t i=0;i<np;i++)
1101 { 1169 {
1102 theTable->AddPoint((*x)[i],(*PAC)[i],(* 1170 theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1103 } 1171 }
1104 1172
1105 if (fVerboseLevel > 2) << 1173 if (verboseLevel > 2)
1106 { 1174 {
1107 G4cout << "**************************** << 1175 G4cout << "*************************************************************************" <<
1108 G4endl; 1176 G4endl;
1109 G4cout << "Sampling table for Penelope 1177 G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
1110 theTable->DumpTable(); 1178 theTable->DumpTable();
1111 } << 1179 }
1112 fSamplingTable->insert(std::make_pair(mat,t << 1180 samplingTable->insert(std::make_pair(mat,theTable));
1113 1181
>> 1182
1114 //Clean up temporary vectors 1183 //Clean up temporary vectors
1115 delete x; 1184 delete x;
1116 delete a; 1185 delete a;
1117 delete b; 1186 delete b;
1118 delete c; 1187 delete c;
1119 delete err; 1188 delete err;
1120 delete area; 1189 delete area;
1121 delete PAC; 1190 delete PAC;
1122 delete ITTL; 1191 delete ITTL;
1123 delete ITTU; 1192 delete ITTU;
1124 1193
1125 //DONE! 1194 //DONE!
1126 return; 1195 return;
>> 1196
1127 } 1197 }
1128 1198
1129 //....oooOO0OOooo........oooOO0OOooo........o 1199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1130 1200
1131 void G4PenelopeRayleighModel::GetPMaxTable(co 1201 void G4PenelopeRayleighModel::GetPMaxTable(const G4Material* mat)
1132 { 1202 {
1133 if (!fPMaxTable) << 1203
>> 1204 if (!pMaxTable)
1134 { 1205 {
1135 G4cout << "G4PenelopeRayleighModel::Bui 1206 G4cout << "G4PenelopeRayleighModel::BuildPMaxTable" << G4endl;
1136 G4cout << "Going to instanziate the fPM << 1207 G4cout << "Going to instanziate the pMaxTable !" << G4endl;
1137 G4cout << "That should _not_ be here! " 1208 G4cout << "That should _not_ be here! " << G4endl;
1138 fPMaxTable = new std::map<const G4Mater << 1209 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
1139 } 1210 }
1140 //check if the table is already there 1211 //check if the table is already there
1141 if (fPMaxTable->count(mat)) << 1212 if (pMaxTable->count(mat))
1142 return; 1213 return;
1143 1214
1144 //otherwise build it 1215 //otherwise build it
1145 if (!fSamplingTable) << 1216 if (!samplingTable)
1146 { 1217 {
1147 G4Exception("G4PenelopeRayleighModel::G 1218 G4Exception("G4PenelopeRayleighModel::GetPMaxTable()",
1148 "em2052",FatalException, 1219 "em2052",FatalException,
1149 "SamplingTable is not properly instanti 1220 "SamplingTable is not properly instantiated");
1150 return; 1221 return;
1151 } 1222 }
1152 1223
1153 //This should not be: the sampling table is 1224 //This should not be: the sampling table is built before the p-table
1154 if (!fSamplingTable->count(mat)) << 1225 if (!samplingTable->count(mat))
1155 { 1226 {
1156 G4ExceptionDescription ed; 1227 G4ExceptionDescription ed;
1157 ed << "Sampling table for material " < 1228 ed << "Sampling table for material " << mat->GetName() << " not found";
1158 G4Exception("G4PenelopeRayleighModel:: 1229 G4Exception("G4PenelopeRayleighModel::GetPMaxTable()",
1159 "em2052",FatalException, 1230 "em2052",FatalException,
1160 ed); 1231 ed);
1161 return; 1232 return;
1162 } 1233 }
1163 1234
1164 G4PenelopeSamplingData *theTable = fSamplin << 1235 G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second;
1165 std::size_t tablePoints = theTable->GetNumb << 1236 size_t tablePoints = theTable->GetNumberOfStoredPoints();
1166 1237
1167 std::size_t nOfEnergyPoints = fLogEnergyGri << 1238 size_t nOfEnergyPoints = logEnergyGridPMax.size();
1168 G4PhysicsFreeVector* theVec = new G4Physics 1239 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
1169 1240
1170 const std::size_t nip = 51; //hard-coded in << 1241 const size_t nip = 51; //hard-coded in Penelope
1171 1242
1172 for (std::size_t ie=0;ie<fLogEnergyGridPMax << 1243 for (size_t ie=0;ie<logEnergyGridPMax.size();ie++)
1173 { 1244 {
1174 G4double energy = G4Exp(fLogEnergyGridP << 1245 G4double energy = std::exp(logEnergyGridPMax[ie]);
1175 G4double Qm = 2.0*energy/electron_mass_ 1246 G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
1176 G4double Qm2 = Qm*Qm; 1247 G4double Qm2 = Qm*Qm;
1177 G4double firstQ2 = theTable->GetX(0); 1248 G4double firstQ2 = theTable->GetX(0);
1178 G4double lastQ2 = theTable->GetX(tableP 1249 G4double lastQ2 = theTable->GetX(tablePoints-1);
1179 G4double thePMax = 0; 1250 G4double thePMax = 0;
1180 << 1251
1181 if (Qm2 > firstQ2) 1252 if (Qm2 > firstQ2)
1182 { 1253 {
1183 if (Qm2 < lastQ2) 1254 if (Qm2 < lastQ2)
1184 { 1255 {
1185 //bisection to look for the index of 1256 //bisection to look for the index of Qm
1186 std::size_t lowerBound = 0; << 1257 size_t lowerBound = 0;
1187 std::size_t upperBound = tablePoints- << 1258 size_t upperBound = tablePoints-1;
1188 while (lowerBound <= upperBound) 1259 while (lowerBound <= upperBound)
1189 { 1260 {
1190 std::size_t midBin = (lowerBound + uppe << 1261 size_t midBin = (lowerBound + upperBound)/2;
1191 if( Qm2 < theTable->GetX(midBin)) 1262 if( Qm2 < theTable->GetX(midBin))
1192 { upperBound = midBin-1; } 1263 { upperBound = midBin-1; }
1193 else 1264 else
1194 { lowerBound = midBin+1; } 1265 { lowerBound = midBin+1; }
1195 } 1266 }
1196 //upperBound is the output (but also 1267 //upperBound is the output (but also lowerBounf --> should be the same!)
1197 G4double Q1 = theTable->GetX(upperBou 1268 G4double Q1 = theTable->GetX(upperBound);
1198 G4double Q2 = Qm2; 1269 G4double Q2 = Qm2;
1199 G4double DQ = (Q2-Q1)/((G4double)(nip 1270 G4double DQ = (Q2-Q1)/((G4double)(nip-1));
1200 G4double theA = theTable->GetA(upperB 1271 G4double theA = theTable->GetA(upperBound);
1201 G4double theB = theTable->GetB(upperB 1272 G4double theB = theTable->GetB(upperBound);
1202 G4double thePAC = theTable->GetPAC(up 1273 G4double thePAC = theTable->GetPAC(upperBound);
1203 G4DataVector* fun = new G4DataVector( 1274 G4DataVector* fun = new G4DataVector();
1204 for (std::size_t k=0;k<nip;++k) << 1275 for (size_t k=0;k<nip;k++)
1205 { 1276 {
1206 G4double qi = Q1 + k*DQ; 1277 G4double qi = Q1 + k*DQ;
1207 G4double tau = (qi-Q1)/ 1278 G4double tau = (qi-Q1)/
1208 (theTable->GetX(upperBound+1)-Q1); 1279 (theTable->GetX(upperBound+1)-Q1);
1209 G4double con1 = 2.0*theB*tau; << 1280 G4double con1 = 2.0*theB*tau;
1210 G4double ci = 1.0+theA+theB; 1281 G4double ci = 1.0+theA+theB;
1211 G4double con2 = ci-theA*tau; 1282 G4double con2 = ci-theA*tau;
1212 G4double etap = 0; 1283 G4double etap = 0;
1213 if (std::fabs(con1) > 1.0e-16*std::fabs 1284 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1214 etap = con2*(1.0-std::sqrt(1.0-2.0*ta 1285 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1215 else 1286 else
1216 etap = tau/con2; 1287 etap = tau/con2;
1217 G4double theFun = (theTable->GetPAC(upp 1288 G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
1218 (1.0+(theA+theB*etap)*etap)*(1.0+(the 1289 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1219 ((1.0-theB*etap*etap)*ci*(theTable->G 1290 ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
1220 fun->push_back(theFun); 1291 fun->push_back(theFun);
1221 } 1292 }
1222 //Now intergrate numerically the fun 1293 //Now intergrate numerically the fun Cavalieri-Simpson's method
1223 G4DataVector* sum = new G4DataVector; 1294 G4DataVector* sum = new G4DataVector;
1224 G4double CONS = DQ*(1./12.); 1295 G4double CONS = DQ*(1./12.);
1225 G4double HCONS = 0.5*CONS; 1296 G4double HCONS = 0.5*CONS;
1226 sum->push_back(0.); 1297 sum->push_back(0.);
1227 G4double secondPoint = (*sum)[0] + << 1298 G4double secondPoint = (*sum)[0] +
1228 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C 1299 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1229 sum->push_back(secondPoint); 1300 sum->push_back(secondPoint);
1230 for (std::size_t hh=2;hh<nip-1;++hh) << 1301 for (size_t hh=2;hh<nip-1;hh++)
1231 { 1302 {
1232 G4double previous = (*sum)[hh-1]; 1303 G4double previous = (*sum)[hh-1];
1233 G4double next = previous+(13.0*((*fun)[ 1304 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1234 (*fun)[hh+1]-(*fun)[hh-2])*HCON 1305 (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1235 sum->push_back(next); 1306 sum->push_back(next);
1236 } 1307 }
1237 G4double last = (*sum)[nip-2]+(5.0*(* 1308 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1238 (*fun)[nip-3])*CONS; 1309 (*fun)[nip-3])*CONS;
1239 sum->push_back(last); << 1310 sum->push_back(last);
1240 thePMax = thePAC + (*sum)[sum->size() 1311 thePMax = thePAC + (*sum)[sum->size()-1]; //last point
1241 delete fun; 1312 delete fun;
1242 delete sum; 1313 delete sum;
1243 } 1314 }
1244 else 1315 else
1245 { 1316 {
1246 thePMax = 1.0; 1317 thePMax = 1.0;
1247 } << 1318 }
1248 } 1319 }
1249 else 1320 else
1250 { 1321 {
1251 thePMax = theTable->GetPAC(0); 1322 thePMax = theTable->GetPAC(0);
1252 } 1323 }
1253 1324
1254 //Write number in the table 1325 //Write number in the table
1255 theVec->PutValue(ie,energy,thePMax); 1326 theVec->PutValue(ie,energy,thePMax);
1256 } 1327 }
1257 << 1328
1258 fPMaxTable->insert(std::make_pair(mat,theVe << 1329 pMaxTable->insert(std::make_pair(mat,theVec));
1259 return; 1330 return;
>> 1331
1260 } 1332 }
1261 1333
1262 //....oooOO0OOooo........oooOO0OOooo........o 1334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1263 1335
1264 void G4PenelopeRayleighModel::DumpFormFactorT 1336 void G4PenelopeRayleighModel::DumpFormFactorTable(const G4Material* mat)
1265 { 1337 {
1266 G4cout << "******************************** 1338 G4cout << "*****************************************************************" << G4endl;
1267 G4cout << "G4PenelopeRayleighModel: Form Fa 1339 G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
1268 //try to use the same format as Penelope-Fo 1340 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1269 G4cout << "Q/(m_e*c) F(Q) 1341 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1270 G4cout << "******************************** 1342 G4cout << "*****************************************************************" << G4endl;
1271 if (!fLogFormFactorTable->count(mat)) << 1343 if (!logFormFactorTable->count(mat))
1272 BuildFormFactorTable(mat); 1344 BuildFormFactorTable(mat);
1273 << 1345
1274 G4PhysicsFreeVector* theVec = fLogFormFacto << 1346 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
1275 for (std::size_t i=0;i<theVec->GetVectorLen << 1347 for (size_t i=0;i<theVec->GetVectorLength();i++)
1276 { 1348 {
1277 G4double logQ2 = theVec->GetLowEdgeEner 1349 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1278 G4double Q = G4Exp(0.5*logQ2); << 1350 G4double Q = std::exp(0.5*logQ2);
1279 G4double logF2 = (*theVec)[i]; 1351 G4double logF2 = (*theVec)[i];
1280 G4double F = G4Exp(0.5*logF2); << 1352 G4double F = std::exp(0.5*logF2);
1281 G4cout << Q << " " << F << 1353 G4cout << Q << " " << F << G4endl;
1282 } 1354 }
1283 //DONE 1355 //DONE
1284 return; 1356 return;
1285 } 1357 }
1286 1358
1287 //....oooOO0OOooo........oooOO0OOooo........o 1359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1288 1360
1289 void G4PenelopeRayleighModel::SetParticle(con 1361 void G4PenelopeRayleighModel::SetParticle(const G4ParticleDefinition* p)
1290 { 1362 {
1291 if(!fParticle) { 1363 if(!fParticle) {
1292 fParticle = p; << 1364 fParticle = p;
1293 } 1365 }
1294 } 1366 }
1295 1367