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