Rolling & Scalar IntegrationDriver: count only 'good' steps towards step number limit

RollingIntegrationDriver: count only 'good' steps towards step number limit
 - Co-works with changes in ScalarIntegrationDriver.

ScalarIntegrationDriver: changed while exit condition re number of iterations
- Co-works with change in RollingIntegrationDriver that counts 'good' steps only
- Suppressed old code for very small steps. [ Cleaned up if( true ) & else .]
   That 'branch' was not used in order to be compatible with vector code.

Agreement shown in test case 3 of testVectorIntegrationDriver

Small changes in aux classes: BaseRkIntegrationDriver, TrialsStats, FormattedReporter
 - BaseRkIntegrationDriver:  number of components: int -> unsigned & check it
 - FormattedReporter: small refinement in ReportOneLane - use separate ostream, then send to cout
 - TrialStats: Added method PrintSummary.  Fixed names in some printouts.

magneticfield/inc/Geant/BaseRkIntegrationDriver.h

@@ -54,11 +54,11 @@ public:
   template <typename T>
   using Vector3D = vecgeom::Vector3D<T>;
 
-  BaseRkIntegrationDriver(double     hminimum, // same
-                          T_Stepper *pStepper,
-                          double     maxRelativeEpsilon,
-                          int        numberOfComponents = 6,
-                          int        statsVerbosity = 1);
+  BaseRkIntegrationDriver(double        hminimum, // same
+                          T_Stepper   * pStepper,
+                          double        maxRelativeEpsilon,
+                          unsigned int  numberOfComponents = 6,
+                          int           statsVerbosity = 1);
 
   virtual ~BaseRkIntegrationDriver();
 
@@ -213,7 +213,7 @@ template <class T_Stepper, unsigned int Nvar>
 BaseRkIntegrationDriver<T_Stepper, Nvar>::BaseRkIntegrationDriver(double     hminimum,
                                                                   T_Stepper *pStepper,
                                                                   double     maxRelativeEpsilon,
-                                                                  int        numComponents,
+                                                                  unsigned int numComponents,
                                                                   int        statisticsVerbose)
     :
       FlexIntegrationDriver( maxRelativeEpsilon ),
