Attachment #303 for problem #1696

View | Details | Raw Unified | Return to problem 1696
Collapse All | Expand All

Lines 16-22 Link Here

(-)History (+9 lines)
16	----------------------------------------------------------	16	----------------------------------------------------------
17	* Reverse chronological order (last date on top), please *	17	* Reverse chronological order (last date on top), please *
18	----------------------------------------------------------	18	----------------------------------------------------------
		19	Jan 16, 2015 J.Apostolakis - field-V10-01-02 (fixed 2)
		20	--------------------------
		21	- Revert uninteded development in G4MagIntegratorDriver
19		22
		23	Jan 15, 2015 J.Apostolakis - field-V10-01-00, 01 (fixed)
		24	--------------------------
		25	- Small refinements in G4FieldTrack - access rest mass, avoid div by 0
		26	Avoid division by zero for Unit direction in FieldTrack
		27	Added method to get Rest Mass
		28
20	Nov 03, 2014 G.Cosmo - field-V10-00-04	29	Nov 03, 2014 G.Cosmo - field-V10-00-04
21	--------------------	30	--------------------
22	- Moved constructors and simple methods to in line in G4ChargeState and	31	- Moved constructors and simple methods to in line in G4ChargeState and

Lines 71-77 Link Here

(-)include/G4FieldTrack.icc (-1 / +2 lines)
71	SixVector[4] = pMomentum.y();	71	SixVector[4] = pMomentum.y();
72	SixVector[5] = pMomentum.z();	72	SixVector[5] = pMomentum.z();
73		73
74	fMomentumDir = pMomentum.unit();	74	fMomentumDir = (pMomentum.mag2() > 0.0) ?
		75	pMomentum.unit() : G4ThreeVector( 0.0, 0.0, 0.0 );
75		76
76	fDistanceAlongCurve= s_curve;	77	fDistanceAlongCurve= s_curve;
77		78

Lines 123-128 Link Here

(-)include/G4FieldTrack.hh (+1 lines)
123	inline G4double GetProperTimeOfFlight() const;	123	inline G4double GetProperTimeOfFlight() const;
124	inline G4double GetKineticEnergy() const;	124	inline G4double GetKineticEnergy() const;
125	inline G4double GetCharge() const;	125	inline G4double GetCharge() const;
		126	inline G4double GetRestMass() const { return fRestMass_c2; }
126	// Accessors.	127	// Accessors.
127		128
128	inline void SetPosition(G4ThreeVector nPos);	129	inline void SetPosition(G4ThreeVector nPos);

Lines 42-48 Link Here

(-)src/G4FieldTrack.cc (-1 / +2 lines)
42	<< G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag(); // mom magnitude	42	<< G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag(); // mom magnitude
43	os << " Ekin= " << SixVec.fKineticEnergy ;	43	os << " Ekin= " << SixVec.fKineticEnergy ;
44	os << " m0= " << SixVec.fRestMass_c2;	44	os << " m0= " << SixVec.fRestMass_c2;
45	os << " Pdir= " << SixVec.fMomentumDir.mag();	45	os << " Pdir= " << SixVec.fMomentumDir.mag();
		46	os << " PolV= " << SixVec.GetPolarization();
46	os << " l= " << SixVec.GetCurveLength();	47	os << " l= " << SixVec.GetCurveLength();
47	os << " t_lab= " << SixVec.fLabTimeOfFlight;	48	os << " t_lab= " << SixVec.fLabTimeOfFlight;
48	os << " t_proper= " << SixVec.fProperTimeOfFlight ;	49	os << " t_proper= " << SixVec.fProperTimeOfFlight ;




  static G4ThreadLocal G4int tot_no_trials=0; 
  const G4int max_trials=100; 

  const G4double max_reduction = 0.1 ; // Largest_reduction_factor for step size
  const G4double local_Pow_shrink = 0.5*GetPshrnk(); 

  // We choose to use a linear approximation of the Taylor series
  //       ( 1 + x ) ^ p =  1 + p * x - 0.5 * p * (1-p) x^2  + O( x^3 ) 
  // if the second order term is smaller than
  const G4double acc_limit= 0.01;
  
  // This limit is reached when
  //      0.5 * p * (1-p) x^2 = acc_limit
  // ie
  //            p * (1-p) x^2 = acc_limit * 2.0

  // NEW constants -- here for now only
  // The following value is the value of 'x' at which this limit is reached:
  const G4double  limit_fac = 0.5 * local_Pow_shrink * (1.0 - local_Pow_shrink) ;
  //  linear_limit_fac
  // To avoid the division here by limit_fac, we check below that
  //            limit_fac * x^2 < acc_limt 

  const G4double err_limit = std::pow( max_reduction / GetSafety(), 2/GetPshrnk() ); 
  // This is constant when the value of the 'shrink' power is constant
  //  TODO: These should become data members ( this one is expensive to calculate! )

  G4ThreeVector Spin(y[9],y[10],y[11]);
  G4double   spin_mag2 =Spin.mag2() ;
  G4bool     hasSpin= (spin_mag2 > 0.0); 

    }

    if ( errmax_sq <= 1.0 )  { break; } // Step succeeded. 

    // Step failed; compute the size of retrial Step.
    htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );

    if (htemp >= 0.1*h)  { h = htemp; }  // Truncation error too large,
    else  { h = 0.1*h; }                 // reduce stepsize, but no more
    G4double hnew = 0.0; 
    G4double x = errmax_sq - 1.0; 
    if (errmax_sq >= err_limit)  { 
       // if( limit_fac * x * x < acc_limt ) {
       //   hnew = 1.0 + local_Pow_shrink * x; 
       // }else{
          hnew = GetSafety()*h* std::pow( errmax_sq, local_Pow_shrink); // 0.5*GetPshrnk() );
       // }
    }  
    else  { hnew = max_reduction*h; } // Truncation error was too large,
                                      // reduce stepsize, but no more than x10
