Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PenelopeRayleighModelMI.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// Author: Luciano Pandola and Gianfranco Paterno
27//
28// -------------------------------------------------------------------
29// History:
30// 03 Dec 2009 L. Pandola 1st implementation
31// 25 May 2011 L. Pandola Renamed (make v2008 as default Penelope)
32// 27 Sep 2013 L. Pandola Migration to MT paradigm
33// 20 Aug 2017 G. Paterno Molecular Interference implementation
34// 24 Mar 2019 G. Paterno Improved Molecular Interference implementation
35// 20 Jun 2020 G. Paterno Read qext separately and leave original atomic form factors
36// 27 Aug 2020 G. Paterno Further improvement of MI implementation
37//
38// -------------------------------------------------------------------
39// Class description:
40// Low Energy Electromagnetic Physics, Rayleigh Scattering
41// with the model from Penelope, version 2008
42// -------------------------------------------------------------------
43//
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46#include "G4PenelopeRayleighModelMI.hh"
47
48#include "G4PenelopeSamplingData.hh"
49#include "G4ParticleDefinition.hh"
50#include "G4MaterialCutsCouple.hh"
51#include "G4ProductionCutsTable.hh"
52#include "G4DynamicParticle.hh"
53#include "G4ElementTable.hh"
54#include "G4Element.hh"
55#include "G4PhysicsFreeVector.hh"
56#include "G4AutoLock.hh"
57#include "G4Exp.hh"
58#include "G4ExtendedMaterial.hh"
59#include "G4CrystalExtension.hh"
60#include "G4EmParameters.hh"
61
62#include "G4PhysicalConstants.hh"
63#include "G4SystemOfUnits.hh"
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66const G4int G4PenelopeRayleighModelMI::fMaxZ;
67G4PhysicsFreeVector* G4PenelopeRayleighModelMI::fLogAtomicCrossSection[] = {nullptr};
68G4PhysicsFreeVector* G4PenelopeRayleighModelMI::fAtomicFormFactor[] = {nullptr};
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4PenelopeRayleighModelMI::G4PenelopeRayleighModelMI(const G4ParticleDefinition* part,
73 const G4String& nam) :
74 G4VEmModel(nam),
75 fParticleChange(nullptr),fParticle(nullptr),fLogFormFactorTable(nullptr),fPMaxTable(nullptr),
76 fSamplingTable(nullptr),fMolInterferenceData(nullptr),fAngularFunction(nullptr), fKnownMaterials(nullptr),
77 fIsInitialised(false),fLocalTable(false),fIsMIActive(true)
78{
79 fIntrinsicLowEnergyLimit = 100.0*eV;
80 fIntrinsicHighEnergyLimit = 100.0*GeV;
81 //SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
82 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
83
84 if (part) SetParticle(part);
85
86 fVerboseLevel = 0;
87 // Verbosity scale:
88 // 0 = nothing
89 // 1 = warning for energy non-conservation
90 // 2 = details of energy budget
91 // 3 = calculation of FF and CS, file openings, sampling of atoms
92 // 4 = entering in methods
93
94 //build the energy grid. It is the same for all materials
95 G4double logenergy = G4Log(fIntrinsicLowEnergyLimit/2.);
96 G4double logmaxenergy = G4Log(1.5*fIntrinsicHighEnergyLimit);
97 //finer grid below 160 keV
98 G4double logtransitionenergy = G4Log(160*keV);
99 G4double logfactor1 = G4Log(10.)/250.;
100 G4double logfactor2 = logfactor1*10;
101 fLogEnergyGridPMax.push_back(logenergy);
102 do {
103 if (logenergy < logtransitionenergy)
104 logenergy += logfactor1;
105 else
106 logenergy += logfactor2;
107 fLogEnergyGridPMax.push_back(logenergy);
108 } while (logenergy < logmaxenergy);
109}
110
111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112
113G4PenelopeRayleighModelMI::~G4PenelopeRayleighModelMI()
114{
115 if (IsMaster() || fLocalTable) {
116
117 for(G4int i=0; i<=fMaxZ; ++i)
118 {
119 if(fLogAtomicCrossSection[i])
120 {
121 delete fLogAtomicCrossSection[i];
122 fLogAtomicCrossSection[i] = nullptr;
123 }
124 if(fAtomicFormFactor[i])
125 {
126 delete fAtomicFormFactor[i];
127 fAtomicFormFactor[i] = nullptr;
128 }
129 }
130 if (fMolInterferenceData) {
131 for (auto& item : (*fMolInterferenceData))
132 if (item.second) delete item.second;
133 delete fMolInterferenceData;
134 fMolInterferenceData = nullptr;
135 }
136 if (fKnownMaterials)
137 {
138 delete fKnownMaterials;
139 fKnownMaterials = nullptr;
140 }
141 if (fAngularFunction)
142 {
143 delete fAngularFunction;
144 fAngularFunction = nullptr;
145 }
146 ClearTables();
147 }
148}
149
150//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151
152void G4PenelopeRayleighModelMI::ClearTables()
153{
154 if (fLogFormFactorTable) {
155 for (auto& item : (*fLogFormFactorTable))
156 if (item.second) delete item.second;
157 delete fLogFormFactorTable;
158 fLogFormFactorTable = nullptr; //zero explicitly
159 }
160
161 if (fPMaxTable) {
162 for (auto& item : (*fPMaxTable))
163 if (item.second) delete item.second;
164 delete fPMaxTable;
165 fPMaxTable = nullptr; //zero explicitly
166 }
167
168 if (fSamplingTable) {
169 for (auto& item : (*fSamplingTable))
170 if (item.second) delete item.second;
171 delete fSamplingTable;
172 fSamplingTable = nullptr; //zero explicitly
173 }
174
175 return;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4PenelopeRayleighModelMI::Initialise(const G4ParticleDefinition* part,
181 const G4DataVector& )
182{
183 if (fVerboseLevel > 3)
184 G4cout << "Calling G4PenelopeRayleighModelMI::Initialise()" << G4endl;
185
186 SetParticle(part);
187
188 if (fVerboseLevel)
189 G4cout << "# Molecular Interference is " << (fIsMIActive ? "ON" : "OFF") << " #" << G4endl;
190
191 //Only the master model creates/fills/destroys the tables
192 if (IsMaster() && part == fParticle) {
193 //clear tables depending on materials, not the atomic ones
194 ClearTables();
195
196 //Use here the highest verbosity, from G4EmParameter or local
197 G4int globVerb = G4EmParameters::Instance()->Verbose();
198 if (globVerb > fVerboseLevel)
199 {
200 fVerboseLevel = globVerb;
201 if (fVerboseLevel)
202 G4cout << "Verbosity level of G4PenelopeRayleighModelMI set to " << fVerboseLevel <<
203 " from G4EmParameters()" << G4endl;
204 }
205 if (fVerboseLevel > 3)
206 G4cout << "Calling G4PenelopeRayleighModelMI::Initialise() [master]" << G4endl;
207
208 //Load the list of known materials and the DCS integration grid
209 if (fIsMIActive)
210 {
211 if (!fKnownMaterials)
212 fKnownMaterials = new std::map<G4String,G4String>;
213 if (!fKnownMaterials->size())
214 LoadKnownMIFFMaterials();
215 if (!fAngularFunction)
216 {
217 //Create and fill once
218 fAngularFunction = new G4PhysicsFreeVector(fNtheta);
219 CalculateThetaAndAngFun(); //angular funtion for DCS integration
220 }
221 }
222 if (fIsMIActive && !fMolInterferenceData)
223 fMolInterferenceData = new std::map<G4String,G4PhysicsFreeVector*>;
224 if (!fLogFormFactorTable)
225 fLogFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
226 if (!fPMaxTable)
227 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
228 if (!fSamplingTable)
229 fSamplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
230
231 //loop on the used materials
232 G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
233
234 for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i) {
235 const G4Material* material =
236 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
237 const G4ElementVector* theElementVector = material->GetElementVector();
238
239 for (std::size_t j=0;j<material->GetNumberOfElements();++j) {
240 G4int iZ = theElementVector->at(j)->GetZasInt();
241 //read data files only in the master
242 if (!fLogAtomicCrossSection[iZ])
243 ReadDataFile(iZ);
244 }
245
246 //1) Read MI form factors
247 if (fIsMIActive && !fMolInterferenceData->count(material->GetName()))
248 ReadMolInterferenceData(material->GetName());
249
250 //2) If the table has not been built for the material, do it!
251 if (!fLogFormFactorTable->count(material))
252 BuildFormFactorTable(material);
253
254 //3) retrieve or build the sampling table
255 if (!(fSamplingTable->count(material)))
256 InitializeSamplingAlgorithm(material);
257
258 //4) retrieve or build the pMax data
259 if (!fPMaxTable->count(material))
260 GetPMaxTable(material);
261 }
262
263 if (fVerboseLevel > 1) {
264 G4cout << G4endl << "Penelope Rayleigh model v2008 is initialized" << G4endl
265 << "Energy range: "
266 << LowEnergyLimit() / keV << " keV - "
267 << HighEnergyLimit() / GeV << " GeV"
268 << G4endl;
269 }
270 }
271
272 if (fIsInitialised)
273 return;
274 fParticleChange = GetParticleChangeForGamma();
275 fIsInitialised = true;
276}
277
278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279
280void G4PenelopeRayleighModelMI::InitialiseLocal(const G4ParticleDefinition* part,
281 G4VEmModel *masterModel)
282{
283 if (fVerboseLevel > 3)
284 G4cout << "Calling G4PenelopeRayleighModelMI::InitialiseLocal()" << G4endl;
285
286 //Check that particle matches: one might have multiple master models
287 //(e.g. for e+ and e-)
288 if (part == fParticle) {
289
290 //Get the const table pointers from the master to the workers
291 const G4PenelopeRayleighModelMI* theModel =
292 static_cast<G4PenelopeRayleighModelMI*> (masterModel);
293
294 //Copy pointers to the data tables
295 for(G4int i=0; i<=fMaxZ; ++i)
296 {
297 fLogAtomicCrossSection[i] = theModel->fLogAtomicCrossSection[i];
298 fAtomicFormFactor[i] = theModel->fAtomicFormFactor[i];
299 }
300 fMolInterferenceData = theModel->fMolInterferenceData;
301 fLogFormFactorTable = theModel->fLogFormFactorTable;
302 fPMaxTable = theModel->fPMaxTable;
303 fSamplingTable = theModel->fSamplingTable;
304 fKnownMaterials = theModel->fKnownMaterials;
305 fAngularFunction = theModel->fAngularFunction;
306
307 //Copy the G4DataVector with the grid
308 fLogQSquareGrid = theModel->fLogQSquareGrid;
309
310 //Same verbosity for all workers, as the master
311 fVerboseLevel = theModel->fVerboseLevel;
312 }
313 return;
314}
315
316//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
317
318namespace {G4Mutex PenelopeRayleighModelMutex = G4MUTEX_INITIALIZER;}
319G4double G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
320 G4double energy,
321 G4double Z,
322 G4double,
323 G4double,
324 G4double)
325{
326 //Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97
327 //tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel
328 //et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.
329
330 if (fVerboseLevel > 3)
331 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModelMI" << G4endl;
332
333 G4int iZ = G4int(Z);
334 if (!fLogAtomicCrossSection[iZ]) {
335 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
336 //not filled up. This can happen in a UnitTest or via G4EmCalculator
337 if (fVerboseLevel > 0) {
338 //Issue a G4Exception (warning) only in verbose mode
339 G4ExceptionDescription ed;
340 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
341 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
342 G4Exception("G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom()",
343 "em2040",JustWarning,ed);
344 }
345
346 //protect file reading via autolock
347 G4AutoLock lock(&PenelopeRayleighModelMutex);
348 ReadDataFile(iZ);
349 lock.unlock();
350 }
351
352 G4double cross = 0;
353 G4PhysicsFreeVector* atom = fLogAtomicCrossSection[iZ];
354 if (!atom) {
355 G4ExceptionDescription ed;
356 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
357 G4Exception("G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom()",
358 "em2041",FatalException,ed);
359 return 0;
360 }
361
362 G4double logene = G4Log(energy);
363 G4double logXS = atom->Value(logene);
364 cross = G4Exp(logXS);
365
366 if (fVerboseLevel > 2) {
367 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z="
368 << Z << " = " << cross/barn << " barn" << G4endl;
369 }
370 return cross;
371}
372
373//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374
375void G4PenelopeRayleighModelMI::CalculateThetaAndAngFun()
376{
377 G4double theta = 0;
378 for(G4int k=0; k<fNtheta; k++) {
379 theta += fDTheta;
380 G4double value = (1+std::cos(theta)*std::cos(theta))*std::sin(theta);
381 fAngularFunction->PutValue(k,theta,value);
382 if (fVerboseLevel > 3)
383 G4cout << "theta[" << k << "]: " << fAngularFunction->Energy(k)
384 << ", angFun[" << k << "]: " << (*fAngularFunction)[k] << G4endl;
385 }
386}
387
388//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389
390G4double G4PenelopeRayleighModelMI::CalculateQSquared(G4double angle, G4double energy)
391{
392 G4double lambda,x,q,q2 = 0;
393
394 lambda = hbarc*twopi/energy;
395 x = 1./lambda*std::sin(angle/2.);
396 q = 2.*h_Planck*x/(electron_mass_c2/c_light);
397
398 if (fVerboseLevel > 3) {
399 G4cout << "E: " << energy/keV << " keV, lambda: " << lambda/nm << " nm"
400 << ", x: " << x*nm << ", q: " << q << G4endl;
401 }
402 q2 = std::pow(q,2);
403 return q2;
404}
405
406//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
407
408//Overriding of parent's (G4VEmModel) method
409G4double G4PenelopeRayleighModelMI::CrossSectionPerVolume(const G4Material* material,
410 const G4ParticleDefinition* p,
411 G4double energy,
412 G4double,
413 G4double)
414{
415 //check if we are in a Unit Test (only for the first time)
416 static G4bool amInAUnitTest = false;
417 if (G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize() == 0 && !amInAUnitTest)
418 {
419 amInAUnitTest = true;
420 G4ExceptionDescription ed;
421 ed << "The ProductionCuts table is empty " << G4endl;
422 ed << "This should happen only in Unit Tests" << G4endl;
423 G4Exception("G4PenelopeRayleighModelMI::CrossSectionPerVolume()",
424 "em2019",JustWarning,ed);
425 }
426 //If the material does not have a MIFF, continue with the old-style calculation
427 G4String matname = material->GetName();
428 if (amInAUnitTest)
429 {
430 //No need to lock, as this is always called in a master
431 const G4ElementVector* theElementVector = material->GetElementVector();
432 //protect file reading via autolock
433 for (std::size_t j=0;j<material->GetNumberOfElements();++j) {
434 G4int iZ = theElementVector->at(j)->GetZasInt();
435 if (!fLogAtomicCrossSection[iZ]) {
436 ReadDataFile(iZ);
437 }
438 }
439 if (fIsMIActive)
440 ReadMolInterferenceData(matname);
441 if (!fLogFormFactorTable->count(material))
442 BuildFormFactorTable(material);
443 if (!(fSamplingTable->count(material)))
444 InitializeSamplingAlgorithm(material);
445 if (!fPMaxTable->count(material))
446 GetPMaxTable(material);
447 }
448 G4bool useMIFF = fIsMIActive && (fMolInterferenceData->count(matname) || matname.find("MedMat") != std::string::npos);
449 if (!useMIFF)
450 {
451 if (fVerboseLevel > 2)
452 G4cout << "Rayleigh CS of: " << matname << " calculated through CSperAtom!" << G4endl;
453 return G4VEmModel::CrossSectionPerVolume(material,p,energy);
454 }
455
456 // If we are here, it means that we have to integrate the cross section
457 if (fVerboseLevel > 2)
458 G4cout << "Rayleigh CS of: " << matname
459 << " calculated through integration of the DCS" << G4endl;
460
461 G4double cs = 0;
462
463 //force null cross-section if below the low-energy edge of the table
464 if (energy < LowEnergyLimit())
465 return cs;
466
467 //if the material is a CRYSTAL, forbid this process
468 if (material->IsExtended() && material->GetName() != "CustomMat") {
469 G4ExtendedMaterial* extendedmaterial = (G4ExtendedMaterial*)material;
470 G4CrystalExtension* crystalExtension = (G4CrystalExtension*)extendedmaterial->RetrieveExtension("crystal");
471 if (crystalExtension != 0) {
472 G4cout << "The material has a crystalline structure, a dedicated diffraction model is used!" << G4endl;
473 return 0;
474 }
475 }
476
477 //GET MATERIAL INFORMATION
478 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
479 std::size_t nElements = material->GetNumberOfElements();
480 const G4ElementVector* elementVector = material->GetElementVector();
481 const G4double* fractionVector = material->GetFractionVector();
482
483 //Stoichiometric factors
484 std::vector<G4double> *StoichiometricFactors = new std::vector<G4double>;
485 for (std::size_t i=0;i<nElements;++i) {
486 G4double fraction = fractionVector[i];
487 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
488 StoichiometricFactors->push_back(fraction/atomicWeigth);
489 }
490 G4double MaxStoichiometricFactor = 0.;
491 for (std::size_t i=0;i<nElements;++i) {
492 if ((*StoichiometricFactors)[i] > MaxStoichiometricFactor)
493 MaxStoichiometricFactor = (*StoichiometricFactors)[i];
494 }
495 for (std::size_t i=0;i<nElements;++i) {
496 if (MaxStoichiometricFactor > 0.)
497 (*StoichiometricFactors)[i] /= MaxStoichiometricFactor;
498 }
499
500 //Equivalent atoms per molecule
501 G4double atPerMol = 0.;
502 for (std::size_t i=0;i<nElements;++i)
503 atPerMol += (*StoichiometricFactors)[i];
504 G4double moleculeDensity = 0.;
505 if (atPerMol != 0.) moleculeDensity = atomDensity/atPerMol;
506
507 if (fVerboseLevel > 2)
508 G4cout << "Material " << material->GetName() << " has " << atPerMol << " atoms "
509 << "per molecule and " << moleculeDensity/(cm*cm*cm) << " molecule/cm3" << G4endl;
510
511 //Equivalent molecular weight (dimensionless)
512 G4double MolWeight = 0.;
513 for (std::size_t i=0;i<nElements;++i)
514 MolWeight += (*StoichiometricFactors)[i]*(*elementVector)[i]->GetA()/(g/mole);
515
516 if (fVerboseLevel > 2)
517 G4cout << "Molecular weight of " << matname << ": " << MolWeight << " g/mol" << G4endl;
518
519 G4double IntegrandFun[fNtheta];
520 for (G4int k=0; k<fNtheta; k++) {
521 G4double theta = fAngularFunction->Energy(k); //the x-value is called "Energy", but is an angle
522 G4double F2 = GetFSquared(material,CalculateQSquared(theta,energy));
523 IntegrandFun[k] = (*fAngularFunction)[k]*F2;
524 }
525
526 G4double constant = pi*classic_electr_radius*classic_electr_radius;
527 cs = constant*IntegrateFun(IntegrandFun,fNtheta,fDTheta);
528
529 //Now cs is the cross section per molecule, let's calculate the cross section per volume
530 G4double csvolume = cs*moleculeDensity;
531
532 //print CS and mfp
533 if (fVerboseLevel > 2)
534 G4cout << "Rayleigh CS of " << matname << " at " << energy/keV
535 << " keV: " << cs/barn << " barn"
536 << ", mean free path: " << 1./csvolume/mm << " mm" << G4endl;
537
538 delete StoichiometricFactors;
539 //return CS
540 return csvolume;
541}
542
543//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
544
545void G4PenelopeRayleighModelMI::BuildFormFactorTable(const G4Material* material)
546{
547 if (fVerboseLevel > 3)
548 G4cout << "Calling BuildFormFactorTable() of G4PenelopeRayleighModelMI" << G4endl;
549
550 //GET MATERIAL INFORMATION
551 std::size_t nElements = material->GetNumberOfElements();
552 const G4ElementVector* elementVector = material->GetElementVector();
553 const G4double* fractionVector = material->GetFractionVector();
554
555 //Stoichiometric factors
556 std::vector<G4double> *StoichiometricFactors = new std::vector<G4double>;
557 for (std::size_t i=0;i<nElements;++i) {
558 G4double fraction = fractionVector[i];
559 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
560 StoichiometricFactors->push_back(fraction/atomicWeigth);
561 }
562 G4double MaxStoichiometricFactor = 0.;
563 for (std::size_t i=0;i<nElements;++i) {
564 if ((*StoichiometricFactors)[i] > MaxStoichiometricFactor)
565 MaxStoichiometricFactor = (*StoichiometricFactors)[i];
566 }
567 if (MaxStoichiometricFactor<1e-16) {
568 G4ExceptionDescription ed;
569 ed << "Inconsistent data of atomic composition for " << material->GetName() << G4endl;
570 G4Exception("G4PenelopeRayleighModelMI::BuildFormFactorTable()",
571 "em2042",FatalException,ed);
572 }
573 for (std::size_t i=0;i<nElements;++i)
574 (*StoichiometricFactors)[i] /= MaxStoichiometricFactor;
575
576 //Equivalent molecular weight (dimensionless)
577 G4double MolWeight = 0.;
578 for (std::size_t i=0;i<nElements;++i)
579 MolWeight += (*StoichiometricFactors)[i]*(*elementVector)[i]->GetA()/(g/mole);
580
581 //CREATE THE FORM FACTOR TABLE
582 //First, the form factors are retrieved [F/sqrt(W)].
583 //Then, they are squared and multiplied for MolWeight -> F2 [dimensionless].
584 //This makes difference for CS calculation, but not for theta sampling.
585 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(fLogQSquareGrid.size(),
586 /*spline=*/true);
587
588 G4String matname = material->GetName();
589 G4String aFileNameFF = "";
590
591 //retrieve MIdata (fFileNameFF)
592 G4MIData* dataMI = GetMIData(material);
593 if (dataMI)
594 aFileNameFF = dataMI->GetFilenameFF();
595
596 //read the MIFF from a file passed by the user
597 if (fIsMIActive && aFileNameFF != "") {
598 if (fVerboseLevel)
599 G4cout << "Read MIFF for " << matname << " from custom file: " << aFileNameFF << G4endl;
600
601 ReadMolInterferenceData(matname,aFileNameFF);
602 G4PhysicsFreeVector* Fvector = fMolInterferenceData->find(matname)->second;
603
604 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
605 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
606 G4double f = Fvector->Value(q);
607 G4double ff2 = f*f*MolWeight;
608 if (ff2)
609 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
610 }
611 }
612 //retrieve the MIFF from the database or use the IAM
613 else {
614 //medical material: composition of fat, water, bonematrix, mineral
615 if (fIsMIActive && matname.find("MedMat") != std::string::npos) {
616 if (fVerboseLevel)
617 G4cout << "Build MIFF from components for " << matname << G4endl;
618
619 //get the material composition from its name
620 G4int ki, kf=6, ktot=19;
621 G4double comp[4];
622 G4String compstring = matname.substr(kf+1, ktot);
623 for (std::size_t j=0; j<4; ++j) {
624 ki = kf+1;
625 kf = ki+4;
626 compstring = matname.substr(ki, 4);
627 comp[j] = atof(compstring.c_str());
628 if (fVerboseLevel > 2)
629 G4cout << " -- MedMat comp[" << j+1 << "]: " << comp[j] << G4endl;
630 }
631
632 const char* path = G4FindDataDir("G4LEDATA");
633 if (!path) {
634 G4String excep = "G4LEDATA environment variable not set!";
635 G4Exception("G4PenelopeRayleighModelMI::BuildFormFactorTable()",
636 "em0006",FatalException,excep);
637 }
638
639 if (!fMolInterferenceData->count("Fat_MI"))
640 ReadMolInterferenceData("Fat_MI");
641 G4PhysicsFreeVector* fatFF = fMolInterferenceData->find("Fat_MI")->second;
642
643 if (!fMolInterferenceData->count("Water_MI"))
644 ReadMolInterferenceData("Water_MI");
645 G4PhysicsFreeVector* waterFF = fMolInterferenceData->find("Water_MI")->second;
646
647 if (!fMolInterferenceData->count("BoneMatrix_MI"))
648 ReadMolInterferenceData("BoneMatrix_MI");
649 G4PhysicsFreeVector* bonematrixFF = fMolInterferenceData->find("BoneMatrix_MI")->second;
650
651 if (!fMolInterferenceData->count("Mineral_MI"))
652 ReadMolInterferenceData("Mineral_MI");
653 G4PhysicsFreeVector* mineralFF = fMolInterferenceData->find("Mineral_MI")->second;
654
655 //get and combine the molecular form factors with interference effect
656 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
657 G4double ff2 = 0;
658 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
659 G4double ffat = fatFF->Value(q);
660 G4double fwater = waterFF->Value(q);
661 G4double fbonematrix = bonematrixFF->Value(q);
662 G4double fmineral = mineralFF->Value(q);
663 ff2 = comp[0]*ffat*ffat+comp[1]*fwater*fwater+
664 comp[2]*fbonematrix*fbonematrix+comp[3]*fmineral*fmineral;
665 ff2 *= MolWeight;
666 if (ff2) theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
667 }
668 }
669 //other materials with MIFF (from the database)
670 else if (fIsMIActive && fMolInterferenceData->count(matname)) {
671 if (fVerboseLevel)
672 G4cout << "Read MIFF from database " << matname << G4endl;
673 G4PhysicsFreeVector* FF = fMolInterferenceData->find(matname)->second;
674 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
675 G4double ff2 = 0;
676 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
677 G4double f = FF->Value(q);
678 ff2 = f*f*MolWeight;
679 if (ff2) theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
680 }
681 }
682 //IAM
683 else {
684 if (fVerboseLevel)
685 G4cout << "FF of " << matname << " calculated according to the IAM" << G4endl;
686 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
687 G4double ff2 = 0;
688 for (std::size_t i=0;i<nElements;++i) {
689 G4int iZ = (*elementVector)[i]->GetZasInt();
690 G4PhysicsFreeVector* theAtomVec = fAtomicFormFactor[iZ];
691 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
692 G4double f = theAtomVec->Value(q);
693 ff2 += f*f*(*StoichiometricFactors)[i];
694 }
695 if (ff2) theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
696 }
697 }
698 }
699 theFFVec->FillSecondDerivatives();
700 fLogFormFactorTable->insert(std::make_pair(material,theFFVec));
701
702 if (fVerboseLevel > 3)
703 DumpFormFactorTable(material);
704 delete StoichiometricFactors;
705
706 return;
707}
708
709//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
710
711void G4PenelopeRayleighModelMI::SampleSecondaries(std::vector<G4DynamicParticle*>*,
712 const G4MaterialCutsCouple* couple,
713 const G4DynamicParticle* aDynamicGamma,
714 G4double,
715 G4double)
716{
717 // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
718 // from the Penelope2008 model. The scattering angle is sampled from the atomic
719 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
720 // anomalous scattering effects. The Form Factor F(Q) function which appears in the
721 // analytical cross section is retrieved via the method GetFSquared(); atomic data
722 // are tabulated for F(Q). Form factor for compounds is calculated according to
723 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
724 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
725 // for each material and managed by G4PenelopeSamplingData objects.
726 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
727 // increases with energy. For E=100 keV the efficiency is 100% and 86% for
728 // hydrogen and uranium, respectively.
729
730 if (fVerboseLevel > 3)
731 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModelMI" << G4endl;
732
733 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
734
735 if (photonEnergy0 <= fIntrinsicLowEnergyLimit) {
736 fParticleChange->ProposeTrackStatus(fStopAndKill);
737 fParticleChange->SetProposedKineticEnergy(0.);
738 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
739 return;
740 }
741
742 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
743 const G4Material* theMat = couple->GetMaterial();
744
745 //1) Verify if tables are ready
746 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
747 //not invoked
748 if (!fPMaxTable || !fSamplingTable || !fLogFormFactorTable) {
749 //create a **thread-local** version of the table. Used only for G4EmCalculator and
750 //Unit Tests
751 fLocalTable = true;
752 if (!fLogFormFactorTable)
753 fLogFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
754 if (!fPMaxTable)
755 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
756 if (!fSamplingTable)
757 fSamplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
758 if (fIsMIActive && !fMolInterferenceData)
759 fMolInterferenceData = new std::map<G4String,G4PhysicsFreeVector*>;
760 }
761
762 if (!fSamplingTable->count(theMat)) {
763 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
764 //not filled up. This can happen in a UnitTest
765 if (fVerboseLevel > 0) {
766 //Issue a G4Exception (warning) only in verbose mode
767 G4ExceptionDescription ed;
768 ed << "Unable to find the fSamplingTable data for " <<
769 theMat->GetName() << G4endl;
770 ed << "This can happen only in Unit Tests" << G4endl;
771 G4Exception("G4PenelopeRayleighModelMI::SampleSecondaries()",
772 "em2019",JustWarning,ed);
773 }
774 const G4ElementVector* theElementVector = theMat->GetElementVector();
775 //protect file reading via autolock
776 G4AutoLock lock(&PenelopeRayleighModelMutex);
777 for (std::size_t j=0;j<theMat->GetNumberOfElements();++j) {
778 G4int iZ = theElementVector->at(j)->GetZasInt();
779 if (!fLogAtomicCrossSection[iZ]) {
780 lock.lock();
781 ReadDataFile(iZ);
782 lock.unlock();
783 }
784 }
785 lock.lock();
786
787 //1) If the table has not been built for the material, do it!
788 if (!fLogFormFactorTable->count(theMat))
789 BuildFormFactorTable(theMat);
790
791 //2) retrieve or build the sampling table
792 if (!(fSamplingTable->count(theMat)))
793 InitializeSamplingAlgorithm(theMat);
794
795 //3) retrieve or build the pMax data
796 if (!fPMaxTable->count(theMat))
797 GetPMaxTable(theMat);
798
799 lock.unlock();
800 }
801
802 //Ok, restart the job
803 G4PenelopeSamplingData* theDataTable = fSamplingTable->find(theMat)->second;
804 G4PhysicsFreeVector* thePMax = fPMaxTable->find(theMat)->second;
805 G4double cosTheta = 1.0;
806
807 //OK, ready to go!
808 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
809
810 if (qmax < 1e-10) { //very low momentum transfer
811 G4bool loopAgain=false;
812 do {
813 loopAgain = false;
814 cosTheta = 1.0-2.0*G4UniformRand();
815 G4double G = 0.5*(1+cosTheta*cosTheta);
816 if (G4UniformRand()>G)
817 loopAgain = true;
818 } while(loopAgain);
819 }
820 else { //larger momentum transfer
821 std::size_t nData = theDataTable->GetNumberOfStoredPoints();
822 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
823 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
824
825 G4bool loopAgain = false;
826 G4double MaxPValue = thePMax->Value(photonEnergy0);
827 G4double xx=0;
828
829 //Sampling by rejection method. The rejection function is
830 //G = 0.5*(1+cos^2(theta))
831 do {
832 loopAgain = false;
833 G4double RandomMax = G4UniformRand()*MaxPValue;
834 xx = theDataTable->SampleValue(RandomMax);
835
836 //xx is a random value of q^2 in (0,q2max),sampled according to
837 //F(Q^2) via the RITA algorithm
838 if (xx > q2max)
839 loopAgain = true;
840 cosTheta = 1.0-2.0*xx/q2max;
841 G4double G = 0.5*(1+cosTheta*cosTheta);
842 if (G4UniformRand()>G)
843 loopAgain = true;
844 } while(loopAgain);
845 }
846
847 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
848
849 //Scattered photon angles. ( Z - axis along the parent photon)
850 G4double phi = twopi * G4UniformRand() ;
851 G4double dirX = sinTheta*std::cos(phi);
852 G4double dirY = sinTheta*std::sin(phi);
853 G4double dirZ = cosTheta;
854
855 //Update G4VParticleChange for the scattered photon
856 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
857
858 photonDirection1.rotateUz(photonDirection0);
859 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
860 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
861
862 return;
863}
864
865//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
866
867void G4PenelopeRayleighModelMI::ReadDataFile(const G4int Z)
868{
869 if (fVerboseLevel > 2) {
870 G4cout << "G4PenelopeRayleighModelMI::ReadDataFile()" << G4endl;
871 G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
872 }
873
874 const char* path = G4FindDataDir("G4LEDATA");
875 if (!path) {
876 G4String excep = "G4LEDATA environment variable not set!";
877 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
878 "em0006",FatalException,excep);
879 return;
880 }
881
882 //Read first the cross section file (all the files have 250 points)
883 std::ostringstream ost;
884 if (Z>9)
885 ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
886 else
887 ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
888 std::ifstream file(ost.str().c_str());
889
890 if (!file.is_open()) {
891 G4String excep = "Data file " + G4String(ost.str()) + " not found!";
892 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
893 "em0003",FatalException,excep);
894 }
895
896 G4int readZ = 0;
897 std::size_t nPoints = 0;
898 file >> readZ >> nPoints;
899
900 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) {
901 G4ExceptionDescription ed;
902 ed << "Corrupted data file for Z=" << Z << G4endl;
903 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
904 "em0005",FatalException,ed);
905 return;
906 }
907
908 fLogAtomicCrossSection[Z] = new G4PhysicsFreeVector((std::size_t)nPoints);
909 G4double ene=0,f1=0,f2=0,xs=0;
910 for (std::size_t i=0;i<nPoints;++i) {
911 file >> ene >> f1 >> f2 >> xs;
912 //dimensional quantities
913 ene *= eV;
914 xs *= cm2;
915 fLogAtomicCrossSection[Z]->PutValue(i,G4Log(ene),G4Log(xs));
916 if (file.eof() && i != (nPoints-1)) { //file ended too early
917 G4ExceptionDescription ed ;
918 ed << "Corrupted data file for Z=" << Z << G4endl;
919 ed << "Found less than " << nPoints << " entries" << G4endl;
920 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
921 "em0005",FatalException,ed);
922 }
923 }
924 file.close();
925
926 //Then, read the extended momentum transfer file
927 std::ostringstream ost2;
928 ost2 << path << "/penelope/rayleigh/MIFF/qext.dat";
929 file.open(ost2.str().c_str());
930
931 if (!file.is_open()) {
932 G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
933 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
934 "em0003",FatalException,excep);
935 }
936 G4bool fillQGrid = false;
937 if (!fLogQSquareGrid.size()) {
938 fillQGrid = true;
939 nPoints = 1142;
940 }
941 G4double qext = 0;
942 if (fillQGrid) { //fLogQSquareGrid filled only one time
943 for (std::size_t i=0;i<nPoints;++i) {
944 file >> qext;
945 fLogQSquareGrid.push_back(2.0*G4Log(qext));
946 }
947 }
948 file.close();
949
950 //Finally, read the form factor file
951 std::ostringstream ost3;
952 if (Z>9)
953 ost3 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
954 else
955 ost3 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
956 file.open(ost3.str().c_str());
957
958 if (!file.is_open()) {
959 G4String excep = "Data file " + G4String(ost3.str()) + " not found!";
960 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
961 "em0003",FatalException,excep);
962 }
963
964 file >> readZ >> nPoints;
965
966 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) {
967 G4ExceptionDescription ed;
968 ed << "Corrupted data file for Z=" << Z << G4endl;
969 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
970 "em0005",FatalException,ed);
971 return;
972 }
973
974 fAtomicFormFactor[Z] = new G4PhysicsFreeVector((std::size_t)nPoints);
975 G4double q=0,ff=0,incoh=0;
976 for (std::size_t i=0;i<nPoints;++i) {
977 file >> q >> ff >> incoh; //q and ff are dimensionless (q is in units of (m_e*c))
978 fAtomicFormFactor[Z]->PutValue(i,q,ff);
979 if (file.eof() && i != (nPoints-1)) { //file ended too early
980 G4ExceptionDescription ed;
981 ed << "Corrupted data file for Z=" << Z << G4endl;
982 ed << "Found less than " << nPoints << " entries" << G4endl;
983 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
984 "em0005",FatalException,ed);
985 }
986 }
987 file.close();
988 return;
989}
990
991//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
992
993void G4PenelopeRayleighModelMI::ReadMolInterferenceData(const G4String& matname,
994 const G4String& FFfilename)
995{
996 if (fVerboseLevel > 2) {
997 G4cout << "G4PenelopeRayleighModelMI::ReadMolInterferenceData() for material " <<
998 matname << G4endl;
999 }
1000 G4bool isLocalFile = (FFfilename != "NULL");
1001
1002 const char* path = G4FindDataDir("G4LEDATA");
1003 if (!path) {
1004 G4String excep = "G4LEDATA environment variable not set!";
1005 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1006 "em0006",FatalException,excep);
1007 return;
1008 }
1009
1010 if (!(fKnownMaterials->count(matname)) && !isLocalFile) //material not found
1011 return;
1012
1013 G4String aFileName = (isLocalFile) ? FFfilename : fKnownMaterials->find(matname)->second;
1014
1015 //if the material has a MIFF, read it from the database
1016 if (aFileName != "") {
1017 if (fVerboseLevel > 2)
1018 G4cout << "ReadMolInterferenceData(). Read material: " << matname << ", filename: " <<
1019 aFileName << " " <<
1020 (isLocalFile ? "(local)" : "(database)") << G4endl;
1021 std::ifstream file;
1022 std::ostringstream ostIMFF;
1023 if (isLocalFile)
1024 ostIMFF << aFileName;
1025 else
1026 ostIMFF << path << "/penelope/rayleigh/MIFF/" << aFileName;
1027 file.open(ostIMFF.str().c_str());
1028
1029 if (!file.is_open()) {
1030 G4String excep = "Data file " + G4String(ostIMFF.str()) + " not found!";
1031 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1032 "em1031",FatalException,excep);
1033 return;
1034 }
1035
1036 //check the number of points in the file
1037 std::size_t nPoints = 0;
1038 G4double x=0,y=0;
1039 while (file.good()) {
1040 file >> x >> y;
1041 nPoints++;
1042 }
1043 file.close();
1044 nPoints--;
1045 if (fVerboseLevel > 3)
1046 G4cout << "Number of nPoints: " << nPoints << G4endl;
1047
1048 //read the file
1049 file.open(ostIMFF.str().c_str());
1050 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((std::size_t)nPoints);
1051 G4double q=0,ff=0;
1052 for (std::size_t i=0; i<nPoints; ++i) {
1053 file >> q >> ff;
1054
1055 //q and ff are dimensionless (q is in units of (m_e*c))
1056 theFFVec->PutValue(i,q,ff);
1057
1058 //file ended too early
1059 if (file.eof() && i != (nPoints-1)) {
1060 G4ExceptionDescription ed;
1061 ed << "Corrupted data file" << G4endl;
1062 ed << "Found less than " << nPoints << " entries" << G4endl;
1063 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1064 "em1005",FatalException,ed);
1065 }
1066 }
1067 if (!fMolInterferenceData) {
1068 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1069 "em2145",FatalException,
1070 "Unable to allocate the Molecular Interference data table");
1071 delete theFFVec;
1072 return;
1073 }
1074 file.close();
1075 fMolInterferenceData->insert(std::make_pair(matname,theFFVec));
1076 }
1077 return;
1078}
1079
1080//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1081
1082G4double G4PenelopeRayleighModelMI::GetFSquared(const G4Material* mat, const G4double QSquared)
1083{
1084 G4double f2 = 0;
1085 //Input value QSquared could be zero: protect the log() below against
1086 //the FPE exception
1087
1088 //If Q<1e-10, set Q to 1e-10
1089 G4double logQSquared = (QSquared>1e-10) ? G4Log(QSquared) : -23.;
1090 //last value of the table
1091 G4double maxlogQ2 = fLogQSquareGrid[fLogQSquareGrid.size()-1];
1092
1093 //now it should be all right
1094 G4PhysicsFreeVector* theVec = fLogFormFactorTable->find(mat)->second;
1095
1096 if (!theVec) {
1097 G4ExceptionDescription ed;
1098 ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
1099 G4Exception("G4PenelopeRayleighModelMI::GetFSquared()",
1100 "em2046",FatalException,ed);
1101 return 0;
1102 }
1103
1104 if (logQSquared < -20) { //Q < 1e-9
1105 G4double logf2 = (*theVec)[0]; //first value of the table
1106 f2 = G4Exp(logf2);
1107 }
1108 else if (logQSquared > maxlogQ2)
1109 f2 =0;
1110 else {
1111 //log(Q^2) vs. log(F^2)
1112 G4double logf2 = theVec->Value(logQSquared);
1113 f2 = G4Exp(logf2);
1114 }
1115
1116 if (fVerboseLevel > 3) {
1117 G4cout << "G4PenelopeRayleighModelMI::GetFSquared() in " << mat->GetName() << G4endl;
1118 G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c)); F^2 = " << f2 << G4endl;
1119 }
1120 return f2;
1121}
1122
1123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1124
1125void G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm(const G4Material* mat)
1126{
1127 G4double q2min = 0;
1128 G4double q2max = 0;
1129 const std::size_t np = 150; //hard-coded in Penelope
1130 for (std::size_t i=1;i<fLogQSquareGrid.size();++i)
1131 {
1132 G4double Q2 = G4Exp(fLogQSquareGrid[i]);
1133 if (GetFSquared(mat,Q2) > 1e-35)
1134 {
1135 q2max = G4Exp(fLogQSquareGrid[i-1]);
1136 }
1137 }
1138 std::size_t nReducedPoints = np/4;
1139
1140 //check for errors
1141 if (np < 16)
1142 {
1143 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1144 "em2047",FatalException,
1145 "Too few points to initialize the sampling algorithm");
1146 }
1147 if (q2min > (q2max-1e-10))
1148 {
1149 G4cout << "q2min= " << q2min << " q2max= " << q2max << G4endl;
1150 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1151 "em2048",FatalException,
1152 "Too narrow grid to initialize the sampling algorithm");
1153 }
1154
1155 //This is subroutine RITAI0 of Penelope
1156 //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material*
1157
1158 //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
1159 G4DataVector* x = new G4DataVector();
1160
1161 /*******************************************************************************
1162 Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
1163 ********************************************************************************/
1164 std::size_t NUNIF = std::min(std::max(((std::size_t)8),nReducedPoints),np/2);
1165 const G4int nip = 51; //hard-coded in Penelope
1166
1167 G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
1168 x->push_back(q2min);
1169 for (std::size_t i=1;i<NUNIF-1;++i)
1170 {
1171 G4double app = q2min + i*dx;
1172 x->push_back(app); //increase
1173 }
1174 x->push_back(q2max);
1175
1176 if (fVerboseLevel> 3)
1177 G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
1178
1179 //Allocate temporary storage vectors
1180 G4DataVector* area = new G4DataVector();
1181 G4DataVector* a = new G4DataVector();
1182 G4DataVector* b = new G4DataVector();
1183 G4DataVector* c = new G4DataVector();
1184 G4DataVector* err = new G4DataVector();
1185
1186 for (std::size_t i=0;i<NUNIF-1;++i) //build all points but the last
1187 {
1188 //Temporary vectors for this loop
1189 G4DataVector* pdfi = new G4DataVector();
1190 G4DataVector* pdfih = new G4DataVector();
1191 G4DataVector* sumi = new G4DataVector();
1192
1193 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
1194 G4double pdfmax = 0;
1195 for (G4int k=0;k<nip;k++)
1196 {
1197 G4double xik = (*x)[i]+k*dxi;
1198 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1199 pdfi->push_back(pdfk);
1200 pdfmax = std::max(pdfmax,pdfk);
1201 if (k < (nip-1))
1202 {
1203 G4double xih = xik + 0.5*dxi;
1204 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1205 pdfih->push_back(pdfIK);
1206 pdfmax = std::max(pdfmax,pdfIK);
1207 }
1208 }
1209
1210 //Simpson's integration
1211 G4double cons = dxi*0.5*(1./3.);
1212 sumi->push_back(0.);
1213 for (G4int k=1;k<nip;k++)
1214 {
1215 G4double previous = (*sumi)[k-1];
1216 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1217 sumi->push_back(next);
1218 }
1219
1220 G4double lastIntegral = (*sumi)[sumi->size()-1];
1221 area->push_back(lastIntegral);
1222 //Normalize cumulative function
1223 G4double factor = 1.0/lastIntegral;
1224 for (std::size_t k=0;k<sumi->size();++k)
1225 (*sumi)[k] *= factor;
1226
1227 //When the PDF vanishes at one of the interval end points, its value is modified
1228 if ((*pdfi)[0] < 1e-35)
1229 (*pdfi)[0] = 1e-5*pdfmax;
1230 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1231 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1232
1233 G4double pli = (*pdfi)[0]*factor;
1234 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1235 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
1236 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
1237 G4double C_temp = 1.0+A_temp+B_temp;
1238 if (C_temp < 1e-35)
1239 {
1240 a->push_back(0.);
1241 b->push_back(0.);
1242 c->push_back(1.);
1243 }
1244 else
1245 {
1246 a->push_back(A_temp);
1247 b->push_back(B_temp);
1248 c->push_back(C_temp);
1249 }
1250
1251 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1252 //and the true pdf, extended over the interval (X(I),X(I+1))
1253 G4int icase = 1; //loop code
1254 G4bool reLoop = false;
1255 err->push_back(0.);
1256 do
1257 {
1258 reLoop = false;
1259 (*err)[i] = 0.; //zero variable
1260 for (G4int k=0;k<nip;k++)
1261 {
1262 G4double rr = (*sumi)[k];
1263 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1264 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1265 if (k == 0 || k == nip-1)
1266 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1267 else
1268 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1269 }
1270 (*err)[i] *= dxi;
1271
1272 //If err(I) is too large, the pdf is approximated by a uniform distribution
1273 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1274 {
1275 (*b)[i] = 0;
1276 (*a)[i] = 0;
1277 (*c)[i] = 1.;
1278 icase = 2;
1279 reLoop = true;
1280 }
1281 }while(reLoop);
1282
1283 delete pdfi;
1284 delete pdfih;
1285 delete sumi;
1286 } //end of first loop over i
1287
1288 //Now assign last point
1289 (*x)[x->size()-1] = q2max;
1290 a->push_back(0.);
1291 b->push_back(0.);
1292 c->push_back(0.);
1293 err->push_back(0.);
1294 area->push_back(0.);
1295
1296 if (x->size() != NUNIF || a->size() != NUNIF ||
1297 err->size() != NUNIF || area->size() != NUNIF)
1298 {
1299 G4ExceptionDescription ed;
1300 ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
1301 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1302 "em2049",FatalException,ed);
1303 }
1304
1305 /*******************************************************************************
1306 New grid points are added by halving the sub-intervals with the largest absolute error
1307 This is done up to np=150 points in the grid
1308 ********************************************************************************/
1309 do
1310 {
1311 G4double maxError = 0.0;
1312 std::size_t iErrMax = 0;
1313 for (std::size_t i=0;i<err->size()-2;++i)
1314 {
1315 //maxError is the lagest of the interval errors err[i]
1316 if ((*err)[i] > maxError)
1317 {
1318 maxError = (*err)[i];
1319 iErrMax = i;
1320 }
1321 }
1322
1323 //OK, now I have to insert one new point in the position iErrMax
1324 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
1325
1326 x->insert(x->begin()+iErrMax+1,newx);
1327 //Add place-holders in the other vectors
1328 area->insert(area->begin()+iErrMax+1,0.);
1329 a->insert(a->begin()+iErrMax+1,0.);
1330 b->insert(b->begin()+iErrMax+1,0.);
1331 c->insert(c->begin()+iErrMax+1,0.);
1332 err->insert(err->begin()+iErrMax+1,0.);
1333
1334 //Now calculate the other parameters
1335 for (std::size_t i=iErrMax;i<=iErrMax+1;++i)
1336 {
1337 //Temporary vectors for this loop
1338 G4DataVector* pdfi = new G4DataVector();
1339 G4DataVector* pdfih = new G4DataVector();
1340 G4DataVector* sumi = new G4DataVector();
1341
1342 G4double dxLocal = (*x)[i+1]-(*x)[i];
1343 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
1344 G4double pdfmax = 0;
1345 for (G4int k=0;k<nip;k++)
1346 {
1347 G4double xik = (*x)[i]+k*dxi;
1348 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1349 pdfi->push_back(pdfk);
1350 pdfmax = std::max(pdfmax,pdfk);
1351 if (k < (nip-1))
1352 {
1353 G4double xih = xik + 0.5*dxi;
1354 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1355 pdfih->push_back(pdfIK);
1356 pdfmax = std::max(pdfmax,pdfIK);
1357 }
1358 }
1359
1360 //Simpson's integration
1361 G4double cons = dxi*0.5*(1./3.);
1362 sumi->push_back(0.);
1363 for (G4int k=1;k<nip;k++)
1364 {
1365 G4double previous = (*sumi)[k-1];
1366 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1367 sumi->push_back(next);
1368 }
1369 G4double lastIntegral = (*sumi)[sumi->size()-1];
1370 (*area)[i] = lastIntegral;
1371
1372 //Normalize cumulative function
1373 G4double factor = 1.0/lastIntegral;
1374 for (std::size_t k=0;k<sumi->size();++k)
1375 (*sumi)[k] *= factor;
1376
1377 //When the PDF vanishes at one of the interval end points, its value is modified
1378 if ((*pdfi)[0] < 1e-35)
1379 (*pdfi)[0] = 1e-5*pdfmax;
1380 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1381 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1382
1383 G4double pli = (*pdfi)[0]*factor;
1384 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1385 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
1386 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
1387 G4double C_temp = 1.0+A_temp+B_temp;
1388 if (C_temp < 1e-35)
1389 {
1390 (*a)[i]= 0.;
1391 (*b)[i] = 0.;
1392 (*c)[i] = 1;
1393 }
1394 else
1395 {
1396 (*a)[i]= A_temp;
1397 (*b)[i] = B_temp;
1398 (*c)[i] = C_temp;
1399 }
1400 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1401 //and the true pdf, extended over the interval (X(I),X(I+1))
1402 G4int icase = 1; //loop code
1403 G4bool reLoop = false;
1404 do
1405 {
1406 reLoop = false;
1407 (*err)[i] = 0.; //zero variable
1408 for (G4int k=0;k<nip;k++)
1409 {
1410 G4double rr = (*sumi)[k];
1411 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1412 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1413 if (k == 0 || k == nip-1)
1414 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1415 else
1416 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1417 }
1418 (*err)[i] *= dxi;
1419
1420 //If err(I) is too large, the pdf is approximated by a uniform distribution
1421 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1422 {
1423 (*b)[i] = 0;
1424 (*a)[i] = 0;
1425 (*c)[i] = 1.;
1426 icase = 2;
1427 reLoop = true;
1428 }
1429 }while(reLoop);
1430 delete pdfi;
1431 delete pdfih;
1432 delete sumi;
1433 }
1434 }while(x->size() < np);
1435
1436 if (x->size() != np || a->size() != np ||
1437 err->size() != np || area->size() != np)
1438 {
1439 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1440 "em2050",FatalException,
1441 "Problem in building the extended Table for Sampling: array dimensions do not match ");
1442 }
1443
1444 /*******************************************************************************
1445 Renormalization
1446 ********************************************************************************/
1447 G4double ws = 0;
1448 for (std::size_t i=0;i<np-1;++i)
1449 ws += (*area)[i];
1450 ws = 1.0/ws;
1451 G4double errMax = 0;
1452 for (std::size_t i=0;i<np-1;++i)
1453 {
1454 (*area)[i] *= ws;
1455 (*err)[i] *= ws;
1456 errMax = std::max(errMax,(*err)[i]);
1457 }
1458
1459 //Vector with the normalized cumulative distribution
1460 G4DataVector* PAC = new G4DataVector();
1461 PAC->push_back(0.);
1462 for (std::size_t i=0;i<np-1;++i)
1463 {
1464 G4double previous = (*PAC)[i];
1465 PAC->push_back(previous+(*area)[i]);
1466 }
1467 (*PAC)[PAC->size()-1] = 1.;
1468
1469 /*******************************************************************************
1470 Pre-calculated limits for the initial binary search for subsequent sampling
1471 ********************************************************************************/
1472 std::vector<std::size_t> *ITTL = new std::vector<std::size_t>;
1473 std::vector<std::size_t> *ITTU = new std::vector<std::size_t>;
1474
1475 //Just create place-holders
1476 for (std::size_t i=0;i<np;++i)
1477 {
1478 ITTL->push_back(0);
1479 ITTU->push_back(0);
1480 }
1481
1482 G4double bin = 1.0/(np-1);
1483 (*ITTL)[0]=0;
1484 for (std::size_t i=1;i<(np-1);++i)
1485 {
1486 G4double ptst = i*bin;
1487 G4bool found = false;
1488 for (std::size_t j=(*ITTL)[i-1];j<np && !found;++j)
1489 {
1490 if ((*PAC)[j] > ptst)
1491 {
1492 (*ITTL)[i] = j-1;
1493 (*ITTU)[i-1] = j;
1494 found = true;
1495 }
1496 }
1497 }
1498 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1499 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1500 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1501
1502 if (ITTU->size() != np || ITTU->size() != np)
1503 {
1504 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1505 "em2051",FatalException,
1506 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1507 }
1508
1509 /********************************************************************************
1510 Copy tables
1511 ********************************************************************************/
1512 G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
1513 for (std::size_t i=0;i<np;++i)
1514 {
1515 theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1516 }
1517 if (fVerboseLevel > 2)
1518 {
1519 G4cout << "*************************************************************************" <<
1520 G4endl;
1521 G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
1522 theTable->DumpTable();
1523 }
1524 fSamplingTable->insert(std::make_pair(mat,theTable));
1525
1526 //Clean up temporary vectors
1527 delete x;
1528 delete a;
1529 delete b;
1530 delete c;
1531 delete err;
1532 delete area;
1533 delete PAC;
1534 delete ITTL;
1535 delete ITTU;
1536
1537 //DONE!
1538 return;
1539}
1540
1541//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1542
1543void G4PenelopeRayleighModelMI::GetPMaxTable(const G4Material* mat)
1544{
1545 if (!fPMaxTable)
1546 {
1547 G4cout << "G4PenelopeRayleighModelMI::BuildPMaxTable" << G4endl;
1548 G4cout << "Going to instanziate the fPMaxTable !" << G4endl;
1549 G4cout << "That should _not_ be here! " << G4endl;
1550 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
1551 }
1552 //check if the table is already there
1553 if (fPMaxTable->count(mat))
1554 return;
1555
1556 //otherwise build it
1557 if (!fSamplingTable)
1558 {
1559 G4Exception("G4PenelopeRayleighModelMI::GetPMaxTable()",
1560 "em2052",FatalException,
1561 "SamplingTable is not properly instantiated");
1562 return;
1563 }
1564
1565 //This should not be: the sampling table is built before the p-table
1566 if (!fSamplingTable->count(mat))
1567 {
1568 G4ExceptionDescription ed;
1569 ed << "Sampling table for material " << mat->GetName() << " not found";
1570 G4Exception("G4PenelopeRayleighModelMI::GetPMaxTable()",
1571 "em2052",FatalException,
1572 ed);
1573 return;
1574 }
1575
1576 G4PenelopeSamplingData *theTable = fSamplingTable->find(mat)->second;
1577 std::size_t tablePoints = theTable->GetNumberOfStoredPoints();
1578 std::size_t nOfEnergyPoints = fLogEnergyGridPMax.size();
1579 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
1580
1581 const std::size_t nip = 51; //hard-coded in Penelope
1582
1583 for (std::size_t ie=0;ie<fLogEnergyGridPMax.size();++ie)
1584 {
1585 G4double energy = G4Exp(fLogEnergyGridPMax[ie]);
1586 G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
1587 G4double Qm2 = Qm*Qm;
1588 G4double firstQ2 = theTable->GetX(0);
1589 G4double lastQ2 = theTable->GetX(tablePoints-1);
1590 G4double thePMax = 0;
1591
1592 if (Qm2 > firstQ2)
1593 {
1594 if (Qm2 < lastQ2)
1595 {
1596 //bisection to look for the index of Qm
1597 std::size_t lowerBound = 0;
1598 std::size_t upperBound = tablePoints-1;
1599 while (lowerBound <= upperBound)
1600 {
1601 std::size_t midBin = (lowerBound + upperBound)/2;
1602 if( Qm2 < theTable->GetX(midBin))
1603 { upperBound = midBin-1; }
1604 else
1605 { lowerBound = midBin+1; }
1606 }
1607 //upperBound is the output (but also lowerBounf --> should be the same!)
1608 G4double Q1 = theTable->GetX(upperBound);
1609 G4double Q2 = Qm2;
1610 G4double DQ = (Q2-Q1)/((G4double)(nip-1));
1611 G4double theA = theTable->GetA(upperBound);
1612 G4double theB = theTable->GetB(upperBound);
1613 G4double thePAC = theTable->GetPAC(upperBound);
1614 G4DataVector* fun = new G4DataVector();
1615 for (std::size_t k=0;k<nip;++k)
1616 {
1617 G4double qi = Q1 + k*DQ;
1618 G4double tau = (qi-Q1)/
1619 (theTable->GetX(upperBound+1)-Q1);
1620 G4double con1 = 2.0*theB*tau;
1621 G4double ci = 1.0+theA+theB;
1622 G4double con2 = ci-theA*tau;
1623 G4double etap = 0;
1624 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1625 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1626 else
1627 etap = tau/con2;
1628 G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
1629 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1630 ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
1631 fun->push_back(theFun);
1632 }
1633 //Now intergrate numerically the fun Cavalieri-Simpson's method
1634 G4DataVector* sum = new G4DataVector;
1635 G4double CONS = DQ*(1./12.);
1636 G4double HCONS = 0.5*CONS;
1637 sum->push_back(0.);
1638 G4double secondPoint = (*sum)[0] +
1639 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1640 sum->push_back(secondPoint);
1641 for (std::size_t hh=2;hh<nip-1;++hh)
1642 {
1643 G4double previous = (*sum)[hh-1];
1644 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1645 (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1646 sum->push_back(next);
1647 }
1648 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1649 (*fun)[nip-3])*CONS;
1650 sum->push_back(last);
1651 thePMax = thePAC + (*sum)[sum->size()-1]; //last point
1652 delete fun;
1653 delete sum;
1654 }
1655 else
1656 {
1657 thePMax = 1.0;
1658 }
1659 }
1660 else
1661 {
1662 thePMax = theTable->GetPAC(0);
1663 }
1664
1665 //Write number in the table
1666 theVec->PutValue(ie,energy,thePMax);
1667 }
1668
1669 fPMaxTable->insert(std::make_pair(mat,theVec));
1670 return;
1671}
1672
1673//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1674
1675void G4PenelopeRayleighModelMI::DumpFormFactorTable(const G4Material* mat)
1676{
1677 G4cout << "*****************************************************************" << G4endl;
1678 G4cout << "G4PenelopeRayleighModelMI: Form Factor Table for " << mat->GetName() << G4endl;
1679 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1680 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1681 G4cout << "*****************************************************************" << G4endl;
1682 if (!fLogFormFactorTable->count(mat))
1683 BuildFormFactorTable(mat);
1684
1685 G4PhysicsFreeVector* theVec = fLogFormFactorTable->find(mat)->second;
1686 for (std::size_t i=0;i<theVec->GetVectorLength();++i)
1687 {
1688 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1689 G4double Q = G4Exp(0.5*logQ2);
1690 G4double logF2 = (*theVec)[i];
1691 G4double F = G4Exp(0.5*logF2);
1692 G4cout << Q << " " << F << G4endl;
1693 }
1694 //DONE
1695 return;
1696}
1697
1698//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1699
1700void G4PenelopeRayleighModelMI::SetParticle(const G4ParticleDefinition* p)
1701{
1702 if(!fParticle) {
1703 fParticle = p;
1704 }
1705}
1706
1707//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1708void G4PenelopeRayleighModelMI::LoadKnownMIFFMaterials()
1709{
1710 fKnownMaterials->insert(std::pair<G4String,G4String>("Fat_MI","FF_fat_Tartari2002.dat"));
1711 fKnownMaterials->insert(std::pair<G4String,G4String>("Water_MI","FF_water_Tartari2002.dat"));
1712 fKnownMaterials->insert(std::pair<G4String,G4String>("BoneMatrix_MI","FF_bonematrix_Tartari2002.dat"));
1713 fKnownMaterials->insert(std::pair<G4String,G4String>("Mineral_MI","FF_mineral_Tartari2002.dat"));
1714 fKnownMaterials->insert(std::pair<G4String,G4String>("adipose_MI","FF_adipose_Poletti2002.dat"));
1715 fKnownMaterials->insert(std::pair<G4String,G4String>("glandular_MI","FF_glandular_Poletti2002.dat"));
1716 fKnownMaterials->insert(std::pair<G4String,G4String>("breast5050_MI","FF_human_breast_Peplow1998.dat"));
1717 fKnownMaterials->insert(std::pair<G4String,G4String>("carcinoma_MI","FF_carcinoma_Kidane1999.dat"));
1718 fKnownMaterials->insert(std::pair<G4String,G4String>("muscle_MI","FF_pork_muscle_Peplow1998.dat"));
1719 fKnownMaterials->insert(std::pair<G4String,G4String>("kidney_MI","FF_pork_kidney_Peplow1998.dat"));
1720 fKnownMaterials->insert(std::pair<G4String,G4String>("liver_MI","FF_pork_liver_Peplow1998.dat"));
1721 fKnownMaterials->insert(std::pair<G4String,G4String>("heart_MI","FF_pork_heart_Peplow1998.dat"));
1722 fKnownMaterials->insert(std::pair<G4String,G4String>("blood_MI","FF_beef_blood_Peplow1998.dat"));
1723 fKnownMaterials->insert(std::pair<G4String,G4String>("grayMatter_MI","FF_gbrain_DeFelici2008.dat"));
1724 fKnownMaterials->insert(std::pair<G4String,G4String>("whiteMatter_MI","FF_wbrain_DeFelici2008.dat"));
1725 fKnownMaterials->insert(std::pair<G4String,G4String>("bone_MI","FF_bone_King2011.dat"));
1726 fKnownMaterials->insert(std::pair<G4String,G4String>("FatLowX_MI","FF_fat_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1727 fKnownMaterials->insert(std::pair<G4String,G4String>("BoneMatrixLowX_MI","FF_bonematrix_Tartari2002_joint_lowXdata.dat"));
1728 fKnownMaterials->insert(std::pair<G4String,G4String>("PMMALowX_MI", "FF_PMMA_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1729 fKnownMaterials->insert(std::pair<G4String,G4String>("dryBoneLowX_MI","FF_drybone_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1730 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS30-70_MI","FF_CIRS30-70_Poletti2002.dat"));
1731 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS50-50_MI","FF_CIRS50-50_Poletti2002.dat"));
1732 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS70-30_MI","FF_CIRS70-30_Poletti2002.dat"));
1733 fKnownMaterials->insert(std::pair<G4String,G4String>("RMI454_MI", "FF_RMI454_Poletti2002.dat"));
1734 fKnownMaterials->insert(std::pair<G4String,G4String>("PMMA_MI","FF_PMMA_Tartari2002.dat"));
1735 fKnownMaterials->insert(std::pair<G4String,G4String>("Lexan_MI","FF_lexan_Peplow1998.dat"));
1736 fKnownMaterials->insert(std::pair<G4String,G4String>("Kapton_MI","FF_kapton_Peplow1998.dat"));
1737 fKnownMaterials->insert(std::pair<G4String,G4String>("Nylon_MI","FF_nylon_Kosanetzky1987.dat"));
1738 fKnownMaterials->insert(std::pair<G4String,G4String>("Polyethylene_MI","FF_polyethylene_Kosanetzky1987.dat"));
1739 fKnownMaterials->insert(std::pair<G4String,G4String>("Polystyrene_MI","FF_polystyrene_Kosanetzky1987.dat"));
1740 fKnownMaterials->insert(std::pair<G4String,G4String>("Formaline_MI","FF_formaline_Peplow1998.dat"));
1741 fKnownMaterials->insert(std::pair<G4String,G4String>("Acetone_MI","FF_acetone_Cozzini2010.dat"));
1742 fKnownMaterials->insert(std::pair<G4String,G4String>("Hperoxide_MI","FF_Hperoxide_Cozzini2010.dat"));
1743}
1744
1745//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1746
1747G4double G4PenelopeRayleighModelMI::IntegrateFun(G4double y[], G4int n, G4double dTheta)
1748{
1749 G4double integral = 0.;
1750 for (G4int k=0; k<n-1; k++) {
1751 integral += (y[k]+y[k+1]);
1752 }
1753 integral *= dTheta*0.5;
1754 return integral;
1755}
1756
1757//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1758G4MIData* G4PenelopeRayleighModelMI::GetMIData(const G4Material* material)
1759{
1760 if (material->IsExtended()) {
1761 G4ExtendedMaterial* aEM = (G4ExtendedMaterial*)material;
1762 G4MIData* dataMI = (G4MIData*)aEM->RetrieveExtension("MI");
1763 return dataMI;
1764 } else {
1765 return nullptr;
1766 }
1767}
G4AutoLock
G4TemplateAutoLock< G4Mutex > G4AutoLock
Definition G4AutoLock.hh:582
G4ElementVector
std::vector< const G4Element * > G4ElementVector
Definition G4ElementVector.hh:33
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:169
G4ParticleMomentum
G4ThreeVector G4ParticleMomentum
Definition G4ParticleMomentum.hh:39
G4MUTEX_INITIALIZER
#define G4MUTEX_INITIALIZER
Definition G4Threading.hh:85
G4Mutex
std::mutex G4Mutex
Definition G4Threading.hh:81
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition ThreeVector.cc:33
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition G4EmParameters.cc:71
G4EmParameters::Verbose
G4int Verbose() const
Definition G4EmParameters.cc:1049
G4ExtendedMaterial::RetrieveExtension
G4VMaterialExtension * RetrieveExtension(const G4String &name)
Definition G4ExtendedMaterial.cc:81
G4MIData::GetFilenameFF
const G4String & GetFilenameFF()
Definition G4MIData.hh:52
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition G4Material.hh:183
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition G4Material.hh:197
G4Material::IsExtended
virtual G4bool IsExtended() const
Definition G4Material.cc:799
G4Material::GetFractionVector
const G4double * GetFractionVector() const
Definition G4Material.hh:186
G4Material::GetNumberOfElements
std::size_t GetNumberOfElements() const
Definition G4Material.hh:180
G4Material::GetName
const G4String & GetName() const
Definition G4Material.hh:171
G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
Definition G4PenelopeRayleighModelMI.cc:319
G4PenelopeRayleighModelMI::InitialiseLocal
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
Definition G4PenelopeRayleighModelMI.cc:280
G4PenelopeRayleighModelMI::SampleSecondaries
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
Definition G4PenelopeRayleighModelMI.cc:711
G4PenelopeRayleighModelMI::CrossSectionPerVolume
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0., G4double maxEnergy=DBL_MAX) override
Definition G4PenelopeRayleighModelMI.cc:409
G4PenelopeRayleighModelMI::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition G4PenelopeRayleighModelMI.cc:180
G4PenelopeRayleighModelMI::DumpFormFactorTable
void DumpFormFactorTable(const G4Material *)
Definition G4PenelopeRayleighModelMI.cc:1675
G4PenelopeRayleighModelMI::G4PenelopeRayleighModelMI
G4PenelopeRayleighModelMI(const G4ParticleDefinition *p=nullptr, const G4String &processName="PenRayleighMI")
Definition G4PenelopeRayleighModelMI.cc:72
G4PenelopeSamplingData::GetA
G4double GetA(size_t index)
Definition G4PenelopeSamplingData.cc:153
G4PenelopeSamplingData::GetPAC
G4double GetPAC(size_t index)
Definition G4PenelopeSamplingData.cc:144
G4PenelopeSamplingData::AddPoint
void AddPoint(G4double x0, G4double pac0, G4double a0, G4double b0, size_t ITTL0, size_t ITTU0)
Definition G4PenelopeSamplingData.cc:97
G4PenelopeSamplingData::GetX
G4double GetX(size_t index)
Definition G4PenelopeSamplingData.cc:135
G4PenelopeSamplingData::SampleValue
G4double SampleValue(G4double rndm)
Definition G4PenelopeSamplingData.cc:171
G4PenelopeSamplingData::GetB
G4double GetB(size_t index)
Definition G4PenelopeSamplingData.cc:162
G4PhysicsFreeVector::PutValue
void PutValue(const std::size_t index, const G4double e, const G4double value)
Definition G4PhysicsFreeVector.cc:180
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(const std::size_t index) const
G4PhysicsVector::Energy
G4double Energy(const std::size_t index) const
G4PhysicsVector::Value
G4double Value(const G4double energy, std::size_t &lastidx) const
G4PhysicsVector::GetVectorLength
std::size_t GetVectorLength() const
G4PhysicsVector::FillSecondDerivatives
void FillSecondDerivatives(const G4SplineType=G4SplineType::Base, const G4double dir1=0.0, const G4double dir2=0.0)
Definition G4PhysicsVector.cc:228
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition G4VEmModel.hh:734
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:121
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition G4VEmModel.hh:643
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition G4VEmModel.hh:713
G4VEmModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition G4VEmModel.cc:178
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition G4VEmModel.hh:636
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185