@@ -233,7 +233,12 @@ BaseRkIntegrationDriver<T_Stepper, Nvar>::BaseRkIntegrationDriver(double     hmi
   // is required. For proper time of flight and spin,  fMinNoVars must be 12
   assert(pStepper != nullptr);
   assert(Nvar <= (unsigned int)numComponents); // Ensure that arrays are large enough for Integr.
-
+  if( Nvar > numComponents ) {
+     std::cerr << " BaseRkIntegrationDriver c-tor:  Incompatibilitye between Nvar= " << Nvar
+               << "  and the number of components " << numComponents << std::endl;
+     exit(1);
+  }
+     
   // fpStepper = pStepper;
 
   SetMaxNoSteps( fMaxStepBase / pStepper->GetIntegratorOrder() );

magneticfield/inc/Geant/FormattedReporter.h

@@ -2,6 +2,10 @@
 #include <cstdlib>
 #include <iomanip>
 
+#include <iostream>
+#include <string>
+#include <sstream>
+
 // #include "Geant/IntegrationDriverConstants.h"  // For ReportOneLane
 
 //  Auxiliary methods - should be encapsulated into a separate helper class
@@ -252,19 +256,20 @@ void
    bool  laneIsDone = vecCore::Get( lanesDone , lanePrint );
    int   prec = 10; // precision
    int   wd = prec + 5;
-   int   oldPrec= cout.precision(prec);
+   // int   oldPrec= cout.precision(prec);
    bool  printSquares = false; // Old version - true
    bool  printValues  = true; 
 
+   std::stringstream ostr;
    // const int trackToPrint = IntegrationDriverConstants::GetInstance()->GetTrackToCheck();   
    
-   cout << std::setw(12) << methodName << " - ReportOneLane : "
+   ostr << std::setw(12) << methodName << " - ReportOneLane : "
         << " trk# "     << setw(3) << trackNum  /* trackToPrint */ << " "
         << " lane: "    << setw(3) << lanePrint << " > "
         << " iter = "   << setw(3) << iter << " #call= " << setw(5) << noCall;
    prec=6;
    wd = prec + 5;
-   cout << std::setprecision( prec )
+   ostr << std::setprecision( prec )
         << " h = "          << setw( wd ) << vecCore::Get( hStep ,       lanePrint )
         << " xStepStart = " << setw( wd ) << vecCore::Get( xStepStart ,  lanePrint )      
         << " Eps/x = "      << setw( wd ) << vecCore::Get( epsPosition , lanePrint );
@@ -276,31 +281,32 @@ void
    int   prec2 = 9; // precision
    int   wd2 = prec2 + 5;
 
-   cout.precision( prec2 );
+   ostr.precision( prec2 );
             
    if( printSquares ) 
-      cout
+      ostr
              << " errSq: x,p = "  << setw( wd2) << errPosLane2  // vecCore::Get( errPosSq ,    lanePrint )
              << " "               << setw( wd2) << errMomLane2  // vecCore::Get( errMomSq ,    lanePrint ) // << " "
              << " errMax^2 = "    << setw( wd2) << errMax2;     // vecCore::Get( errmax_sq ,   lanePrint );
    if( printValues ) 
-      cout
+      ostr
              << " error-x/p = "   << setw( wd2) << sqrt( errPosLane2 ) // Get( errPosSq ,    lanePrint ) )
              << " "               << setw( wd2) << sqrt( errMomLane2 ) // Get( errMomSq ,    lanePrint ) ) // << " "
              << " errMax = "      << setw( wd2) << sqrt( errMax2 )  ;  // Get( errmax_sq ,   lanePrint ) );
 
-   cout << " lane done = "   << laneIsDone;
+   ostr << " lane done = "   << laneIsDone;
 
    if( allDone >= 0 )
-      cout << " allDone = "     << allDone;
+      ostr << " allDone = "     << allDone;
              
-   cout << std::endl;
+   ostr << std::endl;
 
-   if( laneIsDone  ) cout << std::endl;
+   if( laneIsDone  ) ostr << std::endl;
 
-   cout.precision(oldPrec);
-   
-   // if( allDone ) cout << std::endl;
+   // ostr.precision(oldPrec);
+
+   cout << ostr.str(); 
+   // if( allDone ) ostr << std::endl;
    // **> Old mode - that printed all updates - not just ones in which this lane was active.
 }
 

magneticfield/inc/Geant/OldIntegrationDriver.h