#endif

#if 1  // OLD_METHOD
    // Old method of calculation  - always used power
    htemp = GetSafety()*h* std::pow( errmax_sq, local_Pow_shrink); // 0.5*GetPshrnk() );

    if (htemp >= max_reduction*h)  { h = htemp; }  // Truncation error too large,
    else  { h =  max_reduction*h; }                 // reduce stepsize, but no more
                                         // than a factor of 10
#if 1  // NEW_METHOD
    assert( std::fabs(hnew - h) < 0.00001 * std::fabs(hnew+h) ); 
#endif
#endif

#if 1  // NEW_METHOD
    h= hnew; 
#endif

    xnew = x + h;
    if(xnew == x)
    {




     ----------------------------------------------------------
     * Reverse chronological order (last date on top), please *
     ----------------------------------------------------------
January  15, 2015 - J.Apostolakis (geomnav-V10-01-03)
---------------------------------
- G4PathFinder: Pass to equation of motion relevant properties of current
    particle: latest charge state (electric, magnetic dipole moment),
    PDG spin value, momentum and rest mass of current particle
    Part of corrections needed for Bugzilla issue 1696.
	
- Increased maximum number of Navigators from 8 to 16.
    
 rare cases of negative value to sqrt() in equation calculation for rmin/rmax
 intersection.
 Improved warning printout in G4Navigator::GetGlobalExitNormal().
    

Lines 216-222 Link Here

(-)include/G4PathFinder.hh (-1 / +1 lines)
216	G4int fNoActiveNavigators;	216	G4int fNoActiveNavigators;
217	G4bool fNewTrack; // Flag a new track (ensure first step)	217	G4bool fNewTrack; // Flag a new track (ensure first step)
218		218
219	static const G4int fMaxNav = 8; // rename to kMaxNoNav ??	219	static const G4int fMaxNav = 16; // rename to kMaxNoNav ??
220		220
221	// Global state (retained during stepping for one track)	221	// Global state (retained during stepping for one track)
222		222

Lines 164-170 Link Here

(-)include/G4MultiNavigator.hh (-1 / +1 lines)
164	// STATE Information	164	// STATE Information
165		165
166	G4int fNoActiveNavigators;	166	G4int fNoActiveNavigators;
167	static const G4int fMaxNav = 8; // rename to kMaxNoNav ??	167	static const G4int fMaxNav = 16; // rename to kMaxNoNav ??
168	G4VPhysicalVolume* fLastMassWorld;	168	G4VPhysicalVolume* fLastMassWorld;
169		169
170	// Global state (retained during stepping for one track	170	// Global state (retained during stepping for one track

Lines 1155-1160 Link Here

(-)src/G4PathFinder.cc (+9 lines)
1155	G4FieldTrack fieldTrack= initialState;	1155	G4FieldTrack fieldTrack= initialState;
1156	G4ThreeVector startPoint= initialState.GetPosition();	1156	G4ThreeVector startPoint= initialState.GetPosition();
1157		1157
		1158
		1159	G4EquationOfMotion* equationOfMotion =
		1160	(fpFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper())
		1161	->GetEquationOfMotion();
		1162
		1163	equationOfMotion->SetChargeMomentumMass( *(initialState.GetChargeState()),
		1164	initialState.GetMomentum().mag2(),
		1165	initialState.GetRestMass() );
		1166
1158	#ifdef G4DEBUG_PATHFINDER	1167	#ifdef G4DEBUG_PATHFINDER
1159	G4int prc= G4cout.precision(9);	1168	G4int prc= G4cout.precision(9);
1160	if( fVerboseLevel > 2 )	1169	if( fVerboseLevel > 2 )

Lines 61-69 Link Here

(-)src/G4FieldTrackUpdator.cc (-3 / +6 lines)
61	{	61	{
62	const G4DynamicParticle* ptDynamicParticle= trk->GetDynamicParticle();	62	const G4DynamicParticle* ptDynamicParticle= trk->GetDynamicParticle();
63		63
64	// The following properties must be updated ONCE for each new track (at least)	64	// The following properties must be updated 1) for each new track, and
65	ftrk->SetRestMass(ptDynamicParticle->GetMass());	65	ftrk->SetRestMass(ptDynamicParticle->GetMass());
66		66	// 2) Since ion can lose/gain electrons, this must be done at every step
		67
67	ftrk->UpdateState(	68	ftrk->UpdateState(
68	trk->GetPosition(),	69	trk->GetPosition(),
69	trk->GetGlobalTime(),	70	trk->GetGlobalTime(),
Lines 82-88 Link Here
82		83
83	ftrk->SetProperTimeOfFlight(trk->GetProperTime());	84	ftrk->SetProperTimeOfFlight(trk->GetProperTime());
84		85
85	ftrk->SetChargeAndMoments( ptDynamicParticle->GetCharge() );	86	ftrk->SetChargeAndMoments( ptDynamicParticle->GetCharge(),
		87	ptDynamicParticle->GetMagneticMoment());
		88	ftrk->SetPDGSpin( ptDynamicParticle->GetParticleDefinition()->GetPDGSpin() );
86	// The charge can change during tracking	89	// The charge can change during tracking
87	ftrk->SetSpin( ptDynamicParticle->GetPolarization() );	90	ftrk->SetSpin( ptDynamicParticle->GetPolarization() );
88	}	91	}




     * Reverse chronological order (last date on top), please *
     ----------------------------------------------------------

- January 15, 2015 J.Apostolakis (track-V10-01-01)
- Fixed FieldTrackUpdator to pass magnetic moment and PDG spin
  Is part of corrections needed for Bugzilla issue 1696.
     
- August 19, 2014 K. Genser (track-V10-00-04)
- Removed non cost requirement when using G4PhysicsModelCatalog getters and moved 
   implementations to icc file (requires global-V10-00-31)
                                 G4MagIntegratorStepper *pStepper,
                                 G4int                   numComponents,
                                 G4int                   statisticsVerbose)
 : fSmallestFraction( 1.0e-12 ), 
   fNoIntegrationVariables(numComponents), 
   fMinNoVars(12), 
   fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
   fStatisticsVerboseLevel(statisticsVerbose),
   fNoTotalSteps(0),  fNoBadSteps(0), fNoSmallSteps(0),
   fNoInitialSmallSteps(0), 
   fDyerr_max(0.0), fDyerr_mx2(0.0), 
   fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0), 
   fSumH_sm(0.0), fSumH_lg(0.0),
   fVerboseLevel(0)
 // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
 // is required. For proper time of flight and spin,  fMinNoVars must be 12
 RenewStepperAndAdjust( pStepper );
 fMinimumStep= hminimum;
 fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
 fVerboseLevel=2;
 if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
 {
   G4cout << "MagIntDriver version: Accur-Adv: "
          << "invE_nS, QuickAdv-2sqrt with Statistics "
          << " enabled "
          << " disabled "
          << G4endl;
 }
 if( fStatisticsVerboseLevel > 1 )
 {
   PrintStatisticsReport();
 }
                                G4double     hstep,
                                G4double     eps,
                                G4double hinitial )
 // Runge-Kutta driver with adaptive stepsize control. Integrate starting
 // values at y_current over hstep x2 with accuracy eps. 
 // On output ystart is replaced by values at the end of the integration 
 // interval. RightHandSide is the right-hand side of ODE system. 
 // The source is similar to odeint routine from NRC p.721-722 .
 G4int nstp, i, no_warnings=0;
 G4double x, hnext, hdid, h;
 static G4int dbg=1;
 static G4int nStpPr=50;   // For debug printing of long integrations
 G4double ySubStepStart[G4FieldTrack::ncompSVEC];
 G4FieldTrack  yFldTrkStart(y_current);
 G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
 G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC]; 
 G4double  x1, x2;
 G4bool succeeded = true, lastStepSucceeded;
 G4double startCurveLength;
 G4int  noFullIntegr=0, noSmallIntegr = 0 ;
 static G4ThreadLocal G4int  noGoodSteps =0 ;  // Bad = chord > curve-len 
 const  G4int  nvar= fNoVars;
 G4FieldTrack yStartFT(y_current);
 //  Ensure that hstep > 0
 //
 if( hstep <= 0.0 )
 { 
   if(hstep==0.0)
   {
     std::ostringstream message;
     message << "Proposed step is zero; hstep = " << hstep << " !";
     G4Exception("G4MagInt_Driver::AccurateAdvance()", 
                 "GeomField1001", JustWarning, message);
     return succeeded; 
   }
   else
   { 
     std::ostringstream message;
     message << "Invalid run condition." << G4endl
             << "Proposed step is negative; hstep = " << hstep << "." << G4endl
             << "Requested step cannot be negative! Aborting event.";
     G4Exception("G4MagInt_Driver::AccurateAdvance()", 
                 "GeomField0003", EventMustBeAborted, message);
     return false;
   }
 }
 y_current.DumpToArray( ystart );
 startCurveLength= y_current.GetCurveLength();
 x1= startCurveLength; 
 x2= x1 + hstep;
 if ( (hinitial > 0.0) && (hinitial < hstep)
   && (hinitial > perMillion * hstep) )
 {
    h = hinitial;
 }
 else  //  Initial Step size "h" defaults to the full interval
 {
    h = hstep;
 }
 x = x1;
 for (i=0;i<nvar;i++)  { y[i] = ystart[i]; }
 G4bool lastStep= false;
 nstp=1;
 do
 {
   G4ThreeVector StartPos( y[0], y[1], y[2] );
   G4double xSubStepStart= x; 
   for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
   yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
   yFldTrkStart.SetCurveLength(x);
   // Old method - inline call to Equation of Motion
   //   pIntStepper->RightHandSide( y, dydx );
   // New method allows to cache field, or state (eg momentum magnitude)
   pIntStepper->ComputeRightHandSide( y, dydx );
   fNoTotalSteps++;
   // Perform the Integration
   //      
   if( h > fMinimumStep )
   { 
     OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
     //--------------------------------------
     lastStepSucceeded= (hdid == h);   
     if (dbg>2) {
       PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
     }
   }
   else
   {
     G4FieldTrack yFldTrk( G4ThreeVector(0,0,0), 
                           G4ThreeVector(0,0,0), 0., 0., 0., 0. );
     G4double dchord_step, dyerr, dyerr_len;   // What to do with these ?
     yFldTrk.LoadFromArray(y, fNoIntegrationVariables); 
     yFldTrk.SetCurveLength( x );
     QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len ); 
     //-----------------------------------------------------
     yFldTrk.DumpToArray(y);    
     fNoSmallSteps++; 
     if ( dyerr_len > fDyerr_max)  { fDyerr_max= dyerr_len; }
     fDyerrPos_smTot += dyerr_len;
     fSumH_sm += h;  // Length total for 'small' steps
     if (nstp<=1)  { fNoInitialSmallSteps++; }
     if (dbg>1)
     {
       if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
       G4cout << "Another sub-min step, no " << fNoSmallSteps 
              << " of " << fNoTotalSteps << " this time " << nstp << G4endl; 
       PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
       G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h 
              << " epsilon= " << eps << " hstep= " << hstep 
              << " h= " << h << " hmin= " << fMinimumStep << G4endl;
     }
     if( h == 0.0 )
     { 
       G4Exception("G4MagInt_Driver::AccurateAdvance()",
                   "GeomField0003", FatalException,
                   "Integration Step became Zero!"); 
     }
     dyerr = dyerr_len / h;
     hdid= h;
     x += hdid;
     // Compute suggested new step
     hnext= ComputeNewStepSize( dyerr/eps, h);
     // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
     lastStepSucceeded= (dyerr<= eps);
   }
   if (lastStepSucceeded)  { noFullIntegr++; }
   else                    { noSmallIntegr++; }
   G4ThreeVector EndPos( y[0], y[1], y[2] );
   if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
   {
     if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
     G4cout << "MagIntDrv: " ; 
     G4cout << "hdid="  << std::setw(12) << hdid  << " "
            << "hnext=" << std::setw(12) << hnext << " " 
     PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp); 
   }
   // Check the endpoint
   G4double endPointDist= (EndPos-StartPos).mag(); 
   if ( endPointDist >= hdid*(1.+perMillion) )
   {
     fNoBadSteps++;
     // Issue a warning only for gross differences -
     // we understand how small difference occur.
     if ( endPointDist >= hdid*(1.+perThousand) )
     { 
       if (dbg)
       {
         WarnEndPointTooFar ( endPointDist, hdid, eps, dbg ); 
         G4cerr << "  Total steps:  bad " << fNoBadSteps
                << " good " << noGoodSteps << " current h= " << hdid
                << G4endl;
         PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);  
       }
       no_warnings++;
     }
   }
   else
   {
     noGoodSteps ++;
   } 
   //  Avoid numerous small last steps
   if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
   {
     // No more integration -- the next step will not happen
     lastStep = true;  
   }
   else
   {
     // Check the proposed next stepsize
     if(std::fabs(hnext) <= Hmin())
     {
       // If simply a very small interval is being integrated, do not warn
       if( (x < x2 * (1-eps) ) &&        // The last step can be small: OK
           (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
       {
         if(dbg>0)
         {
           WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );  
           PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
         }
         no_warnings++;
       }
       // Make sure that the next step is at least Hmin.
       h = Hmin();
     }
     else
     {
       h = hnext;
     }
     //  Ensure that the next step does not overshoot
     if ( x+h > x2 )
     {                // When stepsize overshoots, decrease it!
       h = x2 - x ;   // Must cope with difficult rounding-error
     }                // issues if hstep << x2
     if ( h == 0.0 )
     {
       // Cannot progress - accept this as last step - by default
       lastStep = true;
       if (dbg>2)
       {
         int prec= G4cout.precision(12); 
         G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
                << G4endl
                << "  Integration step 'h' became "
                << h << " due to roundoff. " << G4endl
                << "  Forcing termination of advance." << G4endl;
         G4cout.precision(prec);
       }          
     }
   }
 } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
    // Have we reached the end ?
    // --> a better test might be x-x2 > an_epsilon
 succeeded=  (x>=x2);  // If it was a "forced" last step
 for (i=0;i<nvar;i++)  { yEnd[i] = y[i]; }
 // Put back the values.
 y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
 y_current.SetCurveLength( x );
 if(nstp > fMaxNoSteps)
 {
   no_warnings++;
   succeeded = false;
   if (dbg)
   {
     WarnTooManyStep( x1, x2, x );  //  Issue WARNING
     PrintStatus( yEnd, x1, y, x, hstep, -nstp);
   }
 }
 if( dbg && no_warnings )
 {
   G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
   PrintStatus( yEnd, x1, y, x, hstep, nstp);
 }
 return succeeded;
                                   G4double h, G4double xDone,
                                   G4int nstp)
 static G4ThreadLocal G4int noWarningsIssued =0;
 const  G4int maxNoWarnings =  10;   // Number of verbose warnings
 std::ostringstream message;
 if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
 {
   message << "The stepsize for the next iteration, " << hnext
           << ", is too small - in Step number " << nstp << "." << G4endl
           << "The minimum for the driver is " << Hmin()  << G4endl
           << "Requested integr. length was " << hstep << " ." << G4endl
           << "The size of this sub-step was " << h     << " ." << G4endl
           << "The integrations has already gone " << xDone;
 }
 else
 {
   message << "Too small 'next' step " << hnext
           << ", step-no: " << nstp << G4endl
           << ", this sub-step: " << h     
           << ",  req_tot_len: " << hstep 
           << ", done: " << xDone << ", min: " << Hmin();
 }
 G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
             JustWarning, message);
 noWarningsIssued++;
                                 G4double x2end, 
                                 G4double xCurrent)
   std::ostringstream message;
   message << "The number of steps used in the Integration driver"
           << " (Runge-Kutta) is too many." << G4endl
           << "Integration of the interval was not completed !" << G4endl
           << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
           << " % fraction of it was done.";
   G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
               JustWarning, message);
                                    G4double   h , 
                                    G4double  eps,
                                    G4int     dbg)
 static G4ThreadLocal G4double maxRelError=0.0;
 G4bool isNewMax, prNewMax;
 isNewMax = endPointDist > (1.0 + maxRelError) * h;
 prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
 if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
 if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) 
         && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
 { 
   static G4ThreadLocal G4int noWarnings = 0;
   std::ostringstream message;
   if( (noWarnings ++ < 10) || (dbg>2) )
   {
     message << "The integration produced an end-point which " << G4endl
             << "is further from the start-point than the curve length."
             << G4endl;
   }
   message << "  Distance of endpoints = " << endPointDist
           << ", curve length = " << h << G4endl
           << "  Difference (curveLen-endpDist)= " << (h - endPointDist)
           << ", relative = " << (h-endPointDist) / h 
           << ", epsilon =  " << eps;
   G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
               JustWarning, message);
 }
                            const G4double dydx[],
                                  G4double& x,         // InOut
                                  G4double htry,
                                  G4double eps_rel_max,
                                  G4double& hdid,      // Out
                                  G4double& hnext )    // Out
 G4double errmax_sq;
 G4double h, htemp, xnew ;
 G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
 h = htry ; // Set stepsize to the initial trial value
 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
 G4double errpos_sq=0.0;    // square of displacement error
 G4double errvel_sq=0.0;    // square of momentum vector difference
 G4double errspin_sq=0.0;   // square of spin vector difference
 G4int iter;
 static G4ThreadLocal G4int tot_no_trials=0; 
 const G4int max_trials=100; 
 G4ThreeVector Spin(y[9],y[10],y[11]);
 G4double   spin_mag2 =Spin.mag2() ;
 G4bool     hasSpin= (spin_mag2 > 0.0); 
 for (iter=0; iter<max_trials ;iter++)
 {
   tot_no_trials++;
   pIntStepper-> Stepper(y,dydx,h,ytemp,yerr); 
   //            *******
   G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep); 
   G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos); 
   // Evaluate accuracy
   //
   errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
   errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
   // Accuracy for momentum
   G4double magvel_sq=  sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
   G4double sumerr_sq =  sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ; 
   if( magvel_sq > 0.0 ) { 
      errvel_sq = sumerr_sq / magvel_sq; 
   }else{
      G4cerr << "** G4MagIntegrationDriver: found case of zero momentum." 
             << " iteration=  " << iter << " h= " << h << G4endl; 
      errvel_sq = sumerr_sq; 
   }
   errvel_sq *= inv_eps_vel_sq;
   errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
   if( hasSpin )
   { 
     // Accuracy for spin
     errspin_sq =  ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
                   /  spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
     errspin_sq *= inv_eps_vel_sq;
     errmax_sq = std::max( errmax_sq, errspin_sq ); 
   }
   if ( errmax_sq <= 1.0 )  { break; } // Step succeeded. 
   // Step failed; compute the size of retrial Step.
   htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
   if (htemp >= 0.1*h)  { h = htemp; }  // Truncation error too large,
   else  { h = 0.1*h; }                 // reduce stepsize, but no more
                                        // than a factor of 10
   xnew = x + h;
   if(xnew == x)
   {
     G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
            << "  Stepsize underflow in Stepper " << G4endl ;
     G4cerr << "  Step's start x=" << x << " and end x= " << xnew 
            << " are equal !! " << G4endl
            <<"  Due to step-size= " << h 
            << " . Note that input step was " << htry << G4endl;
     break;
   }
 }
 // Sum of squares of position error // and momentum dir (underestimated)
 fSumH_lg += h; 
 fDyerrPos_lgTot += errpos_sq;
 fDyerrVel_lgTot += errvel_sq * h * h; 
 // Compute size of next Step
 if (errmax_sq > errcon*errcon)
 { 
   hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
 }
 else
 {
   hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
 }
 x += (hdid = h);
 for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
 return;
                           G4FieldTrack& y_posvel,         // INOUT
                           const G4double     dydx[],  
                                 G4double     hstep,       // In
                                 G4double&    dchord_step,
                                 G4double&    dyerr_pos_sq,
                                 G4double&    dyerr_mom_rel_sq )  
 G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
             FatalException, "Not yet implemented."); 
 // Use the parameters of this method, to please compiler
 dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0]; 
 dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
 return true;
                           G4FieldTrack& y_posvel,         // INOUT
                           const G4double     dydx[],  
                                 G4double     hstep,       // In
                                 G4double&    dchord_step,
                                 G4double&    dyerr )
 G4double dyerr_pos_sq, dyerr_mom_rel_sq;  
 G4double yerr_vec[G4FieldTrack::ncompSVEC],
          yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC]; 
 G4double s_start;
 G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
 static G4ThreadLocal G4int no_call=0; 
 no_call ++; 
 // Move data into array
 y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel 
 s_start = y_posvel.GetCurveLength();
 // Do an Integration Step
 pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ; 
 //            *******
 // Estimate curve-chord distance
 dchord_step= pIntStepper-> DistChord();
 //                         *********
 // Put back the values.  yarrout ==> y_posvel
 y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
 y_posvel.SetCurveLength( s_start + hstep );
 if(fVerboseLevel>2)
 {
   G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
   PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep,  1); 
 }
 // A single measure of the error   
 //      TO-DO :  account for  energy,  spin, ... ? 
 vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
 inv_vel_mag_sq = 1.0 / vel_mag_sq; 
 dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
 dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
 dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
 // Calculate also the change in the momentum squared also ???
 // G4double veloc_square = y_posvel.GetVelocity().mag2();
 // ...
 // The following step cannot be done here because "eps" is not known.
 dyerr_len = std::sqrt( dyerr_len_sq ); 
 dyerr_len_sq /= eps ;
 // Look at the velocity deviation ?
 //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
 // Set suggested new step
 hstep= ComputeNewStepSize( dyerr_len, hstep);
 if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
 {
   dyerr = std::sqrt(dyerr_pos_sq);
 }
 else
 {
   // Scale it to the current step size - for now
   dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
 }
 return true;
                                 G4double     yarrin[],    // In
                           const G4double     dydx[],  
                                 G4double     hstep,       // In
                                 G4double     yarrout[],
                                 G4double&    dchord_step,
                                 G4double&    dyerr )      // In length
 G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
             FatalException, "Not yet implemented.");
 dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
 yarrout[0]= yarrin[0];
                         G4double  errMaxNorm,    // max error  (normalised)
                         G4double  hstepCurrent)  // current step size
 G4double hnew;
 // Compute size of next Step for a failed step
 if(errMaxNorm > 1.0 )
 {
   // Step failed; compute the size of retrial Step.
   hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
 } else if(errMaxNorm > 0.0 ) {
   // Compute size of next Step for a successful step
   hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
 } else {
   // if error estimate is zero (possible) or negative (dubious)
   hnew = max_stepping_increase * hstepCurrent; 
 }
 return hnew;
                         G4double  errMaxNorm,    // max error  (normalised)
                         G4double  hstepCurrent)  // current step size
 G4double hnew;
 // Compute size of next Step for a failed step
 if (errMaxNorm > 1.0 )
 {
   // Step failed; compute the size of retrial Step.
   hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
 
   if (hnew < max_stepping_decrease*hstepCurrent)
   {
     hnew = max_stepping_decrease*hstepCurrent ;
                        // reduce stepsize, but no more
                        // than this factor (value= 1/10)
   }
 }
 else
 {
   // Compute size of next Step for a successful step
   if (errMaxNorm > errcon)
    { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
   else  // No more than a factor of 5 increase
    { hnew = max_stepping_increase * hstepCurrent; }
 }
 return hnew;
                                  G4double          xstart,
                                  const G4double*   CurrentArr, 
                                  G4double          xcurrent,
                                  G4double          requestStep, 
                                  G4int             subStepNo)
 // Potentially add as arguments:  
 //                                 <dydx>           - as Initial Force
 //                                 stepTaken(hdid)  - last step taken
 //                                 nextStep (hnext) - proposal for size
  G4FieldTrack  StartFT(G4ThreeVector(0,0,0),
                G4ThreeVector(0,0,0), 0., 0., 0., 0. );
  G4FieldTrack  CurrentFT (StartFT);
  StartFT.LoadFromArray( StartArr, fNoIntegrationVariables); 
  StartFT.SetCurveLength( xstart);
  CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables); 
  CurrentFT.SetCurveLength( xcurrent );
  PrintStatus(StartFT, CurrentFT, requestStep, subStepNo ); 
                 const G4FieldTrack&  StartFT,
                 const G4FieldTrack&  CurrentFT, 
                 G4double             requestStep, 
                 G4int                subStepNo)
   G4int verboseLevel= fVerboseLevel;
   static G4ThreadLocal G4int noPrecision= 5;
   G4int oldPrec= G4cout.precision(noPrecision);
   // G4cout.setf(ios_base::fixed,ios_base::floatfield);
   const G4ThreeVector StartPosition=       StartFT.GetPosition();
   const G4ThreeVector StartUnitVelocity=   StartFT.GetMomentumDir();
   const G4ThreeVector CurrentPosition=     CurrentFT.GetPosition();
   const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
   G4double  DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
   G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
   G4double subStepSize = step_len;
    
   if( (subStepNo <= 1) || (verboseLevel > 3) )
   {
      subStepNo = - subStepNo;        // To allow printing banner
      G4cout << std::setw( 6)  << " " << std::setw( 25)
             << " G4MagInt_Driver: Current Position  and  Direction" << " "
             << G4endl; 
      G4cout << std::setw( 5) << "Step#" << " "
             << std::setw( 7) << "s-curve" << " "
             << std::setw( 9) << "X(mm)" << " "
             << std::setw( 9) << "Y(mm)" << " "  
             << std::setw( 9) << "Z(mm)" << " "
             << std::setw( 8) << " N_x " << " "
             << std::setw( 8) << " N_y " << " "
             << std::setw( 8) << " N_z " << " "
             << std::setw( 8) << " N^2-1 " << " "
             << std::setw(10) << " N(0).N " << " "
             << std::setw( 7) << "KinEner " << " "
             << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
             << std::setw(12) << "Step-len" << " " 
             << std::setw(12) << "Step-len" << " " 
             << std::setw( 9) << "ReqStep" << " "  
             << G4endl;
   }
   if( (subStepNo <= 0) )
   {
     PrintStat_Aux( StartFT,  requestStep, 0., 
                      0,        0.0,         1.0);
     //*************
   }
   if( verboseLevel <= 3 )
   {
     G4cout.precision(noPrecision);
     PrintStat_Aux( CurrentFT, requestStep, step_len, 
                    subStepNo, subStepSize, DotStartCurrentVeloc );
     //*************
   }
   else // if( verboseLevel > 3 )
   {
      //  Multi-line output
      
      // G4cout << "Current  Position is " << CurrentPosition << G4endl 
      //    << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
      // G4cout << "Step taken was " << step_len  
      //    << " out of PhysicalStep= " <<  requestStep << G4endl;
      // G4cout << "Final safety is: " << safety << G4endl;
      // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
      //        << G4endl << G4endl; 
   }
   G4cout.precision(oldPrec);
                 const G4FieldTrack&  aFieldTrack,
                 G4double             requestStep, 
                 G4double             step_len,
                 G4int                subStepNo,
                 G4double             subStepSize,
                 G4double             dotVeloc_StartCurr)
   const G4ThreeVector Position=      aFieldTrack.GetPosition();
   const G4ThreeVector UnitVelocity=  aFieldTrack.GetMomentumDir();

   if( subStepNo >= 0)
   {
      G4cout << std::setw( 5) << subStepNo << " ";
   }
   else
   {
      G4cout << std::setw( 5) << "Start" << " ";
   }
   G4double curveLen= aFieldTrack.GetCurveLength();
   G4cout << std::setw( 7) << curveLen;
   G4cout << std::setw( 9) << Position.x() << " "
          << std::setw( 9) << Position.y() << " "
          << std::setw( 9) << Position.z() << " "
          << std::setw( 8) << UnitVelocity.x() << " "
          << std::setw( 8) << UnitVelocity.y() << " "
          << std::setw( 8) << UnitVelocity.z() << " ";
   G4int oldprec= G4cout.precision(3);
   G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
   G4cout.precision(6);
   G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
   G4cout.precision(oldprec);
   G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
   G4cout << std::setw(12) << step_len << " ";
   static G4ThreadLocal G4double oldCurveLength= 0.0;
   static G4ThreadLocal G4double oldSubStepLength= 0.0;
   static G4ThreadLocal G4int oldSubStepNo= -1;
   G4double subStep_len=0.0;
   if( curveLen > oldCurveLength )
   {
     subStep_len= curveLen - oldCurveLength;
   }
   else if (subStepNo == oldSubStepNo)
   {
     subStep_len= oldSubStepLength;
   }
   oldCurveLength= curveLen;
   oldSubStepLength= subStep_len;
   G4cout << std::setw(12) << subStep_len << " "; 
   G4cout << std::setw(12) << subStepSize << " "; 
   if( requestStep != -1.0 )
   {
     G4cout << std::setw( 9) << requestStep << " ";
   }
   else
   {
      G4cout << std::setw( 9) << " InitialStep " << " ";
   }
   G4cout << G4endl;
 G4int noPrecBig= 6;
 G4int oldPrec= G4cout.precision(noPrecBig);
 G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
 G4cout << "G4MagInt_Driver: Number of Steps: "
        << " Total= " <<  fNoTotalSteps
        << " Bad= "   <<  fNoBadSteps 
        << " Small= " <<  fNoSmallSteps 
        << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
        << G4endl;
 G4cout << "MID dyerr: " 
        << " maximum= " << fDyerr_max 
        << " Sum small= " << fDyerrPos_smTot 
        << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
        << " vel= " << std::sqrt( fDyerrVel_lgTot )
        << " Total h-distance: small= " << fSumH_sm 
        << " large= " << fSumH_lg
        << G4endl;
 G4int noPrecSmall=4; 
 // Single line precis of statistics ... optional
 G4cout.precision(noPrecSmall);
 G4cout << "MIDnums: " << fMinimumStep
        << "   " << fNoTotalSteps 
        << "  "  <<  fNoSmallSteps
        << "  "  << fNoSmallSteps-fNoInitialSmallSteps
        << "  "  << fNoBadSteps         
        << "   " << fDyerr_max
        << "   " << fDyerr_mx2 
        << "   " << fDyerrPos_smTot 
        << "   " << fSumH_sm
        << "   " << fDyerrPos_lgTot
        << "   " << fDyerrVel_lgTot
        << "   " << fSumH_lg
        << G4endl;
