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