Ignore:
Timestamp:
Feb 16, 2009, 10:14:30 AM (16 years ago)
Author:
garnier
Message:

en test de gl2ps. Problemes de libraries

Location:
trunk/source/geometry/magneticfield/src
Files:
38 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/magneticfield/src/G4CashKarpRKF45.cc

    r850 r921  
    2626//
    2727// $Id: G4CashKarpRKF45.cc,v 1.15 2008/01/11 18:11:44 japost Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth

  • trunk/source/geometry/magneticfield/src/G4ChordFinder.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4ChordFinder.cc,v 1.49 2008/07/15 14:02:06 japost Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4ChordFinder.cc,v 1.51 2008/10/29 14:17:42 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    31 // 25.02.97 John Apostolakis,  design and implimentation
    32 // 05.03.97 V. Grichine , style modification
     31// 25.02.97 - John Apostolakis - Design and implementation
    3332// -------------------------------------------------------------------
     33
     34#include <iomanip>
    3435
    3536#include "G4ChordFinder.hh"
    … …  
    3839#include "G4ClassicalRK4.hh"
    3940
    40 #include <iomanip>
    4141
    4242// ..........................................................................
    … …  
    5454{
    5555  // Simple constructor which does not create equation, ..
    56       // fDeltaChord= fDefaultDeltaChord;
     56
    5757  fIntgrDriver= pIntegrationDriver;
    5858  fAllocatedStepper= false;
    … …  
    6262    // check the values and set the other parameters
    6363}
     64
    6465
    6566// ..........................................................................
    … …  
    8081  //  Construct the Chord Finder
    8182  //  by creating in inverse order the  Driver, the Stepper and EqRhs ...
     83
    8284  G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
    8385  fEquation = pEquation;                           
    … …  
    103105                                     pItsStepper->GetNumberOfVariables() );
    104106}
     107
     108
     109// ......................................................................
     110
     111G4ChordFinder::~G4ChordFinder()
     112{
     113  delete   fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
     114  if( fAllocatedStepper)
     115  {
     116     delete fDriversStepper;
     117  }
     118  delete   fIntgrDriver;
     119
     120  if( fStatsVerbose ) { PrintStatistics(); }
     121}
     122
    105123
    106124// ......................................................................
    … …  
    113131  if( fractNext == -1.0 )   fractNext = 0.98;  // 0.9;
    114132
    115   // fFirstFraction  = 0.999;  // Orig 0.999 A safe value,  range: ~ 0.95 - 0.999
    116   // fMultipleRadius = 15.0;   // For later use, range: ~  2 - 20
    117 
    118   if( fStatsVerbose ) {
     133  // fFirstFraction  = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
     134  // fMultipleRadius = 15.0;  // For later use, range: ~  2 - 20
     135
     136  if( fStatsVerbose )
     137  {
    119138    G4cout << " ChordFnd> Trying to set fractions: "
    120            << " first " << fFirstFraction
    121            << " last " <<  fractLast
    122            << " next " <<  fractNext
    123            << " and multiple " << fMultipleRadius
    124            << G4endl;
     139           << " first " << fFirstFraction
     140           << " last " <<  fractLast
     141           << " next " <<  fractNext
     142           << " and multiple " << fMultipleRadius
     143           << G4endl;
    125144  }
    126145
    127146  if( (fractLast > 0.0) && (fractLast <=1.0) )
    128     { fFractionLast= fractLast; }
    129   else
     147  {
     148    fFractionLast= fractLast;
     149  }
     150  else
     151  {
     152    G4cerr << "G4ChordFinder::SetFractions_Last_Next: Invalid "
     153           << " fraction Last = " << fractLast
     154           << " must be  0 <  fractionLast <= 1 " << G4endl;
     155  }
     156  if( (fractNext > 0.0) && (fractNext <1.0) )
     157  {
     158    fFractionNextEstimate = fractNext;
     159  }
     160  else
     161  {
    130162    G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
    131            << " fraction Last = " << fractLast
    132            << " must be  0 <  fractionLast <= 1 " << G4endl;
    133   if( (fractNext > 0.0) && (fractNext <1.0) )
    134     { fFractionNextEstimate = fractNext; }
    135   else
    136     G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
    137            << " fraction Next = " << fractNext
    138            << " must be  0 <  fractionNext < 1 " << G4endl;
    139 }
     163           << " fraction Next = " << fractNext
     164           << " must be  0 <  fractionNext < 1 " << G4endl;
     165  }
     166}
     167
    140168
    141169// ......................................................................
    142170
    143 G4ChordFinder::~G4ChordFinder()
    144 {
    145   delete   fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
    146   if( fAllocatedStepper)
     171G4double
     172G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
     173                                    G4double      stepMax,
     174                                    G4double      epsStep,
     175                                    const G4ThreeVector latestSafetyOrigin,
     176                                    G4double       latestSafetyRadius )
     177{
     178  G4double stepPossible;
     179  G4double dyErr;
     180  G4FieldTrack yEnd( yCurrent);
     181  G4double  startCurveLen= yCurrent.GetCurveLength();
     182  G4double nextStep;
     183  //            *************
     184  stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
     185                              &nextStep, latestSafetyOrigin, latestSafetyRadius
     186                             );
     187  //            *************
     188
     189  G4bool good_advance;
     190
     191  if ( dyErr < epsStep * stepPossible )
     192  {
     193     // Accept this accuracy.
     194
     195     yCurrent = yEnd;
     196     good_advance = true;
     197  }
     198  else
     199  { 
     200     // Advance more accurately to "end of chord"
     201     //                           ***************
     202     good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible,
     203                                                  epsStep, nextStep);
     204     if ( ! good_advance )
     205     {
     206       // In this case the driver could not do the full distance
     207       stepPossible= yCurrent.GetCurveLength()-startCurveLen;
     208     }
     209  }
     210  return stepPossible;
     211}
     212
     213
     214// ............................................................................
     215
     216G4double
     217G4ChordFinder::FindNextChord( const  G4FieldTrack& yStart,
     218                                     G4double     stepMax,
     219                                     G4FieldTrack&   yEnd, // Endpoint
     220                                     G4double&   dyErrPos, // Error of endpoint
     221                                     G4double    epsStep,
     222                                     G4double*  pStepForAccuracy,
     223                              const  G4ThreeVector, //  latestSafetyOrigin,
     224                                     G4double       //  latestSafetyRadius
     225                                        )
     226{
     227  // Returns Length of Step taken
     228
     229  G4FieldTrack yCurrent=  yStart; 
     230  G4double    stepTrial, stepForAccuracy;
     231  G4double    dydx[G4FieldTrack::ncompSVEC];
     232
     233  //  1.)  Try to "leap" to end of interval
     234  //  2.)  Evaluate if resulting chord gives d_chord that is good enough.
     235  // 2a.)  If d_chord is not good enough, find one that is.
     236 
     237  G4bool    validEndPoint= false;
     238  G4double  dChordStep, lastStepLength; //  stepOfLastGoodChord;
     239
     240  fIntgrDriver-> GetDerivatives( yCurrent, dydx );
     241
     242  G4int     noTrials=0;
     243  const G4double safetyFactor= fFirstFraction; //  0.975 or 0.99 ? was 0.999
     244
     245  stepTrial = std::min( stepMax, safetyFactor*fLastStepEstimate_Unconstrained );
     246
     247  G4double newStepEst_Uncons= 0.0;
     248  do
    147249  {
    148      delete fDriversStepper;
    149   }                                //  fIntgrDriver->pIntStepper;}
    150   delete   fIntgrDriver;
    151 
    152   if( fStatsVerbose ) { PrintStatistics();  }
    153 }
     250     G4double stepForChord; 
     251     yCurrent = yStart;    // Always start from initial point
     252   
     253     //            ************
     254     fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
     255                                 dChordStep, dyErrPos);
     256     //            ************
     257     
     258     //  We check whether the criterion is met here.
     259     validEndPoint = AcceptableMissDist(dChordStep);
     260
     261     lastStepLength = stepTrial;
     262
     263     // This method estimates to step size for a good chord.
     264     stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
     265
     266     if( ! validEndPoint )
     267     {
     268        if( stepTrial<=0.0 )
     269        {
     270          stepTrial = stepForChord;
     271        }
     272        else if (stepForChord <= stepTrial)
     273        {
     274          // Reduce by a fraction, possibly up to 20%
     275          stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
     276        }
     277        else
     278        {
     279          stepTrial *= 0.1;
     280        }
     281     }
     282     noTrials++;
     283  }
     284  while( ! validEndPoint );   // End of do-while  RKD
     285
     286  if( newStepEst_Uncons > 0.0  )
     287  {
     288     fLastStepEstimate_Unconstrained= newStepEst_Uncons;
     289  }
     290
     291  AccumulateStatistics( noTrials );
     292
     293  if( pStepForAccuracy )
     294  {
     295     // Calculate the step size required for accuracy, if it is needed
     296     //
     297     G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
     298     if( dyErr_relative > 1.0 )
     299     {
     300        stepForAccuracy = fIntgrDriver->ComputeNewStepSize( dyErr_relative,
     301                                                            lastStepLength );
     302     }
     303     else
     304     {
     305        stepForAccuracy = 0.0;   // Convention to show step was ok
     306     }
     307     *pStepForAccuracy = stepForAccuracy;
     308  }
     309
     310#ifdef  TEST_CHORD_PRINT
     311  static int dbg=0;
     312  if( dbg )
     313  {
     314    G4cout << "ChordF/FindNextChord:  NoTrials= " << noTrials
     315           << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
     316  }
     317#endif
     318  yEnd=  yCurrent; 
     319  return stepTrial;
     320}
     321
     322
     323// ...........................................................................
     324
     325G4double G4ChordFinder::NewStep(G4double  stepTrialOld,
     326                                G4double  dChordStep, // Curr. dchord achieved
     327                                G4double& stepEstimate_Unconstrained ) 
     328{
     329  // Is called to estimate the next step size, even for successful steps,
     330  // in order to predict an accurate 'chord-sensitive' first step
     331  // which is likely to assist in more performant 'stepping'.
     332
     333  G4double stepTrial;
     334  static G4double lastStepTrial = 1.,  lastDchordStep= 1.;
     335
     336#if 1
     337
     338  if (dChordStep > 0.0)
     339  {
     340    stepEstimate_Unconstrained =
     341                 stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
     342    stepTrial =  fFractionNextEstimate * stepEstimate_Unconstrained;
     343  }
     344  else
     345  {
     346    // Should not update the Unconstrained Step estimate: incorrect!
     347    stepTrial =  stepTrialOld * 2.;
     348  }
     349
     350  if( stepTrial <= 0.001 * stepTrialOld)
     351  {
     352     if ( dChordStep > 1000.0 * fDeltaChord )
     353     {
     354        stepTrial= stepTrialOld * 0.03;   
     355     }
     356     else
     357     {
     358        if ( dChordStep > 100. * fDeltaChord )
     359        {
     360          stepTrial= stepTrialOld * 0.1;   
     361        }
     362        else   // Try halving the length until dChordStep OK
     363        {
     364          stepTrial= stepTrialOld * 0.5;   
     365        }
     366     }
     367  }
     368  else if (stepTrial > 1000.0 * stepTrialOld)
     369  {
     370     stepTrial= 1000.0 * stepTrialOld;
     371  }
     372
     373  if( stepTrial == 0.0 )
     374  {
     375     stepTrial= 0.000001;
     376  }
     377
     378  lastStepTrial = stepTrialOld;
     379  lastDchordStep= dChordStep;
     380
     381#else
     382
     383  if ( dChordStep > 1000. * fDeltaChord )
     384  {
     385        stepTrial= stepTrialOld * 0.03;   
     386  }
     387  else
     388  {
     389     if ( dChordStep > 100. * fDeltaChord )
     390     {
     391        stepTrial= stepTrialOld * 0.1;   
     392     }
     393     else  // Keep halving the length until dChordStep OK
     394     {
     395        stepTrial= stepTrialOld * 0.5;   
     396     }
     397  }
     398
     399#endif
     400
     401  // A more sophisticated chord-finder could figure out a better
     402  // stepTrial, from dChordStep and the required d_geometry
     403  //   e.g.
     404  //      Calculate R, r_helix (eg at orig point)
     405  //      if( stepTrial < 2 pi  R )
     406  //          stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
     407  //      else   
     408  //          ??
     409
     410  return stepTrial;
     411}
     412
     413
     414// ...........................................................................
     415
     416G4FieldTrack
     417G4ChordFinder::ApproxCurvePointS( const G4FieldTrack&  CurveA_PointVelocity,
     418                                  const G4FieldTrack&  CurveB_PointVelocity,
     419                                  const G4FieldTrack&  ApproxCurveV,
     420                                  const G4ThreeVector& CurrentE_Point,
     421                                  const G4ThreeVector& CurrentF_Point,
     422                                  const G4ThreeVector& PointG,
     423                                       G4bool first, G4double eps_step)
     424{
     425  // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
     426  // Use Brent Algorithm (or InvParabolic) when possible.
     427  // Given a starting curve point A (CurveA_PointVelocity), curve point B
     428  // (CurveB_PointVelocity), a point E which is (generally) not on the curve
     429  // and  a point F which is on the curve (first approximation), find new
     430  // point S on the curve closer to point E.
     431  // While advancing towards S utilise 'eps_step' as a measure of the
     432  // relative accuracy of each Step.
     433
     434  G4FieldTrack EndPoint( CurveA_PointVelocity);
     435  G4ThreeVector Point_A=CurveA_PointVelocity.GetPosition();
     436  G4ThreeVector Point_B=CurveB_PointVelocity.GetPosition();
     437  G4double xa,xb,xc,ya,yb,yc;
     438 
     439  // InverseParabolic. AF Intersects (First Part of Curve)
     440
     441  if(first)
     442  {
     443      xa=0.;
     444      ya=(PointG-Point_A).mag();
     445      xb=(Point_A-CurrentF_Point).mag();
     446      yb=-(PointG-CurrentF_Point).mag();
     447      xc=(Point_A-Point_B).mag();
     448      yc=-(CurrentE_Point-Point_B).mag();
     449  }   
     450  else
     451  {
     452     xa=0.;
     453     ya=(Point_A-PointG).mag();
     454     xb=(Point_B-Point_A).mag();
     455     yb=-(PointG-Point_B).mag();
     456     xc=-(Point_A-CurrentF_Point).mag();
     457     yc=-(Point_A-CurrentE_Point).mag();
     458   
     459  }
     460  const G4double tolerance= 1.e-12;
     461  if(ya<=tolerance||std::abs(yc)<=tolerance)
     462  {
     463    ; // What to do for the moment: return the same point as at start
     464      // then PropagatorInField will take care
     465  }
     466  else
     467  {
     468    G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
     469    G4double curve;
     470    if(first)
     471    {
     472      curve=std::abs(EndPoint.GetCurveLength()
     473                    -ApproxCurveV.GetCurveLength());
     474    }
     475    else
     476    {
     477      curve=std::abs(EndPoint.GetCurveLength()
     478                    -CurveB_PointVelocity.GetCurveLength());
     479    }
     480     
     481    if(test_step<=0)    { test_step=0.1*xb; }
     482    if(test_step>=xb)   { test_step=0.5*xb; }
     483    if(test_step>=curve){ test_step=0.5*curve; }
     484
     485    if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
     486    {                          // G4VIntersectionLocator
     487      test_step=0.5*curve;
     488    }
     489
     490    G4bool goodAdvance;
     491    goodAdvance = fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
     492     
     493#ifdef G4DEBUG_FIELD
     494    // Printing Brent and Linear Approximation
     495    //
     496    G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
     497           << test_step << "  EndPoint = " << EndPoint << G4endl;
     498
     499    //  Test Track
     500    //
     501    G4FieldTrack TestTrack( CurveA_PointVelocity);
     502    TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
     503                                   CurveB_PointVelocity,
     504                                   CurrentE_Point, eps_step );
     505    G4cout.precision(14);
     506    G4cout << "G4ChordFinder::BrentApprox = " << EndPoint  << G4endl;
     507    G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
     508#endif
     509  }
     510  return EndPoint;
     511}
     512
     513
     514// ...........................................................................
     515
     516G4FieldTrack G4ChordFinder::
     517ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
     518                   const G4FieldTrack& CurveB_PointVelocity,
     519                   const G4ThreeVector& CurrentE_Point,
     520                         G4double eps_step)
     521{
     522  // If r=|AE|/|AB|, and s=true path lenght (AB)
     523  // return the point that is r*s along the curve!
     524 
     525  G4FieldTrack   Current_PointVelocity = CurveA_PointVelocity;
     526
     527  G4ThreeVector  CurveA_Point= CurveA_PointVelocity.GetPosition();
     528  G4ThreeVector  CurveB_Point= CurveB_PointVelocity.GetPosition();
     529
     530  G4ThreeVector  ChordAB_Vector= CurveB_Point   - CurveA_Point;
     531  G4ThreeVector  ChordAE_Vector= CurrentE_Point - CurveA_Point;
     532
     533  G4double       ABdist= ChordAB_Vector.mag();
     534  G4double  curve_length;  //  A curve length  of AB
     535  G4double  AE_fraction;
     536 
     537  curve_length= CurveB_PointVelocity.GetCurveLength()
     538              - CurveA_PointVelocity.GetCurveLength(); 
     539 
     540  G4double  integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
     541  if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
     542  {
     543#ifdef G4DEBUG_FIELD
     544    G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
     545           << G4endl
     546           << " The two points are further apart than the curve length "
     547           << G4endl
     548           << " Dist = "         << ABdist
     549           << " curve length = " << curve_length
     550           << " relativeDiff = " << (curve_length-ABdist)/ABdist
     551           << G4endl;
     552    if( curve_length < ABdist * (1. - 10*eps_step) )
     553    {
     554      G4cerr << " ERROR: the size of the above difference"
     555             << " exceeds allowed limits.  Aborting." << G4endl;
     556      G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
     557                  FatalException, "Unphysical curve length.");
     558    }
     559#endif
     560    // Take default corrective action: adjust the maximum curve length.
     561    // NOTE: this case only happens for relatively straight paths.
     562    // curve_length = ABdist;
     563  }
     564
     565  G4double  new_st_length;
     566
     567  if ( ABdist > 0.0 )
     568  {
     569     AE_fraction = ChordAE_Vector.mag() / ABdist;
     570  }
     571  else
     572  {
     573     AE_fraction = 0.5;                         // Guess .. ?;
     574#ifdef G4DEBUG_FIELD
     575     G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
     576            << " A and B are the same point!" << G4endl
     577            << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
     578            << G4endl;
     579#endif
     580  }
     581 
     582  if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
     583  {
     584#ifdef G4DEBUG_FIELD
     585    G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
     586           << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
     587           << "   AE_fraction = " <<  AE_fraction << G4endl
     588           << "   Chord AE length = " << ChordAE_Vector.mag() << G4endl
     589           << "   Chord AB length = " << ABdist << G4endl << G4endl;
     590    G4cerr << " OK if this condition occurs after a recalculation of 'B'"
     591           << G4endl << " Otherwise it is an error. " << G4endl ;
     592#endif
     593     // This course can now result if B has been re-evaluated,
     594     // without E being recomputed (1 July 99).
     595     // In this case this is not a "real error" - but it is undesired
     596     // and we cope with it by a default corrective action ...
     597     //
     598     AE_fraction = 0.5;                         // Default value
     599  }
     600
     601  new_st_length= AE_fraction * curve_length;
     602
     603  G4bool good_advance;
     604  if ( AE_fraction > 0.0 )
     605  {
     606     good_advance = fIntgrDriver->AccurateAdvance(Current_PointVelocity,
     607                                                  new_st_length, eps_step );
     608     //
     609     // In this case it does not matter if it cannot advance the full distance
     610  }
     611
     612  // If there was a memory of the step_length actually required at the start
     613  // of the integration Step, this could be re-used ...
     614
     615  G4cout.precision(14);
     616
     617  return Current_PointVelocity;
     618}
     619
     620
     621// ......................................................................
    154622
    155623void
    … …  
    157625{
    158626  // Print Statistics
     627
    159628  G4cout << "G4ChordFinder statistics report: " << G4endl;
    160629  G4cout
    … …  
    171640}
    172641
    173 // ......................................................................
    174 
    175 G4double
    176 G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
    177                                     G4double      stepMax,
    178                                     G4double      epsStep,
    179                                     const G4ThreeVector latestSafetyOrigin,
    180                                     G4double       latestSafetyRadius
    181                                     )
    182 {
    183   G4double stepPossible;
    184   G4double dyErr;
    185   G4FieldTrack yEnd( yCurrent);
    186   G4double  startCurveLen= yCurrent.GetCurveLength();
    187   G4double nextStep;
    188   //            *************
    189   stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep, &nextStep
    190                               , latestSafetyOrigin, latestSafetyRadius
    191                              );
    192   //            *************
    193   //  G4cout<<"Exit Find Next Chord Err= "<<dyErr<<" eps=  "<<epsStep<<"stepPos=  "<<stepPossible<<G4endl;
    194   G4bool good_advance;
    195 
    196   if ( dyErr < epsStep * stepPossible )
    197  
    198       {   //G4cout<<"err comparison = "<<dyErr<<" eps=  "<<epsStep<<G4endl;
    199      // Accept this accuracy.
    200      yCurrent = yEnd;
    201      good_advance = true;
    202   }
    203   else
    204   { 
    205     // G4cout<<"Entering Accurate Advance"<<G4endl;
    206      // Advance more accurately to "end of chord"
    207      //                           ***************
    208      good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep, nextStep);
    209      //                           ***************
    210      if ( ! good_advance ){
    211        // In this case the driver could not do the full distance
    212        stepPossible= yCurrent.GetCurveLength()-startCurveLen;
    213      }
    214   }
    215 
    216 #ifdef G4DEBUG_FIELD
    217   //G4cout << "Exiting FindNextChord Limited with:" << G4endl
    218   //       << "   yCurrent: " << yCurrent<< G4endl;
    219 #endif
    220 
    221   return stepPossible;
    222 }
    223 
    224 // #define  TEST_CHORD_PRINT  1
    225 
    226 // ............................................................................
    227 
    228 G4double
    229 G4ChordFinder::FindNextChord( const  G4FieldTrack& yStart,
    230                                      G4double     stepMax,
    231                                      G4FieldTrack&   yEnd, // Endpoint
    232                                      G4double&   dyErrPos, // Error of endpoint
    233                                      G4double    epsStep,
    234                                      G4double*  pStepForAccuracy,
    235                               const  G4ThreeVector, //  latestSafetyOrigin,
    236                                      G4double       //  latestSafetyRadius
    237                                         )
    238   // Returns Length of Step taken
    239 {
    240   //G4cout<<"Inter Find Next Chord with step="<<stepMax<<G4endl; 
    241   // G4int       stepRKnumber=0;
    242   G4FieldTrack yCurrent=  yStart; 
    243   G4double    stepTrial, stepForAccuracy;
    244   G4double    dydx[G4FieldTrack::ncompSVEC];
    245 
    246   //  1.)  Try to "leap" to end of interval
    247   //  2.)  Evaluate if resulting chord gives d_chord that is good enough.
    248   //     2a.)  If d_chord is not good enough, find one that is.
    249  
    250   G4bool    validEndPoint= false;
    251   G4double  dChordStep, lastStepLength; //  stepOfLastGoodChord;
    252 
    253   fIntgrDriver-> GetDerivatives( yCurrent, dydx )  ;
    254   //for (G4int i=0;i<G4FieldTrack::ncompSVEC;i++){
    255   //  dydx[i]=0.;
    256   //}
    257   G4int     noTrials=0;
    258   const G4double safetyFactor= fFirstFraction; //  0.975 or 0.99 ? was 0.999
    259 
    260   stepTrial = std::min( stepMax,
    261                           safetyFactor * fLastStepEstimate_Unconstrained );
    262 
    263   G4double newStepEst_Uncons= 0.0;
    264   do
    265   {
    266      G4double stepForChord; 
    267      yCurrent = yStart;    // Always start from initial point
    268    
    269      //            ************
    270      fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
    271                                  dChordStep, dyErrPos);
    272      //            ************
    273 
    274      // G4cout<<"AfterQuickAdv step="<<stepTrial<<"  dC="<<dChordStep<<" yErr="<<dyErrPos<<G4endl;
    275      
    276       //  We check whether the criterion is met here.
    277      validEndPoint = AcceptableMissDist(dChordStep);
    278      // if(validEndPoint){G4cout<<"validEndPoint"<<fDeltaChord<<G4endl;}
    279      // else{G4cout<<"No__validEndPoint"<<G4endl;}
    280 
    281 
    282      //&& (dyErrPos < eps) ;
    283      //   validEndPoint = AcceptableMissDist(dChordStep) && (dyErrPos < epsStep) ;
    284 
    285      lastStepLength = stepTrial;
    286 
    287      // This method estimates to step size for a good chord.
    288      stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
    289 
    290      if( ! validEndPoint ) {
    291         if( stepTrial<=0.0 )
    292           stepTrial = stepForChord;
    293         else if (stepForChord <= stepTrial)
    294           // Reduce by a fraction, possibly up to 20%
    295           stepTrial = std::min( stepForChord,
    296                                 fFractionLast * stepTrial);
    297         else
    298           stepTrial *= 0.1;
    299 
    300         // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
    301      }
    302      // #ifdef  TEST_CHORD_PRINT
    303      // TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
    304      // #endif
    305 
    306      noTrials++;
    307   }
    308   while( ! validEndPoint );   // End of do-while  RKD
    309 
    310   if( newStepEst_Uncons > 0.0  ){
    311     fLastStepEstimate_Unconstrained= newStepEst_Uncons;
    312   }
    313 
    314   AccumulateStatistics( noTrials );
    315 
    316   // stepOfLastGoodChord = stepTrial;
    317 
    318   if( pStepForAccuracy ){
    319      // Calculate the step size required for accuracy, if it is needed
    320      G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
    321      if( dyErr_relative > 1.0 ) {
    322         stepForAccuracy =
    323            fIntgrDriver->ComputeNewStepSize( dyErr_relative,
    324                                              lastStepLength );
    325      }else{
    326         stepForAccuracy = 0.0;   // Convention to show step was ok
    327      }
    328      *pStepForAccuracy = stepForAccuracy;
    329   }
    330 
    331 #ifdef  TEST_CHORD_PRINT
    332   static int dbg=0;
    333   if( dbg )
    334     G4cout << "ChordF/FindNextChord:  NoTrials= " << noTrials
    335            << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
    336 #endif
    337   //G4cout << "ChordF/FindNextChord:  NoTrials= " << noTrials
    338   //         << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
    339   yEnd=  yCurrent; 
    340   return stepTrial;
    341 }
    342 
    343 // ----------------------------------------------------------------------------
    344 #if 0         
    345 //   #ifdef G4VERBOSE
    346 if( dbg ) {
    347    G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
    348    G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
    349 }
    350 #endif
    351 // ----------------------------------------------------------------------------
    352642
    353643// ...........................................................................
    354644
    355 G4double G4ChordFinder::NewStep(G4double  stepTrialOld,
    356                                 G4double  dChordStep, // Curr. dchord achieved
    357                                 G4double& stepEstimate_Unconstrained ) 
    358 //
    359 // Is called to estimate the next step size, even for successful steps,
    360 // in order to predict an accurate 'chord-sensitive' first step
    361 // which is likely to assist in more performant 'stepping'.
    362 //
    363        
    364 {
    365   G4double stepTrial;
    366   static G4double lastStepTrial = 1.,  lastDchordStep= 1.;
    367 
    368 #if 1
    369   // const G4double  threshold = 1.21, multiplier = 0.9;
    370   //  0.9 < 1 / std::sqrt(1.21)
    371 
    372   if (dChordStep > 0.0)
    373   {
    374     stepEstimate_Unconstrained = stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
    375     // stepTrial =  0.98 * stepEstimate_Unconstrained;
    376     stepTrial =  fFractionNextEstimate * stepEstimate_Unconstrained;
    377   }
    378   else
    379   {
    380     // Should not update the Unconstrained Step estimate: incorrect!
    381     stepTrial =  stepTrialOld * 2.;
    382   }
    383 
    384   // if ( dChordStep < threshold * fDeltaChord ){
    385   //    stepTrial= stepTrialOld *  multiplier;   
    386   // }
    387   if( stepTrial <= 0.001 * stepTrialOld)
    388   {
    389      if ( dChordStep > 1000.0 * fDeltaChord ){
    390         stepTrial= stepTrialOld * 0.03;   
    391      }else{
    392         if ( dChordStep > 100. * fDeltaChord ){
    393           stepTrial= stepTrialOld * 0.1;   
    394         }else{
    395     // Try halving the length until dChordStep OK
    396           stepTrial= stepTrialOld * 0.5;   
    397         }
    398      }
    399   }else if (stepTrial > 1000.0 * stepTrialOld)
    400   {
    401      stepTrial= 1000.0 * stepTrialOld;
    402   }
    403 
    404   if( stepTrial == 0.0 ){
    405      stepTrial= 0.000001;
    406   }
    407 
    408   lastStepTrial = stepTrialOld;
    409   lastDchordStep= dChordStep;
    410 #else
    411   if ( dChordStep > 1000. * fDeltaChord ){
    412         stepTrial= stepTrialOld * 0.03;   
    413   }else{
    414      if ( dChordStep > 100. * fDeltaChord ){
    415         stepTrial= stepTrialOld * 0.1;   
    416      }else{
    417         // Keep halving the length until dChordStep OK
    418         stepTrial= stepTrialOld * 0.5;   
    419      }
    420   }
    421 #endif
    422 
    423   // A more sophisticated chord-finder could figure out a better
    424   //   stepTrial, from dChordStep and the required d_geometry
    425   //   eg
    426   //      Calculate R, r_helix (eg at orig point)
    427   //      if( stepTrial < 2 pi  R )
    428   //          stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
    429   //      else   
    430   //          ??
    431 
    432   return stepTrial;
    433 }
    434 
    435 //  ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
    436 //  Use Brent Algorithm(or InvParabolic) when it possible.
    437 //  Given a starting curve point A (CurveA_PointVelocity), 
    438 //  curve point B (CurveB_PointVelocity), a point E which is (generally)
    439 //  not on the curve  and  a point F which is on the curve(first approximation)
    440 //  From this information find new point S on the curve closer to point E.
    441 //  While advancing towards S utilise eps_step
    442 //  as a measure of the relative accuracy of each Step.
    443 G4FieldTrack
    444 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
    445                                   const G4FieldTrack& CurveB_PointVelocity,
    446                                   const G4ThreeVector& CurrentE_Point,
    447                                   const G4ThreeVector& CurrentF_Point,
    448                                   const G4ThreeVector& PointG,
    449                                        G4bool first, G4double eps_step)
    450 {
    451 
    452  
    453   G4FieldTrack EndPoint( CurveA_PointVelocity);
    454   G4ThreeVector Point_A=CurveA_PointVelocity.GetPosition();
    455   G4ThreeVector Point_B=CurveB_PointVelocity.GetPosition();
    456   G4double xa,xb,xc,ya,yb,yc;
    457  
    458 //InverseParabolic
    459 //AF Intersects (First Part of Curve)
    460   if(first){
    461       xa=0.;
    462       ya=(PointG-Point_A).mag();
    463       xb=(Point_A-CurrentF_Point).mag();
    464       yb=-(PointG-CurrentF_Point).mag();
    465       xc=(Point_A-Point_B).mag();
    466       yc=-(CurrentE_Point-Point_B).mag();
    467   }   
    468   else{
    469      xa=0.;
    470      ya=(Point_A-PointG).mag();
    471      xb=(Point_B-Point_A).mag();
    472      yb=-(PointG-Point_B).mag();
    473      xc=-(Point_A-CurrentF_Point).mag();
    474      yc=-(Point_A-CurrentE_Point).mag();
    475    
    476   }
    477   const G4double tolerance= 1.e-12;
    478   if(ya<=tolerance||std::abs(yc)<=tolerance){
    479     ; //What to do for the moment return the same point as in begin
    480      //Then PropagatorInField will take care
    481    }
    482    else{
    483          
    484       G4double test_step  =InvParabolic(xa,ya,xb,yb,xc,yc);
    485       G4double curve=std::abs(EndPoint.GetCurveLength()-CurveB_PointVelocity.GetCurveLength());
    486       G4double dist= (EndPoint.GetPosition()-Point_B).mag();
    487       if(test_step<=0) { test_step=0.1*xb;}
    488       if(test_step>=xb){ test_step=0.5*xb;}
    489        
    490 
    491       if(curve*(1.+eps_step)<dist){
    492         test_step=0.5*dist;
    493        }
    494 
    495        G4bool goodAdvance;
    496        goodAdvance=
    497              fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
    498             //            ***************
    499      
    500        #ifdef G4DEBUG_FIELD
    501         G4cout<<"G4ChordFinder:: test-step ShF="<<test_step<<"  EndPoint="<<EndPoint<<G4endl;
    502       //    Test Track
    503        G4FieldTrack TestTrack( CurveA_PointVelocity);
    504        TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
    505                                                   CurveB_PointVelocity,
    506                                                   CurrentE_Point,
    507                                                   eps_step );
    508        G4cout.precision(14);
    509        G4cout<<"G4ChordFinder:: BrentApprox="<<EndPoint<<G4endl;
    510        G4cout<<"G4ChordFinder::LinearApprox="<<TestTrack<<G4endl;
    511        #endif
    512   }
    513 return EndPoint;
    514 } 
    515 
    516 G4FieldTrack
    517 G4ChordFinder::ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
    518                                   const G4FieldTrack& CurveB_PointVelocity,
    519                                   const G4ThreeVector& CurrentE_Point,
    520                                         G4double eps_step)
    521 {
    522   // 1st implementation:
    523   //    if r=|AE|/|AB|, and s=true path lenght (AB)
    524   //    return the point that is r*s along the curve!
    525   /////////////////////////////
    526   //
    527   //2st implementation : Inverse Parabolic Extrapolation by D.C.Williams
    528   //
    529   //    Uses InvParabolic (xa,ya,xb,yb,xc,yc)
    530 
    531   G4FieldTrack    Current_PointVelocity = CurveA_PointVelocity;
    532 
    533   G4ThreeVector  CurveA_Point= CurveA_PointVelocity.GetPosition();
    534   G4ThreeVector  CurveB_Point= CurveB_PointVelocity.GetPosition();
    535 
    536   G4ThreeVector  ChordAB_Vector= CurveB_Point   - CurveA_Point;
    537   G4ThreeVector  ChordAE_Vector= CurrentE_Point - CurveA_Point;
    538 
    539   G4double       ABdist= ChordAB_Vector.mag();
    540   G4double  curve_length;  //  A curve length  of AB
    541   G4double  AE_fraction;
    542  
    543   curve_length= CurveB_PointVelocity.GetCurveLength()
    544               - CurveA_PointVelocity.GetCurveLength(); 
    545 
    546   // const
    547   G4double  integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
    548   if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
    549 #ifdef G4DEBUG_FIELD
    550     G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
    551            << G4endl
    552            << " The two points are further apart than the curve length "
    553            << G4endl
    554            << " Dist = "         << ABdist
    555            << " curve length = " << curve_length
    556            << " relativeDiff = " << (curve_length-ABdist)/ABdist
    557            << G4endl;
    558     if( curve_length < ABdist * (1. - 10*eps_step) ) {
    559       G4cerr << " ERROR: the size of the above difference"
    560              << " exceeds allowed limits.  Aborting." << G4endl;
    561       G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
    562                   FatalException, "Unphysical curve length.");
    563     }
    564 #endif
    565     // Take default corrective action:
    566     //    -->  adjust the maximum curve length.
    567     //  NOTE: this case only happens for relatively straight paths.
    568     curve_length = ABdist;
    569   }
    570 
    571   G4double  new_st_length;
    572 
    573   if ( ABdist > 0.0 ){
    574      AE_fraction = ChordAE_Vector.mag() / ABdist;
    575   }else{
    576      AE_fraction = 0.5;                         // Guess .. ?;
    577 #ifdef G4DEBUG_FIELD
    578      G4cout << "Warning in G4ChordFinder::ApproxCurvePoint:"
    579             << " A and B are the same point!" << G4endl
    580             << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
    581             << G4endl;
    582 #endif
    583   }
    584  
    585   if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) ){
    586 #ifdef G4DEBUG_FIELD
    587     G4cerr << " G4ChordFinder::ApproxCurvePointV - Warning:"
    588            << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
    589            << "   AE_fraction = " <<  AE_fraction << G4endl
    590            << "   Chord AE length = " << ChordAE_Vector.mag() << G4endl
    591            << "   Chord AB length = " << ABdist << G4endl << G4endl;
    592     G4cerr << " OK if this condition occurs after a recalculation of 'B'"
    593            << G4endl << " Otherwise it is an error. " << G4endl ;
    594 #endif
    595      // This course can now result if B has been re-evaluated,
    596      //   without E being recomputed   (1 July 99)
    597      //  In this case this is not a "real error" - but it undesired
    598      //   and we cope with it by a default corrective action ...
    599      AE_fraction = 0.5;                         // Default value
    600   }
    601 
    602   new_st_length= AE_fraction * curve_length;
    603 
    604   G4bool good_advance;
    605   if ( AE_fraction > 0.0 ) {
    606      good_advance =
    607       fIntgrDriver->AccurateAdvance(Current_PointVelocity,
    608                                     new_st_length,
    609                                     eps_step ); // Relative accuracy
    610      // In this case it does not matter if it cannot advance the full distance
    611   }
    612 
    613   // If there was a memory of the step_length actually require at the start
    614   // of the integration Step, this could be re-used ...
    615    G4cout.precision(14);
    616      
    617    //     G4cout<<"G4ChordFinder::LinearApprox="<<Current_PointVelocity<<G4endl;
    618   return Current_PointVelocity;
    619 }
    620 
    621 void
    622 G4ChordFinder::TestChordPrint( G4int    noTrials,
    623                                G4int    lastStepTrial,
    624                                G4double dChordStep,
    625                                G4double nextStepTrial )
     645void G4ChordFinder::TestChordPrint( G4int    noTrials,
     646                                    G4int    lastStepTrial,
     647                                    G4double dChordStep,
     648                                    G4double nextStepTrial )
    626649{
    627650     G4int oldprec= G4cout.precision(5);
    628651     G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
    629652            << " this_step= "       << std::setw(10) << lastStepTrial;
    630      if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
    631              G4cout.precision(8);
    632      }else{  G4cout.precision(6); }
    633      G4cout << " dChordStep=  "     << std::setw(12) << dChordStep;
     653     if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
     654     {
     655       G4cout.precision(8);
     656     }
     657     else
     658     {
     659       G4cout.precision(6);
     660     }
     661     G4cout << " dChordStep=  " << std::setw(12) << dChordStep;
    634662     if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
    635663     else                           { G4cout << " d-"; }

  • trunk/source/geometry/magneticfield/src/G4ChordFinderSaf.cc

    r831 r921  
    116116
    117117G4double
    118 G4ChordFinderSaf::FindNextChord( const  G4FieldTrack  yStart,
     118G4ChordFinderSaf::FindNextChord( const  G4FieldTrack&  yStart,
    119119                                     G4double     stepMax,
    120120                                     G4FieldTrack&   yEnd, // Endpoint

  • trunk/source/geometry/magneticfield/src/G4ClassicalRK4.cc

    r850 r921  
    2626//
    2727// $Id: G4ClassicalRK4.cc,v 1.12 2006/06/29 18:23:37 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4DELPHIMagField.cc

    r850 r921  
    2626//
    2727// $Id: G4DELPHIMagField.cc,v 1.6 2006/06/29 18:23:39 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929// -------------------------------------------------------------------
    3030

  • trunk/source/geometry/magneticfield/src/G4ElectricField.cc

    r850 r921  
    2626//
    2727// $Id: G4ElectricField.cc,v 1.2 2006/06/29 18:23:42 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4ElectroMagneticField.cc

    r850 r921  
    2626//
    2727// $Id: G4ElectroMagneticField.cc,v 1.3 2006/06/29 18:23:44 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4EqEMFieldWithSpin.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithSpin.cc,v 1.2 2008/04/24 12:33:08 tnikitin Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4EqEMFieldWithSpin.cc,v 1.4 2008/11/21 21:17:03 gum Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    4040
    4141#include "G4EqEMFieldWithSpin.hh"
     42#include "G4ElectroMagneticField.hh"
    4243#include "G4ThreeVector.hh"
    4344#include "globals.hh"
    4445
    4546G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField )
    46       : G4EquationOfMotion( emField ) { anomaly = 1.165923e-3; }
     47      : G4EquationOfMotion( emField )
     48{
     49  anomaly = 0.0011659208;
     50}
     51
     52G4EqEMFieldWithSpin::~G4EqEMFieldWithSpin()
     53{
     54}
    4755
    4856void 
    4957G4EqEMFieldWithSpin::SetChargeMomentumMass(G4double particleCharge, // e+ units
    50                                             G4double MomentumXc,
     58                                            G4double MomentumXc,
    5159                                            G4double particleMass)
    5260{
    … …  
    6169   beta  = MomentumXc/E;
    6270   gamma = E/particleMass;
     71
    6372}
    64 
    65 
    6673
    6774void
    6875G4EqEMFieldWithSpin::EvaluateRhsGivenB(const G4double y[],
    69                                         const G4double Field[],
    70                                               G4double dydx[] ) const
     76                                       const G4double Field[],
     77                                             G4double dydx[] ) const
    7178{
    7279
    … …  
    114121
    115122   G4ThreeVector Spin(y[9],y[10],y[11]);
     123
     124   if (Spin.mag() > 0.) Spin = Spin.unit();
     125
    116126   G4ThreeVector dSpin;
    117127

  • trunk/source/geometry/magneticfield/src/G4EqMagElectricField.cc

    r850 r921  
    2626//
    2727// $Id: G4EqMagElectricField.cc,v 1.14 2008/04/24 12:33:35 tnikitin Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4EquationOfMotion.cc

    r850 r921  
    2626//
    2727// $Id: G4EquationOfMotion.cc,v 1.9 2006/06/29 18:23:48 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4ErrorMag_UsualEqRhs.cc

    r850 r921  
    2626//
    2727// $Id: G4ErrorMag_UsualEqRhs.cc,v 1.1 2007/05/16 12:54:02 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4ExactHelixStepper.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4ExactHelixStepper.cc,v 1.8 2007/12/10 16:29:47 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4ExactHelixStepper.cc,v 1.9 2008/10/29 14:34:35 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//  Helix a-la-Explicity Euler: x_1 = x_0 + helix(h)
    … …  
    112112{
    113113  // Implementation : must check whether h/R >  pi  !!
    114   //   If( h/R <  pi)   DistChord=h/2*tan(Ang_curve/4)
     114  //   If( h/R <  pi)   DistChord=h/2*std::tan(Ang_curve/4)
    115115  //   Else             DistChord=R_helix
    116116  //

  • trunk/source/geometry/magneticfield/src/G4ExplicitEuler.cc

    r850 r921  
    2626//
    2727// $Id: G4ExplicitEuler.cc,v 1.8 2006/06/29 18:23:53 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4FieldManager.cc

    r850 r921  
    2626//
    2727// $Id: G4FieldManager.cc,v 1.15 2007/12/07 15:34:10 japost Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4FieldManagerStore.cc

    r850 r921  
    2626//
    2727// $Id: G4FieldManagerStore.cc,v 1.4 2008/01/17 10:56:23 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// G4FieldManagerStore

  • trunk/source/geometry/magneticfield/src/G4FieldTrack.cc

    r850 r921  
    2626//
    2727// $Id: G4FieldTrack.cc,v 1.14 2007/10/03 15:34:42 japost Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4HarmonicPolMagField.cc

    r850 r921  
    2626//
    2727// $Id: G4HarmonicPolMagField.cc,v 1.6 2006/06/29 18:24:00 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4HelixExplicitEuler.cc

    r850 r921  
    2626//
    2727// $Id: G4HelixExplicitEuler.cc,v 1.8 2007/12/10 16:29:49 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4HelixHeum.cc

    r850 r921  
    2626//
    2727// $Id: G4HelixHeum.cc,v 1.6 2006/06/29 18:24:04 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4HelixImplicitEuler.cc

    r850 r921  
    2626//
    2727// $Id: G4HelixImplicitEuler.cc,v 1.6 2006/06/29 18:24:06 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4HelixSimpleRunge.cc

    r850 r921  
    2626//
    2727// $Id: G4HelixSimpleRunge.cc,v 1.7 2006/06/29 18:24:08 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4ImplicitEuler.cc

    r850 r921  
    2626//
    2727// $Id: G4ImplicitEuler.cc,v 1.9 2006/06/29 18:24:11 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4LineCurrentMagField.cc

    r850 r921  
    2525//
    2626// $Id: G4LineCurrentMagField.cc,v 1.6 2006/06/29 18:24:13 gunter Exp $
    27 // GEANT4 tag $Name: HEAD $
     27// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2828// -------------------------------------------------------------------
    2929

  • trunk/source/geometry/magneticfield/src/G4LineSection.cc

    r850 r921  
    2626//
    2727// $Id: G4LineSection.cc,v 1.10 2006/06/29 18:24:16 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4MagErrorStepper.cc

    r850 r921  
    2626//
    2727// $Id: G4MagErrorStepper.cc,v 1.13 2006/06/29 18:24:18 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4MagHelicalStepper.cc

    r850 r921  
    2626//
    2727// $Id: G4MagHelicalStepper.cc,v 1.23 2007/09/05 12:20:17 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4MagIntegratorDriver.cc

    r850 r921  
    2626//
    2727// $Id: G4MagIntegratorDriver.cc,v 1.49 2007/08/17 12:30:33 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4MagIntegratorStepper.cc

    r850 r921  
    2626//
    2727// $Id: G4MagIntegratorStepper.cc,v 1.11 2006/06/29 18:24:34 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4Mag_EqRhs.cc

    r850 r921  
    2626//
    2727// $Id: G4Mag_EqRhs.cc,v 1.11 2006/06/29 18:24:36 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//  This is the standard right-hand side for equation of motion 

  • trunk/source/geometry/magneticfield/src/G4Mag_SpinEqRhs.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4Mag_SpinEqRhs.cc,v 1.12 2006/06/29 18:24:39 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4Mag_SpinEqRhs.cc,v 1.13 2008/11/21 21:18:26 gum Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// This is the standard right-hand side for equation of motion.
    … …  
    4545  : G4Mag_EqRhs( MagField )
    4646{
    47    anomaly = 1.165923e-3;
     47   anomaly = 0.0011659208;
    4848}
    4949
    … …  
    5353G4Mag_SpinEqRhs::SetChargeMomentumMass(G4double particleCharge, // in e+ units
    5454                                       G4double MomentumXc,
    55                                        G4double mass)
     55                                       G4double particleMass)
    5656{
    5757   //  To set fCof_val
    58    G4Mag_EqRhs::SetChargeMomentumMass(particleCharge, MomentumXc, mass);
     58   G4Mag_EqRhs::SetChargeMomentumMass(particleCharge, MomentumXc, particleMass);
    5959
    60    omegac = 0.105658387*GeV/mass * 2.837374841e-3*(rad/cm/kilogauss);
     60   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6161
    6262   ParticleCharge = particleCharge;
    6363
    64    E = std::sqrt(sqr(MomentumXc)+sqr(mass));
     64   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    6565   beta  = MomentumXc/E;
    66    gamma = E/mass;
     66   gamma = E/particleMass;
    6767
    6868}
    … …  
    9696
    9797   G4ThreeVector Spin(y[9],y[10],y[11]);
     98
     99   if (Spin.mag() > 0.) Spin = Spin.unit();
     100
    98101   G4ThreeVector dSpin;
    99102

  • trunk/source/geometry/magneticfield/src/G4Mag_UsualEqRhs.cc

    r850 r921  
    2626//
    2727// $Id: G4Mag_UsualEqRhs.cc,v 1.12 2006/06/29 18:24:42 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4MagneticField.cc

    r850 r921  
    2626//
    2727// $Id: G4MagneticField.cc,v 1.3 2006/06/29 18:24:44 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4QuadrupoleMagField.cc

    r850 r921  
    2626//
    2727// $Id: G4QuadrupoleMagField.cc,v 1.4 2006/06/29 18:24:46 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4RKG3_Stepper.cc

    r850 r921  
    2626//
    2727// $Id: G4RKG3_Stepper.cc,v 1.15 2007/08/21 10:17:41 tnikitin Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4SimpleHeum.cc

    r850 r921  
    2626//
    2727// $Id: G4SimpleHeum.cc,v 1.8 2006/06/29 18:24:51 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//  Simple Heum:

  • trunk/source/geometry/magneticfield/src/G4SimpleRunge.cc

    r850 r921  
    2626//
    2727// $Id: G4SimpleRunge.cc,v 1.10 2006/06/29 18:24:53 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//  Simple Runge:

  • trunk/source/geometry/magneticfield/src/G4UniformElectricField.cc

    r850 r921  
    2626//
    2727// $Id: G4UniformElectricField.cc,v 1.12 2006/06/29 18:24:56 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4UniformMagField.cc

    r850 r921  
    2626//
    2727// $Id: G4UniformMagField.cc,v 1.11 2006/06/29 18:24:58 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
Note: See TracChangeset for help on using the changeset viewer.