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