@@ -1324,13 +1324,13 @@ OldIntegrationDriver< /*Real_v,*/ T_Stepper, Nvar>
 //
 {
    if( partDebug) 
-      std::cout<<"----Storage position (out-arr): "<< indOut << std::endl;
-
-   // int     indOut = indexArr[currIndex]; // might be sent directly to StoreOutput as well
-   // (void)nTracks;
+    std::cout << "----Storage position (out-arr): " << indOut
+              << " (ntracks= " << nTracks << ")" << std::endl;
+   (void) nTracks; // To enable its use in asserts above without warnings (in non-debug builds)
+   
    assert( 0 <= indOut && indOut < nTracks && "Track Index is Out of Range"  );
    assert( 0 <= currIndex && ((unsigned long)currIndex < vecCore::VectorSize<Real_v>() ) && "Lane Index is Out of Range" ) ;
-
+   
    double hOriginal = hstep [indOut];
    
    if (hOriginal >= 0.0)     

magneticfield/inc/Geant/SimpleIntegrationDriver.h

@@ -42,13 +42,13 @@
 #include "Geant/VectorTypes.h" //  Defines geant::Double_v
 #include "Geant/math_wrappers.h"
 
-#define CHECK_ONE_LANE 1
+// #define CHECK_ONE_LANE 1
 //  Define to check a single lane
 
 #define CONST_DEBUG    1
 //  Define to turn 'partDebug' into compile time constant
 
-#define DRIVER_PRINT_PROGRESS   1
+// #define DRIVER_PRINT_PROGRESS   1
 
 #ifdef CHECK_ONE_LANE
 #include "IntegrationDriverConstants.h"
@@ -281,8 +281,8 @@ private:
   
   // ---------------------------------------------------------------
   // Compilation constants
-#ifdef CONST_DEBUG  
-  const bool partDebug  = true;  // false;                 // Enforce debugging output
+#ifdef  CONST_DEBUG  
+  const bool partDebug  = false;                 // Enforce debugging output
 #endif
   const int ncompSVEC   = FieldTrack::NumCompFT; // expect 6, later 8, eventually up to 12
   const bool useOneStep = true;                  //  Algorithm selection - false for KeepStepping
@@ -302,7 +302,7 @@ private:
   // const int fMinNoVars; // Minimum number for TemplateFieldTrack<Real_v>
   // const int fNoVars;    // Full number of variable
 
-  static constexpr int fMaxStepBase = 250;
+  static constexpr int fMaxStepBase = 250;  // Was 250;   2019.05.14  for testing
 
   // static constexpr double fSafetyFactor= 0.9; // -> Failed to compile on clang 9.1 2017.12.05
   static constexpr double kSafetyFactor = 0.9;                                            //     OK ...
@@ -499,8 +499,7 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
   
   static  std::atomic<unsigned int> numCalls(0);
   int  currentCallNo= numCalls++;
-  
-  const Real_v xStart(x);
+
   Real_v xnew;
   Real_v yerr[Nvar], ytemp[Nvar];
 
@@ -517,8 +516,13 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
 #ifndef CONST_DEBUG
   partDebug = ( laneToCheck != -1 );
 #endif
+#ifdef DRIVER_PRINT_PROGRESS  
+  const Real_v xStart(x);
+#endif
+  
   // bool  printNow = (laneToCheck >= 0) && vecCore::Get( htry, laneToCheck ) > 0.0 ;       
-  if( laneToCheck != -1 ) {
+  // if ( laneToCheck != -1 ) {
+  if ( false ) {
      cout << "* SID:AccAdv Global Vars> trackToPrint = " << trackToPrint << " l2c / laneToCheck = " << laneToCheck;
      cout << "  Args>  h[l2c] = " << vecCore::Get( htry, laneToCheck );     
      cout << endl;
@@ -551,17 +555,21 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
 #ifdef CHECK_ONE_LANE    
     bool  printLane = (laneToCheck >= 0) && vecCore::Get( Active, laneToCheck );
     
-    bool anyProgressLastStep = ! vecCore::MaskEmpty(goodStep) ;    
-    cout << " - iter = " << iter << " anyProgressLastStep = " << anyProgressLastStep << std::endl;       
-    ReportManyRowsOfDoubles("Current X/P",  yStart, 6 );
-    ReportRowOfDoubles("Charge",  charge, 6 );       
-    ReportManyRowsOfDoubles("d/ds [X,P] ", dydx, 6 );    
+    bool anyProgressLastStep = ! vecCore::MaskEmpty(goodStep) ;
+    if( partDebug ) { 
+       cout << " - iter = " << iter << " anyProgressLastStep = " << anyProgressLastStep << std::endl;       
+       ReportManyRowsOfDoubles("Current X/P",  yStart, 6 );
+       ReportRowOfDoubles("Charge",  charge, 6 );       
+       ReportManyRowsOfDoubles("d/ds [X,P] ", dydx, 6 );
+    }
 #endif    
     itersLeft--;
     iter++;
+#ifdef DRIVER_PRINT_PROGRESS    
     if (partDebug) cout << " OneGoodStep - iteration = " << iter << "             "
                         << " ( total iterations = " <<  (tot_no_trials + iter) << " ) "
                         << endl;
+#endif
     // #ifdef STORE_ONCE
     vecCore::MaskedAssign(h, finished, Real_v(0.0)); // Set h = 0.0 for finished lanes -- ensure no change !
                                                      // #endif
@@ -630,20 +638,22 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
 #ifdef DRIVER_PRINT_PROGRESS        
     Real_v hAdvance= vecCore::Blend( goodStep, h, Real_v(0.0) );
 
-    // Real_v xNext = x + hStep; //   Updates both good and bad steps -- else it should be  = x + hAdvance;
     Real_v xNext = x + hAdvance; //   Updates only good steps --  2019.01.07
 
-    cout << " ====================================================================" << endl;
-    ReportRowOfDoubles("charge ***>", charge );    
-    ReportRowOfDoubles("x:   Start  =", x );
-    ReportRowOfDoubles("xNext: End  =", xNext );
-    ReportRowOfDoubles("hStep/tried =", h );    
-    ReportRowOfDoubles("hAdvance/got=", hAdvance );
-    // ReportRowOfDoubles("hdid        =", hdid );
-    ReportRowOfBools<Real_v>("goodStep", goodStep);
+    bool ReportAfterGoodStep = false;
+    if (partDebug && ReportAfterGoodStep) {    
+       cout << " ====================================================================" << endl;
+       ReportRowOfDoubles("charge ***>", charge );    
+       ReportRowOfDoubles("x:   Start  =", x );
+       ReportRowOfDoubles("xNext: End  =", xNext );
+       ReportRowOfDoubles("hStep/tried =", h );    
+       ReportRowOfDoubles("hAdvance/got=", hAdvance );
+       // ReportRowOfDoubles("hdid        =", hdid );
+       ReportRowOfBools<Real_v>("goodStep", goodStep);
     
-    ReportRowOfDoubles("x: Updated =", x );    
-    cout << " ====================================================================" << endl;
+       ReportRowOfDoubles("x: Updated =", x );    
+       cout << " ====================================================================" << endl;
+    }
 #endif 
     
     Bool_v laneDone = (goodStep | finished);
@@ -658,18 +668,18 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
                        epsPosition, errpos_sq, errmom_sq, errmax_sq, laneDone,
                        allDone, iter, tot_no_trials, laneToCheck, trackToPrint,
                        "SimpleID" );  // "SimpleIntDrv" );
+#ifdef DRIVER_PRINT_PROGRESS       
        std::cout << " Track to check " << trackToPrint << " laneToCheck = " << laneToCheck << std::endl;       
        ReportRowOfSquareRoots("ErrPos", errpos_sq );
        ReportRowOfSquareRoots("ErrMom", errmom_sq );
        ReportRowOfSquareRoots("ErrMax", errmax_sq );
        ReportManyRowsOfDoubles("yErr", yerr, Nvar);
        
-       if( 1 ) {
+       if( partDebug ) {
           std::cout << "SID: Status after stepper call ----------------------------------------------" << std::endl;
           FormattedReporter::FullReport(yStart, charge, dydx, h /*hStep*/, ytemp /*yStepEnd*/,
                                         yerr, errmax_sq, Active, goodStep );
-       }
-       if( 1 ) {
+
           // ReportRowOfSquareRoots("|err-p|", yerr[3]*yerr[3] + yerr[4]*yerr[4] + yerr[5]*yerr[5] );
           // ReportRowOfDoubles("up = SumErr^2", sumerr_sq );
           // ReportRowOfDoubles("dwn= magMom^2+e", magmom_sq + tinyValue );
@@ -677,6 +687,7 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
           // ReportRowOfDoubles("ErrMom^2", errmom_sq );
           // ReportRowOfSquareRoots("ErrMom", errmom_sq );
        }
+#endif
     }
     // End debug code                 -------------  2019.02.27
 #endif          
@@ -713,7 +724,7 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
     stepSizeUnderflow = Active && (xnew == x);
 
 #ifdef CHECK_ONE_LANE    
-    if( printLane ) {
+    if( printLane && partDebug ) {
       std::cout << "SID: After chance to break. (Not allDone.) Remaining lanes represent work: " << std::endl;
       ReportRowOfBools<Real_v>("laneDone", laneDone);
       ReportRowOfBools<Real_v>("Active", Active);
@@ -750,7 +761,7 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
     }
     // Idea 1.5 :  Use only one store for goodStep & underflow lanes (if continuing)
     if (!vecCore::MaskEmpty(stepSizeUnderflow || goodStep)) {
-#ifdef DRIVER_PRINT_PROGRESS       
+#ifdef DRIVER_PRINT_PROGRESS   
       if (partDebug) cout << "Store and Report Stored lanes - v1.5 allDone - good or Underflow." << endl;
 #endif
       //*************
@@ -819,11 +830,13 @@ void SimpleIntegrationDriver<T_Stepper, Nvar>::OneGoodStep(const Real_v   yStart
   Real_v errStretch =
     vecCore::Blend( underThresh, Real_v(fMaxSteppingIncrease), errStretch1raw );
 
-  ReportRowOfDoubles("errMaxSq(Final)=", errmax_sqFinal );
-  // ReportRowOfBools<Real_v>("underThresh", underThresh);
-  std::cout << " errcon = " << errcon << std::endl;
-  std::cout << " Power-Grow = " << 0.5 * GetPowerGrow() << std::endl;
-  if( errcon == 0.0 ) { std::cerr << " Simple IntDrv> ERROR: errcon is Zero.  Value = " << errcon << std::endl; } 
+#ifdef DRIVER_PRINT_PROGRESS
+  if( partDebug ) {
+     ReportRowOfDoubles("errMaxSq(Final)=", errmax_sqFinal );
+     // ReportRowOfBools<Real_v>("underThresh", underThresh);
+     // std::cout << " errcon = " << errcon << " Power-Grow = " << 0.5 * GetPowerGrow() << std::endl;
+  if( errcon == 0.0 ) { std::cerr << " Simple IntDrv> ERROR: errcon is Zero.  Value = " << errcon << std::endl; }
+#endif
   
   hnext = errStretch * h;
 
@@ -1097,12 +1110,16 @@ int SimpleIntegrationDriver</*Real_v,*/ T_Stepper, Nvar>::InitializeLanes(
       const int trackToPrint = IntegrationDriverConstants::GetInstance()->GetTrackToCheck();      
       // if( j == trackToPrint ) { laneToCheck = slot; }
       if( j == trackToPrint ) {
+         // int oldLane = laneToCheck;
+#ifdef DRIVER_PRINT_PROGRESS           
          // if( laneToCheck == -1 ) ...
          std::cout << "SID::InitLanes> found lane = " << slot
-                   << " (was = " << laneToCheck << " ) "
-                   << " for trackToPrint = " << trackToPrint << std::endl;
+                   << " (was = " << laneToCheck  // << oldLane
+                   
+                   << " )   for trackToPrint = " << trackToPrint << std::endl;
+#endif
          laneToCheck = slot;
-         // std::cout << "SID::InitLanes> found lane = " << laneToCheck << " for trackToPrint = " << trackToPrint << std::endl;         
+;         
       }
 #endif
       ++slot;
@@ -1243,7 +1260,6 @@ void SimpleIntegrationDriver< T_Stepper, Nvar>::
               << indOut
               // << " (ntracks= " << nTracks << ")"
               << std::endl;
-
   (void)nTracks; // Use the value in case on non-Debug builds - avoids compilation warning
 
   assert(0 <= indOut && indOut < nTracks && "Track Index is Out of Range");

magneticfield/inc/Geant/TrialStats.h

@@ -30,8 +30,8 @@ class TrialStats {
    unsigned int GetNoKeepSteps( unsigned int i ) const { return vecCore::Get( fNumKeepSteps, i ); } 
    
    //  Reset the count of selected lanes - after updating totals for lane
-   void fSumAndReset(vecCore::Mask_v<Real_v> resetLane );
-   void fSumAndResetLane(unsigned int i);
+   void SumAndReset(vecCore::Mask_v<Real_v> resetLane );
+   void SumAndResetLane(unsigned int i);
    // void fSumAndReset();
    
    //  Reset the count of selected lanes - no update of totals
@@ -48,15 +48,21 @@ class TrialStats {
    vecCore::Index_v<Real_v> GetSumBadSteps()  const { return fSumBadSteps;  }
    vecCore::Index_v<Real_v> GetSumKeepSteps() const { return fSumKeepSteps; } 
    
-   void PrintSums() const;
-   //  Prints only the Sum statistics
    void PrintStats( bool full = false ) const;
    //  Prints the current statistics - and optionally the Sums
+
+   void PrintSums() const;
+   //  Prints only the Sum statistics - per Lane
+   void PrintSummary() const;
+   //  Print Grand Totals
+   
    // void PrintAllStats() const;
    
    void PrintLaneStats( unsigned int lane ) const;
    
-  public:
+ private:
+   // STATE
+   // 1 - current counters
    vecCore::Index_v<Real_v> fNumTotalSteps= vecCore::Index_v<Real_v>(0);
    vecCore::Index_v<Real_v> fNumGoodSteps = vecCore::Index_v<Real_v>(0);   
    vecCore::Index_v<Real_v> fNumBadSteps  = vecCore::Index_v<Real_v>(0);
@@ -67,8 +73,9 @@ class TrialStats {
    // Extension idea: count 'initial' bad steps - ones at start of integration
    //  2019.05.15
    // vecCore::Index_v<Real_v> fNumInitialBadSteps = vecCore::Index_v<Real_v>(0);   
-   
+
    // Additional capability:  sums
+   // 2 - Sum counters
    vecCore::Index_v<Real_v> fSumTotalSteps= vecCore::Index_v<Real_v>(0);
    vecCore::Index_v<Real_v> fSumBadSteps  = vecCore::Index_v<Real_v>(0);
    vecCore::Index_v<Real_v> fSumGoodSteps = vecCore::Index_v<Real_v>(0);
@@ -89,11 +96,16 @@ void TrialStats<Real_v>::Update(vecCore::Mask_v<Real_v> active,
 template <typename Real_v>
 void TrialStats<Real_v>::ResetLane(unsigned int i)
 {
+   fSumTotalSteps +=  fNumTotalSteps;
+   fSumGoodSteps  +=  fNumGoodSteps;
+   fSumBadSteps   +=  fNumBadSteps;
+   fSumKeepSteps  +=  fNumKeepSteps;
    assert( i <= vecCore::VectorSize<Real_v>() );
+
    vecCore::Set( fNumTotalSteps, i, 0 );
    vecCore::Set( fNumBadSteps,   i, 0 );
    vecCore::Set( fNumGoodSteps,  i, 0 );
-   vecCore::Set( fNumKeepSteps,  i, 0 );   
+   vecCore::Set( fNumKeepSteps,  i, 0 );
 }
 
 template <typename Real_v>
@@ -106,7 +118,7 @@ void TrialStats<Real_v>::Reset(vecCore::Mask_v<Real_v> changeMask )
 }
 
 template <typename Real_v>
-void TrialStats<Real_v>::fSumAndReset(vecCore::Mask_v<Real_v> changeMask )
+void TrialStats<Real_v>::SumAndReset(vecCore::Mask_v<Real_v> changeMask )
 {
    vecCore::MaskedAssign( fSumTotalSteps,  changeMask, fSumTotalSteps + fNumTotalSteps );
    vecCore::MaskedAssign( fSumBadSteps,    changeMask, fSumBadSteps   + fNumBadSteps );
@@ -117,14 +129,14 @@ void TrialStats<Real_v>::fSumAndReset(vecCore::Mask_v<Real_v> changeMask )
 }
 
 template <typename Real_v>
-void TrialStats<Real_v>::fSumAndResetLane(unsigned int i)
+void TrialStats<Real_v>::SumAndResetLane(unsigned int i)
 {
    assert( i <= vecCore::VectorSize<Real_v>() );
    vecCore::Mask_v<Real_v> changeMask= vecCore::Mask_v<Real_v>(false);
 
    Set( changeMask, i, true );
    
-   fSumAndReset( changeMask );
+   SumAndReset( changeMask );
 }
 
 template <typename Real_v>
@@ -154,10 +166,27 @@ void TrialStats<Real_v>::FullReset()
 template <typename Real_v>
 void TrialStats<Real_v>::PrintSums() const
 {
-   FormattedReporter::ReportRowOfInts<Real_v>(  "fSum Total Steps",  fSumTotalSteps );
-   FormattedReporter::ReportRowOfInts<Real_v>(  "fSum Good  Steps",  fSumGoodSteps  );
-   FormattedReporter::ReportRowOfInts<Real_v>(  "fSum Bad   Steps",  fSumBadSteps  );
-   FormattedReporter::ReportRowOfInts<Real_v>(  "fSum Keep  Steps",  fSumKeepSteps );
+   FormattedReporter::ReportRowOfInts<Real_v>(  "Sum Total Steps",  fSumTotalSteps );
+   FormattedReporter::ReportRowOfInts<Real_v>(  "Sum Good  Steps",  fSumGoodSteps  );
+   FormattedReporter::ReportRowOfInts<Real_v>(  "Sum Bad   Steps",  fSumBadSteps  );
+   FormattedReporter::ReportRowOfInts<Real_v>(  "Sum Keep  Steps",  fSumKeepSteps );
+}
+
+template <typename Real_v>
+void TrialStats<Real_v>::PrintSummary() const
+{
+   constexpr unsigned int VecSize = vecCore::VectorSize<Real_v>();   
+   unsigned long  GrandTotalSteps(0UL), GrandTotalGood(0UL), GrandTotalBad(0UL), GrandTotalKeep(0UL);
+   for (unsigned int i = 0; i < VecSize; ++i) {
+      GrandTotalSteps += vecCore::Get( fSumTotalSteps, i );
+      GrandTotalGood  += vecCore::Get( fSumGoodSteps , i );
+      GrandTotalBad   += vecCore::Get( fSumBadSteps,   i );
+      GrandTotalKeep  += vecCore::Get( fSumKeepSteps,  i );
+   }
+   std::cout << "GrandSums:  Total step = " << fGrandSumTotalSteps
+      " Good  Steps= " <<  fGrandSumGoodSteps 
+      " Bad   Steps= " <<  fGrandSumBadSteps  
+      " Keep  Steps= " <<  fGrandSumKeepSteps << G4endl;
 }
 
 template <typename Real_v>
@@ -170,7 +199,7 @@ void TrialStats<Real_v>::PrintStats( bool full ) const
 
    if( full )
    {
-      PrintSums();
+      // PrintSums();
    }
 }
 

magneticfield/src/ScalarIntegrationDriver.cxx

@@ -71,6 +71,7 @@ ScalarIntegrationDriver::ScalarIntegrationDriver(double    hminimum,
 #ifdef GVFLD_STATS
   statsEnabled = true;
 #endif
+  fPrintDerived= false;
 
   fErrZeroStepCount = 0; // Counter for reporting zero step
 
@@ -230,45 +231,35 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
   }
 
   bool lastStep = false;
-  nstp          = 1;
+  nstp          = 0;
 
   // fpStepper->InitializeCharge( charge );   // May be added
   // OLD: fpStepper->GetEquationOfMotion()->InitializeCharge( charge );
 
-  double StartPosAr[3];
-  // ThreeVector StartPos( y[0], y[1], y[2] );
+
   do {
 #ifdef GVFLD_STATS
     gHistStepsLin->Fill(h);
     if (h > 0) gHistStepsLog->Fill(log(h));
 #endif
 
-    // StartPos= ThreeVector( y[0], y[1], y[2] );
-    StartPosAr[0] = y[0];
-    StartPosAr[1] = y[1];
-    StartPosAr[2] = y[2];
+    ThreeVector StartPos( y[0], y[1], y[2] );    
+    // double StartPosAr[3] = { y[0], y[1], y[2] };
 #ifdef GUDEBUG_FIELD
     double xSubStepStart = x;
     for (i = 0; i < nvar; i++) {
       ySubStepStart[i] = y[i];
     }
-// yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
-// yFldTrkStart.SetCurveLength(x);
 #endif
 
-    // Old method - inline call to Equation of Motion
-    //   fpStepper->RightHandSide( y, dydx );
-    // New method allows to cache field, or state (eg momentum magnitude)
-    // fpStepper->ComputeRightHandSide( y, charge, dydx );
-
-    // Back to simple, old method   - JA. 16 Oct 2015
     fpStepper->RightHandSideVIS(y, charge, dydx); // TODO: change to inline
     fNoTotalSteps++;
 
-// #ifdef CHECK_DYDX    
-    std::cout << "ScalarDriver::AccurateAdv>  RHS called with q= " << charge
-              << " at Position = " << y[0] << " y= " << y[1] << " z= " << y[2]
-              << " with Momentum = " << y[3] << " y= " << y[4] << " z= " << y[5] << " ";
+#ifdef CHECK_DYDX
+    if( partDebug )    
+       std::cout << "ScalarDriver::AccurateAdv>  RHS called with q= " << charge
+                 << " at Position = " << y[0] << " y= " << y[1] << " z= " << y[2]
+                 << " with Momentum = " << y[3] << " y= " << y[4] << " z= " << y[5] << " ";
           
     vecgeom::Vector3D<double> Bfield;
     double dydxAgn[ScalarFieldTrack::ncompSVEC];
@@ -276,11 +267,14 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
     using geant::units::tesla;
 
     equationPtr->EvaluateRhsReturnB(y, charge, dydxAgn, Bfield);
+
+    if( partDebug )
+    {     
           // (const double y[], double dydx[], double charge,
-    std::cout << " from B-field,  Bx= " << Bfield.x() / tesla << " By= " << Bfield.y() / tesla
-              << " Bz= " << Bfield.z() / tesla << " ";
-    std::cout << " gives Derivs dydx= :  x = " << dydx[0] << " y = " << dydx[1] << " z = " << dydx[2]        
-              << " px= " << dydx[3] << " py= " << dydx[4] << " pz= " << dydx[5] << std::endl;
+       std::cout << " from B-field,  Bx= " << Bfield.x() / tesla << " By= " << Bfield.y() / tesla
+                 << " Bz= " << Bfield.z() / tesla << " ";
+       std::cout << " gives Derivs dydx= :  x = " << dydx[0] << " y = " << dydx[1] << " z = " << dydx[2]        
+                 << " px= " << dydx[3] << " py= " << dydx[4] << " pz= " << dydx[5] << std::endl;
 /***                                                                                   
     std::cout << " Cross check: ";                                                                               
     std::cout << "   Derivs  dydx= :  x = " << dydx[0] << " y = " << dydx[1] << " z = " << dydx[2]        
@@ -288,80 +282,27 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
     std::cout << "   D/again dydx= :  x = " << dydxAgn[0] << " y = " << dydxAgn[1] << " z = " << dydxAgn[2]        
               << " px= " << dydxAgn[3] << " py= " << dydxAgn[4] << " pz= " << dydxAgn[5] << std::endl;
 ***/
-    std::cout << "   -diff   dydx= :  x = " << dydx[0] - dydxAgn[0] << " y = " << dydx[1] - dydxAgn[1]
-              << " z = " << dydx[2] - dydxAgn[2]
-              << " px= " << dydx[3] - dydxAgn[3]
-              << " py= " << dydx[4]-dydxAgn[4] << " pz= " << dydx[5]-dydxAgn[5] << std::endl;
-// #endif // of CHECK_DYDX
+       std::cout << "   -diff   dydx= :  x = " << dydx[0] - dydxAgn[0] << " y = " << dydx[1] - dydxAgn[1]
+                 << " z = " << dydx[2] - dydxAgn[2]
+                 << " px= " << dydx[3] - dydxAgn[3]
+                 << " py= " << dydx[4]-dydxAgn[4] << " pz= " << dydx[5]-dydxAgn[5] << std::endl;
+    }
+#endif // of CHECK_DYDX
 
-    // Perform the Integration
+    // if ( h > fMinimumStep)
+    // Changed 2019.03.04 : Only 'regular' steps with Runge-Kutta and error control are used.
+    //                      This is in order to emulate Vector mode exactly.
+    // I.e. suppressed 'quick advance' that accepted any estimated error in those small steps
+    
+    // Perform one step of Integration - with error control
     //
-    if ( true )
-    // if ( h > fMinimumStep) //  Changed 2019.03.04 => emulate Vector mode
-    {
-      // std::cout << "Calling       OneGoodStep " << std::endl;
-      OneGoodStep(y, charge, dydx, x, h, fEpsilonRelMax, hdid, hnext);
-      // std::cout << "Returned from OneGoodStep" << std::endl;
-
-      //--------------------------------------
-      lastStepSucceeded = (hdid == h);