| 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: G4PenelopeGammaConversionModel.cc,v 1.7 2010-11-25 09:45:13 pandola Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
26 // 28 //
27 // Author: Luciano Pandola 29 // Author: Luciano Pandola
28 // 30 //
29 // History: 31 // History:
30 // -------- 32 // --------
31 // 13 Jan 2010 L Pandola First implementa << 33 // 06 Oct 2008 L Pandola Migration from process to model
32 // 24 May 2011 L Pandola Renamed (make v2 << 34 // 17 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
33 // 18 Sep 2013 L Pandola Migration to MT << 35 // - apply internal high-energy limit only in constructor
34 // data and create << 36 // - do not apply low-energy limit (default is 0)
>> 37 // - do not apply production threshold on level of the model
>> 38 // 19 May 2009 L Pandola Explicitely set to zero pointers deleted in
>> 39 // Initialise(), since they might be checked later on
35 // 40 //
36 41
37 #include "G4PenelopeGammaConversionModel.hh" 42 #include "G4PenelopeGammaConversionModel.hh"
38 #include "G4PhysicalConstants.hh" <<
39 #include "G4SystemOfUnits.hh" <<
40 #include "G4ParticleDefinition.hh" 43 #include "G4ParticleDefinition.hh"
41 #include "G4MaterialCutsCouple.hh" 44 #include "G4MaterialCutsCouple.hh"
42 #include "G4ProductionCutsTable.hh" 45 #include "G4ProductionCutsTable.hh"
43 #include "G4DynamicParticle.hh" 46 #include "G4DynamicParticle.hh"
44 #include "G4Element.hh" 47 #include "G4Element.hh"
45 #include "G4Gamma.hh" 48 #include "G4Gamma.hh"
46 #include "G4Electron.hh" 49 #include "G4Electron.hh"
47 #include "G4Positron.hh" 50 #include "G4Positron.hh"
48 #include "G4PhysicsFreeVector.hh" << 51 #include "G4CrossSectionHandler.hh"
49 #include "G4MaterialCutsCouple.hh" << 52 #include "G4VEMDataSet.hh"
50 #include "G4AutoLock.hh" <<
51 #include "G4Exp.hh" <<
52 53
53 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 const G4int G4PenelopeGammaConversionModel::fM <<
55 G4PhysicsFreeVector* G4PenelopeGammaConversion <<
56 G4double G4PenelopeGammaConversionModel::fAtom <<
57 1.2281e+02,7.3167e+01,6.9 <<
58 6.4696e+01,6.1228e+01,5.7 <<
59 5.0787e+01,4.7851e+01,4.6 <<
60 4.4503e+01,4.3815e+01,4.3 <<
61 4.1586e+01,4.0953e+01,4.0 <<
62 3.9756e+01,3.9144e+01,3.8 <<
63 3.7174e+01,3.6663e+01,3.5 <<
64 3.4688e+01,3.4197e+01,3.3 <<
65 3.3068e+01,3.2740e+01,3.2 <<
66 3.1884e+01,3.1622e+01,3.1 <<
67 3.0950e+01,3.0758e+01,3.0 <<
68 3.0097e+01,2.9832e+01,2.9 <<
69 2.9247e+01,2.9085e+01,2.8 <<
70 2.8580e+01,2.8442e+01,2.8 <<
71 2.7973e+01,2.7819e+01,2.7 <<
72 2.7285e+01,2.7093e+01,2.6 <<
73 2.6516e+01,2.6304e+01,2.6 <<
74 2.5730e+01,2.5577e+01,2.5 <<
75 2.5100e+01,2.4941e+01,2.4 <<
76 2.4506e+01,2.4391e+01,2.4 <<
77 2.4039e+01,2.3922e+01,2.3 <<
78 2.3621e+01,2.3523e+01,2.3 <<
79 2.3238e+01,2.3139e+01,2.3 <<
80 2.2833e+01,2.2694e+01,2.2 <<
81 2.2446e+01,2.2358e+01,2.2 <<
82 55
83 //....oooOO0OOooo........oooOO0OOooo........oo <<
84 56
85 G4PenelopeGammaConversionModel::G4PenelopeGamm << 57 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition*,
86 const G4String& nam) << 58 const G4String& nam)
87 :G4VEmModel(nam),fParticleChange(nullptr),fP << 59 :G4VEmModel(nam),fTheScreeningRadii(0),crossSectionHandler(0),isInitialised(false)
88 fEffectiveCharge(nullptr),fMaterialInvScree <<
89 fScreeningFunction(nullptr),fIsInitialised( <<
90 { 60 {
91 fIntrinsicLowEnergyLimit = 2.0*electron_mass 61 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
92 fIntrinsicHighEnergyLimit = 100.0*GeV; 62 fIntrinsicHighEnergyLimit = 100.0*GeV;
93 fSmallEnergy = 1.1*MeV; 63 fSmallEnergy = 1.1*MeV;
94 64
95 if (part) <<
96 SetParticle(part); <<
97 <<
98 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 65 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
99 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 66 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
100 // 67 //
101 fVerboseLevel= 0; << 68 verboseLevel= 0;
102 // Verbosity scale: 69 // Verbosity scale:
103 // 0 = nothing << 70 // 0 = nothing
104 // 1 = warning for energy non-conservation << 71 // 1 = warning for energy non-conservation
105 // 2 = details of energy budget 72 // 2 = details of energy budget
106 // 3 = calculation of cross sections, file o 73 // 3 = calculation of cross sections, file openings, sampling of atoms
107 // 4 = entering in methods 74 // 4 = entering in methods
108 } 75 }
109 76
110 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 78
112 G4PenelopeGammaConversionModel::~G4PenelopeGam 79 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
113 { << 80 {
114 //Delete shared tables, they exist only in t << 81 if (crossSectionHandler) delete crossSectionHandler;
115 if (IsMaster() || fLocalTable) << 82 if (fTheScreeningRadii) delete fTheScreeningRadii;
116 { <<
117 for(G4int i=0; i<=fMaxZ; ++i) <<
118 { <<
119 if(fLogAtomicCrossSection[i]) { <<
120 delete fLogAtomicCrossSection[i]; <<
121 fLogAtomicCrossSection[i] = nullptr; <<
122 } <<
123 } <<
124 if (fEffectiveCharge) <<
125 delete fEffectiveCharge; <<
126 if (fMaterialInvScreeningRadius) <<
127 delete fMaterialInvScreeningRadius; <<
128 if (fScreeningFunction) <<
129 delete fScreeningFunction; <<
130 } <<
131 } 83 }
132 84
133 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134 86
135 void G4PenelopeGammaConversionModel::Initialis << 87 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition*,
136 const G4DataVector&) << 88 const G4DataVector& )
137 { 89 {
138 if (fVerboseLevel > 3) << 90 if (verboseLevel > 3)
139 G4cout << "Calling G4PenelopeGammaConvers 91 G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl;
140 92
141 SetParticle(part); << 93 //Delete the old cross section handler, if necessary
142 << 94 if (crossSectionHandler)
143 //Only the master model creates/fills/destro <<
144 if (IsMaster() && part == fParticle) <<
145 { 95 {
146 //delete old material data... << 96 crossSectionHandler->Clear();
147 if (fEffectiveCharge) << 97 delete crossSectionHandler;
148 { << 98 crossSectionHandler = 0;
149 delete fEffectiveCharge; << 99 }
150 fEffectiveCharge = nullptr; << 100
151 } << 101 //Re-initialize cross section handler
152 if (fMaterialInvScreeningRadius) << 102 crossSectionHandler = new G4CrossSectionHandler();
153 { << 103 crossSectionHandler->Initialise(0,fIntrinsicLowEnergyLimit,HighEnergyLimit(),400);
154 delete fMaterialInvScreeningRadius; << 104 crossSectionHandler->Clear();
155 fMaterialInvScreeningRadius = nullptr; << 105 G4String crossSectionFile = "penelope/pp-cs-pen-";
156 } << 106 crossSectionHandler->LoadData(crossSectionFile);
157 if (fScreeningFunction) << 107 //This is used to retrieve cross section values later on
158 { << 108 G4VEMDataSet* emdata =
159 delete fScreeningFunction; << 109 crossSectionHandler->BuildMeanFreePathForMaterials();
160 fScreeningFunction = nullptr; << 110 //The method BuildMeanFreePathForMaterials() is required here only to force
161 } << 111 //the building of an internal table: the output pointer can be deleted
162 //and create new ones << 112 delete emdata;
163 fEffectiveCharge = new std::map<const G4 << 113
164 fMaterialInvScreeningRadius = new std::m << 114 if (verboseLevel > 2)
165 fScreeningFunction = new std::map<const << 115 G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl;
166 << 116
167 G4ProductionCutsTable* theCoupleTable = << 117 if (verboseLevel > 0) {
168 G4ProductionCutsTable::GetProductionCutsTabl << 118 G4cout << "Penelope Gamma Conversion model is initialized " << G4endl
169 << 119 << "Energy range: "
170 for (G4int i=0;i<(G4int)theCoupleTable-> << 120 << LowEnergyLimit() / MeV << " MeV - "
171 { << 121 << HighEnergyLimit() / GeV << " GeV"
172 const G4Material* material = << 122 << G4endl;
173 theCoupleTable->GetMaterialCutsCouple(i) << 123 }
174 const G4ElementVector* theElementVector = <<
175 <<
176 for (std::size_t j=0;j<material->GetNumber <<
177 { <<
178 G4int iZ = theElementVector->at(j)->Ge <<
179 //read data files only in the master <<
180 if (iZ <= fMaxZ && !fLogAtomicCrossSe <<
181 ReadDataFile(iZ); <<
182 } <<
183 124
184 //check if material data are available << 125 if(isInitialised) return;
185 if (!fEffectiveCharge->count(material)) <<
186 InitializeScreeningFunctions(material); <<
187 } <<
188 if (fVerboseLevel > 0) { <<
189 G4cout << "Penelope Gamma Conversion model v <<
190 << "Energy range: " <<
191 << LowEnergyLimit() / MeV << " MeV - <<
192 << HighEnergyLimit() / GeV << " GeV" <<
193 << G4endl; <<
194 } <<
195 } <<
196 if(fIsInitialised) return; <<
197 fParticleChange = GetParticleChangeForGamma( 126 fParticleChange = GetParticleChangeForGamma();
198 fIsInitialised = true; << 127 isInitialised = true;
199 } 128 }
200 129
201 //....oooOO0OOooo........oooOO0OOooo........oo 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 131
203 void G4PenelopeGammaConversionModel::Initialis <<
204 G4VEmModel *masterModel) <<
205 { <<
206 if (fVerboseLevel > 3) <<
207 G4cout << "Calling G4PenelopeGammaConvers <<
208 // <<
209 //Check that particle matches: one might hav <<
210 //for e+ and e-). <<
211 // <<
212 if (part == fParticle) <<
213 { <<
214 //Get the const table pointers from the <<
215 const G4PenelopeGammaConversionModel* th <<
216 static_cast<G4PenelopeGammaConversionModel*> <<
217 <<
218 //Copy pointers to the data tables <<
219 fEffectiveCharge = theModel->fEffectiveC <<
220 fMaterialInvScreeningRadius = theModel-> <<
221 fScreeningFunction = theModel->fScreenin <<
222 for(G4int i=0; i<=fMaxZ; ++i) <<
223 fLogAtomicCrossSection[i] = theModel->fLogAt <<
224 <<
225 //Same verbosity for all workers, as the <<
226 fVerboseLevel = theModel->fVerboseLevel; <<
227 } <<
228 <<
229 return; <<
230 } <<
231 <<
232 //....oooOO0OOooo........oooOO0OOooo........oo <<
233 namespace { G4Mutex PenelopeGammaConversionMo <<
234 <<
235 G4double G4PenelopeGammaConversionModel::Compu 132 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
236 const G4ParticleDefinition << 133 const G4ParticleDefinition*,
237 G4double energy, << 134 G4double energy,
238 G4double Z, G4double, << 135 G4double Z, G4double,
239 G4double, G4double) << 136 G4double, G4double)
240 { 137 {
241 // 138 //
242 // Penelope model v2008. << 139 // Penelope model.
243 // Cross section (including triplet producti << 140 // Cross section (including triplet production) read from database and managed
244 // through the G4CrossSectionHandler utility 141 // through the G4CrossSectionHandler utility. Cross section data are from
245 // M.J. Berger and J.H. Hubbel (XCOM), Repor 142 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
246 // 143 //
>> 144
>> 145 if (verboseLevel > 3)
>> 146 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
>> 147
>> 148 G4int iZ = (G4int) Z;
>> 149 // if (!crossSectionHandler) //VI: should not be checked in run time
>> 150 // {
>> 151 // G4cout << "G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom" << G4endl;
>> 152 // G4cout << "The cross section handler is not correctly initialized" << G4endl;
>> 153 // G4Exception();
>> 154 // }
>> 155 G4double cs = crossSectionHandler->FindValue(iZ,energy);
247 156
248 if (energy < fIntrinsicLowEnergyLimit) << 157 if (verboseLevel > 2)
249 return 0; << 158 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
250 <<
251 G4int iZ = G4int(Z); <<
252 <<
253 if (!fLogAtomicCrossSection[iZ]) <<
254 { <<
255 //If we are here, it means that Initial <<
256 //not filled up. This can happen in a U <<
257 if (fVerboseLevel > 0) <<
258 { <<
259 //Issue a G4Exception (warning) only in v <<
260 G4ExceptionDescription ed; <<
261 ed << "Unable to retrieve the cross secti <<
262 ed << "This can happen only in Unit Tests <<
263 G4Exception("G4PenelopeGammaConversionMod <<
264 "em2018",JustWarning,ed); <<
265 } <<
266 //protect file reading via autolock <<
267 G4AutoLock lock(&PenelopeGammaConversio <<
268 ReadDataFile(iZ); <<
269 lock.unlock(); <<
270 fLocalTable = true; <<
271 } <<
272 G4double cs = 0; <<
273 G4double logene = G4Log(energy); <<
274 G4PhysicsFreeVector* theVec = fLogAtomicCros <<
275 G4double logXS = theVec->Value(logene); <<
276 cs = G4Exp(logXS); <<
277 <<
278 if (fVerboseLevel > 2) <<
279 G4cout << "Gamma conversion cross section <<
280 " = " << cs/barn << " barn" << G4endl; 159 " = " << cs/barn << " barn" << G4endl;
281 return cs; 160 return cs;
282 } 161 }
283 162
284 //....oooOO0OOooo........oooOO0OOooo........oo 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285 164
286 void << 165 void
287 G4PenelopeGammaConversionModel::SampleSecondar 166 G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
288 const G4MaterialCutsCouple* coup 167 const G4MaterialCutsCouple* couple,
289 const G4DynamicParticle* aDynami 168 const G4DynamicParticle* aDynamicGamma,
290 G4double, 169 G4double,
291 G4double) 170 G4double)
292 { 171 {
293 // 172 //
294 // Penelope model v2008. << 173 // Penelope model.
295 // Final state is sampled according to the B 174 // Final state is sampled according to the Bethe-Heitler model with Coulomb
296 // corrections, according to the semi-empiri 175 // corrections, according to the semi-empirical model of
297 // J. Baro' et al., Radiat. Phys. Chem. 44 176 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
298 // 177 //
299 // The model uses the high energy Coulomb co << 178 // The model uses the high energy Coulomb correction from
300 // H. Davies et al., Phys. Rev. 93 (1954) 7 179 // H. Davies et al., Phys. Rev. 93 (1954) 788
301 // and atomic screening radii tabulated from << 180 // and atomic screening radii tabulated from
302 // J.H. Hubbel et al., J. Phys. Chem. Ref. 181 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
303 // for Z= 1 to 92. << 182 // for Z= 1 to 92. This managed in this model by the method
>> 183 // GetScreeningRadius().
304 // 184 //
305 if (fVerboseLevel > 3) << 185 if (verboseLevel > 3)
306 G4cout << "Calling SamplingSecondaries() o 186 G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
307 187
308 G4double photonEnergy = aDynamicGamma->GetKi 188 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
309 189
310 // Always kill primary 190 // Always kill primary
311 fParticleChange->ProposeTrackStatus(fStopAnd 191 fParticleChange->ProposeTrackStatus(fStopAndKill);
312 fParticleChange->SetProposedKineticEnergy(0. 192 fParticleChange->SetProposedKineticEnergy(0.);
313 193
314 if (photonEnergy <= fIntrinsicLowEnergyLimit 194 if (photonEnergy <= fIntrinsicLowEnergyLimit)
315 { 195 {
316 fParticleChange->ProposeLocalEnergyDepos 196 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
317 return ; 197 return ;
318 } 198 }
319 199
320 G4ParticleMomentum photonDirection = aDynami 200 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
321 const G4Material* mat = couple->GetMaterial( <<
322 201
323 //Either Initialize() was not called, or we << 202 G4double eps ;
324 //not invoked << 203 G4double eki = electron_mass_c2 / photonEnergy ;
325 if (!fEffectiveCharge) <<
326 { <<
327 //create a **thread-local** version of t <<
328 //Unit Tests <<
329 fLocalTable = true; <<
330 fEffectiveCharge = new std::map<const G4 <<
331 fMaterialInvScreeningRadius = new std::m <<
332 fScreeningFunction = new std::map<const <<
333 } <<
334 204
335 if (!fEffectiveCharge->count(mat)) << 205 // Do it fast if photon energy < 1.1 MeV
>> 206 if (photonEnergy < fSmallEnergy )
336 { 207 {
337 //If we are here, it means that Initiali << 208 eps = eki + (1-2*eki) * G4UniformRand();
338 //not filled up. This can happen in a Un <<
339 if (fVerboseLevel > 0) <<
340 { <<
341 //Issue a G4Exception (warning) only in ve <<
342 G4ExceptionDescription ed; <<
343 ed << "Unable to allocate the EffectiveCha <<
344 mat->GetName() << G4endl; <<
345 ed << "This can happen only in Unit Tests" <<
346 G4Exception("G4PenelopeGammaConversionMode <<
347 "em2019",JustWarning,ed); <<
348 } <<
349 //protect file reading via autolock <<
350 G4AutoLock lock(&PenelopeGammaConversion <<
351 InitializeScreeningFunctions(mat); <<
352 lock.unlock(); <<
353 } 209 }
354 <<
355 // eps is the fraction of the photon energy <<
356 G4double eps = 0; <<
357 G4double eki = electron_mass_c2/photonEnergy <<
358 <<
359 //Do it fast for photon energy < 1.1 MeV (cl <<
360 if (photonEnergy < fSmallEnergy) <<
361 eps = eki + (1.0-2.0*eki)*G4UniformRand(); <<
362 else 210 else
363 { 211 {
364 //Complete calculation << 212 // Select randomly one element in the current material
365 G4double effC = fEffectiveCharge->find(m << 213 if (verboseLevel > 2)
366 G4double alz = effC*fine_structure_const << 214 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
367 G4double T = std::sqrt(2.0*eki); << 215 //use crossSectionHandler instead of G4EmElementSelector because in this case
368 G4double F00=(-1.774-1.210e1*alz+1.118e1 << 216 //the dimension of the table is equal to the dimension of the database
369 +(8.523+7.326e1*alz-4.441e1*alz*alz)* << 217 //(less interpolation errors)
370 -(1.352e1+1.211e2*alz-9.641e1*alz*alz << 218 G4int Z_int = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
371 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T << 219 if (verboseLevel > 2)
372 << 220 G4cout << "Selected Z = " << Z_int << G4endl;
373 G4double F0b = fScreeningFunction->find( << 221
374 G4double g0 = F0b + F00; << 222 //Low energy and Coulomb corrections
375 G4double invRad = fMaterialInvScreeningR << 223 G4double Z=(G4double) Z_int;
376 G4double bmin = 4.0*eki/invRad; << 224 G4double ZAlpha = Z*fine_structure_const;
377 std::pair<G4double,G4double> scree = Ge << 225 G4double ScreenRadius = GetScreeningRadius(Z);
378 G4double g1 = scree.first; << 226 G4double funct1=0,g0=0;
379 G4double g2 = scree.second; << 227 G4double g1min=0,g2min=0;
380 G4double g1min = g1+g0; << 228 funct1 = 4.0*std::log(ScreenRadius);
381 G4double g2min = g2+g0; << 229 g0 = funct1-4*CoulombCorrection(ZAlpha)+LowEnergyCorrection(ZAlpha,eki);
382 G4double xr = 0.5-eki; << 230 G4double bmin = 2*eki*ScreenRadius;
383 G4double a1 = 2.*g1min*xr*xr/3.; << 231 std::vector<G4double> ScreenFunctionValues = ScreenFunction(bmin);
384 G4double p1 = a1/(a1+g2min); << 232 if (ScreenFunctionValues.size() != 2)
>> 233 {
>> 234 G4cout << "G4PenelopeGammaConversionModel::SampleSecondaries" << G4endl;
>> 235 G4cout << "ScreenFunction did not return 2 values! Something wrong! " << G4endl;
>> 236 G4Exception();
>> 237 }
>> 238 g1min=g0+ScreenFunctionValues[0];
>> 239 g2min=g0+ScreenFunctionValues[1];
>> 240 G4double xr,a1,p1;
>> 241 xr=0.5-eki;
>> 242 a1=(2.0/3.0)*g1min*xr*xr;
>> 243 p1=a1/(a1+g2min);
385 244
386 G4bool loopAgain = false; <<
387 //Random sampling of eps 245 //Random sampling of eps
>> 246 G4double rand1,rand2,rand3,b;
>> 247 G4double g1;
>> 248
388 do{ 249 do{
389 loopAgain = false; << 250 rand1 = G4UniformRand();
390 if (G4UniformRand() <= p1) << 251 if (rand1 < p1) {
391 { << 252 rand2 = 2.0*G4UniformRand()-1.0;
392 G4double ru2m1 = 2.0*G4UniformRand()-1. << 253 if (rand2 < 0) {
393 if (ru2m1 < 0) << 254 eps = 0.5 - xr*std::pow(std::abs(rand2),(1./3.));
394 eps = 0.5-xr*std::pow(std::abs(ru2m1), <<
395 else <<
396 eps = 0.5+xr*std::pow(ru2m1,1./3.); <<
397 G4double B = eki/(invRad*eps*(1.0-eps)); <<
398 scree = GetScreeningFunctions(B); <<
399 g1 = scree.first; <<
400 g1 = std::max(g1+g0,0.); <<
401 if (G4UniformRand()*g1min > g1) <<
402 loopAgain = true; <<
403 } 255 }
404 else << 256 else
>> 257 {
>> 258 eps = 0.5 + xr*std::pow(rand2,(1./3.));
>> 259 }
>> 260 b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
>> 261 std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b);
>> 262 g1 = g0+ScreenFunctionSampling[0];
>> 263 if (g1 < 0) g1=0;
>> 264 rand3 = G4UniformRand()*g1min;
>> 265 }
>> 266 else
405 { 267 {
406 eps = eki+2.0*xr*G4UniformRand(); 268 eps = eki+2.0*xr*G4UniformRand();
407 G4double B = eki/(invRad*eps*(1.0-eps)); << 269 b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
408 scree = GetScreeningFunctions(B); << 270 std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b);
409 g2 = scree.second; << 271 g1 = g0+ScreenFunctionSampling[1];
410 g2 = std::max(g2+g0,0.); << 272 if (g1 < 0) g1=0;
411 if (G4UniformRand()*g2min > g2) << 273 rand3 = G4UniformRand()*g2min;
412 loopAgain = true; << 274 }
413 } << 275 } while (rand3>g1);
414 }while(loopAgain); << 276 } //End of eps sampling
415 } <<
416 if (fVerboseLevel > 4) <<
417 G4cout << "Sampled eps = " << eps << G4end <<
418 277
419 G4double electronTotEnergy = eps*photonEnerg 278 G4double electronTotEnergy = eps*photonEnergy;
420 G4double positronTotEnergy = (1.0-eps)*photo 279 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
421 << 280
422 // Scattered electron (positron) angles. ( Z 281 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
423 282
424 //electron kinematics 283 //electron kinematics
425 G4double electronKineEnergy = std::max(0.,el << 284 G4double costheta_el,costheta_po;
426 G4double costheta_el = G4UniformRand()*2.0-1 << 285 G4double phi_el,phi_po;
>> 286 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
>> 287 costheta_el = G4UniformRand()*2.0-1.0;
427 G4double kk = std::sqrt(electronKineEnergy*( 288 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
428 costheta_el = (costheta_el*electronTotEnergy 289 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
429 G4double phi_el = twopi * G4UniformRand() ; << 290 phi_el = twopi * G4UniformRand() ;
430 G4double dirX_el = std::sqrt(1.-costheta_el* 291 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
431 G4double dirY_el = std::sqrt(1.-costheta_el* 292 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
432 G4double dirZ_el = costheta_el; 293 G4double dirZ_el = costheta_el;
433 294
434 //positron kinematics 295 //positron kinematics
435 G4double positronKineEnergy = std::max(0.,po 296 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
436 G4double costheta_po = G4UniformRand()*2.0-1 << 297 costheta_po = G4UniformRand()*2.0-1.0;
437 kk = std::sqrt(positronKineEnergy*(positronK 298 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
438 costheta_po = (costheta_po*positronTotEnergy 299 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
439 G4double phi_po = twopi * G4UniformRand() ; << 300 phi_po = twopi * G4UniformRand() ;
440 G4double dirX_po = std::sqrt(1.-costheta_po* 301 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
441 G4double dirY_po = std::sqrt(1.-costheta_po* 302 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
442 G4double dirZ_po = costheta_po; 303 G4double dirZ_po = costheta_po;
443 304
444 // Kinematics of the created pair: 305 // Kinematics of the created pair:
445 // the electron and positron are assumed to 306 // the electron and positron are assumed to have a symetric angular
446 // distribution with respect to the Z axis a 307 // distribution with respect to the Z axis along the parent photon
447 G4double localEnergyDeposit = 0. ; 308 G4double localEnergyDeposit = 0. ;
448 309
>> 310 //Generate explicitely the electron in the pair, only if it is > threshold
>> 311 //VI: applying cut here provides inconsistency
>> 312
449 if (electronKineEnergy > 0.0) 313 if (electronKineEnergy > 0.0)
450 { 314 {
451 G4ThreeVector electronDirection ( dirX_e 315 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
452 electronDirection.rotateUz(photonDirecti 316 electronDirection.rotateUz(photonDirection);
453 G4DynamicParticle* electron = new G4Dyna 317 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
454 electronDirection, 318 electronDirection,
455 electronKineEnergy); 319 electronKineEnergy);
456 fvect->push_back(electron); 320 fvect->push_back(electron);
457 } 321 }
458 else 322 else
459 { 323 {
460 localEnergyDeposit += electronKineEnergy 324 localEnergyDeposit += electronKineEnergy;
461 electronKineEnergy = 0; 325 electronKineEnergy = 0;
462 } 326 }
463 327
464 //Generate the positron. Real particle in an 328 //Generate the positron. Real particle in any case, because it will annihilate. If below
465 //threshold, produce it at rest 329 //threshold, produce it at rest
>> 330 // VI: here there was a bug - positron and electron cuts are different
466 if (positronKineEnergy < 0.0) 331 if (positronKineEnergy < 0.0)
467 { 332 {
468 localEnergyDeposit += positronKineEnergy 333 localEnergyDeposit += positronKineEnergy;
469 positronKineEnergy = 0; //produce it at 334 positronKineEnergy = 0; //produce it at rest
470 } 335 }
471 G4ThreeVector positronDirection(dirX_po,dirY 336 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
472 positronDirection.rotateUz(photonDirection); 337 positronDirection.rotateUz(photonDirection);
473 G4DynamicParticle* positron = new G4DynamicP 338 G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
474 positronDirection, positronK 339 positronDirection, positronKineEnergy);
475 fvect->push_back(positron); 340 fvect->push_back(positron);
476 341
477 //Add rest of energy to the local energy dep 342 //Add rest of energy to the local energy deposit
478 fParticleChange->ProposeLocalEnergyDeposit(l 343 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
479 << 344
480 if (fVerboseLevel > 1) << 345 if (verboseLevel > 1)
481 { 346 {
482 G4cout << "----------------------------- 347 G4cout << "-----------------------------------------------------------" << G4endl;
483 G4cout << "Energy balance from G4Penelop 348 G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
484 G4cout << "Incoming photon energy: " << 349 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
485 G4cout << "----------------------------- 350 G4cout << "-----------------------------------------------------------" << G4endl;
486 if (electronKineEnergy) 351 if (electronKineEnergy)
487 G4cout << "Electron (explicitly produced) " << 352 G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV"
488 << G4endl; 353 << G4endl;
489 if (positronKineEnergy) 354 if (positronKineEnergy)
490 G4cout << "Positron (not at rest) " << posit 355 G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
491 G4cout << "Rest masses of e+/- " << 2.0* 356 G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
492 if (localEnergyDeposit) 357 if (localEnergyDeposit)
493 G4cout << "Local energy deposit " << localEn 358 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
494 G4cout << "Total final state: " << (elec 359 G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
495 localEnergyDeposit+2.0*electron_ma 360 localEnergyDeposit+2.0*electron_mass_c2)/keV <<
496 " keV" << G4endl; 361 " keV" << G4endl;
497 G4cout << "----------------------------- 362 G4cout << "-----------------------------------------------------------" << G4endl;
498 } 363 }
499 if (fVerboseLevel > 0) << 364 if (verboseLevel > 0)
500 { 365 {
501 G4double energyDiff = std::fabs(electron 366 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
502 localEnergyDeposit+2.0*electron_ 367 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
503 if (energyDiff > 0.05*keV) 368 if (energyDiff > 0.05*keV)
504 G4cout << "Warning from G4PenelopeGammaConve << 369 G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: "
505 << (electronKineEnergy+positronKineEn 370 << (electronKineEnergy+positronKineEnergy+
506 localEnergyDeposit+2.0*electron_mass_c2 << 371 localEnergyDeposit+2.0*electron_mass_c2)/keV
507 << " keV (final) vs. " << photonEnerg 372 << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
508 } << 373 }
509 } 374 }
510 375
511 //....oooOO0OOooo........oooOO0OOooo........oo 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
512 377
513 void G4PenelopeGammaConversionModel::ReadDataF << 378 std::vector<G4double> G4PenelopeGammaConversionModel::ScreenFunction(G4double b)
514 { 379 {
515 if (!IsMaster()) << 380 std::vector<G4double> result;
516 //Should not be here! << 381 result.clear();
517 G4Exception("G4PenelopeGammaConversionMode << 382 G4double bsquare=b*b;
518 "em0100",FatalException,"Worker thread in << 383 G4double a0,f1,f2;
519 << 384 f1=2.0-2*std::log(1+bsquare);
520 if (fVerboseLevel > 2) << 385 f2=f1-(2.0/3.0);
521 { << 386 if (b < 1.0e-10)
522 G4cout << "G4PenelopeGammaConversionMode <<
523 G4cout << "Going to read Gamma Conversio <<
524 } <<
525 <<
526 const char* path = G4FindDataDir("G4LEDATA <<
527 if(!path) <<
528 { 387 {
529 G4String excep = << 388 f1=f1-twopi*b;
530 "G4PenelopeGammaConversionModel - G4LEDATA e <<
531 G4Exception("G4PenelopeGammaConversionMo <<
532 "em0006",FatalException,excep); <<
533 return; <<
534 } 389 }
535 <<
536 /* <<
537 Read the cross section file <<
538 */ <<
539 std::ostringstream ost; <<
540 if (Z>9) <<
541 ost << path << "/penelope/pairproduction/p <<
542 else 390 else
543 ost << path << "/penelope/pairproduction/p << 391 {
544 std::ifstream file(ost.str().c_str()); << 392 a0 = 4*b*std::atan(1.0/b);
545 if (!file.is_open()) << 393 f1 = f1 - a0;
546 { << 394 f2 = f2+2*bsquare*(4.0-a0-3*std::log((1+bsquare)/bsquare));
547 G4String excep = "G4PenelopeGammaConvers <<
548 G4String(ost.str()) + " not found!"; <<
549 G4Exception("G4PenelopeGammaConversionMo <<
550 "em0003",FatalException,excep); <<
551 } <<
552 <<
553 //I have to know in advance how many points <<
554 //to initialize the G4PhysicsFreeVector() <<
555 std::size_t ndata=0; <<
556 G4String line; <<
557 while( getline(file, line) ) <<
558 ndata++; <<
559 ndata -= 1; //remove one header line <<
560 <<
561 file.clear(); <<
562 file.close(); <<
563 file.open(ost.str().c_str()); <<
564 G4int readZ =0; <<
565 file >> readZ; <<
566 <<
567 if (fVerboseLevel > 3) <<
568 G4cout << "Element Z=" << Z << G4endl; <<
569 <<
570 //check the right file is opened. <<
571 if (readZ != Z) <<
572 { <<
573 G4ExceptionDescription ed; <<
574 ed << "Corrupted data file for Z=" << Z <<
575 G4Exception("G4PenelopeGammaConversionMo <<
576 "em0005",FatalException,ed); <<
577 } <<
578 <<
579 fLogAtomicCrossSection[Z] = new G4PhysicsFre <<
580 G4double ene=0,xs=0; <<
581 for (std::size_t i=0;i<ndata;++i) <<
582 { <<
583 file >> ene >> xs; <<
584 //dimensional quantities <<
585 ene *= eV; <<
586 xs *= barn; <<
587 if (xs < 1e-40*cm2) //protection against <<
588 xs = 1e-40*cm2; <<
589 fLogAtomicCrossSection[Z]->PutValue(i,G4 <<
590 } 395 }
591 file.close(); << 396 result.push_back(0.5*(3*f1-f2));
592 << 397 result.push_back(0.25*(3*f1+f2));
593 return; << 398 return result;
594 } 399 }
595 400
596 //....oooOO0OOooo........oooOO0OOooo........oo 401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
597 402
598 void G4PenelopeGammaConversionModel::Initializ << 403 G4double G4PenelopeGammaConversionModel::GetScreeningRadius(G4double Z)
599 { 404 {
600 // This is subroutine GPPa0 of Penelope << 405 G4double result = 0;
601 // << 406 G4bool foundElement = false;
602 // 1) calculate the effective Z for the purp << 407 G4int iZ = (G4int) Z;
603 // << 408 if (!fTheScreeningRadii)
604 G4double zeff = 0; << 409 fTheScreeningRadii = new std::map<G4int,G4double>;
605 G4int intZ = 0; << 410
606 G4int nElements = (G4int)material->GetNumber << 411 if (fTheScreeningRadii->count(iZ))
607 const G4ElementVector* elementVector = mater << 412 {
608 << 413 //The element is already loaded: just return it
609 //avoid calculations if only one building el << 414 result = fTheScreeningRadii->find(iZ)->second;
610 if (nElements == 1) << 415 return result;
611 { <<
612 zeff = (*elementVector)[0]->GetZ(); <<
613 intZ = (G4int) zeff; <<
614 } 416 }
615 else // many elements...let's do the calcula << 417 else //retrieve all from file
616 { 418 {
617 const G4double* fractionVector = materia << 419 char* path = getenv("G4LEDATA");
618 << 420 if (!path)
619 G4double atot = 0; <<
620 for (G4int i=0;i<nElements;i++) <<
621 { 421 {
622 G4double Zelement = (*elementVector)[i]->G << 422 G4String excep = "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
623 G4double Aelement = (*elementVector)[i]->G << 423 G4Exception(excep);
624 atot += Aelement*fractionVector[i]; << 424 return result;
625 zeff += Zelement*Aelement*fractionVector[i << 425 }
>> 426 G4String pathString(path);
>> 427 G4String pathFile = pathString + "/penelope/pp-pen.dat";
>> 428 std::ifstream file(pathFile);
>> 429
>> 430 if (!(file.is_open()))
>> 431 {
>> 432 G4String excep = "G4PenelopeGammaConversionModel - data file " + pathFile + "not found!";
>> 433 G4Exception(excep);
>> 434 }
>> 435 G4int k;
>> 436 G4double a1,a2;
>> 437 while(!file.eof()) {
>> 438 file >> k >> a1 >> a2;
>> 439 fTheScreeningRadii->insert(std::make_pair(k,a1));
>> 440 if ((G4double) k == Z)
>> 441 {
>> 442 result = a1;
>> 443 foundElement = true;
>> 444 }
>> 445 }
>> 446 file.close();
>> 447 if (verboseLevel > 2)
>> 448 G4cout << "Read file pp-pen.dat" << G4endl;
>> 449 if (foundElement)
>> 450 return result;
>> 451 else
>> 452 {
>> 453 G4String excep = "G4PenelopeGammaConversionModel - Screening Radius for not found in the data file";
>> 454 G4Exception(excep);
>> 455 return 0;
626 } 456 }
627 atot /= material->GetTotNbOfAtomsPerVolu <<
628 zeff /= (material->GetTotNbOfAtomsPerVol <<
629 <<
630 intZ = (G4int) (zeff+0.25); <<
631 if (intZ <= 0) <<
632 intZ = 1; <<
633 if (intZ > fMaxZ) <<
634 intZ = fMaxZ; <<
635 } <<
636 <<
637 if (fEffectiveCharge) <<
638 fEffectiveCharge->insert(std::make_pair(ma <<
639 <<
640 // <<
641 // 2) Calculate Coulomb Correction <<
642 // <<
643 G4double alz = fine_structure_const*zeff; <<
644 G4double alzSquared = alz*alz; <<
645 G4double fc = alzSquared*(0.202059-alzSquar <<
646 (0.03693-alzSquared* <<
647 (0.00835-alzSquared*(0.00201-alzSq <<
648 (0.00049-alzSquared* <<
649 (0.00012-alzSquared*0.00003))) <<
650 +1.0/(alzSquared+1.0)); <<
651 // <<
652 // 3) Screening functions and low-energy cor <<
653 // <<
654 G4double matRadius = 2.0/ fAtomicScreeningRa <<
655 if (fMaterialInvScreeningRadius) <<
656 fMaterialInvScreeningRadius->insert(std::m <<
657 <<
658 std::pair<G4double,G4double> myPair(0,0); <<
659 G4double f0a = 4.0*G4Log(fAtomicScreeningRad <<
660 G4double f0b = f0a - 4.0*fc; <<
661 myPair.first = f0a; <<
662 myPair.second = f0b; <<
663 <<
664 if (fScreeningFunction) <<
665 fScreeningFunction->insert(std::make_pair( <<
666 <<
667 if (fVerboseLevel > 2) <<
668 { <<
669 G4cout << "Average Z for material " << m <<
670 zeff << G4endl; <<
671 G4cout << "Effective radius for material <<
672 fAtomicScreeningRadius[intZ] << " m_e*c/hbar <<
673 matRadius << G4endl; <<
674 G4cout << "Screening parameters F0 for m <<
675 f0a << "," << f0b << G4endl; <<
676 } 457 }
677 return; <<
678 } 458 }
679 459
680 //....oooOO0OOooo........oooOO0OOooo........oo 460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
681 461
682 std::pair<G4double,G4double> << 462 G4double G4PenelopeGammaConversionModel::CoulombCorrection(G4double a)
683 G4PenelopeGammaConversionModel::GetScreeningFu <<
684 { 463 {
685 // This is subroutine SCHIFF of Penelope << 464 G4double fc=0;
686 // << 465 G4double b[7] = {0.202059,-0.03693,0.00835,-0.00201,0.00049,-0.00012,0.00003};
687 // Screening Functions F1(B) and F2(B) in th << 466 G4double aSquared = a*a;
688 // section for pair production << 467 G4double aFourth = aSquared*aSquared;
689 // << 468 G4double aEighth = aFourth*aFourth;
690 std::pair<G4double,G4double> result(0.,0.); << 469
691 G4double BSquared = B*B; << 470 fc = ((1.0/(1.0+a*a))+b[0]+b[1]*aSquared+b[2]*aFourth+b[3]*(aSquared*aFourth)+
692 G4double f1 = 2.0-2.0*G4Log(1.0+BSquared); << 471 b[4]*aEighth+b[5]*(aEighth*aSquared)+b[6]*(aEighth*aFourth));
693 G4double f2 = f1 - 6.66666666e-1; // (-2/3) << 472 fc=aSquared*fc;
694 if (B < 1.0e-10) << 473 return fc;
695 f1 = f1-twopi*B; <<
696 else <<
697 { <<
698 G4double a0 = 4.0*B*std::atan(1./B); <<
699 f1 = f1 - a0; <<
700 f2 += 2.0*BSquared*(4.0-a0-3.0*G4Log((1. <<
701 } <<
702 G4double g1 = 0.5*(3.0*f1-f2); <<
703 G4double g2 = 0.25*(3.0*f1+f2); <<
704 <<
705 result.first = g1; <<
706 result.second = g2; <<
707 <<
708 return result; <<
709 } 474 }
710 475
711 //....oooOO0OOooo........oooOO0OOooo........oo << 476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
712 477
713 void G4PenelopeGammaConversionModel::SetPartic << 478 G4double G4PenelopeGammaConversionModel::LowEnergyCorrection(G4double a,G4double eki)
714 { 479 {
715 if(!fParticle) { << 480 G4double f0=0,t=0;
716 fParticle = p; << 481 G4double b[12] = {-1.744,-12.10,11.18,8.523,73.26,-41.41,-13.52,-121.1,94.41,8.946,62.05,-63.41};
717 } << 482 t=std::sqrt(2.0*eki);
>> 483 G4double tSq = t*t;
>> 484 f0=(b[0]+b[1]*a+b[2]*a*a)*t+(b[3]+b[4]*a+b[5]*a*a)*(tSq)+(b[6]+b[7]*a+b[8]*a*a)*(tSq*t)+
>> 485 (b[9]+b[10]*a+b[11]*a*a)*(tSq*tSq);
>> 486 return f0;
>> 487
718 } 488 }
719 489