ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ErrorSurfaceTrajState.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ErrorSurfaceTrajState.cc
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 //
27 // ------------------------------------------------------------
28 // GEANT 4 class implementation file
29 // ------------------------------------------------------------
30 //
31 
32 #include "G4ErrorSurfaceTrajState.hh"
33 #include "G4ErrorPropagatorData.hh"
34 
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4Field.hh"
38 #include "G4FieldManager.hh"
39 #include "G4TransportationManager.hh"
40 
41 #include "G4ErrorMatrix.hh"
42 
43 #include <iomanip>
44 
45 //------------------------------------------------------------------------
46 G4ErrorSurfaceTrajState::
47 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
48  const G4Vector3D& mom, const G4Vector3D& vecU,
49  const G4Vector3D& vecV, const G4ErrorTrajErr& errmat)
50  : G4ErrorTrajState( partType, pos, mom, errmat )
51 {
52  Init();
53  fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, vecU, vecV );
54 }
55 
56 
57 //------------------------------------------------------------------------
58 G4ErrorSurfaceTrajState::
59 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
60  const G4Vector3D& mom, const G4Plane3D& plane,
61  const G4ErrorTrajErr& errmat )
62  : G4ErrorTrajState( partType, pos, mom, errmat )
63 {
64  Init();
65  fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, plane );
66 
67 }
68 
69 
70 //------------------------------------------------------------------------
71 G4ErrorSurfaceTrajState::
72 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Plane3D& plane )
73  : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
74  tpSC.GetMomentum() )
75 {
76  // fParticleType = tpSC.GetParticleType();
77  // fPosition = tpSC.GetPosition();
78  // fMomentum = tpSC.GetMomentum();
79  fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, plane );
80  Init();
81 
82  //----- Get the error matrix in SC coordinates
83  BuildErrorMatrix( tpSC, GetVectorV(), GetVectorW() );
84 }
85 
86 
87 //------------------------------------------------------------------------
88 G4ErrorSurfaceTrajState::
89 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Vector3D& vecU,
90  const G4Vector3D& vecV , G4ErrorMatrix &transfM)
91  : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
92  tpSC.GetMomentum() )
93 {
94  Init(); // needed to define charge sign
95  fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, vecU, vecV );
96  //----- Get the error matrix in SC coordinates
97  transfM= BuildErrorMatrix( tpSC, vecU, vecV );
98 }
99 
100 
101 //------------------------------------------------------------------------
102 G4ErrorMatrix G4ErrorSurfaceTrajState::
103 BuildErrorMatrix( G4ErrorFreeTrajState& tpSC, const G4Vector3D&,
104  const G4Vector3D& )
105 {
106  G4double sclambda = tpSC.GetParameters().GetLambda();
107  G4double scphi = tpSC.GetParameters().GetPhi();
108  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ){
109  sclambda *= -1;
110  scphi += CLHEP::pi;
111  }
112  G4double cosLambda = std::cos( sclambda );
113  G4double sinLambda = std::sin( sclambda );
114  G4double sinPhi = std::sin( scphi );
115  G4double cosPhi = std::cos( scphi );
116 
117 #ifdef G4EVERBOSE
118  if( iverbose >= 4) G4cout << " PM " << fMomentum.mag() << " pLambda " << sclambda << " pPhi " << scphi << G4endl;
119 #endif
120 
121  G4ThreeVector vTN( cosLambda*cosPhi, cosLambda*sinPhi,sinLambda );
122  G4ThreeVector vUN( -sinPhi, cosPhi, 0. );
123  G4ThreeVector vVN( -vTN.z()*vUN.y(),vTN.z()*vUN.x(), cosLambda );
124 
125 #ifdef G4EVERBOSE
126  if( iverbose >= 4) G4cout << " SC2SD: vTN " << vTN << " vUN " << vUN << " vVN " << vVN << G4endl;
127 #endif
128  G4double UJ = vUN*GetVectorV();
129  G4double UK = vUN*GetVectorW();
130  G4double VJ = vVN*GetVectorV();
131  G4double VK = vVN*GetVectorW();
132 
133 
134  //--- Get transformation first
135  G4ErrorMatrix transfM(5, 5, 0 );
136  //--- Get magnetic field
137  const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
138 
139  G4Vector3D vectorU = GetVectorV().cross( GetVectorW() );
140  G4double T1R = 1. / ( vTN * vectorU );
141 
142 #ifdef G4EVERBOSE
143  if( iverbose >= 4) G4cout << "surf vectors:U,V,W " << vectorU << " " << GetVectorV() << " " << GetVectorW() << " T1R:"<<T1R<<G4endl;
144 #endif
145 
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  return transfM; // want to use trasnfM for the reverse transform
205 }
206 
207 
208 //------------------------------------------------------------------------
209 void G4ErrorSurfaceTrajState::Init()
210 {
211  theTSType = G4eTS_OS;
212  BuildCharge();
213 }
214 
215 
216 //------------------------------------------------------------------------
217 void G4ErrorSurfaceTrajState::Dump( std::ostream& out ) const
218 {
219  out << *this;
220 }
221 
222 
223 //------------------------------------------------------------------------
224 std::ostream& operator<<(std::ostream& out, const G4ErrorSurfaceTrajState& ts)
225 {
226  std::ios::fmtflags oldFlags = out.flags();
227  out.setf(std::ios::fixed,std::ios::floatfield);
228 
229  ts.DumpPosMomError( out );
230 
231  out << " G4ErrorSurfaceTrajState: Params: " << ts.fTrajParam << G4endl;
232  out.flags(oldFlags);
233  return out;
234 }