diff --git a/magneticfield/inc/Geant/BaseRkIntegrationDriver.h b/magneticfield/inc/Geant/BaseRkIntegrationDriver.h
index 94e0d07251722e8aba4957516dc01744ef385ca5..9633da37790cdfb91499f9ce086bea1a6c44bf48 100644
--- a/magneticfield/inc/Geant/BaseRkIntegrationDriver.h
+++ b/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() );
diff --git a/magneticfield/inc/Geant/FormattedReporter.h b/magneticfield/inc/Geant/FormattedReporter.h
index acb05810cc473131fe8fc2a35d38e1d8e18f2ebb..71da8b0e171fb01ba032183fde0c8bda277d0cde 100644
--- a/magneticfield/inc/Geant/FormattedReporter.h
+++ b/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.
}
diff --git a/magneticfield/inc/Geant/OldIntegrationDriver.h b/magneticfield/inc/Geant/OldIntegrationDriver.h
index 560af2abc6d90f6b8b3b0ed7f14e1f3e8565b40b..a1fab6edebe4bffb4a42ac49edfe05f1102eb8c2 100644
--- a/magneticfield/inc/Geant/OldIntegrationDriver.h
+++ b/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)
diff --git a/magneticfield/inc/Geant/RollingIntegrationDriver.h b/magneticfield/inc/Geant/RollingIntegrationDriver.h
index 0510b82f6b6b94c4bb718ab1cb721415dcba3ed0..38eeeb8ade24275b2c34e14a39ffe70b82da1d02 100644
--- a/magneticfield/inc/Geant/RollingIntegrationDriver.h
+++ b/magneticfield/inc/Geant/RollingIntegrationDriver.h
@@ -53,7 +53,7 @@
#include "Geant/AuxVecMethods.h"
// ---- Flags for debugging and diagnostics only - off in benchmarks!
-#define DRIVER_PRINT_PROGRESS 1
+// #define DRIVER_PRINT_PROGRESS 1
#define DRIVER_DIAGNOSTICS 1
#include "ErrorEstimatorSixVec.h"
@@ -238,7 +238,7 @@ protected:
fErrcon = Math::Pow(fMaxSteppingIncrease / fSafetyFactor, 1.0 / GetPowerGrow() );
fShrinkThresh = Math::Pow(fMaxSteppingDecrease / fSafetyFactor, 1.0 / GetPowerShrink() );
- std::cout << "SimpleIntegrationDriverComputAndSetErrcon(): fErrcon = " << fErrcon
+ std::cout << "RollingIntegrationDriverComputAndSetErrcon(): fErrcon = " << fErrcon
<< " from: maxStepIncrease = " << fMaxSteppingIncrease
<< " fSafetyFactor = " << fSafetyFactor
<< " power-grow = " << GetPowerGrow() << std::endl;
@@ -368,7 +368,7 @@ private:
// ---------------------------------------------------------------
// Compilation constants
#ifdef CONST_DEBUG
- const bool partDebug = true; // Enforce debugging output
+ const bool partDebug = false; // true; // 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
@@ -638,7 +638,7 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
using vecCore::Get;
-#ifdef DRIVER_PRINT_PROGRESS
+#if defined(DRIVER_PRINT_PROGRESS) || defined(CHECK_ONE_LANE)
using FormattedReporter::ReportRowOfDoubles;
using FormattedReporter::ReportRowOfSquareRoots;
using FormattedReporter::ReportManyRowsOfDoubles;
@@ -649,11 +649,12 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
using PrintDriverProgress::ReportConditionLanes;
#endif
+#ifdef DRIVER_PRINT_PROGRESS
+ const Real_v xStart= x; // Remember
static std::atomic<unsigned int> numCallsRidKS(0);
+ int currentCallNo= numCallsRidKS++;
-#ifdef DRIVER_PRINT_PROGRESS
if (partDebug) { cout << "\n" << endl; }
- int currentCallNo= numCallsRidKS++;
#endif
Real_v yerr[Nvar];
@@ -661,8 +662,6 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
Real_v yStepEnd[Nvar]; // Could reuse yFinal for this ?? 2018.01.08
Real_v hStep = htry; // Set stepsize to the initial trial value
- const Real_v xStart= x; // Remember
-
// Store values for last good step - for use when current step is not good !
// Real_v hStepLastGood(0.0), errmaxSqLastGood(-1.0);
@@ -697,18 +696,19 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
#ifdef CHECK_ONE_LANE
const int trackToPrint = IntegrationDriverConstants::GetInstance()->GetTrackToCheck();
#ifndef CONST_DEBUG
- partDebug = ( laneToCheck != -1 );
+ partDebug = false;
+ // partDebug = ( laneToCheck != -1 );
#endif
// bool printNow = (laneToCheck >= 0) && vecCore::Get( htry, laneToCheck ) > 0.0 ;
- if( laneToCheck != -1 ) {
+ if( laneToCheck != -1 && partDebug ) {
cout << "* RID:AccAdv Global Vars> trackToPrint = " << trackToPrint << " l2c / laneToCheck = " << laneToCheck;
cout << " Args> h[l2c] = " << vecCore::Get( htry, laneToCheck );
cout << endl;
}
#endif
-#ifdef DRIVER_PRINT_PROGRESS
- if( partDebug ) {
+#ifdef DRIVER_PRINT_PROGRESS
+ if( false ) { // partDebug ) {
std::cout << "RID: KeepStepping called. Call number " << currentCallNo << std::endl;
cout << " -- Input Lane Step statistics." << std::endl;
@@ -756,23 +756,34 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
{
fpStepper ->RightHandSideInl( yStepStart, charge, dydx ); // TODO: change to inline
//---------****************-----------------------
- cout << "RID> Recalculated dy/ds ";
- if( !anyProgressLastStep ) { cout << " (likely unnecessary) "; }
- else { cout << " (believed needed) "; }
+ if( partDebug ) {
+ cout << "RID> Recalculated dy/ds ";
+ if( !anyProgressLastStep ) { cout << " (likely unnecessary) "; }
+ else { cout << " (believed needed) "; }
+ }
} else {
- if( iter > 0 ) { cout << "RID> Kept values of dy/ds from previous step "; }
+ if( partDebug )
+ if( iter > 0 ) { cout << "RID> Kept values of dy/ds from previous step "; }
}
- cout << " - iter = " << iter << " anyProgressLastStep = " << anyProgressLastStep << std::endl;
- ReportManyRowsOfDoubles("Current X/P", yStepStart, 6 );
- ReportRowOfDoubles("Charge", charge, 6 );
- ReportManyRowsOfDoubles("d/ds [X,P] ", dydx, 6 );
+
+#ifdef DRIVER_PRINT_PROGRESS
+ if( partDebug ) {
+ cout << " - iter = " << iter << " anyProgressLastStep = " << anyProgressLastStep << std::endl;
+ ReportManyRowsOfDoubles("Current X/P", yStepStart, 6 );
+ ReportRowOfDoubles("Charge", charge, 6 );
+ ReportManyRowsOfDoubles("d/ds [X,P] ", dydx, 6 );
+ }
+#endif
itersLeft--;
iter++;
+
+#ifdef DRIVER_PRINT_PROGRESS
if (partDebug) cout << " KeepStepping - iteration = " << iter << " "
<< " ( total iterations = " << (tot_no_trials + iter) << " ) "
<< endl;
-
+#endif
+
fpStepper->StepWithErrorEstimate(yStepStart, dydx, charge, hStep, yStepEnd, yerr); // CAREFUL -> changes for others ?
//******************************
@@ -801,17 +812,26 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
#else
// errmax_sq = fErrorEstimator.EstimateError( yerr, hStep, magmomInit_sq ); // Older ver: yStepEnd );
errmax_sq = fErrorEstimator.EstimateError( yerr, hStep, magmom_sq);
- // errmax_sq = fErrorEstimator.EstimateError( yerr, hStep, yStepEnd );
+ // errmax_sq = fErroorEstimator.EstimateError( yerr, hStep, yStepEnd );
#endif
goodStep = Active && (errmax_sq <= 1.0);
progressedLane = progressedLane || goodStep;
- laneStats.Update( Active, goodStep );
+ laneStats.Update( Active, goodStep );
+
+#ifdef DRIVER_PRINT_PROGRESS
+ // Checks 2019.03.21
+ if( false ) { // partDebug ) {
+ cout << " -- Updated Lane Step statistics - Active = " << Active << " goodStep = " << goodStep << std::endl;
+ laneStats.PrintStats( true );
+ ReportRowOfBools<Real_v>("Active(bef)", Active);
+ ReportRowOfBools<Real_v>("goodStep", goodStep);
+ }
+#endif
// vecCore::MaskedAssign( hStepLastGood, goodStep, hStep ); // if( goodStep[i] ) hStepLastGood[i] = hStep[i];
// vecCore::MaskedAssign( errmaxSqLastGood, goodStep, errmax_sq ); // if( goodStep[i] ) errmaxSqLastGood[i] = errmax_sq[i];
- // ****> Moved below ==> only used if falling through the while loop
// Update hdid for successful steps
Real_v hAdvance= vecCore::Blend( goodStep, hStep, Real_v(0.0) );
@@ -819,10 +839,12 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
Real_v xNext = x + hAdvance; // Updates only good steps -- 2019.01.07
+ // Bool_v laneStepNumOver = ( laneStats.GetNoGoodSteps() >= Index_v( GetMaxNoSteps() ) ) ; // Too many iterations
+
// fSmallestFraction
Bool_v endingIntegration = ( xEnd - xNext < geant::units::perMillion * htry ) // Care for very small diff
- || ( laneStats.GetNoSteps() > Index_v( GetMaxNoSteps() ) ) ; // Too many iterations
-
+ || ( laneStats.GetNoGoodSteps() >= Index_v( GetMaxNoSteps() ) ) ; // Too many iterations
+
Bool_v laneDone = ( endingIntegration || finishedLane ); // WAS = (goodStep || finished);
integrationDone = integrationDone || endingIntegration;
@@ -842,7 +864,8 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
ReportRowOfBools<Real_v>("progressedLane", progressedLane);
ReportRowOfBools<Real_v>("endingIntegration", endingIntegration);
ReportRowOfBools<Real_v>("laneDone", laneDone);
- ReportRowOfInts<Real_v>("numStep", laneStats.GetNoSteps() ); // Was nstp);
+ ReportRowOfInts<Real_v>("numStep/Good", laneStats.GetNoGoodSteps() );
+ ReportRowOfInts<Real_v>("numStep", laneStats.GetNoSteps() ); // Was nstp);
cout << " ====================================================================" << endl;
}
#endif
@@ -878,14 +901,16 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
// So must do all the work here!
#endif
-#ifdef DRIVER_PRINT_PROGRESS
- std::cout << "Loop end condition: " << " enoughFinished= " << enoughFinished
- << " composed using : "
- << " someDone= " << someDone
- << " allDone= " << allDone
- << " alreadySomeDone = " << alreadySomeDone << " (from last iteration) "
- << " mustFinishAllStill = " << mustFinishAllStill
- << std::endl;
+#ifdef DRIVER_PRINT_PROGRESS
+ if( partDebug ) {
+ std::cout << "Loop end condition: " << " enoughFinished= " << enoughFinished
+ << " composed using : "
+ << " someDone= " << someDone
+ << " allDone= " << allDone
+ << " alreadySomeDone = " << alreadySomeDone << " (from last iteration) "
+ << " mustFinishAllStill = " << mustFinishAllStill
+ << std::endl;
+ }
#endif
// alreadySomeDone = someDone; // ---> Moved to just before end of while loop
@@ -899,13 +924,13 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
ReportOneLane ( hStep, x, epsPosition, errpos_sq, errmom_sq, errmax_sq, laneDone,
allDone, iter, tot_no_trials, laneToCheck, trackToPrint,
"RollingID" );
- 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 << " Track to check " << trackToPrint << " laneToCheck = " << laneToCheck << std::endl;
+ ReportRowOfSquareRoots("ErrPos", errpos_sq );
+ ReportRowOfSquareRoots("ErrMom", errmom_sq );
+ ReportRowOfSquareRoots("ErrMax", errmax_sq );
+ ReportManyRowsOfDoubles("yErr", yerr, Nvar);
+
std::cout << "RID: Status after stepper call ----------------------------------------------" << std::endl;
ReportManyRowsOfDoubles("RealStartX/P", yStart, 6 );
FormattedReporter::FullReport(yStepStart, charge, dydx, hStep, yStepEnd, yerr, errmax_sq, Active, goodStep);
@@ -940,7 +965,7 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
}
#endif
- if( mustFinishAllStill ){
+ if( mustFinishAllStill && partDebug ){
std::cout << "RollingID> Now in 'FinishAll' mode ******** " << std::endl;
}
@@ -953,32 +978,13 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
// std::cout << "** Exiting loop - Report of status -------------------------------------- " << std::endl;
// FormattedReporter::FullReport(yStepStart, charge, dydx, hStep, yStepEnd, yerr, errmax_sq, Active, goodStep);
-
- // Version 7.x - part 1
- // cout << " Store Good Values calls in *** break *** (v7 part 1): 2019.05.13 17:45" << std::endl;
- // StoreGoodValues(yStepEnd, hStep, errmaxSqBreakOut, progressedLane, yFinal, hFinal, errmax_sqFinal);
- // StoreGoodValues(yStepEnd, hStep, errmax_sq, progressedLane, yFinal, hFinal, errmax_sqFinal);
- // StoreGoodValues(yStepEnd, hStepLastActive, errmaxSqLastActive, progressedLane, yFinal, hFinal, errmax_sqFinal);
-
- // cout << " Store Good Values calls in *** break *** (v7.1 part 1): 2019.05.13 17:45" << std::endl;
- // StoreGoodValues(yStepEnd, hStepLastActive, errmaxSqLastActive, progressedLane && goodStep, yFinal, hFinal, errmax_sqFinal);
- // StoreGoodValues(yStepStart, hStepLastActive, errmaxSqLastActive, progressedLane && !goodStep, yFinal, hFinal, errmax_sqFinal);
-
- // cout << " Store Good Values calls in *** break *** (v7.3 part 1): 2019.05.20 19:25" << std::endl;
- // StoreGoodValues(yStepEnd, hStep, errmax_sq, progressedLane && goodStep, yFinal, hFinal, errmax_sqFinal);
- // StoreGoodValues(yStepStart, hStepLastActive, errmaxSqLastActive, progressedLane && !goodStep, yFinal, hFinal, errmax_sqFinal);
-
- // cout << " Store Good Values calls in *** break *** (v7.4 part 1): 2019.05.20 19:40" << std::endl;
- // for( unsigned int i=0; i<Nvar; i++)
- // vecCore::MaskedAssign( yStepStart[i], goodStep, yStepEnd[i] );
- // StoreGoodValues(yStepStart, hStepLastActive, errmaxSqLastActive, Bool_v(true), yFinal, hFinal, errmax_sqFinal);
-
+ if (partDebug)
+ cout << " Store Good Values calls in *** break *** (v7.5 part 1): 2019.05.20 19:45" << std::endl;
- cout << " Store Good Values calls in *** break *** (v7.5 part 1): 2019.05.20 19:45" << std::endl;
StoreGoodValues( yStepEnd, hStepLastActive, errmaxSqLastActive, goodStep, yFinal, hFinal, errmax_sqFinal);
Bool_v workingOnIt = !goodStep && ( lastActive || progressedLane ) ;
if( ! vecCore::MaskEmpty( workingOnIt ) ) {
- ReportRowOfBools<Real_v>("workingOnIt", workingOnIt);
+ // ReportRowOfBools<Real_v>("workingOnIt", workingOnIt);
StoreGoodValues( yStepStart, hStepLastActive, errmaxSqLastActive, workingOnIt, yFinal, hFinal, errmax_sqFinal);
}
break;
@@ -997,8 +1003,8 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
Bool_v stepSizeUnderflow = Active && (xnew == x);
-#ifdef CHECK_ONE_LANE
- if( printLane ) {
+#ifdef CHECK_ONE_LANE
+ if( printLane && partDebug ) {
std::cout << "RID: After chance to break. (Not enoughFinished.) Remaining lanes represent work: " << std::endl;
ReportRowOfBools<Real_v>("laneDone", laneDone);
ReportRowOfBools<Real_v>("Active (still)", Active);
@@ -1051,7 +1057,7 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
toContinue = itersLeft > 0 && (!vecCore::MaskFull(finishedLane));
if( !toContinue ){
-#ifdef DRIVER_PRINT_PROGRESS
+#ifdef DRIVER_PRINT_PROGRESS
std::cout << "WARNING> RollingIntegrationDriver::KeepStepping: exiting at the Bottom of while." << std::endl;
// FormattedReporter::FullReport(yStepStart, charge, dydx, hStep, yStepEnd, yerr, errmax_sq, Active, goodStep);
#endif
@@ -1067,7 +1073,7 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
bool fellThrough= !toContinue;
- if( !toContinue )
+ if( !toContinue && partDebug )
std::cout << "WARNING> RollingIntegrationDriver::KeepStepping: exiting at the Bottom of while." << std::endl;
// Only ensure that
@@ -1085,10 +1091,12 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
if( fellThrough )
{
-#ifdef DRIVER_PRINT_PROGRESS
- cout << "RiD::KeepStepping> Fell through while loop" << endl;
- cout << " Store Good Values calls - version 7 (fall through): 2019.05.13 17:45" << std::endl;
- ReportRowOfDoubles("errmaxSqFallThru = ", errmaxSqFallThru );
+#ifdef DRIVER_PRINT_PROGRESS
+ if (partDebug) {
+ cout << "RiD::KeepStepping> Fell through while loop" << endl;
+ cout << " Store Good Values calls - version 7 (fall through): 2019.05.13 17:45" << std::endl;
+ ReportRowOfDoubles("errmaxSqFallThru = ", errmaxSqFallThru );
+ }
#endif
// Try 7: 2019.05.09 17:45
StoreGoodValues(yStepStart, hStepLastActive, errmaxSqLastActive, progressedLane, yFinal, hFinal, errmax_sqFinal);
@@ -1097,13 +1105,15 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
// Ideas -- old
// StoreGoodValues(yStepStart, hStepLastGood, errmaxSqLastGood, progressedLane && !goodStep, yFinal, hFinal, errmax_sqFinal);
} else {
-#ifdef DRIVER_PRINT_PROGRESS
- ReportRowOfDoubles("errmaxSqBreakOut = ", errmaxSqBreakOut );
+#ifdef DRIVER_PRINT_PROGRESS
+ if (partDebug) {
+ ReportRowOfDoubles("errmaxSqBreakOut = ", errmaxSqBreakOut );
+ }
#endif
}
#ifdef DRIVER_PRINT_PROGRESS
- if( 1 ) {
+ if( partDebug ) {
cout << " Check LastGood values vs. hstep-last/BreakOut" << std::endl;
ReportRowOfDoubles("hStep =", hStep );
// ReportRowOfDoubles("hStepLastGood = ", hStepLastGood );
@@ -1111,7 +1121,7 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
// template <class Real_v> using FormattedReporter::ReportRowOfInts<Real_v>;
cout << " Lane Step statistics." << std::endl;
- laneStats.PrintStats(true); // ToDo : false + end-of-job PrintSums();
+ laneStats.PrintStats(false); // ToDo : false + end-of-job PrintSums();
}
#endif
@@ -1188,7 +1198,7 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::
#ifdef DRIVER_PRINT_PROGRESS
bool KPSreport = true;
- if ( 1 /*partDebug*/ && KPSreport)
+ if ( partDebug && KPSreport)
{
ReportRowOfDoubles("x: Updated =", x );
// #ifndef CHECK_STRETCH_FACTOR
@@ -1610,13 +1620,13 @@ void RollingIntegrationDriver</*Real_v,*/ T_Stepper, Nvar>::StoreOutput(const Re
template <class T_Stepper, unsigned int Nvar>
template <class Real_v>
void RollingIntegrationDriver<T_Stepper, Nvar>::AccurateAdvance
-(const FieldTrack yInput[],
- const double hstep[],
- const double charge[],
- // double epsilon, // Can be scalar or varying
- FieldTrack yOutput[],
- bool stillOK[],
- int nTracks) const
+ ( const FieldTrack yInput[],
+ const double hstep[],
+ const double charge[],
+ // double epsilon, // Can be scalar or varying
+ FieldTrack yOutput[],
+ bool stillOK[],
+ int nTracks) const
{
// Built on original AccurateAdvance. Takes buffer stream of nTracks
// Converts them to Vc vectors for processing
@@ -1742,13 +1752,12 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::AccurateAdvance
x = x1;
TrialStats<Real_v> laneStats;
-
// trialStats<double> trackStats[nTracks]; // To store them for more examination ... ?
while ( ! vecCore::MaskFull( laneUnemployed )
&&
( ! vecCore::MaskEmpty( ! laneUnemployed
- && ( laneStats.GetNoSteps() <= Index_v(GetMaxNoSteps()) )
+ && ( laneStats.GetNoGoodSteps() < Index_v(GetMaxNoSteps()) )
&& (x < x2)
&& (!isLastStepLane)
)
@@ -1771,11 +1780,14 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::AccurateAdvance
fpStepper ->RightHandSideInl( y, chargeLane, dydx ); // TODO: change to inline
//---------****************-----------------------
-
- ReportRowOfDoubles("hTry(before)", hTry);
bool lastBatch = ( idNext == nTracks);
-
- std::cout << " Just before call to KeepStepping: lastBatch = " << lastBatch << " idNext = " << idNext << " nTracks= " << nTracks << std::endl;
+
+#ifdef DRIVER_PRINT_PROGRESS
+ if( partDebug ) {
+ ReportRowOfDoubles("hTry(before)", hTry);
+ std::cout << " Just before call to KeepStepping: lastBatch = " << lastBatch << " idNext = " << idNext << " nTracks= " << nTracks << std::endl;
+ }
+#endif
KeepStepping<Real_v>( y,
dydx, // Returns last value !
@@ -1864,9 +1876,9 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::AccurateAdvance
hnext = Min(xRemains, hnext); // Ensure that the next step does not overshoot
#ifdef DRIVER_PRINT_PROGRESS
- if ( 1 /*partDebug*/ ) {
+ if ( partDebug ) {
cout << "AccurateAdvance: hnext = " << hnext << " to replace hTry = " << hTry << endl;
- ReportRowOfDoubles("SID/AccAdv: hNext=", hnext);
+ ReportRowOfDoubles("RID/AccAdv: hNext=", hnext);
}
#endif
hTry = hnext;
@@ -1941,9 +1953,11 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::AccurateAdvance
// Common - whether new work exists or not
-#ifdef DRIVER_PRINT_PROGRESS
- cout << "Reseting step stats for lane " << i << " ";
- laneStats.PrintLaneStats(i);
+#ifdef DRIVER_PRINT_PROGRESS
+ if( partDebug ) {
+ cout << "Reseting step stats for lane " << i << " ";
+ laneStats.PrintLaneStats(i);
+ }
#endif
laneStats.ResetLane(i);
@@ -2105,6 +2119,10 @@ void RollingIntegrationDriver<T_Stepper, Nvar>::AccurateAdvance
} // end of while loop
+#ifdef DRIVER_PRINT_PROGRESS
+ laneStats.PrintSums();
+#endif
+
if( numBadSize > 0) {
cout << "Copying yIn => yOut for " << numBadSize << " 'bad' tracks - i.e. hStep <= 0.0 " << endl;
for (int i = 0; i < nTracks; ++i) {
diff --git a/magneticfield/inc/Geant/SimpleIntegrationDriver.h b/magneticfield/inc/Geant/SimpleIntegrationDriver.h
index 8007196d4830c051f1543d9d645cfe71916b9f3b..ee9541746e4f70269e5a3998816096ad49484cb5 100644
--- a/magneticfield/inc/Geant/SimpleIntegrationDriver.h
+++ b/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");
diff --git a/magneticfield/inc/Geant/TrialStats.h b/magneticfield/inc/Geant/TrialStats.h
index 4f5d4a7495194177fa899d80191ab9b39c08c8a5..7f7cbcd163b5dfd9621acc3513d45a143c391384 100644
--- a/magneticfield/inc/Geant/TrialStats.h
+++ b/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();
}
}
diff --git a/magneticfield/src/ScalarIntegrationDriver.cxx b/magneticfield/src/ScalarIntegrationDriver.cxx
index 3387cb0b4900c469c818cf06ac5d1273fee80e01..dbec336859f77ca7b77fad5703036229352bf45c 100644
--- a/magneticfield/src/ScalarIntegrationDriver.cxx
+++ b/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);
-#ifdef GUDEBUG_FIELD
- if (dbg > 2) {
- PrintStatus(ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
- }
-#endif
- } else {
- ScalarFieldTrack yFldTrk(ThreeVector(0, 0, 0), ThreeVector(0, 0, 0), charge);
- // double dchord_step;
- double dyerr_len_sq, dyerr_mom_rel_sq; // What to do with these ?
- yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
- yFldTrk.SetCurveLength(x);
-
- QuickAdvance(yFldTrk, dydx, h, /*dchord_step,*/ dyerr_len_sq, dyerr_mom_rel_sq);
- //----------------------------------------------------------------
-
- yFldTrk.DumpToArray(y);
-
-#ifdef GVFLD_STATS
- fNoSmallSteps++;
- fDyerrPosMaxSq = std::max(fDyerrPosMaxSq, dyerr_len_sq);
- fDyerrDirMaxSq = std::max(fDyerrDirMaxSq, dyerr_mom_rel_sq);
-
- fDyerrPos_smTot += std::sqrt(dyerr_len_sq);
- fSumH_sm += h; // Length total for 'small' steps
- if (nstp <= 1) {
- fNoInitialSmallSteps++;
- }
-#endif
+ OneGoodStep(y, charge, dydx, x, h, fEpsilonRelMax, hdid, hnext);
+
+ //--------------------------------------
+ lastStepSucceeded = (hdid == h);
#ifdef GUDEBUG_FIELD
- if (dbg > 1) {
- if (fNoSmallSteps < 2) {
- PrintStatus(ySubStepStart, x1, y, x, h, -nstp);
- }
- std::cout << "Another sub-min step, no " << fNoSmallSteps << " of " << fNoTotalSteps << " this time " << nstp
- << std::endl;
- PrintStatus(ySubStepStart, x1, y, x, h, nstp); // Only this
- std::cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h << " epsilon= " << fEpsilonRelMax
- << " hstep= " << hstep << " h= " << h << " hmin= " << fMinimumStep << std::endl;
- }
+ if (dbg > 2) { PrintStatus(ySubStepStart, xSubStepStart, y, x, h, nstp); }
#endif
- if (h == 0.0) {
- std::cerr << "ERROR in G4UIntegationDriver::AccurateAdvance: integration failure. " << std::endl;
- std::cerr << "Integration Step became Zero!" << std::endl;
- exit(1);
- }
- double dyerr_sq = std::min(dyerr_len_sq / (h * h), dyerr_mom_rel_sq);
- double dyerr = std::sqrt(dyerr_sq);
- // dyerr = dyerr_len / h;
- hdid = h;
- x += hdid;
-
- // Compute suggested new step
- hnext = ComputeNewStepSize(dyerr * fInvEpsilonRelMax, h);
-
- // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/fEpsilonRelMax, h);
- lastStepSucceeded = (dyerr <= fEpsilonRelMax);
- }
if (lastStepSucceeded) {
noFullIntegr++;
@@ -373,11 +314,8 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
#ifdef GUDEBUG_FIELD
if ((dbg > 0) && (dbg <= 2) && (nstp > nStpPr)) {
- if (nstp == nStpPr) {
- std::cout << "***** Many steps ****" << std::endl;
- }
- std::cout << "MagIntDrv: ";
- std::cout << "hdid=" << std::setw(12) << hdid << " "
+ if (nstp == nStpPr) { std::cout << "***** Many steps ****" << std::endl; }
+ std::cout << "MagIntDrv: " << "hdid=" << std::setw(12) << hdid << " "
<< "hnext=" << std::setw(12) << hnext << " "
<< "hstep=" << std::setw(12) << hstep << " (requested) " << std::endl;
PrintStatus(ystart, x1, y, x, h, (nstp == nStpPr) ? -nstp : nstp);
@@ -385,10 +323,12 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
#endif
// Check the endpoint
+ /***
const double edx = y[0] - StartPosAr[0];
const double edy = y[1] - StartPosAr[1];
- const double edz = y[2] - StartPosAr[2];
- double endPointDist2 = edx * edx + edy * edy + edz * edz; // (EndPos-StartPos).Mag();
+ const double edz = y[2] - StartPosAr[2]; ***/
+ double endPointDist2 = // edx * edx + edy * edy + edz * edz;
+ (EndPos-StartPos).Mag2();
if (endPointDist2 >= hdid * hdid * (1. + 2. * perMillion)) {
double endPointDist = std::sqrt(endPointDist2);
fNoBadSteps++;
@@ -396,14 +336,17 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
// Issue a warning only for gross differences -
// we understand how small difference occur.
if (endPointDist >= hdid * (1. + perThousand)) {
+ int prdbg= 0;
#ifdef GUDEBUG_FIELD
- if (dbg) {
- WarnEndPointTooFar(endPointDist, hdid, fEpsilonRelMax, dbg);
- std::cerr << " Total steps: bad " << fNoBadSteps << " good " << noGoodSteps << " hdid= " << hdid
- << std::endl;
+ prdbg= dbg;
+#endif
+ WarnEndPointTooFar(endPointDist, hdid, fEpsilonRelMax, prdbg);
+ std::cerr << " Total steps: bad " << fNoBadSteps << " good " << noGoodSteps << " hdid= " << hdid
+ << std::endl;
+ if (prdbg) {
PrintStatus(ystart, x1, y, x, hstep, no_warnings ? nstp : -nstp);
}
-#endif
+// #endif
no_warnings++;
}
} else {
@@ -456,7 +399,9 @@ bool ScalarIntegrationDriver::AccurateAdvance(const ScalarFieldTrack &yInput, do
#endif
}
}
- } while (((nstp++) <= fMaxNoSteps) && (x < x2) && (!lastStep));
+ } while (((++nstp) < fMaxNoSteps )
+ && (x < x2)
+ && (!lastStep));
// Have we reached the end ?
// --> a better test might be x-x2 > an_epsilon
@@ -660,10 +605,9 @@ void ScalarIntegrationDriver::OneGoodStep(double y[], // InOut
errmom_sq *= inv_eps_vel_sq;
errmax_sq = std::max(errpos_sq, errmom_sq); // Square of maximum error
-// #ifdef GUDEBUG_FIELD
+#ifdef GUDEBUG_FIELD
if( fPrintDerived
|| fVerboseLevel>2)
- // if(fVerboseLevel>2)
{
int op= std::cout.precision(8);
std::cout << "ScalarIntDrv: 1-good-step - track-num=" << this->GetTrackNumber() << " iter = " << iter << " " << std::endl
@@ -680,7 +624,7 @@ void ScalarIntegrationDriver::OneGoodStep(double y[], // InOut
std::cout.precision(op);
PrintStatus(y, x, ytemp, x + h, h, iter);
}
-// #endif
+#endif
#ifdef CHECK_ONE_LANE
// std::cout << "ScalarIntDrv> print 'Derived' quantities = " << fPrintDerived << std::endl;
@@ -753,21 +697,25 @@ void ScalarIntegrationDriver::OneGoodStep(double y[], // InOut
#endif
// Compute size of next Step
- if( fPrintDerived ) { std::cout << " errmax / fErrcon = " << std::sqrt( errmax_sq ) / fErrcon << " "; }
+ // if( fPrintDerived ) { std::cout << " errmax / fErrcon = " << std::sqrt( errmax_sq ) / fErrcon << " "; }
+
+ // fPrintDerived= false;
if (errmax_sq > fErrcon * fErrcon) {
double growthFac = GetSafety() * Math::Pow(errmax_sq, 0.5 * GetPowerGrow());
hnext = h * growthFac; // GetSafety() * h * Math::Pow(errmax_sq, 0.5 * GetPowerGrow());
- if( fPrintDerived ) { std::cout << " hnext = " << hnext << " ( formula - grow ) stretch-factor = " << growthFac << " "; }
+ // if( fPrintDerived ) { std::cout << " hnext = " << hnext << " ( formula - grow ) stretch-factor = " << growthFac << " "; }
} else {
hnext = fMaxSteppingIncrease * h; // No more than a factor of 5 increase
- if( fPrintDerived ) { std::cout << " hnext = " << hnext << " (else MaxIncrease ) "; }
+ // if( fPrintDerived ) { std::cout << " hnext = " << hnext << " (else MaxIncrease ) "; }
}
x += (hdid = h);
- if( fPrintDerived ) { std::cout << " [ errCon = " << fErrcon << " ] "
- << " errmax / errcon = " << std::sqrt( errmax_sq ) / fErrcon << " "
- << "xEnd = " << x << " hdid = " << hdid << std::endl; }
+ if( false ) { // fPrintDerived ) {
+ std::cout << " [ errCon = " << fErrcon << " ] "
+ << " errmax / errcon = " << std::sqrt( errmax_sq ) / fErrcon << " "
+ << "xEnd = " << x << " hdid = " << hdid << std::endl;
+ }
for (int k = 0; k < fNoIntegrationVariables; k++) {
y[k] = ytemp[k];