source: trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08RayleighModel.cc@ 1346

Last change on this file since 1346 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 36.8 KB
Line 
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: emlowen-V09-03-54 $
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
51G4Penelope08RayleighModel::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
88G4Penelope08RayleighModel::~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....
109void 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
143void 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
184G4double 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....
232void 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
301void 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
421void 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
544G4double 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
587void 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
1009void 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
1134void 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}
Note: See TracBrowser for help on using the repository browser.