source: trunk/source/error_propagation/src/G4ErrorSurfaceTrajState.cc@ 889

Last change on this file since 889 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 8.3 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4ErrorSurfaceTrajState.cc,v 1.6 2007/06/21 15:04:08 gunter Exp $
27// GEANT4 tag $Name: HEAD $
28//
29// ------------------------------------------------------------
30// GEANT 4 class implementation file
31// ------------------------------------------------------------
32//
33
34#include "G4ErrorSurfaceTrajState.hh"
35#include "G4ErrorPropagatorData.hh"
36
37#include "G4Field.hh"
38#include "G4FieldManager.hh"
39#include "G4TransportationManager.hh"
40
41#include "G4ErrorMatrix.hh"
42
43#include <iomanip>
44
45
46//------------------------------------------------------------------------
47G4ErrorSurfaceTrajState::
48G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
49 const G4Vector3D& mom, const G4Vector3D& vecU,
50 const G4Vector3D& vecV, const G4ErrorTrajErr& errmat)
51 : G4ErrorTrajState( partType, pos, mom, errmat )
52{
53 Init();
54 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, vecU, vecV );
55}
56
57
58//------------------------------------------------------------------------
59G4ErrorSurfaceTrajState::
60G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
61 const G4Vector3D& mom, const G4Plane3D& plane,
62 const G4ErrorTrajErr& errmat )
63 : G4ErrorTrajState( partType, pos, mom, errmat )
64{
65 Init();
66 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, plane );
67
68}
69
70
71//------------------------------------------------------------------------
72G4ErrorSurfaceTrajState::
73G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Plane3D& plane )
74 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
75 tpSC.GetMomentum() )
76{
77 // fParticleType = tpSC.GetParticleType();
78 // fPosition = tpSC.GetPosition();
79 // fMomentum = tpSC.GetMomentum();
80 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, plane );
81 Init();
82
83 //----- Get the error matrix in SC coordinates
84 BuildErrorMatrix( tpSC, GetVectorV(), GetVectorW() );
85}
86
87
88//------------------------------------------------------------------------
89G4ErrorSurfaceTrajState::
90G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Vector3D& vecU,
91 const G4Vector3D& vecV )
92 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
93 tpSC.GetMomentum() )
94{
95 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, vecU, vecV );
96 theTSType = G4eTS_OS;
97 //----- Get the error matrix in SC coordinates
98 BuildErrorMatrix( tpSC, vecU, vecV );
99}
100
101
102//------------------------------------------------------------------------
103void G4ErrorSurfaceTrajState::
104BuildErrorMatrix( G4ErrorFreeTrajState& tpSC, const G4Vector3D&,
105 const G4Vector3D& )
106{
107 G4double sclambda = tpSC.GetParameters().GetLambda();
108 G4double scphi = tpSC.GetParameters().GetPhi();
109 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ){
110 sclambda *= -1;
111 scphi += CLHEP::pi;
112 }
113 G4double cosLambda = std::cos( sclambda );
114 G4double sinLambda = std::sin( sclambda );
115 G4double sinPhi = std::sin( scphi );
116 G4double cosPhi = std::cos( scphi );
117
118#ifdef G4EVERBOSE
119 if( iverbose >= 4) G4cout << " PM " << fMomentum.mag() << " pLambda " << sclambda << " pPhi " << scphi << G4endl;
120#endif
121
122 G4ThreeVector vTN( cosLambda*cosPhi, cosLambda*sinPhi,sinLambda );
123 G4ThreeVector vUN( -sinPhi, cosPhi, 0. );
124 G4ThreeVector vVN( -vTN.z()*vUN.y(),vTN.z()*vUN.x(), cosLambda );
125
126#ifdef G4EVERBOSE
127 if( iverbose >= 4) G4cout << " SC2SD: vTN " << vTN << " vUN " << vUN << " vVN " << vVN << G4endl;
128#endif
129 G4double UJ = vUN*GetVectorV();
130 G4double UK = vUN*GetVectorW();
131 G4double VJ = vVN*GetVectorV();
132 G4double VK = vVN*GetVectorW();
133
134
135 //--- Get transformation first
136 G4ErrorMatrix transfM(5, 5, 0 );
137 //--- Get magnetic field
138 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
139
140 G4Vector3D vectorU = GetVectorV().cross( GetVectorW() );
141
142#ifdef G4EVERBOSE
143 if( iverbose >= 4) G4cout << "vectors " << vectorU << " " << GetVectorV() << " " << GetVectorW() << G4endl;
144#endif
145 G4double T1R = 1. / ( vTN * vectorU );
146
147 if( fCharge != 0 && field ) {
148 G4double pos[3]; pos[0] = fPosition.x()*cm; pos[1] = fPosition.y()*cm; pos[2] = fPosition.z()*cm;
149 G4double Hd[3];
150 field->GetFieldValue( pos, Hd );
151 G4ThreeVector H = G4ThreeVector( Hd[0], Hd[1], Hd[2] ) / tesla *10.; //in kilogauus
152 G4double magH = H.mag();
153 G4double invP = 1./(fMomentum.mag()/GeV);
154 G4double magHM = magH * invP;
155 if( magH != 0. ) {
156 G4double magHM2 = fCharge / magH;
157 G4double Q = -magHM * c_light/ (km/ns);
158#ifdef G4EVERBOSE
159 if( iverbose >= 4) G4cout << GeV << " Q " << Q << " magHM " << magHM << " c_light/(km/ns) " << c_light/(km/ns) << G4endl;
160#endif
161
162 G4double sinz = -H*vUN * magHM2;
163 G4double cosz = H*vVN * magHM2;
164 G4double T3R = Q * std::pow(T1R,3);
165 G4double UI = vUN * vectorU;
166 G4double VI = vVN * vectorU;
167#ifdef G4EVERBOSE
168 if( iverbose >= 4) {
169 G4cout << " T1R " << T1R << " T3R " << T3R << G4endl;
170 G4cout << " UI " << UI << " VI " << VI << " vectorU " << vectorU << G4endl;
171 G4cout << " UJ " << UJ << " VJ " << VJ << G4endl;
172 G4cout << " UK " << UK << " VK " << VK << G4endl;
173 }
174#endif
175
176 transfM[1][3] = -UI*( VK*cosz-UK*sinz)*T3R;
177 transfM[1][4] = -VI*( VK*cosz-UK*sinz)*T3R;
178 transfM[2][3] = UI*( VJ*cosz-UJ*sinz)*T3R;
179 transfM[2][4] = VI*( VJ*cosz-UJ*sinz)*T3R;
180 }
181 }
182
183 G4double T2R = T1R * T1R;
184 transfM[0][0] = 1.;
185 transfM[1][1] = -UK*T2R;
186 transfM[1][2] = VK*cosLambda*T2R;
187 transfM[2][1] = UJ*T2R;
188 transfM[2][2] = -VJ*cosLambda*T2R;
189 transfM[3][3] = VK*T1R;
190 transfM[3][4] = -UK*T1R;
191 transfM[4][3] = -VJ*T1R;
192 transfM[4][4] = UJ*T1R;
193
194#ifdef G4EVERBOSE
195 if( iverbose >= 4) G4cout << " SC2SD transf matrix A= " << transfM << G4endl;
196#endif
197 fError = G4ErrorTrajErr( tpSC.GetError().similarity( transfM ) );
198
199#ifdef G4EVERBOSE
200 if( iverbose >= 1) G4cout << "G4EP: error matrix SC2SD " << fError << G4endl;
201 if( iverbose >= 4) G4cout << "G4ErrorSurfaceTrajState from SC " << *this << G4endl;
202#endif
203
204}
205
206
207//------------------------------------------------------------------------
208void G4ErrorSurfaceTrajState::Init()
209{
210 theTSType = G4eTS_OS;
211 BuildCharge();
212}
213
214
215//------------------------------------------------------------------------
216void G4ErrorSurfaceTrajState::Dump( std::ostream& out ) const
217{
218 out << *this;
219}
220
221
222//------------------------------------------------------------------------
223std::ostream& operator<<(std::ostream& out, const G4ErrorSurfaceTrajState& ts)
224{
225 out.setf(std::ios::fixed,std::ios::floatfield);
226
227 ts.DumpPosMomError( out );
228
229 out << " G4ErrorSurfaceTrajState: Params: " << ts.fTrajParam << G4endl;
230 return out;
231}
Note: See TracBrowser for help on using the repository browser.