| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT 4 class implementation file
29 // -------------------------------------------
30 //
31
32 #include "G4ErrorSurfaceTrajState.hh"
33 #include "G4ErrorPropagatorData.hh"
34
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4Field.hh"
38 #include "G4FieldManager.hh"
39 #include "G4TransportationManager.hh"
40
41 #include "G4ErrorMatrix.hh"
42
43 #include <iomanip>
44
45 //--------------------------------------------
46 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta
47 const G4String& partType, const G4Point3D& p
48 const G4Vector3D& vecU, const G4Vector3D& ve
49 : G4ErrorTrajState(partType, pos, mom, errma
50 {
51 Init();
52 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo
53 }
54
55 //--------------------------------------------
56 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta
57
58
59
60
61 : G4ErrorTrajState(partType, pos, mom, errma
62 {
63 Init();
64 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo
65 }
66
67 //--------------------------------------------
68 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta
69
70 : G4ErrorTrajState(tpSC.GetParticleType(), t
71 tpSC.GetMomentum())
72 {
73 // fParticleType = tpSC.GetParticleType();
74 // fPosition = tpSC.GetPosition();
75 // fMomentum = tpSC.GetMomentum();
76 fTrajParam = G4ErrorSurfaceTrajParam(fPositi
77 Init();
78
79 //----- Get the error matrix in SC coordinat
80 BuildErrorMatrix(tpSC, GetVectorV(), GetVect
81 }
82
83 //--------------------------------------------
84 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta
85
86
87
88 : G4ErrorTrajState(tpSC.GetParticleType(), t
89 tpSC.GetMomentum())
90 {
91 Init(); // needed to define charge sign
92 fTrajParam = G4ErrorSurfaceTrajParam(fPositi
93 //----- Get the error matrix in SC coordinat
94 transfM = BuildErrorMatrix(tpSC, vecU, vecV)
95 }
96
97 //--------------------------------------------
98 G4ErrorMatrix G4ErrorSurfaceTrajState::BuildEr
99 G4ErrorFreeTrajState& tpSC, const G4Vector3D
100 {
101 G4double sclambda = tpSC.GetParameters().Get
102 G4double scphi = tpSC.GetParameters().Get
103 if(G4ErrorPropagatorData::GetErrorPropagator
104 G4ErrorMode_PropBackwards)
105 {
106 sclambda *= -1;
107 scphi += CLHEP::pi;
108 }
109 G4double cosLambda = std::cos(sclambda);
110 G4double sinLambda = std::sin(sclambda);
111 G4double sinPhi = std::sin(scphi);
112 G4double cosPhi = std::cos(scphi);
113
114 #ifdef G4EVERBOSE
115 if(iverbose >= 4)
116 G4cout << " PM " << fMomentum.mag() << " p
117 << scphi << G4endl;
118 #endif
119
120 G4ThreeVector vTN(cosLambda * cosPhi, cosLam
121 G4ThreeVector vUN(-sinPhi, cosPhi, 0.);
122 G4ThreeVector vVN(-vTN.z() * vUN.y(), vTN.z(
123
124 #ifdef G4EVERBOSE
125 if(iverbose >= 4)
126 G4cout << " SC2SD: vTN " << vTN << " vUN "
127 << G4endl;
128 #endif
129 G4double UJ = vUN * GetVectorV();
130 G4double UK = vUN * GetVectorW();
131 G4double VJ = vVN * GetVectorV();
132 G4double VK = vVN * GetVectorW();
133
134 //--- Get transformation first
135 G4ErrorMatrix transfM(5, 5, 0);
136 //--- Get magnetic field
137 const G4Field* field = G4TransportationManag
138 ->GetFieldManager()
139 ->GetDetectorField(
140
141 G4Vector3D vectorU = GetVectorV().cross(GetV
142 G4double T1R = 1. / (vTN * vectorU);
143
144 #ifdef G4EVERBOSE
145 if(iverbose >= 4)
146 G4cout << "surf vectors:U,V,W " << vectorU
147 << GetVectorW() << " T1R:" << T1R
148 #endif
149
150 if(fCharge != 0 && field)
151 {
152 G4double pos[3];
153 pos[0] = fPosition.x() * cm;
154 pos[1] = fPosition.y() * cm;
155 pos[2] = fPosition.z() * cm;
156 G4double Hd[3];
157 field->GetFieldValue(pos, Hd);
158 G4ThreeVector H =
159 G4ThreeVector(Hd[0], Hd[1], Hd[2]) / tes
160 G4double magH = H.mag();
161 G4double invP = 1. / (fMomentum.mag() / G
162 G4double magHM = magH * invP;
163 if(magH != 0.)
164 {
165 G4double magHM2 = fCharge / magH;
166 G4double Q = -magHM * c_light / (km
167 #ifdef G4EVERBOSE
168 if(iverbose >= 4)
169 G4cout << GeV << " Q " << Q << " magHM
170 << c_light / (km / ns) << G4end
171 #endif
172
173 G4double sinz = -H * vUN * magHM2;
174 G4double cosz = H * vVN * magHM2;
175 G4double T3R = Q * std::pow(T1R, 3);
176 G4double UI = vUN * vectorU;
177 G4double VI = vVN * vectorU;
178 #ifdef G4EVERBOSE
179 if(iverbose >= 4)
180 {
181 G4cout << " T1R " << T1R << " T3R " <<
182 G4cout << " UI " << UI << " VI " << VI
183 << G4endl;
184 G4cout << " UJ " << UJ << " VJ " << VJ
185 G4cout << " UK " << UK << " VK " << VK
186 }
187 #endif
188
189 transfM[1][3] = -UI * (VK * cosz - UK *
190 transfM[1][4] = -VI * (VK * cosz - UK *
191 transfM[2][3] = UI * (VJ * cosz - UJ * s
192 transfM[2][4] = VI * (VJ * cosz - UJ * s
193 }
194 }
195
196 G4double T2R = T1R * T1R;
197 transfM[0][0] = 1.;
198 transfM[1][1] = -UK * T2R;
199 transfM[1][2] = VK * cosLambda * T2R;
200 transfM[2][1] = UJ * T2R;
201 transfM[2][2] = -VJ * cosLambda * T2R;
202 transfM[3][3] = VK * T1R;
203 transfM[3][4] = -UK * T1R;
204 transfM[4][3] = -VJ * T1R;
205 transfM[4][4] = UJ * T1R;
206
207 #ifdef G4EVERBOSE
208 if(iverbose >= 4)
209 G4cout << " SC2SD transf matrix A= " << tr
210 #endif
211 fError = G4ErrorTrajErr(tpSC.GetError().simi
212
213 #ifdef G4EVERBOSE
214 if(iverbose >= 1)
215 G4cout << "G4EP: error matrix SC2SD " << f
216 if(iverbose >= 4)
217 G4cout << "G4ErrorSurfaceTrajState from SC
218 #endif
219
220 return transfM; // want to use trasnfM for
221 }
222
223 //--------------------------------------------
224 void G4ErrorSurfaceTrajState::Init()
225 {
226 theTSType = G4eTS_OS;
227 BuildCharge();
228 }
229
230 //--------------------------------------------
231 void G4ErrorSurfaceTrajState::Dump(std::ostrea
232
233 //--------------------------------------------
234 std::ostream& operator<<(std::ostream& out, co
235 {
236 std::ios::fmtflags oldFlags = out.flags();
237 out.setf(std::ios::fixed, std::ios::floatfie
238
239 ts.DumpPosMomError(out);
240
241 out << " G4ErrorSurfaceTrajState: Params: "
242 out.flags(oldFlags);
243 return out;
244 }
245