Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4MagIntegratorDriver.cc
이 파일의 문서화 페이지로 가기
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 110753 2018-06-12 15:44:03Z 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 // ---------------------------------------------------------
51 
52 // Constructor
53 //
55  G4MagIntegratorStepper *pStepper,
56  G4int numComponents,
57  G4int statisticsVerbose)
58  : fSmallestFraction( 1.0e-12 ),
59  fNoIntegrationVariables(numComponents),
60  fMinNoVars(12),
61  fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
62  fStatisticsVerboseLevel(statisticsVerbose),
63  fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
64  fNoInitialSmallSteps(0), fNoCalls(0),
65  fDyerr_max(0.0), fDyerr_mx2(0.0),
66  fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
67  fSumH_sm(0.0), fSumH_lg(0.0),
68  fVerboseLevel(0)
69 {
70  // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
71  // is required. For proper time of flight and spin, fMinNoVars must be 12
72 
73  RenewStepperAndAdjust( pStepper );
74  fMinimumStep= hminimum;
75 
76  // The (default) maximum number of steps is Base
77  // divided by the order of Stepper
78  //
79  fMaxStepBase = 250; // Was 5000
80 
82 #ifdef G4DEBUG_FIELD
83  fVerboseLevel=2;
84 #endif
85 
86  if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
87  {
88  G4cout << "MagIntDriver version: Accur-Adv: "
89  << "invE_nS, QuickAdv-2sqrt with Statistics "
90 #ifdef G4FLD_STATS
91  << " enabled "
92 #else
93  << " disabled "
94 #endif
95  << G4endl;
96  }
97 }
98 
99 // ---------------------------------------------------------
100 
101 // Destructor
102 //
104 {
105  if( fStatisticsVerboseLevel > 1 )
106  {
108  }
109 }
110 
111 // ---------------------------------------------------------
112 
113 G4bool
115  G4double hstep,
116  G4double eps,
117  G4double hinitial )
118 {
119  // Runge-Kutta driver with adaptive stepsize control. Integrate starting
120  // values at y_current over hstep x2 with accuracy eps.
121  // On output ystart is replaced by values at the end of the integration
122  // interval. RightHandSide is the right-hand side of ODE system.
123  // The source is similar to odeint routine from NRC p.721-722 .
124 
125  G4int nstp, i, no_warnings=0;
126  G4double x, hnext, hdid, h;
127 
128 #ifdef G4DEBUG_FIELD
129  static G4int dbg=1;
130  static G4int nStpPr=50; // For debug printing of long integrations
131  G4double ySubStepStart[G4FieldTrack::ncompSVEC];
132  G4FieldTrack yFldTrkStart(y_current);
133 #endif
134 
137  G4double x1, x2;
138  G4bool succeeded = true, lastStepSucceeded;
139 
140  G4double startCurveLength;
141 
142  G4int noFullIntegr=0, noSmallIntegr = 0 ;
143  static G4ThreadLocal G4int noGoodSteps =0 ; // Bad = chord > curve-len
144  const G4int nvar= fNoVars;
145 
146  G4FieldTrack yStartFT(y_current);
147 
148  // Ensure that hstep > 0
149  //
150  if( hstep <= 0.0 )
151  {
152  if(hstep==0.0)
153  {
154  std::ostringstream message;
155  message << "Proposed step is zero; hstep = " << hstep << " !";
156  G4Exception("G4MagInt_Driver::AccurateAdvance()",
157  "GeomField1001", JustWarning, message);
158  return succeeded;
159  }
160  else
161  {
162  std::ostringstream message;
163  message << "Invalid run condition." << G4endl
164  << "Proposed step is negative; hstep = " << hstep << "." << G4endl
165  << "Requested step cannot be negative! Aborting event.";
166  G4Exception("G4MagInt_Driver::AccurateAdvance()",
167  "GeomField0003", EventMustBeAborted, message);
168  return false;
169  }
170  }
171 
172  y_current.DumpToArray( ystart );
173 
174  startCurveLength= y_current.GetCurveLength();
175  x1= startCurveLength;
176  x2= x1 + hstep;
177 
178  if ( (hinitial > 0.0) && (hinitial < hstep)
179  && (hinitial > perMillion * hstep) )
180  {
181  h = hinitial;
182  }
183  else // Initial Step size "h" defaults to the full interval
184  {
185  h = hstep;
186  }
187 
188  x = x1;
189 
190  for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
191 
192  G4bool lastStep= false;
193  nstp=1;
194 
195  do
196  {
197  G4ThreeVector StartPos( y[0], y[1], y[2] );
198 
199 #ifdef G4DEBUG_FIELD
200  G4double xSubStepStart= x;
201  for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
202  yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
203  yFldTrkStart.SetCurveLength(x);
204 #endif
205 
206  // Old method - inline call to Equation of Motion
207  // pIntStepper->RightHandSide( y, dydx );
208  // New method allows to cache field, or state (eg momentum magnitude)
209  pIntStepper->ComputeRightHandSide( y, dydx );
210  fNoTotalSteps++;
211 
212  // Perform the Integration
213  //
214  if( h > fMinimumStep )
215  {
216  OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
217  //--------------------------------------
218  lastStepSucceeded= (hdid == h);
219 #ifdef G4DEBUG_FIELD
220  if (dbg>2)
221  {
222  PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
223  }
224 #endif
225  }
226  else
227  {
228  G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
229  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
230  G4double dchord_step, dyerr, dyerr_len; // What to do with these ?
232  yFldTrk.SetCurveLength( x );
233 
234  QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
235  //-----------------------------------------------------
236 
237  yFldTrk.DumpToArray(y);
238 
239 #ifdef G4FLD_STATS
240  fNoSmallSteps++;
241  if ( dyerr_len > fDyerr_max) { fDyerr_max= dyerr_len; }
242  fDyerrPos_smTot += dyerr_len;
243  fSumH_sm += h; // Length total for 'small' steps
244  if (nstp<=1) { fNoInitialSmallSteps++; }
245 #endif
246 #ifdef G4DEBUG_FIELD
247  if (dbg>1)
248  {
249  if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
250  G4cout << "Another sub-min step, no " << fNoSmallSteps
251  << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
252  PrintStatus( ySubStepStart, x1, y, x, h, nstp); // Only this
253  G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
254  << " epsilon= " << eps << " hstep= " << hstep
255  << " h= " << h << " hmin= " << fMinimumStep << G4endl;
256  }
257 #endif
258  if( h == 0.0 )
259  {
260  G4Exception("G4MagInt_Driver::AccurateAdvance()",
261  "GeomField0003", FatalException,
262  "Integration Step became Zero!");
263  }
264  dyerr = dyerr_len / h;
265  hdid= h;
266  x += hdid;
267 
268  // Compute suggested new step
269  hnext= ComputeNewStepSize( dyerr/eps, h);
270 
271  // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
272  lastStepSucceeded= (dyerr<= eps);
273  }
274 
275  if (lastStepSucceeded) { noFullIntegr++; }
276  else { noSmallIntegr++; }
277 
278  G4ThreeVector EndPos( y[0], y[1], y[2] );
279 
280 #ifdef G4DEBUG_FIELD
281  if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
282  {
283  if( nstp==nStpPr ) { G4cout << "***** Many steps ****" << G4endl; }
284  G4cout << "MagIntDrv: " ;
285  G4cout << "hdid=" << std::setw(12) << hdid << " "
286  << "hnext=" << std::setw(12) << hnext << " "
287  << "hstep=" << std::setw(12) << hstep << " (requested) "
288  << G4endl;
289  PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
290  }
291 #endif
292 
293  // Check the endpoint
294  G4double endPointDist= (EndPos-StartPos).mag();
295  if ( endPointDist >= hdid*(1.+perMillion) )
296  {
297  fNoBadSteps++;
298 
299  // Issue a warning only for gross differences -
300  // we understand how small difference occur.
301  if ( endPointDist >= hdid*(1.+perThousand) )
302  {
303 #ifdef G4DEBUG_FIELD
304  if (dbg)
305  {
306  WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
307  G4cerr << " Total steps: bad " << fNoBadSteps
308  << " good " << noGoodSteps << " current h= " << hdid
309  << G4endl;
310  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
311  }
312 #endif
313  no_warnings++;
314  }
315  }
316  else
317  {
318  noGoodSteps ++;
319  }
320 // #endif
321 
322  // Avoid numerous small last steps
323  if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
324  {
325  // No more integration -- the next step will not happen
326  lastStep = true;
327  }
328  else
329  {
330  // Check the proposed next stepsize
331  if(std::fabs(hnext) <= Hmin())
332  {
333 #ifdef G4DEBUG_FIELD
334  // If simply a very small interval is being integrated, do not warn
335  if( (x < x2 * (1-eps) ) && // The last step can be small: OK
336  (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
337  {
338  if(dbg>0)
339  {
340  WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
341  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
342  }
343  no_warnings++;
344  }
345 #endif
346  // Make sure that the next step is at least Hmin.
347  h = Hmin();
348  }
349  else
350  {
351  h = hnext;
352  }
353 
354  // Ensure that the next step does not overshoot
355  if ( x+h > x2 )
356  { // When stepsize overshoots, decrease it!
357  h = x2 - x ; // Must cope with difficult rounding-error
358  } // issues if hstep << x2
359 
360  if ( h == 0.0 )
361  {
362  // Cannot progress - accept this as last step - by default
363  lastStep = true;
364 #ifdef G4DEBUG_FIELD
365  if (dbg>2)
366  {
367  int prec= G4cout.precision(12);
368  G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
369  << G4endl
370  << " Integration step 'h' became "
371  << h << " due to roundoff. " << G4endl
372  << " Calculated as difference of x2= "<< x2 << " and x=" << x
373  << " Forcing termination of advance." << G4endl;
374  G4cout.precision(prec);
375  }
376 #endif
377  }
378  }
379  } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
380  // Loop checking, 07.10.2016, J. Apostolakis
381 
382  // Have we reached the end ?
383  // --> a better test might be x-x2 > an_epsilon
384 
385  succeeded= (x>=x2); // If it was a "forced" last step
386 
387  for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
388 
389  // Put back the values.
390  y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
391  y_current.SetCurveLength( x );
392 
393  if(nstp > fMaxNoSteps)
394  {
395  no_warnings++;
396  succeeded = false;
397 #ifdef G4DEBUG_FIELD
398  if (dbg)
399  {
400  WarnTooManyStep( x1, x2, x ); // Issue WARNING
401  PrintStatus( yEnd, x1, y, x, hstep, -nstp);
402  }
403 #endif
404  }
405 
406 #ifdef G4DEBUG_FIELD
407  if( dbg && no_warnings )
408  {
409  G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp << G4endl;
410  PrintStatus( yEnd, x1, y, x, hstep, nstp);
411  }
412 #endif
413 
414  return succeeded;
415 } // end of AccurateAdvance ...........................
416 
417 // ---------------------------------------------------------
418 
419 void
421  G4double h, G4double xDone,
422  G4int nstp)
423 {
424  static G4ThreadLocal G4int noWarningsIssued =0;
425  const G4int maxNoWarnings = 10; // Number of verbose warnings
426  std::ostringstream message;
427  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
428  {
429  message << "The stepsize for the next iteration, " << hnext
430  << ", is too small - in Step number " << nstp << "." << G4endl
431  << "The minimum for the driver is " << Hmin() << G4endl
432  << "Requested integr. length was " << hstep << " ." << G4endl
433  << "The size of this sub-step was " << h << " ." << G4endl
434  << "The integrations has already gone " << xDone;
435  }
436  else
437  {
438  message << "Too small 'next' step " << hnext
439  << ", step-no: " << nstp << G4endl
440  << ", this sub-step: " << h
441  << ", req_tot_len: " << hstep
442  << ", done: " << xDone << ", min: " << Hmin();
443  }
444  G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
445  JustWarning, message);
446  noWarningsIssued++;
447 }
448 
449 // ---------------------------------------------------------
450 
451 void
453  G4double x2end,
454  G4double xCurrent)
455 {
456  std::ostringstream message;
457  message << "The number of steps used in the Integration driver"
458  << " (Runge-Kutta) is too many." << G4endl
459  << "Integration of the interval was not completed !" << G4endl
460  << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
461  << " % fraction of it was done.";
462  G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
463  JustWarning, message);
464 }
465 
466 // ---------------------------------------------------------
467 
468 void
470  G4double h ,
471  G4double eps,
472  G4int dbg)
473 {
474  static G4ThreadLocal G4double maxRelError=0.0;
475  G4bool isNewMax, prNewMax;
476 
477  isNewMax = endPointDist > (1.0 + maxRelError) * h;
478  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
479  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
480 
482  && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
483  {
484  static G4ThreadLocal G4int noWarnings = 0;
485  std::ostringstream message;
486  if( (noWarnings ++ < 10) || (dbg>2) )
487  {
488  message << "The integration produced an end-point which " << G4endl
489  << "is further from the start-point than the curve length."
490  << G4endl;
491  }
492  message << " Distance of endpoints = " << endPointDist
493  << ", curve length = " << h << G4endl
494  << " Difference (curveLen-endpDist)= " << (h - endPointDist)
495  << ", relative = " << (h-endPointDist) / h
496  << ", epsilon = " << eps;
497  G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
498  JustWarning, message);
499  }
500 }
501 
502 // ---------------------------------------------------------
503 
504 void
506  const G4double dydx[],
507  G4double& x, // InOut
508  G4double htry,
509  G4double eps_rel_max,
510  G4double& hdid, // Out
511  G4double& hnext ) // Out
512 
513 // Driver for one Runge-Kutta Step with monitoring of local truncation error
514 // to ensure accuracy and adjust stepsize. Input are dependent variable
515 // array y[0,...,5] and its derivative dydx[0,...,5] at the
516 // starting value of the independent variable x . Also input are stepsize
517 // to be attempted htry, and the required accuracy eps. On output y and x
518 // are replaced by their new values, hdid is the stepsize that was actually
519 // accomplished, and hnext is the estimated next stepsize.
520 // This is similar to the function rkqs from the book:
521 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
522 // Edition, by William H. Press, Saul A. Teukolsky, William T.
523 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
524 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
525 
526 {
527  G4double errmax_sq;
528  G4double h, htemp, xnew ;
529 
531 
532  h = htry ; // Set stepsize to the initial trial value
533 
534  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
535 
536  G4double errpos_sq=0.0; // square of displacement error
537  G4double errvel_sq=0.0; // square of momentum vector difference
538  G4double errspin_sq=0.0; // square of spin vector difference
539 
540  G4int iter;
541 
542  static G4ThreadLocal G4int tot_no_trials=0;
543  const G4int max_trials=100;
544 
545  G4ThreeVector Spin(y[9],y[10],y[11]);
546  G4double spin_mag2 =Spin.mag2() ;
547  G4bool hasSpin= (spin_mag2 > 0.0);
548 
549  for (iter=0; iter<max_trials ;iter++)
550  {
551  tot_no_trials++;
552  pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
553  // *******
554  G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
555  G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
556 
557  // Evaluate accuracy
558  //
559  errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
560  errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
561 
562  // Accuracy for momentum
563  G4double magvel_sq= sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
564  G4double sumerr_sq = sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ;
565  if( magvel_sq > 0.0 )
566  {
567  errvel_sq = sumerr_sq / magvel_sq;
568  }
569  else
570  {
571  std::ostringstream message;
572  message << "Found case of zero momentum." << G4endl
573  << "- iteration= " << iter << "; h= " << h;
574  G4Exception("G4MagInt_Driver::OneGoodStep()",
575  "GeomField1001", JustWarning, message);
576  errvel_sq = sumerr_sq;
577  }
578  errvel_sq *= inv_eps_vel_sq;
579  errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
580 
581  if( hasSpin )
582  {
583  // Accuracy for spin
584  errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
585  / spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
586  errspin_sq *= inv_eps_vel_sq;
587  errmax_sq = std::max( errmax_sq, errspin_sq );
588  }
589 
590  if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
591 
592  // Step failed; compute the size of retrial Step.
593  htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
594 
595  if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
596  else { h = 0.1*h; } // reduce stepsize, but no more
597  // than a factor of 10
598  xnew = x + h;
599  if(xnew == x)
600  {
601  std::ostringstream message;
602  message << "Stepsize underflow in Stepper !" << G4endl
603  << "- Step's start x=" << x << " and end x= " << xnew
604  << " are equal !! " << G4endl
605  << " Due to step-size= " << h
606  << ". Note that input step was " << htry;
607  G4Exception("G4MagInt_Driver::OneGoodStep()",
608  "GeomField1001", JustWarning, message);
609  break;
610  }
611  }
612 
613  // Compute size of next Step
614  if (errmax_sq > errcon*errcon)
615  {
616  hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
617  }
618  else
619  {
620  hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
621  }
622  x += (hdid = h);
623 
624  for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
625 
626  return;
627 }
628 
629 //----------------------------------------------------------------------
630 
631 // QuickAdvance just tries one Step - it does not ensure accuracy
632 //
634  G4FieldTrack& y_posvel, // INOUT
635  const G4double dydx[],
636  G4double hstep, // In
637  G4double& dchord_step,
638  G4double& dyerr_pos_sq,
639  G4double& dyerr_mom_rel_sq )
640 {
641  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
642  FatalException, "Not yet implemented.");
643 
644  // Use the parameters of this method, to please compiler
645  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
646  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
647  return true;
648 }
649 
650 //----------------------------------------------------------------------
651 
653  G4FieldTrack& y_posvel, // INOUT
654  const G4double dydx[],
655  G4double hstep, // In
656  G4double& dchord_step,
657  G4double& dyerr )
658 {
659  G4double dyerr_pos_sq, dyerr_mom_rel_sq;
662  G4double s_start;
663  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
664 
665  static G4ThreadLocal G4int no_call=0;
666  no_call ++;
667 
668  // Move data into array
669  y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
670  s_start = y_posvel.GetCurveLength();
671 
672  // Do an Integration Step
673  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
674 
675  // Estimate curve-chord distance
676  dchord_step= pIntStepper-> DistChord();
677 
678  // Put back the values. yarrout ==> y_posvel
679  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
680  y_posvel.SetCurveLength( s_start + hstep );
681 
682 #ifdef G4DEBUG_FIELD
683  if(fVerboseLevel>2)
684  {
685  G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
686  PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
687  }
688 #endif
689 
690  // A single measure of the error
691  // TO-DO : account for energy, spin, ... ?
692  vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
693  inv_vel_mag_sq = 1.0 / vel_mag_sq;
694  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
695  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
696  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
697 
698  // Calculate also the change in the momentum squared also ???
699  // G4double veloc_square = y_posvel.GetVelocity().mag2();
700  // ...
701 
702 #ifdef RETURN_A_NEW_STEP_LENGTH
703  // The following step cannot be done here because "eps" is not known.
704  dyerr_len = std::sqrt( dyerr_len_sq );
705  dyerr_len_sq /= eps ;
706 
707  // Look at the velocity deviation ?
708  // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
709 
710  // Set suggested new step
711  hstep= ComputeNewStepSize( dyerr_len, hstep);
712 #endif
713 
714  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
715  {
716  dyerr = std::sqrt(dyerr_pos_sq);
717  }
718  else
719  {
720  // Scale it to the current step size - for now
721  dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
722  }
723 
724  return true;
725 }
726 
727 // --------------------------------------------------------------------------
728 
729 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
731  G4double yarrin[], // In
732  const G4double dydx[],
733  G4double hstep, // In
734  G4double yarrout[],
735  G4double& dchord_step,
736  G4double& dyerr ) // In length
737 {
738  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
739  FatalException, "Not yet implemented.");
740  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
741  yarrout[0]= yarrin[0];
742 }
743 #endif
744 
745 // --------------------------------------------------------------------------
746 
747 // This method computes new step sizes - but does not limit changes to
748 // within certain factors
749 //
751 ComputeNewStepSize(G4double errMaxNorm, // max error (normalised)
752  G4double hstepCurrent) // current step size
753 {
755 
756  // Compute size of next Step for a failed step
757  if(errMaxNorm > 1.0 )
758  {
759  // Step failed; compute the size of retrial Step.
760  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
761  }
762  else if(errMaxNorm > 0.0 )
763  {
764  // Compute size of next Step for a successful step
765  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
766  }
767  else
768  {
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
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 
842  const G4FieldTrack& StartFT,
843  const G4FieldTrack& CurrentFT,
844  G4double requestStep,
845  G4int subStepNo)
846 {
847  G4int verboseLevel= fVerboseLevel;
848  const 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  if( verboseLevel <= 3 )
894  {
895  G4cout.precision(noPrecision);
896  PrintStat_Aux( CurrentFT, requestStep, step_len,
897  subStepNo, subStepSize, DotStartCurrentVeloc );
898  }
899 
900  G4cout.precision(oldPrec);
901 }
902 
903 // ---------------------------------------------------------------------------
904 
906  const G4FieldTrack& aFieldTrack,
907  G4double requestStep,
908  G4double step_len,
909  G4int subStepNo,
910  G4double subStepSize,
911  G4double dotVeloc_StartCurr)
912 {
913  const G4ThreeVector Position= aFieldTrack.GetPosition();
914  const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
915 
916  if( subStepNo >= 0)
917  {
918  G4cout << std::setw( 5) << subStepNo << " ";
919  }
920  else
921  {
922  G4cout << std::setw( 5) << "Start" << " ";
923  }
924  G4double curveLen= aFieldTrack.GetCurveLength();
925  G4cout << std::setw( 7) << curveLen;
926  G4cout << std::setw( 9) << Position.x() << " "
927  << std::setw( 9) << Position.y() << " "
928  << std::setw( 9) << Position.z() << " "
929  << std::setw( 8) << UnitVelocity.x() << " "
930  << std::setw( 8) << UnitVelocity.y() << " "
931  << std::setw( 8) << UnitVelocity.z() << " ";
932  G4int oldprec= G4cout.precision(3);
933  G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
934  G4cout.precision(6);
935  G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
936  G4cout.precision(oldprec);
937  G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
938  G4cout << std::setw(12) << step_len << " ";
939 
940  static G4ThreadLocal G4double oldCurveLength= 0.0;
941  static G4ThreadLocal G4double oldSubStepLength= 0.0;
942  static G4ThreadLocal G4int oldSubStepNo= -1;
943 
944  G4double subStep_len=0.0;
945  if( curveLen > oldCurveLength )
946  {
947  subStep_len= curveLen - oldCurveLength;
948  }
949  else if (subStepNo == oldSubStepNo)
950  {
951  subStep_len= oldSubStepLength;
952  }
953  oldCurveLength= curveLen;
954  oldSubStepLength= subStep_len;
955 
956  G4cout << std::setw(12) << subStep_len << " ";
957  G4cout << std::setw(12) << subStepSize << " ";
958  if( requestStep != -1.0 )
959  {
960  G4cout << std::setw( 9) << requestStep << " ";
961  }
962  else
963  {
964  G4cout << std::setw( 9) << " InitialStep " << " ";
965  }
966  G4cout << G4endl;
967 }
968 
969 // ---------------------------------------------------------------------------
970 
972 {
973  G4int noPrecBig= 6;
974  G4int oldPrec= G4cout.precision(noPrecBig);
975 
976  G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
977  G4cout << "G4MagInt_Driver: Number of Steps: "
978  << " Total= " << fNoTotalSteps
979  << " Bad= " << fNoBadSteps
980  << " Small= " << fNoSmallSteps
981  << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
982  << G4endl;
983  G4cout.precision(oldPrec);
984 }
985 
986 // ---------------------------------------------------------------------------
987 
989 {
990  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
991  {
992  fSmallestFraction= newFraction;
993  }
994  else
995  {
996  std::ostringstream message;
997  message << "Smallest Fraction not changed. " << G4endl
998  << " Proposed value was " << newFraction << G4endl
999  << " Value must be between 1.e-8 and 1.e-16";
1000  G4Exception("G4MagInt_Driver::SetSmallestFraction()",
1001  "GeomField1001", JustWarning, message);
1002  }
1003 }
1004 
1005 void G4MagInt_Driver::
1006 GetDerivatives(const G4FieldTrack& y_curr, G4double* dydx) const
1007 {
1009  y_curr.DumpToArray(ytemp);
1010  GetStepper()->RightHandSide(ytemp, dydx);
1011 }
1012 
1014 {
1015  return pIntStepper->GetEquationOfMotion();
1016 }
1017 
1019 {
1020  pIntStepper->SetEquationOfMotion(equation);
1021 }
1022 
1024 {
1025  return pIntStepper;
1026 }
1027 
1029 {
1030  return pIntStepper;
1031 }
1032 
1033 void G4MagInt_Driver::
1035 {
1036  pIntStepper = pItsStepper;
1037  ReSetParameters();
1038 }
Float_t x
Definition: compare.C:6
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0) override
static const double prec
Definition: RanecuEngine.cc:58
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
CLHEP::Hep3Vector G4ThreeVector
virtual G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr) override
G4EquationOfMotion * GetEquationOfMotion()
void SetEquationOfMotion(G4EquationOfMotion *newEquation)
Float_t x1[n_points_granero]
Definition: compare.C:5
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
void DumpToArray(G4double valArr[ncompSVEC]) const
double z() const
G4double GetPshrnk() const
static constexpr double perMillion
Definition: G4SIunits.hh:334
double dot(const Hep3Vector &) const
G4double G4MagInt_Driver::G4double hstepCurrent G4double hnew
void SetSmallestFraction(G4double val)
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent) override
#define G4ThreadLocal
Definition: tls.hh:69
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
const G4int fNoIntegrationVariables
G4MagInt_Driver(G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
unsigned long fNoSmallSteps
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
G4ThreeVector GetPosition() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double GetSafety() const
virtual void GetDerivatives(const G4FieldTrack &y_curr, G4double dydx[]) const override
G4double GetKineticEnergy() const
G4MagIntegratorStepper * pIntStepper
void SetCurveLength(G4double nCurve_s)
void RightHandSide(const double y[], double dydx[]) const
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
double mag2() const
virtual void SetEquationOfMotion(G4EquationOfMotion *equation) override
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
G4GLOB_DLL std::ostream G4cerr
virtual G4EquationOfMotion * GetEquationOfMotion() override
G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, G4double hstepCurrent)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
G4double GetPgrow() const
virtual ~G4MagInt_Driver() override
const G4ThreeVector & GetMomentumDir() const
void ReSetParameters(G4double new_safety=0.9)
G4GLOB_DLL std::ostream G4cout
double x() const
static G4GeometryTolerance * GetInstance()
T sqr(const T &x)
Definition: templates.hh:145
double y() const
virtual G4int IntegratorOrder() const =0
static const G4double eps
unsigned long fNoInitialSmallSteps
unsigned long fNoTotalSteps
Float_t x2[n_points_geant4]
Definition: compare.C:26
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
G4double Hmin() const
virtual void RenewStepperAndAdjust(G4MagIntegratorStepper *pItsStepper) override
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
virtual const G4MagIntegratorStepper * GetStepper() const override
static constexpr double perThousand
Definition: G4SIunits.hh:333
G4double GetSurfaceTolerance() const
G4double GetCurveLength() const