Geant4-11
G4ChordFinder.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// G4ChordFinder implementation
27//
28// Author: J.Apostolakis - Design and implementation - 25.02.1997
29// -------------------------------------------------------------------
30
31#include <iomanip>
32
33#include "G4ChordFinder.hh"
34#include "G4SystemOfUnits.hh"
35#include "G4MagneticField.hh"
36#include "G4Mag_UsualEqRhs.hh"
38// #include "G4ClassicalRK4.hh"
39// #include "G4CashKarpRKF45.hh"
40// #include "G4NystromRK4.hh"
41// #include "G4BogackiShampine23.hh"
42// #include "G4BogackiShampine45.hh"
43
44#include "G4DormandPrince745.hh"
45
46// New templated stepper(s) -- avoid virtual calls to equation rhs
47#include "G4TDormandPrince45.hh"
48
49// FSAL type driver / steppers -----
52#include "G4RK547FEq1.hh"
53// #include "G4RK547FEq2.hh"
54// #include "G4RK547FEq3.hh"
55// #include "G4FSALBogackiShampine45.hh"
56// #include "G4FSALDormandPrince745.hh"
57
58// Templated type drivers -----
61
62#include "G4HelixHeum.hh"
64
66
67#include <cassert>
68
69static G4bool gVerboseCtor = false; // true;
70
71// ..........................................................................
72
74 : fDefaultDeltaChord(0.25 * mm), fIntgrDriver(pIntegrationDriver)
75{
76 // Simple constructor -- it does not create equation
77 if( gVerboseCtor )
78 G4cout << "G4ChordFinder: Simple constructor -- it uses pre-existing driver." << G4endl;
79
80 fDeltaChord = fDefaultDeltaChord; // Parameters
81}
82
83
84// ..........................................................................
85
87 G4double stepMinimum,
88 G4MagIntegratorStepper* pItsStepper,
89 G4int stepperDriverId )
90 : fDefaultDeltaChord(0.25 * mm)
91{
92 // Construct the Chord Finder
93 // by creating in inverse order the Driver, the Stepper and EqRhs ...
94 constexpr G4int nVar6 = 6; // Components integrated in Usual Equation of motion
95
96 fDeltaChord = fDefaultDeltaChord; // Parameters
97
98 // stepperDriverId = 2;
99
100 G4bool useFSALstepper= (stepperDriverId == 1);
101 G4bool useTemplatedStepper= (stepperDriverId == 2);
102 G4bool useRegularStepper = (stepperDriverId == 3);
103 // G4bool useBFieldDriver = !useRegularStepper && !useFSALstepper && !useTemplatedStepper;
104
105 // G4bool useRegularStepper = !stepperDriverId != 3) && !useFSALStepper && !useTemplatedStepper;
106 // If it's not 0, 1 or 2 then 'BFieldDriver' which combines DoPri5 (short) and helix is used.
107
108 using EquationType = G4Mag_UsualEqRhs;
109
110 using TemplatedStepperType =
111 G4TDormandPrince45<EquationType,nVar6>; // 5th order embedded method. High efficiency.
112 const char* TemplatedStepperName =
113 "G4TDormandPrince745 (templated Dormand-Prince45, aka DoPri5): 5th/4th Order 7-stage embedded";
114
115 using RegularStepperType =
116 G4DormandPrince745; // 5th order embedded method. High efficiency.
117 // G4ClassicalRK4; // The old default
118 // G4CashKarpRKF45; // First embedded method in G4
119 // G4BogackiShampine45; // High efficiency 5th order embedded method
120 // G4NystromRK4; // Nystrom stepper 4th order
121 // G4RK547FEq1; // or 2 or 3
122 const char* RegularStepperName =
123 "G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded";
124 // "BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)";
125 // "Nystrom stepper 4th order";
126
127 using NewFsalStepperType = G4DormandPrince745; // Now works -- 2020.10.08
128 // Was G4RK547FEq1; // or 2 or 3
129 const char* NewFSALStepperName =
130 "G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1.";
131
132#ifdef G4DEBUG_FIELD
133 static G4bool verboseDebug = true;
134 if( verboseDebug )
135 {
136 G4cout << "G4ChordFinder 2nd Constructor called. " << G4endl;
137 G4cout << " Arguments: " << G4endl
138 << " - min step = " << stepMinimum << G4endl
139 << " - stepper ptr provided : "
140 << ( pItsStepper==nullptr ? " no " : " yes " ) << G4endl;
141 if( pItsStepper==nullptr )
142 G4cout << " - stepper/driver Id = " << stepperDriverId << " i.e. "
143 << " useFSAL = " << useFSALstepper
144 << " , useTemplated = " << useTemplatedStepper
145 << " , useRegular = " << useRegularStepper
146 << " , useFSAL = " << useFSALstepper
147 << G4endl;
148 }
149#endif
150
151 // useHigherStepper = forceHigherEffiencyStepper || useHigherStepper;
152
153 EquationType* pEquation = new G4Mag_UsualEqRhs(theMagField);
154 fEquation = pEquation;
155
156 // G4MagIntegratorStepper* regularStepper = nullptr;
157 // G4VFSALIntegrationStepper* fsalStepper = nullptr; // for FSAL steppers only
158 // G4MagIntegratorStepper* oldFSALStepper = nullptr;
159
160 G4bool errorInStepperCreation = false;
161
162 std::ostringstream message; // In case of failure, load with description !
163
164 if( pItsStepper != nullptr )
165 {
166 if( gVerboseCtor )
167 G4cout << " G4ChordFinder: Creating G4IntegrationDriver<G4MagIntegratorStepper> with "
168 << " stepMinimum = " << stepMinimum
169 << " numVar= " << pItsStepper->GetNumberOfVariables() << G4endl;
170 // Stepper type is not known - so must use base class G4MagIntegratorStepper
172 stepMinimum, pItsStepper, pItsStepper->GetNumberOfVariables());
173
174 // -- Older:
175 // G4cout << " G4ChordFinder: Creating G4MagInt_Driver with " ...
176 // Type is not known - so must use old class
177 // fIntgrDriver = new G4MagInt_Driver( stepMinimum, pItsStepper,
178 // pItsStepper->GetNumberOfVariables());
179 }
180 else if ( useTemplatedStepper )
181 {
182 if( gVerboseCtor )
183 G4cout << " G4ChordFinder: Creating Templated Stepper of type> "
184 << TemplatedStepperName << G4endl;
185 // RegularStepperType* regularStepper = nullptr; // To check the exception
186 auto templatedStepper = new TemplatedStepperType(pEquation);
187 // *** ******************
188 //
189 // Alternative - for G4NystromRK4:
190 // = new G4NystromRK4(pEquation, 0.1*mm );
191 fRegularStepperOwned = templatedStepper;
192 if( templatedStepper == nullptr )
193 {
194 message << "Templated Stepper instantiation FAILED." << G4endl;
195 message << "G4ChordFinder: Attempted to instantiate "
196 << TemplatedStepperName << " type stepper " << G4endl;
197 errorInStepperCreation = true;
198 }
199 else
200 {
202 stepMinimum, templatedStepper, nVar6 );
203 if( gVerboseCtor )
204 G4cout << " G4ChordFinder: Using G4IntegrationDriver. " << G4endl;
205 }
206
207 }
208 else if ( useRegularStepper ) // Plain stepper -- not double ...
209 {
210 auto regularStepper = new RegularStepperType(pEquation);
211 // *** ******************
212 fRegularStepperOwned = regularStepper;
213
214 if( gVerboseCtor )
215 G4cout << " G4ChordFinder: Creating Driver for regular stepper.";
216
217 if( regularStepper == nullptr )
218 {
219 message << "Regular Stepper instantiation FAILED." << G4endl;
220 message << "G4ChordFinder: Attempted to instantiate "
221 << RegularStepperName << " type stepper " << G4endl;
222 errorInStepperCreation = true;
223 }
224 else
225 {
226 auto dp5= dynamic_cast<G4DormandPrince745*>(regularStepper);
227 if( dp5 ) {
229 stepMinimum, dp5, nVar6 );
230 if( gVerboseCtor )
231 G4cout << " Using InterpolationDriver<DoPri5> " << G4endl;
232 } else {
234 stepMinimum, regularStepper, nVar6 );
235 if( gVerboseCtor )
236 G4cout << " Using IntegrationDriver<DoPri5> " << G4endl;
237 }
238 }
239 }
240 else if ( !useFSALstepper )
241 {
242 auto regularStepper = new G4DormandPrince745(pEquation);
243 // *** ******************
244 //
245 fRegularStepperOwned = regularStepper;
246
247 {
248 using SmallStepDriver = G4InterpolationDriver<G4DormandPrince745>;
249 using LargeStepDriver = G4IntegrationDriver<G4HelixHeum>;
250
251 fLongStepper = std::unique_ptr<G4HelixHeum>(new G4HelixHeum(pEquation));
252
254 std::unique_ptr<SmallStepDriver>(new SmallStepDriver(stepMinimum,
255 regularStepper, regularStepper->GetNumberOfVariables())),
256 std::unique_ptr<LargeStepDriver>(new LargeStepDriver(stepMinimum,
257 fLongStepper.get(), regularStepper->GetNumberOfVariables())) );
258
259 if( fIntgrDriver == nullptr)
260 {
261 message << "Using G4BFieldIntegrationDriver with "
262 << RegularStepperName << " type stepper " << G4endl;
263 message << "Driver instantiation FAILED." << G4endl;
264 G4Exception("G4ChordFinder::G4ChordFinder()",
265 "GeomField1001", JustWarning, message);
266 }
267 }
268 }
269 else
270 {
271 auto fsalStepper= new NewFsalStepperType(pEquation);
272 // *** ******************
273 fNewFSALStepperOwned = fsalStepper;
274
275 if( fsalStepper == nullptr )
276 {
277 message << "Stepper instantiation FAILED." << G4endl;
278 message << "Attempted to instantiate "
279 << NewFSALStepperName << " type stepper " << G4endl;
280 G4Exception("G4ChordFinder::G4ChordFinder()",
281 "GeomField1001", JustWarning, message);
282 errorInStepperCreation = true;
283 }
284 else
285 {
286 fIntgrDriver = new
287 G4FSALIntegrationDriver<NewFsalStepperType>(stepMinimum, fsalStepper,
288 fsalStepper->GetNumberOfVariables() );
289 // ==== Create the driver which knows the class type
290
291 if( fIntgrDriver == nullptr )
292 {
293 message << "Using G4FSALIntegrationDriver with stepper type: "
294 << NewFSALStepperName << G4endl;
295 message << "Integration Driver instantiation FAILED." << G4endl;
296 G4Exception("G4ChordFinder::G4ChordFinder()",
297 "GeomField1001", JustWarning, message);
298 }
299 }
300 }
301
302 // -- Main work is now done
303
304 // Now check that no error occured, and report it if one did.
305
306 // To test failure to create driver
307 // delete fIntgrDriver;
308 // fIntgrDriver = nullptr;
309
310 // Detect and report Error conditions
311 //
312 if( errorInStepperCreation || (fIntgrDriver == nullptr ))
313 {
314 std::ostringstream errmsg;
315
316 if( errorInStepperCreation )
317 {
318 errmsg << "ERROR> Failure to create Stepper object." << G4endl
319 << " --------------------------------" << G4endl;
320 }
321 if (fIntgrDriver == nullptr )
322 {
323 errmsg << "ERROR> Failure to create Integration-Driver object."
324 << G4endl
325 << " -------------------------------------------"
326 << G4endl;
327 }
328 const std::string BoolName[2]= { "False", "True" };
329 errmsg << " Configuration: (constructor arguments) " << G4endl
330 << " provided Stepper = " << pItsStepper << G4endl
331 << " stepper/driver Id = " << stepperDriverId << " i.e. "
332 << " useTemplated = " << BoolName[useTemplatedStepper]
333 << " useRegular = " << BoolName[useRegularStepper]
334 << " useFSAL = " << BoolName[useFSALstepper]
335 << " using combo BField Driver = " <<
336 BoolName[ ! (useFSALstepper||useTemplatedStepper
337 || useRegularStepper ) ]
338 << G4endl;
339 errmsg << message.str();
340 errmsg << "Aborting.";
341 G4Exception("G4ChordFinder::G4ChordFinder() - constructor 2",
342 "GeomField0003", FatalException, errmsg);
343 }
344
345 assert( ( pItsStepper != nullptr )
346 || ( fRegularStepperOwned != nullptr )
347 || ( fNewFSALStepperOwned != nullptr )
348 );
349 assert( fIntgrDriver != nullptr );
350}
351
352
353// ......................................................................
354
356{
357 delete fEquation;
360 delete fCachedField;
361 delete fIntgrDriver;
362}
363
364// ...........................................................................
365
367G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
368 const G4FieldTrack& CurveB_PointVelocity,
369 const G4FieldTrack& ApproxCurveV,
370 const G4ThreeVector& CurrentE_Point,
371 const G4ThreeVector& CurrentF_Point,
372 const G4ThreeVector& PointG,
373 G4bool first, G4double eps_step)
374{
375 // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
376 // Use Brent Algorithm (or InvParabolic) when possible.
377 // Given a starting curve point A (CurveA_PointVelocity), curve point B
378 // (CurveB_PointVelocity), a point E which is (generally) not on the curve
379 // and a point F which is on the curve (first approximation), find new
380 // point S on the curve closer to point E.
381 // While advancing towards S utilise 'eps_step' as a measure of the
382 // relative accuracy of each Step.
383
384 G4FieldTrack EndPoint(CurveA_PointVelocity);
385 if(!first) { EndPoint = ApproxCurveV; }
386
387 G4ThreeVector Point_A,Point_B;
388 Point_A=CurveA_PointVelocity.GetPosition();
389 Point_B=CurveB_PointVelocity.GetPosition();
390
391 G4double xa,xb,xc,ya,yb,yc;
392
393 // InverseParabolic. AF Intersects (First Part of Curve)
394
395 if(first)
396 {
397 xa=0.;
398 ya=(PointG-Point_A).mag();
399 xb=(Point_A-CurrentF_Point).mag();
400 yb=-(PointG-CurrentF_Point).mag();
401 xc=(Point_A-Point_B).mag();
402 yc=-(CurrentE_Point-Point_B).mag();
403 }
404 else
405 {
406 xa=0.;
407 ya=(Point_A-CurrentE_Point).mag();
408 xb=(Point_A-CurrentF_Point).mag();
409 yb=(PointG-CurrentF_Point).mag();
410 xc=(Point_A-Point_B).mag();
411 yc=-(Point_B-PointG).mag();
412 if(xb==0.)
413 {
414 EndPoint = ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
415 CurrentE_Point, eps_step);
416 return EndPoint;
417 }
418 }
419
420 const G4double tolerance = 1.e-12;
421 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
422 {
423 ; // What to do for the moment: return the same point as at start
424 // then PropagatorInField will take care
425 }
426 else
427 {
428 G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
429 G4double curve;
430 if(first)
431 {
432 curve=std::abs(EndPoint.GetCurveLength()
433 -ApproxCurveV.GetCurveLength());
434 }
435 else
436 {
437 test_step = test_step - xb;
438 curve=std::abs(EndPoint.GetCurveLength()
439 -CurveB_PointVelocity.GetCurveLength());
440 xb = (CurrentF_Point-Point_B).mag();
441 }
442
443 if(test_step<=0) { test_step=0.1*xb; }
444 if(test_step>=xb) { test_step=0.5*xb; }
445 if(test_step>=curve){ test_step=0.5*curve; }
446
447 if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
448 { // G4VIntersectionLocator
449 test_step=0.5*curve;
450 }
451
452 fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
453
454#ifdef G4DEBUG_FIELD
455 // Printing Brent and Linear Approximation
456 //
457 G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
458 << test_step << " EndPoint = " << EndPoint << G4endl;
459
460 // Test Track
461 //
462 G4FieldTrack TestTrack( CurveA_PointVelocity);
463 TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
464 CurveB_PointVelocity,
465 CurrentE_Point, eps_step );
466 G4cout.precision(14);
467 G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
468 G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
469#endif
470 }
471 return EndPoint;
472}
473
474
475// ...........................................................................
476
478ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
479 const G4FieldTrack& CurveB_PointVelocity,
480 const G4ThreeVector& CurrentE_Point,
481 G4double eps_step)
482{
483 // If r=|AE|/|AB|, and s=true path lenght (AB)
484 // return the point that is r*s along the curve!
485
486 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
487
488 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
489 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
490
491 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
492 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
493
494 G4double ABdist= ChordAB_Vector.mag();
495 G4double curve_length; // A curve length of AB
496 G4double AE_fraction;
497
498 curve_length= CurveB_PointVelocity.GetCurveLength()
499 - CurveA_PointVelocity.GetCurveLength();
500
501 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
502 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
503 {
504#ifdef G4DEBUG_FIELD
505 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
506 << G4endl
507 << " The two points are further apart than the curve length "
508 << G4endl
509 << " Dist = " << ABdist
510 << " curve length = " << curve_length
511 << " relativeDiff = " << (curve_length-ABdist)/ABdist
512 << G4endl;
513 if( curve_length < ABdist * (1. - 10*eps_step) )
514 {
515 std::ostringstream message;
516 message << "Unphysical curve length." << G4endl
517 << "The size of the above difference exceeds allowed limits."
518 << G4endl
519 << "Aborting.";
520 G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
521 FatalException, message);
522 }
523#endif
524 // Take default corrective action: adjust the maximum curve length.
525 // NOTE: this case only happens for relatively straight paths.
526 // curve_length = ABdist;
527 }
528
529 G4double new_st_length;
530
531 if ( ABdist > 0.0 )
532 {
533 AE_fraction = ChordAE_Vector.mag() / ABdist;
534 }
535 else
536 {
537 AE_fraction = 0.5; // Guess .. ?;
538#ifdef G4DEBUG_FIELD
539 G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
540 << " A and B are the same point!" << G4endl
541 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
542 << G4endl;
543#endif
544 }
545
546 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
547 {
548#ifdef G4DEBUG_FIELD
549 G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
550 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
551 << " AE_fraction = " << AE_fraction << G4endl
552 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
553 << " Chord AB length = " << ABdist << G4endl << G4endl;
554 G4cerr << " OK if this condition occurs after a recalculation of 'B'"
555 << G4endl << " Otherwise it is an error. " << G4endl ;
556#endif
557 // This course can now result if B has been re-evaluated,
558 // without E being recomputed (1 July 99).
559 // In this case this is not a "real error" - but it is undesired
560 // and we cope with it by a default corrective action ...
561 //
562 AE_fraction = 0.5; // Default value
563 }
564
565 new_st_length = AE_fraction * curve_length;
566
567 if ( AE_fraction > 0.0 )
568 {
569 fIntgrDriver->AccurateAdvance(Current_PointVelocity,
570 new_st_length, eps_step );
571 //
572 // In this case it does not matter if it cannot advance the full distance
573 }
574
575 // If there was a memory of the step_length actually required at the start
576 // of the integration Step, this could be re-used ...
577
578 G4cout.precision(14);
579
580 return Current_PointVelocity;
581}
582
583// ...........................................................................
584
585std::ostream& operator<<( std::ostream& os, const G4ChordFinder& cf)
586{
587 // Dumping the state of G4ChordFinder
588 os << "State of G4ChordFinder : " << std::endl;
589 os << " delta_chord = " << cf.fDeltaChord;
590 os << " Default d_c = " << cf.fDefaultDeltaChord;
591
592 os << " stats-verbose = " << cf.fStatsVerbose;
593
594 return os;
595}
static G4bool gVerboseCtor
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
static constexpr double perMillion
Definition: G4SIunits.hh:327
static constexpr double mm
Definition: G4SIunits.hh:95
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
double mag() const
G4EquationOfMotion * fEquation
G4double InvParabolic(const G4double xa, const G4double ya, const G4double xb, const G4double yb, const G4double xc, const G4double yc)
std::unique_ptr< G4HelixHeum > fLongStepper
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector &currentEPoint, const G4ThreeVector &currentFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
G4MagIntegratorStepper * fNewFSALStepperOwned
virtual ~G4ChordFinder()
const G4double fDefaultDeltaChord
G4MagIntegratorStepper * fRegularStepperOwned
G4CachedMagneticField * fCachedField
G4double fDeltaChord
G4VIntegrationDriver * fIntgrDriver
G4ChordFinder(G4VIntegrationDriver *pIntegrationDriver)
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
G4int GetNumberOfVariables() const
virtual G4bool AccurateAdvance(G4FieldTrack &track, G4double hstep, G4double eps, G4double hinitial=0)=0
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
T max(const T t1, const T t2)
brief Return the largest of the two arguments