1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: G4Penelope08RayleighModel.cc,v 1.3 2010/07/28 07:09:16 pandola Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-04-ref-00 $ |
---|
28 | // |
---|
29 | // Author: Luciano Pandola |
---|
30 | // |
---|
31 | // History: |
---|
32 | // -------- |
---|
33 | // 03 Dec 2009 L Pandola First implementation |
---|
34 | // |
---|
35 | |
---|
36 | #include "G4Penelope08RayleighModel.hh" |
---|
37 | #include "G4PenelopeSamplingData.hh" |
---|
38 | #include "G4ParticleDefinition.hh" |
---|
39 | #include "G4MaterialCutsCouple.hh" |
---|
40 | #include "G4ProductionCutsTable.hh" |
---|
41 | #include "G4DynamicParticle.hh" |
---|
42 | #include "G4PhysicsTable.hh" |
---|
43 | #include "G4ElementTable.hh" |
---|
44 | #include "G4Element.hh" |
---|
45 | #include "G4PhysicsFreeVector.hh" |
---|
46 | |
---|
47 | |
---|
48 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
49 | |
---|
50 | |
---|
51 | G4Penelope08RayleighModel::G4Penelope08RayleighModel(const G4ParticleDefinition*, |
---|
52 | const G4String& nam) |
---|
53 | :G4VEmModel(nam),isInitialised(false),logAtomicCrossSection(0), |
---|
54 | atomicFormFactor(0),logFormFactorTable(0),pMaxTable(0),samplingTable(0) |
---|
55 | { |
---|
56 | fIntrinsicLowEnergyLimit = 100.0*eV; |
---|
57 | fIntrinsicHighEnergyLimit = 100.0*GeV; |
---|
58 | // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); |
---|
59 | SetHighEnergyLimit(fIntrinsicHighEnergyLimit); |
---|
60 | // |
---|
61 | verboseLevel= 0; |
---|
62 | // Verbosity scale: |
---|
63 | // 0 = nothing |
---|
64 | // 1 = warning for energy non-conservation |
---|
65 | // 2 = details of energy budget |
---|
66 | // 3 = calculation of cross sections, file openings, sampling of atoms |
---|
67 | // 4 = entering in methods |
---|
68 | |
---|
69 | //build the energy grid. It is the same for all materials |
---|
70 | G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.); |
---|
71 | G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit); |
---|
72 | //finer grid below 160 keV |
---|
73 | G4double logtransitionenergy = std::log(160*keV); |
---|
74 | G4double logfactor1 = std::log(10.)/250.; |
---|
75 | G4double logfactor2 = logfactor1*10; |
---|
76 | logEnergyGridPMax.push_back(logenergy); |
---|
77 | do{ |
---|
78 | if (logenergy < logtransitionenergy) |
---|
79 | logenergy += logfactor1; |
---|
80 | else |
---|
81 | logenergy += logfactor2; |
---|
82 | logEnergyGridPMax.push_back(logenergy); |
---|
83 | }while (logenergy < logmaxenergy); |
---|
84 | } |
---|
85 | |
---|
86 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
87 | |
---|
88 | G4Penelope08RayleighModel::~G4Penelope08RayleighModel() |
---|
89 | { |
---|
90 | std::map <const G4int,G4PhysicsFreeVector*>::iterator i; |
---|
91 | if (logAtomicCrossSection) |
---|
92 | { |
---|
93 | for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++) |
---|
94 | if (i->second) delete i->second; |
---|
95 | delete logAtomicCrossSection; |
---|
96 | } |
---|
97 | |
---|
98 | if (atomicFormFactor) |
---|
99 | { |
---|
100 | for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++) |
---|
101 | if (i->second) delete i->second; |
---|
102 | delete atomicFormFactor; |
---|
103 | } |
---|
104 | |
---|
105 | ClearTables(); |
---|
106 | } |
---|
107 | |
---|
108 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
109 | void G4Penelope08RayleighModel::ClearTables() |
---|
110 | { |
---|
111 | std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i; |
---|
112 | |
---|
113 | if (logFormFactorTable) |
---|
114 | { |
---|
115 | for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++) |
---|
116 | if (i->second) delete i->second; |
---|
117 | delete logFormFactorTable; |
---|
118 | logFormFactorTable = 0; //zero explicitely |
---|
119 | } |
---|
120 | |
---|
121 | if (pMaxTable) |
---|
122 | { |
---|
123 | for (i=pMaxTable->begin(); i != pMaxTable->end(); i++) |
---|
124 | if (i->second) delete i->second; |
---|
125 | delete pMaxTable; |
---|
126 | pMaxTable = 0; //zero explicitely |
---|
127 | } |
---|
128 | |
---|
129 | std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii; |
---|
130 | if (samplingTable) |
---|
131 | { |
---|
132 | for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++) |
---|
133 | if (ii->second) delete ii->second; |
---|
134 | delete samplingTable; |
---|
135 | samplingTable = 0; //zero explicitely |
---|
136 | } |
---|
137 | |
---|
138 | return; |
---|
139 | } |
---|
140 | |
---|
141 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
142 | |
---|
143 | void G4Penelope08RayleighModel::Initialise(const G4ParticleDefinition* , |
---|
144 | const G4DataVector& ) |
---|
145 | { |
---|
146 | if (verboseLevel > 3) |
---|
147 | G4cout << "Calling G4Penelope08RayleighModel::Initialise()" << G4endl; |
---|
148 | |
---|
149 | //clear tables depending on materials, not the atomic ones |
---|
150 | ClearTables(); |
---|
151 | |
---|
152 | //create new tables |
---|
153 | // |
---|
154 | // logAtomicCrossSection and atomicFormFactor are created only once, |
---|
155 | // since they are never cleared |
---|
156 | if (!logAtomicCrossSection) |
---|
157 | logAtomicCrossSection = new std::map<const G4int,G4PhysicsFreeVector*>; |
---|
158 | if (!atomicFormFactor) |
---|
159 | atomicFormFactor = new std::map<const G4int,G4PhysicsFreeVector*>; |
---|
160 | |
---|
161 | if (!logFormFactorTable) |
---|
162 | logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
163 | if (!pMaxTable) |
---|
164 | pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
165 | if (!samplingTable) |
---|
166 | samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>; |
---|
167 | |
---|
168 | |
---|
169 | if (verboseLevel > 0) { |
---|
170 | G4cout << "Penelope08 Rayleigh model is initialized " << G4endl |
---|
171 | << "Energy range: " |
---|
172 | << LowEnergyLimit() / keV << " keV - " |
---|
173 | << HighEnergyLimit() / GeV << " GeV" |
---|
174 | << G4endl; |
---|
175 | } |
---|
176 | |
---|
177 | if(isInitialised) return; |
---|
178 | fParticleChange = GetParticleChangeForGamma(); |
---|
179 | isInitialised = true; |
---|
180 | } |
---|
181 | |
---|
182 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
183 | |
---|
184 | G4double G4Penelope08RayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, |
---|
185 | G4double energy, |
---|
186 | G4double Z, |
---|
187 | G4double, |
---|
188 | G4double, |
---|
189 | G4double) |
---|
190 | { |
---|
191 | // Cross section of Rayleigh scattering in Penelope2008 is calculated by the EPDL97 |
---|
192 | // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel |
---|
193 | // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615. |
---|
194 | |
---|
195 | if (verboseLevel > 3) |
---|
196 | G4cout << "Calling CrossSectionPerAtom() of G4Penelope08RayleighModel" << G4endl; |
---|
197 | |
---|
198 | G4int iZ = (G4int) Z; |
---|
199 | |
---|
200 | //read data files |
---|
201 | if (!logAtomicCrossSection->count(iZ)) |
---|
202 | ReadDataFile(iZ); |
---|
203 | //now it should be ok |
---|
204 | if (!logAtomicCrossSection->count(iZ)) |
---|
205 | { |
---|
206 | G4cout << "Problem in G4Penelope08RayleighModel::ComputeCrossSectionPerAtom" |
---|
207 | << G4endl; |
---|
208 | G4Exception(); |
---|
209 | } |
---|
210 | |
---|
211 | G4double cross = 0; |
---|
212 | |
---|
213 | G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second; |
---|
214 | if (!atom) |
---|
215 | { |
---|
216 | G4cout << "Problem in G4Penelope08RayleighModel::ComputeCrossSectionPerAtom" |
---|
217 | << G4endl; |
---|
218 | G4Exception(); |
---|
219 | } |
---|
220 | G4double logene = std::log(energy); |
---|
221 | G4double logXS = atom->Value(logene); |
---|
222 | cross = std::exp(logXS); |
---|
223 | |
---|
224 | if (verboseLevel > 2) |
---|
225 | G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z << |
---|
226 | " = " << cross/barn << " barn" << G4endl; |
---|
227 | return cross; |
---|
228 | } |
---|
229 | |
---|
230 | |
---|
231 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
232 | void G4Penelope08RayleighModel::BuildFormFactorTable(const G4Material* material) |
---|
233 | { |
---|
234 | |
---|
235 | /* |
---|
236 | 1) get composition and equivalent molecular density |
---|
237 | */ |
---|
238 | |
---|
239 | G4int nElements = material->GetNumberOfElements(); |
---|
240 | const G4ElementVector* elementVector = material->GetElementVector(); |
---|
241 | const G4double* fractionVector = material->GetFractionVector(); |
---|
242 | |
---|
243 | std::vector<G4double> *StechiometricFactors = new std::vector<G4double>; |
---|
244 | for (G4int i=0;i<nElements;i++) |
---|
245 | { |
---|
246 | G4double fraction = fractionVector[i]; |
---|
247 | G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole); |
---|
248 | StechiometricFactors->push_back(fraction/atomicWeigth); |
---|
249 | } |
---|
250 | //Find max |
---|
251 | G4double MaxStechiometricFactor = 0.; |
---|
252 | for (G4int i=0;i<nElements;i++) |
---|
253 | { |
---|
254 | if ((*StechiometricFactors)[i] > MaxStechiometricFactor) |
---|
255 | MaxStechiometricFactor = (*StechiometricFactors)[i]; |
---|
256 | } |
---|
257 | if (MaxStechiometricFactor<1e-16) |
---|
258 | { |
---|
259 | G4cout << "G4Penelope08RayleighModel::BuildFormFactorTable" << G4endl; |
---|
260 | G4cout << "Problem with the mass composition of " << material->GetName() << G4endl; |
---|
261 | G4Exception(); |
---|
262 | } |
---|
263 | //Normalize |
---|
264 | for (G4int i=0;i<nElements;i++) |
---|
265 | (*StechiometricFactors)[i] /= MaxStechiometricFactor; |
---|
266 | |
---|
267 | // Equivalent atoms per molecule |
---|
268 | G4double atomsPerMolecule = 0; |
---|
269 | for (G4int i=0;i<nElements;i++) |
---|
270 | atomsPerMolecule += (*StechiometricFactors)[i]; |
---|
271 | |
---|
272 | /* |
---|
273 | CREATE THE FORM FACTOR TABLE |
---|
274 | */ |
---|
275 | G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size()); |
---|
276 | theFFVec->SetSpline(true); |
---|
277 | |
---|
278 | for (size_t k=0;k<logQSquareGrid.size();k++) |
---|
279 | { |
---|
280 | G4double ff2 = 0; //squared form factor |
---|
281 | for (G4int i=0;i<nElements;i++) |
---|
282 | { |
---|
283 | G4int iZ = (G4int) (*elementVector)[i]->GetZ(); |
---|
284 | G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second; |
---|
285 | G4double f = (*theAtomVec)[k]; //the q-grid is always the same |
---|
286 | ff2 += f*f*(*StechiometricFactors)[i]; |
---|
287 | } |
---|
288 | if (ff2) |
---|
289 | theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2) |
---|
290 | } |
---|
291 | logFormFactorTable->insert(std::make_pair(material,theFFVec)); |
---|
292 | |
---|
293 | delete StechiometricFactors; |
---|
294 | |
---|
295 | return; |
---|
296 | } |
---|
297 | |
---|
298 | |
---|
299 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
300 | |
---|
301 | void G4Penelope08RayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* , |
---|
302 | const G4MaterialCutsCouple* couple, |
---|
303 | const G4DynamicParticle* aDynamicGamma, |
---|
304 | G4double, |
---|
305 | G4double) |
---|
306 | { |
---|
307 | // Sampling of the Rayleigh final state (namely, scattering angle of the photon) |
---|
308 | // from the Penelope2008 model. The scattering angle is sampled from the atomic |
---|
309 | // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding |
---|
310 | // anomalous scattering effects. The Form Factor F(Q) function which appears in the |
---|
311 | // analytical cross section is retrieved via the method GetFSquared(); atomic data |
---|
312 | // are tabulated for F(Q). Form factor for compounds is calculated according to |
---|
313 | // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse |
---|
314 | // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once |
---|
315 | // for each material and managed by G4PenelopeSamplingData objects. |
---|
316 | // The sampling algorithm (rejection method) has efficiency 67% at low energy, and |
---|
317 | // increases with energy. For E=100 keV the efficiency is 100% and 86% for |
---|
318 | // hydrogen and uranium, respectively. |
---|
319 | |
---|
320 | if (verboseLevel > 3) |
---|
321 | G4cout << "Calling SamplingSecondaries() of G4Penelope08RayleighModel" << G4endl; |
---|
322 | |
---|
323 | G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); |
---|
324 | |
---|
325 | if (photonEnergy0 <= fIntrinsicLowEnergyLimit) |
---|
326 | { |
---|
327 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
328 | fParticleChange->SetProposedKineticEnergy(0.); |
---|
329 | fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); |
---|
330 | return ; |
---|
331 | } |
---|
332 | |
---|
333 | G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); |
---|
334 | |
---|
335 | const G4Material* theMat = couple->GetMaterial(); |
---|
336 | |
---|
337 | //1) Verify if tables are ready |
---|
338 | if (!pMaxTable || !samplingTable) |
---|
339 | { |
---|
340 | G4cout << " G4Penelope08RayleighModel::SampleSecondaries" << G4endl; |
---|
341 | G4cout << "Looks like the model is _not_ properly initialized" << G4endl; |
---|
342 | G4Exception(); |
---|
343 | } |
---|
344 | |
---|
345 | //2) retrieve or build the sampling table |
---|
346 | if (!(samplingTable->count(theMat))) |
---|
347 | InitializeSamplingAlgorithm(theMat); |
---|
348 | G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second; |
---|
349 | |
---|
350 | //3) retrieve or build the pMax data |
---|
351 | if (!pMaxTable->count(theMat)) |
---|
352 | GetPMaxTable(theMat); |
---|
353 | G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second; |
---|
354 | |
---|
355 | G4double cosTheta = 1.0; |
---|
356 | |
---|
357 | //OK, ready to go! |
---|
358 | G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now |
---|
359 | |
---|
360 | if (qmax < 1e-10) //very low momentum transfer |
---|
361 | { |
---|
362 | G4bool loopAgain=false; |
---|
363 | do |
---|
364 | { |
---|
365 | loopAgain = false; |
---|
366 | cosTheta = 1.0-2.0*G4UniformRand(); |
---|
367 | G4double G = 0.5*(1+cosTheta*cosTheta); |
---|
368 | if (G4UniformRand()>G) |
---|
369 | loopAgain = true; |
---|
370 | }while(loopAgain); |
---|
371 | } |
---|
372 | else //larger momentum transfer |
---|
373 | { |
---|
374 | size_t nData = theDataTable->GetNumberOfStoredPoints(); |
---|
375 | G4double LastQ2inTheTable = theDataTable->GetX(nData-1); |
---|
376 | G4double q2max = std::min(qmax*qmax,LastQ2inTheTable); |
---|
377 | |
---|
378 | G4bool loopAgain = false; |
---|
379 | G4double MaxPValue = thePMax->Value(photonEnergy0); |
---|
380 | G4double xx=0; |
---|
381 | |
---|
382 | //Sampling by rejection method. The rejection function is |
---|
383 | //G = 0.5*(1+cos^2(theta)) |
---|
384 | // |
---|
385 | do{ |
---|
386 | loopAgain = false; |
---|
387 | G4double RandomMax = G4UniformRand()*MaxPValue; |
---|
388 | xx = theDataTable->SampleValue(RandomMax); |
---|
389 | //xx is a random value of q^2 in (0,q2max),sampled according to |
---|
390 | //F(Q^2) via the RITA algorithm |
---|
391 | if (xx > q2max) |
---|
392 | loopAgain = true; |
---|
393 | cosTheta = 1.0-2.0*xx/q2max; |
---|
394 | G4double G = 0.5*(1+cosTheta*cosTheta); |
---|
395 | if (G4UniformRand()>G) |
---|
396 | loopAgain = true; |
---|
397 | }while(loopAgain); |
---|
398 | } |
---|
399 | |
---|
400 | G4double sinTheta = std::sqrt(1-cosTheta*cosTheta); |
---|
401 | |
---|
402 | // Scattered photon angles. ( Z - axis along the parent photon) |
---|
403 | G4double phi = twopi * G4UniformRand() ; |
---|
404 | G4double dirX = sinTheta*std::cos(phi); |
---|
405 | G4double dirY = sinTheta*std::sin(phi); |
---|
406 | G4double dirZ = cosTheta; |
---|
407 | |
---|
408 | // Update G4VParticleChange for the scattered photon |
---|
409 | G4ThreeVector photonDirection1(dirX, dirY, dirZ); |
---|
410 | |
---|
411 | photonDirection1.rotateUz(photonDirection0); |
---|
412 | fParticleChange->ProposeMomentumDirection(photonDirection1) ; |
---|
413 | fParticleChange->SetProposedKineticEnergy(photonEnergy0) ; |
---|
414 | |
---|
415 | return; |
---|
416 | } |
---|
417 | |
---|
418 | |
---|
419 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
420 | |
---|
421 | void G4Penelope08RayleighModel::ReadDataFile(const G4int Z) |
---|
422 | { |
---|
423 | if (verboseLevel > 2) |
---|
424 | { |
---|
425 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
426 | G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl; |
---|
427 | } |
---|
428 | |
---|
429 | char* path = getenv("G4LEDATA"); |
---|
430 | if (!path) |
---|
431 | { |
---|
432 | G4String excep = "G4Penelope08RayleighModel - G4LEDATA environment variable not set!"; |
---|
433 | G4Exception(excep); |
---|
434 | } |
---|
435 | |
---|
436 | /* |
---|
437 | Read first the cross section file |
---|
438 | */ |
---|
439 | std::ostringstream ost; |
---|
440 | if (Z>9) |
---|
441 | ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08"; |
---|
442 | else |
---|
443 | ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08"; |
---|
444 | std::ifstream file(ost.str().c_str()); |
---|
445 | if (!file.is_open()) |
---|
446 | { |
---|
447 | G4String excep = "G4Penelope08RayleighModel - data file " + G4String(ost.str()) + " not found!"; |
---|
448 | G4Exception(excep); |
---|
449 | } |
---|
450 | G4int readZ =0; |
---|
451 | size_t nPoints= 0; |
---|
452 | file >> readZ >> nPoints; |
---|
453 | //check the right file is opened. |
---|
454 | if (readZ != Z || nPoints <= 0) |
---|
455 | { |
---|
456 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
457 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
458 | G4Exception(); |
---|
459 | } |
---|
460 | G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints); |
---|
461 | G4double ene=0,f1=0,f2=0,xs=0; |
---|
462 | for (size_t i=0;i<nPoints;i++) |
---|
463 | { |
---|
464 | file >> ene >> f1 >> f2 >> xs; |
---|
465 | //dimensional quantities |
---|
466 | ene *= eV; |
---|
467 | xs *= cm2; |
---|
468 | theVec->PutValue(i,std::log(ene),std::log(xs)); |
---|
469 | if (file.eof() && i != (nPoints-1)) //file ended too early |
---|
470 | { |
---|
471 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
472 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
473 | G4cout << "Found less than " << nPoints << "entries " <<G4endl; |
---|
474 | G4Exception(); |
---|
475 | } |
---|
476 | } |
---|
477 | if (!logAtomicCrossSection) |
---|
478 | { |
---|
479 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
480 | G4cout << "Problem with allocation of logAtomicCrossSection data table " << G4endl; |
---|
481 | G4Exception(); |
---|
482 | } |
---|
483 | logAtomicCrossSection->insert(std::make_pair(Z,theVec)); |
---|
484 | file.close(); |
---|
485 | |
---|
486 | /* |
---|
487 | Then read the form factor file |
---|
488 | */ |
---|
489 | std::ostringstream ost2; |
---|
490 | if (Z>9) |
---|
491 | ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08"; |
---|
492 | else |
---|
493 | ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08"; |
---|
494 | file.open(ost2.str().c_str()); |
---|
495 | if (!file.is_open()) |
---|
496 | { |
---|
497 | G4String excep = "G4Penelope08RayleighModel - data file " + G4String(ost2.str()) + " not found!"; |
---|
498 | G4Exception(excep); |
---|
499 | } |
---|
500 | file >> readZ >> nPoints; |
---|
501 | //check the right file is opened. |
---|
502 | if (readZ != Z || nPoints <= 0) |
---|
503 | { |
---|
504 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
505 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
506 | G4Exception(); |
---|
507 | } |
---|
508 | G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints); |
---|
509 | G4double q=0,ff=0,incoh=0; |
---|
510 | G4bool fillQGrid = false; |
---|
511 | //fill this vector only the first time. |
---|
512 | if (!logQSquareGrid.size()) |
---|
513 | fillQGrid = true; |
---|
514 | for (size_t i=0;i<nPoints;i++) |
---|
515 | { |
---|
516 | file >> q >> ff >> incoh; |
---|
517 | //q and ff are dimensionless (q is in units of (m_e*c) |
---|
518 | theFFVec->PutValue(i,q,ff); |
---|
519 | if (fillQGrid) |
---|
520 | { |
---|
521 | logQSquareGrid.push_back(2.0*std::log(q)); |
---|
522 | } |
---|
523 | if (file.eof() && i != (nPoints-1)) //file ended too early |
---|
524 | { |
---|
525 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
526 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
527 | G4cout << "Found less than " << nPoints << "entries " <<G4endl; |
---|
528 | G4Exception(); |
---|
529 | } |
---|
530 | } |
---|
531 | if (!atomicFormFactor) |
---|
532 | { |
---|
533 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
534 | G4cout << "Problem with allocation of atomicFormFactor data table " << G4endl; |
---|
535 | G4Exception(); |
---|
536 | } |
---|
537 | atomicFormFactor->insert(std::make_pair(Z,theFFVec)); |
---|
538 | file.close(); |
---|
539 | return; |
---|
540 | } |
---|
541 | |
---|
542 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
543 | |
---|
544 | G4double G4Penelope08RayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared) |
---|
545 | { |
---|
546 | G4double f2 = 0; |
---|
547 | G4double logQSquared = std::log(QSquared); |
---|
548 | //last value of the table |
---|
549 | G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1]; |
---|
550 | //If the table has not been built for the material, do it! |
---|
551 | if (!logFormFactorTable->count(mat)) |
---|
552 | BuildFormFactorTable(mat); |
---|
553 | |
---|
554 | //now it should be all right |
---|
555 | G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second; |
---|
556 | |
---|
557 | if (!theVec) |
---|
558 | { |
---|
559 | G4cout << " G4Penelope08RayleighModel::GetFSquared()" << G4endl; |
---|
560 | G4cout << "Unable to retrieve F squared table for " << mat->GetName(); |
---|
561 | G4Exception(); |
---|
562 | } |
---|
563 | if (logQSquared < -20) // Q < 1e-9 |
---|
564 | { |
---|
565 | G4double logf2 = (*theVec)[0]; //first value of the table |
---|
566 | f2 = std::exp(logf2); |
---|
567 | } |
---|
568 | else if (logQSquared > maxlogQ2) |
---|
569 | f2 =0; |
---|
570 | else |
---|
571 | { |
---|
572 | //log(Q^2) vs. log(F^2) |
---|
573 | G4double logf2 = theVec->Value(logQSquared); |
---|
574 | f2 = std::exp(logf2); |
---|
575 | |
---|
576 | } |
---|
577 | if (verboseLevel > 3) |
---|
578 | { |
---|
579 | G4cout << "G4Penelope08RayleighModel::GetFSquared() in " << mat->GetName() << G4endl; |
---|
580 | G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl; |
---|
581 | } |
---|
582 | return f2; |
---|
583 | } |
---|
584 | |
---|
585 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
586 | |
---|
587 | void G4Penelope08RayleighModel::InitializeSamplingAlgorithm(const G4Material* mat) |
---|
588 | { |
---|
589 | G4double q2min = 0; |
---|
590 | G4double q2max = 0; |
---|
591 | const size_t np = 150; //hard-coded in Penelope |
---|
592 | for (size_t i=1;i<logQSquareGrid.size();i++) |
---|
593 | { |
---|
594 | G4double Q2 = std::exp(logQSquareGrid[i]); |
---|
595 | if (GetFSquared(mat,Q2) > 1e-35) |
---|
596 | { |
---|
597 | q2max = std::exp(logQSquareGrid[i-1]); |
---|
598 | } |
---|
599 | } |
---|
600 | |
---|
601 | size_t nReducedPoints = np/4; |
---|
602 | |
---|
603 | //check for errors |
---|
604 | if (np < 16) |
---|
605 | { |
---|
606 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
607 | G4cout << "Too few points to initialize the sampling algorithm" << G4endl; |
---|
608 | G4Exception(); |
---|
609 | } |
---|
610 | if (q2min > (q2max-1e-10)) |
---|
611 | { |
---|
612 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
613 | G4cout << "Too narrow grid to initialize the sampling algorithm" << G4endl; |
---|
614 | G4Exception(); |
---|
615 | } |
---|
616 | |
---|
617 | //This is subroutine RITAI0 of Penelope |
---|
618 | //Create an object of type G4Penelope08RayleighSamplingData --> store in a map::Material* |
---|
619 | |
---|
620 | //temporary vectors --> Then everything is stored in G4PenelopeSamplingData |
---|
621 | G4DataVector* x = new G4DataVector(); |
---|
622 | |
---|
623 | /******************************************************************************* |
---|
624 | Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max |
---|
625 | ********************************************************************************/ |
---|
626 | size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2); |
---|
627 | const G4int nip = 51; //hard-coded in Penelope |
---|
628 | |
---|
629 | G4double dx = (q2max-q2min)/((G4double) NUNIF-1); |
---|
630 | x->push_back(q2min); |
---|
631 | for (size_t i=1;i<NUNIF-1;i++) |
---|
632 | { |
---|
633 | G4double app = q2min + i*dx; |
---|
634 | x->push_back(app); //increase |
---|
635 | } |
---|
636 | x->push_back(q2max); |
---|
637 | |
---|
638 | if (verboseLevel> 3) |
---|
639 | G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl; |
---|
640 | |
---|
641 | //Allocate temporary storage vectors |
---|
642 | G4DataVector* area = new G4DataVector(); |
---|
643 | G4DataVector* a = new G4DataVector(); |
---|
644 | G4DataVector* b = new G4DataVector(); |
---|
645 | G4DataVector* c = new G4DataVector(); |
---|
646 | G4DataVector* err = new G4DataVector(); |
---|
647 | |
---|
648 | for (size_t i=0;i<NUNIF-1;i++) //build all points but the last |
---|
649 | { |
---|
650 | //Temporary vectors for this loop |
---|
651 | G4DataVector* pdfi = new G4DataVector(); |
---|
652 | G4DataVector* pdfih = new G4DataVector(); |
---|
653 | G4DataVector* sumi = new G4DataVector(); |
---|
654 | |
---|
655 | G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1)); |
---|
656 | G4double pdfmax = 0; |
---|
657 | for (G4int k=0;k<nip;k++) |
---|
658 | { |
---|
659 | G4double xik = (*x)[i]+k*dxi; |
---|
660 | G4double pdfk = std::max(GetFSquared(mat,xik),0.); |
---|
661 | pdfi->push_back(pdfk); |
---|
662 | pdfmax = std::max(pdfmax,pdfk); |
---|
663 | if (k < (nip-1)) |
---|
664 | { |
---|
665 | G4double xih = xik + 0.5*dxi; |
---|
666 | G4double pdfIK = std::max(GetFSquared(mat,xih),0.); |
---|
667 | pdfih->push_back(pdfIK); |
---|
668 | pdfmax = std::max(pdfmax,pdfIK); |
---|
669 | } |
---|
670 | } |
---|
671 | |
---|
672 | //Simpson's integration |
---|
673 | G4double cons = dxi*0.5*(1./3.); |
---|
674 | sumi->push_back(0.); |
---|
675 | for (G4int k=1;k<nip;k++) |
---|
676 | { |
---|
677 | G4double previous = (*sumi)[k-1]; |
---|
678 | G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]); |
---|
679 | sumi->push_back(next); |
---|
680 | } |
---|
681 | |
---|
682 | G4double lastIntegral = (*sumi)[sumi->size()-1]; |
---|
683 | area->push_back(lastIntegral); |
---|
684 | //Normalize cumulative function |
---|
685 | G4double factor = 1.0/lastIntegral; |
---|
686 | for (size_t k=0;k<sumi->size();k++) |
---|
687 | (*sumi)[k] *= factor; |
---|
688 | |
---|
689 | //When the PDF vanishes at one of the interval end points, its value is modified |
---|
690 | if ((*pdfi)[0] < 1e-35) |
---|
691 | (*pdfi)[0] = 1e-5*pdfmax; |
---|
692 | if ((*pdfi)[pdfi->size()-1] < 1e-35) |
---|
693 | (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; |
---|
694 | |
---|
695 | G4double pli = (*pdfi)[0]*factor; |
---|
696 | G4double pui = (*pdfi)[pdfi->size()-1]*factor; |
---|
697 | G4double B_temp = 1.0-1.0/(pli*pui*dx*dx); |
---|
698 | G4double A_temp = (1.0/(pli*dx))-1.0-B_temp; |
---|
699 | G4double C_temp = 1.0+A_temp+B_temp; |
---|
700 | if (C_temp < 1e-35) |
---|
701 | { |
---|
702 | a->push_back(0.); |
---|
703 | b->push_back(0.); |
---|
704 | c->push_back(1.); |
---|
705 | } |
---|
706 | else |
---|
707 | { |
---|
708 | a->push_back(A_temp); |
---|
709 | b->push_back(B_temp); |
---|
710 | c->push_back(C_temp); |
---|
711 | } |
---|
712 | |
---|
713 | //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation |
---|
714 | //and the true pdf, extended over the interval (X(I),X(I+1)) |
---|
715 | G4int icase = 1; //loop code |
---|
716 | G4bool reLoop = false; |
---|
717 | err->push_back(0.); |
---|
718 | do |
---|
719 | { |
---|
720 | reLoop = false; |
---|
721 | (*err)[i] = 0.; //zero variable |
---|
722 | for (G4int k=0;k<nip;k++) |
---|
723 | { |
---|
724 | G4double rr = (*sumi)[k]; |
---|
725 | G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/ |
---|
726 | ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i])); |
---|
727 | if (k == 0 || k == nip-1) |
---|
728 | (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]); |
---|
729 | else |
---|
730 | (*err)[i] += std::fabs(pap-(*pdfi)[k]); |
---|
731 | } |
---|
732 | (*err)[i] *= dxi; |
---|
733 | |
---|
734 | //If err(I) is too large, the pdf is approximated by a uniform distribution |
---|
735 | if ((*err)[i] > 0.1*(*area)[i] && icase == 1) |
---|
736 | { |
---|
737 | (*b)[i] = 0; |
---|
738 | (*a)[i] = 0; |
---|
739 | (*c)[i] = 1.; |
---|
740 | icase = 2; |
---|
741 | reLoop = true; |
---|
742 | } |
---|
743 | }while(reLoop); |
---|
744 | |
---|
745 | if (pdfi) delete pdfi; |
---|
746 | if (pdfih) delete pdfih; |
---|
747 | if (sumi) delete sumi; |
---|
748 | } //end of first loop over i |
---|
749 | |
---|
750 | //Now assign last point |
---|
751 | (*x)[x->size()-1] = q2max; |
---|
752 | a->push_back(0.); |
---|
753 | b->push_back(0.); |
---|
754 | c->push_back(0.); |
---|
755 | err->push_back(0.); |
---|
756 | area->push_back(0.); |
---|
757 | |
---|
758 | if (x->size() != NUNIF || a->size() != NUNIF || |
---|
759 | err->size() != NUNIF || area->size() != NUNIF) |
---|
760 | { |
---|
761 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
762 | G4cout << "Problem in building the Table for Sampling: array dimensions do not match " << G4endl; |
---|
763 | G4Exception(); |
---|
764 | } |
---|
765 | |
---|
766 | /******************************************************************************* |
---|
767 | New grid points are added by halving the sub-intervals with the largest absolute error |
---|
768 | This is done up to np=150 points in the grid |
---|
769 | ********************************************************************************/ |
---|
770 | do |
---|
771 | { |
---|
772 | G4double maxError = 0.0; |
---|
773 | size_t iErrMax = 0; |
---|
774 | for (size_t i=0;i<err->size()-2;i++) |
---|
775 | { |
---|
776 | //maxError is the lagest of the interval errors err[i] |
---|
777 | if ((*err)[i] > maxError) |
---|
778 | { |
---|
779 | maxError = (*err)[i]; |
---|
780 | iErrMax = i; |
---|
781 | } |
---|
782 | } |
---|
783 | |
---|
784 | //OK, now I have to insert one new point in the position iErrMax |
---|
785 | G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]); |
---|
786 | |
---|
787 | x->insert(x->begin()+iErrMax+1,newx); |
---|
788 | //Add place-holders in the other vectors |
---|
789 | area->insert(area->begin()+iErrMax+1,0.); |
---|
790 | a->insert(a->begin()+iErrMax+1,0.); |
---|
791 | b->insert(b->begin()+iErrMax+1,0.); |
---|
792 | c->insert(c->begin()+iErrMax+1,0.); |
---|
793 | err->insert(err->begin()+iErrMax+1,0.); |
---|
794 | |
---|
795 | //Now calculate the other parameters |
---|
796 | for (size_t i=iErrMax;i<=iErrMax+1;i++) |
---|
797 | { |
---|
798 | //Temporary vectors for this loop |
---|
799 | G4DataVector* pdfi = new G4DataVector(); |
---|
800 | G4DataVector* pdfih = new G4DataVector(); |
---|
801 | G4DataVector* sumi = new G4DataVector(); |
---|
802 | |
---|
803 | G4double dx = (*x)[i+1]-(*x)[i]; |
---|
804 | G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1)); |
---|
805 | G4double pdfmax = 0; |
---|
806 | for (G4int k=0;k<nip;k++) |
---|
807 | { |
---|
808 | G4double xik = (*x)[i]+k*dxi; |
---|
809 | G4double pdfk = std::max(GetFSquared(mat,xik),0.); |
---|
810 | pdfi->push_back(pdfk); |
---|
811 | pdfmax = std::max(pdfmax,pdfk); |
---|
812 | if (k < (nip-1)) |
---|
813 | { |
---|
814 | G4double xih = xik + 0.5*dxi; |
---|
815 | G4double pdfIK = std::max(GetFSquared(mat,xih),0.); |
---|
816 | pdfih->push_back(pdfIK); |
---|
817 | pdfmax = std::max(pdfmax,pdfIK); |
---|
818 | } |
---|
819 | } |
---|
820 | |
---|
821 | //Simpson's integration |
---|
822 | G4double cons = dxi*0.5*(1./3.); |
---|
823 | sumi->push_back(0.); |
---|
824 | for (G4int k=1;k<nip;k++) |
---|
825 | { |
---|
826 | G4double previous = (*sumi)[k-1]; |
---|
827 | G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]); |
---|
828 | sumi->push_back(next); |
---|
829 | } |
---|
830 | G4double lastIntegral = (*sumi)[sumi->size()-1]; |
---|
831 | (*area)[i] = lastIntegral; |
---|
832 | |
---|
833 | //Normalize cumulative function |
---|
834 | G4double factor = 1.0/lastIntegral; |
---|
835 | for (size_t k=0;k<sumi->size();k++) |
---|
836 | (*sumi)[k] *= factor; |
---|
837 | |
---|
838 | //When the PDF vanishes at one of the interval end points, its value is modified |
---|
839 | if ((*pdfi)[0] < 1e-35) |
---|
840 | (*pdfi)[0] = 1e-5*pdfmax; |
---|
841 | if ((*pdfi)[pdfi->size()-1] < 1e-35) |
---|
842 | (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; |
---|
843 | |
---|
844 | G4double pli = (*pdfi)[0]*factor; |
---|
845 | G4double pui = (*pdfi)[pdfi->size()-1]*factor; |
---|
846 | G4double B_temp = 1.0-1.0/(pli*pui*dx*dx); |
---|
847 | G4double A_temp = (1.0/(pli*dx))-1.0-B_temp; |
---|
848 | G4double C_temp = 1.0+A_temp+B_temp; |
---|
849 | if (C_temp < 1e-35) |
---|
850 | { |
---|
851 | (*a)[i]= 0.; |
---|
852 | (*b)[i] = 0.; |
---|
853 | (*c)[i] = 0; |
---|
854 | } |
---|
855 | else |
---|
856 | { |
---|
857 | (*a)[i]= A_temp; |
---|
858 | (*b)[i] = B_temp; |
---|
859 | (*c)[i] = C_temp; |
---|
860 | } |
---|
861 | //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation |
---|
862 | //and the true pdf, extended over the interval (X(I),X(I+1)) |
---|
863 | G4int icase = 1; //loop code |
---|
864 | G4bool reLoop = false; |
---|
865 | do |
---|
866 | { |
---|
867 | reLoop = false; |
---|
868 | (*err)[i] = 0.; //zero variable |
---|
869 | for (G4int k=0;k<nip;k++) |
---|
870 | { |
---|
871 | G4double rr = (*sumi)[k]; |
---|
872 | G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/ |
---|
873 | ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i])); |
---|
874 | if (k == 0 || k == nip-1) |
---|
875 | (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]); |
---|
876 | else |
---|
877 | (*err)[i] += std::fabs(pap-(*pdfi)[k]); |
---|
878 | } |
---|
879 | (*err)[i] *= dxi; |
---|
880 | |
---|
881 | //If err(I) is too large, the pdf is approximated by a uniform distribution |
---|
882 | if ((*err)[i] > 0.1*(*area)[i] && icase == 1) |
---|
883 | { |
---|
884 | (*b)[i] = 0; |
---|
885 | (*a)[i] = 0; |
---|
886 | (*c)[i] = 1.; |
---|
887 | icase = 2; |
---|
888 | reLoop = true; |
---|
889 | } |
---|
890 | }while(reLoop); |
---|
891 | if (pdfi) delete pdfi; |
---|
892 | if (pdfih) delete pdfih; |
---|
893 | if (sumi) delete sumi; |
---|
894 | } |
---|
895 | }while(x->size() < np); |
---|
896 | |
---|
897 | if (x->size() != np || a->size() != np || |
---|
898 | err->size() != np || area->size() != np) |
---|
899 | { |
---|
900 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
901 | G4cout << "Problem in building the extended Table for Sampling: array dimensions do not match " << G4endl; |
---|
902 | G4Exception(); |
---|
903 | } |
---|
904 | |
---|
905 | /******************************************************************************* |
---|
906 | Renormalization |
---|
907 | ********************************************************************************/ |
---|
908 | G4double ws = 0; |
---|
909 | for (size_t i=0;i<np-1;i++) |
---|
910 | ws += (*area)[i]; |
---|
911 | ws = 1.0/ws; |
---|
912 | G4double errMax = 0; |
---|
913 | for (size_t i=0;i<np-1;i++) |
---|
914 | { |
---|
915 | (*area)[i] *= ws; |
---|
916 | (*err)[i] *= ws; |
---|
917 | errMax = std::max(errMax,(*err)[i]); |
---|
918 | } |
---|
919 | |
---|
920 | //Vector with the normalized cumulative distribution |
---|
921 | G4DataVector* PAC = new G4DataVector(); |
---|
922 | PAC->push_back(0.); |
---|
923 | for (size_t i=0;i<np-1;i++) |
---|
924 | { |
---|
925 | G4double previous = (*PAC)[i]; |
---|
926 | PAC->push_back(previous+(*area)[i]); |
---|
927 | } |
---|
928 | (*PAC)[PAC->size()-1] = 1.; |
---|
929 | |
---|
930 | /******************************************************************************* |
---|
931 | Pre-calculated limits for the initial binary search for subsequent sampling |
---|
932 | ********************************************************************************/ |
---|
933 | |
---|
934 | //G4DataVector* ITTL = new G4DataVector(); |
---|
935 | std::vector<size_t> *ITTL = new std::vector<size_t>; |
---|
936 | std::vector<size_t> *ITTU = new std::vector<size_t>; |
---|
937 | |
---|
938 | //Just create place-holders |
---|
939 | for (size_t i=0;i<np;i++) |
---|
940 | { |
---|
941 | ITTL->push_back(0); |
---|
942 | ITTU->push_back(0); |
---|
943 | } |
---|
944 | |
---|
945 | G4double bin = 1.0/(np-1); |
---|
946 | (*ITTL)[0]=0; |
---|
947 | for (size_t i=1;i<(np-1);i++) |
---|
948 | { |
---|
949 | G4double ptst = i*bin; |
---|
950 | G4bool found = false; |
---|
951 | for (size_t j=(*ITTL)[i-1];j<np && !found;j++) |
---|
952 | { |
---|
953 | if ((*PAC)[j] > ptst) |
---|
954 | { |
---|
955 | (*ITTL)[i] = j-1; |
---|
956 | (*ITTU)[i-1] = j; |
---|
957 | found = true; |
---|
958 | } |
---|
959 | } |
---|
960 | } |
---|
961 | (*ITTU)[ITTU->size()-2] = ITTU->size()-1; |
---|
962 | (*ITTU)[ITTU->size()-1] = ITTU->size()-1; |
---|
963 | (*ITTL)[ITTL->size()-1] = ITTU->size()-2; |
---|
964 | |
---|
965 | if (ITTU->size() != np || ITTU->size() != np) |
---|
966 | { |
---|
967 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
968 | G4cout << "Problem in building the Limit Tables for Sampling: array dimensions do not match " << G4endl; |
---|
969 | G4Exception(); |
---|
970 | } |
---|
971 | |
---|
972 | |
---|
973 | /******************************************************************************** |
---|
974 | Copy tables |
---|
975 | ********************************************************************************/ |
---|
976 | G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np); |
---|
977 | for (size_t i=0;i<np;i++) |
---|
978 | { |
---|
979 | theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]); |
---|
980 | } |
---|
981 | |
---|
982 | if (verboseLevel > 2) |
---|
983 | { |
---|
984 | G4cout << "*************************************************************************" << G4endl; |
---|
985 | G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl; |
---|
986 | theTable->DumpTable(); |
---|
987 | } |
---|
988 | samplingTable->insert(std::make_pair(mat,theTable)); |
---|
989 | |
---|
990 | |
---|
991 | //Clean up temporary vectors |
---|
992 | if (x) delete x; |
---|
993 | if (a) delete a; |
---|
994 | if (b) delete b; |
---|
995 | if (c) delete c; |
---|
996 | if (err) delete err; |
---|
997 | if (area) delete area; |
---|
998 | if (PAC) delete PAC; |
---|
999 | if (ITTL) delete ITTL; |
---|
1000 | if (ITTU) delete ITTU; |
---|
1001 | |
---|
1002 | //DONE! |
---|
1003 | return; |
---|
1004 | |
---|
1005 | } |
---|
1006 | |
---|
1007 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
1008 | |
---|
1009 | void G4Penelope08RayleighModel::GetPMaxTable(const G4Material* mat) |
---|
1010 | { |
---|
1011 | if (!pMaxTable) |
---|
1012 | { |
---|
1013 | G4cout << "G4Penelope08RayleighModel::BuildPMaxTable" << G4endl; |
---|
1014 | G4cout << "Going to instanziate the pMaxTable !" << G4endl; |
---|
1015 | G4cout << "That should _not_ be here! " << G4endl; |
---|
1016 | pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
1017 | } |
---|
1018 | //check if the table is already there |
---|
1019 | if (pMaxTable->count(mat)) |
---|
1020 | return; |
---|
1021 | |
---|
1022 | //otherwise build it |
---|
1023 | if (!samplingTable) |
---|
1024 | { |
---|
1025 | G4cout << "G4Penelope08RayleighModel::BuildPMaxTable" << G4endl; |
---|
1026 | G4cout << "WARNING: samplingTable is not properly instantiated" << G4endl; |
---|
1027 | G4Exception(); |
---|
1028 | } |
---|
1029 | |
---|
1030 | if (!samplingTable->count(mat)) |
---|
1031 | InitializeSamplingAlgorithm(mat); |
---|
1032 | |
---|
1033 | G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second; |
---|
1034 | size_t tablePoints = theTable->GetNumberOfStoredPoints(); |
---|
1035 | |
---|
1036 | size_t nOfEnergyPoints = logEnergyGridPMax.size(); |
---|
1037 | G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints); |
---|
1038 | |
---|
1039 | const size_t nip = 51; //hard-coded in Penelope |
---|
1040 | |
---|
1041 | for (size_t ie=0;ie<logEnergyGridPMax.size();ie++) |
---|
1042 | { |
---|
1043 | G4double energy = std::exp(logEnergyGridPMax[ie]); |
---|
1044 | G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now |
---|
1045 | G4double Qm2 = Qm*Qm; |
---|
1046 | G4double firstQ2 = theTable->GetX(0); |
---|
1047 | G4double lastQ2 = theTable->GetX(tablePoints-1); |
---|
1048 | G4double thePMax = 0; |
---|
1049 | |
---|
1050 | if (Qm2 > firstQ2) |
---|
1051 | { |
---|
1052 | if (Qm2 < lastQ2) |
---|
1053 | { |
---|
1054 | //bisection to look for the index of Qm |
---|
1055 | size_t lowerBound = 0; |
---|
1056 | size_t upperBound = tablePoints-1; |
---|
1057 | while (lowerBound <= upperBound) |
---|
1058 | { |
---|
1059 | size_t midBin = (lowerBound + upperBound)/2; |
---|
1060 | if( Qm2 < theTable->GetX(midBin)) |
---|
1061 | { upperBound = midBin-1; } |
---|
1062 | else |
---|
1063 | { lowerBound = midBin+1; } |
---|
1064 | } |
---|
1065 | //upperBound is the output (but also lowerBounf --> should be the same!) |
---|
1066 | G4double Q1 = theTable->GetX(upperBound); |
---|
1067 | G4double Q2 = Qm2; |
---|
1068 | G4double DQ = (Q2-Q1)/((G4double)(nip-1)); |
---|
1069 | G4double theA = theTable->GetA(upperBound); |
---|
1070 | G4double theB = theTable->GetB(upperBound); |
---|
1071 | G4double thePAC = theTable->GetPAC(upperBound); |
---|
1072 | G4DataVector* fun = new G4DataVector(); |
---|
1073 | for (size_t k=0;k<nip;k++) |
---|
1074 | { |
---|
1075 | G4double qi = Q1 + k*DQ; |
---|
1076 | G4double tau = (qi-Q1)/ |
---|
1077 | (theTable->GetX(upperBound+1)-Q1); |
---|
1078 | G4double con1 = 2.0*theB*tau; |
---|
1079 | G4double ci = 1.0+theA+theB; |
---|
1080 | G4double con2 = ci-theA*tau; |
---|
1081 | G4double etap = 0; |
---|
1082 | if (std::fabs(con1) > 1.0e-16*std::fabs(con2)) |
---|
1083 | etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1; |
---|
1084 | else |
---|
1085 | etap = tau/con2; |
---|
1086 | G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)* |
---|
1087 | (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/ |
---|
1088 | ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1)); |
---|
1089 | fun->push_back(theFun); |
---|
1090 | } |
---|
1091 | //Now intergrate numerically the fun Cavalieri-Simpson's method |
---|
1092 | G4DataVector* sum = new G4DataVector; |
---|
1093 | G4double CONS = DQ*(1./12.); |
---|
1094 | G4double HCONS = 0.5*CONS; |
---|
1095 | sum->push_back(0.); |
---|
1096 | G4double secondPoint = (*sum)[0] + |
---|
1097 | (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS; |
---|
1098 | sum->push_back(secondPoint); |
---|
1099 | for (size_t hh=2;hh<nip-1;hh++) |
---|
1100 | { |
---|
1101 | G4double previous = (*sum)[hh-1]; |
---|
1102 | G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])- |
---|
1103 | (*fun)[hh+1]-(*fun)[hh-2])*HCONS; |
---|
1104 | sum->push_back(next); |
---|
1105 | } |
---|
1106 | G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]- |
---|
1107 | (*fun)[nip-3])*CONS; |
---|
1108 | sum->push_back(last); |
---|
1109 | thePMax = thePAC + (*sum)[sum->size()-1]; //last point |
---|
1110 | if (fun) delete fun; |
---|
1111 | if (sum) delete sum; |
---|
1112 | } |
---|
1113 | else |
---|
1114 | { |
---|
1115 | thePMax = 1.0; |
---|
1116 | } |
---|
1117 | } |
---|
1118 | else |
---|
1119 | { |
---|
1120 | thePMax = theTable->GetPAC(0); |
---|
1121 | } |
---|
1122 | |
---|
1123 | //Write number in the table |
---|
1124 | theVec->PutValue(ie,energy,thePMax); |
---|
1125 | } |
---|
1126 | |
---|
1127 | pMaxTable->insert(std::make_pair(mat,theVec)); |
---|
1128 | return; |
---|
1129 | |
---|
1130 | } |
---|
1131 | |
---|
1132 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
1133 | |
---|
1134 | void G4Penelope08RayleighModel::DumpFormFactorTable(const G4Material* mat) |
---|
1135 | { |
---|
1136 | G4cout << "*****************************************************************" << G4endl; |
---|
1137 | G4cout << "G4Penelope08RayleighModel: Form Factor Table for " << mat->GetName() << G4endl; |
---|
1138 | //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F |
---|
1139 | G4cout << "Q/(m_e*c) F(Q) " << G4endl; |
---|
1140 | G4cout << "*****************************************************************" << G4endl; |
---|
1141 | if (!logFormFactorTable->count(mat)) |
---|
1142 | BuildFormFactorTable(mat); |
---|
1143 | |
---|
1144 | G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second; |
---|
1145 | for (size_t i=0;i<theVec->GetVectorLength();i++) |
---|
1146 | { |
---|
1147 | G4double logQ2 = theVec->GetLowEdgeEnergy(i); |
---|
1148 | G4double Q = std::exp(0.5*logQ2); |
---|
1149 | G4double logF2 = (*theVec)[i]; |
---|
1150 | G4double F = std::exp(0.5*logF2); |
---|
1151 | G4cout << Q << " " << F << G4endl; |
---|
1152 | } |
---|
1153 | //DONE |
---|
1154 | return; |
---|
1155 | } |
---|