Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 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::G4ErrorSurfaceTrajSta 47 const G4String& partType, const G4Point3D& p 48 const G4Vector3D& vecU, const G4Vector3D& ve 49 : G4ErrorTrajState(partType, pos, mom, errma 50 { 51 Init(); 52 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo 53 } 54 55 //-------------------------------------------- 56 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta 57 58 59 60 61 : G4ErrorTrajState(partType, pos, mom, errma 62 { 63 Init(); 64 fTrajParam = G4ErrorSurfaceTrajParam(pos, mo 65 } 66 67 //-------------------------------------------- 68 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta 69 70 : G4ErrorTrajState(tpSC.GetParticleType(), t 71 tpSC.GetMomentum()) 72 { 73 // fParticleType = tpSC.GetParticleType(); 74 // fPosition = tpSC.GetPosition(); 75 // fMomentum = tpSC.GetMomentum(); 76 fTrajParam = G4ErrorSurfaceTrajParam(fPositi 77 Init(); 78 79 //----- Get the error matrix in SC coordinat 80 BuildErrorMatrix(tpSC, GetVectorV(), GetVect 81 } 82 83 //-------------------------------------------- 84 G4ErrorSurfaceTrajState::G4ErrorSurfaceTrajSta 85 86 87 88 : G4ErrorTrajState(tpSC.GetParticleType(), t 89 tpSC.GetMomentum()) 90 { 91 Init(); // needed to define charge sign 92 fTrajParam = G4ErrorSurfaceTrajParam(fPositi 93 //----- Get the error matrix in SC coordinat 94 transfM = BuildErrorMatrix(tpSC, vecU, vecV) 95 } 96 97 //-------------------------------------------- 98 G4ErrorMatrix G4ErrorSurfaceTrajState::BuildEr 99 G4ErrorFreeTrajState& tpSC, const G4Vector3D 100 { 101 G4double sclambda = tpSC.GetParameters().Get 102 G4double scphi = tpSC.GetParameters().Get 103 if(G4ErrorPropagatorData::GetErrorPropagator 104 G4ErrorMode_PropBackwards) 105 { 106 sclambda *= -1; 107 scphi += CLHEP::pi; 108 } 109 G4double cosLambda = std::cos(sclambda); 110 G4double sinLambda = std::sin(sclambda); 111 G4double sinPhi = std::sin(scphi); 112 G4double cosPhi = std::cos(scphi); 113 114 #ifdef G4EVERBOSE 115 if(iverbose >= 4) 116 G4cout << " PM " << fMomentum.mag() << " p 117 << scphi << G4endl; 118 #endif 119 120 G4ThreeVector vTN(cosLambda * cosPhi, cosLam 121 G4ThreeVector vUN(-sinPhi, cosPhi, 0.); 122 G4ThreeVector vVN(-vTN.z() * vUN.y(), vTN.z( 123 124 #ifdef G4EVERBOSE 125 if(iverbose >= 4) 126 G4cout << " SC2SD: vTN " << vTN << " vUN " 127 << 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 //--- Get transformation first 135 G4ErrorMatrix transfM(5, 5, 0); 136 //--- Get magnetic field 137 const G4Field* field = G4TransportationManag 138 ->GetFieldManager() 139 ->GetDetectorField( 140 141 G4Vector3D vectorU = GetVectorV().cross(GetV 142 G4double T1R = 1. / (vTN * vectorU); 143 144 #ifdef G4EVERBOSE 145 if(iverbose >= 4) 146 G4cout << "surf vectors:U,V,W " << vectorU 147 << GetVectorW() << " T1R:" << T1R 148 #endif 149 150 if(fCharge != 0 && field) 151 { 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]; 157 field->GetFieldValue(pos, Hd); 158 G4ThreeVector H = 159 G4ThreeVector(Hd[0], Hd[1], Hd[2]) / tes 160 G4double magH = H.mag(); 161 G4double invP = 1. / (fMomentum.mag() / G 162 G4double magHM = magH * invP; 163 if(magH != 0.) 164 { 165 G4double magHM2 = fCharge / magH; 166 G4double Q = -magHM * c_light / (km 167 #ifdef G4EVERBOSE 168 if(iverbose >= 4) 169 G4cout << GeV << " Q " << Q << " magHM 170 << c_light / (km / ns) << G4end 171 #endif 172 173 G4double sinz = -H * vUN * magHM2; 174 G4double cosz = H * vVN * magHM2; 175 G4double T3R = Q * std::pow(T1R, 3); 176 G4double UI = vUN * vectorU; 177 G4double VI = vVN * vectorU; 178 #ifdef G4EVERBOSE 179 if(iverbose >= 4) 180 { 181 G4cout << " T1R " << T1R << " T3R " << 182 G4cout << " UI " << UI << " VI " << VI 183 << G4endl; 184 G4cout << " UJ " << UJ << " VJ " << VJ 185 G4cout << " UK " << UK << " VK " << VK 186 } 187 #endif 188 189 transfM[1][3] = -UI * (VK * cosz - UK * 190 transfM[1][4] = -VI * (VK * cosz - UK * 191 transfM[2][3] = UI * (VJ * cosz - UJ * s 192 transfM[2][4] = VI * (VJ * cosz - UJ * s 193 } 194 } 195 196 G4double T2R = T1R * T1R; 197 transfM[0][0] = 1.; 198 transfM[1][1] = -UK * T2R; 199 transfM[1][2] = VK * cosLambda * T2R; 200 transfM[2][1] = UJ * T2R; 201 transfM[2][2] = -VJ * cosLambda * T2R; 202 transfM[3][3] = VK * T1R; 203 transfM[3][4] = -UK * T1R; 204 transfM[4][3] = -VJ * T1R; 205 transfM[4][4] = UJ * T1R; 206 207 #ifdef G4EVERBOSE 208 if(iverbose >= 4) 209 G4cout << " SC2SD transf matrix A= " << tr 210 #endif 211 fError = G4ErrorTrajErr(tpSC.GetError().simi 212 213 #ifdef G4EVERBOSE 214 if(iverbose >= 1) 215 G4cout << "G4EP: error matrix SC2SD " << f 216 if(iverbose >= 4) 217 G4cout << "G4ErrorSurfaceTrajState from SC 218 #endif 219 220 return transfM; // want to use trasnfM for 221 } 222 223 //-------------------------------------------- 224 void G4ErrorSurfaceTrajState::Init() 225 { 226 theTSType = G4eTS_OS; 227 BuildCharge(); 228 } 229 230 //-------------------------------------------- 231 void G4ErrorSurfaceTrajState::Dump(std::ostrea 232 233 //-------------------------------------------- 234 std::ostream& operator<<(std::ostream& out, co 235 { 236 std::ios::fmtflags oldFlags = out.flags(); 237 out.setf(std::ios::fixed, std::ios::floatfie 238 239 ts.DumpPosMomError(out); 240 241 out << " G4ErrorSurfaceTrajState: Params: " 242 out.flags(oldFlags); 243 return out; 244 } 245