Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4RKG3_Stepper.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/magneticfield/src/G4RKG3_Stepper.cc (Version 11.3.0) and /geometry/magneticfield/src/G4RKG3_Stepper.cc (Version 5.2.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // G4RKG3_Stepper implementation               << 
 27 //                                                 23 //
 28 // Created: J.Apostolakis, V.Grichine - 30.01. <<  24 // $Id: G4RKG3_Stepper.cc,v 1.7 2003/04/02 08:52:21 gcosmo Exp $
 29 // ------------------------------------------- <<  25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
 30                                                <<  26 //
 31 #include "G4RKG3_Stepper.hh"                       27 #include "G4RKG3_Stepper.hh"
                                                   >>  28 #include "G4ThreeVector.hh"
 32 #include "G4LineSection.hh"                        29 #include "G4LineSection.hh"
 33 #include "G4Mag_EqRhs.hh"                      << 
 34                                                << 
 35 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs* Eq << 
 36   : G4MagIntegratorStepper(EqRhs,6)            << 
 37 {                                              << 
 38 }                                              << 
 39                                                    30 
 40 G4RKG3_Stepper::~G4RKG3_Stepper() = default;   <<  31 void G4RKG3_Stepper::Stepper(  const G4double  yInput[7],
 41                                                <<  32                                const G4double dydx[7],
 42 void G4RKG3_Stepper::Stepper( const G4double y <<  33                                      G4double Step,
 43                               const G4double d <<  34                                      G4double yOut[7],
 44                                     G4double S <<  35                                      G4double yErr[])
 45                                     G4double y << 
 46                                     G4double y << 
 47 {                                                  36 {
 48    G4double  B[3];                                 37    G4double  B[3];
                                                   >>  38    //   G4double  yderiv[6];
                                                   >>  39    //   G4double  alpha2, beta2;
 49    G4int nvar = 6 ;                                40    G4int nvar = 6 ;
 50    G4double  by15 = 1. / 15. ; // was  0.06666 <<  41    //   G4double beTemp2, beta2=0;
 51                                                    42 
 52    G4double yTemp[8], dydxTemp[6], yIn[8];     <<  43    G4int i;
                                                   >>  44    G4double  by15 = 1. / 15. ; // was  0.066666666 ;
                                                   >>  45    G4double yTemp[7], dydxTemp[6], yIn[7] ;
                                                   >>  46    //  Saving yInput because yInput and yOut can be aliases for same array
                                                   >>  47    for(i=0;i<nvar;i++) yIn[i]=yInput[i];
 53                                                    48 
 54    // Saving yInput because yInput and yOut ca << 
 55    //                                          << 
 56    for(G4int i=0; i<nvar; ++i)                 << 
 57    {                                           << 
 58      yIn[i]=yInput[i];                         << 
 59    }                                           << 
 60    yIn[6] = yInput[6];                         << 
 61    yIn[7] = yInput[7];                         << 
 62    G4double h = Step * 0.5;                        49    G4double h = Step * 0.5; 
 63    hStep = Step;                               << 
 64      // Do two half steps                      << 
 65                                                    50 
 66    StepNoErr(yIn, dydx,h, yTemp,B) ;           <<  51    // Do two half steps
 67                                                << 
 68    // Store Bfld for DistChord Calculation     << 
 69    //                                          << 
 70    for(auto i=0; i<3; ++i)                     << 
 71    {                                           << 
 72      BfldIn[i] = B[i];                         << 
 73    }                                           << 
 74    // RightHandSide(yTemp,dydxTemp) ;          << 
 75                                                    52 
                                                   >>  53 
                                                   >>  54    // To obtain B1 ...
                                                   >>  55    // GetEquationOfMotion()->GetFieldValue(yIn,B);
                                                   >>  56    // G4RKG3_Stepper::StepWithEst(yIn, dydx, Step, yOut,alpha2, beta2, B1, B2 );
                                                   >>  57 
                                                   >>  58    StepNoErr(yIn, dydx,h, yTemp,B) ;
                                                   >>  59                                      //   RightHandSide(yTemp,dydxTemp) ;
 76    GetEquationOfMotion()->EvaluateRhsGivenB(yT     60    GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;  
 77    StepNoErr(yTemp,dydxTemp,h,yOut,B);         <<  61    StepNoErr(yTemp,dydxTemp,h,yOut,B);              // ,beTemp2) ;
 78                                                <<  62                                                      //   beta2 += beTemp2;
                                                   >>  63                                                      //   beta2 *= 0.5;
                                                   >>  64 
 79    // Store midpoint, chord calculation            65    // Store midpoint, chord calculation
 80                                                    66                                  
 81    fyMidPoint = G4ThreeVector(yTemp[0],  yTemp <<  67    fyMidPoint = G4ThreeVector( yTemp[0],  yTemp[1],  yTemp[2]); 
 82                                                    68 
 83    // Do a full Step                               69    // Do a full Step
 84   //                                           <<  70 
 85    h *= 2 ;                                        71    h *= 2 ;
 86    StepNoErr(yIn,dydx,h,yTemp,B);              <<  72    StepNoErr(yIn,dydx,h,yTemp,B); // ,beTemp2) ;
 87    for(G4int i=0; i<nvar; ++i)                 <<  73    for(i=0;i<nvar;i++)
 88    {                                               74    {
 89       yErr[i] = yOut[i] - yTemp[i] ;               75       yErr[i] = yOut[i] - yTemp[i] ;
 90       yOut[i] += yErr[i]*by15 ;          // Pr     76       yOut[i] += yErr[i]*by15 ;          // Provides 5th order of accuracy
 91    }                                               77    }
 92                                                    78 
 93    // Store values for DistChord method        <<  79 //    for(i=0;i<ncomp;i++)
 94    //                                          <<  80 //    {
 95    fyInitial = G4ThreeVector( yIn[0],   yIn[1] <<  81 //       fyInitial[i]  = yIn[i]; 
 96    fpInitial = G4ThreeVector( yIn[3],   yIn[4] <<  82 //       fyFinal[i]    = yOut[i]; 
                                                   >>  83 //    }
                                                   >>  84 
                                                   >>  85    fyInitial = G4ThreeVector( yIn[0],   yIn[1],   yIn[2]); 
 97    fyFinal   = G4ThreeVector( yOut[0],  yOut[1     86    fyFinal   = G4ThreeVector( yOut[0],  yOut[1],  yOut[2]); 
                                                   >>  87    //   beta2 += beTemp2 ;
                                                   >>  88    //   beta2 *= 0.5 ;   
                                                   >>  89    // NormaliseTangentVector( yOut );  // Deleted
                                                   >>  90    return ;
 98 }                                                  91 }
 99                                                    92 
100 // -------------------------------------------     93 // ---------------------------------------------------------------------------
101                                                    94 
102 // Integrator for RK from G3 with evaluation o     95 // Integrator for RK from G3 with evaluation of error in solution and delta
103 // geometry based on naive similarity with the     96 // geometry based on naive similarity with the case of uniform magnetic field.
104 // B1[3] is input  and is the first magnetic f     97 // B1[3] is input  and is the first magnetic field values
105 // B2[3] is output and is the final magnetic f     98 // B2[3] is output and is the final magnetic field values.
106 //                                             <<  99 
107 void G4RKG3_Stepper::StepWithEst( const G4doub    100 void G4RKG3_Stepper::StepWithEst( const G4double*,
108                                   const G4doub    101                                   const G4double*,
109                                         G4doub    102                                         G4double,
110                                         G4doub    103                                         G4double*,
111                                         G4doub    104                                         G4double&,
112                                         G4doub    105                                         G4double&,
113                                   const G4doub    106                                   const G4double*,
114                                         G4doub    107                                         G4double* )
115                                                   108    
116 {                                                 109 {
117   G4Exception("G4RKG3_Stepper::StepWithEst()", << 110   G4Exception(" ERROR - G4RKG3_Stepper::StepWithEst(): method no longer used.");
118               FatalException, "Method no longe << 111   return ;
119 }                                                 112 }
120                                                   113 
121 // -------------------------------------------    114 // -----------------------------------------------------------------
122                                                   115 
123 // Integrator RK Stepper from G3 with only two    116 // Integrator RK Stepper from G3 with only two field evaluation per Step. 
124 // It is used in propagation initial Step by s    117 // It is used in propagation initial Step by small substeps after solution 
125 // error and delta geometry considerations. B[    118 // error and delta geometry considerations. B[3] is magnetic field which 
126 // is passed from substep to substep.             119 // is passed from substep to substep.
127 //                                             << 120 
128 void G4RKG3_Stepper::StepNoErr(const G4double  << 121 void G4RKG3_Stepper::StepNoErr(const G4double tIn[7],
129                                const G4double  << 122                                const G4double dydx[7],
130                                      G4double     123                                      G4double Step,
131                                      G4double  << 124                                      G4double tOut[7],
132                                      G4double  << 125                                      G4double B[3]      )     // const
133                                                   126    
134 {                                              << 127 {
135                                                << 128   //  Copy and edit the routine above, to delete alpha2, beta2, ...
136    // Copy and edit the routine above, to dele << 129    G4double K1[7],K2[7],K3[7],K4[7] ;
137    //                                          << 130    G4double tTemp[7], yderiv[6] ;
138    G4double K1[7], K2[7], K3[7], K4[7];        << 131    G4int i ;
139    G4double tTemp[8]={0.0}, yderiv[6]={0.0};   << 132 
140                                                << 133 #ifdef END_CODE_G3STEPPER
141    // Need Momentum value to give correct valu << 134    G4Exception(" G4RKG3_Stepper::StepNoErr(): method to be no longer used.");
142    // equation. Integration on unit velocity,  << 135 #else
143                                                << 136    // GetEquationOfMotion()->EvaluateRhsReturnB(tIn,dydx,B1) ;
144    G4double mom, inverse_mom;                  << 137    
145    const G4double c1=0.5, c2=0.125, c3=1./6.;  << 138    for(i=0;i<3;i++)
146                                                << 
147    // Correction for momentum not a velocity   << 
148    // Need the protection !!! must be not zero << 
149    //                                          << 
150    mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4] << 
151    inverse_mom = 1./mom;                       << 
152    for(auto i=0; i<3; ++i)                     << 
153    {                                              139    {
154       K1[i] = Step * dydx[i+3]*inverse_mom;    << 140       K1[i] = Step * dydx[i+3];
155       tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*in << 141       tTemp[i] = tIn[i] + Step*(0.5*tIn[i+3] + 0.125*K1[i]) ;
156       tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ;   << 142       tTemp[i+3] = tIn[i+3] + 0.5*K1[i] ;
157    }                                              143    }
158                                                << 
159    GetEquationOfMotion()->EvaluateRhsReturnB(t    144    GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
160                                                << 145      //  Calculates yderiv & returns B too!
161    for(auto i=0; i<3; ++i)                     << 146 
                                                   >> 147    for(i=0;i<3;i++)
162    {                                              148    {
163       K2[i] = Step * yderiv[i+3]*inverse_mom;  << 149       K2[i] = Step * yderiv[i+3];
164       tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ;   << 150       tTemp[i+3] = tIn[i+3] + 0.5*K2[i] ;
165    }                                              151    }
166                                                << 152 
167    // Given B, calculate yderiv !              << 153    //  Given B, calculate yderiv !
168    //                                          << 
169    GetEquationOfMotion()->EvaluateRhsGivenB(tT    154    GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;  
170                                                << 155    
171    for(auto i=0; i<3; ++i)                     << 156    for(i=0;i<3;i++)
172    {                                              157    {
173       K3[i] = Step * yderiv[i+3]*inverse_mom;  << 158       K3[i] = Step * yderiv[i+3];
174       tTemp[i] = tIn[i] + Step*(tIn[i+3]*inver << 159       tTemp[i] = tIn[i] + Step*(tIn[i+3] + 0.5*K3[i]) ;
175       tTemp[i+3] = tIn[i+3] + K3[i]*mom ;      << 160       tTemp[i+3] = tIn[i+3] + K3[i] ;
176    }                                              161    }
177                                                   162 
178    // Calculates y-deriv(atives) & returns B t << 163    //  Calculates y-deriv(atives) & returns B too!
179    //                                          << 
180    GetEquationOfMotion()->EvaluateRhsReturnB(t    164    GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;  
181                                                   165 
182    for(auto i=0; i<3; ++i)        // Output tr << 166    for(i=0;i<3;i++)        // Output trajectory vector
183    {                                              167    {
184       K4[i] = Step * yderiv[i+3]*inverse_mom;  << 168       K4[i] = Step * yderiv[i+3];
185       tOut[i] = tIn[i] + Step*(tIn[i+3]*invers << 169       tOut[i] = tIn[i] + Step*(tIn[i+3] + (K1[i] + K2[i] + K3[i])/6.0) ;
186       tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2 << 170       tOut[i+3] = tIn[i+3] + (K1[i] + 2*K2[i] + 2*K3[i] +K4[i])/6.0 ;
187    }                                              171    }
188    tOut[6] = tIn[6];                           << 172    // NormaliseTangentVector( tOut );
189    tOut[7] = tIn[7];                           << 173 #endif
                                                   >> 174    
                                                   >> 175    return ;
190 }                                                 176 }
191                                                   177 
192 // -------------------------------------------    178 // ---------------------------------------------------------------------------
193                                                << 
194 G4double G4RKG3_Stepper::DistChord() const     << 
195 {                                              << 
196    // Soon: must check whether h/R > 2 pi  !!  << 
197    // Method below is good only for < 2 pi     << 
198                                                   179 
199    G4double distChord,distLine;                << 180 G4double G4RKG3_Stepper::DistChord()   const 
200                                                << 181 {
201    if (fyInitial != fyFinal)                   << 182   // Soon: must check whether h/R > 2 pi  !!
202    {                                           << 183   //  Method below is good only for < 2 pi
203       distLine = G4LineSection::Distline(fyMid << 
204       distChord = distLine;                    << 
205    }                                           << 
206    else                                        << 
207    {                                           << 
208       distChord = (fyMidPoint-fyInitial).mag() << 
209    }                                           << 
210                                                   184 
211    return distChord;                           << 185   return G4LineSection::Distline( fyMidPoint, fyInitial, fyFinal );
                                                   >> 186   // This is a class method that gives distance of Mid 
                                                   >> 187   //  from the Chord between the Initial and Final points.
212 }                                                 188 }
213                                                   189