| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4HelixExplicitEuler implementation << 26 //
>> 27 // $Id: G4HelixExplicitEuler.cc,v 1.7 2007/08/21 09:43:50 tnikitin Exp $
>> 28 // GEANT4 tag $Name: geant4-09-00-patch-01 $
>> 29 //
27 // 30 //
28 // Helix Explicit Euler: x_1 = x_0 + helix(h) 31 // Helix Explicit Euler: x_1 = x_0 + helix(h)
29 // with helix(h) being a helix piece of lengt << 32 // with helix(h) being a helix piece of length h
30 // Most simple approach for solving linear di << 33 // most simple approach for solving linear differential equations.
31 // Take the current derivative and add it to 34 // Take the current derivative and add it to the current position.
32 // 35 //
33 // Author: W.Wander <wwc@mit.edu>, 12.09.1997 << 36 // W.Wander <wwc@mit.edu> 12/09/97
34 // ------------------------------------------- 37 // -------------------------------------------------------------------
35 38
36 #include "G4HelixExplicitEuler.hh" 39 #include "G4HelixExplicitEuler.hh"
37 #include "G4PhysicalConstants.hh" <<
38 #include "G4ThreeVector.hh" 40 #include "G4ThreeVector.hh"
39 41
40 G4HelixExplicitEuler::G4HelixExplicitEuler(G4M <<
41 : G4MagHelicalStepper(EqRhs) <<
42 { <<
43 } <<
44 <<
45 G4HelixExplicitEuler::~G4HelixExplicitEuler() <<
46 42
47 void G4HelixExplicitEuler::Stepper( const G4do << 43 void G4HelixExplicitEuler::Stepper( const G4double yInput[7],
48 const G4do << 44 const G4double*,
49 G4do << 45 G4double Step,
50 G4do << 46 G4double yOut[7],
51 G4do << 47 G4double yErr[])
>> 48
52 { 49 {
53 // Estimation of the Stepping Angle << 50
54 // << 51 //Estimation of the Stepping Angle
>> 52
55 G4ThreeVector Bfld; 53 G4ThreeVector Bfld;
56 MagFieldEvaluate(yInput, Bfld); 54 MagFieldEvaluate(yInput, Bfld);
57 55
58 const G4int nvar = 6 ; 56 const G4int nvar = 6 ;
59 G4double yTemp[8], yIn[8] ; << 57 G4int i;
>> 58 G4double yTemp[7], yIn[7] ;
60 G4ThreeVector Bfld_midpoint; 59 G4ThreeVector Bfld_midpoint;
61 << 60 // Saving yInput because yInput and yOut can be aliases for same array
62 // Saving yInput because yInput and yOut can << 61 for(i=0;i<nvar;i++) yIn[i]=yInput[i];
63 // <<
64 for(G4int i=0; i<nvar; ++i) <<
65 { <<
66 yIn[i] = yInput[i]; <<
67 } <<
68 62
69 G4double h = Step * 0.5; << 63 G4double h = Step * 0.5;
70 64
71 // Do full step and two half steps << 65 // Do full step and two half steps
72 // << 66 G4double yTemp2[7];
73 G4double yTemp2[7]; << 67 AdvanceHelix(yIn, Bfld, h, yTemp2,yTemp);
74 AdvanceHelix(yIn, Bfld, h, yTemp2,yTemp); << 68 MagFieldEvaluate(yTemp2, Bfld_midpoint) ;
75 MagFieldEvaluate(yTemp2, Bfld_midpoint) ; << 69 AdvanceHelix(yTemp2, Bfld_midpoint, h, yOut);
76 AdvanceHelix(yTemp2, Bfld_midpoint, h, yOut) << 70
77 SetAngCurve(GetAngCurve() * 2); << 71 // Error estimation
>> 72 for(i=0;i<nvar;i++) {
>> 73 yErr[i] = yOut[i] - yTemp[i] ;
>> 74 }
78 75
79 // Error estimation <<
80 // <<
81 for(G4int i=0; i<nvar; ++i) <<
82 { <<
83 yErr[i] = yOut[i] - yTemp[i]; <<
84 } <<
85 } 76 }
86 77
87 G4double G4HelixExplicitEuler::DistChord() c 78 G4double G4HelixExplicitEuler::DistChord() const
88 { 79 {
89 // Implementation : must check whether h/R > 80 // Implementation : must check whether h/R > 2 pi !!
90 // If( h/R < pi) use G4LineSection::DistL 81 // If( h/R < pi) use G4LineSection::DistLine
91 // Else DistChord=R_helix 82 // Else DistChord=R_helix
92 // 83 //
93 G4double distChord; 84 G4double distChord;
94 G4double Ang_curve=GetAngCurve(); 85 G4double Ang_curve=GetAngCurve();
95 86
96 << 87
97 if(Ang_curve<=pi) << 88 if(Ang_curve<=pi){
98 { << 89 distChord=GetRadHelix()*(1-cos(0.5*Ang_curve));
99 distChord=GetRadHelix()*(1-std::cos(0.5*An << 90 }
100 } << 91 else
101 else if(Ang_curve<twopi) << 92 if(Ang_curve<twopi){
102 { << 93 distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
103 distChord=GetRadHelix()*(1+std::cos(0.5*(t << 94 }
104 } << 95 else{
105 else << 96 distChord=2.*GetRadHelix();
106 { << 97 }
107 distChord=2.*GetRadHelix(); <<
108 } <<
109 98
110 return distChord; 99 return distChord;
>> 100
111 } 101 }
112 << 102 void
113 void G4HelixExplicitEuler::DumbStepper( const << 103 G4HelixExplicitEuler::DumbStepper( const G4double yIn[],
114 << 104 G4ThreeVector Bfld,
115 << 105 G4double h,
116 << 106 G4double yOut[])
117 { 107 {
118 AdvanceHelix(yIn, Bfld, h, yOut); << 108
>> 109 AdvanceHelix(yIn, Bfld, h, yOut);
>> 110
119 } 111 }
120 112