| 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 // 26 //
27 // ------------------------------------------- 27 // ------------------------------------------------------------
28 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file
29 // ------------------------------------------- 29 // ------------------------------------------------------------
30 // 30 //
31 31
32 #include <iomanip> 32 #include <iomanip>
33 33
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "G4Field.hh" 36 #include "G4Field.hh"
37 #include "G4FieldManager.hh" 37 #include "G4FieldManager.hh"
38 #include "G4TransportationManager.hh" 38 #include "G4TransportationManager.hh"
39 #include "G4GeometryTolerance.hh" 39 #include "G4GeometryTolerance.hh"
40 #include "G4Material.hh" 40 #include "G4Material.hh"
41 #include "G4ErrorPropagatorData.hh" 41 #include "G4ErrorPropagatorData.hh"
42 #include "G4ErrorFreeTrajState.hh" 42 #include "G4ErrorFreeTrajState.hh"
43 #include "G4ErrorFreeTrajParam.hh" 43 #include "G4ErrorFreeTrajParam.hh"
44 #include "G4ErrorSurfaceTrajState.hh" 44 #include "G4ErrorSurfaceTrajState.hh"
45 #include "G4ErrorMatrix.hh" 45 #include "G4ErrorMatrix.hh"
46 46
47 //-------------------------------------------- 47 //------------------------------------------------------------------------
48 G4ErrorFreeTrajState::G4ErrorFreeTrajState(con 48 G4ErrorFreeTrajState::G4ErrorFreeTrajState(const G4String& partName,
49 con 49 const G4Point3D& pos,
50 con 50 const G4Vector3D& mom,
51 con 51 const G4ErrorTrajErr& errmat)
52 : G4ErrorTrajState(partName, pos, mom, errma 52 : G4ErrorTrajState(partName, pos, mom, errmat)
53 { 53 {
54 fTrajParam = G4ErrorFreeTrajParam(pos, mom); 54 fTrajParam = G4ErrorFreeTrajParam(pos, mom);
55 Init(); 55 Init();
56 } 56 }
57 57
58 //-------------------------------------------- 58 //------------------------------------------------------------------------
59 G4ErrorFreeTrajState::G4ErrorFreeTrajState(con 59 G4ErrorFreeTrajState::G4ErrorFreeTrajState(const G4ErrorSurfaceTrajState& tpSD)
60 : G4ErrorTrajState(tpSD.GetParticleType(), t 60 : G4ErrorTrajState(tpSD.GetParticleType(), tpSD.GetPosition(),
61 tpSD.GetMomentum()) 61 tpSD.GetMomentum())
62 { 62 {
63 // G4ThreeVector planeNormal = tpSD.GetPlane 63 // G4ThreeVector planeNormal = tpSD.GetPlaneNormal();
64 // G4double fPt = tpSD.GetMomentum()*planeNo 64 // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to
65 // plane G4ErrorSurfaceTrajParam tpSDparam = 65 // plane G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
66 // G4ThreeVector Psc = fPt * planeNormal + 66 // G4ThreeVector Psc = fPt * planeNormal +
67 // tpSDparam.GetPU()*tpSDparam.GetVectorU() 67 // tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW();
68 68
69 fTrajParam = G4ErrorFreeTrajParam(fPosition, 69 fTrajParam = G4ErrorFreeTrajParam(fPosition, fMomentum);
70 Init(); 70 Init();
71 71
72 //----- Get the error matrix in SC coordinat 72 //----- Get the error matrix in SC coordinates
73 G4ErrorSurfaceTrajParam tpSDparam = tpSD.Get 73 G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
74 G4double mom = fMomentu 74 G4double mom = fMomentum.mag();
75 G4double mom2 = fMomentu 75 G4double mom2 = fMomentum.mag2();
76 G4double TVW1 = 76 G4double TVW1 =
77 std::sqrt(mom2 / (mom2 + tpSDparam.GetPV() 77 std::sqrt(mom2 / (mom2 + tpSDparam.GetPV() * tpSDparam.GetPV() +
78 tpSDparam.GetPW() * tpSD 78 tpSDparam.GetPW() * tpSDparam.GetPW()));
79 G4ThreeVector vTVW(TVW1, tpSDparam.GetPV() / 79 G4ThreeVector vTVW(TVW1, tpSDparam.GetPV() / mom * TVW1,
80 tpSDparam.GetPW() / mom * 80 tpSDparam.GetPW() / mom * TVW1);
81 G4Vector3D vectorU = tpSDparam.GetVectorV(). 81 G4Vector3D vectorU = tpSDparam.GetVectorV().cross(tpSDparam.GetVectorW());
82 G4Vector3D vTN = vTVW.x() * vectorU + vT 82 G4Vector3D vTN = vTVW.x() * vectorU + vTVW.y() * tpSDparam.GetVectorV() +
83 vTVW.z() * tpSDparam.GetVec 83 vTVW.z() * tpSDparam.GetVectorW();
84 84
85 #ifdef G4EVERBOSE 85 #ifdef G4EVERBOSE
86 if(iverbose >= 5) 86 if(iverbose >= 5)
87 { 87 {
88 G4double pc2 = std::asin(vTN.z()); 88 G4double pc2 = std::asin(vTN.z());
89 G4double pc3 = std::atan(vTN.y() / vTN.x() 89 G4double pc3 = std::atan(vTN.y() / vTN.x());
90 90
91 G4cout << " CHECK: pc2 " << pc2 << " = " < 91 G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda()
92 << " diff " << pc2 - GetParameters( 92 << " diff " << pc2 - GetParameters().GetLambda() << G4endl;
93 G4cout << " CHECK: pc3 " << pc3 << " = " < 93 G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi()
94 << " diff " << pc3 - GetParameters( 94 << " diff " << pc3 - GetParameters().GetPhi() << G4endl;
95 } 95 }
96 #endif 96 #endif
97 97
98 //--- Get the unit vectors perp to P 98 //--- Get the unit vectors perp to P
99 G4double cosl = std::cos(GetParameters().Get 99 G4double cosl = std::cos(GetParameters().GetLambda());
100 if(cosl < 1.E-30) 100 if(cosl < 1.E-30)
101 cosl = 1.E-30; 101 cosl = 1.E-30;
102 G4double cosl1 = 1. / cosl; 102 G4double cosl1 = 1. / cosl;
103 G4Vector3D vUN(-vTN.y() * cosl1, vTN.x() * c 103 G4Vector3D vUN(-vTN.y() * cosl1, vTN.x() * cosl1, 0.);
104 G4Vector3D vVN(-vTN.z() * vUN.y(), vTN.z() * 104 G4Vector3D vVN(-vTN.z() * vUN.y(), vTN.z() * vUN.x(), cosl);
105 105
106 G4Vector3D vUperp = G4Vector3D(-fMomentum.y( 106 G4Vector3D vUperp = G4Vector3D(-fMomentum.y(), fMomentum.x(), 0.);
107 G4Vector3D vVperp = vUperp.cross(fMomentum); 107 G4Vector3D vVperp = vUperp.cross(fMomentum);
108 vUperp *= 1. / vUperp.mag(); 108 vUperp *= 1. / vUperp.mag();
109 vVperp *= 1. / vVperp.mag(); 109 vVperp *= 1. / vVperp.mag();
110 110
111 #ifdef G4EVERBOSE 111 #ifdef G4EVERBOSE
112 if(iverbose >= 5) 112 if(iverbose >= 5)
113 { 113 {
114 G4cout << " CHECK: vUN " << vUN << " = " < 114 G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff "
115 << (vUN - vUperp).mag() << G4endl; 115 << (vUN - vUperp).mag() << G4endl;
116 G4cout << " CHECK: vVN " << vVN << " = " < 116 G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff "
117 << (vVN - vVperp).mag() << G4endl; 117 << (vVN - vVperp).mag() << G4endl;
118 } 118 }
119 #endif 119 #endif
120 120
121 // get the dot products of vectors perpendic 121 // get the dot products of vectors perpendicular to direction and vector
122 // defining SD plane 122 // defining SD plane
123 G4double dUU = vUperp * tpSD.GetVectorV(); 123 G4double dUU = vUperp * tpSD.GetVectorV();
124 G4double dUV = vUperp * tpSD.GetVectorW(); 124 G4double dUV = vUperp * tpSD.GetVectorW();
125 G4double dVU = vVperp * tpSD.GetVectorV(); 125 G4double dVU = vVperp * tpSD.GetVectorV();
126 G4double dVV = vVperp * tpSD.GetVectorW(); 126 G4double dVV = vVperp * tpSD.GetVectorW();
127 127
128 //--- Get transformation first 128 //--- Get transformation first
129 G4ErrorMatrix transfM(5, 5, 1); 129 G4ErrorMatrix transfM(5, 5, 1);
130 //--- Get magnetic field 130 //--- Get magnetic field
131 const G4Field* field = G4TransportationManag 131 const G4Field* field = G4TransportationManager::GetTransportationManager()
132 ->GetFieldManager() 132 ->GetFieldManager()
133 ->GetDetectorField( 133 ->GetDetectorField();
134 G4ThreeVector dir = fTrajParam.GetDirecti 134 G4ThreeVector dir = fTrajParam.GetDirection();
135 G4double invCosTheta = 1. / std::cos(dir.the 135 G4double invCosTheta = 1. / std::cos(dir.theta());
136 G4cout << " dir=" << dir << " invCosTheta " 136 G4cout << " dir=" << dir << " invCosTheta " << invCosTheta << G4endl;
137 137
138 if(fCharge != 0 && field) 138 if(fCharge != 0 && field)
139 { 139 {
140 G4double pos1[3]; 140 G4double pos1[3];
141 pos1[0] = fPosition.x() * cm; 141 pos1[0] = fPosition.x() * cm;
142 pos1[1] = fPosition.y() * cm; 142 pos1[1] = fPosition.y() * cm;
143 pos1[2] = fPosition.z() * cm; 143 pos1[2] = fPosition.z() * cm;
144 G4double h1[3]; 144 G4double h1[3];
145 field->GetFieldValue(pos1, h1); 145 field->GetFieldValue(pos1, h1);
146 G4ThreeVector HPre = G4ThreeVector(h1[0], 146 G4ThreeVector HPre = G4ThreeVector(h1[0], h1[1], h1[2]) / tesla * 10.;
147 G4double magHPre = HPre.mag(); 147 G4double magHPre = HPre.mag();
148 G4double invP = 1. / fMomentum.mag(); 148 G4double invP = 1. / fMomentum.mag();
149 G4double magHPreM = magHPre * invP; 149 G4double magHPreM = magHPre * invP;
150 if(magHPre != 0.) 150 if(magHPre != 0.)
151 { 151 {
152 G4double magHPreM2 = fCharge / magHPre; 152 G4double magHPreM2 = fCharge / magHPre;
153 153
154 G4double Q = -magHPreM * c_light; 154 G4double Q = -magHPreM * c_light;
155 G4double sinz = -HPre * vUperp * magHPre 155 G4double sinz = -HPre * vUperp * magHPreM2;
156 G4double cosz = HPre * vVperp * magHPreM 156 G4double cosz = HPre * vVperp * magHPreM2;
157 157
158 transfM[1][3] = -Q * dir.y() * sinz; 158 transfM[1][3] = -Q * dir.y() * sinz;
159 transfM[1][4] = -Q * dir.z() * sinz; 159 transfM[1][4] = -Q * dir.z() * sinz;
160 transfM[2][3] = -Q * dir.y() * cosz * in 160 transfM[2][3] = -Q * dir.y() * cosz * invCosTheta;
161 transfM[2][4] = -Q * dir.z() * cosz * in 161 transfM[2][4] = -Q * dir.z() * cosz * invCosTheta;
162 } 162 }
163 } 163 }
164 164
165 transfM[0][0] = 1.; 165 transfM[0][0] = 1.;
166 transfM[1][1] = dir.x() * dVU; 166 transfM[1][1] = dir.x() * dVU;
167 transfM[1][2] = dir.x() * dVV; 167 transfM[1][2] = dir.x() * dVV;
168 transfM[2][1] = dir.x() * dUU * invCosTheta; 168 transfM[2][1] = dir.x() * dUU * invCosTheta;
169 transfM[2][2] = dir.x() * dUV * invCosTheta; 169 transfM[2][2] = dir.x() * dUV * invCosTheta;
170 transfM[3][3] = dUU; 170 transfM[3][3] = dUU;
171 transfM[3][4] = dUV; 171 transfM[3][4] = dUV;
172 transfM[4][3] = dVU; 172 transfM[4][3] = dVU;
173 transfM[4][4] = dVV; 173 transfM[4][4] = dVV;
174 174
175 fError = G4ErrorTrajErr(tpSD.GetError().simi 175 fError = G4ErrorTrajErr(tpSD.GetError().similarity(transfM));
176 176
177 #ifdef G4EVERBOSE 177 #ifdef G4EVERBOSE
178 if(iverbose >= 1) 178 if(iverbose >= 1)
179 G4cout << "error matrix SD2SC " << fError 179 G4cout << "error matrix SD2SC " << fError << G4endl;
180 if(iverbose >= 4) 180 if(iverbose >= 4)
181 G4cout << "G4ErrorFreeTrajState from SD " 181 G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl;
182 #endif 182 #endif
183 } 183 }
184 184
185 //-------------------------------------------- 185 //------------------------------------------------------------------------
186 void G4ErrorFreeTrajState::Init() 186 void G4ErrorFreeTrajState::Init()
187 { 187 {
188 theTSType = G4eTS_FREE; 188 theTSType = G4eTS_FREE;
189 BuildCharge(); 189 BuildCharge();
190 theTransfMat = G4ErrorMatrix(5, 5, 0); 190 theTransfMat = G4ErrorMatrix(5, 5, 0);
191 theFirstStep = true; 191 theFirstStep = true;
192 } 192 }
193 193
194 //-------------------------------------------- 194 //------------------------------------------------------------------------
195 void G4ErrorFreeTrajState::Dump(std::ostream& 195 void G4ErrorFreeTrajState::Dump(std::ostream& out) const { out << *this; }
196 196
197 //-------------------------------------------- 197 //------------------------------------------------------------------------
198 G4int G4ErrorFreeTrajState::Update(const G4Tra 198 G4int G4ErrorFreeTrajState::Update(const G4Track* aTrack)
199 { 199 {
200 G4int ierr = 0; 200 G4int ierr = 0;
201 fTrajParam.Update(aTrack); 201 fTrajParam.Update(aTrack);
202 UpdatePosMom(aTrack->GetPosition(), aTrack-> 202 UpdatePosMom(aTrack->GetPosition(), aTrack->GetMomentum());
203 return ierr; 203 return ierr;
204 } 204 }
205 205
206 //-------------------------------------------- 206 //------------------------------------------------------------------------
207 std::ostream& operator<<(std::ostream& out, co 207 std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts)
208 { 208 {
209 std::ios::fmtflags orig_flags = out.flags(); 209 std::ios::fmtflags orig_flags = out.flags();
210 210
211 out.setf(std::ios::fixed, std::ios::floatfie 211 out.setf(std::ios::fixed, std::ios::floatfield);
212 212
213 ts.DumpPosMomError(out); 213 ts.DumpPosMomError(out);
214 214
215 out << " G4ErrorFreeTrajState: Params: " << 215 out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
216 216
217 out.flags(orig_flags); 217 out.flags(orig_flags);
218 218
219 return out; 219 return out;
220 } 220 }
221 221
222 //-------------------------------------------- 222 //------------------------------------------------------------------------
223 G4int G4ErrorFreeTrajState::PropagateError(con 223 G4int G4ErrorFreeTrajState::PropagateError(const G4Track* aTrack)
224 { 224 {
225 G4double stepLengthCm = aTrack->GetStep()->G 225 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
226 if(G4ErrorPropagatorData::GetErrorPropagator 226 if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() ==
227 G4ErrorStage_Deflation) 227 G4ErrorStage_Deflation)
228 stepLengthCm *= -1.; 228 stepLengthCm *= -1.;
229 229
230 G4double kCarTolerance = 230 G4double kCarTolerance =
231 G4GeometryTolerance::GetInstance()->GetSur 231 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
232 232
233 if(std::fabs(stepLengthCm) <= kCarTolerance 233 if(std::fabs(stepLengthCm) <= kCarTolerance / cm)
234 return 0; 234 return 0;
235 235
236 #ifdef G4EVERBOSE 236 #ifdef G4EVERBOSE
237 if(iverbose >= 2) 237 if(iverbose >= 2)
238 G4cout << " G4ErrorFreeTrajState::Propaga 238 G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl;
239 G4cout << "G4EP: iverbose=" << iverbose << G 239 G4cout << "G4EP: iverbose=" << iverbose << G4endl;
240 #endif 240 #endif
241 241
242 // * *** ERROR PROPAGATION ON A HELIX ASSUMI 242 // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
243 G4Point3D vposPost = aTrack->GetPosition() / 243 G4Point3D vposPost = aTrack->GetPosition() / cm;
244 G4Vector3D vpPost = aTrack->GetMomentum() / 244 G4Vector3D vpPost = aTrack->GetMomentum() / GeV;
245 // G4Point3D vposPre = fPosition/cm; 245 // G4Point3D vposPre = fPosition/cm;
246 // G4Vector3D vpPre = fMomentum/GeV; 246 // G4Vector3D vpPre = fMomentum/GeV;
247 G4Point3D vposPre = aTrack->GetStep()->GetPr 247 G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition() / cm;
248 G4Vector3D vpPre = aTrack->GetStep()->GetPr 248 G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum() / GeV;
249 // correct to avoid propagation along Z 249 // correct to avoid propagation along Z
250 if(vpPre.mag() == vpPre.z()) 250 if(vpPre.mag() == vpPre.z())
251 vpPre.setX(1.E-6 * MeV); 251 vpPre.setX(1.E-6 * MeV);
252 if(vpPost.mag() == vpPost.z()) 252 if(vpPost.mag() == vpPost.z())
253 vpPost.setX(1.E-6 * MeV); 253 vpPost.setX(1.E-6 * MeV);
254 254
255 G4double pPre = vpPre.mag(); 255 G4double pPre = vpPre.mag();
256 G4double pPost = vpPost.mag(); 256 G4double pPost = vpPost.mag();
257 #ifdef G4EVERBOSE 257 #ifdef G4EVERBOSE
258 if(iverbose >= 2) 258 if(iverbose >= 2)
259 { 259 {
260 G4cout << "G4EP: vposPre " << vposPre << G 260 G4cout << "G4EP: vposPre " << vposPre << G4endl << "G4EP: vposPost "
261 << vposPost << G4endl; 261 << vposPost << G4endl;
262 G4cout << "G4EP: vpPre " << vpPre << G4end 262 G4cout << "G4EP: vpPre " << vpPre << G4endl << "G4EP: vpPost " << vpPost
263 << G4endl; 263 << G4endl;
264 G4cout << " err start step " << fError << 264 G4cout << " err start step " << fError << G4endl;
265 G4cout << "G4EP: stepLengthCm " << stepLen 265 G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
266 } 266 }
267 #endif 267 #endif
268 268
269 if(pPre == 0. || pPost == 0) 269 if(pPre == 0. || pPost == 0)
270 return 2; 270 return 2;
271 G4double pInvPre = 1. / pPre; 271 G4double pInvPre = 1. / pPre;
272 G4double pInvPost = 1. / pPost; 272 G4double pInvPost = 1. / pPost;
273 G4double deltaPInv = pInvPost - pInvPre; 273 G4double deltaPInv = pInvPost - pInvPre;
274 if(iverbose >= 2) 274 if(iverbose >= 2)
275 G4cout << "G4EP: pInvPre" << pInvPre << " 275 G4cout << "G4EP: pInvPre" << pInvPre << " pInvPost:" << pInvPost
276 << " deltaPInv:" << deltaPInv << G 276 << " deltaPInv:" << deltaPInv << G4endl;
277 277
278 G4Vector3D vpPreNorm = vpPre * pInvPre; 278 G4Vector3D vpPreNorm = vpPre * pInvPre;
279 G4Vector3D vpPostNorm = vpPost * pInvPost; 279 G4Vector3D vpPostNorm = vpPost * pInvPost;
280 if(iverbose >= 2) 280 if(iverbose >= 2)
281 G4cout << "G4EP: vpPreNorm " << vpPreNorm 281 G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm
282 << G4endl; 282 << G4endl;
283 // return if propagation along Z?? 283 // return if propagation along Z??
284 if(1. - std::fabs(vpPreNorm.z()) < kCarToler 284 if(1. - std::fabs(vpPreNorm.z()) < kCarTolerance)
285 return 4; 285 return 4;
286 if(1. - std::fabs(vpPostNorm.z()) < kCarTole 286 if(1. - std::fabs(vpPostNorm.z()) < kCarTolerance)
287 return 4; 287 return 4;
288 G4double sinpPre = 288 G4double sinpPre =
289 std::sin(vpPreNorm.theta()); // cosine pe 289 std::sin(vpPreNorm.theta()); // cosine perpendicular to pPre = sine pPre
290 G4double sinpPost = 290 G4double sinpPost =
291 std::sin(vpPostNorm.theta()); // cosine p 291 std::sin(vpPostNorm.theta()); // cosine perpendicular to pPost = sine pPost
292 G4double sinpPostInv = 1. / std::sin(vpPostN 292 G4double sinpPostInv = 1. / std::sin(vpPostNorm.theta());
293 293
294 #ifdef G4EVERBOSE 294 #ifdef G4EVERBOSE
295 if(iverbose >= 2) 295 if(iverbose >= 2)
296 G4cout << "G4EP: cosl " << sinpPre << " co 296 G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
297 #endif 297 #endif
298 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN 298 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
299 //* *** NEUTRAL PARTICLE OR FIELDFREE REGION 299 //* *** NEUTRAL PARTICLE OR FIELDFREE REGION
300 G4ErrorMatrix transf(5, 5, 0); 300 G4ErrorMatrix transf(5, 5, 0);
301 301
302 transf[3][2] = stepLengthCm * sinpPost; 302 transf[3][2] = stepLengthCm * sinpPost;
303 transf[4][1] = stepLengthCm; 303 transf[4][1] = stepLengthCm;
304 for(auto ii = 0; ii < 5; ++ii) << 304 for(size_t ii = 0; ii < 5; ii++)
305 { 305 {
306 transf[ii][ii] = 1.; 306 transf[ii][ii] = 1.;
307 } 307 }
308 #ifdef G4EVERBOSE 308 #ifdef G4EVERBOSE
309 if(iverbose >= 2) 309 if(iverbose >= 2)
310 { 310 {
311 G4cout << "G4EP: transf matrix neutral " < 311 G4cout << "G4EP: transf matrix neutral " << transf;
312 } 312 }
313 #endif 313 #endif
314 314
315 // charge X propagation direction 315 // charge X propagation direction
316 G4double charge = aTrack->GetDynamicParticle 316 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
317 if(G4ErrorPropagatorData::GetErrorPropagator 317 if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() ==
318 G4ErrorMode_PropBackwards) 318 G4ErrorMode_PropBackwards)
319 { 319 {
320 charge *= -1.; 320 charge *= -1.;
321 } 321 }
322 // G4cout << " charge " << charge << G4endl 322 // G4cout << " charge " << charge << G4endl;
323 // t check if particle has charge 323 // t check if particle has charge
324 // t if( charge == 0 ) goto 45; 324 // t if( charge == 0 ) goto 45;
325 // check if the magnetic field is = 0. 325 // check if the magnetic field is = 0.
326 326
327 // position is from geant4, it is assumed to 327 // position is from geant4, it is assumed to be in mm (for debugging,
328 // eventually it will not be transformed) it 328 // eventually it will not be transformed) it is assumed vposPre[] is in cm and
329 // pos1[] is in mm. 329 // pos1[] is in mm.
330 G4double pos1[3]; 330 G4double pos1[3];
331 pos1[0] = vposPre.x() * cm; 331 pos1[0] = vposPre.x() * cm;
332 pos1[1] = vposPre.y() * cm; 332 pos1[1] = vposPre.y() * cm;
333 pos1[2] = vposPre.z() * cm; 333 pos1[2] = vposPre.z() * cm;
334 G4double pos2[3]; 334 G4double pos2[3];
335 pos2[0] = vposPost.x() * cm; 335 pos2[0] = vposPost.x() * cm;
336 pos2[1] = vposPost.y() * cm; 336 pos2[1] = vposPost.y() * cm;
337 pos2[2] = vposPost.z() * cm; 337 pos2[2] = vposPost.z() * cm;
338 G4double h1[3], h2[3]; 338 G4double h1[3], h2[3];
339 339
340 const G4Field* field = G4TransportationManag 340 const G4Field* field = G4TransportationManager::GetTransportationManager()
341 ->GetFieldManager() 341 ->GetFieldManager()
342 ->GetDetectorField( 342 ->GetDetectorField();
343 if(!field) 343 if(!field)
344 return 0; // goto 45 344 return 0; // goto 45
345 345
346 // calculate transformation except it NEUTRA 346 // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
347 if(charge != 0. && field) 347 if(charge != 0. && field)
348 { 348 {
349 field->GetFieldValue(pos1, 349 field->GetFieldValue(pos1,
350 h1); // here pos1[], 350 h1); // here pos1[], pos2[] are in mm, not changed
351 field->GetFieldValue(pos2, h2); 351 field->GetFieldValue(pos2, h2);
352 G4ThreeVector HPre = 352 G4ThreeVector HPre =
353 G4ThreeVector(h1[0], h1[1], h1[2]) / tes 353 G4ThreeVector(h1[0], h1[1], h1[2]) / tesla *
354 10.; // 10. is to get same dimensions a 354 10.; // 10. is to get same dimensions as GEANT3 (kilogauss)
355 G4ThreeVector HPost = G4ThreeVector(h2[0], 355 G4ThreeVector HPost = G4ThreeVector(h2[0], h2[1], h2[2]) / tesla * 10.;
356 G4double magHPre = HPre.mag(); 356 G4double magHPre = HPre.mag();
357 G4double magHPost = HPost.mag(); 357 G4double magHPost = HPost.mag();
358 #ifdef G4EVERBOSE 358 #ifdef G4EVERBOSE
359 if(iverbose >= 2) 359 if(iverbose >= 2)
360 { 360 {
361 G4cout << "G4EP: h1 = " << h1[0] << ", " 361 G4cout << "G4EP: h1 = " << h1[0] << ", " << h1[1] << ", " << h1[2]
362 << G4endl; 362 << G4endl;
363 G4cout << "G4EP: pos1/mm = " << pos1[0] 363 G4cout << "G4EP: pos1/mm = " << pos1[0] << ", " << pos1[1] << ", "
364 << pos1[2] << G4endl; 364 << pos1[2] << G4endl;
365 G4cout << "G4EP: pos2/mm = " << pos2[0] 365 G4cout << "G4EP: pos2/mm = " << pos2[0] << ", " << pos2[1] << ", "
366 << pos2[2] << G4endl; 366 << pos2[2] << G4endl;
367 G4cout << "G4EP: B-filed in KGauss HPre 367 G4cout << "G4EP: B-filed in KGauss HPre " << HPre << G4endl
368 << "G4EP: in KGauss HPost " << HP 368 << "G4EP: in KGauss HPost " << HPost << G4endl;
369 } 369 }
370 #endif 370 #endif
371 371
372 if(magHPre + magHPost != 0.) 372 if(magHPre + magHPost != 0.)
373 { 373 {
374 //* *** CHECK WHETHER H*ALFA/P IS TOO DI 374 //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
375 G4double gam; 375 G4double gam;
376 if(magHPost != 0.) 376 if(magHPost != 0.)
377 { 377 {
378 gam = HPost * vpPostNorm / magHPost; 378 gam = HPost * vpPostNorm / magHPost;
379 } 379 }
380 else 380 else
381 { 381 {
382 gam = HPre * vpPreNorm / magHPre; 382 gam = HPre * vpPreNorm / magHPre;
383 } 383 }
384 384
385 // G4eMagneticLimitsProcess will limit t 385 // G4eMagneticLimitsProcess will limit the step, but based on an straight
386 // line trajectory 386 // line trajectory
387 G4double alphaSqr = 1. - gam * gam; 387 G4double alphaSqr = 1. - gam * gam;
388 G4double diffHSqr = (HPre * pInvPre - H 388 G4double diffHSqr = (HPre * pInvPre - HPost * pInvPost).mag2();
389 G4double delhp6Sqr = 300. * 300.; 389 G4double delhp6Sqr = 300. * 300.;
390 #ifdef G4EVERBOSE 390 #ifdef G4EVERBOSE
391 if(iverbose >= 2) 391 if(iverbose >= 2)
392 { 392 {
393 G4cout << " G4EP: gam " << gam << " al 393 G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr
394 << " diffHSqr " << diffHSqr << 394 << " diffHSqr " << diffHSqr << G4endl;
395 G4cout << " alpha= " << std::sqrt(alph 395 G4cout << " alpha= " << std::sqrt(alphaSqr) << G4endl;
396 } 396 }
397 #endif 397 #endif
398 if(diffHSqr * alphaSqr > delhp6Sqr) 398 if(diffHSqr * alphaSqr > delhp6Sqr)
399 return 3; 399 return 3;
400 400
401 //* *** DEFINE AVERAGE MAGNETIC FIELD AN 401 //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
402 G4double pInvAver = 1. / (pInvPre + pInv 402 G4double pInvAver = 1. / (pInvPre + pInvPost);
403 G4double CFACT8 = 2.997925E-4; 403 G4double CFACT8 = 2.997925E-4;
404 // G4double HAver 404 // G4double HAver
405 G4ThreeVector vHAverNorm((HPre * pInvPre 405 G4ThreeVector vHAverNorm((HPre * pInvPre + HPost * pInvPost) * pInvAver *
406 charge * CFACT8 406 charge * CFACT8);
407 G4double HAver = vHAverNorm.mag(); 407 G4double HAver = vHAverNorm.mag();
408 G4double invHAver = 1. / HAver; 408 G4double invHAver = 1. / HAver;
409 vHAverNorm *= invHAver; 409 vHAverNorm *= invHAver;
410 #ifdef G4EVERBOSE 410 #ifdef G4EVERBOSE
411 if(iverbose >= 2) 411 if(iverbose >= 2)
412 G4cout << " G4EP: HaverNorm " << vHAve 412 G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver
413 << " charge " << charge << G4en 413 << " charge " << charge << G4endl;
414 #endif 414 #endif
415 415
416 G4double pAver = (pPre + pPost) * 416 G4double pAver = (pPre + pPost) * 0.5;
417 G4double QAver = -HAver / pAver; 417 G4double QAver = -HAver / pAver;
418 G4double thetaAver = QAver * stepLeng 418 G4double thetaAver = QAver * stepLengthCm;
419 G4double sinThetaAver = std::sin(thetaAv 419 G4double sinThetaAver = std::sin(thetaAver);
420 G4double cosThetaAver = std::cos(thetaAv 420 G4double cosThetaAver = std::cos(thetaAver);
421 G4double gamma = vHAverNorm * vpP 421 G4double gamma = vHAverNorm * vpPostNorm;
422 G4ThreeVector AN2 = vHAverNorm.cross 422 G4ThreeVector AN2 = vHAverNorm.cross(vpPostNorm);
423 423
424 #ifdef G4EVERBOSE 424 #ifdef G4EVERBOSE
425 if(iverbose >= 2) 425 if(iverbose >= 2)
426 G4cout << " G4EP: AN2 " << AN2 << " g 426 G4cout << " G4EP: AN2 " << AN2 << " gamma:" << gamma
427 << " theta=" << thetaAver << G 427 << " theta=" << thetaAver << G4endl;
428 #endif 428 #endif
429 G4double AU = 1. / vpPreNorm.perp(); 429 G4double AU = 1. / vpPreNorm.perp();
430 // t G4ThreeVector vU( vpPreNorm.cross( 430 // t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
431 G4ThreeVector vUPre(-AU * vpPreNorm.y(), 431 G4ThreeVector vUPre(-AU * vpPreNorm.y(), AU * vpPreNorm.x(), 0.);
432 G4ThreeVector vVPre(-vpPreNorm.z() * vUP 432 G4ThreeVector vVPre(-vpPreNorm.z() * vUPre.y(), vpPreNorm.z() * vUPre.x(),
433 vpPreNorm.x() * vUPr 433 vpPreNorm.x() * vUPre.y() -
434 vpPreNorm.y() * vU 434 vpPreNorm.y() * vUPre.x());
435 435
436 // 436 //
437 AU = 1. / vpPostNorm.perp(); 437 AU = 1. / vpPostNorm.perp();
438 // t G4ThreeVector vU( vpPostNorm.cross 438 // t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU
439 // ); 439 // );
440 G4ThreeVector vUPost(-AU * vpPostNorm.y( 440 G4ThreeVector vUPost(-AU * vpPostNorm.y(), AU * vpPostNorm.x(), 0.);
441 G4ThreeVector vVPost( 441 G4ThreeVector vVPost(
442 -vpPostNorm.z() * vUPost.y(), vpPostNo 442 -vpPostNorm.z() * vUPost.y(), vpPostNorm.z() * vUPost.x(),
443 vpPostNorm.x() * vUPost.y() - vpPostNo 443 vpPostNorm.x() * vUPost.y() - vpPostNorm.y() * vUPost.x());
444 #ifdef G4EVERBOSE 444 #ifdef G4EVERBOSE
445 G4cout << " vpPostNorm " << vpPostNorm < 445 G4cout << " vpPostNorm " << vpPostNorm << G4endl;
446 if(iverbose >= 2) 446 if(iverbose >= 2)
447 G4cout << " G4EP: AU " << AU << " vUPr 447 G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre
448 << " vUPost " << vUPost << " vV 448 << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
449 #endif 449 #endif
450 G4Point3D deltaPos(vposPre - vposPost); 450 G4Point3D deltaPos(vposPre - vposPost);
451 451
452 // * *** COMPLETE TRANSFORMATION MATRIX 452 // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
453 // * *** FIELD GRADIENT PERPENDICULAR TO 453 // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
454 // * *** TAKEN INTO ACCOUNT 454 // * *** TAKEN INTO ACCOUNT
455 455
456 G4double QP = QAver * pAver; // = -HAve 456 G4double QP = QAver * pAver; // = -HAver
457 #ifdef G4EVERBOSE 457 #ifdef G4EVERBOSE
458 if(iverbose >= 2) 458 if(iverbose >= 2)
459 G4cout << " G4EP: QP " << QP << " QAve 459 G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver
460 << G4endl; 460 << G4endl;
461 #endif 461 #endif
462 G4double ANV = 462 G4double ANV =
463 -(vHAverNorm.x() * vUPost.x() + vHAver 463 -(vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y());
464 G4double ANU = 464 G4double ANU =
465 (vHAverNorm.x() * vVPost.x() + vHAverN 465 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
466 vHAverNorm.z() * vVPost.z()); 466 vHAverNorm.z() * vVPost.z());
467 G4double OMcosThetaAver = 1. - cosThetaA 467 G4double OMcosThetaAver = 1. - cosThetaAver;
468 #ifdef G4EVERBOSE 468 #ifdef G4EVERBOSE
469 if(iverbose >= 2) 469 if(iverbose >= 2)
470 G4cout << "G4EP: OMcosThetaAver " << O 470 G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver "
471 << cosThetaAver << " thetaAver 471 << cosThetaAver << " thetaAver " << thetaAver << " QAver "
472 << QAver << " stepLengthCm " << 472 << QAver << " stepLengthCm " << stepLengthCm << G4endl;
473 #endif 473 #endif
474 G4double TMSINT = thetaAver - sinThetaAv 474 G4double TMSINT = thetaAver - sinThetaAver;
475 #ifdef G4EVERBOSE 475 #ifdef G4EVERBOSE
476 if(iverbose >= 2) 476 if(iverbose >= 2)
477 G4cout << " G4EP: ANV " << ANV << " AN 477 G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
478 #endif 478 #endif
479 479
480 G4ThreeVector vHUPre( 480 G4ThreeVector vHUPre(
481 -vHAverNorm.z() * vUPre.y(), vHAverNor 481 -vHAverNorm.z() * vUPre.y(), vHAverNorm.z() * vUPre.x(),
482 vHAverNorm.x() * vUPre.y() - vHAverNor 482 vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x());
483 #ifdef G4EVERBOSE 483 #ifdef G4EVERBOSE
484 // if( iverbose >= 2 ) G4cout << "G4E 484 // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " "
485 // << vHAverNorm.z() << " " << vUPre. 485 // << vHAverNorm.z() << " " << vUPre.y() << G4endl;
486 #endif 486 #endif
487 G4ThreeVector vHVPre( 487 G4ThreeVector vHVPre(
488 vHAverNorm.y() * vVPre.z() - vHAverNor 488 vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
489 vHAverNorm.z() * vVPre.x() - vHAverNor 489 vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
490 vHAverNorm.x() * vVPre.y() - vHAverNor 490 vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x());
491 #ifdef G4EVERBOSE 491 #ifdef G4EVERBOSE
492 if(iverbose >= 2) 492 if(iverbose >= 2)
493 G4cout << " G4EP: HUPre " << vHUPre << 493 G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
494 #endif 494 #endif
495 495
496 //------------------- COMPUTE MATRIX 496 //------------------- COMPUTE MATRIX
497 //---------- 1/P 497 //---------- 1/P
498 498
499 transf[0][0] = 499 transf[0][0] =
500 1. - 500 1. -
501 deltaPInv * pAver * 501 deltaPInv * pAver *
502 (1. + (vpPostNorm.x() * deltaPos.x() 502 (1. + (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
503 vpPostNorm.z() * deltaPos.z() 503 vpPostNorm.z() * deltaPos.z()) /
504 stepLengthCm) + 504 stepLengthCm) +
505 2. * deltaPInv * pAver; 505 2. * deltaPInv * pAver;
506 506
507 transf[0][1] = 507 transf[0][1] =
508 -deltaPInv / thetaAver * 508 -deltaPInv / thetaAver *
509 (TMSINT * gamma * 509 (TMSINT * gamma *
510 (vHAverNorm.x() * vVPre.x() + vHAve 510 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
511 vHAverNorm.z() * vVPre.z()) + 511 vHAverNorm.z() * vVPre.z()) +
512 sinThetaAver * 512 sinThetaAver *
513 (vVPre.x() * vpPostNorm.x() + vVPre 513 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
514 vVPre.z() * vpPostNorm.z()) + 514 vVPre.z() * vpPostNorm.z()) +
515 OMcosThetaAver * 515 OMcosThetaAver *
516 (vHVPre.x() * vpPostNorm.x() + vHVP 516 (vHVPre.x() * vpPostNorm.x() + vHVPre.y() * vpPostNorm.y() +
517 vHVPre.z() * vpPostNorm.z())); 517 vHVPre.z() * vpPostNorm.z()));
518 518
519 transf[0][2] = 519 transf[0][2] =
520 -sinpPre * deltaPInv / thetaAver * 520 -sinpPre * deltaPInv / thetaAver *
521 (TMSINT * gamma * 521 (TMSINT * gamma *
522 (vHAverNorm.x() * vUPre.x() + vHAve 522 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) +
523 sinThetaAver * 523 sinThetaAver *
524 (vUPre.x() * vpPostNorm.x() + vUPre 524 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
525 OMcosThetaAver * 525 OMcosThetaAver *
526 (vHUPre.x() * vpPostNorm.x() + vHUP 526 (vHUPre.x() * vpPostNorm.x() + vHUPre.y() * vpPostNorm.y() +
527 vHUPre.z() * vpPostNorm.z())); 527 vHUPre.z() * vpPostNorm.z()));
528 528
529 transf[0][3] = -deltaPInv / stepLengthCm 529 transf[0][3] = -deltaPInv / stepLengthCm *
530 (vUPre.x() * vpPostNorm.x 530 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y());
531 531
532 transf[0][4] = -deltaPInv / stepLengthCm 532 transf[0][4] = -deltaPInv / stepLengthCm *
533 (vVPre.x() * vpPostNorm.x 533 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
534 vVPre.z() * vpPostNorm.z 534 vVPre.z() * vpPostNorm.z());
535 535
536 // *** Lambda 536 // *** Lambda
537 transf[1][0] = 537 transf[1][0] =
538 -QP * ANV * 538 -QP * ANV *
539 (vpPostNorm.x() * deltaPos.x() + vpPos 539 (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
540 vpPostNorm.z() * deltaPos.z()) * 540 vpPostNorm.z() * deltaPos.z()) *
541 (1. + deltaPInv * pAver); 541 (1. + deltaPInv * pAver);
542 #ifdef G4EVERBOSE 542 #ifdef G4EVERBOSE
543 if(iverbose >= 3) 543 if(iverbose >= 3)
544 G4cout << "ctransf10= " << transf[1][0 544 G4cout << "ctransf10= " << transf[1][0] << " " << -QP << " " << ANV
545 << " " << vpPostNorm.x() << " " 545 << " " << vpPostNorm.x() << " " << deltaPos.x() << " "
546 << vpPostNorm.y() << " " << del 546 << vpPostNorm.y() << " " << deltaPos.y() << " " << vpPostNorm.z()
547 << " " << deltaPos.z() << " " < 547 << " " << deltaPos.z() << " " << deltaPInv << " " << pAver
548 << G4endl; 548 << G4endl;
549 #endif 549 #endif
550 550
551 transf[1][1] = 551 transf[1][1] =
552 cosThetaAver * (vVPre.x() * vVPost.x() 552 cosThetaAver * (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
553 vVPre.z() * vVPost.z() 553 vVPre.z() * vVPost.z()) +
554 sinThetaAver * (vHVPre.x() * vVPost.x( 554 sinThetaAver * (vHVPre.x() * vVPost.x() + vHVPre.y() * vVPost.y() +
555 vHVPre.z() * vVPost.z( 555 vHVPre.z() * vVPost.z()) +
556 OMcosThetaAver * 556 OMcosThetaAver *
557 (vHAverNorm.x() * vVPre.x() + vHAver 557 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
558 vHAverNorm.z() * vVPre.z()) * 558 vHAverNorm.z() * vVPre.z()) *
559 (vHAverNorm.x() * vVPost.x() + vHAve 559 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
560 vHAverNorm.z() * vVPost.z()) + 560 vHAverNorm.z() * vVPost.z()) +
561 ANV * (-sinThetaAver * 561 ANV * (-sinThetaAver *
562 (vVPre.x() * vpPostNorm.x() + 562 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
563 vVPre.z() * vpPostNorm.z()) 563 vVPre.z() * vpPostNorm.z()) +
564 OMcosThetaAver * (vVPre.x() * A 564 OMcosThetaAver * (vVPre.x() * AN2.x() + vVPre.y() * AN2.y() +
565 vVPre.z() * A 565 vVPre.z() * AN2.z()) -
566 TMSINT * gamma * 566 TMSINT * gamma *
567 (vHAverNorm.x() * vVPre.x() + 567 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
568 vHAverNorm.z() * vVPre.z())) 568 vHAverNorm.z() * vVPre.z()));
569 569
570 transf[1][2] = 570 transf[1][2] =
571 cosThetaAver * (vUPre.x() * vVPost.x() 571 cosThetaAver * (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y()) +
572 sinThetaAver * (vHUPre.x() * vVPost.x( 572 sinThetaAver * (vHUPre.x() * vVPost.x() + vHUPre.y() * vVPost.y() +
573 vHUPre.z() * vVPost.z( 573 vHUPre.z() * vVPost.z()) +
574 OMcosThetaAver * 574 OMcosThetaAver *
575 (vHAverNorm.x() * vUPre.x() + vHAver 575 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) *
576 (vHAverNorm.x() * vVPost.x() + vHAve 576 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
577 vHAverNorm.z() * vVPost.z()) + 577 vHAverNorm.z() * vVPost.z()) +
578 ANV * (-sinThetaAver * 578 ANV * (-sinThetaAver *
579 (vUPre.x() * vpPostNorm.x() + 579 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
580 OMcosThetaAver * (vUPre.x() * A 580 OMcosThetaAver * (vUPre.x() * AN2.x() + vUPre.y() * AN2.y()) -
581 TMSINT * gamma * 581 TMSINT * gamma *
582 (vHAverNorm.x() * vUPre.x() + 582 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()));
583 transf[1][2] = sinpPre * transf[1][2]; 583 transf[1][2] = sinpPre * transf[1][2];
584 584
585 transf[1][3] = -QAver * ANV * 585 transf[1][3] = -QAver * ANV *
586 (vUPre.x() * vpPostNorm.x 586 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y());
587 587
588 transf[1][4] = -QAver * ANV * 588 transf[1][4] = -QAver * ANV *
589 (vVPre.x() * vpPostNorm.x 589 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
590 vVPre.z() * vpPostNorm.z 590 vVPre.z() * vpPostNorm.z());
591 591
592 // *** Phi 592 // *** Phi
593 593
594 transf[2][0] = 594 transf[2][0] =
595 -QP * ANU * 595 -QP * ANU *
596 (vpPostNorm.x() * deltaPos.x() + vpPos 596 (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
597 vpPostNorm.z() * deltaPos.z()) * 597 vpPostNorm.z() * deltaPos.z()) *
598 sinpPostInv * (1. + deltaPInv * pAver) 598 sinpPostInv * (1. + deltaPInv * pAver);
599 #ifdef G4EVERBOSE 599 #ifdef G4EVERBOSE
600 if(iverbose >= 3) 600 if(iverbose >= 3)
601 G4cout << "ctransf20= " << transf[2][0 601 G4cout << "ctransf20= " << transf[2][0] << " " << -QP << " " << ANU
602 << " " << vpPostNorm.x() << " " 602 << " " << vpPostNorm.x() << " " << deltaPos.x() << " "
603 << vpPostNorm.y() << " " << del 603 << vpPostNorm.y() << " " << deltaPos.y() << " " << vpPostNorm.z()
604 << " " << deltaPos.z() << " " < 604 << " " << deltaPos.z() << " " << sinpPostInv << " " << deltaPInv
605 << " " << pAver << G4endl; 605 << " " << pAver << G4endl;
606 #endif 606 #endif
607 transf[2][1] = 607 transf[2][1] =
608 cosThetaAver * (vVPre.x() * vUPost.x() 608 cosThetaAver * (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y()) +
609 sinThetaAver * (vHVPre.x() * vUPost.x( 609 sinThetaAver * (vHVPre.x() * vUPost.x() + vHVPre.y() * vUPost.y()) +
610 OMcosThetaAver * 610 OMcosThetaAver *
611 (vHAverNorm.x() * vVPre.x() + vHAver 611 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
612 vHAverNorm.z() * vVPre.z()) * 612 vHAverNorm.z() * vVPre.z()) *
613 (vHAverNorm.x() * vUPost.x() + vHAve 613 (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) +
614 ANU * (-sinThetaAver * 614 ANU * (-sinThetaAver *
615 (vVPre.x() * vpPostNorm.x() + 615 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
616 vVPre.z() * vpPostNorm.z()) 616 vVPre.z() * vpPostNorm.z()) +
617 OMcosThetaAver * (vVPre.x() * A 617 OMcosThetaAver * (vVPre.x() * AN2.x() + vVPre.y() * AN2.y() +
618 vVPre.z() * A 618 vVPre.z() * AN2.z()) -
619 TMSINT * gamma * 619 TMSINT * gamma *
620 (vHAverNorm.x() * vVPre.x() + 620 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
621 vHAverNorm.z() * vVPre.z())) 621 vHAverNorm.z() * vVPre.z()));
622 transf[2][1] = sinpPostInv * transf[2][1 622 transf[2][1] = sinpPostInv * transf[2][1];
623 623
624 transf[2][2] = 624 transf[2][2] =
625 cosThetaAver * (vUPre.x() * vUPost.x() 625 cosThetaAver * (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y()) +
626 sinThetaAver * (vHUPre.x() * vUPost.x( 626 sinThetaAver * (vHUPre.x() * vUPost.x() + vHUPre.y() * vUPost.y()) +
627 OMcosThetaAver * 627 OMcosThetaAver *
628 (vHAverNorm.x() * vUPre.x() + vHAver 628 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) *
629 (vHAverNorm.x() * vUPost.x() + vHAve 629 (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) +
630 ANU * (-sinThetaAver * 630 ANU * (-sinThetaAver *
631 (vUPre.x() * vpPostNorm.x() + 631 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
632 OMcosThetaAver * (vUPre.x() * A 632 OMcosThetaAver * (vUPre.x() * AN2.x() + vUPre.y() * AN2.y()) -
633 TMSINT * gamma * 633 TMSINT * gamma *
634 (vHAverNorm.x() * vUPre.x() + 634 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()));
635 transf[2][2] = sinpPostInv * sinpPre * t 635 transf[2][2] = sinpPostInv * sinpPre * transf[2][2];
636 636
637 transf[2][3] = -QAver * ANU * 637 transf[2][3] = -QAver * ANU *
638 (vUPre.x() * vpPostNorm.x 638 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) *
639 sinpPostInv; 639 sinpPostInv;
640 #ifdef G4EVERBOSE 640 #ifdef G4EVERBOSE
641 if(iverbose >= 3) 641 if(iverbose >= 3)
642 G4cout << "ctransf23= " << transf[2][3 642 G4cout << "ctransf23= " << transf[2][3] << " " << -QAver << " " << ANU
643 << " " << vUPre.x() << " " << v 643 << " " << vUPre.x() << " " << vpPostNorm.x() << " " << vUPre.y()
644 << " " << vpPostNorm.y() << " " 644 << " " << vpPostNorm.y() << " " << sinpPostInv << G4endl;
645 #endif 645 #endif
646 646
647 transf[2][4] = -QAver * ANU * 647 transf[2][4] = -QAver * ANU *
648 (vVPre.x() * vpPostNorm.x 648 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
649 vVPre.z() * vpPostNorm.z 649 vVPre.z() * vpPostNorm.z()) *
650 sinpPostInv; 650 sinpPostInv;
651 651
652 // *** Yt 652 // *** Yt
653 653
654 transf[3][0] = pAver * 654 transf[3][0] = pAver *
655 (vUPost.x() * deltaPos.x( 655 (vUPost.x() * deltaPos.x() + vUPost.y() * deltaPos.y()) *
656 (1. + deltaPInv * pAver); 656 (1. + deltaPInv * pAver);
657 #ifdef G4EVERBOSE 657 #ifdef G4EVERBOSE
658 if(iverbose >= 3) 658 if(iverbose >= 3)
659 G4cout << "ctransf30= " << transf[3][0 659 G4cout << "ctransf30= " << transf[3][0] << " " << pAver << " "
660 << vUPost.x() << " " << deltaPo 660 << vUPost.x() << " " << deltaPos.x() << " " << vUPost.y() << " "
661 << deltaPos.y() << " " << delta 661 << deltaPos.y() << " " << deltaPInv << " " << pAver << G4endl;
662 #endif 662 #endif
663 663
664 transf[3][1] = 664 transf[3][1] =
665 (sinThetaAver * (vVPre.x() * vUPost.x( 665 (sinThetaAver * (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y()) +
666 OMcosThetaAver * (vHVPre.x() * vUPost 666 OMcosThetaAver * (vHVPre.x() * vUPost.x() + vHVPre.y() * vUPost.y()) +
667 TMSINT * (vHAverNorm.x() * vUPost.x() 667 TMSINT * (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) *
668 (vHAverNorm.x() * vVPre.x() + vHAve 668 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
669 vHAverNorm.z() * vVPre.z())) / 669 vHAverNorm.z() * vVPre.z())) /
670 QAver; 670 QAver;
671 671
672 transf[3][2] = 672 transf[3][2] =
673 (sinThetaAver * (vUPre.x() * vUPost.x( 673 (sinThetaAver * (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y()) +
674 OMcosThetaAver * (vHUPre.x() * vUPost 674 OMcosThetaAver * (vHUPre.x() * vUPost.x() + vHUPre.y() * vUPost.y()) +
675 TMSINT * (vHAverNorm.x() * vUPost.x() 675 TMSINT * (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) *
676 (vHAverNorm.x() * vUPre.x() + vHAve 676 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y())) *
677 sinpPre / QAver; 677 sinpPre / QAver;
678 #ifdef G4EVERBOSE 678 #ifdef G4EVERBOSE
679 if(iverbose >= 3) 679 if(iverbose >= 3)
680 G4cout << "ctransf32= " << transf[3][2 680 G4cout << "ctransf32= " << transf[3][2] << " " << sinThetaAver << " "
681 << vUPre.x() << " " << vUPost.x 681 << vUPre.x() << " " << vUPost.x() << " " << vUPre.y() << " "
682 << vUPost.y() << " " << OMcosTh 682 << vUPost.y() << " " << OMcosThetaAver << " " << vHUPre.x()
683 << " " << vUPost.x() << " " << 683 << " " << vUPost.x() << " " << vHUPre.y() << " " << vUPost.y()
684 << " " << TMSINT << " " << vHAv 684 << " " << TMSINT << " " << vHAverNorm.x() << " " << vUPost.x()
685 << " " << vHAverNorm.y() << " " 685 << " " << vHAverNorm.y() << " " << vUPost.y() << " "
686 << vHAverNorm.x() << " " << vUP 686 << vHAverNorm.x() << " " << vUPre.x() << " " << vHAverNorm.y()
687 << " " << vUPre.y() << " " << s 687 << " " << vUPre.y() << " " << sinpPre << " " << QAver << G4endl;
688 #endif 688 #endif
689 689
690 transf[3][3] = (vUPre.x() * vUPost.x() + 690 transf[3][3] = (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y());
691 691
692 transf[3][4] = (vVPre.x() * vUPost.x() + 692 transf[3][4] = (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y());
693 693
694 // *** Zt 694 // *** Zt
695 transf[4][0] = pAver * 695 transf[4][0] = pAver *
696 (vVPost.x() * deltaPos.x( 696 (vVPost.x() * deltaPos.x() + vVPost.y() * deltaPos.y() +
697 vVPost.z() * deltaPos.z( 697 vVPost.z() * deltaPos.z()) *
698 (1. + deltaPInv * pAver); 698 (1. + deltaPInv * pAver);
699 699
700 transf[4][1] = 700 transf[4][1] =
701 (sinThetaAver * (vVPre.x() * vVPost.x( 701 (sinThetaAver * (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
702 vVPre.z() * vVPost.z( 702 vVPre.z() * vVPost.z()) +
703 OMcosThetaAver * (vHVPre.x() * vVPost 703 OMcosThetaAver * (vHVPre.x() * vVPost.x() + vHVPre.y() * vVPost.y() +
704 vHVPre.z() * vVPost 704 vHVPre.z() * vVPost.z()) +
705 TMSINT * 705 TMSINT *
706 (vHAverNorm.x() * vVPost.x() + vHAv 706 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
707 vHAverNorm.z() * vVPost.z()) * 707 vHAverNorm.z() * vVPost.z()) *
708 (vHAverNorm.x() * vVPre.x() + vHAve 708 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
709 vHAverNorm.z() * vVPre.z())) / 709 vHAverNorm.z() * vVPre.z())) /
710 QAver; 710 QAver;
711 #ifdef G4EVERBOSE 711 #ifdef G4EVERBOSE
712 if(iverbose >= 3) 712 if(iverbose >= 3)
713 G4cout << "ctransf41= " << transf[4][1 713 G4cout << "ctransf41= " << transf[4][1] << " " << sinThetaAver << " "
714 << OMcosThetaAver << " " << TMS 714 << OMcosThetaAver << " " << TMSINT << " " << vVPre << " "
715 << vVPost << " " << vHVPre << " 715 << vVPost << " " << vHVPre << " " << vHAverNorm << " " << QAver
716 << G4endl; 716 << G4endl;
717 #endif 717 #endif
718 718
719 transf[4][2] = 719 transf[4][2] =
720 (sinThetaAver * (vUPre.x() * vVPost.x( 720 (sinThetaAver * (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y()) +
721 OMcosThetaAver * (vHUPre.x() * vVPost 721 OMcosThetaAver * (vHUPre.x() * vVPost.x() + vHUPre.y() * vVPost.y() +
722 vHUPre.z() * vVPost 722 vHUPre.z() * vVPost.z()) +
723 TMSINT * 723 TMSINT *
724 (vHAverNorm.x() * vVPost.x() + vHAv 724 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
725 vHAverNorm.z() * vVPost.z()) * 725 vHAverNorm.z() * vVPost.z()) *
726 (vHAverNorm.x() * vUPre.x() + vHAve 726 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y())) *
727 sinpPre / QAver; 727 sinpPre / QAver;
728 728
729 transf[4][3] = (vUPre.x() * vVPost.x() + 729 transf[4][3] = (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y());
730 730
731 transf[4][4] = (vVPre.x() * vVPost.x() + 731 transf[4][4] = (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
732 vVPre.z() * vVPost.z()); 732 vVPre.z() * vVPost.z());
733 // if(iverbose >= 3) G4cout <<"ctransf 733 // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<<
734 // vVPre.x() <<" "<<vVPost.x() <<" "< 734 // vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<"
735 // "<< vVPre.z() <<" "<< vVPost.z() << 735 // "<< vVPre.z() <<" "<< vVPost.z() << G4endl;
736 736
737 #ifdef G4EVERBOSE 737 #ifdef G4EVERBOSE
738 if(iverbose >= 1) 738 if(iverbose >= 1)
739 G4cout << "G4EP: transf matrix compute 739 G4cout << "G4EP: transf matrix computed " << transf << G4endl;
740 #endif 740 #endif
741 /* for( G4int ii=0;ii<5;ii++){ 741 /* for( G4int ii=0;ii<5;ii++){
742 for( G4int jj=0;jj<5;jj++){ 742 for( G4int jj=0;jj<5;jj++){
743 G4cout << transf[ii][jj] << " "; 743 G4cout << transf[ii][jj] << " ";
744 } 744 }
745 G4cout << G4endl; 745 G4cout << G4endl;
746 } */ 746 } */
747 } 747 }
748 } 748 }
749 // end of calculate transformation except it 749 // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE
750 // REGION 750 // REGION
751 /* if( iverbose >= 1 ) G4cout << "G4EP: tra 751 /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized "
752 << theFirstStep << G4endl; if( theFirstStep 752 << theFirstStep << G4endl; if( theFirstStep ) { theTransfMat = transf;
753 theFirstStep = false; 753 theFirstStep = false;
754 }else{ 754 }else{
755 theTransfMat = theTransfMat * transf; 755 theTransfMat = theTransfMat * transf;
756 if( iverbose >= 1 ) G4cout << "G4EP: trans 756 if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" <<
757 theTransfMat << G4endl; 757 theTransfMat << G4endl;
758 } 758 }
759 */ 759 */
760 theTransfMat = std::move(transf); << 760 theTransfMat = transf;
761 #ifdef G4EVERBOSE 761 #ifdef G4EVERBOSE
762 if(iverbose >= 1) 762 if(iverbose >= 1)
763 G4cout << "G4EP: error matrix before trans 763 G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
764 if(iverbose >= 2) 764 if(iverbose >= 2)
765 G4cout << " tf * err " << theTransfMat * f 765 G4cout << " tf * err " << theTransfMat * fError << G4endl
766 << " transf matrix " << theTransfMa 766 << " transf matrix " << theTransfMat.T() << G4endl;
767 #endif 767 #endif
768 768
769 fError = fError.similarity(theTransfMat).T() 769 fError = fError.similarity(theTransfMat).T();
770 //- fError = transf * fError * transf.T() 770 //- fError = transf * fError * transf.T();
771 #ifdef G4EVERBOSE 771 #ifdef G4EVERBOSE
772 if(iverbose >= 1) 772 if(iverbose >= 1)
773 G4cout << "G4EP: error matrix propagated " 773 G4cout << "G4EP: error matrix propagated " << fError << G4endl;
774 #endif 774 #endif
775 775
776 //? S = B*S*BT S.similarity(B) 776 //? S = B*S*BT S.similarity(B)
777 //? R = S 777 //? R = S
778 // not needed * *** TRANSFORM ERROR MATRIX F 778 // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL
779 // VARIABLES; 779 // VARIABLES;
780 780
781 PropagateErrorMSC(aTrack); 781 PropagateErrorMSC(aTrack);
782 782
783 PropagateErrorIoni(aTrack); 783 PropagateErrorIoni(aTrack);
784 784
785 return 0; 785 return 0;
786 } 786 }
787 787
788 //-------------------------------------------- 788 //------------------------------------------------------------------------
789 G4int G4ErrorFreeTrajState::PropagateErrorMSC( 789 G4int G4ErrorFreeTrajState::PropagateErrorMSC(const G4Track* aTrack)
790 { 790 {
791 G4ThreeVector vpPre = aTrack->GetMomentum( 791 G4ThreeVector vpPre = aTrack->GetMomentum() / GeV;
792 G4double pPre = vpPre.mag(); 792 G4double pPre = vpPre.mag();
793 G4double pBeta = pPre * pPre / (aTrac 793 G4double pBeta = pPre * pPre / (aTrack->GetTotalEnergy() / GeV);
794 G4double stepLengthCm = aTrack->GetStep()->G 794 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
795 795
796 G4Material* mate = aTrack->GetVolume()->GetL 796 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
797 G4double effZ, effA; 797 G4double effZ, effA;
798 CalculateEffectiveZandA(mate, effZ, effA); 798 CalculateEffectiveZandA(mate, effZ, effA);
799 799
800 #ifdef G4EVERBOSE 800 #ifdef G4EVERBOSE
801 if(iverbose >= 4) 801 if(iverbose >= 4)
802 G4cout << "material " 802 G4cout << "material "
803 << mate->GetName() 803 << mate->GetName()
804 //<< " " << mate->GetZ() << " " << 804 //<< " " << mate->GetZ() << " " << mate->GetA()
805 << " effZ:" << effZ << " effA:" << 805 << " effZ:" << effZ << " effA:" << effA
806 << " dens(g/mole):" << mate->GetDen 806 << " dens(g/mole):" << mate->GetDensity() / g * mole
807 << " Radlen/cm:" << mate->GetRadlen 807 << " Radlen/cm:" << mate->GetRadlen() / cm << " nuclLen/cm"
808 << mate->GetNuclearInterLength() / 808 << mate->GetNuclearInterLength() / cm << G4endl;
809 #endif 809 #endif
810 810
811 G4double RI = stepLengthCm / (mate->GetRadle 811 G4double RI = stepLengthCm / (mate->GetRadlen() / cm);
812 #ifdef G4EVERBOSE 812 #ifdef G4EVERBOSE
813 if(iverbose >= 4) 813 if(iverbose >= 4)
814 G4cout << std::setprecision(6) << std::set 814 G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI=X/X0 " << RI
815 << " stepLengthCm " << stepLengthCm 815 << " stepLengthCm " << stepLengthCm << " radlen/cm "
816 << (mate->GetRadlen() / cm) << " RI 816 << (mate->GetRadlen() / cm) << " RI*1.e10:" << RI * 1.e10 << G4endl;
817 #endif 817 #endif
818 G4double charge = aTrack->GetDynamicParticle 818 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
819 G4double DD = 1.8496E-4 * RI * (charge / 819 G4double DD = 1.8496E-4 * RI * (charge / pBeta * charge / pBeta);
820 #ifdef G4EVERBOSE 820 #ifdef G4EVERBOSE
821 if(iverbose >= 3) 821 if(iverbose >= 3)
822 G4cout << "G4EP:MSC: D*1E6= " << DD * 1.E6 822 G4cout << "G4EP:MSC: D*1E6= " << DD * 1.E6 << " pBeta " << pBeta << G4endl;
823 #endif 823 #endif
824 G4double S1 = DD * stepLengthCm * stepLength 824 G4double S1 = DD * stepLengthCm * stepLengthCm / 3.;
825 G4double S2 = DD; 825 G4double S2 = DD;
826 G4double S3 = DD * stepLengthCm / 2.; 826 G4double S3 = DD * stepLengthCm / 2.;
827 827
828 G4double CLA = 828 G4double CLA =
829 std::sqrt(vpPre.x() * vpPre.x() + vpPre.y( 829 std::sqrt(vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y()) / pPre;
830 #ifdef G4EVERBOSE 830 #ifdef G4EVERBOSE
831 if(iverbose >= 2) 831 if(iverbose >= 2)
832 G4cout << std::setw(6) << "G4EP:MSC: RI " 832 G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 "
833 << S2 << " S3 " << S3 << " CLA " << 833 << S2 << " S3 " << S3 << " CLA " << CLA << G4endl;
834 #endif 834 #endif
835 fError[1][1] += S2; 835 fError[1][1] += S2;
836 fError[1][4] -= S3; 836 fError[1][4] -= S3;
837 fError[2][2] += S2 / CLA / CLA; 837 fError[2][2] += S2 / CLA / CLA;
838 fError[2][3] += S3 / CLA; 838 fError[2][3] += S3 / CLA;
839 fError[3][3] += S1; 839 fError[3][3] += S1;
840 fError[4][4] += S1; 840 fError[4][4] += S1;
841 841
842 #ifdef G4EVERBOSE 842 #ifdef G4EVERBOSE
843 if(iverbose >= 2) 843 if(iverbose >= 2)
844 G4cout << "G4EP:MSC: error matrix propagat 844 G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl;
845 #endif 845 #endif
846 846
847 return 0; 847 return 0;
848 } 848 }
849 849
850 //-------------------------------------------- 850 //------------------------------------------------------------------------
851 void G4ErrorFreeTrajState::CalculateEffectiveZ 851 void G4ErrorFreeTrajState::CalculateEffectiveZandA(const G4Material* mate,
852 852 G4double& effZ,
853 853 G4double& effA)
854 { 854 {
855 effZ = 0.; 855 effZ = 0.;
856 effA = 0.; 856 effA = 0.;
857 auto nelem = mate->GetNumberOfElements(); << 857 G4int ii, nelem = mate->GetNumberOfElements();
858 const G4double* fracVec = mate->GetFractionV 858 const G4double* fracVec = mate->GetFractionVector();
859 for(G4int ii = 0; ii < (G4int)nelem; ++ii) << 859 for(ii = 0; ii < nelem; ii++)
860 { 860 {
861 effZ += mate->GetElement(ii)->GetZ() * fra 861 effZ += mate->GetElement(ii)->GetZ() * fracVec[ii];
862 effA += mate->GetElement(ii)->GetA() * fra 862 effA += mate->GetElement(ii)->GetA() * fracVec[ii] / g * mole;
863 } 863 }
864 } 864 }
865 865
866 //-------------------------------------------- 866 //------------------------------------------------------------------------
867 G4int G4ErrorFreeTrajState::PropagateErrorIoni 867 G4int G4ErrorFreeTrajState::PropagateErrorIoni(const G4Track* aTrack)
868 { 868 {
869 G4double stepLengthCm = aTrack->GetStep()->G 869 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
870 #ifdef G4EVERBOSE 870 #ifdef G4EVERBOSE
871 G4double DEDX2; 871 G4double DEDX2;
872 if(stepLengthCm < 1.E-7) 872 if(stepLengthCm < 1.E-7)
873 { 873 {
874 DEDX2 = 0.; 874 DEDX2 = 0.;
875 } 875 }
876 #endif 876 #endif
877 // * Calculate xi factor (KeV). 877 // * Calculate xi factor (KeV).
878 G4Material* mate = aTrack->GetVolume()->GetL 878 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
879 G4double effZ, effA; 879 G4double effZ, effA;
880 CalculateEffectiveZandA(mate, effZ, effA); 880 CalculateEffectiveZandA(mate, effZ, effA);
881 881
882 G4double Etot = aTrack->GetTotalEnergy() / 882 G4double Etot = aTrack->GetTotalEnergy() / GeV;
883 G4double beta = aTrack->GetMomentum().mag() 883 G4double beta = aTrack->GetMomentum().mag() / GeV / Etot;
884 G4double mass = aTrack->GetDynamicParticle( 884 G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV;
885 G4double gamma = Etot / mass; 885 G4double gamma = Etot / mass;
886 886
887 // * Calculate xi factor (keV). 887 // * Calculate xi factor (keV).
888 G4double XI = 153.5 * effZ * stepLengthCm * 888 G4double XI = 153.5 * effZ * stepLengthCm * (mate->GetDensity() / mg * mole) /
889 (effA * beta * beta); 889 (effA * beta * beta);
890 890
891 #ifdef G4EVERBOSE 891 #ifdef G4EVERBOSE
892 if(iverbose >= 2) 892 if(iverbose >= 2)
893 { 893 {
894 G4cout << "G4EP:IONI: XI/keV " << XI << " 894 G4cout << "G4EP:IONI: XI/keV " << XI << " beta " << beta << " gamma "
895 << gamma << G4endl; 895 << gamma << G4endl;
896 G4cout << " density " << (mate->GetDensity 896 G4cout << " density " << (mate->GetDensity() / mg * mole) << " effA "
897 << effA << " step " << stepLengthCm 897 << effA << " step " << stepLengthCm << G4endl;
898 } 898 }
899 #endif 899 #endif
900 // * Maximum energy transfer to atomic e 900 // * Maximum energy transfer to atomic electron (KeV).
901 G4double eta = beta * gamma; 901 G4double eta = beta * gamma;
902 G4double etasq = eta * eta; 902 G4double etasq = eta * eta;
903 G4double eMass = 0.51099906 / GeV; 903 G4double eMass = 0.51099906 / GeV;
904 G4double massRatio = eMass / mass; 904 G4double massRatio = eMass / mass;
905 G4double F1 = 2 * eMass * etasq; 905 G4double F1 = 2 * eMass * etasq;
906 G4double F2 = 1. + 2. * massRatio * g 906 G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio;
907 G4double Emax = 1.E+6 * F1 / F2; // no 907 G4double Emax = 1.E+6 * F1 / F2; // now in keV
908 908
909 // * *** and now sigma**2 in GeV 909 // * *** and now sigma**2 in GeV
910 G4double dedxSq = 910 G4double dedxSq =
911 XI * Emax * (1. - (beta * beta / 2.)) * 1. 911 XI * Emax * (1. - (beta * beta / 2.)) * 1.E-12; // now in GeV^2
912 /*The above formula for var(1/p) good for d 912 /*The above formula for var(1/p) good for dens scatterers. However, for MIPS
913 passing through a gas it leads to overesti 913 passing through a gas it leads to overestimation. Further more for incident
914 electrons the Emax is almost equal to inci 914 electrons the Emax is almost equal to incident energy. This leads to
915 k=Xi/Emax as small as e-6 and gradually 915 k=Xi/Emax as small as e-6 and gradually the cov matrix explodes.
916 916
917 http://www2.pv.infn.it/~rotondi/kalman_1.p 917 http://www2.pv.infn.it/~rotondi/kalman_1.pdf
918 918
919 Since I do not have enough info at the mom 919 Since I do not have enough info at the moment to implement Landau &
920 sub-Landau models for k=Xi/Emax <0.01 I'll 920 sub-Landau models for k=Xi/Emax <0.01 I'll saturate k at this value for now
921 */ 921 */
922 922
923 if(XI / Emax < 0.01) 923 if(XI / Emax < 0.01)
924 dedxSq *= 924 dedxSq *=
925 XI / Emax * 100; // Quench for low Elos 925 XI / Emax * 100; // Quench for low Elos, see above: newVar=odVar *k/0.01
926 926
927 #ifdef G4EVERBOSE 927 #ifdef G4EVERBOSE
928 if(iverbose >= 2) 928 if(iverbose >= 2)
929 G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << d 929 G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << dedxSq << " emass/GeV: " << eMass
930 << " Emax/keV: " << Emax << " k=Xi 930 << " Emax/keV: " << Emax << " k=Xi/Emax=" << XI / Emax << G4endl;
931 931
932 #endif 932 #endif
933 933
934 G4double pPre6 = 934 G4double pPre6 =
935 (aTrack->GetStep()->GetPreStepPoint()->Get 935 (aTrack->GetStep()->GetPreStepPoint()->GetMomentum() / GeV).mag();
936 pPre6 = std::pow(pPre6, 6); 936 pPre6 = std::pow(pPre6, 6);
937 // Apply it to error 937 // Apply it to error
938 fError[0][0] += Etot * Etot * dedxSq / pPre6 938 fError[0][0] += Etot * Etot * dedxSq / pPre6;
939 #ifdef G4EVERBOSE 939 #ifdef G4EVERBOSE
940 if(iverbose >= 2) 940 if(iverbose >= 2)
941 G4cout << "G4:IONI Etot/GeV: " << Etot << 941 G4cout << "G4:IONI Etot/GeV: " << Etot << " err_dedx^2/GeV^2: " << dedxSq
942 << " p^6: " << pPre6 << G4endl; 942 << " p^6: " << pPre6 << G4endl;
943 if(iverbose >= 2) 943 if(iverbose >= 2)
944 G4cout << "G4EP:IONI: error2_from_ionisati 944 G4cout << "G4EP:IONI: error2_from_ionisation "
945 << (Etot * Etot * dedxSq) / pPre6 < 945 << (Etot * Etot * dedxSq) / pPre6 << G4endl;
946 #endif 946 #endif
947 947
948 return 0; 948 return 0;
949 } 949 }
950 950