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