Geant4  10.02.p02
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 97613 2016-06-06 12:24:51Z 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 {
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 = G4Exp(logXS);
339 
340  if (verboseLevel > 2)
341  {
342  G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z <<
343  " = " << cross/barn << " barn" << G4endl;
344  }
345 
346  return cross;
347 }
348 
349 
350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 {
353 
354  /*
355  1) get composition and equivalent molecular density
356  */
357 
358  G4int nElements = material->GetNumberOfElements();
359  const G4ElementVector* elementVector = material->GetElementVector();
360  const G4double* fractionVector = material->GetFractionVector();
361 
362  std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
363  for (G4int i=0;i<nElements;i++)
364  {
365  G4double fraction = fractionVector[i];
366  G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
367  StechiometricFactors->push_back(fraction/atomicWeigth);
368  }
369  //Find max
370  G4double MaxStechiometricFactor = 0.;
371  for (G4int i=0;i<nElements;i++)
372  {
373  if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
374  MaxStechiometricFactor = (*StechiometricFactors)[i];
375  }
376  if (MaxStechiometricFactor<1e-16)
377  {
379  ed << "Inconsistent data of atomic composition for " <<
380  material->GetName() << G4endl;
381  G4Exception("G4PenelopeRayleighModel::BuildFormFactorTable()",
382  "em2042",FatalException,ed);
383  }
384  //Normalize
385  for (G4int i=0;i<nElements;i++)
386  (*StechiometricFactors)[i] /= MaxStechiometricFactor;
387 
388  // Equivalent atoms per molecule
389  G4double atomsPerMolecule = 0;
390  for (G4int i=0;i<nElements;i++)
391  atomsPerMolecule += (*StechiometricFactors)[i];
392 
393  /*
394  CREATE THE FORM FACTOR TABLE
395  */
397  theFFVec->SetSpline(true);
398 
399  for (size_t k=0;k<logQSquareGrid.size();k++)
400  {
401  G4double ff2 = 0; //squared form factor
402  for (G4int i=0;i<nElements;i++)
403  {
404  G4int iZ = (G4int) (*elementVector)[i]->GetZ();
405  G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second;
406  G4double f = (*theAtomVec)[k]; //the q-grid is always the same
407  ff2 += f*f*(*StechiometricFactors)[i];
408  }
409  if (ff2)
410  theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2)
411  }
412  logFormFactorTable->insert(std::make_pair(material,theFFVec));
413 
414  delete StechiometricFactors;
415 
416  return;
417 }
418 
419 
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
421 
422 void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
423  const G4MaterialCutsCouple* couple,
424  const G4DynamicParticle* aDynamicGamma,
425  G4double,
426  G4double)
427 {
428  // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
429  // from the Penelope2008 model. The scattering angle is sampled from the atomic
430  // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
431  // anomalous scattering effects. The Form Factor F(Q) function which appears in the
432  // analytical cross section is retrieved via the method GetFSquared(); atomic data
433  // are tabulated for F(Q). Form factor for compounds is calculated according to
434  // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
435  // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
436  // for each material and managed by G4PenelopeSamplingData objects.
437  // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
438  // increases with energy. For E=100 keV the efficiency is 100% and 86% for
439  // hydrogen and uranium, respectively.
440 
441  if (verboseLevel > 3)
442  G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
443 
444  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
445 
446  if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
447  {
451  return ;
452  }
453 
454  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
455 
456  const G4Material* theMat = couple->GetMaterial();
457 
458 
459  //1) Verify if tables are ready
460  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
461  //not invoked
464  {
465  //create a **thread-local** version of the table. Used only for G4EmCalculator and
466  //Unit Tests
467  fLocalTable = true;
469  logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
470  if (!atomicFormFactor)
471  atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
472  if (!logFormFactorTable)
473  logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
474  if (!pMaxTable)
475  pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
476  if (!samplingTable)
477  samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
478  }
479 
480  if (!samplingTable->count(theMat))
481  {
482  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
483  //not filled up. This can happen in a UnitTest
484  if (verboseLevel > 0)
485  {
486  //Issue a G4Exception (warning) only in verbose mode
488  ed << "Unable to find the samplingTable data for " <<
489  theMat->GetName() << G4endl;
490  ed << "This can happen only in Unit Tests" << G4endl;
491  G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
492  "em2019",JustWarning,ed);
493  }
494  const G4ElementVector* theElementVector = theMat->GetElementVector();
495  //protect file reading via autolock
496  G4AutoLock lock(&PenelopeRayleighModelMutex);
497  for (size_t j=0;j<theMat->GetNumberOfElements();j++)
498  {
499  G4int iZ = (G4int) theElementVector->at(j)->GetZ();
500  if (!logAtomicCrossSection->count(iZ))
501  {
502  lock.lock();
503  ReadDataFile(iZ);
504  lock.unlock();
505  }
506  }
507  lock.lock();
508  //1) If the table has not been built for the material, do it!
509  if (!logFormFactorTable->count(theMat))
510  BuildFormFactorTable(theMat);
511 
512  //2) retrieve or build the sampling table
513  if (!(samplingTable->count(theMat)))
515 
516  //3) retrieve or build the pMax data
517  if (!pMaxTable->count(theMat))
518  GetPMaxTable(theMat);
519  lock.unlock();
520  }
521 
522  //Ok, restart the job
523 
524  G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
525  G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;
526 
527  G4double cosTheta = 1.0;
528 
529  //OK, ready to go!
530  G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
531 
532  if (qmax < 1e-10) //very low momentum transfer
533  {
534  G4bool loopAgain=false;
535  do
536  {
537  loopAgain = false;
538  cosTheta = 1.0-2.0*G4UniformRand();
539  G4double G = 0.5*(1+cosTheta*cosTheta);
540  if (G4UniformRand()>G)
541  loopAgain = true;
542  }while(loopAgain);
543  }
544  else //larger momentum transfer
545  {
546  size_t nData = theDataTable->GetNumberOfStoredPoints();
547  G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
548  G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
549 
550  G4bool loopAgain = false;
551  G4double MaxPValue = thePMax->Value(photonEnergy0);
552  G4double xx=0;
553 
554  //Sampling by rejection method. The rejection function is
555  //G = 0.5*(1+cos^2(theta))
556  //
557  do{
558  loopAgain = false;
559  G4double RandomMax = G4UniformRand()*MaxPValue;
560  xx = theDataTable->SampleValue(RandomMax);
561  //xx is a random value of q^2 in (0,q2max),sampled according to
562  //F(Q^2) via the RITA algorithm
563  if (xx > q2max)
564  loopAgain = true;
565  cosTheta = 1.0-2.0*xx/q2max;
566  G4double G = 0.5*(1+cosTheta*cosTheta);
567  if (G4UniformRand()>G)
568  loopAgain = true;
569  }while(loopAgain);
570  }
571 
572  G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
573 
574  // Scattered photon angles. ( Z - axis along the parent photon)
575  G4double phi = twopi * G4UniformRand() ;
576  G4double dirX = sinTheta*std::cos(phi);
577  G4double dirY = sinTheta*std::sin(phi);
578  G4double dirZ = cosTheta;
579 
580  // Update G4VParticleChange for the scattered photon
581  G4ThreeVector photonDirection1(dirX, dirY, dirZ);
582 
583  photonDirection1.rotateUz(photonDirection0);
584  fParticleChange->ProposeMomentumDirection(photonDirection1) ;
585  fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
586 
587  return;
588 }
589 
590 
591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
592 
594 {
595 
596  if (verboseLevel > 2)
597  {
598  G4cout << "G4PenelopeRayleighModel::ReadDataFile()" << G4endl;
599  G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
600  }
601 
602  char* path = getenv("G4LEDATA");
603  if (!path)
604  {
605  G4String excep = "G4LEDATA environment variable not set!";
606  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
607  "em0006",FatalException,excep);
608  return;
609  }
610 
611  /*
612  Read first the cross section file
613  */
614  std::ostringstream ost;
615  if (Z>9)
616  ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
617  else
618  ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
619  std::ifstream file(ost.str().c_str());
620  if (!file.is_open())
621  {
622  G4String excep = "Data file " + G4String(ost.str()) + " not found!";
623  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
624  "em0003",FatalException,excep);
625  }
626  G4int readZ =0;
627  size_t nPoints= 0;
628  file >> readZ >> nPoints;
629  //check the right file is opened.
630  if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
631  {
633  ed << "Corrupted data file for Z=" << Z << G4endl;
634  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
635  "em0005",FatalException,ed);
636  return;
637  }
638  G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints);
639  G4double ene=0,f1=0,f2=0,xs=0;
640  for (size_t i=0;i<nPoints;i++)
641  {
642  file >> ene >> f1 >> f2 >> xs;
643  //dimensional quantities
644  ene *= eV;
645  xs *= cm2;
646  theVec->PutValue(i,std::log(ene),std::log(xs));
647  if (file.eof() && i != (nPoints-1)) //file ended too early
648  {
650  ed << "Corrupted data file for Z=" << Z << G4endl;
651  ed << "Found less than " << nPoints << "entries " <<G4endl;
652  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
653  "em0005",FatalException,ed);
654  }
655  }
657  {
658  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
659  "em2044",FatalException,"Unable to allocate the atomic cross section table");
660  delete theVec;
661  return;
662  }
663  logAtomicCrossSection->insert(std::make_pair(Z,theVec));
664  file.close();
665 
666  /*
667  Then read the form factor file
668  */
669  std::ostringstream ost2;
670  if (Z>9)
671  ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
672  else
673  ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
674  file.open(ost2.str().c_str());
675  if (!file.is_open())
676  {
677  G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
678  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
679  "em0003",FatalException,excep);
680  }
681  file >> readZ >> nPoints;
682  //check the right file is opened.
683  if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
684  {
686  ed << "Corrupted data file for Z=" << Z << G4endl;
687  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
688  "em0005",FatalException,ed);
689  return;
690  }
691  G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
692  G4double q=0,ff=0,incoh=0;
693  G4bool fillQGrid = false;
694  //fill this vector only the first time.
695  if (!logQSquareGrid.size())
696  fillQGrid = true;
697  for (size_t i=0;i<nPoints;i++)
698  {
699  file >> q >> ff >> incoh;
700  //q and ff are dimensionless (q is in units of (m_e*c)
701  theFFVec->PutValue(i,q,ff);
702  if (fillQGrid)
703  {
704  logQSquareGrid.push_back(2.0*std::log(q));
705  }
706  if (file.eof() && i != (nPoints-1)) //file ended too early
707  {
709  ed << "Corrupted data file for Z=" << Z << G4endl;
710  ed << "Found less than " << nPoints << "entries " <<G4endl;
711  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
712  "em0005",FatalException,ed);
713  }
714  }
715  if (!atomicFormFactor)
716  {
717  G4Exception("G4PenelopeRayleighModel::ReadDataFile()",
718  "em2045",FatalException,
719  "Unable to allocate the atomicFormFactor data table");
720  delete theFFVec;
721  return;
722  }
723  atomicFormFactor->insert(std::make_pair(Z,theFFVec));
724  file.close();
725  return;
726 }
727 
728 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
729 
731 {
732  G4double f2 = 0;
733  //Input value QSquared could be zero: protect the log() below against
734  //the FPE exception
735  //If Q<1e-10, set Q to 1e-10
736  G4double logQSquared = (QSquared>1e-10) ? std::log(QSquared) : -23.;
737  //last value of the table
738  G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1];
739 
740 
741  //now it should be all right
742  G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
743 
744  if (!theVec)
745  {
747  ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
748  G4Exception("G4PenelopeRayleighModel::GetFSquared()",
749  "em2046",FatalException,ed);
750  return 0;
751  }
752  if (logQSquared < -20) // Q < 1e-9
753  {
754  G4double logf2 = (*theVec)[0]; //first value of the table
755  f2 = G4Exp(logf2);
756  }
757  else if (logQSquared > maxlogQ2)
758  f2 =0;
759  else
760  {
761  //log(Q^2) vs. log(F^2)
762  G4double logf2 = theVec->Value(logQSquared);
763  f2 = G4Exp(logf2);
764 
765  }
766  if (verboseLevel > 3)
767  {
768  G4cout << "G4PenelopeRayleighModel::GetFSquared() in " << mat->GetName() << G4endl;
769  G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl;
770  }
771  return f2;
772 }
773 
774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
775 
777 {
778 
779  G4double q2min = 0;
780  G4double q2max = 0;
781  const size_t np = 150; //hard-coded in Penelope
782  //G4cout << "Init N= " << logQSquareGrid.size() << G4endl;
783  for (size_t i=1;i<logQSquareGrid.size();i++)
784  {
785  G4double Q2 = G4Exp(logQSquareGrid[i]);
786  if (GetFSquared(mat,Q2) > 1e-35)
787  {
788  q2max = G4Exp(logQSquareGrid[i-1]);
789  }
790  //G4cout << "Q2= " << Q2 << " q2max= " << q2max << G4endl;
791  }
792 
793  size_t nReducedPoints = np/4;
794 
795  //check for errors
796  if (np < 16)
797  {
798  G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
799  "em2047",FatalException,
800  "Too few points to initialize the sampling algorithm");
801  }
802  if (q2min > (q2max-1e-10))
803  {
804  G4cout << "q2min= " << q2min << " q2max= " << q2max << G4endl;
805  G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
806  "em2048",FatalException,
807  "Too narrow grid to initialize the sampling algorithm");
808  }
809 
810  //This is subroutine RITAI0 of Penelope
811  //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material*
812 
813  //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
814  G4DataVector* x = new G4DataVector();
815 
816  /*******************************************************************************
817  Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
818  ********************************************************************************/
819  size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2);
820  const G4int nip = 51; //hard-coded in Penelope
821 
822  G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
823  x->push_back(q2min);
824  for (size_t i=1;i<NUNIF-1;i++)
825  {
826  G4double app = q2min + i*dx;
827  x->push_back(app); //increase
828  }
829  x->push_back(q2max);
830 
831  if (verboseLevel> 3)
832  G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
833 
834  //Allocate temporary storage vectors
835  G4DataVector* area = new G4DataVector();
836  G4DataVector* a = new G4DataVector();
837  G4DataVector* b = new G4DataVector();
838  G4DataVector* c = new G4DataVector();
839  G4DataVector* err = new G4DataVector();
840 
841  for (size_t i=0;i<NUNIF-1;i++) //build all points but the last
842  {
843  //Temporary vectors for this loop
844  G4DataVector* pdfi = new G4DataVector();
845  G4DataVector* pdfih = new G4DataVector();
846  G4DataVector* sumi = new G4DataVector();
847 
848  G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
849  G4double pdfmax = 0;
850  for (G4int k=0;k<nip;k++)
851  {
852  G4double xik = (*x)[i]+k*dxi;
853  G4double pdfk = std::max(GetFSquared(mat,xik),0.);
854  pdfi->push_back(pdfk);
855  pdfmax = std::max(pdfmax,pdfk);
856  if (k < (nip-1))
857  {
858  G4double xih = xik + 0.5*dxi;
859  G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
860  pdfih->push_back(pdfIK);
861  pdfmax = std::max(pdfmax,pdfIK);
862  }
863  }
864 
865  //Simpson's integration
866  G4double cons = dxi*0.5*(1./3.);
867  sumi->push_back(0.);
868  for (G4int k=1;k<nip;k++)
869  {
870  G4double previous = (*sumi)[k-1];
871  G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
872  sumi->push_back(next);
873  }
874 
875  G4double lastIntegral = (*sumi)[sumi->size()-1];
876  area->push_back(lastIntegral);
877  //Normalize cumulative function
878  G4double factor = 1.0/lastIntegral;
879  for (size_t k=0;k<sumi->size();k++)
880  (*sumi)[k] *= factor;
881 
882  //When the PDF vanishes at one of the interval end points, its value is modified
883  if ((*pdfi)[0] < 1e-35)
884  (*pdfi)[0] = 1e-5*pdfmax;
885  if ((*pdfi)[pdfi->size()-1] < 1e-35)
886  (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
887 
888  G4double pli = (*pdfi)[0]*factor;
889  G4double pui = (*pdfi)[pdfi->size()-1]*factor;
890  G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
891  G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
892  G4double C_temp = 1.0+A_temp+B_temp;
893  if (C_temp < 1e-35)
894  {
895  a->push_back(0.);
896  b->push_back(0.);
897  c->push_back(1.);
898  }
899  else
900  {
901  a->push_back(A_temp);
902  b->push_back(B_temp);
903  c->push_back(C_temp);
904  }
905 
906  //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
907  //and the true pdf, extended over the interval (X(I),X(I+1))
908  G4int icase = 1; //loop code
909  G4bool reLoop = false;
910  err->push_back(0.);
911  do
912  {
913  reLoop = false;
914  (*err)[i] = 0.; //zero variable
915  for (G4int k=0;k<nip;k++)
916  {
917  G4double rr = (*sumi)[k];
918  G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
919  ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
920  if (k == 0 || k == nip-1)
921  (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
922  else
923  (*err)[i] += std::fabs(pap-(*pdfi)[k]);
924  }
925  (*err)[i] *= dxi;
926 
927  //If err(I) is too large, the pdf is approximated by a uniform distribution
928  if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
929  {
930  (*b)[i] = 0;
931  (*a)[i] = 0;
932  (*c)[i] = 1.;
933  icase = 2;
934  reLoop = true;
935  }
936  }while(reLoop);
937 
938  delete pdfi;
939  delete pdfih;
940  delete sumi;
941  } //end of first loop over i
942 
943  //Now assign last point
944  (*x)[x->size()-1] = q2max;
945  a->push_back(0.);
946  b->push_back(0.);
947  c->push_back(0.);
948  err->push_back(0.);
949  area->push_back(0.);
950 
951  if (x->size() != NUNIF || a->size() != NUNIF ||
952  err->size() != NUNIF || area->size() != NUNIF)
953  {
955  ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
956  G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
957  "em2049",FatalException,ed);
958  }
959 
960  /*******************************************************************************
961  New grid points are added by halving the sub-intervals with the largest absolute error
962  This is done up to np=150 points in the grid
963  ********************************************************************************/
964  do
965  {
966  G4double maxError = 0.0;
967  size_t iErrMax = 0;
968  for (size_t i=0;i<err->size()-2;i++)
969  {
970  //maxError is the lagest of the interval errors err[i]
971  if ((*err)[i] > maxError)
972  {
973  maxError = (*err)[i];
974  iErrMax = i;
975  }
976  }
977 
978  //OK, now I have to insert one new point in the position iErrMax
979  G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
980 
981  x->insert(x->begin()+iErrMax+1,newx);
982  //Add place-holders in the other vectors
983  area->insert(area->begin()+iErrMax+1,0.);
984  a->insert(a->begin()+iErrMax+1,0.);
985  b->insert(b->begin()+iErrMax+1,0.);
986  c->insert(c->begin()+iErrMax+1,0.);
987  err->insert(err->begin()+iErrMax+1,0.);
988 
989  //Now calculate the other parameters
990  for (size_t i=iErrMax;i<=iErrMax+1;i++)
991  {
992  //Temporary vectors for this loop
993  G4DataVector* pdfi = new G4DataVector();
994  G4DataVector* pdfih = new G4DataVector();
995  G4DataVector* sumi = new G4DataVector();
996 
997  G4double dxLocal = (*x)[i+1]-(*x)[i];
998  G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
999  G4double pdfmax = 0;
1000  for (G4int k=0;k<nip;k++)
1001  {
1002  G4double xik = (*x)[i]+k*dxi;
1003  G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1004  pdfi->push_back(pdfk);
1005  pdfmax = std::max(pdfmax,pdfk);
1006  if (k < (nip-1))
1007  {
1008  G4double xih = xik + 0.5*dxi;
1009  G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1010  pdfih->push_back(pdfIK);
1011  pdfmax = std::max(pdfmax,pdfIK);
1012  }
1013  }
1014 
1015  //Simpson's integration
1016  G4double cons = dxi*0.5*(1./3.);
1017  sumi->push_back(0.);
1018  for (G4int k=1;k<nip;k++)
1019  {
1020  G4double previous = (*sumi)[k-1];
1021  G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1022  sumi->push_back(next);
1023  }
1024  G4double lastIntegral = (*sumi)[sumi->size()-1];
1025  (*area)[i] = lastIntegral;
1026 
1027  //Normalize cumulative function
1028  G4double factor = 1.0/lastIntegral;
1029  for (size_t k=0;k<sumi->size();k++)
1030  (*sumi)[k] *= factor;
1031 
1032  //When the PDF vanishes at one of the interval end points, its value is modified
1033  if ((*pdfi)[0] < 1e-35)
1034  (*pdfi)[0] = 1e-5*pdfmax;
1035  if ((*pdfi)[pdfi->size()-1] < 1e-35)
1036  (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1037 
1038  G4double pli = (*pdfi)[0]*factor;
1039  G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1040  G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
1041  G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
1042  G4double C_temp = 1.0+A_temp+B_temp;
1043  if (C_temp < 1e-35)
1044  {
1045  (*a)[i]= 0.;
1046  (*b)[i] = 0.;
1047  (*c)[i] = 1;
1048  }
1049  else
1050  {
1051  (*a)[i]= A_temp;
1052  (*b)[i] = B_temp;
1053  (*c)[i] = C_temp;
1054  }
1055  //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1056  //and the true pdf, extended over the interval (X(I),X(I+1))
1057  G4int icase = 1; //loop code
1058  G4bool reLoop = false;
1059  do
1060  {
1061  reLoop = false;
1062  (*err)[i] = 0.; //zero variable
1063  for (G4int k=0;k<nip;k++)
1064  {
1065  G4double rr = (*sumi)[k];
1066  G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1067  ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1068  if (k == 0 || k == nip-1)
1069  (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1070  else
1071  (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1072  }
1073  (*err)[i] *= dxi;
1074 
1075  //If err(I) is too large, the pdf is approximated by a uniform distribution
1076  if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1077  {
1078  (*b)[i] = 0;
1079  (*a)[i] = 0;
1080  (*c)[i] = 1.;
1081  icase = 2;
1082  reLoop = true;
1083  }
1084  }while(reLoop);
1085  delete pdfi;
1086  delete pdfih;
1087  delete sumi;
1088  }
1089  }while(x->size() < np);
1090 
1091  if (x->size() != np || a->size() != np ||
1092  err->size() != np || area->size() != np)
1093  {
1094  G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1095  "em2050",FatalException,
1096  "Problem in building the extended Table for Sampling: array dimensions do not match ");
1097  }
1098 
1099  /*******************************************************************************
1100  Renormalization
1101  ********************************************************************************/
1102  G4double ws = 0;
1103  for (size_t i=0;i<np-1;i++)
1104  ws += (*area)[i];
1105  ws = 1.0/ws;
1106  G4double errMax = 0;
1107  for (size_t i=0;i<np-1;i++)
1108  {
1109  (*area)[i] *= ws;
1110  (*err)[i] *= ws;
1111  errMax = std::max(errMax,(*err)[i]);
1112  }
1113 
1114  //Vector with the normalized cumulative distribution
1115  G4DataVector* PAC = new G4DataVector();
1116  PAC->push_back(0.);
1117  for (size_t i=0;i<np-1;i++)
1118  {
1119  G4double previous = (*PAC)[i];
1120  PAC->push_back(previous+(*area)[i]);
1121  }
1122  (*PAC)[PAC->size()-1] = 1.;
1123 
1124  /*******************************************************************************
1125  Pre-calculated limits for the initial binary search for subsequent sampling
1126  ********************************************************************************/
1127 
1128  //G4DataVector* ITTL = new G4DataVector();
1129  std::vector<size_t> *ITTL = new std::vector<size_t>;
1130  std::vector<size_t> *ITTU = new std::vector<size_t>;
1131 
1132  //Just create place-holders
1133  for (size_t i=0;i<np;i++)
1134  {
1135  ITTL->push_back(0);
1136  ITTU->push_back(0);
1137  }
1138 
1139  G4double bin = 1.0/(np-1);
1140  (*ITTL)[0]=0;
1141  for (size_t i=1;i<(np-1);i++)
1142  {
1143  G4double ptst = i*bin;
1144  G4bool found = false;
1145  for (size_t j=(*ITTL)[i-1];j<np && !found;j++)
1146  {
1147  if ((*PAC)[j] > ptst)
1148  {
1149  (*ITTL)[i] = j-1;
1150  (*ITTU)[i-1] = j;
1151  found = true;
1152  }
1153  }
1154  }
1155  (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1156  (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1157  (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1158 
1159  if (ITTU->size() != np || ITTU->size() != np)
1160  {
1161  G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1162  "em2051",FatalException,
1163  "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1164  }
1165 
1166 
1167  /********************************************************************************
1168  Copy tables
1169  ********************************************************************************/
1170  G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
1171  for (size_t i=0;i<np;i++)
1172  {
1173  theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1174  }
1175 
1176  if (verboseLevel > 2)
1177  {
1178  G4cout << "*************************************************************************" <<
1179  G4endl;
1180  G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
1181  theTable->DumpTable();
1182  }
1183  samplingTable->insert(std::make_pair(mat,theTable));
1184 
1185 
1186  //Clean up temporary vectors
1187  delete x;
1188  delete a;
1189  delete b;
1190  delete c;
1191  delete err;
1192  delete area;
1193  delete PAC;
1194  delete ITTL;
1195  delete ITTU;
1196 
1197  //DONE!
1198  return;
1199 
1200 }
1201 
1202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1203 
1205 {
1206 
1207  if (!pMaxTable)
1208  {
1209  G4cout << "G4PenelopeRayleighModel::BuildPMaxTable" << G4endl;
1210  G4cout << "Going to instanziate the pMaxTable !" << G4endl;
1211  G4cout << "That should _not_ be here! " << G4endl;
1212  pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
1213  }
1214  //check if the table is already there
1215  if (pMaxTable->count(mat))
1216  return;
1217 
1218  //otherwise build it
1219  if (!samplingTable)
1220  {
1221  G4Exception("G4PenelopeRayleighModel::GetPMaxTable()",
1222  "em2052",FatalException,
1223  "SamplingTable is not properly instantiated");
1224  return;
1225  }
1226 
1227  //This should not be: the sampling table is built before the p-table
1228  if (!samplingTable->count(mat))
1229  {
1231  ed << "Sampling table for material " << mat->GetName() << " not found";
1232  G4Exception("G4PenelopeRayleighModel::GetPMaxTable()",
1233  "em2052",FatalException,
1234  ed);
1235  return;
1236  }
1237 
1238  G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second;
1239  size_t tablePoints = theTable->GetNumberOfStoredPoints();
1240 
1241  size_t nOfEnergyPoints = logEnergyGridPMax.size();
1242  G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
1243 
1244  const size_t nip = 51; //hard-coded in Penelope
1245 
1246  for (size_t ie=0;ie<logEnergyGridPMax.size();ie++)
1247  {
1249  G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
1250  G4double Qm2 = Qm*Qm;
1251  G4double firstQ2 = theTable->GetX(0);
1252  G4double lastQ2 = theTable->GetX(tablePoints-1);
1253  G4double thePMax = 0;
1254 
1255  if (Qm2 > firstQ2)
1256  {
1257  if (Qm2 < lastQ2)
1258  {
1259  //bisection to look for the index of Qm
1260  size_t lowerBound = 0;
1261  size_t upperBound = tablePoints-1;
1262  while (lowerBound <= upperBound)
1263  {
1264  size_t midBin = (lowerBound + upperBound)/2;
1265  if( Qm2 < theTable->GetX(midBin))
1266  { upperBound = midBin-1; }
1267  else
1268  { lowerBound = midBin+1; }
1269  }
1270  //upperBound is the output (but also lowerBounf --> should be the same!)
1271  G4double Q1 = theTable->GetX(upperBound);
1272  G4double Q2 = Qm2;
1273  G4double DQ = (Q2-Q1)/((G4double)(nip-1));
1274  G4double theA = theTable->GetA(upperBound);
1275  G4double theB = theTable->GetB(upperBound);
1276  G4double thePAC = theTable->GetPAC(upperBound);
1277  G4DataVector* fun = new G4DataVector();
1278  for (size_t k=0;k<nip;k++)
1279  {
1280  G4double qi = Q1 + k*DQ;
1281  G4double tau = (qi-Q1)/
1282  (theTable->GetX(upperBound+1)-Q1);
1283  G4double con1 = 2.0*theB*tau;
1284  G4double ci = 1.0+theA+theB;
1285  G4double con2 = ci-theA*tau;
1286  G4double etap = 0;
1287  if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1288  etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1289  else
1290  etap = tau/con2;
1291  G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
1292  (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1293  ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
1294  fun->push_back(theFun);
1295  }
1296  //Now intergrate numerically the fun Cavalieri-Simpson's method
1297  G4DataVector* sum = new G4DataVector;
1298  G4double CONS = DQ*(1./12.);
1299  G4double HCONS = 0.5*CONS;
1300  sum->push_back(0.);
1301  G4double secondPoint = (*sum)[0] +
1302  (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1303  sum->push_back(secondPoint);
1304  for (size_t hh=2;hh<nip-1;hh++)
1305  {
1306  G4double previous = (*sum)[hh-1];
1307  G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1308  (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1309  sum->push_back(next);
1310  }
1311  G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1312  (*fun)[nip-3])*CONS;
1313  sum->push_back(last);
1314  thePMax = thePAC + (*sum)[sum->size()-1]; //last point
1315  delete fun;
1316  delete sum;
1317  }
1318  else
1319  {
1320  thePMax = 1.0;
1321  }
1322  }
1323  else
1324  {
1325  thePMax = theTable->GetPAC(0);
1326  }
1327 
1328  //Write number in the table
1329  theVec->PutValue(ie,energy,thePMax);
1330  }
1331 
1332  pMaxTable->insert(std::make_pair(mat,theVec));
1333  return;
1334 
1335 }
1336 
1337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1338 
1340 {
1341  G4cout << "*****************************************************************" << G4endl;
1342  G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
1343  //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1344  G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1345  G4cout << "*****************************************************************" << G4endl;
1346  if (!logFormFactorTable->count(mat))
1347  BuildFormFactorTable(mat);
1348 
1349  G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
1350  for (size_t i=0;i<theVec->GetVectorLength();i++)
1351  {
1352  G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1353  G4double Q = G4Exp(0.5*logQ2);
1354  G4double logF2 = (*theVec)[i];
1355  G4double F = G4Exp(0.5*logF2);
1356  G4cout << Q << " " << F << G4endl;
1357  }
1358  //DONE
1359  return;
1360 }
1361 
1362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1363 
1365 {
1366  if(!fParticle) {
1367  fParticle = p;
1368  }
1369 }
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
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
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 *)