G4cout.precision(oldPrec);

 if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
 {
   fSmallestFraction= newFraction;
 }
 else
 { 
   G4cerr << "Warning: SmallestFraction not changed. " << G4endl
          << "  Proposed value was " << newFraction << G4endl
          << "  Value must be between 1.e-8 and 1.e-16" << G4endl;
 }

Return to problem 1696

Lines 552-581 Link Here

(-)src/G4MagIntegratorDriver.cc (-54 / +4 lines)
552	static G4ThreadLocal G4int tot_no_trials=0;	552	static G4ThreadLocal G4int tot_no_trials=0;
553	const G4int max_trials=100;	553	const G4int max_trials=100;
554		554
555	const G4double max_reduction = 0.1 ; // Largest_reduction_factor for step size
556	const G4double local_Pow_shrink = 0.5*GetPshrnk();
557
558	// We choose to use a linear approximation of the Taylor series
559	// ( 1 + x ) ^ p = 1 + p * x - 0.5 * p * (1-p) x^2 + O( x^3 )
560	// if the second order term is smaller than
561	const G4double acc_limit= 0.01;
562
563	// This limit is reached when
564	// 0.5 * p * (1-p) x^2 = acc_limit
565	// ie
566	// p * (1-p) x^2 = acc_limit * 2.0
567
568	// NEW constants -- here for now only
569	// The following value is the value of 'x' at which this limit is reached:
570	const G4double limit_fac = 0.5 * local_Pow_shrink * (1.0 - local_Pow_shrink) ;
571	// linear_limit_fac
572	// To avoid the division here by limit_fac, we check below that
573	// limit_fac * x^2 < acc_limt
574
575	const G4double err_limit = std::pow( max_reduction / GetSafety(), 2/GetPshrnk() );
576	// This is constant when the value of the 'shrink' power is constant
577	// TODO: These should become data members ( this one is expensive to calculate! )
578
579	G4ThreeVector Spin(y[9],y[10],y[11]);	555	G4ThreeVector Spin(y[9],y[10],y[11]);
580	G4double spin_mag2 =Spin.mag2() ;	556	G4double spin_mag2 =Spin.mag2() ;
581	G4bool hasSpin= (spin_mag2 > 0.0);	557	G4bool hasSpin= (spin_mag2 > 0.0);
Lines 616-654 Link Here
616	}	592	}
617		593
618	if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.	594	if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
		595
