| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4PenelopeRayleighModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-04 $
26 // 28 //
27 // Author: Luciano Pandola 29 // Author: Luciano Pandola
28 // 30 //
29 // History: 31 // History:
30 // -------- 32 // --------
31 // 03 Dec 2009 L Pandola First implementa << 33 // 14 Oct 2008 L Pandola Migration from process to model
32 // 25 May 2011 L.Pandola Renamed (make v2 <<
33 // 19 Sep 2013 L.Pandola Migration to MT <<
34 // 34 //
35 35
36 #include "G4PenelopeRayleighModel.hh" 36 #include "G4PenelopeRayleighModel.hh"
37 #include "G4PhysicalConstants.hh" <<
38 #include "G4SystemOfUnits.hh" <<
39 #include "G4PenelopeSamplingData.hh" <<
40 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh"
41 #include "G4MaterialCutsCouple.hh" 38 #include "G4MaterialCutsCouple.hh"
42 #include "G4ProductionCutsTable.hh" 39 #include "G4ProductionCutsTable.hh"
43 #include "G4DynamicParticle.hh" 40 #include "G4DynamicParticle.hh"
44 #include "G4PhysicsTable.hh" 41 #include "G4PhysicsTable.hh"
45 #include "G4ElementTable.hh" 42 #include "G4ElementTable.hh"
46 #include "G4Element.hh" 43 #include "G4Element.hh"
47 #include "G4PhysicsFreeVector.hh" << 44 #include "G4PenelopeIntegrator.hh"
48 #include "G4AutoLock.hh" <<
49 #include "G4Exp.hh" <<
50 45
51 //....oooOO0OOooo........oooOO0OOooo........oo <<
52 <<
53 const G4int G4PenelopeRayleighModel::fMaxZ; <<
54 G4PhysicsFreeVector* G4PenelopeRayleighModel:: <<
55 G4PhysicsFreeVector* G4PenelopeRayleighModel:: <<
56 46
57 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 48
59 G4PenelopeRayleighModel::G4PenelopeRayleighMod << 49
60 const G4String& nam) << 50 G4PenelopeRayleighModel::G4PenelopeRayleighModel(const G4ParticleDefinition*,
61 :G4VEmModel(nam),fParticleChange(nullptr),fP << 51 const G4String& nam)
62 fLogFormFactorTable(nullptr),fPMaxTable(nul << 52 :G4VEmModel(nam),samplingFunction_x(0),samplingFunction_xNoLog(0),
63 fIsInitialised(false),fLocalTable(false) << 53 theMaterial(0),
>> 54 isInitialised(false)
64 { 55 {
65 fIntrinsicLowEnergyLimit = 100.0*eV; 56 fIntrinsicLowEnergyLimit = 100.0*eV;
66 fIntrinsicHighEnergyLimit = 100.0*GeV; 57 fIntrinsicHighEnergyLimit = 100.0*GeV;
67 // SetLowEnergyLimit(fIntrinsicLowEnergyLim << 58 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
68 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 59 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
69 <<
70 if (part) <<
71 SetParticle(part); <<
72 <<
73 // 60 //
74 fVerboseLevel= 0; << 61 verboseLevel= 0;
75 // Verbosity scale: 62 // Verbosity scale:
76 // 0 = nothing << 63 // 0 = nothing
77 // 1 = warning for energy non-conservation << 64 // 1 = warning for energy non-conservation
78 // 2 = details of energy budget 65 // 2 = details of energy budget
79 // 3 = calculation of cross sections, file o 66 // 3 = calculation of cross sections, file openings, sampling of atoms
80 // 4 = entering in methods 67 // 4 = entering in methods
81 68
82 //build the energy grid. It is the same for << 69 PrepareConstants();
83 G4double logenergy = G4Log(fIntrinsicLowEner <<
84 G4double logmaxenergy = G4Log(1.5*fIntrinsic <<
85 //finer grid below 160 keV <<
86 G4double logtransitionenergy = G4Log(160*keV <<
87 G4double logfactor1 = G4Log(10.)/250.; <<
88 G4double logfactor2 = logfactor1*10; <<
89 fLogEnergyGridPMax.push_back(logenergy); <<
90 do{ <<
91 if (logenergy < logtransitionenergy) <<
92 logenergy += logfactor1; <<
93 else <<
94 logenergy += logfactor2; <<
95 fLogEnergyGridPMax.push_back(logenergy); <<
96 }while (logenergy < logmaxenergy); <<
97 } 70 }
98 71
99 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 73
101 G4PenelopeRayleighModel::~G4PenelopeRayleighMo 74 G4PenelopeRayleighModel::~G4PenelopeRayleighModel()
102 { 75 {
103 if (IsMaster() || fLocalTable) << 76 std::map <const G4Material*,G4DataVector*>::iterator i;
104 { << 77 for(i=SamplingTable.begin(); i != SamplingTable.end(); i++) {
105 << 78 delete (*i).second;
106 for(G4int i=0; i<=fMaxZ; ++i) << 79 }
107 { << 80 if (samplingFunction_x) delete samplingFunction_x;
108 if(fLogAtomicCrossSection[i]) << 81 if (samplingFunction_xNoLog) delete samplingFunction_xNoLog;
109 { <<
110 delete fLogAtomicCrossSection[i]; <<
111 fLogAtomicCrossSection[i] = nullptr; <<
112 } <<
113 if(fAtomicFormFactor[i]) <<
114 { <<
115 delete fAtomicFormFactor[i]; <<
116 fAtomicFormFactor[i] = nullptr; <<
117 } <<
118 } <<
119 ClearTables(); <<
120 } <<
121 } <<
122 <<
123 //....oooOO0OOooo........oooOO0OOooo........oo <<
124 void G4PenelopeRayleighModel::ClearTables() <<
125 { <<
126 if (fLogFormFactorTable) <<
127 { <<
128 for (auto& item : (*fLogFormFactorTable <<
129 if (item.second) delete item.second; <<
130 delete fLogFormFactorTable; <<
131 fLogFormFactorTable = nullptr; //zero e <<
132 } <<
133 if (fPMaxTable) <<
134 { <<
135 for (auto& item : (*fPMaxTable)) <<
136 if (item.second) delete item.second; <<
137 delete fPMaxTable; <<
138 fPMaxTable = nullptr; //zero explicitly <<
139 } <<
140 if (fSamplingTable) <<
141 { <<
142 for (auto& item : (*fSamplingTable)) <<
143 if (item.second) delete item.second; <<
144 delete fSamplingTable; <<
145 fSamplingTable = nullptr; //zero explic <<
146 } <<
147 return; <<
148 } 82 }
149 83
150 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 85
152 void G4PenelopeRayleighModel::Initialise(const << 86 void G4PenelopeRayleighModel::Initialise(const G4ParticleDefinition* ,
153 const G4DataVector& ) 87 const G4DataVector& )
154 { 88 {
155 if (fVerboseLevel > 3) << 89 if (verboseLevel > 3)
156 G4cout << "Calling G4PenelopeRayleighModel 90 G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl;
157 91
158 SetParticle(part); << 92 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
159 <<
160 //Only the master model creates/fills/destro <<
161 if (IsMaster() && part == fParticle) <<
162 { 93 {
163 //clear tables depending on materials, n << 94 G4cout << "G4PenelopeRayleighModel: low energy limit increased from " <<
164 ClearTables(); << 95 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
165 << 96 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
166 if (fVerboseLevel > 3) << 97 }
167 G4cout << "Calling G4PenelopeRayleighModel:: << 98 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
168 << 99 {
169 //create new tables << 100 G4cout << "G4PenelopeRayleighModel: high energy limit decreased from " <<
170 if (!fLogFormFactorTable) << 101 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
171 fLogFormFactorTable = new std::map<const G4M << 102 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
172 if (!fPMaxTable) << 103 }
173 fPMaxTable = new std::map<const G4Material*, <<
174 if (!fSamplingTable) <<
175 fSamplingTable = new std::map<const G4Materi <<
176 <<
177 G4ProductionCutsTable* theCoupleTable = <<
178 G4ProductionCutsTable::GetProductionCutsTabl <<
179 104
180 for (G4int i=0;i<(G4int)theCoupleTable-> << 105 G4cout << "Penelope Rayleigh model is initialized " << G4endl
181 { << 106 << "Energy range: "
182 const G4Material* material = << 107 << LowEnergyLimit() / keV << " keV - "
183 theCoupleTable->GetMaterialCutsCouple(i) << 108 << HighEnergyLimit() / GeV << " GeV"
184 const G4ElementVector* theElementVector = << 109 << G4endl;
185 110
186 for (std::size_t j=0;j<material->GetNumber << 111 if(isInitialised) return;
187 { <<
188 G4int iZ = theElementVector->at(j)->Ge <<
189 //read data files only in the master <<
190 if (!fLogAtomicCrossSection[iZ]) <<
191 ReadDataFile(iZ); <<
192 } <<
193 <<
194 //1) If the table has not been built for t <<
195 if (!fLogFormFactorTable->count(material)) <<
196 BuildFormFactorTable(material); <<
197 <<
198 //2) retrieve or build the sampling table <<
199 if (!(fSamplingTable->count(material))) <<
200 InitializeSamplingAlgorithm(material); <<
201 <<
202 //3) retrieve or build the pMax data <<
203 if (!fPMaxTable->count(material)) <<
204 GetPMaxTable(material); <<
205 } <<
206 112
207 if (fVerboseLevel > 1) { << 113 if(pParticleChange)
208 G4cout << "Penelope Rayleigh model v2008 is << 114 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
209 << "Energy range: " << 115 else
210 << LowEnergyLimit() / keV << " keV - << 116 fParticleChange = new G4ParticleChangeForGamma();
211 << HighEnergyLimit() / GeV << " GeV" << 117 isInitialised = true;
212 << G4endl; <<
213 } <<
214 } <<
215 <<
216 if(fIsInitialised) return; <<
217 fParticleChange = GetParticleChangeForGamma( <<
218 fIsInitialised = true; <<
219 } 118 }
220 119
221 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 121
223 void G4PenelopeRayleighModel::InitialiseLocal( << 122 G4double G4PenelopeRayleighModel::CrossSectionPerVolume(const G4Material* material,
224 G4VEmModel *masterModel) << 123 const G4ParticleDefinition* p,
>> 124 G4double ekin,
>> 125 G4double,
>> 126 G4double)
225 { 127 {
226 if (fVerboseLevel > 3) << 128 // Penelope model to calculate the Rayleigh scattering inverse mean
227 G4cout << "Calling G4PenelopeRayleighMode << 129 // free path.
228 // 130 //
229 //Check that particle matches: one might hav << 131 // The basic method is from
230 //for e+ and e-). << 132 // M. Born, Atomic physics, Ed. Blackie and Sons (1969)
>> 133 // using numerical approximations developed in
>> 134 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531
>> 135 // Data used for the form factor used in the calculation (numerical integral of
>> 136 // dSigma/dOmega) are been derived by fitting the atomic forn factor tables
>> 137 // tabulated in
>> 138 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 4 (1975) 471; erratum ibid.
>> 139 // 6 (1977) 615.
>> 140 // The numerical integration of the differential cross section dSigma/dOmega,
>> 141 // which is implemented in the DifferentialCrossSection() method, is performed
>> 142 // by 20-point Gaussian method (managed by G4PenelopeIntegrator).
231 // 143 //
232 if (part == fParticle) <<
233 { <<
234 //Get the const table pointers from the <<
235 const G4PenelopeRayleighModel* theModel <<
236 static_cast<G4PenelopeRayleighModel*> (maste <<
237 <<
238 //Copy pointers to the data tables <<
239 for(G4int i=0; i<=fMaxZ; ++i) <<
240 { <<
241 fLogAtomicCrossSection[i] = theModel->fLog <<
242 fAtomicFormFactor[i] = theModel->fAtomicFo <<
243 } <<
244 fLogFormFactorTable = theModel->fLogForm <<
245 fPMaxTable = theModel->fPMaxTable; <<
246 fSamplingTable = theModel->fSamplingTabl <<
247 <<
248 //copy the G4DataVector with the grid <<
249 fLogQSquareGrid = theModel->fLogQSquareG <<
250 <<
251 //Same verbosity for all workers, as the <<
252 fVerboseLevel = theModel->fVerboseLevel; <<
253 } <<
254 <<
255 return; <<
256 } <<
257 <<
258 <<
259 //....oooOO0OOooo........oooOO0OOooo........oo <<
260 namespace { G4Mutex PenelopeRayleighModelMute <<
261 G4double G4PenelopeRayleighModel::ComputeCross <<
262 G4double energy, <<
263 G4double Z, <<
264 G4double, <<
265 G4double, <<
266 G4double) <<
267 { <<
268 // Cross section of Rayleigh scattering in P <<
269 // tabulation, Cuellen et al. (1997), with n <<
270 // et al. J. Phys. Chem. Ref. Data 4 (1975) <<
271 <<
272 if (fVerboseLevel > 3) <<
273 G4cout << "Calling CrossSectionPerAtom() o <<
274 <<
275 G4int iZ = G4int(Z); <<
276 <<
277 if (!fLogAtomicCrossSection[iZ]) <<
278 { <<
279 //If we are here, it means that Initial <<
280 //not filled up. This can happen in a U <<
281 if (fVerboseLevel > 0) <<
282 { <<
283 //Issue a G4Exception (warning) only in ve <<
284 G4ExceptionDescription ed; <<
285 ed << "Unable to retrieve the cross sectio <<
286 ed << "This can happen only in Unit Tests <<
287 G4Exception("G4PenelopeRayleighModel::Comp <<
288 "em2040",JustWarning,ed); <<
289 } <<
290 //protect file reading via autolock <<
291 G4AutoLock lock(&PenelopeRayleighModelM <<
292 ReadDataFile(iZ); <<
293 lock.unlock(); <<
294 } <<
295 <<
296 G4double cross = 0; <<
297 G4PhysicsFreeVector* atom = fLogAtomicCross <<
298 if (!atom) <<
299 { <<
300 G4ExceptionDescription ed; <<
301 ed << "Unable to find Z=" << iZ << " in <<
302 G4Exception("G4PenelopeRayleighModel::C <<
303 "em2041",FatalException,ed); <<
304 return 0; <<
305 } <<
306 G4double logene = G4Log(energy); <<
307 G4double logXS = atom->Value(logene); <<
308 cross = G4Exp(logXS); <<
309 144
310 if (fVerboseLevel > 2) << 145 if (verboseLevel > 3)
311 { << 146 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeRayleighModel" << G4endl;
312 G4cout << "Rayleigh cross section at " << 147 SetupForMaterial(p, material, ekin);
313 " = " << cross/barn << " barn" << G4endl; << 148
314 } << 149 //Assign local variable "material" to private member "theMaterial", because
315 return cross; << 150 //this information is necessary to calculate the cross section
316 } << 151 theMaterial = material;
317 152
318 << 153 G4int nElements = material->GetNumberOfElements();
319 //....oooOO0OOooo........oooOO0OOooo........oo <<
320 void G4PenelopeRayleighModel::BuildFormFactorT <<
321 { <<
322 /* <<
323 1) get composition and equivalent molecula <<
324 */ <<
325 std::size_t nElements = material->GetNumberO <<
326 const G4ElementVector* elementVector = mater 154 const G4ElementVector* elementVector = material->GetElementVector();
327 const G4double* fractionVector = material->G <<
328 155
329 std::vector<G4double> *StechiometricFactors << 156 G4int maxZ=0;
330 for (std::size_t i=0;i<nElements;++i) << 157 for (G4int i=0; i<nElements; i++)
331 { <<
332 G4double fraction = fractionVector[i]; <<
333 G4double atomicWeigth = (*elementVector) <<
334 StechiometricFactors->push_back(fraction <<
335 } <<
336 //Find max <<
337 G4double MaxStechiometricFactor = 0.; <<
338 for (std::size_t i=0;i<nElements;++i) <<
339 { 158 {
340 if ((*StechiometricFactors)[i] > MaxStec << 159 G4int Z = (G4int) (*elementVector)[i]->GetZ();
341 MaxStechiometricFactor = (*Stechiometr << 160 if (Z>maxZ)
>> 161 maxZ = Z;
342 } 162 }
343 if (MaxStechiometricFactor<1e-16) << 163
>> 164 G4double ec=std::min(ekin,0.5*maxZ);
>> 165 factorE = 849.3315*(ec/electron_mass_c2)*(ec/electron_mass_c2);
>> 166 G4double cs=0;
>> 167 //
>> 168 //Integrate the Differential Cross Section dSigma/dCosTheta between -1 and 1.
>> 169 //
>> 170 G4PenelopeIntegrator<G4PenelopeRayleighModel,G4double(G4PenelopeRayleighModel::*)(G4double)>
>> 171 theIntegrator;
>> 172 cs =
>> 173 theIntegrator.Calculate(this,&G4PenelopeRayleighModel::DifferentialCrossSection,
>> 174 -1.0,0.90,1e-06);
>> 175 cs += theIntegrator.Calculate(this,&G4PenelopeRayleighModel::DifferentialCrossSection,
>> 176 0.90,0.9999999,1e-06);
>> 177 cs = cs*(ec/ekin)*(ec/ekin)*pi*classic_electr_radius*classic_electr_radius;
>> 178 //
>> 179 // Here cs represents the cross section per molecule for materials defined with chemical
>> 180 // formulas and the average cross section per atom in compounds (defined with the mass
>> 181 // fraction)
>> 182 //
>> 183 const G4int* stechiometric = material->GetAtomsVector();
>> 184 //This is the total density of atoms in the material
>> 185 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
>> 186 G4double moleculeDensity = 0;
>> 187
>> 188 //Default case: the material is a compound. In this case, cs is the average cross section
>> 189 // _per_atom_ and one has to multiply for the atom density.
>> 190 G4double cross = atomDensity*cs;
>> 191 G4bool isAMolecule = false;
>> 192
>> 193 //Alternative case: the material is a molecule. In this case cs is the cross section
>> 194 // _per_molecule_ and one has to multiply for the molecule density
>> 195 if (stechiometric)
344 { 196 {
345 G4ExceptionDescription ed; << 197 //Calculate the total number of atoms per molecule
346 ed << "Inconsistent data of atomic compo << 198 G4int atomsPerMolecule = 0;
347 material->GetName() << G4endl; << 199 for (G4int k=0;k<nElements;k++)
348 G4Exception("G4PenelopeRayleighModel::Bu << 200 atomsPerMolecule += stechiometric[k];
349 "em2042",FatalException,ed); << 201 if (atomsPerMolecule)
>> 202 {
>> 203 isAMolecule = true;
>> 204 if (verboseLevel > 3)
>> 205 {
>> 206 G4cout << "Material " << material->GetName() << " is a molecule composed by " <<
>> 207 atomsPerMolecule << " atoms" << G4endl;
>> 208 }
>> 209 moleculeDensity = atomDensity/((G4double) atomsPerMolecule);
>> 210 cross = cs*moleculeDensity;
>> 211 }
350 } 212 }
351 //Normalize <<
352 for (std::size_t i=0;i<nElements;++i) <<
353 (*StechiometricFactors)[i] /= MaxStechiom <<
354 <<
355 /* <<
356 CREATE THE FORM FACTOR TABLE <<
357 */ <<
358 G4PhysicsFreeVector* theFFVec = new G4Physic <<
359 213
360 for (std::size_t k=0;k<fLogQSquareGrid.size( << 214 if (verboseLevel > 2)
361 { 215 {
362 G4double ff2 = 0; //squared form factor << 216 if (isAMolecule)
363 for (std::size_t i=0;i<nElements;++i) <<
364 { 217 {
365 G4int iZ = (*elementVector)[i]->GetZasInt( << 218 G4cout << "Rayleigh cross section at " << ekin/keV << " keV for material " <<
366 G4PhysicsFreeVector* theAtomVec = fAtomicF << 219 material->GetName() << " (molecule) = " << cs/barn << " barn/molecule." << G4endl;
367 G4double f = (*theAtomVec)[k]; //the q-gri << 220 G4cout << "Mean free path: " << (1./cross)/mm << " mm" << G4endl;
368 ff2 += f*f*(*StechiometricFactors)[i]; << 221 }
>> 222 else
>> 223 {
>> 224 G4cout << "Rayleigh cross section at " << ekin/keV << " keV for material " <<
>> 225 material->GetName() << " (compound) = " << cs/barn << " barn/atom." << G4endl;
>> 226 G4cout << "Mean free path: " << (1./cross)/mm << " mm" << G4endl;
369 } 227 }
370 if (ff2) <<
371 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo <<
372 } 228 }
373 theFFVec->FillSecondDerivatives(); //vector << 229 return cross;
374 fLogFormFactorTable->insert(std::make_pair(m <<
375 <<
376 delete StechiometricFactors; <<
377 return; <<
378 } 230 }
379 231
380 //....oooOO0OOooo........oooOO0OOooo........oo 232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 233
382 void G4PenelopeRayleighModel::SampleSecondarie 234 void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
383 const G4MaterialCutsCouple* couple << 235 const G4MaterialCutsCouple* couple,
384 const G4DynamicParticle* aDynamicG << 236 const G4DynamicParticle* aDynamicGamma,
385 G4double, << 237 G4double,
386 G4double) << 238 G4double)
387 { 239 {
388 // Sampling of the Rayleigh final state (nam << 240 // Penelope model to sample the Rayleigh scattering final state.
389 // from the Penelope2008 model. The scatteri << 241 //
390 // cross section dOmega/d(cosTheta) from Bor << 242 // The angular deflection of the scattered photon is sampled according to the
391 // anomalous scattering effects. The Form Fa << 243 // differential cross section dSigma/dOmega used for the numerical integration,
392 // analytical cross section is retrieved via << 244 // and implemented in the DifferentialCrossSection() method. See comments in
393 // are tabulated for F(Q). Form factor for c << 245 // method CrossSectionPerVolume() for more details on the original references
394 // the additivity rule. The sampling from th << 246 // of the model.
395 // Transform with Aliasing (RITA) algorithm; << 247 //
396 // for each material and managed by G4Penelo <<
397 // The sampling algorithm (rejection method) <<
398 // increases with energy. For E=100 keV the <<
399 // hydrogen and uranium, respectively. <<
400 248
401 if (fVerboseLevel > 3) << 249 if (verboseLevel > 3)
402 G4cout << "Calling SamplingSecondaries() o 250 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl;
403 << 251
404 G4double photonEnergy0 = aDynamicGamma->GetK 252 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
405 << 253
406 if (photonEnergy0 <= fIntrinsicLowEnergyLimi << 254 if (photonEnergy0 <= LowEnergyLimit())
407 { << 255 {
408 fParticleChange->ProposeTrackStatus(fSto 256 fParticleChange->ProposeTrackStatus(fStopAndKill);
409 fParticleChange->SetProposedKineticEnerg 257 fParticleChange->SetProposedKineticEnergy(0.);
410 fParticleChange->ProposeLocalEnergyDepos 258 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
411 return ; 259 return ;
412 } << 260 }
413 261
414 G4ParticleMomentum photonDirection0 = aDynam 262 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
415 << 263
416 const G4Material* theMat = couple->GetMateri << 264 // Sampling inizialitation (build internal table)
417 << 265 theMaterial = couple->GetMaterial();
418 //1) Verify if tables are ready << 266 InitialiseSampling();
419 //Either Initialize() was not called, or we << 267
420 //not invoked << 268 G4DataVector* samplingFunction_y = SamplingTable.find(theMaterial)->second;
421 if (!fPMaxTable || !fSamplingTable || !fLogF << 269
422 { << 270 // Sample the angle of the scattered photon
423 //create a **thread-local** version of t << 271 G4double x2max = 2.0*std::log(41.2148*photonEnergy0/electron_mass_c2);
424 //Unit Tests << 272 G4int jm=0;
425 fLocalTable = true; << 273 G4int asize = samplingFunction_x->size();
426 if (!fLogFormFactorTable) << 274 if (x2max < (*samplingFunction_x)[1])
427 fLogFormFactorTable = new std::map<const G4M << 275 jm=0;
428 if (!fPMaxTable) << 276 else if(x2max>(*samplingFunction_x)[asize-2])
429 fPMaxTable = new std::map<const G4Material*, << 277 jm=asize-2;
430 if (!fSamplingTable) << 278 else
431 fSamplingTable = new std::map<const G4Materi << 279 {
432 } << 280 //spacing in the logTable
433 << 281 G4double logScalingFactor = (*samplingFunction_x)[1]-(*samplingFunction_x)[0];
434 if (!fSamplingTable->count(theMat)) << 282 jm=(G4int) ((x2max-(*samplingFunction_x)[0])/logScalingFactor);
435 { << 283 }
436 //If we are here, it means that Initiali << 284
437 //not filled up. This can happen in a Un << 285 G4double rumax = (*samplingFunction_y)[jm]+((*samplingFunction_y)[jm+1]-(*samplingFunction_y)[jm])*
438 if (fVerboseLevel > 0) << 286 (x2max-(*samplingFunction_x)[jm])/((*samplingFunction_x)[jm+1]-(*samplingFunction_x)[jm]);
439 { << 287 G4double cosTheta=0;
440 //Issue a G4Exception (warning) only in ve << 288 G4double rejectionValue = 0;
441 G4ExceptionDescription ed; << 289 do{
442 ed << "Unable to find the fSamplingTable d << 290 G4double ru = rumax + std::log(G4UniformRand());
443 theMat->GetName() << G4endl; << 291 G4int j=0;
444 ed << "This can happen only in Unit Tests" << 292 G4int ju=jm+1;
445 G4Exception("G4PenelopeRayleighModel::Samp << 293 do{
446 "em2019",JustWarning,ed); << 294 G4int jt=(j+ju)/2; //bipartition
447 } << 295 if (ru > (*samplingFunction_y)[jt])
448 const G4ElementVector* theElementVector << 296 j=jt;
449 //protect file reading via autolock << 297 else
450 G4AutoLock lock(&PenelopeRayleighModelMu << 298 ju=jt;
451 for (std::size_t j=0;j<theMat->GetNumber << 299 }while ((ju-j)>1);
452 { << 300 G4double x2rat = 0;
453 G4int iZ = theElementVector->at(j)->GetZas << 301 G4double denomin = (*samplingFunction_y)[j+1]-(*samplingFunction_y)[j];
454 if (!fLogAtomicCrossSection[iZ]) << 302 if (denomin > 1e-12)
455 { << 303 {
456 lock.lock(); << 304 x2rat = (*samplingFunction_x)[j]+(((*samplingFunction_x)[j+1]-(*samplingFunction_x)[j])*
457 ReadDataFile(iZ); << 305 (ru-(*samplingFunction_y)[j])/denomin)-x2max;
458 lock.unlock(); << 306 }
459 } << 307 else
460 } << 308 {
461 lock.lock(); << 309 x2rat = (*samplingFunction_x)[j]-x2max;
462 //1) If the table has not been built for << 310 }
463 if (!fLogFormFactorTable->count(theMat)) << 311 cosTheta = 1.0-2.0*std::exp(x2rat);
464 BuildFormFactorTable(theMat); << 312 rejectionValue = 0.5*(1.0+cosTheta*cosTheta);
465 << 313 }while (G4UniformRand() > rejectionValue);
466 //2) retrieve or build the sampling tabl <<
467 if (!(fSamplingTable->count(theMat))) <<
468 InitializeSamplingAlgorithm(theMat); <<
469 <<
470 //3) retrieve or build the pMax data <<
471 if (!fPMaxTable->count(theMat)) <<
472 GetPMaxTable(theMat); <<
473 lock.unlock(); <<
474 } <<
475 <<
476 //Ok, restart the job <<
477 G4PenelopeSamplingData* theDataTable = fSamp <<
478 G4PhysicsFreeVector* thePMax = fPMaxTable->f <<
479 <<
480 G4double cosTheta = 1.0; <<
481 <<
482 //OK, ready to go! <<
483 G4double qmax = 2.0*photonEnergy0/electron_m <<
484 <<
485 if (qmax < 1e-10) //very low momentum transf <<
486 { <<
487 G4bool loopAgain=false; <<
488 do <<
489 { <<
490 loopAgain = false; <<
491 cosTheta = 1.0-2.0*G4UniformRand(); <<
492 G4double G = 0.5*(1+cosTheta*cosTheta); <<
493 if (G4UniformRand()>G) <<
494 loopAgain = true; <<
495 }while(loopAgain); <<
496 } <<
497 else //larger momentum transfer <<
498 { <<
499 std::size_t nData = theDataTable->GetNum <<
500 G4double LastQ2inTheTable = theDataTable <<
501 G4double q2max = std::min(qmax*qmax,Last <<
502 <<
503 G4bool loopAgain = false; <<
504 G4double MaxPValue = thePMax->Value(phot <<
505 G4double xx=0; <<
506 <<
507 //Sampling by rejection method. The reje <<
508 //G = 0.5*(1+cos^2(theta)) <<
509 // <<
510 do{ <<
511 loopAgain = false; <<
512 G4double RandomMax = G4UniformRand()*MaxPVal <<
513 xx = theDataTable->SampleValue(RandomMax); <<
514 //xx is a random value of q^2 in (0,q2max),s <<
515 //F(Q^2) via the RITA algorithm <<
516 if (xx > q2max) <<
517 loopAgain = true; <<
518 cosTheta = 1.0-2.0*xx/q2max; <<
519 G4double G = 0.5*(1+cosTheta*cosTheta); <<
520 if (G4UniformRand()>G) <<
521 loopAgain = true; <<
522 }while(loopAgain); <<
523 } <<
524 314
525 G4double sinTheta = std::sqrt(1-cosTheta*cos 315 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
526 << 316
527 // Scattered photon angles. ( Z - axis along 317 // Scattered photon angles. ( Z - axis along the parent photon)
528 G4double phi = twopi * G4UniformRand() ; 318 G4double phi = twopi * G4UniformRand() ;
529 G4double dirX = sinTheta*std::cos(phi); 319 G4double dirX = sinTheta*std::cos(phi);
530 G4double dirY = sinTheta*std::sin(phi); 320 G4double dirY = sinTheta*std::sin(phi);
531 G4double dirZ = cosTheta; 321 G4double dirZ = cosTheta;
532 << 322
533 // Update G4VParticleChange for the scattere << 323 // Update G4VParticleChange for the scattered photon
534 G4ThreeVector photonDirection1(dirX, dirY, d 324 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
535 325
536 photonDirection1.rotateUz(photonDirection0); 326 photonDirection1.rotateUz(photonDirection0);
537 fParticleChange->ProposeMomentumDirection(ph 327 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
538 fParticleChange->SetProposedKineticEnergy(ph 328 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
539 << 329
540 return; 330 return;
541 } 331 }
542 332
543 <<
544 //....oooOO0OOooo........oooOO0OOooo........oo 333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
545 334
546 void G4PenelopeRayleighModel::ReadDataFile(con << 335 G4double G4PenelopeRayleighModel::DifferentialCrossSection(G4double cosTheta)
547 { 336 {
548 if (fVerboseLevel > 2) << 337 //Differential cross section for Rayleigh scattering
549 { << 338 G4double x2 = factorE*(1-cosTheta);
550 G4cout << "G4PenelopeRayleighModel::Read << 339 G4double gradx = (1+cosTheta*cosTheta)*MolecularFormFactor(x2);
551 G4cout << "Going to read Rayleigh data f << 340 return gradx;
552 } << 341 }
553 const char* path = G4FindDataDir("G4LEDATA << 342
554 if(!path) << 343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
555 { << 344
556 G4String excep = "G4LEDATA environment v << 345 G4double G4PenelopeRayleighModel::MolecularFormFactor(G4double x2)
557 G4Exception("G4PenelopeRayleighModel::Re << 346 {
558 "em0006",FatalException,excep); << 347 //Squared molecular form factor (additivity rule)
559 return; << 348 const G4int ntot=95;
560 } << 349 G4double RA1[ntot] = {0.0e0, 3.9265e+0, 4.3100e+1, 5.2757e+1, 2.5021e+1,
561 << 350 1.2211e+1, 9.3229e+0, 3.2455e+0, 2.4197e+0, 1.5985e+0,
562 /* << 351 3.0926e+1, 1.5315e+1, 7.7061e+0, 3.9493e+0, 2.2042e+0,
563 Read first the cross section file << 352 1.9453e+1, 1.9354e+1, 8.0374e+0, 8.3779e+1, 5.7370e+1,
564 */ << 353 5.2310e+1, 4.7514e+1, 4.3785e+1, 4.2048e+1, 3.6972e+1,
565 std::ostringstream ost; << 354 3.3849e+1, 3.1609e+1, 2.8763e+1, 2.7217e+1, 2.4263e+1,
566 if (Z>9) << 355 2.2403e+1, 1.8606e+1, 1.5143e+1, 1.4226e+1, 1.1792e+1,
567 ost << path << "/penelope/rayleigh/pdgra" << 356 9.7574e+0, 1.2796e+1, 1.2854e+1, 1.2368e+1, 1.0208e+1,
568 else << 357 8.2823e+0, 7.4677e+0, 7.6028e+0, 6.1090e+0, 5.5346e+0,
569 ost << path << "/penelope/rayleigh/pdgra0" << 358 4.2340e+0, 4.0444e+0, 4.2905e+0, 4.7950e+0, 5.1112e+0,
570 std::ifstream file(ost.str().c_str()); << 359 5.2407e+0, 5.2153e+0, 5.1639e+0, 4.8814e+0, 5.8054e+0,
571 if (!file.is_open()) << 360 6.6724e+0, 6.5104e+0, 6.3364e+0, 6.2889e+0, 6.3028e+0,
572 { << 361 6.3853e+0, 6.3475e+0, 6.5779e+0, 6.8486e+0, 7.0993e+0,
573 G4String excep = "Data file " + G4String << 362 7.6122e+0, 7.9681e+0, 8.3481e+0, 6.3875e+0, 8.0042e+0,
574 G4Exception("G4PenelopeRayleighModel::Re << 363 8.0820e+0, 7.6940e+0, 7.1927e+0, 6.6751e+0, 6.1623e+0,
575 "em0003",FatalException,excep); << 364 5.8335e+0, 5.5599e+0, 4.6551e+0, 4.4327e+0, 4.7601e+0,
576 } << 365 5.2872e+0, 5.6084e+0, 5.7680e+0, 5.8041e+0, 5.7566e+0,
577 G4int readZ =0; << 366 5.6541e+0, 6.3932e+0, 6.9313e+0, 7.0027e+0, 6.8796e+0,
578 std::size_t nPoints= 0; << 367 6.4739e+0, 6.2405e+0, 6.0081e+0, 5.5708e+0, 5.3680e+0};
579 file >> readZ >> nPoints; << 368
580 //check the right file is opened. << 369
581 if (readZ != Z || nPoints <= 0 || nPoints >= << 370 G4double RA2[ntot] = {0.0e0, 1.3426e-1, 9.4875e+1,-1.0896e+2,-4.5494e+1,
582 { << 371 -1.9572e+1,-1.2382e+1,-3.6827e+0,-2.4542e+0,-1.4453e+0,
583 G4ExceptionDescription ed; << 372 1.3401e+2, 7.9717e+1, 6.2164e+1, 4.0300e+1, 3.1682e+1,
584 ed << "Corrupted data file for Z=" << Z << 373 -1.3639e+1,-1.5950e+1,-5.1523e+0, 1.8351e+2, 1.2205e+2,
585 G4Exception("G4PenelopeRayleighModel::Re << 374 1.0007e+2, 8.5632e+1, 7.9145e+1, 6.3675e+1, 6.2954e+1,
586 "em0005",FatalException,ed); << 375 5.6601e+1, 5.4171e+1, 4.8752e+1, 3.8062e+1, 3.9933e+1,
587 return; << 376 4.8343e+1, 4.2137e+1, 3.4617e+1, 2.9430e+1, 2.4010e+1,
588 } << 377 1.9744e+1, 4.0009e+1, 5.1614e+1, 5.0456e+1, 3.9088e+1,
589 << 378 2.6824e+1, 2.2953e+1, 2.4773e+1, 1.6893e+1, 1.4548e+1,
590 fLogAtomicCrossSection[Z] = new G4PhysicsFre << 379 9.7226e+0, 1.0192e+1, 1.1153e+1, 1.3188e+1, 1.4733e+1,
591 G4double ene=0,f1=0,f2=0,xs=0; << 380 1.5644e+1, 1.5939e+1, 1.5923e+1, 1.5254e+1, 2.0748e+1,
592 for (std::size_t i=0;i<nPoints;++i) << 381 2.6901e+1, 2.7032e+1, 2.4938e+1, 2.1528e+1, 2.0362e+1,
593 { << 382 1.9474e+1, 1.8238e+1, 1.7898e+1, 1.9174e+1, 1.9023e+1,
594 file >> ene >> f1 >> f2 >> xs; << 383 1.8194e+1, 1.8504e+1, 1.8955e+1, 1.4276e+1, 1.7558e+1,
595 //dimensional quantities << 384 1.8651e+1, 1.7984e+1, 1.6793e+1, 1.5469e+1, 1.4143e+1,
596 ene *= eV; << 385 1.3149e+1, 1.2255e+1, 9.2352e+0, 8.6067e+0, 9.7460e+0,
597 xs *= cm2; << 386 1.1749e+1, 1.3281e+1, 1.4326e+1, 1.4920e+1, 1.5157e+1,
598 fLogAtomicCrossSection[Z]->PutValue(i,G4 << 387 1.5131e+1, 1.9489e+1, 2.3649e+1, 2.4686e+1, 2.4760e+1,
599 if (file.eof() && i != (nPoints-1)) //fi << 388 2.1519e+1, 2.0099e+1, 1.8746e+1, 1.5943e+1, 1.4880e+1};
600 { << 389
601 G4ExceptionDescription ed ; << 390 G4double RA3[ntot] = {0.0e0, 2.2648e+0, 1.0579e+3, 8.6177e+2, 2.4422e+2,
602 ed << "Corrupted data file for Z=" << Z << << 391 7.8788e+1, 3.8293e+1, 1.2564e+1, 6.9091e+0, 3.7926e+0,
603 ed << "Found less than " << nPoints << "en << 392 0.0000e+0, 0.0000e+0, 1.6759e-9, 1.3026e+1, 3.0569e+0,
604 G4Exception("G4PenelopeRayleighModel::Read << 393 1.5521e+2, 1.2815e+2, 4.7378e+1, 9.2802e+2, 4.7508e+2,
605 "em0005",FatalException,ed); << 394 3.6612e+2, 2.7582e+2, 2.1008e+2, 1.5903e+2, 1.2322e+2,
606 } << 395 9.2898e+1, 7.1345e+1, 5.1651e+1, 3.8474e+1, 2.7410e+1,
607 } << 396 1.9126e+1, 1.0889e+1, 5.3479e+0, 8.2223e+0, 5.0837e+0,
608 file.close(); << 397 2.8905e+0, 2.7457e+0, 6.7082e-1, 0.0000e+0, 0.0000e+0,
609 << 398 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0,
610 /* << 399 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0,
611 Then read the form factor file << 400 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0,
612 */ << 401 0.0000e+0, 0.0000e+0, 0.0000e+0, 1.7264e-1, 2.7322e-1,
613 std::ostringstream ost2; << 402 3.9444e-1, 4.5648e-1, 6.2286e-1, 7.2468e-1, 8.4296e-1,
614 if (Z>9) << 403 1.1698e+0, 1.2994e+0, 1.4295e+0, 0.0000e+0, 8.1570e-1,
615 ost2 << path << "/penelope/rayleigh/pdaff" << 404 6.9349e-1, 4.9536e-1, 3.1211e-1, 1.5931e-1, 2.9512e-2,
616 else << 405 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0,
617 ost2 << path << "/penelope/rayleigh/pdaff0 << 406 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0,
618 file.open(ost2.str().c_str()); << 407 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0,
619 if (!file.is_open()) << 408 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0, 0.0000e+0};
620 { << 409
621 G4String excep = "Data file " + G4String << 410 G4double RA4[ntot] = {1.1055e1,6.3519e0,4.7367e+1, 3.9402e+1, 2.2896e+1,
622 G4Exception("G4PenelopeRayleighModel::Re << 411 1.3979e+1, 1.0766e+1, 6.5252e+0, 5.1631e+0, 4.0524e+0,
623 "em0003",FatalException,excep); << 412 2.7145e+1, 1.8724e+1, 1.4782e+1, 1.1608e+1, 9.7750e+0,
624 } << 413 1.6170e+1, 1.5249e+1, 9.1916e+0, 5.4499e+1, 4.1381e+1,
625 file >> readZ >> nPoints; << 414 3.7395e+1, 3.3815e+1, 3.1135e+1, 2.8273e+1, 2.6140e+1,
626 //check the right file is opened. << 415 2.3948e+1, 2.2406e+1, 2.0484e+1, 1.8453e+1, 1.7386e+1,
627 if (readZ != Z || nPoints <= 0 || nPoints >= << 416 1.7301e+1, 1.5388e+1, 1.3411e+1, 1.2668e+1, 1.1133e+1,
628 { << 417 9.8081e+0, 1.3031e+1, 1.4143e+1, 1.3815e+1, 1.2077e+1,
629 G4ExceptionDescription ed; << 418 1.0033e+1, 9.2549e+0, 9.5338e+0, 7.9076e+0, 7.3263e+0,
630 ed << "Corrupted data file for Z=" << Z << 419 5.9996e+0, 6.0087e+0, 6.2660e+0, 6.7914e+0, 7.1501e+0,
631 G4Exception("G4PenelopeRayleighModel::Re << 420 7.3367e+0, 7.3729e+0, 7.3508e+0, 7.1465e+0, 8.2731e+0,
632 "em0005",FatalException,ed); << 421 9.3745e+0, 9.3508e+0, 8.9897e+0, 8.4566e+0, 8.2690e+0,
633 return; << 422 8.1398e+0, 7.9183e+0, 7.9123e+0, 8.1677e+0, 8.1871e+0,
634 } << 423 8.1766e+0, 8.2881e+0, 8.4227e+0, 7.0273e+0, 8.0002e+0,
635 fAtomicFormFactor[Z] = new G4PhysicsFreeVect << 424 8.1440e+0, 7.9104e+0, 7.5685e+0, 7.1970e+0, 6.8184e+0,
636 G4double q=0,ff=0,incoh=0; << 425 6.5469e+0, 6.3056e+0, 5.4844e+0, 5.2832e+0, 5.5889e+0,
637 G4bool fillQGrid = false; << 426 6.0919e+0, 6.4340e+0, 6.6426e+0, 6.7428e+0, 6.7636e+0,
638 //fill this vector only the first time. << 427 6.7281e+0, 7.5729e+0, 8.2808e+0, 8.4400e+0, 8.4220e+0,
639 if (!fLogQSquareGrid.size()) << 428 7.8662e+0, 7.5993e+0, 7.3353e+0, 6.7829e+0, 6.5520e+0};
640 fillQGrid = true; << 429
641 for (std::size_t i=0;i<nPoints;++i) << 430 G4double RA5[ntot] = {0.0e0, 4.9828e+0, 5.5674e+1, 3.0902e+1, 1.1496e+1,
642 { << 431 4.8936e+0, 2.5506e+0, 1.2236e+0, 7.4698e-1, 4.7042e-1,
643 file >> q >> ff >> incoh; << 432 4.7809e+0, 4.6315e+0, 4.3677e+0, 4.9269e+0, 2.6033e+0,
644 //q and ff are dimensionless (q is in un << 433 9.6229e+0, 7.2592e+0, 4.1634e+0, 1.3999e+1, 8.6975e+0,
645 fAtomicFormFactor[Z]->PutValue(i,q,ff); << 434 6.9630e+0, 5.4681e+0, 4.2653e+0, 3.2848e+0, 2.7354e+0,
646 if (fillQGrid) << 435 2.1617e+0, 1.7030e+0, 1.2826e+0, 9.7080e-1, 7.2227e-1,
647 { << 436 5.0874e-1, 3.1402e-1, 1.6360e-1, 3.2918e-1, 2.3570e-1,
648 fLogQSquareGrid.push_back(2.0*G4Log(q)); << 437 1.5868e-1, 1.5146e-1, 9.7662e-2, 7.3151e-2, 6.4206e-2,
649 } << 438 4.8945e-2, 4.3189e-2, 4.4368e-2, 3.3976e-2, 3.0466e-2,
650 if (file.eof() && i != (nPoints-1)) //fi << 439 2.4477e-2, 3.7202e-2, 3.7093e-2, 3.8161e-2, 3.8576e-2,
651 { << 440 3.8403e-2, 3.7806e-2, 3.4958e-2, 3.6029e-2, 4.3087e-2,
652 G4ExceptionDescription ed; << 441 4.7069e-2, 4.6452e-2, 4.2486e-2, 4.1517e-2, 4.1691e-2,
653 ed << "Corrupted data file for Z=" << Z << << 442 4.2813e-2, 4.2294e-2, 4.5287e-2, 4.8462e-2, 4.9726e-2,
654 ed << "Found less than " << nPoints << "en << 443 5.5097e-2, 5.6568e-2, 5.8069e-2, 1.2270e-2, 3.8006e-2,
655 G4Exception("G4PenelopeRayleighModel::Read << 444 3.5048e-2, 3.0050e-2, 2.5069e-2, 2.0485e-2, 1.6151e-2,
656 "em0005",FatalException,ed); << 445 1.4631e-2, 1.4034e-2, 1.1978e-2, 1.1522e-2, 1.2375e-2,
657 } << 446 1.3805e-2, 1.4954e-2, 1.5832e-2, 1.6467e-2, 1.6896e-2,
658 } << 447 1.7166e-2, 1.9954e-2, 2.2497e-2, 2.1942e-2, 2.1965e-2,
659 file.close(); << 448 2.0005e-2, 1.8927e-2, 1.8167e-2, 1.6314e-2, 1.5522e-2};
660 return; << 449
661 } << 450 G4double x=std::sqrt(x2);
662 << 451 G4double gradx1=0.0;
663 //....oooOO0OOooo........oooOO0OOooo........oo << 452
664 << 453 G4int nElements = theMaterial->GetNumberOfElements();
665 G4double G4PenelopeRayleighModel::GetFSquared( << 454 const G4ElementVector* elementVector = theMaterial->GetElementVector();
666 { << 455 const G4int* stechiometric = theMaterial->GetAtomsVector();
667 G4double f2 = 0; << 456 const G4double* vector_of_atoms = theMaterial->GetVecNbOfAtomsPerVolume();
668 //Input value QSquared could be zero: protec << 457 const G4double tot_atoms = theMaterial->GetTotNbOfAtomsPerVolume();
669 //the FPE exception << 458 for (G4int i=0;i<nElements;i++)
670 //If Q<1e-10, set Q to 1e-10 << 459 {
671 G4double logQSquared = (QSquared>1e-10) ? G4 << 460 G4int Z = (G4int) (*elementVector)[i]->GetZ();
672 //last value of the table << 461 if (Z>ntot) Z=95;
673 G4double maxlogQ2 = fLogQSquareGrid[fLogQSqu << 462 G4double denomin = 1.+x2*(RA4[Z-1]+x2*RA5[Z-1]);
674 << 463 G4double fa=Z*(1+x2*(RA1[Z-1]+x*(RA2[Z-1]+x*RA3[Z-1])))/(denomin*denomin);
675 //now it should be all right << 464 if (Z>10 && fa>2.0)
676 G4PhysicsFreeVector* theVec = fLogFormFactor << 465 {
677 << 466 G4double k1=0.3125;
678 if (!theVec) << 467 G4double k2=2.426311e-02;
679 { << 468 G4double Pa=(Z-k1)*fine_structure_const;
680 G4ExceptionDescription ed; << 469 G4double Pg=std::sqrt(1-(Pa*Pa));
681 ed << "Unable to retrieve F squared tabl << 470 G4double Pq=k2*x/Pa;
682 G4Exception("G4PenelopeRayleighModel::Ge << 471 G4double fb=std::sin(2.0*Pg*std::atan(Pq))/(Pg*Pq*std::pow((1+Pq*Pq),Pg));
683 "em2046",FatalException,ed); << 472 fa=std::max(fa,fb);
684 return 0; <<
685 } <<
686 if (logQSquared < -20) // Q < 1e-9 <<
687 { <<
688 G4double logf2 = (*theVec)[0]; //first v <<
689 f2 = G4Exp(logf2); <<
690 } <<
691 else if (logQSquared > maxlogQ2) <<
692 f2 =0; <<
693 else <<
694 { <<
695 //log(Q^2) vs. log(F^2) <<
696 G4double logf2 = theVec->Value(logQSquar <<
697 f2 = G4Exp(logf2); <<
698 <<
699 } <<
700 if (fVerboseLevel > 3) <<
701 { <<
702 G4cout << "G4PenelopeRayleighModel::GetF <<
703 G4cout << "Q^2 = " << QSquared << " (un <<
704 } <<
705 return f2; <<
706 } <<
707 <<
708 //....oooOO0OOooo........oooOO0OOooo........oo <<
709 <<
710 void G4PenelopeRayleighModel::InitializeSampli <<
711 { <<
712 G4double q2min = 0; <<
713 G4double q2max = 0; <<
714 const std::size_t np = 150; //hard-coded in <<
715 //G4cout << "Init N= " << fLogQSquareGrid.si <<
716 for (std::size_t i=1;i<fLogQSquareGrid.size( <<
717 { <<
718 G4double Q2 = G4Exp(fLogQSquareGrid[i]); <<
719 if (GetFSquared(mat,Q2) > 1e-35) <<
720 { <<
721 q2max = G4Exp(fLogQSquareGrid[i-1]); <<
722 } <<
723 //G4cout << "Q2= " << Q2 << " q2max= " < <<
724 } <<
725 <<
726 std::size_t nReducedPoints = np/4; <<
727 <<
728 //check for errors <<
729 if (np < 16) <<
730 { <<
731 G4Exception("G4PenelopeRayleighModel::In <<
732 "em2047",FatalException, <<
733 "Too few points to initialize the sampli <<
734 } <<
735 if (q2min > (q2max-1e-10)) <<
736 { <<
737 G4cout << "q2min= " << q2min << " q2max= <<
738 G4Exception("G4PenelopeRayleighModel::In <<
739 "em2048",FatalException, <<
740 "Too narrow grid to initialize the sampl <<
741 } <<
742 <<
743 //This is subroutine RITAI0 of Penelope <<
744 //Create an object of type G4PenelopeRayleig <<
745 <<
746 //temporary vectors --> Then everything is s <<
747 G4DataVector* x = new G4DataVector(); <<
748 <<
749 /******************************************* <<
750 Start with a grid of NUNIF points uniforml <<
751 ******************************************** <<
752 std::size_t NUNIF = std::min(std::max(((std: <<
753 const G4int nip = 51; //hard-coded in Penelo <<
754 <<
755 G4double dx = (q2max-q2min)/((G4double) NUNI <<
756 x->push_back(q2min); <<
757 for (std::size_t i=1;i<NUNIF-1;++i) <<
758 { <<
759 G4double app = q2min + i*dx; <<
760 x->push_back(app); //increase <<
761 } <<
762 x->push_back(q2max); <<
763 <<
764 if (fVerboseLevel> 3) <<
765 G4cout << "Vector x has " << x->size() << <<
766 <<
767 //Allocate temporary storage vectors <<
768 G4DataVector* area = new G4DataVector(); <<
769 G4DataVector* a = new G4DataVector(); <<
770 G4DataVector* b = new G4DataVector(); <<
771 G4DataVector* c = new G4DataVector(); <<
772 G4DataVector* err = new G4DataVector(); <<
773 <<
774 for (std::size_t i=0;i<NUNIF-1;++i) //build <<
775 { <<
776 //Temporary vectors for this loop <<
777 G4DataVector* pdfi = new G4DataVector(); <<
778 G4DataVector* pdfih = new G4DataVector() <<
779 G4DataVector* sumi = new G4DataVector(); <<
780 <<
781 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4do <<
782 G4double pdfmax = 0; <<
783 for (G4int k=0;k<nip;k++) <<
784 { <<
785 G4double xik = (*x)[i]+k*dxi; <<
786 G4double pdfk = std::max(GetFSquared(mat,x <<
787 pdfi->push_back(pdfk); <<
788 pdfmax = std::max(pdfmax,pdfk); <<
789 if (k < (nip-1)) <<
790 { <<
791 G4double xih = xik + 0.5*dxi; <<
792 G4double pdfIK = std::max(GetFSquared( <<
793 pdfih->push_back(pdfIK); <<
794 pdfmax = std::max(pdfmax,pdfIK); <<
795 } <<
796 } <<
797 <<
798 //Simpson's integration <<
799 G4double cons = dxi*0.5*(1./3.); <<
800 sumi->push_back(0.); <<
801 for (G4int k=1;k<nip;k++) <<
802 { <<
803 G4double previous = (*sumi)[k-1]; <<
804 G4double next = previous + cons*((*pdfi)[k <<
805 sumi->push_back(next); <<
806 } <<
807 <<
808 G4double lastIntegral = (*sumi)[sumi->si <<
809 area->push_back(lastIntegral); <<
810 //Normalize cumulative function <<
811 G4double factor = 1.0/lastIntegral; <<
812 for (std::size_t k=0;k<sumi->size();++k) <<
813 (*sumi)[k] *= factor; <<
814 <<
815 //When the PDF vanishes at one of the in <<
816 if ((*pdfi)[0] < 1e-35) <<
817 (*pdfi)[0] = 1e-5*pdfmax; <<
818 if ((*pdfi)[pdfi->size()-1] < 1e-35) <<
819 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; <<
820 <<
821 G4double pli = (*pdfi)[0]*factor; <<
822 G4double pui = (*pdfi)[pdfi->size()-1]*f <<
823 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx <<
824 G4double A_temp = (1.0/(pli*dx))-1.0-B_t <<
825 G4double C_temp = 1.0+A_temp+B_temp; <<
826 if (C_temp < 1e-35) <<
827 { <<
828 a->push_back(0.); <<
829 b->push_back(0.); <<
830 c->push_back(1.); <<
831 } 473 }
>> 474 if (stechiometric && stechiometric[i]!=0)
>> 475 gradx1 += stechiometric[i]*(fa*fa); //sum on the molecule
832 else 476 else
833 { << 477 gradx1 += (vector_of_atoms[i]/tot_atoms)*(fa*fa); //weighted mean
834 a->push_back(A_temp); <<
835 b->push_back(B_temp); <<
836 c->push_back(C_temp); <<
837 } <<
838 <<
839 //OK, now get ERR(I), the integral of th <<
840 //and the true pdf, extended over the in <<
841 G4int icase = 1; //loop code <<
842 G4bool reLoop = false; <<
843 err->push_back(0.); <<
844 do <<
845 { <<
846 reLoop = false; <<
847 (*err)[i] = 0.; //zero variable <<
848 for (G4int k=0;k<nip;k++) <<
849 { <<
850 G4double rr = (*sumi)[k]; <<
851 G4double pap = (*area)[i]*(1.0+((*a)[i <<
852 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(* <<
853 if (k == 0 || k == nip-1) <<
854 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]) <<
855 else <<
856 (*err)[i] += std::fabs(pap-(*pdfi)[k]); <<
857 } <<
858 (*err)[i] *= dxi; <<
859 <<
860 //If err(I) is too large, the pdf is appro <<
861 if ((*err)[i] > 0.1*(*area)[i] && icase == <<
862 { <<
863 (*b)[i] = 0; <<
864 (*a)[i] = 0; <<
865 (*c)[i] = 1.; <<
866 icase = 2; <<
867 reLoop = true; <<
868 } <<
869 }while(reLoop); <<
870 delete pdfi; <<
871 delete pdfih; <<
872 delete sumi; <<
873 } //end of first loop over i <<
874 <<
875 //Now assign last point <<
876 (*x)[x->size()-1] = q2max; <<
877 a->push_back(0.); <<
878 b->push_back(0.); <<
879 c->push_back(0.); <<
880 err->push_back(0.); <<
881 area->push_back(0.); <<
882 <<
883 if (x->size() != NUNIF || a->size() != NUNIF <<
884 err->size() != NUNIF || area->size() != <<
885 { <<
886 G4ExceptionDescription ed; <<
887 ed << "Problem in building the Table for <<
888 G4Exception("G4PenelopeRayleighModel::In <<
889 "em2049",FatalException,ed); <<
890 } 478 }
891 << 479 return gradx1;
892 /******************************************* <<
893 New grid points are added by halving the su <<
894 This is done up to np=150 points in the grid <<
895 ******************************************** <<
896 do <<
897 { <<
898 G4double maxError = 0.0; <<
899 std::size_t iErrMax = 0; <<
900 for (std::size_t i=0;i<err->size()-2;++i <<
901 { <<
902 //maxError is the lagest of the interval e <<
903 if ((*err)[i] > maxError) <<
904 { <<
905 maxError = (*err)[i]; <<
906 iErrMax = i; <<
907 } <<
908 } <<
909 <<
910 //OK, now I have to insert one new point <<
911 G4double newx = 0.5*((*x)[iErrMax]+(*x)[ <<
912 <<
913 x->insert(x->begin()+iErrMax+1,newx); <<
914 //Add place-holders in the other vectors <<
915 area->insert(area->begin()+iErrMax+1,0.) <<
916 a->insert(a->begin()+iErrMax+1,0.); <<
917 b->insert(b->begin()+iErrMax+1,0.); <<
918 c->insert(c->begin()+iErrMax+1,0.); <<
919 err->insert(err->begin()+iErrMax+1,0.); <<
920 <<
921 //Now calculate the other parameters <<
922 for (std::size_t i=iErrMax;i<=iErrMax+1; <<
923 { <<
924 //Temporary vectors for this loop <<
925 G4DataVector* pdfi = new G4DataVector(); <<
926 G4DataVector* pdfih = new G4DataVector(); <<
927 G4DataVector* sumi = new G4DataVector(); <<
928 <<
929 G4double dxLocal = (*x)[i+1]-(*x)[i]; <<
930 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4doub <<
931 G4double pdfmax = 0; <<
932 for (G4int k=0;k<nip;k++) <<
933 { <<
934 G4double xik = (*x)[i]+k*dxi; <<
935 G4double pdfk = std::max(GetFSquared(m <<
936 pdfi->push_back(pdfk); <<
937 pdfmax = std::max(pdfmax,pdfk); <<
938 if (k < (nip-1)) <<
939 { <<
940 G4double xih = xik + 0.5*dxi; <<
941 G4double pdfIK = std::max(GetFSquared(ma <<
942 pdfih->push_back(pdfIK); <<
943 pdfmax = std::max(pdfmax,pdfIK); <<
944 } <<
945 } <<
946 <<
947 //Simpson's integration <<
948 G4double cons = dxi*0.5*(1./3.); <<
949 sumi->push_back(0.); <<
950 for (G4int k=1;k<nip;k++) <<
951 { <<
952 G4double previous = (*sumi)[k-1]; <<
953 G4double next = previous + cons*((*pdf <<
954 sumi->push_back(next); <<
955 } <<
956 G4double lastIntegral = (*sumi)[sumi->size <<
957 (*area)[i] = lastIntegral; <<
958 <<
959 //Normalize cumulative function <<
960 G4double factor = 1.0/lastIntegral; <<
961 for (std::size_t k=0;k<sumi->size();++k) <<
962 (*sumi)[k] *= factor; <<
963 <<
964 //When the PDF vanishes at one of the inte <<
965 if ((*pdfi)[0] < 1e-35) <<
966 (*pdfi)[0] = 1e-5*pdfmax; <<
967 if ((*pdfi)[pdfi->size()-1] < 1e-35) <<
968 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; <<
969 <<
970 G4double pli = (*pdfi)[0]*factor; <<
971 G4double pui = (*pdfi)[pdfi->size()-1]*fac <<
972 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal <<
973 G4double A_temp = (1.0/(pli*dxLocal))-1.0- <<
974 G4double C_temp = 1.0+A_temp+B_temp; <<
975 if (C_temp < 1e-35) <<
976 { <<
977 (*a)[i]= 0.; <<
978 (*b)[i] = 0.; <<
979 (*c)[i] = 1; <<
980 } <<
981 else <<
982 { <<
983 (*a)[i]= A_temp; <<
984 (*b)[i] = B_temp; <<
985 (*c)[i] = C_temp; <<
986 } <<
987 //OK, now get ERR(I), the integral of the <<
988 //and the true pdf, extended over the inte <<
989 G4int icase = 1; //loop code <<
990 G4bool reLoop = false; <<
991 do <<
992 { <<
993 reLoop = false; <<
994 (*err)[i] = 0.; //zero variable <<
995 for (G4int k=0;k<nip;k++) <<
996 { <<
997 G4double rr = (*sumi)[k]; <<
998 G4double pap = (*area)[i]*(1.0+((*a)[i]+ <<
999 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1 <<
1000 if (k == 0 || k == nip-1) <<
1001 (*err)[i] += 0.5*std::fabs(pap-(*pdfi <<
1002 else <<
1003 (*err)[i] += std::fabs(pap-(*pdfi)[k] <<
1004 } <<
1005 (*err)[i] *= dxi; <<
1006 <<
1007 //If err(I) is too large, the pdf is <<
1008 if ((*err)[i] > 0.1*(*area)[i] && ica <<
1009 { <<
1010 (*b)[i] = 0; <<
1011 (*a)[i] = 0; <<
1012 (*c)[i] = 1.; <<
1013 icase = 2; <<
1014 reLoop = true; <<
1015 } <<
1016 }while(reLoop); <<
1017 delete pdfi; <<
1018 delete pdfih; <<
1019 delete sumi; <<
1020 } <<
1021 }while(x->size() < np); <<
1022 <<
1023 if (x->size() != np || a->size() != np || <<
1024 err->size() != np || area->size() != np <<
1025 { <<
1026 G4Exception("G4PenelopeRayleighModel::I <<
1027 "em2050",FatalException, <<
1028 "Problem in building the extended Table <<
1029 } <<
1030 <<
1031 /****************************************** <<
1032 Renormalization <<
1033 ******************************************* <<
1034 G4double ws = 0; <<
1035 for (std::size_t i=0;i<np-1;++i) <<
1036 ws += (*area)[i]; <<
1037 ws = 1.0/ws; <<
1038 G4double errMax = 0; <<
1039 for (std::size_t i=0;i<np-1;++i) <<
1040 { <<
1041 (*area)[i] *= ws; <<
1042 (*err)[i] *= ws; <<
1043 errMax = std::max(errMax,(*err)[i]); <<
1044 } <<
1045 <<
1046 //Vector with the normalized cumulative dis <<
1047 G4DataVector* PAC = new G4DataVector(); <<
1048 PAC->push_back(0.); <<
1049 for (std::size_t i=0;i<np-1;++i) <<
1050 { <<
1051 G4double previous = (*PAC)[i]; <<
1052 PAC->push_back(previous+(*area)[i]); <<
1053 } <<
1054 (*PAC)[PAC->size()-1] = 1.; <<
1055 <<
1056 /****************************************** <<
1057 Pre-calculated limits for the initial binar <<
1058 ******************************************* <<
1059 std::vector<std::size_t> *ITTL = new std::v <<
1060 std::vector<std::size_t> *ITTU = new std::v <<
1061 <<
1062 //Just create place-holders <<
1063 for (std::size_t i=0;i<np;++i) <<
1064 { <<
1065 ITTL->push_back(0); <<
1066 ITTU->push_back(0); <<
1067 } <<
1068 <<
1069 G4double bin = 1.0/(np-1); <<
1070 (*ITTL)[0]=0; <<
1071 for (std::size_t i=1;i<(np-1);++i) <<
1072 { <<
1073 G4double ptst = i*bin; <<
1074 G4bool found = false; <<
1075 for (std::size_t j=(*ITTL)[i-1];j<np && <<
1076 { <<
1077 if ((*PAC)[j] > ptst) <<
1078 { <<
1079 (*ITTL)[i] = j-1; <<
1080 (*ITTU)[i-1] = j; <<
1081 found = true; <<
1082 } <<
1083 } <<
1084 } <<
1085 (*ITTU)[ITTU->size()-2] = ITTU->size()-1; <<
1086 (*ITTU)[ITTU->size()-1] = ITTU->size()-1; <<
1087 (*ITTL)[ITTL->size()-1] = ITTU->size()-2; <<
1088 <<
1089 if (ITTU->size() != np || ITTU->size() != n <<
1090 { <<
1091 G4Exception("G4PenelopeRayleighModel::I <<
1092 "em2051",FatalException, <<
1093 "Problem in building the Limit Tables f <<
1094 } <<
1095 <<
1096 /****************************************** <<
1097 Copy tables <<
1098 ******************************************* <<
1099 G4PenelopeSamplingData* theTable = new G4Pe <<
1100 for (std::size_t i=0;i<np;++i) <<
1101 { <<
1102 theTable->AddPoint((*x)[i],(*PAC)[i],(* <<
1103 } <<
1104 <<
1105 if (fVerboseLevel > 2) <<
1106 { <<
1107 G4cout << "**************************** <<
1108 G4endl; <<
1109 G4cout << "Sampling table for Penelope <<
1110 theTable->DumpTable(); <<
1111 } <<
1112 fSamplingTable->insert(std::make_pair(mat,t <<
1113 <<
1114 //Clean up temporary vectors <<
1115 delete x; <<
1116 delete a; <<
1117 delete b; <<
1118 delete c; <<
1119 delete err; <<
1120 delete area; <<
1121 delete PAC; <<
1122 delete ITTL; <<
1123 delete ITTU; <<
1124 <<
1125 //DONE! <<
1126 return; <<
1127 } 480 }
1128 481
1129 //....oooOO0OOooo........oooOO0OOooo........o 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1130 483
1131 void G4PenelopeRayleighModel::GetPMaxTable(co << 484 void G4PenelopeRayleighModel::InitialiseSampling()
1132 { 485 {
1133 if (!fPMaxTable) << 486 if (!samplingFunction_x || !samplingFunction_xNoLog)
1134 { <<
1135 G4cout << "G4PenelopeRayleighModel::Bui <<
1136 G4cout << "Going to instanziate the fPM <<
1137 G4cout << "That should _not_ be here! " <<
1138 fPMaxTable = new std::map<const G4Mater <<
1139 } <<
1140 //check if the table is already there <<
1141 if (fPMaxTable->count(mat)) <<
1142 return; <<
1143 <<
1144 //otherwise build it <<
1145 if (!fSamplingTable) <<
1146 { <<
1147 G4Exception("G4PenelopeRayleighModel::G <<
1148 "em2052",FatalException, <<
1149 "SamplingTable is not properly instanti <<
1150 return; <<
1151 } <<
1152 <<
1153 //This should not be: the sampling table is <<
1154 if (!fSamplingTable->count(mat)) <<
1155 { 487 {
1156 G4ExceptionDescription ed; << 488 G4cout << "G4PenelopeRayleighModel::InitialiseSampling(), something wrong" << G4endl;
1157 ed << "Sampling table for material " < << 489 G4cout << "It looks like G4PenelopeRayleighModel::PrepareConstants() has not been called" << G4endl;
1158 G4Exception("G4PenelopeRayleighModel:: << 490 G4Exception();
1159 "em2052",FatalException, <<
1160 ed); <<
1161 return; <<
1162 } 491 }
1163 << 492 if (!SamplingTable.count(theMaterial)) //material not defined yet
1164 G4PenelopeSamplingData *theTable = fSamplin << 493 {
1165 std::size_t tablePoints = theTable->GetNumb << 494 G4double XlowInte = 0.;
1166 << 495 G4double XhighInte=(*samplingFunction_xNoLog)[0];
1167 std::size_t nOfEnergyPoints = fLogEnergyGri << 496 //material not inizialized yet: initialize it now
1168 G4PhysicsFreeVector* theVec = new G4Physics << 497 G4DataVector* samplingFunction_y = new G4DataVector();
1169 << 498 G4PenelopeIntegrator<G4PenelopeRayleighModel,G4double(G4PenelopeRayleighModel::*)(G4double)>
1170 const std::size_t nip = 51; //hard-coded in << 499 theIntegrator;
1171 << 500 G4double sum = theIntegrator.Calculate(this,&G4PenelopeRayleighModel::MolecularFormFactor,
1172 for (std::size_t ie=0;ie<fLogEnergyGridPMax << 501 XlowInte,XhighInte,1e-10);
1173 { << 502 samplingFunction_y->push_back(sum);
1174 G4double energy = G4Exp(fLogEnergyGridP << 503 for (G4int i=1;i<nPoints;i++)
1175 G4double Qm = 2.0*energy/electron_mass_ <<
1176 G4double Qm2 = Qm*Qm; <<
1177 G4double firstQ2 = theTable->GetX(0); <<
1178 G4double lastQ2 = theTable->GetX(tableP <<
1179 G4double thePMax = 0; <<
1180 <<
1181 if (Qm2 > firstQ2) <<
1182 { <<
1183 if (Qm2 < lastQ2) <<
1184 { <<
1185 //bisection to look for the index of <<
1186 std::size_t lowerBound = 0; <<
1187 std::size_t upperBound = tablePoints- <<
1188 while (lowerBound <= upperBound) <<
1189 { <<
1190 std::size_t midBin = (lowerBound + uppe <<
1191 if( Qm2 < theTable->GetX(midBin)) <<
1192 { upperBound = midBin-1; } <<
1193 else <<
1194 { lowerBound = midBin+1; } <<
1195 } <<
1196 //upperBound is the output (but also <<
1197 G4double Q1 = theTable->GetX(upperBou <<
1198 G4double Q2 = Qm2; <<
1199 G4double DQ = (Q2-Q1)/((G4double)(nip <<
1200 G4double theA = theTable->GetA(upperB <<
1201 G4double theB = theTable->GetB(upperB <<
1202 G4double thePAC = theTable->GetPAC(up <<
1203 G4DataVector* fun = new G4DataVector( <<
1204 for (std::size_t k=0;k<nip;++k) <<
1205 { <<
1206 G4double qi = Q1 + k*DQ; <<
1207 G4double tau = (qi-Q1)/ <<
1208 (theTable->GetX(upperBound+1)-Q1); <<
1209 G4double con1 = 2.0*theB*tau; <<
1210 G4double ci = 1.0+theA+theB; <<
1211 G4double con2 = ci-theA*tau; <<
1212 G4double etap = 0; <<
1213 if (std::fabs(con1) > 1.0e-16*std::fabs <<
1214 etap = con2*(1.0-std::sqrt(1.0-2.0*ta <<
1215 else <<
1216 etap = tau/con2; <<
1217 G4double theFun = (theTable->GetPAC(upp <<
1218 (1.0+(theA+theB*etap)*etap)*(1.0+(the <<
1219 ((1.0-theB*etap*etap)*ci*(theTable->G <<
1220 fun->push_back(theFun); <<
1221 } <<
1222 //Now intergrate numerically the fun <<
1223 G4DataVector* sum = new G4DataVector; <<
1224 G4double CONS = DQ*(1./12.); <<
1225 G4double HCONS = 0.5*CONS; <<
1226 sum->push_back(0.); <<
1227 G4double secondPoint = (*sum)[0] + <<
1228 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C <<
1229 sum->push_back(secondPoint); <<
1230 for (std::size_t hh=2;hh<nip-1;++hh) <<
1231 { <<
1232 G4double previous = (*sum)[hh-1]; <<
1233 G4double next = previous+(13.0*((*fun)[ <<
1234 (*fun)[hh+1]-(*fun)[hh-2])*HCON <<
1235 sum->push_back(next); <<
1236 } <<
1237 G4double last = (*sum)[nip-2]+(5.0*(* <<
1238 (*fun)[nip-3])*CONS; <<
1239 sum->push_back(last); <<
1240 thePMax = thePAC + (*sum)[sum->size() <<
1241 delete fun; <<
1242 delete sum; <<
1243 } <<
1244 else <<
1245 { <<
1246 thePMax = 1.0; <<
1247 } <<
1248 } <<
1249 else <<
1250 { 504 {
1251 thePMax = theTable->GetPAC(0); << 505 XlowInte=(*samplingFunction_xNoLog)[i-1];
>> 506 XhighInte=(*samplingFunction_xNoLog)[i];
>> 507 sum += theIntegrator.Calculate(this,
>> 508 &G4PenelopeRayleighModel::MolecularFormFactor,
>> 509 XlowInte,XhighInte,1e-10);
>> 510 samplingFunction_y->push_back(sum);
1252 } 511 }
1253 << 512 for (G4int i=0;i<nPoints;i++)
1254 //Write number in the table << 513 (*samplingFunction_y)[i]=std::log((*samplingFunction_y)[i]);
1255 theVec->PutValue(ie,energy,thePMax); << 514
>> 515 //
>> 516 /*
>> 517 G4String nnn = theMaterial->GetName()+".ntab";
>> 518 std::ofstream file(nnn);
>> 519 for (size_t k=0;k<samplingFunction_x->size();k++)
>> 520 file << (*samplingFunction_x)[k] << " " << (*samplingFunction_y)[k] << G4endl;
>> 521 file.close();
>> 522 */
>> 523 //
>> 524 SamplingTable[theMaterial] = samplingFunction_y;
1256 } 525 }
1257 <<
1258 fPMaxTable->insert(std::make_pair(mat,theVe <<
1259 return; <<
1260 } 526 }
1261 527
1262 //....oooOO0OOooo........oooOO0OOooo........o << 528 void G4PenelopeRayleighModel::PrepareConstants()
1263 <<
1264 void G4PenelopeRayleighModel::DumpFormFactorT <<
1265 { 529 {
1266 G4cout << "******************************** << 530 if (verboseLevel > 3)
1267 G4cout << "G4PenelopeRayleighModel: Form Fa << 531 G4cout << "Calling G4PenelopeRayleighModel::PrepareConstants()" << G4endl;
1268 //try to use the same format as Penelope-Fo << 532 nPoints=241;
1269 G4cout << "Q/(m_e*c) F(Q) << 533 Xlow=1e-04;
1270 G4cout << "******************************** << 534 Xhigh=1e06;
1271 if (!fLogFormFactorTable->count(mat)) << 535 if (samplingFunction_x)
1272 BuildFormFactorTable(mat); << 536 delete samplingFunction_x;
1273 << 537 if (samplingFunction_xNoLog)
1274 G4PhysicsFreeVector* theVec = fLogFormFacto << 538 delete samplingFunction_xNoLog;
1275 for (std::size_t i=0;i<theVec->GetVectorLen << 539
1276 { << 540 samplingFunction_x = new G4DataVector();
1277 G4double logQ2 = theVec->GetLowEdgeEner << 541 samplingFunction_xNoLog = new G4DataVector();
1278 G4double Q = G4Exp(0.5*logQ2); << 542
1279 G4double logF2 = (*theVec)[i]; << 543 G4double scalingFactor = std::pow((Xhigh/Xlow),((1/((G4double)(nPoints-1)))));
1280 G4double F = G4Exp(0.5*logF2); << 544 G4double logScalingFactor=(1/((G4double)(nPoints-1)))*std::log(Xhigh/Xlow);
1281 G4cout << Q << " " << F << << 545 //Logarithmic table between log(Xlow) and log(Xhigh)
>> 546 samplingFunction_x->push_back(std::log(Xlow));
>> 547 //Table between Xlow and Xhigh with log spacement: needed for InitialiseSampling()
>> 548 samplingFunction_xNoLog->push_back(Xlow);
>> 549
>> 550 for (G4int i=1;i<nPoints;i++)
>> 551 {
>> 552 G4double nextx = (*samplingFunction_x)[i-1]+logScalingFactor;
>> 553 samplingFunction_x->push_back(nextx);
>> 554 nextx = (*samplingFunction_xNoLog)[i-1]*scalingFactor;
>> 555 samplingFunction_xNoLog->push_back(nextx);
1282 } 556 }
1283 //DONE <<
1284 return; <<
1285 } <<
1286 <<
1287 //....oooOO0OOooo........oooOO0OOooo........o <<
1288 <<
1289 void G4PenelopeRayleighModel::SetParticle(con <<
1290 { <<
1291 if(!fParticle) { <<
1292 fParticle = p; <<
1293 } <<
1294 } 557 }
1295 558