00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "G4ErrorSurfaceTrajState.hh"
00034 #include "G4ErrorPropagatorData.hh"
00035
00036 #include "G4PhysicalConstants.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4Field.hh"
00039 #include "G4FieldManager.hh"
00040 #include "G4TransportationManager.hh"
00041
00042 #include "G4ErrorMatrix.hh"
00043
00044 #include <iomanip>
00045
00046
00047 G4ErrorSurfaceTrajState::
00048 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
00049 const G4Vector3D& mom, const G4Vector3D& vecU,
00050 const G4Vector3D& vecV, const G4ErrorTrajErr& errmat)
00051 : G4ErrorTrajState( partType, pos, mom, errmat )
00052 {
00053 Init();
00054 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, vecU, vecV );
00055 }
00056
00057
00058
00059 G4ErrorSurfaceTrajState::
00060 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos,
00061 const G4Vector3D& mom, const G4Plane3D& plane,
00062 const G4ErrorTrajErr& errmat )
00063 : G4ErrorTrajState( partType, pos, mom, errmat )
00064 {
00065 Init();
00066 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, plane );
00067
00068 }
00069
00070
00071
00072 G4ErrorSurfaceTrajState::
00073 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Plane3D& plane )
00074 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
00075 tpSC.GetMomentum() )
00076 {
00077
00078
00079
00080 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, plane );
00081 Init();
00082
00083
00084 BuildErrorMatrix( tpSC, GetVectorV(), GetVectorW() );
00085 }
00086
00087
00088
00089 G4ErrorSurfaceTrajState::
00090 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Vector3D& vecU,
00091 const G4Vector3D& vecV , G4ErrorMatrix &transfM)
00092 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(),
00093 tpSC.GetMomentum() )
00094 {
00095 Init();
00096 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, vecU, vecV );
00097
00098 transfM= BuildErrorMatrix( tpSC, vecU, vecV );
00099 }
00100
00101
00102
00103 G4ErrorMatrix G4ErrorSurfaceTrajState::
00104 BuildErrorMatrix( G4ErrorFreeTrajState& tpSC, const G4Vector3D&,
00105 const G4Vector3D& )
00106 {
00107 G4double sclambda = tpSC.GetParameters().GetLambda();
00108 G4double scphi = tpSC.GetParameters().GetPhi();
00109 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ){
00110 sclambda *= -1;
00111 scphi += CLHEP::pi;
00112 }
00113 G4double cosLambda = std::cos( sclambda );
00114 G4double sinLambda = std::sin( sclambda );
00115 G4double sinPhi = std::sin( scphi );
00116 G4double cosPhi = std::cos( scphi );
00117
00118 #ifdef G4EVERBOSE
00119 if( iverbose >= 4) G4cout << " PM " << fMomentum.mag() << " pLambda " << sclambda << " pPhi " << scphi << G4endl;
00120 #endif
00121
00122 G4ThreeVector vTN( cosLambda*cosPhi, cosLambda*sinPhi,sinLambda );
00123 G4ThreeVector vUN( -sinPhi, cosPhi, 0. );
00124 G4ThreeVector vVN( -vTN.z()*vUN.y(),vTN.z()*vUN.x(), cosLambda );
00125
00126 #ifdef G4EVERBOSE
00127 if( iverbose >= 4) G4cout << " SC2SD: vTN " << vTN << " vUN " << vUN << " vVN " << vVN << G4endl;
00128 #endif
00129 G4double UJ = vUN*GetVectorV();
00130 G4double UK = vUN*GetVectorW();
00131 G4double VJ = vVN*GetVectorV();
00132 G4double VK = vVN*GetVectorW();
00133
00134
00135
00136 G4ErrorMatrix transfM(5, 5, 0 );
00137
00138 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
00139
00140 G4Vector3D vectorU = GetVectorV().cross( GetVectorW() );
00141 G4double T1R = 1. / ( vTN * vectorU );
00142
00143 #ifdef G4EVERBOSE
00144 if( iverbose >= 4) G4cout << "surf vectors:U,V,W " << vectorU << " " << GetVectorV() << " " << GetVectorW() << " T1R:"<<T1R<<G4endl;
00145 #endif
00146
00147
00148 if( fCharge != 0 && field ) {
00149 G4double pos[3]; pos[0] = fPosition.x()*cm; pos[1] = fPosition.y()*cm; pos[2] = fPosition.z()*cm;
00150 G4double Hd[3];
00151 field->GetFieldValue( pos, Hd );
00152 G4ThreeVector H = G4ThreeVector( Hd[0], Hd[1], Hd[2] ) / tesla *10.;
00153 G4double magH = H.mag();
00154 G4double invP = 1./(fMomentum.mag()/GeV);
00155 G4double magHM = magH * invP;
00156 if( magH != 0. ) {
00157 G4double magHM2 = fCharge / magH;
00158 G4double Q = -magHM * c_light/ (km/ns);
00159 #ifdef G4EVERBOSE
00160 if( iverbose >= 4) G4cout << GeV << " Q " << Q << " magHM " << magHM << " c_light/(km/ns) " << c_light/(km/ns) << G4endl;
00161 #endif
00162
00163 G4double sinz = -H*vUN * magHM2;
00164 G4double cosz = H*vVN * magHM2;
00165 G4double T3R = Q * std::pow(T1R,3);
00166 G4double UI = vUN * vectorU;
00167 G4double VI = vVN * vectorU;
00168 #ifdef G4EVERBOSE
00169 if( iverbose >= 4) {
00170 G4cout << " T1R " << T1R << " T3R " << T3R << G4endl;
00171 G4cout << " UI " << UI << " VI " << VI << " vectorU " << vectorU << G4endl;
00172 G4cout << " UJ " << UJ << " VJ " << VJ << G4endl;
00173 G4cout << " UK " << UK << " VK " << VK << G4endl;
00174 }
00175 #endif
00176
00177 transfM[1][3] = -UI*( VK*cosz-UK*sinz)*T3R;
00178 transfM[1][4] = -VI*( VK*cosz-UK*sinz)*T3R;
00179 transfM[2][3] = UI*( VJ*cosz-UJ*sinz)*T3R;
00180 transfM[2][4] = VI*( VJ*cosz-UJ*sinz)*T3R;
00181 }
00182 }
00183
00184 G4double T2R = T1R * T1R;
00185 transfM[0][0] = 1.;
00186 transfM[1][1] = -UK*T2R;
00187 transfM[1][2] = VK*cosLambda*T2R;
00188 transfM[2][1] = UJ*T2R;
00189 transfM[2][2] = -VJ*cosLambda*T2R;
00190 transfM[3][3] = VK*T1R;
00191 transfM[3][4] = -UK*T1R;
00192 transfM[4][3] = -VJ*T1R;
00193 transfM[4][4] = UJ*T1R;
00194
00195 #ifdef G4EVERBOSE
00196 if( iverbose >= 4) G4cout << " SC2SD transf matrix A= " << transfM << G4endl;
00197 #endif
00198 fError = G4ErrorTrajErr( tpSC.GetError().similarity( transfM ) );
00199
00200 #ifdef G4EVERBOSE
00201 if( iverbose >= 1) G4cout << "G4EP: error matrix SC2SD " << fError << G4endl;
00202 if( iverbose >= 4) G4cout << "G4ErrorSurfaceTrajState from SC " << *this << G4endl;
00203 #endif
00204
00205 return transfM;
00206 }
00207
00208
00209
00210 void G4ErrorSurfaceTrajState::Init()
00211 {
00212 theTSType = G4eTS_OS;
00213 BuildCharge();
00214 }
00215
00216
00217
00218 void G4ErrorSurfaceTrajState::Dump( std::ostream& out ) const
00219 {
00220 out << *this;
00221 }
00222
00223
00224
00225 std::ostream& operator<<(std::ostream& out, const G4ErrorSurfaceTrajState& ts)
00226 {
00227 std::ios::fmtflags oldFlags = out.flags();
00228 out.setf(std::ios::fixed,std::ios::floatfield);
00229
00230 ts.DumpPosMomError( out );
00231
00232 out << " G4ErrorSurfaceTrajState: Params: " << ts.fTrajParam << G4endl;
00233 out.flags(oldFlags);
00234 return out;
00235 }