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 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------ 28 // GEANT 4 class implementation file << 28 // GEANT 4 class implementation file 29 // ------------------------------------------- 29 // ------------------------------------------------------------ 30 // 30 // 31 31 32 #include "G4ErrorSurfaceTrajState.hh" 32 #include "G4ErrorSurfaceTrajState.hh" 33 #include "G4ErrorPropagatorData.hh" 33 #include "G4ErrorPropagatorData.hh" 34 34 35 #include "G4PhysicalConstants.hh" 35 #include "G4PhysicalConstants.hh" 36 #include "G4SystemOfUnits.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 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 46 G4ErrorSurfaceTrajState:: 47 const G4String& partType, const G4Point3D& p << 47 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos, 48 const G4Vector3D& vecU, const G4Vector3D& ve << 48 const G4Vector3D& mom, const G4Vector3D& vecU, 49 : G4ErrorTrajState(partType, pos, mom, errma << 49 const G4Vector3D& vecV, const G4ErrorTrajErr& errmat) >> 50 : G4ErrorTrajState( partType, pos, mom, errmat ) 50 { 51 { 51 Init(); 52 Init(); 52 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo << 53 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, vecU, vecV ); 53 } 54 } 54 55 >> 56 55 //-------------------------------------------- 57 //------------------------------------------------------------------------ 56 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 58 G4ErrorSurfaceTrajState:: 57 << 59 G4ErrorSurfaceTrajState( const G4String& partType, const G4Point3D& pos, 58 << 60 const G4Vector3D& mom, const G4Plane3D& plane, 59 << 61 const G4ErrorTrajErr& errmat ) 60 << 62 : G4ErrorTrajState( partType, pos, mom, errmat ) 61 : G4ErrorTrajState(partType, pos, mom, errma << 62 { 63 { 63 Init(); 64 Init(); 64 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo << 65 fTrajParam = G4ErrorSurfaceTrajParam( pos, mom, plane ); >> 66 65 } 67 } 66 68 >> 69 67 //-------------------------------------------- 70 //------------------------------------------------------------------------ 68 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 71 G4ErrorSurfaceTrajState:: 69 << 72 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Plane3D& plane ) 70 : G4ErrorTrajState(tpSC.GetParticleType(), t << 73 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(), 71 tpSC.GetMomentum()) << 74 tpSC.GetMomentum() ) 72 { 75 { 73 // fParticleType = tpSC.GetParticleType(); 76 // fParticleType = tpSC.GetParticleType(); 74 // fPosition = tpSC.GetPosition(); 77 // fPosition = tpSC.GetPosition(); 75 // fMomentum = tpSC.GetMomentum(); 78 // fMomentum = tpSC.GetMomentum(); 76 fTrajParam = G4ErrorSurfaceTrajParam(fPositi << 79 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, plane ); 77 Init(); 80 Init(); 78 81 79 //----- Get the error matrix in SC coordinat 82 //----- Get the error matrix in SC coordinates 80 BuildErrorMatrix(tpSC, GetVectorV(), GetVect << 83 BuildErrorMatrix( tpSC, GetVectorV(), GetVectorW() ); 81 } 84 } >> 85 82 86 83 //-------------------------------------------- 87 //------------------------------------------------------------------------ 84 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta << 88 G4ErrorSurfaceTrajState:: 85 << 89 G4ErrorSurfaceTrajState( G4ErrorFreeTrajState& tpSC, const G4Vector3D& vecU, 86 << 90 const G4Vector3D& vecV , G4ErrorMatrix &transfM) 87 << 91 : G4ErrorTrajState( tpSC.GetParticleType(), tpSC.GetPosition(), 88 : G4ErrorTrajState(tpSC.GetParticleType(), t << 92 tpSC.GetMomentum() ) 89 tpSC.GetMomentum()) << 90 { 93 { 91 Init(); // needed to define charge sign << 94 Init(); // needed to define charge sign 92 fTrajParam = G4ErrorSurfaceTrajParam(fPositi << 95 fTrajParam = G4ErrorSurfaceTrajParam( fPosition, fMomentum, vecU, vecV ); 93 //----- Get the error matrix in SC coordinat 96 //----- Get the error matrix in SC coordinates 94 transfM = BuildErrorMatrix(tpSC, vecU, vecV) << 97 transfM= BuildErrorMatrix( tpSC, vecU, vecV ); 95 } 98 } 96 99 >> 100 97 //-------------------------------------------- 101 //------------------------------------------------------------------------ 98 G4ErrorMatrix G4ErrorSurfaceTrajState::BuildEr << 102 G4ErrorMatrix G4ErrorSurfaceTrajState:: 99 G4ErrorFreeTrajState& tpSC, const G4Vector3D << 103 BuildErrorMatrix( G4ErrorFreeTrajState& tpSC, const G4Vector3D&, >> 104 const G4Vector3D& ) 100 { 105 { 101 G4double sclambda = tpSC.GetParameters().Get 106 G4double sclambda = tpSC.GetParameters().GetLambda(); 102 G4double scphi = tpSC.GetParameters().Get << 107 G4double scphi = tpSC.GetParameters().GetPhi(); 103 if(G4ErrorPropagatorData::GetErrorPropagator << 108 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ){ 104 G4ErrorMode_PropBackwards) << 105 { << 106 sclambda *= -1; 109 sclambda *= -1; 107 scphi += CLHEP::pi; 110 scphi += CLHEP::pi; 108 } 111 } 109 G4double cosLambda = std::cos(sclambda); << 112 G4double cosLambda = std::cos( sclambda ); 110 G4double sinLambda = std::sin(sclambda); << 113 G4double sinLambda = std::sin( sclambda ); 111 G4double sinPhi = std::sin(scphi); << 114 G4double sinPhi = std::sin( scphi ); 112 G4double cosPhi = std::cos(scphi); << 115 G4double cosPhi = std::cos( scphi ); 113 116 114 #ifdef G4EVERBOSE 117 #ifdef G4EVERBOSE 115 if(iverbose >= 4) << 118 if( iverbose >= 4) G4cout << " PM " << fMomentum.mag() << " pLambda " << sclambda << " pPhi " << scphi << G4endl; 116 G4cout << " PM " << fMomentum.mag() << " p << 117 << scphi << G4endl; << 118 #endif 119 #endif 119 120 120 G4ThreeVector vTN(cosLambda * cosPhi, cosLam << 121 G4ThreeVector vTN( cosLambda*cosPhi, cosLambda*sinPhi,sinLambda ); 121 G4ThreeVector vUN(-sinPhi, cosPhi, 0.); << 122 G4ThreeVector vUN( -sinPhi, cosPhi, 0. ); 122 G4ThreeVector vVN(-vTN.z() * vUN.y(), vTN.z( << 123 G4ThreeVector vVN( -vTN.z()*vUN.y(),vTN.z()*vUN.x(), cosLambda ); 123 << 124 124 #ifdef G4EVERBOSE 125 #ifdef G4EVERBOSE 125 if(iverbose >= 4) << 126 if( iverbose >= 4) G4cout << " SC2SD: vTN " << vTN << " vUN " << vUN << " vVN " << vVN << G4endl; 126 G4cout << " SC2SD: vTN " << vTN << " vUN " << 127 << G4endl; << 128 #endif 127 #endif 129 G4double UJ = vUN * GetVectorV(); << 128 G4double UJ = vUN*GetVectorV(); 130 G4double UK = vUN * GetVectorW(); << 129 G4double UK = vUN*GetVectorW(); 131 G4double VJ = vVN * GetVectorV(); << 130 G4double VJ = vVN*GetVectorV(); 132 G4double VK = vVN * GetVectorW(); << 131 G4double VK = vVN*GetVectorW(); >> 132 133 133 134 //--- Get transformation first 134 //--- Get transformation first 135 G4ErrorMatrix transfM(5, 5, 0); << 135 G4ErrorMatrix transfM(5, 5, 0 ); 136 //--- Get magnetic field 136 //--- Get magnetic field 137 const G4Field* field = G4TransportationManag << 137 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); 138 ->GetFieldManager() << 139 ->GetDetectorField( << 140 138 141 G4Vector3D vectorU = GetVectorV().cross(GetV << 139 G4Vector3D vectorU = GetVectorV().cross( GetVectorW() ); 142 G4double T1R = 1. / (vTN * vectorU); << 140 G4double T1R = 1. / ( vTN * vectorU ); 143 141 144 #ifdef G4EVERBOSE 142 #ifdef G4EVERBOSE 145 if(iverbose >= 4) << 143 if( iverbose >= 4) G4cout << "surf vectors:U,V,W " << vectorU << " " << GetVectorV() << " " << GetVectorW() << " T1R:"<<T1R<<G4endl; 146 G4cout << "surf vectors:U,V,W " << vectorU << 147 << GetVectorW() << " T1R:" << T1R << 148 #endif 144 #endif 149 145 150 if(fCharge != 0 && field) << 146 151 { << 147 if( fCharge != 0 && field ) { 152 G4double pos[3]; << 148 G4double pos[3]; pos[0] = fPosition.x()*cm; pos[1] = fPosition.y()*cm; pos[2] = fPosition.z()*cm; 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 << 204 return transfM; // want to use trasnfM for the reverse transform 221 } 205 } 222 206 >> 207 223 //-------------------------------------------- 208 //------------------------------------------------------------------------ 224 void G4ErrorSurfaceTrajState::Init() 209 void G4ErrorSurfaceTrajState::Init() 225 { 210 { 226 theTSType = G4eTS_OS; << 211 theTSType = G4eTS_OS; 227 BuildCharge(); << 212 BuildCharge(); 228 } 213 } 229 214 >> 215 230 //-------------------------------------------- 216 //------------------------------------------------------------------------ 231 void G4ErrorSurfaceTrajState::Dump(std::ostrea << 217 void G4ErrorSurfaceTrajState::Dump( std::ostream& out ) const >> 218 { >> 219 out << *this; >> 220 } >> 221 232 222 233 //-------------------------------------------- 223 //------------------------------------------------------------------------ 234 std::ostream& operator<<(std::ostream& out, co 224 std::ostream& operator<<(std::ostream& out, const G4ErrorSurfaceTrajState& ts) 235 { 225 { 236 std::ios::fmtflags oldFlags = out.flags(); 226 std::ios::fmtflags oldFlags = out.flags(); 237 out.setf(std::ios::fixed, std::ios::floatfie << 227 out.setf(std::ios::fixed,std::ios::floatfield); 238 << 228 239 ts.DumpPosMomError(out); << 229 ts.DumpPosMomError( out ); 240 << 230 241 out << " G4ErrorSurfaceTrajState: Params: " 231 out << " G4ErrorSurfaceTrajState: Params: " << ts.fTrajParam << G4endl; 242 out.flags(oldFlags); 232 out.flags(oldFlags); 243 return out; 233 return out; 244 } 234 } 245 235