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