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