Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

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

Legend:

Unmodified
Added
Removed

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4CashKarpRKF45.cc,v 1.15 2008/01/11 18:11:44 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4CashKarpRKF45.cc,v 1.16 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
    … …  
    4949                                 G4int noIntegrationVariables,
    5050                                 G4bool primary)
    51   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
    52 {
    53   // unsigned int noVariables= std::max(numberOfVariables,8); // For Time .. 7+1
     51  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
     52    fLastStepLength(0.), fAuxStepper(0)
     53{
    5454  const G4int numberOfVariables = noIntegrationVariables;
    5555
    … …  
    6969  fMidVector = new G4double[numberOfVariables];
    7070  fMidError =  new G4double[numberOfVariables];
    71   fAuxStepper = 0;   
    72   if( primary )
    73       fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
    74 
     71  if( primary )
     72  {
     73    fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
     74  }
    7575}
    7676

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ConstRK4.cc,v 1.2 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ConstRK4.cc,v 1.5 2010/09/10 15:51:10 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    3838//////////////////////////////////////////////////////////////////
    3939//
    40 // Constructor sets the number of variables (default = 8)
    41 
    42 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numberOfVariables)
    43   : G4MagErrorStepper(EqRhs, numberOfVariables)
    44 {
    45   if(numberOfVariables !=8 )
     40// Constructor sets the number of *State* variables (default = 8)
     41//   The number of variables integrated is always 6
     42
     43G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numStateVariables)
     44  : G4MagErrorStepper(EqRhs, 6, numStateVariables)
     45{
     46  // const G4int numberOfVariables= 6;
     47  if( numStateVariables < 8 )
    4648  {
     49    G4cerr << "ERROR in G4ConstRK4::G4ConstRK4 "
     50           << "   The number of State variables at least 8 " << G4endl;
     51    G4cerr << "   Instead it is  numStateVariables= " << numStateVariables << G4endl;
    4752    G4Exception("G4ConstRK4::G4ConstRK4()", "InvalidSetup", FatalException,
    48                 "Valid only for number of variables=8. Use another Stepper!");
     53        "Valid only for number of state variables of 8 or more. Use another Stepper!");
    4954  }
    50   else
    51   {
    52     fEq=EqRhs;
    53     yMiddle= new G4double[8];
    54     dydxMid= new G4double[8];
    55     yInitial= new G4double[8];
    56     yOneStep= new G4double[8];
    57 
    58     dydxm = new G4double[8];
    59     dydxt = new G4double[8];
    60     yt    = new G4double[8];
    61     Field[0]=0.;Field[1]=0.;Field[2]=0.;
    62   }
     55
     56  fEq = EqRhs;
     57  yMiddle  = new G4double[8];
     58  dydxMid  = new G4double[8];
     59  yInitial = new G4double[8];
     60  yOneStep = new G4double[8];
     61
     62  dydxm = new G4double[8];
     63  dydxt = new G4double[8];
     64  yt    = new G4double[8];
     65  Field[0]=0.; Field[1]=0.; Field[2]=0.;
    6366}
    6467
    … …  
    150153                           G4double yError [] )
    151154{
    152    const G4int nvar = 8 ;
    153    const G4int maxvar= 8;
     155   const G4int nvar = 6;  // number of variables integrated
     156   const G4int maxvar= GetNumberOfStateVariables();
     157
     158   // Correction for Richardson extrapolation
     159   G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
    154160
    155161   G4int i;
    156 
    157    // Correction for Richardson extrapolation
    158    //
    159    G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
    160162   
    161163   // Saving yInput because yInput and yOutput can be aliases for same array
    162 
    163    for (i=0;i<nvar;i++) { yInitial[i]=yInput[i]; }
    164 
    165    yInitial[7]= yInput[7];  // Copy the time in case...even if not really needed
    166    yMiddle[7] = yInput[7];  // Copy the time from initial value
    167    yOneStep[7] = yInput[7]; // As it contributes to final value of yOutput ?
    168    yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
    169    for (i=nvar;i<maxvar;i++) { yOutput[i]=yInput[i]; }
     164   for (i=0;    i<maxvar; i++) { yInitial[i]= yInput[i]; }
     165 
     166   // Must copy the part of the state *not* integrated to the output
     167   for (i=nvar; i<maxvar; i++) { yOutput[i]=  yInput[i]; }
     168
     169   // yInitial[7]= yInput[7];  //  The time is typically needed
     170   yMiddle[7]  = yInput[7];   // Copy the time from initial value
     171   yOneStep[7] = yInput[7];   // As it contributes to final value of yOutput ?
     172   // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
    170173   yError[7] = 0.0;         
    171174

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithEDM.cc,v 1.2 2009/11/06 22:31:54 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4EqEMFieldWithEDM.cc,v 1.3 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    4343
    4444G4EqEMFieldWithEDM::G4EqEMFieldWithEDM(G4ElectroMagneticField *emField )
    45       : G4EquationOfMotion( emField )
     45      : G4EquationOfMotion( emField ), fElectroMagCof(0.), fMassCof(0.),
     46        omegac(0.), anomaly(0.0011659208), eta(0.), pcharge(0.), E(0.),
     47        gamma(0.), beta(0.)
    4648{
    47   anomaly = 0.0011659208;
    48   eta = 0.;
    4949}
    5050
    … …  
    6363   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6464
    65    ParticleCharge = particleCharge;
     65   pcharge = particleCharge;
    6666
    6767   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    … …  
    145145
    146146   G4ThreeVector dSpin
    147      = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
     147     = pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
    148148                               // from Jackson
    149149                               // -uce*Spin.cross(u.cross(EField)) )
    150150                               // but this form has one less operation
    151                                - uce*(u*(Spin*EField) - EField*(Spin*u))
    152                                +eta/2.*(Spin.cross(EField) - ude*(Spin.cross(u))
    153                                         // +Spin.cross(u.cross(Bfield))
    154                                         + (u*(Spin*BField) - BField*(Spin*u))
    155                                        )
    156                              );
    157                              
     151                       - uce*(u*(Spin*EField) - EField*(Spin*u))
     152                       + eta/2.*(Spin.cross(EField) - ude*(Spin.cross(u))
     153                               // +Spin.cross(u.cross(Bfield))
     154                       + (u*(Spin*BField) - BField*(Spin*u)) ) );
     155     
    158156   dydx[ 9] = dSpin.x();
    159157   dydx[10] = dSpin.y();

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithSpin.cc,v 1.8 2009/11/06 22:31:35 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4EqEMFieldWithSpin.cc,v 1.9 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    3333//  30.08.2007 Chris Gong, Peter Gumplinger
    3434//  14.02.2009 Kevin Lynch
    35 //  06.11.2009 Hiromi Iinuma see:
    36 //  http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/161.html
     35//  06.11.2009 Hiromi Iinuma
    3736//
    3837// -------------------------------------------------------------------
    … …  
    4443
    4544G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField )
    46   : G4EquationOfMotion( emField )
     45  : G4EquationOfMotion( emField ), fElectroMagCof(0.), fMassCof(0.),
     46    omegac(0.), anomaly(0.0011659208), pcharge(0.), E(0.), gamma(0.), beta(0.)
    4747{
    48   anomaly = 0.0011659208;
    4948}
    5049
    … …  
    6362   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6463
    65    ParticleCharge = particleCharge;
     64   pcharge = particleCharge;
    6665
    6766   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    … …  
    135134
    136135   G4ThreeVector dSpin
    137      = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
     136     = pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
    138137                               // from Jackson
    139138                               // -uce*Spin.cross(u.cross(EField)) );
    140139                               // but this form has one less operation
    141                                - uce*(u*(Spin*EField) - EField*(Spin*u)) );
     140                      - uce*(u*(Spin*EField) - EField*(Spin*u)) );
    142141
    143142   dydx[ 9] = dSpin.x();

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ExactHelixStepper.cc,v 1.9 2008/10/29 14:34:35 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ExactHelixStepper.cc,v 1.11 2010/07/21 13:46:01 tnikitin Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//  Helix a-la-Explicity Euler: x_1 = x_0 + helix(h)
    … …  
    4343#include "G4LineSection.hh"
    4444
    45 
    4645G4ExactHelixStepper::G4ExactHelixStepper(G4Mag_EqRhs *EqRhs)
    4746  : G4MagHelicalStepper(EqRhs),
    48     fBfieldValue(DBL_MAX, DBL_MAX, DBL_MAX), yInitialEHS(DBL_MAX), yFinalEHS(-DBL_MAX)
    49    
     47    fBfieldValue(DBL_MAX, DBL_MAX, DBL_MAX),
     48    fPtrMagEqOfMot(EqRhs)
    5049{
    51    const G4int nvar = 6 ;
    52    G4int i;
    53    for(i=0;i<nvar;i++)  {
    54      fYInSav[i]= DBL_MAX;
    55    }
    56     fPtrMagEqOfMot=EqRhs;
     50  ;
    5751}
    5852
    … …  
    6155void
    6256G4ExactHelixStepper::Stepper( const G4double yInput[],
    63                               const G4double*,
    64                                     G4double hstep,
    65                                     G4double yOut[],
    66                                     G4double yErr[]      )
     57                              const G4double*,
     58                                    G4double hstep,
     59                                    G4double yOut[],
     60                                    G4double yErr[]      )
    6761{ 
    68    const G4int nvar = 6 ;
     62   const G4int nvar = 6;
    6963
    7064   G4int i;
    71    // G4double      yTemp[7], yIn[7] ;
    7265   G4ThreeVector Bfld_value;
    7366
    74    for(i=0;i<nvar;i++)  {
    75       // yIn[i]=     yInput[i];
    76       fYInSav[i]= yInput[i];
    77    }
    78 
    79    MagFieldEvaluate(yInput, Bfld_value) ;       
    80      
    81    // DumbStepper(yIn, Bfld_value, hstep, yTemp);
     67   MagFieldEvaluate(yInput, Bfld_value);       
    8268   AdvanceHelix(yInput, Bfld_value, hstep, yOut);
    8369
    84    // We are assuming a constant field: helix is exact.
    85    for(i=0;i<nvar;i++) {
    86      yErr[i] = 0.0 ;
    87    }
     70  // We are assuming a constant field: helix is exact
     71  //
     72  for(i=0;i<nvar;i++)
     73  {
     74    yErr[i] = 0.0 ;
     75  }
    8876
    89     yInitialEHS = G4ThreeVector( yInput[0],   yInput[1],   yInput[2]);
    90     yFinalEHS   = G4ThreeVector( yOut[0],  yOut[1],  yOut[2]);
    91     fBfieldValue=Bfld_value;
    92    
     77  fBfieldValue=Bfld_value;
    9378}
    9479
    9580void
    9681G4ExactHelixStepper::DumbStepper( const G4double  yIn[],
    97                                    G4ThreeVector   Bfld,
    98                                    G4double  h,
    99                                    G4double  yOut[])
     82                                        G4ThreeVector   Bfld,
     83                                        G4double  h,
     84                                        G4double  yOut[])
    10085{
    10186  // Assuming a constant field: solution is a helix
     87
    10288  AdvanceHelix(yIn, Bfld, h, yOut);
    10389
    104   G4Exception("G4ExactHelixStepper::DumbStepper should not be called.",
    105               "EHS:NoDumbStepper", FatalException, "Stepper must do all the work." );
     90  G4Exception("G4ExactHelixStepper::DumbStepper",
     91              "EHS:NoDumbStepper", FatalException,
     92              "Should not be called. Stepper must do all the work." );
    10693} 
    10794
    … …  
    10996// ---------------------------------------------------------------------------
    11097
    111 G4double G4ExactHelixStepper::DistChord()   const
     98G4double
     99G4ExactHelixStepper::DistChord() const
    112100{
    113101  // Implementation : must check whether h/R >  pi  !!
    114102  //   If( h/R <  pi)   DistChord=h/2*std::tan(Ang_curve/4)
    115103  //   Else             DistChord=R_helix
    116   //
     104
    117105  G4double distChord;
    118106  G4double Ang_curve=GetAngCurve();
    119107
    120      
    121          if(Ang_curve<=pi){
    122            distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
    123          }
    124          else
    125          if(Ang_curve<twopi){
    126            distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
    127          }
    128          else{
    129           distChord=2.*GetRadHelix(); 
    130          }
    131 
    132    
     108  if (Ang_curve<=pi)
     109  {
     110    distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
     111  }
     112  else if(Ang_curve<twopi)
     113  {
     114    distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
     115  }
     116  else
     117  {
     118    distChord=2.*GetRadHelix(); 
     119  }
    133120
    134121  return distChord;
    135  
    136122}   
    137123

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4FieldTrack.cc,v 1.14 2007/10/03 15:34:42 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4FieldTrack.cc,v 1.15 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3636     const G4double *SixV = SixVec.SixVector;
    3737     os << " ( ";
    38      os << " X= " << SixV[0] << " " << SixV[1] << " " << SixV[2] << " ";  // Position
    39      os << " V= " << SixV[3] << " " << SixV[4] << " " << SixV[5] << " ";  // Momentum
    40      os << " v2= " << G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag();     // mom magnitude
     38     os << " X= " << SixV[0] << " " << SixV[1] << " "
     39                  << SixV[2] << " ";  // Position
     40     os << " V= " << SixV[3] << " " << SixV[4] << " "
     41                  << SixV[5] << " ";  // Momentum
     42     os << " v2= "
     43        << G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag(); // mom magnitude
    4144     os << " mdm= " << SixVec.fMomentumDir.mag();
    4245     os << " l= " << SixVec.GetCurveLength();
    … …  
    5760   fRestMass_c2(restMass_c2),
    5861   fLabTimeOfFlight(LaboratoryTimeOfFlight),
    59    // fProperTimeOfFlight(0.0),
     62   fProperTimeOfFlight(0.),
    6063   // fMomentumDir(pMomentumDirection),
    6164   fChargeState(  charge, magnetic_dipole_moment )
    … …  
    6366  G4double momentum  = std::sqrt(kineticEnergy*kineticEnergy
    6467                            +2.0*restMass_c2*kineticEnergy);
    65 
    6668  G4ThreeVector pMomentum= momentum * pMomentumDirection;
    6769  SetCurvePnt( pPosition, pMomentum, curve_length );
    68   // Sets momentum direction as well.
     70    // Sets momentum direction as well.
    6971
    70   // Set the momentum direction again - keeping value from argument exactly
    7172  fMomentumDir=pMomentumDirection;
     73    // Set the momentum direction again - keeping value from argument exactly
    7274
    7375  InitialiseSpin( Spin );
    74 
    75   // fpChargeState = new G4ChargeState( charge, magnetic_dipole_moment );
    7676}
    7777
    … …  
    107107  else         Spin= *pSpin;
    108108  InitialiseSpin( Spin );
    109 
    110   // fpChargeState = new G4ChargeState( DBL_MAX );  //  charge not yet set !!
    111109}
    112110
    113111G4FieldTrack::G4FieldTrack( char )                  //  Nothing is set !!
    114   : fRestMass_c2(0.0), fLabTimeOfFlight(0.0),
    115    fChargeState( DBL_MAX ) //  charge not set
     112  : fKineticEnergy(0.), fRestMass_c2(0.), fLabTimeOfFlight(0.),
     113    fProperTimeOfFlight(0.), fChargeState( DBL_MAX )
    116114{
    117115  G4ThreeVector Zero(0.0, 0.0, 0.0);
    118116  SetCurvePnt( Zero, Zero, 0.0 );
    119117  InitialiseSpin( Zero );
    120   // fpChargeState = new G4ChargeState( DBL_MAX );   
    121118}
    122119

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagHelicalStepper.cc,v 1.23 2007/09/05 12:20:17 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagHelicalStepper.cc,v 1.24 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    4545
    4646G4MagHelicalStepper::G4MagHelicalStepper(G4Mag_EqRhs *EqRhs)
    47    : G4MagIntegratorStepper(EqRhs, 6)  // integrate over 6 variables only !!
     47   : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
    4848                                       // position & velocity
    49 {
    50   fPtrMagEqOfMot = EqRhs;
     49     fPtrMagEqOfMot(EqRhs), fAngCurve(0.), frCurve(0.), frHelix(0.)
     50{
    5151}
    5252

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagIntegratorDriver.cc,v 1.53 2009/11/05 22:31:43 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagIntegratorDriver.cc,v 1.57 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    7878    fNoInitialSmallSteps(0), fDyerr_max(0.0), fDyerr_mx2(0.0),
    7979    fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
    80     fSumH_sm(0.0),   fSumH_lg(0.0),
     80    fSumH_sm(0.0), fSumH_lg(0.0),
    8181    fStatisticsVerboseLevel(statisticsVerbose)
    8282{ 
    … …  
    208208
    209209#ifdef G4DEBUG_FIELD
     210    G4double xSubStepStart= x;
    210211    for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
    211212    yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
    … …  
    227228      lastStepSucceeded= (hdid == h);   
    228229#ifdef G4DEBUG_FIELD
    229       if (dbg>2)
    230       {
    231         PrintStatus( ySubStepStart, xSubStart, y, x, h,  nstp); // Only
     230      if (dbg>2) {
     231        PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
    232232      }
    233233#endif
    … …  
    291291    {
    292292      if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
     293      G4cout << "MagIntDrv: " ;
    293294      G4cout << "hdid="  << std::setw(12) << hdid  << " "
    294              << "hnext=" << std::setw(12) << hnext << " " << G4endl;
     295             << "hnext=" << std::setw(12) << hnext << " "
     296             << "hstep=" << std::setw(12) << hstep << " (requested) "
     297             << G4endl;
    295298      PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
    296299    }
    … …  
    369372        lastStep = true;
    370373#ifdef G4DEBUG_FIELD
    371         if (dbg)
     374        if (dbg>2)
    372375        {
     376          int prec= G4cout.precision(12);
    373377          G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
    374378                 << G4endl
    375379                 << "  Integration step 'h' became "
    376                  << h << " due to roundoff " << G4endl
     380                 << h << " due to roundoff. " << G4endl
     381                 << " Calculated as difference of x2= "<< x2 << " and x=" << x
    377382                 << "  Forcing termination of advance." << G4endl;
     383          G4cout.precision(prec);
    378384        }         
    379385#endif
    … …  
    581587                 /  ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
    582588      errspin_sq *= inv_eps_vel_sq;
    583     }
     589      errmax_sq = std::max( errmax_sq, errspin_sq );
     590   }
    584591
    585592    if ( errmax_sq <= 1.0 )  { break; } // Step succeeded.

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_EqRhs.cc,v 1.11 2006/06/29 18:24:36 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_EqRhs.cc,v 1.12 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//  This is the standard right-hand side for equation of motion 
    … …  
    4848//
    4949G4Mag_EqRhs::G4Mag_EqRhs( G4MagneticField *magField )
    50    : G4EquationOfMotion(magField)
     50   : G4EquationOfMotion(magField), fCof_val(0.)
    5151{
    5252}

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_SpinEqRhs.cc,v 1.15 2009/03/25 15:29:02 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_SpinEqRhs.cc,v 1.16 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// This is the standard right-hand side for equation of motion.
    … …  
    4343
    4444G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
    45   : G4Mag_EqRhs( MagField )
     45  : G4Mag_EqRhs( MagField ), omegac(0.), anomaly(0.0011659208),
     46    pcharge(0.), E(0.), gamma(0.), beta(0.)
    4647{
    47    anomaly = 0.0011659208;
    4848}
    4949
    … …  
    6262   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6363
    64    ParticleCharge = particleCharge;
     64   pcharge = particleCharge;
    6565
    6666   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    … …  
    101101   G4ThreeVector dSpin;
    102102
    103    dSpin = ParticleCharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
     103   dSpin = pcharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
    104104
    105105   dydx[ 9] = dSpin.x();

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_UsualEqRhs.cc,v 1.12 2006/06/29 18:24:42 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_UsualEqRhs.cc,v 1.13 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    4040
    4141G4Mag_UsualEqRhs::G4Mag_UsualEqRhs( G4MagneticField* MagField )
    42   : G4Mag_EqRhs( MagField ) {}
     42  : G4Mag_EqRhs( MagField ), fInvCurrentMomentumXc(0.) {}
    4343
    4444G4Mag_UsualEqRhs::~G4Mag_UsualEqRhs() {}

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4NystromRK4.cc,v 1.6 2009/12/16 17:49:57 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4NystromRK4.cc,v 1.9 2010/09/10 15:42:09 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// History:
    … …  
    4949    m_cachedMom( false )
    5050{
    51   m_fldPosition[0]  = m_iPoint[1] = m_fPoint[2] = 9.9999999e+99 ;
    52   m_fldPosition[1]  = m_iPoint[1] = m_fPoint[2] = 9.9999999e+99 ;
    53   m_fldPosition[2]  = m_iPoint[1] = m_fPoint[2] = 9.9999999e+99 ;
     51  m_fldPosition[0]  = m_iPoint[0] = m_fPoint[0] = m_mPoint[0] = 9.9999999e+99 ;
     52  m_fldPosition[1]  = m_iPoint[1] = m_fPoint[1] = m_mPoint[1] = 9.9999999e+99 ;
     53  m_fldPosition[2]  = m_iPoint[2] = m_fPoint[2] = m_mPoint[2] = 9.9999999e+99 ;
    5454  m_fldPosition[3]  = -9.9999999e+99;
     55  m_lastField[0] = m_lastField[1] = m_lastField[2] = 0.0;
    5556
    5657  m_magdistance2 = distanceConstField*distanceConstField;
    … …  
    196197G4NystromRK4::ComputeRightHandSide(const G4double P[],G4double dPdS[])
    197198{
    198   getField(P);
     199  G4double P4vec[4]= { P[0], P[1], P[2], P[7] }; // Time is P[7]
     200  getField(P4vec);
    199201  m_mom   = std::sqrt(P[3]*P[3]+P[4]*P[4]+P[5]*P[5])     ;
    200202  m_imom  = 1./m_mom                                ;

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4RKG3_Stepper.cc,v 1.15 2007/08/21 10:17:41 tnikitin Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4RKG3_Stepper.cc,v 1.17 2010/07/23 14:13:49 tnikitin Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3535
    3636G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs *EqRhs)
    37   : G4MagIntegratorStepper(EqRhs,6)
    38 {
    39  
    40   fPtrMagEqOfMot=EqRhs;
    41 
     37  : G4MagIntegratorStepper(EqRhs,6), hStep(0.), fPtrMagEqOfMot(EqRhs)
     38{
    4239}
    4340
    … …  
    4643}
    4744
    48 void G4RKG3_Stepper::Stepper(  const G4double  yInput[7],
    49                                const G4double dydx[7],
     45void G4RKG3_Stepper::Stepper(  const G4double yInput[8],
     46                               const G4double dydx[6],
    5047                                     G4double Step,
    51                                      G4double yOut[7],
     48                                     G4double yOut[8],
    5249                                     G4double yErr[])
    5350{
    … …  
    5754   G4double  by15 = 1. / 15. ; // was  0.066666666 ;
    5855
    59    G4double yTemp[7], dydxTemp[6], yIn[7] ;
     56   G4double yTemp[8], dydxTemp[6], yIn[8] ;
    6057   //  Saving yInput because yInput and yOut can be aliases for same array
    6158   for(i=0;i<nvar;i++) yIn[i]=yInput[i];
    62 
     59   yIn[6] = yInput[6];
     60   yIn[7] = yInput[7];
    6361   G4double h = Step * 0.5;
    6462   hStep=Step;
    … …  
    8280
    8381   h *= 2 ;
    84    StepNoErr(yIn,dydx,h,yTemp,B); // ,beTemp2) ;
     82   StepNoErr(yIn,dydx,h,yTemp,B);
    8583   for(i=0;i<nvar;i++)
    8684   {
    … …  
    127125// is passed from substep to substep.
    128126
    129 void G4RKG3_Stepper::StepNoErr(const G4double tIn[7],
    130                                const G4double dydx[7],
     127void G4RKG3_Stepper::StepNoErr(const G4double tIn[8],
     128                               const G4double dydx[6],
    131129                                     G4double Step,
    132                                      G4double tOut[7],
     130                                     G4double tOut[8],
    133131                                     G4double B[3]      )     // const
    134132   
    … …  
    137135  //  Copy and edit the routine above, to delete alpha2, beta2, ...
    138136   G4double K1[7],K2[7],K3[7],K4[7] ;
    139    G4double tTemp[7], yderiv[6] ;
     137   G4double tTemp[8], yderiv[6] ;
    140138
    141139  // Need Momentum value to give correct values to the coefficients in equation
    … …  
    188186      tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ;
    189187   }
    190  
     188   tOut[6] = tIn[6];
     189   tOut[7] = tIn[7];
    191190   // NormaliseTangentVector( tOut );
    192191 

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4UniformElectricField.cc,v 1.12 2006/06/29 18:24:56 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4UniformElectricField.cc,v 1.13 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    5252                                               G4double vPhi    )
    5353{
    54    if(vField >= 0 &&
    55       vTheta >= 0 && vTheta <= pi &&
    56       vPhi >= 0 && vPhi <= twopi)
    57    {
    58       fFieldComponents[0] = 0.0;
    59       fFieldComponents[1] = 0.0;
    60       fFieldComponents[2] = 0.0;
    61       fFieldComponents[3] = vField*std::sin(vTheta)*std::cos(vPhi) ;
    62       fFieldComponents[4] = vField*std::sin(vTheta)*std::sin(vPhi) ;
    63       fFieldComponents[5] = vField*std::cos(vTheta) ;
    64    }
    65    else
     54   if ( (vField<0) || (vTheta<0) || (vTheta>pi) || (vPhi<0) || (vPhi>twopi) )
    6655   {
    6756      G4Exception("G4UniformElectricField::G4UniformElectricField()",
    6857                  "WrongArgumentValue", FatalException, "Invalid parameters.");
    6958   }
     59
     60   fFieldComponents[0] = 0.0;
     61   fFieldComponents[1] = 0.0;
     62   fFieldComponents[2] = 0.0;
     63   fFieldComponents[3] = vField*std::sin(vTheta)*std::cos(vPhi) ;
     64   fFieldComponents[4] = vField*std::sin(vTheta)*std::sin(vPhi) ;
     65   fFieldComponents[5] = vField*std::cos(vTheta) ;
    7066}
    7167

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

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4UniformMagField.cc,v 1.11 2006/06/29 18:24:58 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4UniformMagField.cc,v 1.12 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    5656                                     G4double vPhi    )
    5757{
    58    if(vField >= 0 &&
    59       vTheta >= 0 && vTheta <= pi &&
    60       vPhi >= 0 && vPhi <= twopi)
    61    {
    62       fFieldComponents[0] = vField*std::sin(vTheta)*std::cos(vPhi) ;
    63       fFieldComponents[1] = vField*std::sin(vTheta)*std::sin(vPhi) ;
    64       fFieldComponents[2] = vField*std::cos(vTheta) ;
    65    }
    66    else
     58   if ( (vField<0) || (vTheta<0) || (vTheta>pi) || (vPhi<0) || (vPhi>twopi) )
    6759   {
    6860      G4Exception("G4UniformMagField::G4UniformMagField()",
    6961                  "WrongArgumentValue", FatalException, "Invalid parameters.") ;
    7062   }
     63   fFieldComponents[0] = vField*std::sin(vTheta)*std::cos(vPhi) ;
     64   fFieldComponents[1] = vField*std::sin(vTheta)*std::sin(vPhi) ;
     65   fFieldComponents[2] = vField*std::cos(vTheta) ;
    7166}
    7267
Note: See TracChangeset for help on using the changeset viewer.