| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Author: Luciano Pandola and Gianfranco Pate
27 //
28 // -------------------------------------------
29 // History:
30 // 03 Dec 2009 L. Pandola 1st implementati
31 // 25 May 2011 L. Pandola Renamed (make v2
32 // 27 Sep 2013 L. Pandola Migration to MT
33 // 20 Aug 2017 G. Paterno Molecular Interf
34 // 24 Mar 2019 G. Paterno Improved Molecul
35 // 20 Jun 2020 G. Paterno Read qext separa
36 // 27 Aug 2020 G. Paterno Further improvem
37 //
38 // -------------------------------------------
39 // Class description:
40 // Low Energy Electromagnetic Physics, Rayleig
41 // with the model from Penelope, version 2008
42 // -------------------------------------------
43 //
44 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
65
66 const G4int G4PenelopeRayleighModelMI::fMaxZ;
67 G4PhysicsFreeVector* G4PenelopeRayleighModelMI
68 G4PhysicsFreeVector* G4PenelopeRayleighModelMI
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 G4PenelopeRayleighModelMI::G4PenelopeRayleighM
73 const G4String& nam) :
74 G4VEmModel(nam),
75 fParticleChange(nullptr),fParticle(nullptr),
76 fSamplingTable(nullptr),fMolInterferenceData
77 fIsInitialised(false),fLocalTable(false),fIs
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 openin
92 // 4 = entering in methods
93
94 //build the energy grid. It is the same for
95 G4double logenergy = G4Log(fIntrinsicLowEner
96 G4double logmaxenergy = G4Log(1.5*fIntrinsic
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........oo
112
113 G4PenelopeRayleighModelMI::~G4PenelopeRayleigh
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........oo
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 expl
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 explicitl
173 }
174
175 return;
176 }
177
178 //....oooOO0OOooo........oooOO0OOooo........oo
179
180 void G4PenelopeRayleighModelMI::Initialise(con
181 const G4DataVector& )
182 {
183 if (fVerboseLevel > 3)
184 G4cout << "Calling G4PenelopeRayleighModel
185
186 SetParticle(part);
187
188 if (fVerboseLevel)
189 G4cout << "# Molecular Interference is " <
190
191 //Only the master model creates/fills/destro
192 if (IsMaster() && part == fParticle) {
193 //clear tables depending on materials, not
194 ClearTables();
195
196 //Use here the highest verbosity, from G4E
197 G4int globVerb = G4EmParameters::Instance(
198 if (globVerb > fVerboseLevel)
199 {
200 fVerboseLevel = globVerb;
201 if (fVerboseLevel)
202 G4cout << "Verbosity level of G4PenelopeRa
203 " from G4EmParameters()" << G4endl;
204 }
205 if (fVerboseLevel > 3)
206 G4cout << "Calling G4PenelopeRayleighMod
207
208 //Load the list of known materials and the
209 if (fIsMIActive)
210 {
211 if (!fKnownMaterials)
212 fKnownMaterials = new std::map<G4String,G4
213 if (!fKnownMaterials->size())
214 LoadKnownMIFFMaterials();
215 if (!fAngularFunction)
216 {
217 //Create and fill once
218 fAngularFunction = new G4PhysicsFreeVect
219 CalculateThetaAndAngFun(); //angular fun
220 }
221 }
222 if (fIsMIActive && !fMolInterferenceData)
223 fMolInterferenceData = new std::map<G4St
224 if (!fLogFormFactorTable)
225 fLogFormFactorTable = new std::map<const
226 if (!fPMaxTable)
227 fPMaxTable = new std::map<const G4Materi
228 if (!fSamplingTable)
229 fSamplingTable = new std::map<const G4Ma
230
231 //loop on the used materials
232 G4ProductionCutsTable* theCoupleTable = G4
233
234 for (G4int i=0;i<(G4int)theCoupleTable->Ge
235 const G4Material* material =
236 theCoupleTable->GetMaterialCutsCouple(i)->Ge
237 const G4ElementVector* theElementVector
238
239 for (std::size_t j=0;j<material->GetNumb
240 G4int iZ = theElementVector->at(j)->GetZasIn
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
248 ReadMolInterferenceData(material->GetName())
249
250 //2) If the table has not been built for
251 if (!fLogFormFactorTable->count(material
252 BuildFormFactorTable(material);
253
254 //3) retrieve or build the sampling tabl
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 m
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........oo
279
280 void G4PenelopeRayleighModelMI::InitialiseLoca
281 G4VEmModel *masterModel)
282 {
283 if (fVerboseLevel > 3)
284 G4cout << "Calling G4PenelopeRayleighMode
285
286 //Check that particle matches: one might hav
287 //(e.g. for e+ and e-)
288 if (part == fParticle) {
289
290 //Get the const table pointers from the ma
291 const G4PenelopeRayleighModelMI* theModel
292 static_cast<G4PenelopeRayleighModelMI*>
293
294 //Copy pointers to the data tables
295 for(G4int i=0; i<=fMaxZ; ++i)
296 {
297 fLogAtomicCrossSection[i] = theModel->fLogAt
298 fAtomicFormFactor[i] = theModel->fAtomicForm
299 }
300 fMolInterferenceData = theModel->fMolInter
301 fLogFormFactorTable = theModel->fLogFormFa
302 fPMaxTable = theModel->fPMaxTable;
303 fSamplingTable = theModel->fSamplingTable;
304 fKnownMaterials = theModel->fKnownMaterial
305 fAngularFunction = theModel->fAngularFunct
306
307 //Copy the G4DataVector with the grid
308 fLogQSquareGrid = theModel->fLogQSquareGri
309
310 //Same verbosity for all workers, as the m
311 fVerboseLevel = theModel->fVerboseLevel;
312 }
313 return;
314 }
315
316 //....oooOO0OOooo........oooOO0OOooo........oo
317
318 namespace {G4Mutex PenelopeRayleighModelMutex
319 G4double G4PenelopeRayleighModelMI::ComputeCro
320 G4double energy,
321 G4double Z,
322 G4double,
323 G4double,
324 G4double)
325 {
326 //Cross section of Rayleigh scattering in Pe
327 //tabulation, Cuellen et al. (1997), with no
328 //et al. J. Phys. Chem. Ref. Data 4 (1975) 4
329
330 if (fVerboseLevel > 3)
331 G4cout << "Calling CrossSectionPerAtom() o
332
333 G4int iZ = G4int(Z);
334 if (!fLogAtomicCrossSection[iZ]) {
335 //If we are here, it means that Initialize
336 //not filled up. This can happen in a Unit
337 if (fVerboseLevel > 0) {
338 //Issue a G4Exception (warning) only in
339 G4ExceptionDescription ed;
340 ed << "Unable to retrieve the cross sect
341 ed << "This can happen only in Unit Test
342 G4Exception("G4PenelopeRayleighModelMI::
343 "em2040",JustWarning,ed);
344 }
345
346 //protect file reading via autolock
347 G4AutoLock lock(&PenelopeRayleighModelMute
348 ReadDataFile(iZ);
349 lock.unlock();
350 }
351
352 G4double cross = 0;
353 G4PhysicsFreeVector* atom = fLogAtomicCrossS
354 if (!atom) {
355 G4ExceptionDescription ed;
356 ed << "Unable to find Z=" << iZ << " in th
357 G4Exception("G4PenelopeRayleighModelMI::Co
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 " <<
368 << Z << " = " << cross/barn << " barn" <<
369 }
370 return cross;
371 }
372
373 //....oooOO0OOooo........oooOO0OOooo........oo
374
375 void G4PenelopeRayleighModelMI::CalculateTheta
376 {
377 G4double theta = 0;
378 for(G4int k=0; k<fNtheta; k++) {
379 theta += fDTheta;
380 G4double value = (1+std::cos(theta)*std::c
381 fAngularFunction->PutValue(k,theta,value);
382 if (fVerboseLevel > 3)
383 G4cout << "theta[" << k << "]: " << fAn
384 << ", angFun[" << k << "]: " << (*fAngu
385 }
386 }
387
388 //....oooOO0OOooo........oooOO0OOooo........oo
389
390 G4double G4PenelopeRayleighModelMI::CalculateQ
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, la
400 << ", x: " << x*nm << ", q: " << q << G4e
401 }
402 q2 = std::pow(q,2);
403 return q2;
404 }
405
406 //....oooOO0OOooo........oooOO0OOooo........oo
407
408 //Overriding of parent's (G4VEmModel) method
409 G4double G4PenelopeRayleighModelMI::CrossSecti
410 const G4ParticleDefinition* p,
411 G4double energy,
412 G4double,
413 G4double)
414 {
415 //check if we are in a Unit Test (only for t
416 static G4bool amInAUnitTest = false;
417 if (G4ProductionCutsTable::GetProductionCuts
418 {
419 amInAUnitTest = true;
420 G4ExceptionDescription ed;
421 ed << "The ProductionCuts table is empty
422 ed << "This should happen only in Unit T
423 G4Exception("G4PenelopeRayleighModelMI::
424 "em2019",JustWarning,ed);
425 }
426 //If the material does not have a MIFF, cont
427 G4String matname = material->GetName();
428 if (amInAUnitTest)
429 {
430 //No need to lock, as this is always cal
431 const G4ElementVector* theElementVector
432 //protect file reading via autolock
433 for (std::size_t j=0;j<material->GetNumb
434 G4int iZ = theElementVector->at(j)->GetZasIn
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 && (fMolInterfe
449 if (!useMIFF)
450 {
451 if (fVerboseLevel > 2)
452 G4cout << "Rayleigh CS of: " << matname << "
453 return G4VEmModel::CrossSectionPerVolume
454 }
455
456 // If we are here, it means that we have to
457 if (fVerboseLevel > 2)
458 G4cout << "Rayleigh CS of: " << matname
459 << " calculated through integration of th
460
461 G4double cs = 0;
462
463 //force null cross-section if below the low-
464 if (energy < LowEnergyLimit())
465 return cs;
466
467 //if the material is a CRYSTAL, forbid this
468 if (material->IsExtended() && material->GetN
469 G4ExtendedMaterial* extendedmaterial = (G4
470 G4CrystalExtension* crystalExtension = (G4
471 if (crystalExtension != 0) {
472 G4cout << "The material has a crystallin
473 return 0;
474 }
475 }
476
477 //GET MATERIAL INFORMATION
478 G4double atomDensity = material->GetTotNbOfA
479 std::size_t nElements = material->GetNumberO
480 const G4ElementVector* elementVector = mater
481 const G4double* fractionVector = material->G
482
483 //Stoichiometric factors
484 std::vector<G4double> *StoichiometricFactors
485 for (std::size_t i=0;i<nElements;++i) {
486 G4double fraction = fractionVector[i];
487 G4double atomicWeigth = (*elementVector)[i
488 StoichiometricFactors->push_back(fraction/
489 }
490 G4double MaxStoichiometricFactor = 0.;
491 for (std::size_t i=0;i<nElements;++i) {
492 if ((*StoichiometricFactors)[i] > MaxStoic
493 MaxStoichiometricFactor = (*Stoichiometr
494 }
495 for (std::size_t i=0;i<nElements;++i) {
496 (*StoichiometricFactors)[i] /= MaxStoichi
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 = atomDe
505
506 if (fVerboseLevel > 2)
507 G4cout << "Material " << material->GetName
508 << "per molecule and " << moleculeDensity
509
510 //Equivalent molecular weight (dimensionless
511 G4double MolWeight = 0.;
512 for (std::size_t i=0;i<nElements;++i)
513 MolWeight += (*StoichiometricFactors)[i]*(
514
515 if (fVerboseLevel > 2)
516 G4cout << "Molecular weight of " << matnam
517
518 G4double IntegrandFun[fNtheta];
519 for (G4int k=0; k<fNtheta; k++) {
520 G4double theta = fAngularFunction->Energy(
521 G4double F2 = GetFSquared(material,Calcula
522 IntegrandFun[k] = (*fAngularFunction)[k]*F
523 }
524
525 G4double constant = pi*classic_electr_radius
526 cs = constant*IntegrateFun(IntegrandFun,fNth
527
528 //Now cs is the cross section per molecule,
529 G4double csvolume = cs*moleculeDensity;
530
531 //print CS and mfp
532 if (fVerboseLevel > 2)
533 G4cout << "Rayleigh CS of " << matname <<
534 << " keV: " << cs/barn << " barn"
535 << ", mean free path: " << 1./csvolume/mm
536
537 delete StoichiometricFactors;
538 //return CS
539 return csvolume;
540 }
541
542 //....oooOO0OOooo........oooOO0OOooo........oo
543
544 void G4PenelopeRayleighModelMI::BuildFormFacto
545 {
546 if (fVerboseLevel > 3)
547 G4cout << "Calling BuildFormFactorTable()
548
549 //GET MATERIAL INFORMATION
550 std::size_t nElements = material->GetNumberO
551 const G4ElementVector* elementVector = mater
552 const G4double* fractionVector = material->G
553
554 //Stoichiometric factors
555 std::vector<G4double> *StoichiometricFactors
556 for (std::size_t i=0;i<nElements;++i) {
557 G4double fraction = fractionVector[i];
558 G4double atomicWeigth = (*elementVector)[i
559 StoichiometricFactors->push_back(fraction/
560 }
561 G4double MaxStoichiometricFactor = 0.;
562 for (std::size_t i=0;i<nElements;++i) {
563 if ((*StoichiometricFactors)[i] > MaxStoic
564 MaxStoichiometricFactor = (*Stoichiometr
565 }
566 if (MaxStoichiometricFactor<1e-16) {
567 G4ExceptionDescription ed;
568 ed << "Inconsistent data of atomic composi
569 G4Exception("G4PenelopeRayleighModelMI::Bu
570 "em2042",FatalException,ed);
571 }
572 for (std::size_t i=0;i<nElements;++i)
573 (*StoichiometricFactors)[i] /= MaxStoichi
574
575 //Equivalent molecular weight (dimensionless
576 G4double MolWeight = 0.;
577 for (std::size_t i=0;i<nElements;++i)
578 MolWeight += (*StoichiometricFactors)[i]*(
579
580 //CREATE THE FORM FACTOR TABLE
581 //First, the form factors are retrieved [F/s
582 //Then, they are squared and multiplied for
583 //This makes difference for CS calculation,
584 G4PhysicsFreeVector* theFFVec = new G4Physic
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 us
596 if (fIsMIActive && aFileNameFF != "") {
597 if (fVerboseLevel)
598 G4cout << "Read MIFF for " << matname <<
599
600 ReadMolInterferenceData(matname,aFileNameF
601 G4PhysicsFreeVector* Fvector = fMolInterf
602
603 for (std::size_t k=0;k<fLogQSquareGrid.siz
604 G4double q = std::pow(G4Exp(fLogQSquareG
605 G4double f = Fvector->Value(q);
606 G4double ff2 = f*f*MolWeight;
607 if (ff2)
608 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo
609 }
610 }
611 //retrieve the MIFF from the database or use
612 else {
613 //medical material: composition of fat, wa
614 if (fIsMIActive && matname.find("MedMat")
615 if (fVerboseLevel)
616 G4cout << "Build MIFF from components for "
617
618 //get the material composition from its
619 G4int ki, kf=6, ktot=19;
620 G4double comp[4];
621 G4String compstring = matname.substr(kf+
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 << "]:
629 }
630
631 const char* path = G4FindDataDir("G4LEDA
632 if (!path) {
633 G4String excep = "G4LEDATA environment varia
634 G4Exception("G4PenelopeRayleighModelMI::Buil
635 "em0006",FatalException,excep);
636 }
637
638 if (!fMolInterferenceData->count("Fat_MI
639 ReadMolInterferenceData("Fat_MI");
640 G4PhysicsFreeVector* fatFF = fMolInterfe
641
642 if (!fMolInterferenceData->count("Water_
643 ReadMolInterferenceData("Water_MI");
644 G4PhysicsFreeVector* waterFF = fMolInter
645
646 if (!fMolInterferenceData->count("BoneMa
647 ReadMolInterferenceData("BoneMatrix_MI");
648 G4PhysicsFreeVector* bonematrixFF = fMol
649
650 if (!fMolInterferenceData->count("Minera
651 ReadMolInterferenceData("Mineral_MI");
652 G4PhysicsFreeVector* mineralFF = fMolInt
653
654 //get and combine the molecular form fac
655 for (std::size_t k=0;k<fLogQSquareGrid.s
656 G4double ff2 = 0;
657 G4double q = std::pow(G4Exp(fLogQSquareGrid[
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*fwate
663 comp[2]*fbonematrix*fbonematrix+comp[3]*fm
664 ff2 *= MolWeight;
665 if (ff2) theFFVec->PutValue(k,fLogQSquareGri
666 }
667 }
668 //other materials with MIFF (from the data
669 else if (fIsMIActive && fMolInterferenceDa
670 if (fVerboseLevel)
671 G4cout << "Read MIFF from database " << matn
672 G4PhysicsFreeVector* FF = fMolInterferen
673 for (std::size_t k=0;k<fLogQSquareGrid.s
674 G4double ff2 = 0;
675 G4double q = std::pow(G4Exp(fLogQSquareGrid[
676 G4double f = FF->Value(q);
677 ff2 = f*f*MolWeight;
678 if (ff2) theFFVec->PutValue(k,fLogQSquareGri
679 }
680 }
681 //IAM
682 else {
683 if (fVerboseLevel)
684 G4cout << "FF of " << matname << " calculate
685 for (std::size_t k=0;k<fLogQSquareGrid.s
686 G4double ff2 = 0;
687 for (std::size_t i=0;i<nElements;++i) {
688 G4int iZ = (*elementVector)[i]->GetZasInt(
689 G4PhysicsFreeVector* theAtomVec = fAtomicF
690 G4double q = std::pow(G4Exp(fLogQSquareGri
691 G4double f = theAtomVec->Value(q);
692 ff2 += f*f*(*StoichiometricFactors)[i];
693 }
694 if (ff2) theFFVec->PutValue(k,fLogQSquareGri
695 }
696 }
697 }
698 theFFVec->FillSecondDerivatives();
699 fLogFormFactorTable->insert(std::make_pair(m
700
701 if (fVerboseLevel > 3)
702 DumpFormFactorTable(material);
703 delete StoichiometricFactors;
704
705 return;
706 }
707
708 //....oooOO0OOooo........oooOO0OOooo........oo
709
710 void G4PenelopeRayleighModelMI::SampleSecondar
711 const G4MaterialCutsCouple* coup
712 const G4DynamicParticle* aDynami
713 G4double,
714 G4double)
715 {
716 // Sampling of the Rayleigh final state (nam
717 // from the Penelope2008 model. The scatteri
718 // cross section dOmega/d(cosTheta) from Bor
719 // anomalous scattering effects. The Form Fa
720 // analytical cross section is retrieved via
721 // are tabulated for F(Q). Form factor for c
722 // the additivity rule. The sampling from th
723 // Transform with Aliasing (RITA) algorithm;
724 // for each material and managed by G4Penelo
725 // The sampling algorithm (rejection method)
726 // increases with energy. For E=100 keV the
727 // hydrogen and uranium, respectively.
728
729 if (fVerboseLevel > 3)
730 G4cout << "Calling SamplingSecondaries() o
731
732 G4double photonEnergy0 = aDynamicGamma->GetK
733
734 if (photonEnergy0 <= fIntrinsicLowEnergyLimi
735 fParticleChange->ProposeTrackStatus(fStopA
736 fParticleChange->SetProposedKineticEnergy(
737 fParticleChange->ProposeLocalEnergyDeposit
738 return;
739 }
740
741 G4ParticleMomentum photonDirection0 = aDynam
742 const G4Material* theMat = couple->GetMateri
743
744 //1) Verify if tables are ready
745 //Either Initialize() was not called, or we
746 //not invoked
747 if (!fPMaxTable || !fSamplingTable || !fLogF
748 //create a **thread-local** version of the
749 //Unit Tests
750 fLocalTable = true;
751 if (!fLogFormFactorTable)
752 fLogFormFactorTable = new std::map<const
753 if (!fPMaxTable)
754 fPMaxTable = new std::map<const G4Materi
755 if (!fSamplingTable)
756 fSamplingTable = new std::map<const G4Ma
757 if (fIsMIActive && !fMolInterferenceData)
758 fMolInterferenceData = new std::map<G4St
759 }
760
761 if (!fSamplingTable->count(theMat)) {
762 //If we are here, it means that Initialize
763 //not filled up. This can happen in a Unit
764 if (fVerboseLevel > 0) {
765 //Issue a G4Exception (warning) only in
766 G4ExceptionDescription ed;
767 ed << "Unable to find the fSamplingTable
768 theMat->GetName() << G4endl;
769 ed << "This can happen only in Unit Test
770 G4Exception("G4PenelopeRayleighModelMI::
771 "em2019",JustWarning,ed);
772 }
773 const G4ElementVector* theElementVector =
774 //protect file reading via autolock
775 G4AutoLock lock(&PenelopeRayleighModelMute
776 for (std::size_t j=0;j<theMat->GetNumberOf
777 G4int iZ = theElementVector->at(j)->GetZ
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 t
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 = fSamp
803 G4PhysicsFreeVector* thePMax = fPMaxTable->f
804 G4double cosTheta = 1.0;
805
806 //OK, ready to go!
807 G4double qmax = 2.0*photonEnergy0/electron_m
808
809 if (qmax < 1e-10) { //very low momentum tran
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->GetNumbe
821 G4double LastQ2inTheTable = theDataTable->
822 G4double q2max = std::min(qmax*qmax,LastQ2
823
824 G4bool loopAgain = false;
825 G4double MaxPValue = thePMax->Value(photon
826 G4double xx=0;
827
828 //Sampling by rejection method. The reject
829 //G = 0.5*(1+cos^2(theta))
830 do {
831 loopAgain = false;
832 G4double RandomMax = G4UniformRand()*Max
833 xx = theDataTable->SampleValue(RandomMax
834
835 //xx is a random value of q^2 in (0,q2ma
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*cos
847
848 //Scattered photon angles. ( Z - axis along
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
855 G4ThreeVector photonDirection1(dirX, dirY, d
856
857 photonDirection1.rotateUz(photonDirection0);
858 fParticleChange->ProposeMomentumDirection(ph
859 fParticleChange->SetProposedKineticEnergy(ph
860
861 return;
862 }
863
864 //....oooOO0OOooo........oooOO0OOooo........oo
865
866 void G4PenelopeRayleighModelMI::ReadDataFile(c
867 {
868 if (fVerboseLevel > 2) {
869 G4cout << "G4PenelopeRayleighModelMI::Read
870 G4cout << "Going to read Rayleigh data fil
871 }
872
873 const char* path = G4FindDataDir("G4LEDATA")
874 if (!path) {
875 G4String excep = "G4LEDATA environment var
876 G4Exception("G4PenelopeRayleighModelMI::Re
877 "em0006",FatalException,excep);
878 return;
879 }
880
881 //Read first the cross section file (all the
882 std::ostringstream ost;
883 if (Z>9)
884 ost << path << "/penelope/rayleigh/pdgra"
885 else
886 ost << path << "/penelope/rayleigh/pdgra0"
887 std::ifstream file(ost.str().c_str());
888
889 if (!file.is_open()) {
890 G4String excep = "Data file " + G4String(o
891 G4Exception("G4PenelopeRayleighModelMI::Re
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 >=
900 G4ExceptionDescription ed;
901 ed << "Corrupted data file for Z=" << Z <<
902 G4Exception("G4PenelopeRayleighModelMI::Re
903 "em0005",FatalException,ed);
904 return;
905 }
906
907 fLogAtomicCrossSection[Z] = new G4PhysicsFre
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,G4Lo
915 if (file.eof() && i != (nPoints-1)) { //fi
916 G4ExceptionDescription ed ;
917 ed << "Corrupted data file for Z=" << Z
918 ed << "Found less than " << nPoints << "
919 G4Exception("G4PenelopeRayleighModelMI::
920 "em0005",FatalException,ed);
921 }
922 }
923 file.close();
924
925 //Then, read the extended momentum transfer
926 std::ostringstream ost2;
927 ost2 << path << "/penelope/rayleigh/MIFF/qex
928 file.open(ost2.str().c_str());
929
930 if (!file.is_open()) {
931 G4String excep = "Data file " + G4String(o
932 G4Exception("G4PenelopeRayleighModelMI::Re
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 o
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"
953 else
954 ost3 << path << "/penelope/rayleigh/pdaff0
955 file.open(ost3.str().c_str());
956
957 if (!file.is_open()) {
958 G4String excep = "Data file " + G4String(o
959 G4Exception("G4PenelopeRayleighModelMI::Re
960 "em0003",FatalException,excep);
961 }
962
963 file >> readZ >> nPoints;
964
965 if (readZ != Z || nPoints <= 0 || nPoints >=
966 G4ExceptionDescription ed;
967 ed << "Corrupted data file for Z=" << Z <<
968 G4Exception("G4PenelopeRayleighModelMI::Re
969 "em0005",FatalException,ed);
970 return;
971 }
972
973 fAtomicFormFactor[Z] = new G4PhysicsFreeVect
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 d
977 fAtomicFormFactor[Z]->PutValue(i,q,ff);
978 if (file.eof() && i != (nPoints-1)) { //fi
979 G4ExceptionDescription ed;
980 ed << "Corrupted data file for Z=" << Z
981 ed << "Found less than " << nPoints << "
982 G4Exception("G4PenelopeRayleighModelMI::
983 "em0005",FatalException,ed);
984 }
985 }
986 file.close();
987 return;
988 }
989
990 //....oooOO0OOooo........oooOO0OOooo........oo
991
992 void G4PenelopeRayleighModelMI::ReadMolInterfe
993 const G4String& FFfilename)
994 {
995 if (fVerboseLevel > 2) {
996 G4cout << "G4PenelopeRayleighModelMI::Read
997 matname << G4endl;
998 }
999 G4bool isLocalFile = (FFfilename != "NULL");
1000
1001 const char* path = G4FindDataDir("G4LEDATA"
1002 if (!path) {
1003 G4String excep = "G4LEDATA environment va
1004 G4Exception("G4PenelopeRayleighModelMI::R
1005 "em0006",FatalException,excep);
1006 return;
1007 }
1008
1009 if (!(fKnownMaterials->count(matname)) && !
1010 return;
1011
1012 G4String aFileName = (isLocalFile) ? FFfile
1013
1014 //if the material has a MIFF, read it from
1015 if (aFileName != "") {
1016 if (fVerboseLevel > 2)
1017 G4cout << "ReadMolInterferenceData(). R
1018 aFileName << " " <<
1019 (isLocalFile ? "(local)" : "(database)") <<
1020 std::ifstream file;
1021 std::ostringstream ostIMFF;
1022 if (isLocalFile)
1023 ostIMFF << aFileName;
1024 else
1025 ostIMFF << path << "/penelope/rayleigh/
1026 file.open(ostIMFF.str().c_str());
1027
1028 if (!file.is_open()) {
1029 G4String excep = "Data file " + G4Strin
1030 G4Exception("G4PenelopeRayleighModelMI:
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: " << nPoi
1046
1047 //read the file
1048 file.open(ostIMFF.str().c_str());
1049 G4PhysicsFreeVector* theFFVec = new G4Phy
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 u
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 << " en
1062 G4Exception("G4PenelopeRayleighModelMI::Rea
1063 "em1005",FatalException,ed);
1064 }
1065 }
1066 if (!fMolInterferenceData) {
1067 G4Exception("G4PenelopeRayleighModelMI:
1068 "em2145",FatalException,
1069 "Unable to allocate the Molecular Inter
1070 delete theFFVec;
1071 return;
1072 }
1073 file.close();
1074 fMolInterferenceData->insert(std::make_pa
1075 }
1076 return;
1077 }
1078
1079 //....oooOO0OOooo........oooOO0OOooo........o
1080
1081 G4double G4PenelopeRayleighModelMI::GetFSquar
1082 {
1083 G4double f2 = 0;
1084 //Input value QSquared could be zero: prote
1085 //the FPE exception
1086
1087 //If Q<1e-10, set Q to 1e-10
1088 G4double logQSquared = (QSquared>1e-10) ? G
1089 //last value of the table
1090 G4double maxlogQ2 = fLogQSquareGrid[fLogQSq
1091
1092 //now it should be all right
1093 G4PhysicsFreeVector* theVec = fLogFormFacto
1094
1095 if (!theVec) {
1096 G4ExceptionDescription ed;
1097 ed << "Unable to retrieve F squared table
1098 G4Exception("G4PenelopeRayleighModelMI::G
1099 "em2046",FatalException,ed);
1100 return 0;
1101 }
1102
1103 if (logQSquared < -20) { //Q < 1e-9
1104 G4double logf2 = (*theVec)[0]; //first va
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(logQSquare
1112 f2 = G4Exp(logf2);
1113 }
1114
1115 if (fVerboseLevel > 3) {
1116 G4cout << "G4PenelopeRayleighModelMI::Get
1117 G4cout << "Q^2 = " << QSquared << " (uni
1118 }
1119 return f2;
1120 }
1121
1122 //....oooOO0OOooo........oooOO0OOooo........o
1123
1124 void G4PenelopeRayleighModelMI::InitializeSam
1125 {
1126 G4double q2min = 0;
1127 G4double q2max = 0;
1128 const std::size_t np = 150; //hard-coded in
1129 for (std::size_t i=1;i<fLogQSquareGrid.size
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:
1143 "em2047",FatalException,
1144 "Too few points to initialize the sampl
1145 }
1146 if (q2min > (q2max-1e-10))
1147 {
1148 G4cout << "q2min= " << q2min << " q2max
1149 G4Exception("G4PenelopeRayleighModelMI:
1150 "em2048",FatalException,
1151 "Too narrow grid to initialize the samp
1152 }
1153
1154 //This is subroutine RITAI0 of Penelope
1155 //Create an object of type G4PenelopeRaylei
1156
1157 //temporary vectors --> Then everything is
1158 G4DataVector* x = new G4DataVector();
1159
1160 /******************************************
1161 Start with a grid of NUNIF points uniform
1162 *******************************************
1163 std::size_t NUNIF = std::min(std::max(((std
1164 const G4int nip = 51; //hard-coded in Penel
1165
1166 G4double dx = (q2max-q2min)/((G4double) NUN
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() <<
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
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])/(G4d
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,
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
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)[
1216 sumi->push_back(next);
1217 }
1218
1219 G4double lastIntegral = (*sumi)[sumi->s
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 i
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]*
1234 G4double B_temp = 1.0-1.0/(pli*pui*dx*d
1235 G4double A_temp = (1.0/(pli*dx))-1.0-B_
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 t
1251 //and the true pdf, extended over the i
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)[
1263 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(
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 appr
1272 if ((*err)[i] > 0.1*(*area)[i] && icase =
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() != NUNI
1296 err->size() != NUNIF || area->size() !=
1297 {
1298 G4ExceptionDescription ed;
1299 ed << "Problem in building the Table fo
1300 G4Exception("G4PenelopeRayleighModelMI:
1301 "em2049",FatalException,ed);
1302 }
1303
1304 /******************************************
1305 New grid points are added by halving the s
1306 This is done up to np=150 points in the gri
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;++
1313 {
1314 //maxError is the lagest of the interval
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 poin
1323 G4double newx = 0.5*((*x)[iErrMax]+(*x)
1324
1325 x->insert(x->begin()+iErrMax+1,newx);
1326 //Add place-holders in the other vector
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
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])/(G4dou
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(
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(m
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*((*pd
1366 sumi->push_back(next);
1367 }
1368 G4double lastIntegral = (*sumi)[sumi->siz
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 int
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]*fa
1384 G4double B_temp = 1.0-1.0/(pli*pui*dxLoca
1385 G4double A_temp = (1.0/(pli*dxLocal))-1.0
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
1400 //and the true pdf, extended over the int
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]
1411 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+
1412 if (k == 0 || k == nip-1)
1413 (*err)[i] += 0.5*std::fabs(pap-(*pdfi
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
1420 if ((*err)[i] > 0.1*(*area)[i] && ica
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:
1439 "em2050",FatalException,
1440 "Problem in building the extended Table
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 dis
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 binar
1470 *******************************************
1471 std::vector<std::size_t> *ITTL = new std::v
1472 std::vector<std::size_t> *ITTU = new std::v
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 &&
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() != n
1502 {
1503 G4Exception("G4PenelopeRayleighModelMI:
1504 "em2051",FatalException,
1505 "Problem in building the Limit Tables f
1506 }
1507
1508 /******************************************
1509 Copy tables
1510 *******************************************
1511 G4PenelopeSamplingData* theTable = new G4Pe
1512 for (std::size_t i=0;i<np;++i)
1513 {
1514 theTable->AddPoint((*x)[i],(*PAC)[i],(*
1515 }
1516 if (fVerboseLevel > 2)
1517 {
1518 G4cout << "****************************
1519 G4endl;
1520 G4cout << "Sampling table for Penelope
1521 theTable->DumpTable();
1522 }
1523 fSamplingTable->insert(std::make_pair(mat,t
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........o
1541
1542 void G4PenelopeRayleighModelMI::GetPMaxTable(
1543 {
1544 if (!fPMaxTable)
1545 {
1546 G4cout << "G4PenelopeRayleighModelMI::B
1547 G4cout << "Going to instanziate the fPM
1548 G4cout << "That should _not_ be here! "
1549 fPMaxTable = new std::map<const G4Mater
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:
1559 "em2052",FatalException,
1560 "SamplingTable is not properly instanti
1561 return;
1562 }
1563
1564 //This should not be: the sampling table is
1565 if (!fSamplingTable->count(mat))
1566 {
1567 G4ExceptionDescription ed;
1568 ed << "Sampling table for material " <<
1569 G4Exception("G4PenelopeRayleighModelMI:
1570 "em2052",FatalException,
1571 ed);
1572 return;
1573 }
1574
1575 G4PenelopeSamplingData *theTable = fSamplin
1576 std::size_t tablePoints = theTable->GetNumb
1577 std::size_t nOfEnergyPoints = fLogEnergyGri
1578 G4PhysicsFreeVector* theVec = new G4Physics
1579
1580 const std::size_t nip = 51; //hard-coded in
1581
1582 for (std::size_t ie=0;ie<fLogEnergyGridPMax
1583 {
1584 G4double energy = G4Exp(fLogEnergyGridP
1585 G4double Qm = 2.0*energy/electron_mass_
1586 G4double Qm2 = Qm*Qm;
1587 G4double firstQ2 = theTable->GetX(0);
1588 G4double lastQ2 = theTable->GetX(tableP
1589 G4double thePMax = 0;
1590
1591 if (Qm2 > firstQ2)
1592 {
1593 if (Qm2 < lastQ2)
1594 {
1595 //bisection to look for the index of
1596 std::size_t lowerBound = 0;
1597 std::size_t upperBound = tablePoints-
1598 while (lowerBound <= upperBound)
1599 {
1600 std::size_t midBin = (lowerBound + uppe
1601 if( Qm2 < theTable->GetX(midBin))
1602 { upperBound = midBin-1; }
1603 else
1604 { lowerBound = midBin+1; }
1605 }
1606 //upperBound is the output (but also
1607 G4double Q1 = theTable->GetX(upperBou
1608 G4double Q2 = Qm2;
1609 G4double DQ = (Q2-Q1)/((G4double)(nip
1610 G4double theA = theTable->GetA(upperB
1611 G4double theB = theTable->GetB(upperB
1612 G4double thePAC = theTable->GetPAC(up
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
1624 etap = con2*(1.0-std::sqrt(1.0-2.0*ta
1625 else
1626 etap = tau/con2;
1627 G4double theFun = (theTable->GetPAC(upp
1628 (1.0+(theA+theB*etap)*etap)*(1.0+(the
1629 ((1.0-theB*etap*etap)*ci*(theTable->G
1630 fun->push_back(theFun);
1631 }
1632 //Now intergrate numerically the fun
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])*C
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)[
1644 (*fun)[hh+1]-(*fun)[hh-2])*HCON
1645 sum->push_back(next);
1646 }
1647 G4double last = (*sum)[nip-2]+(5.0*(*
1648 (*fun)[nip-3])*CONS;
1649 sum->push_back(last);
1650 thePMax = thePAC + (*sum)[sum->size()
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,theVe
1669 return;
1670 }
1671
1672 //....oooOO0OOooo........oooOO0OOooo........o
1673
1674 void G4PenelopeRayleighModelMI::DumpFormFacto
1675 {
1676 G4cout << "********************************
1677 G4cout << "G4PenelopeRayleighModelMI: Form
1678 //try to use the same format as Penelope-Fo
1679 G4cout << "Q/(m_e*c) F(Q)
1680 G4cout << "********************************
1681 if (!fLogFormFactorTable->count(mat))
1682 BuildFormFactorTable(mat);
1683
1684 G4PhysicsFreeVector* theVec = fLogFormFacto
1685 for (std::size_t i=0;i<theVec->GetVectorLen
1686 {
1687 G4double logQ2 = theVec->GetLowEdgeEner
1688 G4double Q = G4Exp(0.5*logQ2);
1689 G4double logF2 = (*theVec)[i];
1690 G4double F = G4Exp(0.5*logF2);
1691 G4cout << Q << " " << F <<
1692 }
1693 //DONE
1694 return;
1695 }
1696
1697 //....oooOO0OOooo........oooOO0OOooo........o
1698
1699 void G4PenelopeRayleighModelMI::SetParticle(c
1700 {
1701 if(!fParticle) {
1702 fParticle = p;
1703 }
1704 }
1705
1706 //....oooOO0OOooo........oooOO0OOooo........o
1707 void G4PenelopeRayleighModelMI::LoadKnownMIFF
1708 {
1709 fKnownMaterials->insert(std::pair<G4String,
1710 fKnownMaterials->insert(std::pair<G4String,
1711 fKnownMaterials->insert(std::pair<G4String,
1712 fKnownMaterials->insert(std::pair<G4String,
1713 fKnownMaterials->insert(std::pair<G4String,
1714 fKnownMaterials->insert(std::pair<G4String,
1715 fKnownMaterials->insert(std::pair<G4String,
1716 fKnownMaterials->insert(std::pair<G4String,
1717 fKnownMaterials->insert(std::pair<G4String,
1718 fKnownMaterials->insert(std::pair<G4String,
1719 fKnownMaterials->insert(std::pair<G4String,
1720 fKnownMaterials->insert(std::pair<G4String,
1721 fKnownMaterials->insert(std::pair<G4String,
1722 fKnownMaterials->insert(std::pair<G4String,
1723 fKnownMaterials->insert(std::pair<G4String,
1724 fKnownMaterials->insert(std::pair<G4String,
1725 fKnownMaterials->insert(std::pair<G4String,
1726 fKnownMaterials->insert(std::pair<G4String,
1727 fKnownMaterials->insert(std::pair<G4String,
1728 fKnownMaterials->insert(std::pair<G4String,
1729 fKnownMaterials->insert(std::pair<G4String,
1730 fKnownMaterials->insert(std::pair<G4String,
1731 fKnownMaterials->insert(std::pair<G4String,
1732 fKnownMaterials->insert(std::pair<G4String,
1733 fKnownMaterials->insert(std::pair<G4String,
1734 fKnownMaterials->insert(std::pair<G4String,
1735 fKnownMaterials->insert(std::pair<G4String,
1736 fKnownMaterials->insert(std::pair<G4String,
1737 fKnownMaterials->insert(std::pair<G4String,
1738 fKnownMaterials->insert(std::pair<G4String,
1739 fKnownMaterials->insert(std::pair<G4String,
1740 fKnownMaterials->insert(std::pair<G4String,
1741 fKnownMaterials->insert(std::pair<G4String,
1742 }
1743
1744 //....oooOO0OOooo........oooOO0OOooo........o
1745
1746 G4double G4PenelopeRayleighModelMI::Integrate
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........o
1757 G4MIData* G4PenelopeRayleighModelMI::GetMIDat
1758 {
1759 if (material->IsExtended()) {
1760 G4ExtendedMaterial* aEM = (G4ExtendedMate
1761 G4MIData* dataMI = (G4MIData*)aEM->Retrie
1762 return dataMI;
1763 } else {
1764 return nullptr;
1765 }
1766 }
1767