619	// Step failed; compute the size of retrial Step.	596	// Step failed; compute the size of retrial Step.
		597	htemp = GetSafety()h std::pow( errmax_sq, 0.5*GetPshrnk() );
620		598
621	#if 1 // NEW_METHOD	599	if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
622	// New method of calculation - avoids power if step is very small	600	else { h = 0.1*h; } // reduce stepsize, but no more
623	G4double hnew = 0.0;
624	G4double x = errmax_sq - 1.0;
625	if (errmax_sq >= err_limit) {
626	// if( limit_fac * x * x < acc_limt ) {
627	// hnew = 1.0 + local_Pow_shrink * x;
628	// }else{
629	hnew = GetSafety()h std::pow( errmax_sq, local_Pow_shrink); // 0.5*GetPshrnk() );
630	// }
631	}
632	else { hnew = max_reduction*h; } // Truncation error was too large,
633	// reduce stepsize, but no more than x10
634	#endif
635
636	#if 1 // OLD_METHOD
637	// Old method of calculation - always used power
638	htemp = GetSafety()h std::pow( errmax_sq, local_Pow_shrink); // 0.5*GetPshrnk() );
639
640	if (htemp >= max_reduction*h) { h = htemp; } // Truncation error too large,
641	else { h = max_reduction*h; } // reduce stepsize, but no more
642	// than a factor of 10	601	// than a factor of 10
643	#if 1 // NEW_METHOD
644	assert( std::fabs(hnew - h) < 0.00001 * std::fabs(hnew+h) );
645	#endif
646	#endif
647
648	#if 1 // NEW_METHOD
649	h= hnew;
650	#endif
651
652	xnew = x + h;	602	xnew = x + h;
653	if(xnew == x)	603	if(xnew == x)
654	{	604	{

Lines 16-25 Link Here

(-)History (-3 / +9 lines)
16	----------------------------------------------------------	16	----------------------------------------------------------
17	* Reverse chronological order (last date on top), please *	17	* Reverse chronological order (last date on top), please *
18	----------------------------------------------------------	18	----------------------------------------------------------
19	--------------------------	19	January 15, 2015 - J.Apostolakis (geomnav-V10-01-03)
20	Added protection in G4ReplicaNavigation::DistanceToOutRad() for potential	20	---------------------------------
		21	- G4PathFinder: Pass to equation of motion relevant properties of current
		22	particle: latest charge state (electric, magnetic dipole moment),
		23	PDG spin value, momentum and rest mass of current particle
		24	Part of corrections needed for Bugzilla issue 1696.
		25
		26	- Increased maximum number of Navigators from 8 to 16.
		27
21	rare cases of negative value to sqrt() in equation calculation for rmin/rmax	28	rare cases of negative value to sqrt() in equation calculation for rmin/rmax
22	intersection.	29	intersection.
23	Improved warning printout in G4Navigator::GetGlobalExitNormal().
24		30