Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4eBremParametrizedModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class file
30//
31//
32// File name: G4eBremParametrizedModel
33//
34// Author: Andreas Schaelicke
35//
36// Creation date: 06.04.2011
37//
38// Modifications:
39//
40// Main References:
41// - based on G4eBremsstrahlungModel and G4eBremsstrahlungRelModel
42// -------------------------------------------------------------------
43//
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
49#include "G4SystemOfUnits.hh"
50#include "G4Electron.hh"
51#include "G4Positron.hh"
52#include "G4Gamma.hh"
53#include "Randomize.hh"
54#include "G4Material.hh"
55#include "G4Element.hh"
56#include "G4ElementVector.hh"
59#include "G4LossTableManager.hh"
60#include "G4ModifiedTsai.hh"
61#include "G4Exp.hh"
62#include "G4Log.hh"
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66const G4double G4eBremParametrizedModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
67 0.5917, 0.7628, 0.8983, 0.9801 };
68const G4double G4eBremParametrizedModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
69 0.1813, 0.1569, 0.1112, 0.0506 };
70
71static const G4double tlow = 1.*CLHEP::MeV;
72
73//
74// GEANT4 internal units.
75//
76static const G4double
77 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
78 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
79 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
80
81static const G4double
82 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
83 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
84 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
85
86static const G4double
87 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
88 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
89 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
90
91static const G4double
92 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
93 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
94 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
95
96using namespace std;
97
99 const G4String& nam)
100 : G4VEmModel(nam),
101 particle(nullptr),
102 fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
103 bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
104 isInitialised(false),
105 isElectron(true)
106{
108
109 minThreshold = 0.1*keV;
110 lowKinEnergy = 10.*MeV;
111 SetLowEnergyLimit(lowKinEnergy);
112
114
116
119
120 InitialiseConstants();
121 if(nullptr != p) { SetParticle(p); }
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126void G4eBremParametrizedModel::InitialiseConstants()
127{
128 facFel = G4Log(184.15);
129 facFinel = G4Log(1194.);
130}
131
132//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137
138void G4eBremParametrizedModel::SetParticle(const G4ParticleDefinition* p)
139{
140 particle = p;
142 if(p == G4Electron::Electron()) { isElectron = true; }
143 else { isElectron = false;}
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
150{
151 return minThreshold;
152}
153
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155
157 const G4Material* mat,
158 G4double kineticEnergy)
159{
160 densityFactor = mat->GetElectronDensity()*fMigdalConstant;
161
162 // calculate threshold for density effect
163 kinEnergy = kineticEnergy;
164 totalEnergy = kineticEnergy + particleMass;
166}
167
168
169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170
172 const G4DataVector& cuts)
173{
174 if(p) { SetParticle(p); }
175
176 lowKinEnergy = LowEnergyLimit();
177
178 currentZ = 0.;
179
180 if(IsMaster()) { InitialiseElementSelectors(p, cuts); }
181
182 if(isInitialised) { return; }
184 isInitialised = true;
185}
186
187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188
190 G4VEmModel* masterModel)
191{
193}
194
195//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196
198 const G4Material* material,
199 const G4ParticleDefinition* p,
200 G4double kineticEnergy,
201 G4double cutEnergy)
202{
203 if(!particle) { SetParticle(p); }
204 if(kineticEnergy < lowKinEnergy) { return 0.0; }
205 G4double cut = std::min(cutEnergy, kineticEnergy);
206 if(cut == 0.0) { return 0.0; }
207
208 SetupForMaterial(particle, material,kineticEnergy);
209
210 const G4ElementVector* theElementVector = material->GetElementVector();
211 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
212
213 G4double dedx = 0.0;
214
215 // loop for elements in the material
216 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
217
218 G4VEmModel::SetCurrentElement((*theElementVector)[i]);
219 SetCurrentElement((*theElementVector)[i]->GetZ());
220
221 dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
222 }
223 dedx *= bremFactor;
224
225 return dedx;
226}
227
228//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229
230G4double G4eBremParametrizedModel::ComputeBremLoss(G4double cut)
231{
232 G4double loss = 0.0;
233
234 // number of intervals and integration step
235 G4double vcut = cut/totalEnergy;
236 G4int n = (G4int)(20*vcut) + 3;
237 G4double delta = vcut/G4double(n);
238
239 G4double e0 = 0.0;
240 G4double xs;
241
242 // integration
243 for(G4int l=0; l<n; l++) {
244
245 for(G4int i=0; i<8; i++) {
246
247 G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
248
249 xs = ComputeDXSectionPerAtom(eg);
250
251 loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
252 }
253 e0 += delta;
254 }
255
256 loss *= delta*totalEnergy;
257
258 return loss;
259}
260
261//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262
264 const G4ParticleDefinition* p,
265 G4double kineticEnergy,
267 G4double cutEnergy,
268 G4double maxEnergy)
269{
270 if(!particle) { SetParticle(p); }
271 if(kineticEnergy < lowKinEnergy) { return 0.0; }
272 G4double cut = std::min(cutEnergy, kineticEnergy);
273 G4double tmax = std::min(maxEnergy, kineticEnergy);
274
275 if(cut >= tmax) { return 0.0; }
276
277 SetCurrentElement(Z);
278
279 G4double cross = ComputeXSectionPerAtom(cut);
280
281 // allow partial integration
282 if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
283
284 cross *= Z*Z*bremFactor;
285
286 return cross;
287}
288
289//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290
291G4double G4eBremParametrizedModel::ComputeXSectionPerAtom(G4double cut)
292{
293 G4double cross = 0.0;
294
295 // number of intervals and integration step
296 G4double vcut = G4Log(cut/totalEnergy);
298 G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
299 // n=1; // integration test
300 G4double delta = (vmax - vcut)/G4double(n);
301
302 G4double e0 = vcut;
303 G4double xs;
304
305 // integration
306 for(G4int l=0; l<n; l++) {
307
308 for(G4int i=0; i<8; i++) {
309
310 G4double eg = G4Exp(e0 + xgi[i]*delta)*totalEnergy;
311
312 xs = ComputeDXSectionPerAtom(eg);
313
314 cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
315 }
316 e0 += delta;
317 }
318
319 cross *= delta;
320
321 return cross;
322}
323
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325
326// compute the value of the screening function 3*PHI1 - PHI2
327
328G4double G4eBremParametrizedModel::ScreenFunction1(G4double ScreenVariable)
329{
330 G4double screenVal;
331
332 if (ScreenVariable > 1.)
333 screenVal = 42.24 - 8.368*G4Log(ScreenVariable+0.952);
334 else
335 screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
336
337 return screenVal;
338}
339
340//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
341
342// compute the value of the screening function 1.5*PHI1 - 0.5*PHI2
343
344G4double G4eBremParametrizedModel::ScreenFunction2(G4double ScreenVariable)
345{
346 G4double screenVal;
347
348 if (ScreenVariable > 1.)
349 screenVal = 42.24 - 8.368*G4Log(ScreenVariable+0.952);
350 else
351 screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
352
353 return screenVal;
354}
355
356//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357
358// Parametrized cross section
359G4double G4eBremParametrizedModel::ComputeParametrizedDXSectionPerAtom(
360 G4double kineticEnergy,
361 G4double gammaEnergy, G4double Z)
362{
363 SetCurrentElement(Z);
364 G4double FZ = lnZ* (4.- 0.55*lnZ);
365 G4double Z3 = z13;
366 G4double ZZ = z13*nist->GetZ13(G4lrint(Z)+1);
367
368 totalEnergy = kineticEnergy + electron_mass_c2;
369
370 // G4double x, epsil, greject, migdal, grejmax, q;
371 G4double epsil, greject;
372 G4double U = G4Log(kineticEnergy/electron_mass_c2);
373 G4double U2 = U*U;
374
375 // precalculated parameters
376 G4double ah, bh;
377
378 if (kineticEnergy > tlow) {
379
380 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
381 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
382 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
383
384 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
385 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
386 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
387
388 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
389 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
390
391 // limit of the screening variable
392 G4double screenfac =
393 136.*electron_mass_c2/(Z3*totalEnergy);
394
395 epsil = gammaEnergy/totalEnergy; // epsil = x*kineticEnergy/totalEnergy;
396 G4double screenvar = screenfac*epsil/(1.0-epsil);
397 G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
398 G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
399
400
401 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.; // 1./(42.392 - FZ);
402
403 std::cout << " yy = "<<epsil<<std::endl;
404 std::cout << " F1/(...) "<<F1/(42.392 - FZ)<<std::endl;
405 std::cout << " F2/(...) "<<F2/(42.392 - FZ)<<std::endl;
406 std::cout << " (42.392 - FZ) " << (42.392 - FZ) <<std::endl;
407
408 } else {
409
410 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
411 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
412 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
413
414 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
415 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
416 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
417
418 ah = al0 + al1*U + al2*U2;
419 bh = bl0 + bl1*U + bl2*U2;
420
421 G4double x=gammaEnergy/kineticEnergy;
422 greject=(1. + x* (ah + bh*x));
423
424 /*
425 // Compute the maximum of the rejection function
426 grejmax = max(1. + xmin* (ah + bh*xmin), 1.+ah+bh);
427 G4double xm = -ah/(2.*bh);
428 if ( xmin < xm && xm < xmax) grejmax = max(grejmax, 1.+ xm* (ah + bh*xm));
429 */
430 }
431
432 return greject;
433}
434
435//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
436
437G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
438{
439
440 if(gammaEnergy < 0.0) { return 0.0; }
441
442 G4double y = gammaEnergy/totalEnergy;
443
444 G4double main=0.;
445 //secondTerm=0.;
446
447 // ** form factors complete screening case **
448 // only valid for high energies (and if LPM suppression does not play a role)
449 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
450 // secondTerm = (1.-y)/12.*(1.+1./currentZ);
451
452 std::cout<<" F1(0) "<<ScreenFunction1(0.) <<std::endl;
453 std::cout<<" F1(0) "<<ScreenFunction2(0.) <<std::endl;
454 std::cout<<"Ekin = "<<kinEnergy<<std::endl;
455 std::cout<<"Z = "<<currentZ<<std::endl;
456 std::cout<<"main = "<<main<<std::endl;
457 std::cout<<" y = "<<y<<std::endl;
458 std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb) <<std::endl;
459
460 G4double main2 = ComputeParametrizedDXSectionPerAtom(kinEnergy,gammaEnergy,currentZ);
461 std::cout<<"main2 = "<<main2<<std::endl;
462 std::cout<<"main2tot = "<<main2 * ( (Fel-fCoulomb) + Finel/currentZ )/(Fel-fCoulomb);
463
464 G4double cross = main2; //main+secondTerm;
465 return cross;
466}
467
468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469
471 std::vector<G4DynamicParticle*>* vdp,
472 const G4MaterialCutsCouple* couple,
473 const G4DynamicParticle* dp,
474 G4double cutEnergy,
475 G4double maxEnergy)
476{
477 G4double kineticEnergy = dp->GetKineticEnergy();
478 if(kineticEnergy < lowKinEnergy) { return; }
479 G4double cut = std::min(cutEnergy, kineticEnergy);
480 G4double emax = std::min(maxEnergy, kineticEnergy);
481 if(cut >= emax) { return; }
482
483 SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
484
485 const G4Element* elm = SelectTargetAtom(couple,particle,kineticEnergy,
486 dp->GetLogKineticEnergy(),cut,emax);
487 SetCurrentElement(elm->GetZ());
488
489 kinEnergy = kineticEnergy;
490 totalEnergy = kineticEnergy + particleMass;
492
493 G4double xmin = G4Log(cut*cut + densityCorr);
494 G4double xmax = G4Log(emax*emax + densityCorr);
495 G4double gammaEnergy, f, x;
496
497 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
498
499 do {
500 x = G4Exp(xmin + rndmEngine->flat()*(xmax - xmin)) - densityCorr;
501 if(x < 0.0) x = 0.0;
502 gammaEnergy = sqrt(x);
503 f = ComputeDXSectionPerAtom(gammaEnergy);
504
505 if ( f > fMax ) {
506 G4cout << "### G4eBremParametrizedModel Warning: Majoranta exceeded! "
507 << f << " > " << fMax
508 << " Egamma(MeV)= " << gammaEnergy
509 << " E(mEV)= " << kineticEnergy
510 << G4endl;
511 }
512
513 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
514 } while (f < fMax*rndmEngine->flat());
515
516 //
517 // angles of the emitted gamma. ( Z - axis along the parent particle)
518 // use general interface
519 //
520 G4ThreeVector gammaDirection =
523 couple->GetMaterial());
524
525 // create G4DynamicParticle object for the Gamma
526 auto gamma = new G4DynamicParticle(theGamma,gammaDirection, gammaEnergy);
527 vdp->push_back(gamma);
528
529 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
530 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
531 - gammaEnergy*gammaDirection).unit();
532
533 // energy of primary
534 G4double finalE = kineticEnergy - gammaEnergy;
535
536 // stop tracking and create new secondary instead of primary
537 if(gammaEnergy > SecondaryThreshold()) {
540 auto el =
542 direction, finalE);
543 vdp->push_back(el);
544
545 // continue tracking
546 } else {
549 }
550}
551
552//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553
554
std::vector< const G4Element * > G4ElementVector
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double G4Log(G4double x)
Definition: G4Log.hh:227
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
virtual double flat()=0
const G4ThreeVector & GetMomentumDirection() const
G4double GetLogKineticEnergy() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4double GetZ() const
Definition: G4Element.hh:131
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:185
size_t GetNumberOfElements() const
Definition: G4Material.hh:181
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:211
G4double GetElectronDensity() const
Definition: G4Material.hh:212
G4double GetZ13(G4double Z) const
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:831
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:600
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:823
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
void SetCurrentElement(const G4Element *)
Definition: G4VEmModel.hh:493
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:753
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:607
const G4Element * SelectTargetAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:577
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:139
G4double SecondaryThreshold() const
Definition: G4VEmModel.hh:669
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:109
void ProposeTrackStatus(G4TrackStatus status)
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
static const G4double wgi[8]
~G4eBremParametrizedModel() override
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *) override
G4ParticleChangeForLoss * fParticleChange
const G4ParticleDefinition * particle
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double cutEnergy, G4double maxEnergy) override
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double tkin, G4double Z, G4double, G4double cutEnergy, G4double maxEnergy=DBL_MAX) override
G4eBremParametrizedModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="eBremParam")
static const G4double xgi[8]
void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double) override
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
T max(const T t1, const T t2)
brief Return the largest of the two arguments
int G4lrint(double ad)
Definition: templates.hh:134