Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MagIntegratorDriver.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 // $Id: G4MagIntegratorDriver.cc 69786 2013-05-15 09:38:51Z gcosmo $
28 //
29 //
30 //
31 // Implementation for class G4MagInt_Driver
32 // Tracking in space dependent magnetic field
33 //
34 // History of major changes:
35 // 8 Nov 01 J. Apostolakis: Respect minimum step in AccurateAdvance
36 // 27 Jul 99 J. Apostolakis: Ensured that AccurateAdvance does not loop
37 // due to very small eps & step size (precision)
38 // 28 Jan 98 W. Wander: Added ability for low order integrators
39 // 7 Oct 96 V. Grichine First version
40 // --------------------------------------------------------------------
41 
42 #include <iomanip>
43 
44 #include "globals.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4GeometryTolerance.hh"
47 #include "G4MagIntegratorDriver.hh"
48 #include "G4FieldTrack.hh"
49 
50 // Stepsize can increase by no more than 5.0
51 // and decrease by no more than 1/10. = 0.1
52 //
53 const G4double G4MagInt_Driver::max_stepping_increase = 5.0;
54 const G4double G4MagInt_Driver::max_stepping_decrease = 0.1;
55 
56 // The (default) maximum number of steps is Base
57 // divided by the order of Stepper
58 //
59 const G4int G4MagInt_Driver::fMaxStepBase = 250; // Was 5000
60 
61 #ifndef G4NO_FIELD_STATISTICS
62 #define G4FLD_STATS 1
63 #endif
64 
65 // ---------------------------------------------------------
66 
67 // Constructor
68 //
69 G4MagInt_Driver::G4MagInt_Driver( G4double hminimum,
70  G4MagIntegratorStepper *pStepper,
71  G4int numComponents,
72  G4int statisticsVerbose)
73  : fSmallestFraction( 1.0e-12 ),
74  fNoIntegrationVariables(numComponents),
75  fMinNoVars(12),
76  fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
77  fStatisticsVerboseLevel(statisticsVerbose),
78  fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
79  fNoInitialSmallSteps(0),
80  fDyerr_max(0.0), fDyerr_mx2(0.0),
81  fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
82  fSumH_sm(0.0), fSumH_lg(0.0),
83  fVerboseLevel(0)
84 {
85  // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
86  // is required. For proper time of flight and spin, fMinNoVars must be 12
87 
88  RenewStepperAndAdjust( pStepper );
89  fMinimumStep= hminimum;
90  fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
91 #ifdef G4DEBUG_FIELD
92  fVerboseLevel=2;
93 #endif
94 
95  if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
96  {
97  G4cout << "MagIntDriver version: Accur-Adv: "
98  << "invE_nS, QuickAdv-2sqrt with Statistics "
99 #ifdef G4FLD_STATS
100  << " enabled "
101 #else
102  << " disabled "
103 #endif
104  << G4endl;
105  }
106 }
107 
108 // ---------------------------------------------------------
109 
110 // Destructor
111 //
112 G4MagInt_Driver::~G4MagInt_Driver()
113 {
114  if( fStatisticsVerboseLevel > 1 )
115  {
116  PrintStatisticsReport();
117  }
118 }
119 
120 // To add much printing for debugging purposes, uncomment the following
121 // and set verbose level to 1 or higher value !
122 // #define G4DEBUG_FIELD 1
123 
124 // ---------------------------------------------------------
125 
126 G4bool
127 G4MagInt_Driver::AccurateAdvance(G4FieldTrack& y_current,
128  G4double hstep,
129  G4double eps,
130  G4double hinitial )
131 {
132  // Runge-Kutta driver with adaptive stepsize control. Integrate starting
133  // values at y_current over hstep x2 with accuracy eps.
134  // On output ystart is replaced by values at the end of the integration
135  // interval. RightHandSide is the right-hand side of ODE system.
136  // The source is similar to odeint routine from NRC p.721-722 .
137 
138  G4int nstp, i, no_warnings=0;
139  G4double x, hnext, hdid, h;
140 
141 #ifdef G4DEBUG_FIELD
142  static G4int dbg=1;
143  static G4int nStpPr=50; // For debug printing of long integrations
144  G4double ySubStepStart[G4FieldTrack::ncompSVEC];
145  G4FieldTrack yFldTrkStart(y_current);
146 #endif
147 
148  G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
149  G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC];
150  G4double x1, x2;
151  G4bool succeeded = true, lastStepSucceeded;
152 
153  G4double startCurveLength;
154 
155  G4int noFullIntegr=0, noSmallIntegr = 0 ;
156  static G4int noGoodSteps =0 ; // Bad = chord > curve-len
157  const G4int nvar= fNoVars;
158 
159  G4FieldTrack yStartFT(y_current);
160 
161  // Ensure that hstep > 0
162  //
163  if( hstep <= 0.0 )
164  {
165  if(hstep==0.0)
166  {
167  std::ostringstream message;
168  message << "Proposed step is zero; hstep = " << hstep << " !";
169  G4Exception("G4MagInt_Driver::AccurateAdvance()",
170  "GeomField1001", JustWarning, message);
171  return succeeded;
172  }
173  else
174  {
175  std::ostringstream message;
176  message << "Invalid run condition." << G4endl
177  << "Proposed step is negative; hstep = " << hstep << "." << G4endl
178  << "Requested step cannot be negative! Aborting event.";
179  G4Exception("G4MagInt_Driver::AccurateAdvance()",
180  "GeomField0003", EventMustBeAborted, message);
181  return false;
182  }
183  }
184 
185  y_current.DumpToArray( ystart );
186 
187  startCurveLength= y_current.GetCurveLength();
188  x1= startCurveLength;
189  x2= x1 + hstep;
190 
191  if ( (hinitial > 0.0) && (hinitial < hstep)
192  && (hinitial > perMillion * hstep) )
193  {
194  h = hinitial;
195  }
196  else // Initial Step size "h" defaults to the full interval
197  {
198  h = hstep;
199  }
200 
201  x = x1;
202 
203  for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
204 
205  G4bool lastStep= false;
206  nstp=1;
207 
208  do
209  {
210  G4ThreeVector StartPos( y[0], y[1], y[2] );
211 
212 #ifdef G4DEBUG_FIELD
213  G4double xSubStepStart= x;
214  for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
215  yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
216  yFldTrkStart.SetCurveLength(x);
217 #endif
218 
219  // Old method - inline call to Equation of Motion
220  // pIntStepper->RightHandSide( y, dydx );
221  // New method allows to cache field, or state (eg momentum magnitude)
222  pIntStepper->ComputeRightHandSide( y, dydx );
223  fNoTotalSteps++;
224 
225  // Perform the Integration
226  //
227  if( h > fMinimumStep )
228  {
229  OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
230  //--------------------------------------
231  lastStepSucceeded= (hdid == h);
232 #ifdef G4DEBUG_FIELD
233  if (dbg>2) {
234  PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
235  }
236 #endif
237  }
238  else
239  {
240  G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
241  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
242  G4double dchord_step, dyerr, dyerr_len; // What to do with these ?
243  yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
244  yFldTrk.SetCurveLength( x );
245 
246  QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
247  //-----------------------------------------------------
248 
249  yFldTrk.DumpToArray(y);
250 
251 #ifdef G4FLD_STATS
252  fNoSmallSteps++;
253  if ( dyerr_len > fDyerr_max) { fDyerr_max= dyerr_len; }
254  fDyerrPos_smTot += dyerr_len;
255  fSumH_sm += h; // Length total for 'small' steps
256  if (nstp<=1) { fNoInitialSmallSteps++; }
257 #endif
258 #ifdef G4DEBUG_FIELD
259  if (dbg>1)
260  {
261  if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
262  G4cout << "Another sub-min step, no " << fNoSmallSteps
263  << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
264  PrintStatus( ySubStepStart, x1, y, x, h, nstp); // Only this
265  G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
266  << " epsilon= " << eps << " hstep= " << hstep
267  << " h= " << h << " hmin= " << fMinimumStep << G4endl;
268  }
269 #endif
270  if( h == 0.0 )
271  {
272  G4Exception("G4MagInt_Driver::AccurateAdvance()",
273  "GeomField0003", FatalException,
274  "Integration Step became Zero!");
275  }
276  dyerr = dyerr_len / h;
277  hdid= h;
278  x += hdid;
279 
280  // Compute suggested new step
281  hnext= ComputeNewStepSize( dyerr/eps, h);
282 
283  // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
284  lastStepSucceeded= (dyerr<= eps);
285  }
286 
287  if (lastStepSucceeded) { noFullIntegr++; }
288  else { noSmallIntegr++; }
289 
290  G4ThreeVector EndPos( y[0], y[1], y[2] );
291 
292 #ifdef G4DEBUG_FIELD
293  if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
294  {
295  if( nstp==nStpPr ) { G4cout << "***** Many steps ****" << G4endl; }
296  G4cout << "MagIntDrv: " ;
297  G4cout << "hdid=" << std::setw(12) << hdid << " "
298  << "hnext=" << std::setw(12) << hnext << " "
299  << "hstep=" << std::setw(12) << hstep << " (requested) "
300  << G4endl;
301  PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
302  }
303 #endif
304 
305  // Check the endpoint
306  G4double endPointDist= (EndPos-StartPos).mag();
307  if ( endPointDist >= hdid*(1.+perMillion) )
308  {
309  fNoBadSteps++;
310 
311  // Issue a warning only for gross differences -
312  // we understand how small difference occur.
313  if ( endPointDist >= hdid*(1.+perThousand) )
314  {
315 #ifdef G4DEBUG_FIELD
316  if (dbg)
317  {
318  WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
319  G4cerr << " Total steps: bad " << fNoBadSteps
320  << " good " << noGoodSteps << " current h= " << hdid
321  << G4endl;
322  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
323  }
324 #endif
325  no_warnings++;
326  }
327  }
328  else
329  {
330  noGoodSteps ++;
331  }
332 // #endif
333 
334  // Avoid numerous small last steps
335  if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
336  {
337  // No more integration -- the next step will not happen
338  lastStep = true;
339  }
340  else
341  {
342  // Check the proposed next stepsize
343  if(std::fabs(hnext) <= Hmin())
344  {
345 #ifdef G4DEBUG_FIELD
346  // If simply a very small interval is being integrated, do not warn
347  if( (x < x2 * (1-eps) ) && // The last step can be small: OK
348  (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
349  {
350  if(dbg>0)
351  {
352  WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
353  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
354  }
355  no_warnings++;
356  }
357 #endif
358  // Make sure that the next step is at least Hmin.
359  h = Hmin();
360  }
361  else
362  {
363  h = hnext;
364  }
365 
366  // Ensure that the next step does not overshoot
367  if ( x+h > x2 )
368  { // When stepsize overshoots, decrease it!
369  h = x2 - x ; // Must cope with difficult rounding-error
370  } // issues if hstep << x2
371 
372  if ( h == 0.0 )
373  {
374  // Cannot progress - accept this as last step - by default
375  lastStep = true;
376 #ifdef G4DEBUG_FIELD
377  if (dbg>2)
378  {
379  int prec= G4cout.precision(12);
380  G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
381  << G4endl
382  << " Integration step 'h' became "
383  << h << " due to roundoff. " << G4endl
384  << " Calculated as difference of x2= "<< x2 << " and x=" << x
385  << " Forcing termination of advance." << G4endl;
386  G4cout.precision(prec);
387  }
388 #endif
389  }
390  }
391  } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
392  // Have we reached the end ?
393  // --> a better test might be x-x2 > an_epsilon
394 
395  succeeded= (x>=x2); // If it was a "forced" last step
396 
397  for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
398 
399  // Put back the values.
400  y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
401  y_current.SetCurveLength( x );
402 
403  if(nstp > fMaxNoSteps)
404  {
405  no_warnings++;
406  succeeded = false;
407 #ifdef G4DEBUG_FIELD
408  if (dbg)
409  {
410  WarnTooManyStep( x1, x2, x ); // Issue WARNING
411  PrintStatus( yEnd, x1, y, x, hstep, -nstp);
412  }
413 #endif
414  }
415 
416 #ifdef G4DEBUG_FIELD
417  if( dbg && no_warnings )
418  {
419  G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
420  PrintStatus( yEnd, x1, y, x, hstep, nstp);
421  }
422 #endif
423 
424  return succeeded;
425 } // end of AccurateAdvance ...........................
426 
427 // ---------------------------------------------------------
428 
429 void
430 G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep,
431  G4double h, G4double xDone,
432  G4int nstp)
433 {
434  static G4int noWarningsIssued =0;
435  const G4int maxNoWarnings = 10; // Number of verbose warnings
436  std::ostringstream message;
437  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
438  {
439  message << "The stepsize for the next iteration, " << hnext
440  << ", is too small - in Step number " << nstp << "." << G4endl
441  << "The minimum for the driver is " << Hmin() << G4endl
442  << "Requested integr. length was " << hstep << " ." << G4endl
443  << "The size of this sub-step was " << h << " ." << G4endl
444  << "The integrations has already gone " << xDone;
445  }
446  else
447  {
448  message << "Too small 'next' step " << hnext
449  << ", step-no: " << nstp << G4endl
450  << ", this sub-step: " << h
451  << ", req_tot_len: " << hstep
452  << ", done: " << xDone << ", min: " << Hmin();
453  }
454  G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
455  JustWarning, message);
456  noWarningsIssued++;
457 }
458 
459 // ---------------------------------------------------------
460 
461 void
462 G4MagInt_Driver::WarnTooManyStep( G4double x1start,
463  G4double x2end,
464  G4double xCurrent)
465 {
466  std::ostringstream message;
467  message << "The number of steps used in the Integration driver"
468  << " (Runge-Kutta) is too many." << G4endl
469  << "Integration of the interval was not completed !" << G4endl
470  << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
471  << " % fraction of it was done.";
472  G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
473  JustWarning, message);
474 }
475 
476 // ---------------------------------------------------------
477 
478 void
479 G4MagInt_Driver::WarnEndPointTooFar (G4double endPointDist,
480  G4double h ,
481  G4double eps,
482  G4int dbg)
483 {
484  static G4double maxRelError=0.0;
485  G4bool isNewMax, prNewMax;
486 
487  isNewMax = endPointDist > (1.0 + maxRelError) * h;
488  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
489  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
490 
491  if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
492  && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
493  {
494  static G4int noWarnings = 0;
495  std::ostringstream message;
496  if( (noWarnings ++ < 10) || (dbg>2) )
497  {
498  message << "The integration produced an end-point which " << G4endl
499  << "is further from the start-point than the curve length."
500  << G4endl;
501  }
502  message << " Distance of endpoints = " << endPointDist
503  << ", curve length = " << h << G4endl
504  << " Difference (curveLen-endpDist)= " << (h - endPointDist)
505  << ", relative = " << (h-endPointDist) / h
506  << ", epsilon = " << eps;
507  G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
508  JustWarning, message);
509  }
510 }
511 
512 // ---------------------------------------------------------
513 
514 void
515 G4MagInt_Driver::OneGoodStep( G4double y[], // InOut
516  const G4double dydx[],
517  G4double& x, // InOut
518  G4double htry,
519  G4double eps_rel_max,
520  G4double& hdid, // Out
521  G4double& hnext ) // Out
522 
523 // Driver for one Runge-Kutta Step with monitoring of local truncation error
524 // to ensure accuracy and adjust stepsize. Input are dependent variable
525 // array y[0,...,5] and its derivative dydx[0,...,5] at the
526 // starting value of the independent variable x . Also input are stepsize
527 // to be attempted htry, and the required accuracy eps. On output y and x
528 // are replaced by their new values, hdid is the stepsize that was actually
529 // accomplished, and hnext is the estimated next stepsize.
530 // This is similar to the function rkqs from the book:
531 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
532 // Edition, by William H. Press, Saul A. Teukolsky, William T.
533 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
534 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
535 
536 {
537  G4double errmax_sq;
538  G4double h, htemp, xnew ;
539 
540  G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
541 
542  h = htry ; // Set stepsize to the initial trial value
543 
544  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
545 
546  G4double errpos_sq=0.0; // square of displacement error
547  G4double errvel_sq=0.0; // square of momentum vector difference
548  G4double errspin_sq=0.0; // square of spin vector difference
549 
550  G4int iter;
551 
552  static G4int tot_no_trials=0;
553  const G4int max_trials=100;
554 
555  G4ThreeVector Spin(y[9],y[10],y[11]);
556  G4bool hasSpin= (Spin.mag2() > 0.0);
557 
558  for (iter=0; iter<max_trials ;iter++)
559  {
560  tot_no_trials++;
561  pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
562  // *******
563  G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
564  G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
565 
566  // Evaluate accuracy
567  //
568  errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
569  errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
570 
571  // Accuracy for momentum
572  errvel_sq = (sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) )
573  / (sqr(y[3]) + sqr(y[4]) + sqr(y[5]) );
574  errvel_sq *= inv_eps_vel_sq;
575  errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
576 
577  if( hasSpin )
578  {
579  // Accuracy for spin
580  errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
581  / ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
582  errspin_sq *= inv_eps_vel_sq;
583  errmax_sq = std::max( errmax_sq, errspin_sq );
584  }
585 
586  if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
587 
588  // Step failed; compute the size of retrial Step.
589  htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
590 
591  if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
592  else { h = 0.1*h; } // reduce stepsize, but no more
593  // than a factor of 10
594  xnew = x + h;
595  if(xnew == x)
596  {
597  G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
598  << " Stepsize underflow in Stepper " << G4endl ;
599  G4cerr << " Step's start x=" << x << " and end x= " << xnew
600  << " are equal !! " << G4endl
601  <<" Due to step-size= " << h
602  << " . Note that input step was " << htry << G4endl;
603  break;
604  }
605  }
606 
607 #ifdef G4FLD_STATS
608  // Sum of squares of position error // and momentum dir (underestimated)
609  fSumH_lg += h;
610  fDyerrPos_lgTot += errpos_sq;
611  fDyerrVel_lgTot += errvel_sq * h * h;
612 #endif
613 
614  // Compute size of next Step
615  if (errmax_sq > errcon*errcon)
616  {
617  hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
618  }
619  else
620  {
621  hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
622  }
623  x += (hdid = h);
624 
625  for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
626 
627  return;
628 } // end of OneGoodStep .............................
629 
630 //----------------------------------------------------------------------
631 
632 // QuickAdvance just tries one Step - it does not ensure accuracy
633 //
634 G4bool G4MagInt_Driver::QuickAdvance(
635  G4FieldTrack& y_posvel, // INOUT
636  const G4double dydx[],
637  G4double hstep, // In
638  G4double& dchord_step,
639  G4double& dyerr_pos_sq,
640  G4double& dyerr_mom_rel_sq )
641 {
642  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
643  FatalException, "Not yet implemented.");
644 
645  // Use the parameters of this method, to please compiler
646  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
647  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
648  return true;
649 }
650 
651 //----------------------------------------------------------------------
652 
653 G4bool G4MagInt_Driver::QuickAdvance(
654  G4FieldTrack& y_posvel, // INOUT
655  const G4double dydx[],
656  G4double hstep, // In
657  G4double& dchord_step,
658  G4double& dyerr )
659 {
660  G4double dyerr_pos_sq, dyerr_mom_rel_sq;
661  G4double yerr_vec[G4FieldTrack::ncompSVEC],
662  yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC];
663  G4double s_start;
664  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
665 
666  static G4int no_call=0;
667  no_call ++;
668 
669  // Move data into array
670  y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
671  s_start = y_posvel.GetCurveLength();
672 
673  // Do an Integration Step
674  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
675  // *******
676 
677  // Estimate curve-chord distance
678  dchord_step= pIntStepper-> DistChord();
679  // *********
680 
681  // Put back the values. yarrout ==> y_posvel
682  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
683  y_posvel.SetCurveLength( s_start + hstep );
684 
685 #ifdef G4DEBUG_FIELD
686  if(fVerboseLevel>2)
687  {
688  G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
689  PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
690  }
691 #endif
692 
693  // A single measure of the error
694  // TO-DO : account for energy, spin, ... ?
695  vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
696  inv_vel_mag_sq = 1.0 / vel_mag_sq;
697  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
698  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
699  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
700 
701  // Calculate also the change in the momentum squared also ???
702  // G4double veloc_square = y_posvel.GetVelocity().mag2();
703  // ...
704 
705 #ifdef RETURN_A_NEW_STEP_LENGTH
706  // The following step cannot be done here because "eps" is not known.
707  dyerr_len = std::sqrt( dyerr_len_sq );
708  dyerr_len_sq /= eps ;
709 
710  // Look at the velocity deviation ?
711  // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
712 
713  // Set suggested new step
714  hstep= ComputeNewStepSize( dyerr_len, hstep);
715 #endif
716 
717  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
718  {
719  dyerr = std::sqrt(dyerr_pos_sq);
720  }
721  else
722  {
723  // Scale it to the current step size - for now
724  dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
725  }
726 
727  return true;
728 }
729 
730 // --------------------------------------------------------------------------
731 
732 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
733 G4bool G4MagInt_Driver::QuickAdvance(
734  G4double yarrin[], // In
735  const G4double dydx[],
736  G4double hstep, // In
737  G4double yarrout[],
738  G4double& dchord_step,
739  G4double& dyerr ) // In length
740 {
741  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
742  FatalException, "Not yet implemented.");
743  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
744  yarrout[0]= yarrin[0];
745 }
746 #endif
747 
748 // --------------------------------------------------------------------------
749 
750 // This method computes new step sizes - but does not limit changes to
751 // within certain factors
752 //
753 G4double
754 G4MagInt_Driver::ComputeNewStepSize(
755  G4double errMaxNorm, // max error (normalised)
756  G4double hstepCurrent) // current step size
757 {
758  G4double hnew;
759 
760  // Compute size of next Step for a failed step
761  if(errMaxNorm > 1.0 )
762  {
763  // Step failed; compute the size of retrial Step.
764  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
765  } else if(errMaxNorm > 0.0 ) {
766  // Compute size of next Step for a successful step
767  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
768  } else {
769  // if error estimate is zero (possible) or negative (dubious)
770  hnew = max_stepping_increase * hstepCurrent;
771  }
772 
773  return hnew;
774 }
775 
776 // ---------------------------------------------------------------------------
777 
778 // This method computes new step sizes limiting changes within certain factors
779 //
780 // It shares its logic with AccurateAdvance.
781 // They are kept separate currently for optimisation.
782 //
783 G4double
784 G4MagInt_Driver::ComputeNewStepSize_WithinLimits(
785  G4double errMaxNorm, // max error (normalised)
786  G4double hstepCurrent) // current step size
787 {
788  G4double hnew;
789 
790  // Compute size of next Step for a failed step
791  if (errMaxNorm > 1.0 )
792  {
793  // Step failed; compute the size of retrial Step.
794  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
795 
796  if (hnew < max_stepping_decrease*hstepCurrent)
797  {
798  hnew = max_stepping_decrease*hstepCurrent ;
799  // reduce stepsize, but no more
800  // than this factor (value= 1/10)
801  }
802  }
803  else
804  {
805  // Compute size of next Step for a successful step
806  if (errMaxNorm > errcon)
807  { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
808  else // No more than a factor of 5 increase
809  { hnew = max_stepping_increase * hstepCurrent; }
810  }
811  return hnew;
812 }
813 
814 // ---------------------------------------------------------------------------
815 
816 void G4MagInt_Driver::PrintStatus( const G4double* StartArr,
817  G4double xstart,
818  const G4double* CurrentArr,
819  G4double xcurrent,
820  G4double requestStep,
821  G4int subStepNo)
822  // Potentially add as arguments:
823  // <dydx> - as Initial Force
824  // stepTaken(hdid) - last step taken
825  // nextStep (hnext) - proposal for size
826 {
827  G4FieldTrack StartFT(G4ThreeVector(0,0,0),
828  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
829  G4FieldTrack CurrentFT (StartFT);
830 
831  StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
832  StartFT.SetCurveLength( xstart);
833  CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
834  CurrentFT.SetCurveLength( xcurrent );
835 
836  PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
837 }
838 
839 // ---------------------------------------------------------------------------
840 
841 void G4MagInt_Driver::PrintStatus(
842  const G4FieldTrack& StartFT,
843  const G4FieldTrack& CurrentFT,
844  G4double requestStep,
845  G4int subStepNo)
846 {
847  G4int verboseLevel= fVerboseLevel;
848  static G4int noPrecision= 5;
849  G4int oldPrec= G4cout.precision(noPrecision);
850  // G4cout.setf(ios_base::fixed,ios_base::floatfield);
851 
852  const G4ThreeVector StartPosition= StartFT.GetPosition();
853  const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
854  const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
855  const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
856 
857  G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
858 
859  G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
860  G4double subStepSize = step_len;
861 
862  if( (subStepNo <= 1) || (verboseLevel > 3) )
863  {
864  subStepNo = - subStepNo; // To allow printing banner
865 
866  G4cout << std::setw( 6) << " " << std::setw( 25)
867  << " G4MagInt_Driver: Current Position and Direction" << " "
868  << G4endl;
869  G4cout << std::setw( 5) << "Step#" << " "
870  << std::setw( 7) << "s-curve" << " "
871  << std::setw( 9) << "X(mm)" << " "
872  << std::setw( 9) << "Y(mm)" << " "
873  << std::setw( 9) << "Z(mm)" << " "
874  << std::setw( 8) << " N_x " << " "
875  << std::setw( 8) << " N_y " << " "
876  << std::setw( 8) << " N_z " << " "
877  << std::setw( 8) << " N^2-1 " << " "
878  << std::setw(10) << " N(0).N " << " "
879  << std::setw( 7) << "KinEner " << " "
880  << std::setw(12) << "Track-l" << " " // Add the Sub-step ??
881  << std::setw(12) << "Step-len" << " "
882  << std::setw(12) << "Step-len" << " "
883  << std::setw( 9) << "ReqStep" << " "
884  << G4endl;
885  }
886 
887  if( (subStepNo <= 0) )
888  {
889  PrintStat_Aux( StartFT, requestStep, 0.,
890  0, 0.0, 1.0);
891  //*************
892  }
893 
894  if( verboseLevel <= 3 )
895  {
896  G4cout.precision(noPrecision);
897  PrintStat_Aux( CurrentFT, requestStep, step_len,
898  subStepNo, subStepSize, DotStartCurrentVeloc );
899  //*************
900  }
901 
902  else // if( verboseLevel > 3 )
903  {
904  // Multi-line output
905 
906  // G4cout << "Current Position is " << CurrentPosition << G4endl
907  // << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
908  // G4cout << "Step taken was " << step_len
909  // << " out of PhysicalStep= " << requestStep << G4endl;
910  // G4cout << "Final safety is: " << safety << G4endl;
911  // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
912  // << G4endl << G4endl;
913  }
914  G4cout.precision(oldPrec);
915 }
916 
917 // ---------------------------------------------------------------------------
918 
919 void G4MagInt_Driver::PrintStat_Aux(
920  const G4FieldTrack& aFieldTrack,
921  G4double requestStep,
922  G4double step_len,
923  G4int subStepNo,
924  G4double subStepSize,
925  G4double dotVeloc_StartCurr)
926 {
927  const G4ThreeVector Position= aFieldTrack.GetPosition();
928  const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
929 
930  if( subStepNo >= 0)
931  {
932  G4cout << std::setw( 5) << subStepNo << " ";
933  }
934  else
935  {
936  G4cout << std::setw( 5) << "Start" << " ";
937  }
938  G4double curveLen= aFieldTrack.GetCurveLength();
939  G4cout << std::setw( 7) << curveLen;
940  G4cout << std::setw( 9) << Position.x() << " "
941  << std::setw( 9) << Position.y() << " "
942  << std::setw( 9) << Position.z() << " "
943  << std::setw( 8) << UnitVelocity.x() << " "
944  << std::setw( 8) << UnitVelocity.y() << " "
945  << std::setw( 8) << UnitVelocity.z() << " ";
946  G4int oldprec= G4cout.precision(3);
947  G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
948  G4cout.precision(6);
949  G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
950  G4cout.precision(oldprec);
951  G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
952  G4cout << std::setw(12) << step_len << " ";
953 
954  static G4double oldCurveLength= 0.0;
955  static G4double oldSubStepLength= 0.0;
956  static G4int oldSubStepNo= -1;
957 
958  G4double subStep_len=0.0;
959  if( curveLen > oldCurveLength )
960  {
961  subStep_len= curveLen - oldCurveLength;
962  }
963  else if (subStepNo == oldSubStepNo)
964  {
965  subStep_len= oldSubStepLength;
966  }
967  oldCurveLength= curveLen;
968  oldSubStepLength= subStep_len;
969 
970  G4cout << std::setw(12) << subStep_len << " ";
971  G4cout << std::setw(12) << subStepSize << " ";
972  if( requestStep != -1.0 )
973  {
974  G4cout << std::setw( 9) << requestStep << " ";
975  }
976  else
977  {
978  G4cout << std::setw( 9) << " InitialStep " << " ";
979  }
980  G4cout << G4endl;
981 }
982 
983 // ---------------------------------------------------------------------------
984 
985 void G4MagInt_Driver::PrintStatisticsReport()
986 {
987  G4int noPrecBig= 6;
988  G4int oldPrec= G4cout.precision(noPrecBig);
989 
990  G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
991  G4cout << "G4MagInt_Driver: Number of Steps: "
992  << " Total= " << fNoTotalSteps
993  << " Bad= " << fNoBadSteps
994  << " Small= " << fNoSmallSteps
995  << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
996  << G4endl;
997 
998 #ifdef G4FLD_STATS
999  G4cout << "MID dyerr: "
1000  << " maximum= " << fDyerr_max
1001  << " Sum small= " << fDyerrPos_smTot
1002  << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
1003  << " vel= " << std::sqrt( fDyerrVel_lgTot )
1004  << " Total h-distance: small= " << fSumH_sm
1005  << " large= " << fSumH_lg
1006  << G4endl;
1007 
1008 #if 0
1009  G4int noPrecSmall=4;
1010  // Single line precis of statistics ... optional
1011  G4cout.precision(noPrecSmall);
1012  G4cout << "MIDnums: " << fMinimumStep
1013  << " " << fNoTotalSteps
1014  << " " << fNoSmallSteps
1015  << " " << fNoSmallSteps-fNoInitialSmallSteps
1016  << " " << fNoBadSteps
1017  << " " << fDyerr_max
1018  << " " << fDyerr_mx2
1019  << " " << fDyerrPos_smTot
1020  << " " << fSumH_sm
1021  << " " << fDyerrPos_lgTot
1022  << " " << fDyerrVel_lgTot
1023  << " " << fSumH_lg
1024  << G4endl;
1025 #endif
1026 #endif
1027 
1028  G4cout.precision(oldPrec);
1029 }
1030 
1031 // ---------------------------------------------------------------------------
1032 
1033 void G4MagInt_Driver::SetSmallestFraction(G4double newFraction)
1034 {
1035  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
1036  {
1037  fSmallestFraction= newFraction;
1038  }
1039  else
1040  {
1041  G4cerr << "Warning: SmallestFraction not changed. " << G4endl
1042  << " Proposed value was " << newFraction << G4endl
1043  << " Value must be between 1.e-8 and 1.e-16" << G4endl;
1044  }
1045 }
Generated on Sat May 25 2013 14:33:10 for Geant4 by   doxygen 1.8.4