| Geant4 Cross Reference |
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
66 const G4int G4PenelopeRayleighModelMI::fMaxZ;
67 G4PhysicsFreeVector* G4PenelopeRayleighModelMI::fLogAtomicCrossSection[] = {nullptr};
68 G4PhysicsFreeVector* G4PenelopeRayleighModelMI::fAtomicFormFactor[] = {nullptr};
69
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72 G4PenelopeRayleighModelMI::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
113 G4PenelopeRayleighModelMI::~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
152 void 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
180 void 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
280 void 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
318 namespace {G4Mutex PenelopeRayleighModelMutex = G4MUTEX_INITIALIZER;}
319 G4double 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
375 void 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
390 G4double 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
409 G4double 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 (*StoichiometricFactors)[i] /= MaxStoichiometricFactor;
497 }
498
499 //Equivalent atoms per molecule
500 G4double atPerMol = 0.;
501 for (std::size_t i=0;i<nElements;++i)
502 atPerMol += (*StoichiometricFactors)[i];
503 G4double moleculeDensity = 0.;
504 if (atPerMol != 0.) moleculeDensity = atomDensity/atPerMol;
505
506 if (fVerboseLevel > 2)
507 G4cout << "Material " << material->GetName() << " has " << atPerMol << " atoms "
508 << "per molecule and " << moleculeDensity/(cm*cm*cm) << " molecule/cm3" << G4endl;
509
510 //Equivalent molecular weight (dimensionless)
511 G4double MolWeight = 0.;
512 for (std::size_t i=0;i<nElements;++i)
513 MolWeight += (*StoichiometricFactors)[i]*(*elementVector)[i]->GetA()/(g/mole);
514
515 if (fVerboseLevel > 2)
516 G4cout << "Molecular weight of " << matname << ": " << MolWeight << " g/mol" << G4endl;
517
518 G4double IntegrandFun[fNtheta];
519 for (G4int k=0; k<fNtheta; k++) {
520 G4double theta = fAngularFunction->Energy(k); //the x-value is called "Energy", but is an angle
521 G4double F2 = GetFSquared(material,CalculateQSquared(theta,energy));
522 IntegrandFun[k] = (*fAngularFunction)[k]*F2;
523 }
524
525 G4double constant = pi*classic_electr_radius*classic_electr_radius;
526 cs = constant*IntegrateFun(IntegrandFun,fNtheta,fDTheta);
527
528 //Now cs is the cross section per molecule, let's calculate the cross section per volume
529 G4double csvolume = cs*moleculeDensity;
530
531 //print CS and mfp
532 if (fVerboseLevel > 2)
533 G4cout << "Rayleigh CS of " << matname << " at " << energy/keV
534 << " keV: " << cs/barn << " barn"
535 << ", mean free path: " << 1./csvolume/mm << " mm" << G4endl;
536
537 delete StoichiometricFactors;
538 //return CS
539 return csvolume;
540 }
541
542 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
543
544 void G4PenelopeRayleighModelMI::BuildFormFactorTable(const G4Material* material)
545 {
546 if (fVerboseLevel > 3)
547 G4cout << "Calling BuildFormFactorTable() of G4PenelopeRayleighModelMI" << G4endl;
548
549 //GET MATERIAL INFORMATION
550 std::size_t nElements = material->GetNumberOfElements();
551 const G4ElementVector* elementVector = material->GetElementVector();
552 const G4double* fractionVector = material->GetFractionVector();
553
554 //Stoichiometric factors
555 std::vector<G4double> *StoichiometricFactors = new std::vector<G4double>;
556 for (std::size_t i=0;i<nElements;++i) {
557 G4double fraction = fractionVector[i];
558 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
559 StoichiometricFactors->push_back(fraction/atomicWeigth);
560 }
561 G4double MaxStoichiometricFactor = 0.;
562 for (std::size_t i=0;i<nElements;++i) {
563 if ((*StoichiometricFactors)[i] > MaxStoichiometricFactor)
564 MaxStoichiometricFactor = (*StoichiometricFactors)[i];
565 }
566 if (MaxStoichiometricFactor<1e-16) {
567 G4ExceptionDescription ed;
568 ed << "Inconsistent data of atomic composition for " << material->GetName() << G4endl;
569 G4Exception("G4PenelopeRayleighModelMI::BuildFormFactorTable()",
570 "em2042",FatalException,ed);
571 }
572 for (std::size_t i=0;i<nElements;++i)
573 (*StoichiometricFactors)[i] /= MaxStoichiometricFactor;
574
575 //Equivalent molecular weight (dimensionless)
576 G4double MolWeight = 0.;
577 for (std::size_t i=0;i<nElements;++i)
578 MolWeight += (*StoichiometricFactors)[i]*(*elementVector)[i]->GetA()/(g/mole);
579
580 //CREATE THE FORM FACTOR TABLE
581 //First, the form factors are retrieved [F/sqrt(W)].
582 //Then, they are squared and multiplied for MolWeight -> F2 [dimensionless].
583 //This makes difference for CS calculation, but not for theta sampling.
584 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(fLogQSquareGrid.size(),
585 /*spline=*/true);
586
587 G4String matname = material->GetName();
588 G4String aFileNameFF = "";
589
590 //retrieve MIdata (fFileNameFF)
591 G4MIData* dataMI = GetMIData(material);
592 if (dataMI)
593 aFileNameFF = dataMI->GetFilenameFF();
594
595 //read the MIFF from a file passed by the user
596 if (fIsMIActive && aFileNameFF != "") {
597 if (fVerboseLevel)
598 G4cout << "Read MIFF for " << matname << " from custom file: " << aFileNameFF << G4endl;
599
600 ReadMolInterferenceData(matname,aFileNameFF);
601 G4PhysicsFreeVector* Fvector = fMolInterferenceData->find(matname)->second;
602
603 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
604 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
605 G4double f = Fvector->Value(q);
606 G4double ff2 = f*f*MolWeight;
607 if (ff2)
608 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
609 }
610 }
611 //retrieve the MIFF from the database or use the IAM
612 else {
613 //medical material: composition of fat, water, bonematrix, mineral
614 if (fIsMIActive && matname.find("MedMat") != std::string::npos) {
615 if (fVerboseLevel)
616 G4cout << "Build MIFF from components for " << matname << G4endl;
617
618 //get the material composition from its name
619 G4int ki, kf=6, ktot=19;
620 G4double comp[4];
621 G4String compstring = matname.substr(kf+1, ktot);
622 for (std::size_t j=0; j<4; ++j) {
623 ki = kf+1;
624 kf = ki+4;
625 compstring = matname.substr(ki, 4);
626 comp[j] = atof(compstring.c_str());
627 if (fVerboseLevel > 2)
628 G4cout << " -- MedMat comp[" << j+1 << "]: " << comp[j] << G4endl;
629 }
630
631 const char* path = G4FindDataDir("G4LEDATA");
632 if (!path) {
633 G4String excep = "G4LEDATA environment variable not set!";
634 G4Exception("G4PenelopeRayleighModelMI::BuildFormFactorTable()",
635 "em0006",FatalException,excep);
636 }
637
638 if (!fMolInterferenceData->count("Fat_MI"))
639 ReadMolInterferenceData("Fat_MI");
640 G4PhysicsFreeVector* fatFF = fMolInterferenceData->find("Fat_MI")->second;
641
642 if (!fMolInterferenceData->count("Water_MI"))
643 ReadMolInterferenceData("Water_MI");
644 G4PhysicsFreeVector* waterFF = fMolInterferenceData->find("Water_MI")->second;
645
646 if (!fMolInterferenceData->count("BoneMatrix_MI"))
647 ReadMolInterferenceData("BoneMatrix_MI");
648 G4PhysicsFreeVector* bonematrixFF = fMolInterferenceData->find("BoneMatrix_MI")->second;
649
650 if (!fMolInterferenceData->count("Mineral_MI"))
651 ReadMolInterferenceData("Mineral_MI");
652 G4PhysicsFreeVector* mineralFF = fMolInterferenceData->find("Mineral_MI")->second;
653
654 //get and combine the molecular form factors with interference effect
655 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
656 G4double ff2 = 0;
657 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
658 G4double ffat = fatFF->Value(q);
659 G4double fwater = waterFF->Value(q);
660 G4double fbonematrix = bonematrixFF->Value(q);
661 G4double fmineral = mineralFF->Value(q);
662 ff2 = comp[0]*ffat*ffat+comp[1]*fwater*fwater+
663 comp[2]*fbonematrix*fbonematrix+comp[3]*fmineral*fmineral;
664 ff2 *= MolWeight;
665 if (ff2) theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
666 }
667 }
668 //other materials with MIFF (from the database)
669 else if (fIsMIActive && fMolInterferenceData->count(matname)) {
670 if (fVerboseLevel)
671 G4cout << "Read MIFF from database " << matname << G4endl;
672 G4PhysicsFreeVector* FF = fMolInterferenceData->find(matname)->second;
673 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
674 G4double ff2 = 0;
675 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
676 G4double f = FF->Value(q);
677 ff2 = f*f*MolWeight;
678 if (ff2) theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
679 }
680 }
681 //IAM
682 else {
683 if (fVerboseLevel)
684 G4cout << "FF of " << matname << " calculated according to the IAM" << G4endl;
685 for (std::size_t k=0;k<fLogQSquareGrid.size();++k) {
686 G4double ff2 = 0;
687 for (std::size_t i=0;i<nElements;++i) {
688 G4int iZ = (*elementVector)[i]->GetZasInt();
689 G4PhysicsFreeVector* theAtomVec = fAtomicFormFactor[iZ];
690 G4double q = std::pow(G4Exp(fLogQSquareGrid[k]),0.5);
691 G4double f = theAtomVec->Value(q);
692 ff2 += f*f*(*StoichiometricFactors)[i];
693 }
694 if (ff2) theFFVec->PutValue(k,fLogQSquareGrid[k],G4Log(ff2));
695 }
696 }
697 }
698 theFFVec->FillSecondDerivatives();
699 fLogFormFactorTable->insert(std::make_pair(material,theFFVec));
700
701 if (fVerboseLevel > 3)
702 DumpFormFactorTable(material);
703 delete StoichiometricFactors;
704
705 return;
706 }
707
708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
709
710 void G4PenelopeRayleighModelMI::SampleSecondaries(std::vector<G4DynamicParticle*>*,
711 const G4MaterialCutsCouple* couple,
712 const G4DynamicParticle* aDynamicGamma,
713 G4double,
714 G4double)
715 {
716 // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
717 // from the Penelope2008 model. The scattering angle is sampled from the atomic
718 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
719 // anomalous scattering effects. The Form Factor F(Q) function which appears in the
720 // analytical cross section is retrieved via the method GetFSquared(); atomic data
721 // are tabulated for F(Q). Form factor for compounds is calculated according to
722 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
723 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
724 // for each material and managed by G4PenelopeSamplingData objects.
725 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
726 // increases with energy. For E=100 keV the efficiency is 100% and 86% for
727 // hydrogen and uranium, respectively.
728
729 if (fVerboseLevel > 3)
730 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModelMI" << G4endl;
731
732 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
733
734 if (photonEnergy0 <= fIntrinsicLowEnergyLimit) {
735 fParticleChange->ProposeTrackStatus(fStopAndKill);
736 fParticleChange->SetProposedKineticEnergy(0.);
737 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
738 return;
739 }
740
741 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
742 const G4Material* theMat = couple->GetMaterial();
743
744 //1) Verify if tables are ready
745 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
746 //not invoked
747 if (!fPMaxTable || !fSamplingTable || !fLogFormFactorTable) {
748 //create a **thread-local** version of the table. Used only for G4EmCalculator and
749 //Unit Tests
750 fLocalTable = true;
751 if (!fLogFormFactorTable)
752 fLogFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
753 if (!fPMaxTable)
754 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
755 if (!fSamplingTable)
756 fSamplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
757 if (fIsMIActive && !fMolInterferenceData)
758 fMolInterferenceData = new std::map<G4String,G4PhysicsFreeVector*>;
759 }
760
761 if (!fSamplingTable->count(theMat)) {
762 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
763 //not filled up. This can happen in a UnitTest
764 if (fVerboseLevel > 0) {
765 //Issue a G4Exception (warning) only in verbose mode
766 G4ExceptionDescription ed;
767 ed << "Unable to find the fSamplingTable data for " <<
768 theMat->GetName() << G4endl;
769 ed << "This can happen only in Unit Tests" << G4endl;
770 G4Exception("G4PenelopeRayleighModelMI::SampleSecondaries()",
771 "em2019",JustWarning,ed);
772 }
773 const G4ElementVector* theElementVector = theMat->GetElementVector();
774 //protect file reading via autolock
775 G4AutoLock lock(&PenelopeRayleighModelMutex);
776 for (std::size_t j=0;j<theMat->GetNumberOfElements();++j) {
777 G4int iZ = theElementVector->at(j)->GetZasInt();
778 if (!fLogAtomicCrossSection[iZ]) {
779 lock.lock();
780 ReadDataFile(iZ);
781 lock.unlock();
782 }
783 }
784 lock.lock();
785
786 //1) If the table has not been built for the material, do it!
787 if (!fLogFormFactorTable->count(theMat))
788 BuildFormFactorTable(theMat);
789
790 //2) retrieve or build the sampling table
791 if (!(fSamplingTable->count(theMat)))
792 InitializeSamplingAlgorithm(theMat);
793
794 //3) retrieve or build the pMax data
795 if (!fPMaxTable->count(theMat))
796 GetPMaxTable(theMat);
797
798 lock.unlock();
799 }
800
801 //Ok, restart the job
802 G4PenelopeSamplingData* theDataTable = fSamplingTable->find(theMat)->second;
803 G4PhysicsFreeVector* thePMax = fPMaxTable->find(theMat)->second;
804 G4double cosTheta = 1.0;
805
806 //OK, ready to go!
807 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
808
809 if (qmax < 1e-10) { //very low momentum transfer
810 G4bool loopAgain=false;
811 do {
812 loopAgain = false;
813 cosTheta = 1.0-2.0*G4UniformRand();
814 G4double G = 0.5*(1+cosTheta*cosTheta);
815 if (G4UniformRand()>G)
816 loopAgain = true;
817 } while(loopAgain);
818 }
819 else { //larger momentum transfer
820 std::size_t nData = theDataTable->GetNumberOfStoredPoints();
821 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
822 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
823
824 G4bool loopAgain = false;
825 G4double MaxPValue = thePMax->Value(photonEnergy0);
826 G4double xx=0;
827
828 //Sampling by rejection method. The rejection function is
829 //G = 0.5*(1+cos^2(theta))
830 do {
831 loopAgain = false;
832 G4double RandomMax = G4UniformRand()*MaxPValue;
833 xx = theDataTable->SampleValue(RandomMax);
834
835 //xx is a random value of q^2 in (0,q2max),sampled according to
836 //F(Q^2) via the RITA algorithm
837 if (xx > q2max)
838 loopAgain = true;
839 cosTheta = 1.0-2.0*xx/q2max;
840 G4double G = 0.5*(1+cosTheta*cosTheta);
841 if (G4UniformRand()>G)
842 loopAgain = true;
843 } while(loopAgain);
844 }
845
846 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
847
848 //Scattered photon angles. ( Z - axis along the parent photon)
849 G4double phi = twopi * G4UniformRand() ;
850 G4double dirX = sinTheta*std::cos(phi);
851 G4double dirY = sinTheta*std::sin(phi);
852 G4double dirZ = cosTheta;
853
854 //Update G4VParticleChange for the scattered photon
855 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
856
857 photonDirection1.rotateUz(photonDirection0);
858 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
859 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
860
861 return;
862 }
863
864 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
865
866 void G4PenelopeRayleighModelMI::ReadDataFile(const G4int Z)
867 {
868 if (fVerboseLevel > 2) {
869 G4cout << "G4PenelopeRayleighModelMI::ReadDataFile()" << G4endl;
870 G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
871 }
872
873 const char* path = G4FindDataDir("G4LEDATA");
874 if (!path) {
875 G4String excep = "G4LEDATA environment variable not set!";
876 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
877 "em0006",FatalException,excep);
878 return;
879 }
880
881 //Read first the cross section file (all the files have 250 points)
882 std::ostringstream ost;
883 if (Z>9)
884 ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
885 else
886 ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
887 std::ifstream file(ost.str().c_str());
888
889 if (!file.is_open()) {
890 G4String excep = "Data file " + G4String(ost.str()) + " not found!";
891 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
892 "em0003",FatalException,excep);
893 }
894
895 G4int readZ = 0;
896 std::size_t nPoints = 0;
897 file >> readZ >> nPoints;
898
899 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) {
900 G4ExceptionDescription ed;
901 ed << "Corrupted data file for Z=" << Z << G4endl;
902 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
903 "em0005",FatalException,ed);
904 return;
905 }
906
907 fLogAtomicCrossSection[Z] = new G4PhysicsFreeVector((std::size_t)nPoints);
908 G4double ene=0,f1=0,f2=0,xs=0;
909 for (std::size_t i=0;i<nPoints;++i) {
910 file >> ene >> f1 >> f2 >> xs;
911 //dimensional quantities
912 ene *= eV;
913 xs *= cm2;
914 fLogAtomicCrossSection[Z]->PutValue(i,G4Log(ene),G4Log(xs));
915 if (file.eof() && i != (nPoints-1)) { //file ended too early
916 G4ExceptionDescription ed ;
917 ed << "Corrupted data file for Z=" << Z << G4endl;
918 ed << "Found less than " << nPoints << " entries" << G4endl;
919 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
920 "em0005",FatalException,ed);
921 }
922 }
923 file.close();
924
925 //Then, read the extended momentum transfer file
926 std::ostringstream ost2;
927 ost2 << path << "/penelope/rayleigh/MIFF/qext.dat";
928 file.open(ost2.str().c_str());
929
930 if (!file.is_open()) {
931 G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
932 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
933 "em0003",FatalException,excep);
934 }
935 G4bool fillQGrid = false;
936 if (!fLogQSquareGrid.size()) {
937 fillQGrid = true;
938 nPoints = 1142;
939 }
940 G4double qext = 0;
941 if (fillQGrid) { //fLogQSquareGrid filled only one time
942 for (std::size_t i=0;i<nPoints;++i) {
943 file >> qext;
944 fLogQSquareGrid.push_back(2.0*G4Log(qext));
945 }
946 }
947 file.close();
948
949 //Finally, read the form factor file
950 std::ostringstream ost3;
951 if (Z>9)
952 ost3 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
953 else
954 ost3 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
955 file.open(ost3.str().c_str());
956
957 if (!file.is_open()) {
958 G4String excep = "Data file " + G4String(ost3.str()) + " not found!";
959 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
960 "em0003",FatalException,excep);
961 }
962
963 file >> readZ >> nPoints;
964
965 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) {
966 G4ExceptionDescription ed;
967 ed << "Corrupted data file for Z=" << Z << G4endl;
968 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
969 "em0005",FatalException,ed);
970 return;
971 }
972
973 fAtomicFormFactor[Z] = new G4PhysicsFreeVector((std::size_t)nPoints);
974 G4double q=0,ff=0,incoh=0;
975 for (std::size_t i=0;i<nPoints;++i) {
976 file >> q >> ff >> incoh; //q and ff are dimensionless (q is in units of (m_e*c))
977 fAtomicFormFactor[Z]->PutValue(i,q,ff);
978 if (file.eof() && i != (nPoints-1)) { //file ended too early
979 G4ExceptionDescription ed;
980 ed << "Corrupted data file for Z=" << Z << G4endl;
981 ed << "Found less than " << nPoints << " entries" << G4endl;
982 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
983 "em0005",FatalException,ed);
984 }
985 }
986 file.close();
987 return;
988 }
989
990 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
991
992 void G4PenelopeRayleighModelMI::ReadMolInterferenceData(const G4String& matname,
993 const G4String& FFfilename)
994 {
995 if (fVerboseLevel > 2) {
996 G4cout << "G4PenelopeRayleighModelMI::ReadMolInterferenceData() for material " <<
997 matname << G4endl;
998 }
999 G4bool isLocalFile = (FFfilename != "NULL");
1000
1001 const char* path = G4FindDataDir("G4LEDATA");
1002 if (!path) {
1003 G4String excep = "G4LEDATA environment variable not set!";
1004 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1005 "em0006",FatalException,excep);
1006 return;
1007 }
1008
1009 if (!(fKnownMaterials->count(matname)) && !isLocalFile) //material not found
1010 return;
1011
1012 G4String aFileName = (isLocalFile) ? FFfilename : fKnownMaterials->find(matname)->second;
1013
1014 //if the material has a MIFF, read it from the database
1015 if (aFileName != "") {
1016 if (fVerboseLevel > 2)
1017 G4cout << "ReadMolInterferenceData(). Read material: " << matname << ", filename: " <<
1018 aFileName << " " <<
1019 (isLocalFile ? "(local)" : "(database)") << G4endl;
1020 std::ifstream file;
1021 std::ostringstream ostIMFF;
1022 if (isLocalFile)
1023 ostIMFF << aFileName;
1024 else
1025 ostIMFF << path << "/penelope/rayleigh/MIFF/" << aFileName;
1026 file.open(ostIMFF.str().c_str());
1027
1028 if (!file.is_open()) {
1029 G4String excep = "Data file " + G4String(ostIMFF.str()) + " not found!";
1030 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1031 "em1031",FatalException,excep);
1032 return;
1033 }
1034
1035 //check the number of points in the file
1036 std::size_t nPoints = 0;
1037 G4double x=0,y=0;
1038 while (file.good()) {
1039 file >> x >> y;
1040 nPoints++;
1041 }
1042 file.close();
1043 nPoints--;
1044 if (fVerboseLevel > 3)
1045 G4cout << "Number of nPoints: " << nPoints << G4endl;
1046
1047 //read the file
1048 file.open(ostIMFF.str().c_str());
1049 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((std::size_t)nPoints);
1050 G4double q=0,ff=0;
1051 for (std::size_t i=0; i<nPoints; ++i) {
1052 file >> q >> ff;
1053
1054 //q and ff are dimensionless (q is in units of (m_e*c))
1055 theFFVec->PutValue(i,q,ff);
1056
1057 //file ended too early
1058 if (file.eof() && i != (nPoints-1)) {
1059 G4ExceptionDescription ed;
1060 ed << "Corrupted data file" << G4endl;
1061 ed << "Found less than " << nPoints << " entries" << G4endl;
1062 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1063 "em1005",FatalException,ed);
1064 }
1065 }
1066 if (!fMolInterferenceData) {
1067 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1068 "em2145",FatalException,
1069 "Unable to allocate the Molecular Interference data table");
1070 delete theFFVec;
1071 return;
1072 }
1073 file.close();
1074 fMolInterferenceData->insert(std::make_pair(matname,theFFVec));
1075 }
1076 return;
1077 }
1078
1079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1080
1081 G4double G4PenelopeRayleighModelMI::GetFSquared(const G4Material* mat, const G4double QSquared)
1082 {
1083 G4double f2 = 0;
1084 //Input value QSquared could be zero: protect the log() below against
1085 //the FPE exception
1086
1087 //If Q<1e-10, set Q to 1e-10
1088 G4double logQSquared = (QSquared>1e-10) ? G4Log(QSquared) : -23.;
1089 //last value of the table
1090 G4double maxlogQ2 = fLogQSquareGrid[fLogQSquareGrid.size()-1];
1091
1092 //now it should be all right
1093 G4PhysicsFreeVector* theVec = fLogFormFactorTable->find(mat)->second;
1094
1095 if (!theVec) {
1096 G4ExceptionDescription ed;
1097 ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
1098 G4Exception("G4PenelopeRayleighModelMI::GetFSquared()",
1099 "em2046",FatalException,ed);
1100 return 0;
1101 }
1102
1103 if (logQSquared < -20) { //Q < 1e-9
1104 G4double logf2 = (*theVec)[0]; //first value of the table
1105 f2 = G4Exp(logf2);
1106 }
1107 else if (logQSquared > maxlogQ2)
1108 f2 =0;
1109 else {
1110 //log(Q^2) vs. log(F^2)
1111 G4double logf2 = theVec->Value(logQSquared);
1112 f2 = G4Exp(logf2);
1113 }
1114
1115 if (fVerboseLevel > 3) {
1116 G4cout << "G4PenelopeRayleighModelMI::GetFSquared() in " << mat->GetName() << G4endl;
1117 G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c)); F^2 = " << f2 << G4endl;
1118 }
1119 return f2;
1120 }
1121
1122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1123
1124 void G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm(const G4Material* mat)
1125 {
1126 G4double q2min = 0;
1127 G4double q2max = 0;
1128 const std::size_t np = 150; //hard-coded in Penelope
1129 for (std::size_t i=1;i<fLogQSquareGrid.size();++i)
1130 {
1131 G4double Q2 = G4Exp(fLogQSquareGrid[i]);
1132 if (GetFSquared(mat,Q2) > 1e-35)
1133 {
1134 q2max = G4Exp(fLogQSquareGrid[i-1]);
1135 }
1136 }
1137 std::size_t nReducedPoints = np/4;
1138
1139 //check for errors
1140 if (np < 16)
1141 {
1142 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1143 "em2047",FatalException,
1144 "Too few points to initialize the sampling algorithm");
1145 }
1146 if (q2min > (q2max-1e-10))
1147 {
1148 G4cout << "q2min= " << q2min << " q2max= " << q2max << G4endl;
1149 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1150 "em2048",FatalException,
1151 "Too narrow grid to initialize the sampling algorithm");
1152 }
1153
1154 //This is subroutine RITAI0 of Penelope
1155 //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material*
1156
1157 //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
1158 G4DataVector* x = new G4DataVector();
1159
1160 /*******************************************************************************
1161 Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
1162 ********************************************************************************/
1163 std::size_t NUNIF = std::min(std::max(((std::size_t)8),nReducedPoints),np/2);
1164 const G4int nip = 51; //hard-coded in Penelope
1165
1166 G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
1167 x->push_back(q2min);
1168 for (std::size_t i=1;i<NUNIF-1;++i)
1169 {
1170 G4double app = q2min + i*dx;
1171 x->push_back(app); //increase
1172 }
1173 x->push_back(q2max);
1174
1175 if (fVerboseLevel> 3)
1176 G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
1177
1178 //Allocate temporary storage vectors
1179 G4DataVector* area = new G4DataVector();
1180 G4DataVector* a = new G4DataVector();
1181 G4DataVector* b = new G4DataVector();
1182 G4DataVector* c = new G4DataVector();
1183 G4DataVector* err = new G4DataVector();
1184
1185 for (std::size_t i=0;i<NUNIF-1;++i) //build all points but the last
1186 {
1187 //Temporary vectors for this loop
1188 G4DataVector* pdfi = new G4DataVector();
1189 G4DataVector* pdfih = new G4DataVector();
1190 G4DataVector* sumi = new G4DataVector();
1191
1192 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
1193 G4double pdfmax = 0;
1194 for (G4int k=0;k<nip;k++)
1195 {
1196 G4double xik = (*x)[i]+k*dxi;
1197 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1198 pdfi->push_back(pdfk);
1199 pdfmax = std::max(pdfmax,pdfk);
1200 if (k < (nip-1))
1201 {
1202 G4double xih = xik + 0.5*dxi;
1203 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1204 pdfih->push_back(pdfIK);
1205 pdfmax = std::max(pdfmax,pdfIK);
1206 }
1207 }
1208
1209 //Simpson's integration
1210 G4double cons = dxi*0.5*(1./3.);
1211 sumi->push_back(0.);
1212 for (G4int k=1;k<nip;k++)
1213 {
1214 G4double previous = (*sumi)[k-1];
1215 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1216 sumi->push_back(next);
1217 }
1218
1219 G4double lastIntegral = (*sumi)[sumi->size()-1];
1220 area->push_back(lastIntegral);
1221 //Normalize cumulative function
1222 G4double factor = 1.0/lastIntegral;
1223 for (std::size_t k=0;k<sumi->size();++k)
1224 (*sumi)[k] *= factor;
1225
1226 //When the PDF vanishes at one of the interval end points, its value is modified
1227 if ((*pdfi)[0] < 1e-35)
1228 (*pdfi)[0] = 1e-5*pdfmax;
1229 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1230 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1231
1232 G4double pli = (*pdfi)[0]*factor;
1233 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1234 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
1235 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
1236 G4double C_temp = 1.0+A_temp+B_temp;
1237 if (C_temp < 1e-35)
1238 {
1239 a->push_back(0.);
1240 b->push_back(0.);
1241 c->push_back(1.);
1242 }
1243 else
1244 {
1245 a->push_back(A_temp);
1246 b->push_back(B_temp);
1247 c->push_back(C_temp);
1248 }
1249
1250 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1251 //and the true pdf, extended over the interval (X(I),X(I+1))
1252 G4int icase = 1; //loop code
1253 G4bool reLoop = false;
1254 err->push_back(0.);
1255 do
1256 {
1257 reLoop = false;
1258 (*err)[i] = 0.; //zero variable
1259 for (G4int k=0;k<nip;k++)
1260 {
1261 G4double rr = (*sumi)[k];
1262 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1263 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1264 if (k == 0 || k == nip-1)
1265 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1266 else
1267 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1268 }
1269 (*err)[i] *= dxi;
1270
1271 //If err(I) is too large, the pdf is approximated by a uniform distribution
1272 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1273 {
1274 (*b)[i] = 0;
1275 (*a)[i] = 0;
1276 (*c)[i] = 1.;
1277 icase = 2;
1278 reLoop = true;
1279 }
1280 }while(reLoop);
1281
1282 delete pdfi;
1283 delete pdfih;
1284 delete sumi;
1285 } //end of first loop over i
1286
1287 //Now assign last point
1288 (*x)[x->size()-1] = q2max;
1289 a->push_back(0.);
1290 b->push_back(0.);
1291 c->push_back(0.);
1292 err->push_back(0.);
1293 area->push_back(0.);
1294
1295 if (x->size() != NUNIF || a->size() != NUNIF ||
1296 err->size() != NUNIF || area->size() != NUNIF)
1297 {
1298 G4ExceptionDescription ed;
1299 ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
1300 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1301 "em2049",FatalException,ed);
1302 }
1303
1304 /*******************************************************************************
1305 New grid points are added by halving the sub-intervals with the largest absolute error
1306 This is done up to np=150 points in the grid
1307 ********************************************************************************/
1308 do
1309 {
1310 G4double maxError = 0.0;
1311 std::size_t iErrMax = 0;
1312 for (std::size_t i=0;i<err->size()-2;++i)
1313 {
1314 //maxError is the lagest of the interval errors err[i]
1315 if ((*err)[i] > maxError)
1316 {
1317 maxError = (*err)[i];
1318 iErrMax = i;
1319 }
1320 }
1321
1322 //OK, now I have to insert one new point in the position iErrMax
1323 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
1324
1325 x->insert(x->begin()+iErrMax+1,newx);
1326 //Add place-holders in the other vectors
1327 area->insert(area->begin()+iErrMax+1,0.);
1328 a->insert(a->begin()+iErrMax+1,0.);
1329 b->insert(b->begin()+iErrMax+1,0.);
1330 c->insert(c->begin()+iErrMax+1,0.);
1331 err->insert(err->begin()+iErrMax+1,0.);
1332
1333 //Now calculate the other parameters
1334 for (std::size_t i=iErrMax;i<=iErrMax+1;++i)
1335 {
1336 //Temporary vectors for this loop
1337 G4DataVector* pdfi = new G4DataVector();
1338 G4DataVector* pdfih = new G4DataVector();
1339 G4DataVector* sumi = new G4DataVector();
1340
1341 G4double dxLocal = (*x)[i+1]-(*x)[i];
1342 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
1343 G4double pdfmax = 0;
1344 for (G4int k=0;k<nip;k++)
1345 {
1346 G4double xik = (*x)[i]+k*dxi;
1347 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1348 pdfi->push_back(pdfk);
1349 pdfmax = std::max(pdfmax,pdfk);
1350 if (k < (nip-1))
1351 {
1352 G4double xih = xik + 0.5*dxi;
1353 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1354 pdfih->push_back(pdfIK);
1355 pdfmax = std::max(pdfmax,pdfIK);
1356 }
1357 }
1358
1359 //Simpson's integration
1360 G4double cons = dxi*0.5*(1./3.);
1361 sumi->push_back(0.);
1362 for (G4int k=1;k<nip;k++)
1363 {
1364 G4double previous = (*sumi)[k-1];
1365 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1366 sumi->push_back(next);
1367 }
1368 G4double lastIntegral = (*sumi)[sumi->size()-1];
1369 (*area)[i] = lastIntegral;
1370
1371 //Normalize cumulative function
1372 G4double factor = 1.0/lastIntegral;
1373 for (std::size_t k=0;k<sumi->size();++k)
1374 (*sumi)[k] *= factor;
1375
1376 //When the PDF vanishes at one of the interval end points, its value is modified
1377 if ((*pdfi)[0] < 1e-35)
1378 (*pdfi)[0] = 1e-5*pdfmax;
1379 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1380 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1381
1382 G4double pli = (*pdfi)[0]*factor;
1383 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1384 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
1385 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
1386 G4double C_temp = 1.0+A_temp+B_temp;
1387 if (C_temp < 1e-35)
1388 {
1389 (*a)[i]= 0.;
1390 (*b)[i] = 0.;
1391 (*c)[i] = 1;
1392 }
1393 else
1394 {
1395 (*a)[i]= A_temp;
1396 (*b)[i] = B_temp;
1397 (*c)[i] = C_temp;
1398 }
1399 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1400 //and the true pdf, extended over the interval (X(I),X(I+1))
1401 G4int icase = 1; //loop code
1402 G4bool reLoop = false;
1403 do
1404 {
1405 reLoop = false;
1406 (*err)[i] = 0.; //zero variable
1407 for (G4int k=0;k<nip;k++)
1408 {
1409 G4double rr = (*sumi)[k];
1410 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1411 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1412 if (k == 0 || k == nip-1)
1413 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1414 else
1415 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1416 }
1417 (*err)[i] *= dxi;
1418
1419 //If err(I) is too large, the pdf is approximated by a uniform distribution
1420 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1421 {
1422 (*b)[i] = 0;
1423 (*a)[i] = 0;
1424 (*c)[i] = 1.;
1425 icase = 2;
1426 reLoop = true;
1427 }
1428 }while(reLoop);
1429 delete pdfi;
1430 delete pdfih;
1431 delete sumi;
1432 }
1433 }while(x->size() < np);
1434
1435 if (x->size() != np || a->size() != np ||
1436 err->size() != np || area->size() != np)
1437 {
1438 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1439 "em2050",FatalException,
1440 "Problem in building the extended Table for Sampling: array dimensions do not match ");
1441 }
1442
1443 /*******************************************************************************
1444 Renormalization
1445 ********************************************************************************/
1446 G4double ws = 0;
1447 for (std::size_t i=0;i<np-1;++i)
1448 ws += (*area)[i];
1449 ws = 1.0/ws;
1450 G4double errMax = 0;
1451 for (std::size_t i=0;i<np-1;++i)
1452 {
1453 (*area)[i] *= ws;
1454 (*err)[i] *= ws;
1455 errMax = std::max(errMax,(*err)[i]);
1456 }
1457
1458 //Vector with the normalized cumulative distribution
1459 G4DataVector* PAC = new G4DataVector();
1460 PAC->push_back(0.);
1461 for (std::size_t i=0;i<np-1;++i)
1462 {
1463 G4double previous = (*PAC)[i];
1464 PAC->push_back(previous+(*area)[i]);
1465 }
1466 (*PAC)[PAC->size()-1] = 1.;
1467
1468 /*******************************************************************************
1469 Pre-calculated limits for the initial binary search for subsequent sampling
1470 ********************************************************************************/
1471 std::vector<std::size_t> *ITTL = new std::vector<std::size_t>;
1472 std::vector<std::size_t> *ITTU = new std::vector<std::size_t>;
1473
1474 //Just create place-holders
1475 for (std::size_t i=0;i<np;++i)
1476 {
1477 ITTL->push_back(0);
1478 ITTU->push_back(0);
1479 }
1480
1481 G4double bin = 1.0/(np-1);
1482 (*ITTL)[0]=0;
1483 for (std::size_t i=1;i<(np-1);++i)
1484 {
1485 G4double ptst = i*bin;
1486 G4bool found = false;
1487 for (std::size_t j=(*ITTL)[i-1];j<np && !found;++j)
1488 {
1489 if ((*PAC)[j] > ptst)
1490 {
1491 (*ITTL)[i] = j-1;
1492 (*ITTU)[i-1] = j;
1493 found = true;
1494 }
1495 }
1496 }
1497 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1498 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1499 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1500
1501 if (ITTU->size() != np || ITTU->size() != np)
1502 {
1503 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1504 "em2051",FatalException,
1505 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1506 }
1507
1508 /********************************************************************************
1509 Copy tables
1510 ********************************************************************************/
1511 G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
1512 for (std::size_t i=0;i<np;++i)
1513 {
1514 theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1515 }
1516 if (fVerboseLevel > 2)
1517 {
1518 G4cout << "*************************************************************************" <<
1519 G4endl;
1520 G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
1521 theTable->DumpTable();
1522 }
1523 fSamplingTable->insert(std::make_pair(mat,theTable));
1524
1525 //Clean up temporary vectors
1526 delete x;
1527 delete a;
1528 delete b;
1529 delete c;
1530 delete err;
1531 delete area;
1532 delete PAC;
1533 delete ITTL;
1534 delete ITTU;
1535
1536 //DONE!
1537 return;
1538 }
1539
1540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1541
1542 void G4PenelopeRayleighModelMI::GetPMaxTable(const G4Material* mat)
1543 {
1544 if (!fPMaxTable)
1545 {
1546 G4cout << "G4PenelopeRayleighModelMI::BuildPMaxTable" << G4endl;
1547 G4cout << "Going to instanziate the fPMaxTable !" << G4endl;
1548 G4cout << "That should _not_ be here! " << G4endl;
1549 fPMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
1550 }
1551 //check if the table is already there
1552 if (fPMaxTable->count(mat))
1553 return;
1554
1555 //otherwise build it
1556 if (!fSamplingTable)
1557 {
1558 G4Exception("G4PenelopeRayleighModelMI::GetPMaxTable()",
1559 "em2052",FatalException,
1560 "SamplingTable is not properly instantiated");
1561 return;
1562 }
1563
1564 //This should not be: the sampling table is built before the p-table
1565 if (!fSamplingTable->count(mat))
1566 {
1567 G4ExceptionDescription ed;
1568 ed << "Sampling table for material " << mat->GetName() << " not found";
1569 G4Exception("G4PenelopeRayleighModelMI::GetPMaxTable()",
1570 "em2052",FatalException,
1571 ed);
1572 return;
1573 }
1574
1575 G4PenelopeSamplingData *theTable = fSamplingTable->find(mat)->second;
1576 std::size_t tablePoints = theTable->GetNumberOfStoredPoints();
1577 std::size_t nOfEnergyPoints = fLogEnergyGridPMax.size();
1578 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
1579
1580 const std::size_t nip = 51; //hard-coded in Penelope
1581
1582 for (std::size_t ie=0;ie<fLogEnergyGridPMax.size();++ie)
1583 {
1584 G4double energy = G4Exp(fLogEnergyGridPMax[ie]);
1585 G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
1586 G4double Qm2 = Qm*Qm;
1587 G4double firstQ2 = theTable->GetX(0);
1588 G4double lastQ2 = theTable->GetX(tablePoints-1);
1589 G4double thePMax = 0;
1590
1591 if (Qm2 > firstQ2)
1592 {
1593 if (Qm2 < lastQ2)
1594 {
1595 //bisection to look for the index of Qm
1596 std::size_t lowerBound = 0;
1597 std::size_t upperBound = tablePoints-1;
1598 while (lowerBound <= upperBound)
1599 {
1600 std::size_t midBin = (lowerBound + upperBound)/2;
1601 if( Qm2 < theTable->GetX(midBin))
1602 { upperBound = midBin-1; }
1603 else
1604 { lowerBound = midBin+1; }
1605 }
1606 //upperBound is the output (but also lowerBounf --> should be the same!)
1607 G4double Q1 = theTable->GetX(upperBound);
1608 G4double Q2 = Qm2;
1609 G4double DQ = (Q2-Q1)/((G4double)(nip-1));
1610 G4double theA = theTable->GetA(upperBound);
1611 G4double theB = theTable->GetB(upperBound);
1612 G4double thePAC = theTable->GetPAC(upperBound);
1613 G4DataVector* fun = new G4DataVector();
1614 for (std::size_t k=0;k<nip;++k)
1615 {
1616 G4double qi = Q1 + k*DQ;
1617 G4double tau = (qi-Q1)/
1618 (theTable->GetX(upperBound+1)-Q1);
1619 G4double con1 = 2.0*theB*tau;
1620 G4double ci = 1.0+theA+theB;
1621 G4double con2 = ci-theA*tau;
1622 G4double etap = 0;
1623 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1624 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1625 else
1626 etap = tau/con2;
1627 G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
1628 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1629 ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
1630 fun->push_back(theFun);
1631 }
1632 //Now intergrate numerically the fun Cavalieri-Simpson's method
1633 G4DataVector* sum = new G4DataVector;
1634 G4double CONS = DQ*(1./12.);
1635 G4double HCONS = 0.5*CONS;
1636 sum->push_back(0.);
1637 G4double secondPoint = (*sum)[0] +
1638 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1639 sum->push_back(secondPoint);
1640 for (std::size_t hh=2;hh<nip-1;++hh)
1641 {
1642 G4double previous = (*sum)[hh-1];
1643 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1644 (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1645 sum->push_back(next);
1646 }
1647 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1648 (*fun)[nip-3])*CONS;
1649 sum->push_back(last);
1650 thePMax = thePAC + (*sum)[sum->size()-1]; //last point
1651 delete fun;
1652 delete sum;
1653 }
1654 else
1655 {
1656 thePMax = 1.0;
1657 }
1658 }
1659 else
1660 {
1661 thePMax = theTable->GetPAC(0);
1662 }
1663
1664 //Write number in the table
1665 theVec->PutValue(ie,energy,thePMax);
1666 }
1667
1668 fPMaxTable->insert(std::make_pair(mat,theVec));
1669 return;
1670 }
1671
1672 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1673
1674 void G4PenelopeRayleighModelMI::DumpFormFactorTable(const G4Material* mat)
1675 {
1676 G4cout << "*****************************************************************" << G4endl;
1677 G4cout << "G4PenelopeRayleighModelMI: Form Factor Table for " << mat->GetName() << G4endl;
1678 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1679 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1680 G4cout << "*****************************************************************" << G4endl;
1681 if (!fLogFormFactorTable->count(mat))
1682 BuildFormFactorTable(mat);
1683
1684 G4PhysicsFreeVector* theVec = fLogFormFactorTable->find(mat)->second;
1685 for (std::size_t i=0;i<theVec->GetVectorLength();++i)
1686 {
1687 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1688 G4double Q = G4Exp(0.5*logQ2);
1689 G4double logF2 = (*theVec)[i];
1690 G4double F = G4Exp(0.5*logF2);
1691 G4cout << Q << " " << F << G4endl;
1692 }
1693 //DONE
1694 return;
1695 }
1696
1697 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1698
1699 void G4PenelopeRayleighModelMI::SetParticle(const G4ParticleDefinition* p)
1700 {
1701 if(!fParticle) {
1702 fParticle = p;
1703 }
1704 }
1705
1706 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1707 void G4PenelopeRayleighModelMI::LoadKnownMIFFMaterials()
1708 {
1709 fKnownMaterials->insert(std::pair<G4String,G4String>("Fat_MI","FF_fat_Tartari2002.dat"));
1710 fKnownMaterials->insert(std::pair<G4String,G4String>("Water_MI","FF_water_Tartari2002.dat"));
1711 fKnownMaterials->insert(std::pair<G4String,G4String>("BoneMatrix_MI","FF_bonematrix_Tartari2002.dat"));
1712 fKnownMaterials->insert(std::pair<G4String,G4String>("Mineral_MI","FF_mineral_Tartari2002.dat"));
1713 fKnownMaterials->insert(std::pair<G4String,G4String>("adipose_MI","FF_adipose_Poletti2002.dat"));
1714 fKnownMaterials->insert(std::pair<G4String,G4String>("glandular_MI","FF_glandular_Poletti2002.dat"));
1715 fKnownMaterials->insert(std::pair<G4String,G4String>("breast5050_MI","FF_human_breast_Peplow1998.dat"));
1716 fKnownMaterials->insert(std::pair<G4String,G4String>("carcinoma_MI","FF_carcinoma_Kidane1999.dat"));
1717 fKnownMaterials->insert(std::pair<G4String,G4String>("muscle_MI","FF_pork_muscle_Peplow1998.dat"));
1718 fKnownMaterials->insert(std::pair<G4String,G4String>("kidney_MI","FF_pork_kidney_Peplow1998.dat"));
1719 fKnownMaterials->insert(std::pair<G4String,G4String>("liver_MI","FF_pork_liver_Peplow1998.dat"));
1720 fKnownMaterials->insert(std::pair<G4String,G4String>("heart_MI","FF_pork_heart_Peplow1998.dat"));
1721 fKnownMaterials->insert(std::pair<G4String,G4String>("blood_MI","FF_beef_blood_Peplow1998.dat"));
1722 fKnownMaterials->insert(std::pair<G4String,G4String>("grayMatter_MI","FF_gbrain_DeFelici2008.dat"));
1723 fKnownMaterials->insert(std::pair<G4String,G4String>("whiteMatter_MI","FF_wbrain_DeFelici2008.dat"));
1724 fKnownMaterials->insert(std::pair<G4String,G4String>("bone_MI","FF_bone_King2011.dat"));
1725 fKnownMaterials->insert(std::pair<G4String,G4String>("FatLowX_MI","FF_fat_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1726 fKnownMaterials->insert(std::pair<G4String,G4String>("BoneMatrixLowX_MI","FF_bonematrix_Tartari2002_joint_lowXdata.dat"));
1727 fKnownMaterials->insert(std::pair<G4String,G4String>("PMMALowX_MI", "FF_PMMA_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1728 fKnownMaterials->insert(std::pair<G4String,G4String>("dryBoneLowX_MI","FF_drybone_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1729 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS30-70_MI","FF_CIRS30-70_Poletti2002.dat"));
1730 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS50-50_MI","FF_CIRS50-50_Poletti2002.dat"));
1731 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS70-30_MI","FF_CIRS70-30_Poletti2002.dat"));
1732 fKnownMaterials->insert(std::pair<G4String,G4String>("RMI454_MI", "FF_RMI454_Poletti2002.dat"));
1733 fKnownMaterials->insert(std::pair<G4String,G4String>("PMMA_MI","FF_PMMA_Tartari2002.dat"));
1734 fKnownMaterials->insert(std::pair<G4String,G4String>("Lexan_MI","FF_lexan_Peplow1998.dat"));
1735 fKnownMaterials->insert(std::pair<G4String,G4String>("Kapton_MI","FF_kapton_Peplow1998.dat"));
1736 fKnownMaterials->insert(std::pair<G4String,G4String>("Nylon_MI","FF_nylon_Kosanetzky1987.dat"));
1737 fKnownMaterials->insert(std::pair<G4String,G4String>("Polyethylene_MI","FF_polyethylene_Kosanetzky1987.dat"));
1738 fKnownMaterials->insert(std::pair<G4String,G4String>("Polystyrene_MI","FF_polystyrene_Kosanetzky1987.dat"));
1739 fKnownMaterials->insert(std::pair<G4String,G4String>("Formaline_MI","FF_formaline_Peplow1998.dat"));
1740 fKnownMaterials->insert(std::pair<G4String,G4String>("Acetone_MI","FF_acetone_Cozzini2010.dat"));
1741 fKnownMaterials->insert(std::pair<G4String,G4String>("Hperoxide_MI","FF_Hperoxide_Cozzini2010.dat"));
1742 }
1743
1744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1745
1746 G4double G4PenelopeRayleighModelMI::IntegrateFun(G4double y[], G4int n, G4double dTheta)
1747 {
1748 G4double integral = 0.;
1749 for (G4int k=0; k<n-1; k++) {
1750 integral += (y[k]+y[k+1]);
1751 }
1752 integral *= dTheta*0.5;
1753 return integral;
1754 }
1755
1756 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1757 G4MIData* G4PenelopeRayleighModelMI::GetMIData(const G4Material* material)
1758 {
1759 if (material->IsExtended()) {
1760 G4ExtendedMaterial* aEM = (G4ExtendedMaterial*)material;
1761 G4MIData* dataMI = (G4MIData*)aEM->RetrieveExtension("MI");
1762 return dataMI;
1763 } else {
1764 return nullptr;
1765 }
1766 }
1767