ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HelixMixedStepper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4HelixMixedStepper.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // class G4HelixMixedStepper
27 //
28 // Class description:
29 //
30 // G4HelixMixedStepper split the Method used for Integration in two:
31 //
32 // If Stepping Angle ( h / R_curve) < pi/3
33 // use Stepper for small step(ClassicalRK4 by default)
34 // Else use HelixExplicitEuler Stepper
35 //
36 // Created: T.Nikitina, CERN - 18.05.2007, derived from G4ExactHelicalStepper
37 // -------------------------------------------------------------------------
38 
39 #include "G4HelixMixedStepper.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4ClassicalRK4.hh"
42 #include "G4CashKarpRKF45.hh"
43 #include "G4SimpleRunge.hh"
44 #include "G4HelixImplicitEuler.hh"
45 #include "G4HelixExplicitEuler.hh"
46 #include "G4HelixSimpleRunge.hh"
47 #include "G4ExactHelixStepper.hh"
48 #include "G4ExplicitEuler.hh"
49 #include "G4ImplicitEuler.hh"
50 #include "G4SimpleHeum.hh"
51 #include "G4RKG3_Stepper.hh"
52 #include "G4NystromRK4.hh"
53 
54 // Additional potential steppers
55 #include "G4DormandPrince745.hh"
56 #include "G4BogackiShampine23.hh"
57 #include "G4BogackiShampine45.hh"
58 #include "G4TsitourasRK45.hh"
59 
60 #include "G4ThreeVector.hh"
61 #include "G4LineSection.hh"
62 
63 // ---------------------------------------------------------------------------
66  G4int stepperNumber,
67  G4double angleThreshold)
68  : G4MagHelicalStepper(EqRhs)
69 {
70  if( angleThreshold < 0.0 )
71  {
72  fAngle_threshold = (1.0/3.0)*pi;
73  }
74  else
75  {
76  fAngle_threshold = angleThreshold;
77  }
78 
79  if(stepperNumber<0)
80  {
81  // stepperNumber = 4; // Default is RK4 (original)
82  stepperNumber = 745; // Default is DormandPrince745 (ie DoPri5)
83  // stepperNumber = 8; // Default is CashKarp
84  }
85 
86  fStepperNumber = stepperNumber; // Store the choice
88 }
89 
90 // ---------------------------------------------------------------------------
92 {
93  delete fRK4Stepper;
94  if (fVerbose>0) { PrintCalls(); }
95 }
96 
97 // ---------------------------------------------------------------------------
98 void G4HelixMixedStepper::Stepper( const G4double yInput[7],
99  const G4double dydx[7],
100  G4double Step,
101  G4double yOut[7],
102  G4double yErr[])
103 {
104  // Estimation of the Stepping Angle
105  //
106  G4ThreeVector Bfld;
107  MagFieldEvaluate(yInput, Bfld);
108 
109  G4double Bmag = Bfld.mag();
110  const G4double* pIn = yInput+3;
111  G4ThreeVector initVelocity = G4ThreeVector( pIn[0], pIn[1], pIn[2] );
112  G4double velocityVal = initVelocity.mag();
113 
114  const G4double R_1 = std::abs(GetInverseCurve(velocityVal,Bmag));
115  // curv = inverse Radius
116  G4double Ang_curve = R_1 * Step;
117  // SetAngCurve(Ang_curve);
118  // SetCurve(std::abs(1/R_1));
119 
120  if(Ang_curve < fAngle_threshold)
121  {
122  ++fNumCallsRK4;
123  fRK4Stepper->Stepper(yInput,dydx,Step,yOut,yErr);
124  }
125  else
126  {
127  constexpr G4int nvar = 6 ;
128  constexpr G4int nvarMax = 8 ;
129  G4double yTemp[nvarMax], yIn[nvarMax], yTemp2[nvarMax];
130  G4ThreeVector Bfld_midpoint;
131 
132  SetAngCurve(Ang_curve);
133  SetCurve(std::abs(1.0/R_1));
134  ++fNumCallsHelix;
135 
136  // Saving yInput because yInput and yOut can be aliases for same array
137  //
138  for(G4int i=0; i<nvar; ++i)
139  {
140  yIn[i]=yInput[i];
141  }
142 
143  G4double halfS = Step * 0.5;
144 
145  // 1. Do first half step and full step
146  //
147  AdvanceHelix(yIn, Bfld, halfS, yTemp, yTemp2); // yTemp2 for s=2*h (halfS)
148 
149  MagFieldEvaluate(yTemp, Bfld_midpoint) ;
150 
151  // 2. Do second half step - with revised field
152  // NOTE: Could avoid this call if 'Bfld_midpoint == Bfld'
153  // or diff 'almost' zero
154  //
155  AdvanceHelix(yTemp, Bfld_midpoint, halfS, yOut);
156  // Not requesting y at s=2*h (halfS)
157 
158  // 3. Estimate the integration error
159  // should be (nearly) zero if Bfield= constant
160  //
161  for(G4int i=0; i<nvar; ++i)
162  {
163  yErr[i] = yOut[i] - yTemp2[i];
164  }
165  }
166 }
167 
168 // ---------------------------------------------------------------------------
170  G4ThreeVector Bfld,
171  G4double h,
172  G4double yOut[] )
173 {
174  AdvanceHelix(yIn, Bfld, h, yOut);
175 }
176 
177 // ---------------------------------------------------------------------------
179 {
180  // Implementation : must check whether h/R > 2 pi !!
181  // If( h/R < pi) use G4LineSection::DistLine
182  // Else DistChord=R_helix
183  //
184  G4double distChord;
185  G4double Ang_curve=GetAngCurve();
186 
187  if(Ang_curve<=pi)
188  {
189  distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
190  }
191  else
192  {
193  if(Ang_curve<twopi)
194  {
195  distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
196  }
197  else
198  {
199  distChord=2.*GetRadHelix();
200  }
201  }
202 
203  return distChord;
204 }
205 
206 // ---------------------------------------------------------------------------
208 {
209  G4cout << "In HelixMixedStepper::Number of calls to smallStepStepper = "
210  << fNumCallsRK4
211  << " and Number of calls to Helix = " << fNumCallsHelix << G4endl;
212 }
213 
214 // ---------------------------------------------------------------------------
217 {
218  G4MagIntegratorStepper* pStepper;
219  if (fVerbose>0) G4cout << " G4HelixMixedStepper: ";
220  switch ( StepperNumber )
221  {
222  // Robust, classic method
223  case 4:
224  pStepper = new G4ClassicalRK4( pE );
225  if (fVerbose>0) G4cout << "G4ClassicalRK4";
226  break;
227 
228  // Steppers with embedded estimation of error
229  case 8:
230  pStepper = new G4CashKarpRKF45( pE );
231  if (fVerbose>0) G4cout << "G4CashKarpRKF45";
232  break;
233  case 13:
234  pStepper = new G4NystromRK4( pE );
235  if (fVerbose>0) G4cout << "G4NystromRK4";
236  break;
237 
238  // Lowest order RK Stepper - experimental
239  case 1:
240  pStepper = new G4ImplicitEuler( pE );
241  if (fVerbose>0) G4cout << "G4ImplicitEuler";
242  break;
243 
244  // Lower order RK Steppers - ok overall, good for uneven fields
245  case 2:
246  pStepper = new G4SimpleRunge( pE );
247  if (fVerbose>0) G4cout << "G4SimpleRunge";
248  break;
249  case 3:
250  pStepper = new G4SimpleHeum( pE );
251  if (fVerbose>0) G4cout << "G4SimpleHeum";
252  break;
253  case 23:
254  pStepper = new G4BogackiShampine23( pE );
255  if (fVerbose>0) G4cout << "G4BogackiShampine23";
256  break;
257 
258  // Higher order RK Steppers
259  // for smoother fields and high accuracy requirements
260  case 45:
261  pStepper = new G4BogackiShampine45( pE );
262  if (fVerbose>0) G4cout << "G4BogackiShampine45";
263  break;
264  case 145:
265  pStepper = new G4TsitourasRK45( pE );
266  if (fVerbose>0) G4cout << "G4TsitourasRK45";
267  break;
268  case 745:
269  pStepper = new G4DormandPrince745( pE );
270  if (fVerbose>0) G4cout << "G4DormandPrince745";
271  break;
272 
273  // Helical Steppers
274  case 6:
275  pStepper = new G4HelixImplicitEuler( pE );
276  if (fVerbose>0) G4cout << "G4HelixImplicitEuler";
277  break;
278  case 7:
279  pStepper = new G4HelixSimpleRunge( pE );
280  if (fVerbose>0) G4cout << "G4HelixSimpleRunge";
281  break;
282  case 5:
283  pStepper = new G4HelixExplicitEuler( pE );
284  if (fVerbose>0) G4cout << "G4HelixExplicitEuler";
285  break; // Since Helix Explicit is used for long steps,
286  // this is useful only to measure overhead.
287  // Exact Helix - likely good only for cases of
288  // i) uniform field (potentially over small distances)
289  // ii) segmented uniform field (maybe)
290  case 9:
291  pStepper = new G4ExactHelixStepper( pE );
292  if (fVerbose>0) G4cout << "G4ExactHelixStepper";
293  break;
294  case 10:
295  pStepper = new G4RKG3_Stepper( pE );
296  if (fVerbose>0) G4cout << "G4RKG3_Stepper";
297  break;
298 
299  // Low Order Steppers - not good except for very weak fields
300  case 11:
301  pStepper = new G4ExplicitEuler( pE );
302  if (fVerbose>0) G4cout << "G4ExplicitEuler";
303  break;
304  case 12:
305  pStepper = new G4ImplicitEuler( pE );
306  if (fVerbose>0) G4cout << "G4ImplicitEuler";
307  break;
308 
309  case 0:
310  case -1:
311  default:
312  pStepper = new G4DormandPrince745( pE ); // Was G4ClassicalRK4( pE );
313  if (fVerbose>0) G4cout << "G4DormandPrince745 (Default)";
314  break;
315  }
316 
317  if(fVerbose>0)
318  {
319  G4cout << " chosen as stepper for small steps in G4HelixMixedStepper."
320  << G4endl;
321  }
322 
323  return pStepper;
324 }