| 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: G4ErrorSurfaceTrajState.cc,v 1.6 2007-06-21 15:04:08 gunter Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
26 // 28 //
27 // ------------------------------------------- 29 // ------------------------------------------------------------
28 // GEANT 4 class implementation file << 30 // GEANT 4 class implementation file
29 // ------------------------------------------- 31 // ------------------------------------------------------------
30 // 32 //
31 33
32 #include "G4ErrorSurfaceTrajState.hh" 34 #include "G4ErrorSurfaceTrajState.hh"
33 #include "G4ErrorPropagatorData.hh" 35 #include "G4ErrorPropagatorData.hh"
34 36
35 #include "G4PhysicalConstants.hh" <<
36 #include "G4SystemOfUnits.hh" <<
37 #include "G4Field.hh" 37 #include "G4Field.hh"
38 #include "G4FieldManager.hh" 38 #include "G4FieldManager.hh"
39 #include "G4TransportationManager.hh" 39 #include "G4TransportationManager.hh"
40 40
41 #include "G4ErrorMatrix.hh" 41 #include "G4ErrorMatrix.hh"
42 42
43 #include <iomanip> 43 #include <iomanip>
44 44
>> 45
45 //-------------------------------------------- 46 //------------------------------------------------------------------------
46 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 47 G4ErrorSurfaceTrajState::
47 const G4String& partType, const G4Point3D& p << 48 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
48 const G4Vector3D& vecU, const G4Vector3D& ve << 49 const G4Vector3D& mom, const G4Vector3D& vecU,
49 : G4ErrorTrajState(partType, pos, mom, errma << 50 const G4Vector3D& vecV, const G4ErrorTrajErr& errmat)
>> 51 : G4ErrorTrajState( partType, pos, mom, errmat )
50 { 52 {
51 Init(); 53 Init();
52 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo << 54 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, vecU, vecV );
53 } 55 }
54 56
>> 57
55 //-------------------------------------------- 58 //------------------------------------------------------------------------
56 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 59 G4ErrorSurfaceTrajState::
57 << 60 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
58 << 61 const G4Vector3D& mom, const G4Plane3D& plane,
59 << 62 const G4ErrorTrajErr& errmat )
60 << 63 : G4ErrorTrajState( partType, pos, mom, errmat )
61 : G4ErrorTrajState(partType, pos, mom, errma <<
62 { 64 {
63 Init(); 65 Init();
64 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo << 66 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, plane );
>> 67
65 } 68 }
66 69
>> 70
67 //-------------------------------------------- 71 //------------------------------------------------------------------------
68 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 72 G4ErrorSurfaceTrajState::
69 << 73 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Plane3D& plane )
70 : G4ErrorTrajState(tpSC.GetParticleType(), t << 74 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
71 tpSC.GetMomentum()) << 75 tpSC.GetMomentum() )
72 { 76 {
73 // fParticleType = tpSC.GetParticleType(); 77 // fParticleType = tpSC.GetParticleType();
74 // fPosition = tpSC.GetPosition(); 78 // fPosition = tpSC.GetPosition();
75 // fMomentum = tpSC.GetMomentum(); 79 // fMomentum = tpSC.GetMomentum();
76 fTrajParam = G4ErrorSurfaceTrajParam(fPositi << 80 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, plane );
77 Init(); 81 Init();
78 82
79 //----- Get the error matrix in SC coordinat 83 //----- Get the error matrix in SC coordinates
80 BuildErrorMatrix(tpSC, GetVectorV(), GetVect << 84 BuildErrorMatrix( tpSC, GetVectorV(), GetVectorW() );
81 } 85 }
82 86
>> 87
83 //-------------------------------------------- 88 //------------------------------------------------------------------------
84 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 89 G4ErrorSurfaceTrajState::
85 << 90 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Vector3D& vecU,
86 << 91 const G4Vector3D& vecV )
87 << 92 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
88 : G4ErrorTrajState(tpSC.GetParticleType(), t << 93 tpSC.GetMomentum() )
89 tpSC.GetMomentum()) <<
90 { 94 {
91 Init(); // needed to define charge sign << 95 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, vecU, vecV );
92 fTrajParam = G4ErrorSurfaceTrajParam(fPositi << 96 theTSType = G4eTS_OS;
93 //----- Get the error matrix in SC coordinat 97 //----- Get the error matrix in SC coordinates
94 transfM = BuildErrorMatrix(tpSC, vecU, vecV) << 98 BuildErrorMatrix( tpSC, vecU, vecV );
95 } 99 }
96 100
>> 101
97 //-------------------------------------------- 102 //------------------------------------------------------------------------
98 G4ErrorMatrix G4ErrorSurfaceTrajState::BuildEr << 103 void G4ErrorSurfaceTrajState::
99 G4ErrorFreeTrajState& tpSC, const G4Vector3D << 104 BuildErrorMatrix( G4ErrorFreeTrajState& tpSC, const G4Vector3D&,
>> 105 const G4Vector3D& )
100 { 106 {
101 G4double sclambda = tpSC.GetParameters().Get 107 G4double sclambda = tpSC.GetParameters().GetLambda();
102 G4double scphi = tpSC.GetParameters().Get << 108 G4double scphi = tpSC.GetParameters().GetPhi();
103 if(G4ErrorPropagatorData::GetErrorPropagator << 109 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ){
104 G4ErrorMode_PropBackwards) <<
105 { <<
106 sclambda *= -1; 110 sclambda *= -1;
107 scphi += CLHEP::pi; 111 scphi += CLHEP::pi;
108 } 112 }
109 G4double cosLambda = std::cos(sclambda); << 113 G4double cosLambda = std::cos( sclambda );
110 G4double sinLambda = std::sin(sclambda); << 114 G4double sinLambda = std::sin( sclambda );
111 G4double sinPhi = std::sin(scphi); << 115 G4double sinPhi = std::sin( scphi );
112 G4double cosPhi = std::cos(scphi); << 116 G4double cosPhi = std::cos( scphi );
113 117
114 #ifdef G4EVERBOSE 118 #ifdef G4EVERBOSE
115 if(iverbose >= 4) << 119 if( iverbose >= 4) G4cout << " PM " << fMomentum.mag() << " pLambda " << sclambda << " pPhi " << scphi << G4endl;
116 G4cout << " PM " << fMomentum.mag() << " p <<
117 << scphi << G4endl; <<
118 #endif 120 #endif
119 121
120 G4ThreeVector vTN(cosLambda * cosPhi, cosLam << 122 G4ThreeVector vTN( cosLambda*cosPhi, cosLambda*sinPhi,sinLambda );
121 G4ThreeVector vUN(-sinPhi, cosPhi, 0.); << 123 G4ThreeVector vUN( -sinPhi, cosPhi, 0. );
122 G4ThreeVector vVN(-vTN.z() * vUN.y(), vTN.z( << 124 G4ThreeVector vVN( -vTN.z()*vUN.y(),vTN.z()*vUN.x(), cosLambda );
123 << 125
124 #ifdef G4EVERBOSE 126 #ifdef G4EVERBOSE
125 if(iverbose >= 4) << 127 if( iverbose >= 4) G4cout << " SC2SD: vTN " << vTN << " vUN " << vUN << " vVN " << vVN << G4endl;
126 G4cout << " SC2SD: vTN " << vTN << " vUN " <<
127 << G4endl; <<
128 #endif 128 #endif
129 G4double UJ = vUN * GetVectorV(); << 129 G4double UJ = vUN*GetVectorV();
130 G4double UK = vUN * GetVectorW(); << 130 G4double UK = vUN*GetVectorW();
131 G4double VJ = vVN * GetVectorV(); << 131 G4double VJ = vVN*GetVectorV();
132 G4double VK = vVN * GetVectorW(); << 132 G4double VK = vVN*GetVectorW();
>> 133
133 134
134 //--- Get transformation first 135 //--- Get transformation first
135 G4ErrorMatrix transfM(5, 5, 0); << 136 G4ErrorMatrix transfM(5, 5, 0 );
136 //--- Get magnetic field 137 //--- Get magnetic field
137 const G4Field* field = G4TransportationManag << 138 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
138 ->GetFieldManager() <<
139 ->GetDetectorField( <<
140 139
141 G4Vector3D vectorU = GetVectorV().cross(GetV << 140 G4Vector3D vectorU = GetVectorV().cross( GetVectorW() );
142 G4double T1R = 1. / (vTN * vectorU); <<
143 141
144 #ifdef G4EVERBOSE 142 #ifdef G4EVERBOSE
145 if(iverbose >= 4) << 143 if( iverbose >= 4) G4cout << "vectors " << vectorU << " " << GetVectorV() << " " << GetVectorW() << G4endl;
146 G4cout << "surf vectors:U,V,W " << vectorU <<
147 << GetVectorW() << " T1R:" << T1R <<
148 #endif 144 #endif
>> 145 G4double T1R = 1. / ( vTN * vectorU );
149 146
150 if(fCharge != 0 && field) << 147 if( fCharge != 0 && field ) {
151 { << 148 G4double pos[3]; pos[0] = fPosition.x()*cm; pos[1] = fPosition.y()*cm; pos[2] = fPosition.z()*cm;
152 G4double pos[3]; <<
153 pos[0] = fPosition.x() * cm; <<
154 pos[1] = fPosition.y() * cm; <<
155 pos[2] = fPosition.z() * cm; <<
156 G4double Hd[3]; 149 G4double Hd[3];
157 field->GetFieldValue(pos, Hd); << 150 field->GetFieldValue( pos, Hd );
158 G4ThreeVector H = << 151 G4ThreeVector H = G4ThreeVector( Hd[0], Hd[1], Hd[2] ) / tesla *10.; //in kilogauus
159 G4ThreeVector(Hd[0], Hd[1], Hd[2]) / tes << 152 G4double magH = H.mag();
160 G4double magH = H.mag(); << 153 G4double invP = 1./(fMomentum.mag()/GeV);
161 G4double invP = 1. / (fMomentum.mag() / G <<
162 G4double magHM = magH * invP; 154 G4double magHM = magH * invP;
163 if(magH != 0.) << 155 if( magH != 0. ) {
164 { <<
165 G4double magHM2 = fCharge / magH; 156 G4double magHM2 = fCharge / magH;
166 G4double Q = -magHM * c_light / (km << 157 G4double Q = -magHM * c_light/ (km/ns);
167 #ifdef G4EVERBOSE 158 #ifdef G4EVERBOSE
168 if(iverbose >= 4) << 159 if( iverbose >= 4) G4cout << GeV << " Q " << Q << " magHM " << magHM << " c_light/(km/ns) " << c_light/(km/ns) << G4endl;
169 G4cout << GeV << " Q " << Q << " magHM <<
170 << c_light / (km / ns) << G4end <<
171 #endif 160 #endif
172 161
173 G4double sinz = -H * vUN * magHM2; << 162 G4double sinz = -H*vUN * magHM2;
174 G4double cosz = H * vVN * magHM2; << 163 G4double cosz = H*vVN * magHM2;
175 G4double T3R = Q * std::pow(T1R, 3); << 164 G4double T3R = Q * std::pow(T1R,3);
176 G4double UI = vUN * vectorU; << 165 G4double UI = vUN * vectorU;
177 G4double VI = vVN * vectorU; << 166 G4double VI = vVN * vectorU;
178 #ifdef G4EVERBOSE 167 #ifdef G4EVERBOSE
179 if(iverbose >= 4) << 168 if( iverbose >= 4) {
180 { <<
181 G4cout << " T1R " << T1R << " T3R " << 169 G4cout << " T1R " << T1R << " T3R " << T3R << G4endl;
182 G4cout << " UI " << UI << " VI " << VI << 170 G4cout << " UI " << UI << " VI " << VI << " vectorU " << vectorU << G4endl;
183 << G4endl; <<
184 G4cout << " UJ " << UJ << " VJ " << VJ 171 G4cout << " UJ " << UJ << " VJ " << VJ << G4endl;
185 G4cout << " UK " << UK << " VK " << VK 172 G4cout << " UK " << UK << " VK " << VK << G4endl;
186 } 173 }
187 #endif 174 #endif
188 175
189 transfM[1][3] = -UI * (VK * cosz - UK * << 176 transfM[1][3] = -UI*( VK*cosz-UK*sinz)*T3R;
190 transfM[1][4] = -VI * (VK * cosz - UK * << 177 transfM[1][4] = -VI*( VK*cosz-UK*sinz)*T3R;
191 transfM[2][3] = UI * (VJ * cosz - UJ * s << 178 transfM[2][3] = UI*( VJ*cosz-UJ*sinz)*T3R;
192 transfM[2][4] = VI * (VJ * cosz - UJ * s << 179 transfM[2][4] = VI*( VJ*cosz-UJ*sinz)*T3R;
193 } 180 }
194 } 181 }
195 182
196 G4double T2R = T1R * T1R; << 183 G4double T2R = T1R * T1R;
197 transfM[0][0] = 1.; 184 transfM[0][0] = 1.;
198 transfM[1][1] = -UK * T2R; << 185 transfM[1][1] = -UK*T2R;
199 transfM[1][2] = VK * cosLambda * T2R; << 186 transfM[1][2] = VK*cosLambda*T2R;
200 transfM[2][1] = UJ * T2R; << 187 transfM[2][1] = UJ*T2R;
201 transfM[2][2] = -VJ * cosLambda * T2R; << 188 transfM[2][2] = -VJ*cosLambda*T2R;
202 transfM[3][3] = VK * T1R; << 189 transfM[3][3] = VK*T1R;
203 transfM[3][4] = -UK * T1R; << 190 transfM[3][4] = -UK*T1R;
204 transfM[4][3] = -VJ * T1R; << 191 transfM[4][3] = -VJ*T1R;
205 transfM[4][4] = UJ * T1R; << 192 transfM[4][4] = UJ*T1R;
206 193
207 #ifdef G4EVERBOSE 194 #ifdef G4EVERBOSE
208 if(iverbose >= 4) << 195 if( iverbose >= 4) G4cout << " SC2SD transf matrix A= " << transfM << G4endl;
209 G4cout << " SC2SD transf matrix A= " << tr <<
210 #endif 196 #endif
211 fError = G4ErrorTrajErr(tpSC.GetError().simi << 197 fError = G4ErrorTrajErr( tpSC.GetError().similarity( transfM ) );
212 198
213 #ifdef G4EVERBOSE 199 #ifdef G4EVERBOSE
214 if(iverbose >= 1) << 200 if( iverbose >= 1) G4cout << "G4EP: error matrix SC2SD " << fError << G4endl;
215 G4cout << "G4EP: error matrix SC2SD " << f << 201 if( iverbose >= 4) G4cout << "G4ErrorSurfaceTrajState from SC " << *this << G4endl;
216 if(iverbose >= 4) <<
217 G4cout << "G4ErrorSurfaceTrajState from SC <<
218 #endif 202 #endif
219 203
220 return transfM; // want to use trasnfM for <<
221 } 204 }
222 205
>> 206
223 //-------------------------------------------- 207 //------------------------------------------------------------------------
224 void G4ErrorSurfaceTrajState::Init() 208 void G4ErrorSurfaceTrajState::Init()
225 { 209 {
226 theTSType = G4eTS_OS; << 210 theTSType = G4eTS_OS;
227 BuildCharge(); << 211 BuildCharge();
228 } 212 }
229 213
>> 214
230 //-------------------------------------------- 215 //------------------------------------------------------------------------
231 void G4ErrorSurfaceTrajState::Dump(std::ostrea << 216 void G4ErrorSurfaceTrajState::Dump( std::ostream& out ) const
>> 217 {
>> 218 out << *this;
>> 219 }
>> 220
232 221
233 //-------------------------------------------- 222 //------------------------------------------------------------------------
234 std::ostream& operator<<(std::ostream& out, co 223 std::ostream& operator<<(std::ostream& out, const G4ErrorSurfaceTrajState& ts)
235 { 224 {
236 std::ios::fmtflags oldFlags = out.flags(); 225 std::ios::fmtflags oldFlags = out.flags();
237 out.setf(std::ios::fixed, std::ios::floatfie << 226 out.setf(std::ios::fixed,std::ios::floatfield);
238 << 227
239 ts.DumpPosMomError(out); << 228 ts.DumpPosMomError( out );
240 << 229
241 out << " G4ErrorSurfaceTrajState: Params: " 230 out << " G4ErrorSurfaceTrajState: Params: " << ts.fTrajParam << G4endl;
242 out.flags(oldFlags); 231 out.flags(oldFlags);
243 return out; 232 return out;
244 } 233 }
245 234