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