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 // $Id: G4ErrorFreeTrajState.cc,v 1.8 2009/05/14 13:53:06 arce Exp $ >> 27 // GEANT4 tag $Name: geant4-09-03 $ 26 // 28 // 27 // ------------------------------------------- 29 // ------------------------------------------------------------ 28 // GEANT 4 class implementation file << 30 // GEANT 4 class implementation file 29 // ------------------------------------------- 31 // ------------------------------------------------------------ 30 // 32 // >> 33 #include "G4ErrorFreeTrajState.hh" >> 34 #include "G4ErrorFreeTrajParam.hh" >> 35 #include "G4ErrorSurfaceTrajState.hh" 31 36 >> 37 #include "G4ErrorMatrix.hh" 32 #include <iomanip> 38 #include <iomanip> 33 39 34 #include "G4PhysicalConstants.hh" << 35 #include "G4SystemOfUnits.hh" << 36 #include "G4Field.hh" 40 #include "G4Field.hh" 37 #include "G4FieldManager.hh" 41 #include "G4FieldManager.hh" 38 #include "G4TransportationManager.hh" 42 #include "G4TransportationManager.hh" 39 #include "G4GeometryTolerance.hh" 43 #include "G4GeometryTolerance.hh" 40 #include "G4Material.hh" 44 #include "G4Material.hh" 41 #include "G4ErrorPropagatorData.hh" 45 #include "G4ErrorPropagatorData.hh" 42 #include "G4ErrorFreeTrajState.hh" << 43 #include "G4ErrorFreeTrajParam.hh" << 44 #include "G4ErrorSurfaceTrajState.hh" << 45 #include "G4ErrorMatrix.hh" << 46 46 47 //-------------------------------------------- 47 //------------------------------------------------------------------------ 48 G4ErrorFreeTrajState::G4ErrorFreeTrajState(con << 48 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partType, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partType, pos, mom, errmat ) 49 con << 50 con << 51 con << 52 : G4ErrorTrajState(partName, pos, mom, errma << 53 { 49 { 54 fTrajParam = G4ErrorFreeTrajParam(pos, mom); << 50 fTrajParam = G4ErrorFreeTrajParam( pos, mom ); 55 Init(); 51 Init(); 56 } 52 } 57 53 >> 54 58 //-------------------------------------------- 55 //------------------------------------------------------------------------ 59 G4ErrorFreeTrajState::G4ErrorFreeTrajState(con << 56 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() ) 60 : G4ErrorTrajState(tpSD.GetParticleType(), t << 61 tpSD.GetMomentum()) << 62 { 57 { 63 // G4ThreeVector planeNormal = tpSD.GetPlane << 58 // G4ThreeVector planeNormal = tpSD.GetPlaneNormal(); 64 // G4double fPt = tpSD.GetMomentum()*planeNo << 59 // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to plane 65 // plane G4ErrorSurfaceTrajParam tpSDparam = << 60 // G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters(); 66 // G4ThreeVector Psc = fPt * planeNormal + << 61 // G4ThreeVector Psc = fPt * planeNormal + tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW(); 67 // tpSDparam.GetPU()*tpSDparam.GetVectorU() << 68 62 69 fTrajParam = G4ErrorFreeTrajParam(fPosition, << 63 fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum ); 70 Init(); 64 Init(); 71 65 72 //----- Get the error matrix in SC coordinat 66 //----- Get the error matrix in SC coordinates 73 G4ErrorSurfaceTrajParam tpSDparam = tpSD.Get 67 G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters(); 74 G4double mom = fMomentu << 68 G4double mom = fMomentum.mag(); 75 G4double mom2 = fMomentu << 69 G4double mom2 = fMomentum.mag2(); 76 G4double TVW1 = << 70 G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPV()*tpSDparam.GetPV()) ); 77 std::sqrt(mom2 / (mom2 + tpSDparam.GetPV() << 71 G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 ); 78 tpSDparam.GetPW() * tpSD << 72 G4Vector3D vectorU = tpSDparam.GetVectorV().cross( tpSDparam.GetVectorW() ); 79 G4ThreeVector vTVW(TVW1, tpSDparam.GetPV() / << 73 G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW(); 80 tpSDparam.GetPW() / mom * << 74 81 G4Vector3D vectorU = tpSDparam.GetVectorV(). << 75 #ifdef G4EVERBOSE 82 G4Vector3D vTN = vTVW.x() * vectorU + vT << 76 if( iverbose >= 5){ 83 vTVW.z() * tpSDparam.GetVec << 77 G4double pc2 = std::asin( vTN.z() ); 84 << 78 G4double pc3 = std::atan (vTN.y()/vTN.x()); 85 #ifdef G4EVERBOSE << 79 86 if(iverbose >= 5) << 80 G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() << " diff " << pc2-GetParameters().GetLambda() << G4endl; 87 { << 81 G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() << " diff " << pc3-GetParameters().GetPhi() << G4endl; 88 G4double pc2 = std::asin(vTN.z()); << 82 } 89 G4double pc3 = std::atan(vTN.y() / vTN.x() << 83 #endif 90 << 84 91 G4cout << " CHECK: pc2 " << pc2 << " = " < << 85 //--- Get the unit vectors perp to P 92 << " diff " << pc2 - GetParameters( << 86 G4double cosl = std::cos( GetParameters().GetLambda() ); 93 G4cout << " CHECK: pc3 " << pc3 << " = " < << 87 if (cosl < 1.E-30) cosl = 1.E-30; 94 << " diff " << pc3 - GetParameters( << 88 G4double cosl1 = 1./cosl; 95 } << 89 G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. ); 96 #endif << 90 G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl ); 97 << 91 98 //--- Get the unit vectors perp to P << 92 G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.); 99 G4double cosl = std::cos(GetParameters().Get << 93 G4Vector3D vVperp = vUperp.cross( fMomentum ); 100 if(cosl < 1.E-30) << 94 vUperp *= 1./vUperp.mag(); 101 cosl = 1.E-30; << 95 vVperp *= 1./vVperp.mag(); 102 G4double cosl1 = 1. / cosl; << 96 103 G4Vector3D vUN(-vTN.y() * cosl1, vTN.x() * c << 97 #ifdef G4EVERBOSE 104 G4Vector3D vVN(-vTN.z() * vUN.y(), vTN.z() * << 98 if( iverbose >= 5){ 105 << 99 G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff " << (vUN-vUperp).mag() << G4endl; 106 G4Vector3D vUperp = G4Vector3D(-fMomentum.y( << 100 G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff " << (vVN-vVperp).mag() << G4endl; 107 G4Vector3D vVperp = vUperp.cross(fMomentum); << 101 } 108 vUperp *= 1. / vUperp.mag(); << 109 vVperp *= 1. / vVperp.mag(); << 110 << 111 #ifdef G4EVERBOSE << 112 if(iverbose >= 5) << 113 { << 114 G4cout << " CHECK: vUN " << vUN << " = " < << 115 << (vUN - vUperp).mag() << G4endl; << 116 G4cout << " CHECK: vVN " << vVN << " = " < << 117 << (vVN - vVperp).mag() << G4endl; << 118 } << 119 #endif 102 #endif 120 103 121 // get the dot products of vectors perpendic << 104 //get the dot products of vectors perpendicular to direction and vector defining SD plane 122 // defining SD plane << 123 G4double dUU = vUperp * tpSD.GetVectorV(); 105 G4double dUU = vUperp * tpSD.GetVectorV(); 124 G4double dUV = vUperp * tpSD.GetVectorW(); 106 G4double dUV = vUperp * tpSD.GetVectorW(); 125 G4double dVU = vVperp * tpSD.GetVectorV(); 107 G4double dVU = vVperp * tpSD.GetVectorV(); 126 G4double dVV = vVperp * tpSD.GetVectorW(); 108 G4double dVV = vVperp * tpSD.GetVectorW(); 127 109 >> 110 128 //--- Get transformation first 111 //--- Get transformation first 129 G4ErrorMatrix transfM(5, 5, 1); << 112 G4ErrorMatrix transfM(5, 5, 1 ); 130 //--- Get magnetic field 113 //--- Get magnetic field 131 const G4Field* field = G4TransportationManag << 114 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); 132 ->GetFieldManager() << 115 G4ThreeVector dir = fTrajParam.GetDirection(); 133 ->GetDetectorField( << 116 G4double invCosTheta = 1./std::cos( dir.theta() ); 134 G4ThreeVector dir = fTrajParam.GetDirecti << 117 135 G4double invCosTheta = 1. / std::cos(dir.the << 118 if( fCharge != 0 136 G4cout << " dir=" << dir << " invCosTheta " << 119 && field ) { 137 << 120 G4double pos1[3]; pos1[0] = fPosition.x()*cm; pos1[1] = fPosition.y()*cm; pos1[2] = fPosition.z()*cm; 138 if(fCharge != 0 && field) << 139 { << 140 G4double pos1[3]; << 141 pos1[0] = fPosition.x() * cm; << 142 pos1[1] = fPosition.y() * cm; << 143 pos1[2] = fPosition.z() * cm; << 144 G4double h1[3]; 121 G4double h1[3]; 145 field->GetFieldValue(pos1, h1); << 122 field->GetFieldValue( pos1, h1 ); 146 G4ThreeVector HPre = G4ThreeVector(h1[0], << 123 G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; 147 G4double magHPre = HPre.mag(); << 124 G4double magHPre = HPre.mag(); 148 G4double invP = 1. / fMomentum.mag(); << 125 G4double invP = 1./fMomentum.mag(); 149 G4double magHPreM = magHPre * invP; << 126 G4double magHPreM = magHPre * invP; 150 if(magHPre != 0.) << 127 if( magHPre != 0. ) { 151 { << 152 G4double magHPreM2 = fCharge / magHPre; 128 G4double magHPreM2 = fCharge / magHPre; 153 129 154 G4double Q = -magHPreM * c_light; << 130 G4double Q = -magHPreM * c_light; 155 G4double sinz = -HPre * vUperp * magHPre << 131 G4double sinz = -HPre*vUperp * magHPreM2; 156 G4double cosz = HPre * vVperp * magHPreM << 132 G4double cosz = HPre*vVperp * magHPreM2; 157 << 133 158 transfM[1][3] = -Q * dir.y() * sinz; << 134 transfM[1][3] = -Q*dir.y()*sinz; 159 transfM[1][4] = -Q * dir.z() * sinz; << 135 transfM[1][4] = -Q*dir.z()*sinz; 160 transfM[2][3] = -Q * dir.y() * cosz * in << 136 transfM[2][3] = -Q*dir.y()*cosz*invCosTheta; 161 transfM[2][4] = -Q * dir.z() * cosz * in << 137 transfM[2][4] = -Q*dir.z()*cosz*invCosTheta; 162 } 138 } 163 } 139 } 164 140 165 transfM[0][0] = 1.; 141 transfM[0][0] = 1.; 166 transfM[1][1] = dir.x() * dVU; << 142 transfM[1][1] = dir.x()*dVU; 167 transfM[1][2] = dir.x() * dVV; << 143 transfM[1][2] = dir.x()*dVV; 168 transfM[2][1] = dir.x() * dUU * invCosTheta; << 144 transfM[2][1] = dir.x()*dUU*invCosTheta; 169 transfM[2][2] = dir.x() * dUV * invCosTheta; << 145 transfM[2][2] = dir.x()*dUV*invCosTheta; 170 transfM[3][3] = dUU; 146 transfM[3][3] = dUU; 171 transfM[3][4] = dUV; 147 transfM[3][4] = dUV; 172 transfM[4][3] = dVU; 148 transfM[4][3] = dVU; 173 transfM[4][4] = dVV; 149 transfM[4][4] = dVV; 174 150 175 fError = G4ErrorTrajErr(tpSD.GetError().simi << 151 fError = G4ErrorTrajErr( tpSD.GetError().similarity( transfM ) ); 176 152 177 #ifdef G4EVERBOSE 153 #ifdef G4EVERBOSE 178 if(iverbose >= 1) << 154 if( iverbose >= 1) G4cout << "error matrix SD2SC " << fError << G4endl; 179 G4cout << "error matrix SD2SC " << fError << 155 if( iverbose >= 4) G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl; 180 if(iverbose >= 4) << 181 G4cout << "G4ErrorFreeTrajState from SD " << 182 #endif 156 #endif 183 } 157 } 184 158 >> 159 185 //-------------------------------------------- 160 //------------------------------------------------------------------------ 186 void G4ErrorFreeTrajState::Init() 161 void G4ErrorFreeTrajState::Init() 187 { 162 { 188 theTSType = G4eTS_FREE; 163 theTSType = G4eTS_FREE; 189 BuildCharge(); 164 BuildCharge(); 190 theTransfMat = G4ErrorMatrix(5, 5, 0); << 165 theTransfMat = G4ErrorMatrix(5,5,0); 191 theFirstStep = true; << 166 //- theFirstStep = true; 192 } 167 } 193 168 194 //-------------------------------------------- 169 //------------------------------------------------------------------------ 195 void G4ErrorFreeTrajState::Dump(std::ostream& << 170 void G4ErrorFreeTrajState::Dump( std::ostream& out ) const >> 171 { >> 172 out << *this; >> 173 } 196 174 197 //-------------------------------------------- 175 //------------------------------------------------------------------------ 198 G4int G4ErrorFreeTrajState::Update(const G4Tra << 176 G4int G4ErrorFreeTrajState::Update( const G4Track* aTrack ) 199 { 177 { 200 G4int ierr = 0; 178 G4int ierr = 0; 201 fTrajParam.Update(aTrack); << 179 fTrajParam.Update( aTrack ); 202 UpdatePosMom(aTrack->GetPosition(), aTrack-> << 180 UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() ); 203 return ierr; 181 return ierr; >> 182 204 } 183 } 205 184 >> 185 206 //-------------------------------------------- 186 //------------------------------------------------------------------------ 207 std::ostream& operator<<(std::ostream& out, co 187 std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts) 208 { 188 { 209 std::ios::fmtflags orig_flags = out.flags(); << 189 out.setf(std::ios::fixed,std::ios::floatfield); 210 << 211 out.setf(std::ios::fixed, std::ios::floatfie << 212 << 213 ts.DumpPosMomError(out); << 214 190 >> 191 >> 192 ts.DumpPosMomError( out ); >> 193 215 out << " G4ErrorFreeTrajState: Params: " << 194 out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl; 216 195 217 out.flags(orig_flags); << 218 << 219 return out; 196 return out; >> 197 220 } 198 } 221 199 >> 200 222 //-------------------------------------------- 201 //------------------------------------------------------------------------ 223 G4int G4ErrorFreeTrajState::PropagateError(con << 202 G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack ) 224 { 203 { 225 G4double stepLengthCm = aTrack->GetStep()->G << 204 G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; 226 if(G4ErrorPropagatorData::GetErrorPropagator << 205 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() == G4ErrorStage_Deflation ) stepLengthCm *= -1.; 227 G4ErrorStage_Deflation) << 228 stepLengthCm *= -1.; << 229 << 230 G4double kCarTolerance = << 231 G4GeometryTolerance::GetInstance()->GetSur << 232 206 233 if(std::fabs(stepLengthCm) <= kCarTolerance << 207 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 234 return 0; << 235 208 >> 209 if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0; >> 210 236 #ifdef G4EVERBOSE 211 #ifdef G4EVERBOSE 237 if(iverbose >= 2) << 212 if( iverbose >= 2 )G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl; 238 G4cout << " G4ErrorFreeTrajState::Propaga << 239 G4cout << "G4EP: iverbose=" << iverbose << G << 240 #endif 213 #endif 241 214 242 // * *** ERROR PROPAGATION ON A HELIX ASSUMI 215 // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES 243 G4Point3D vposPost = aTrack->GetPosition() / << 216 G4Point3D vposPost = aTrack->GetPosition()/cm; 244 G4Vector3D vpPost = aTrack->GetMomentum() / << 217 G4Vector3D vpPost = aTrack->GetMomentum()/GeV; 245 // G4Point3D vposPre = fPosition/cm; 218 // G4Point3D vposPre = fPosition/cm; 246 // G4Vector3D vpPre = fMomentum/GeV; 219 // G4Vector3D vpPre = fMomentum/GeV; 247 G4Point3D vposPre = aTrack->GetStep()->GetPr << 220 G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm; 248 G4Vector3D vpPre = aTrack->GetStep()->GetPr << 221 G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV; 249 // correct to avoid propagation along Z << 222 //correct to avoid propagation along Z 250 if(vpPre.mag() == vpPre.z()) << 223 if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV ); 251 vpPre.setX(1.E-6 * MeV); << 224 if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV ); 252 if(vpPost.mag() == vpPost.z()) << 253 vpPost.setX(1.E-6 * MeV); << 254 225 255 G4double pPre = vpPre.mag(); << 226 G4double pPre = vpPre.mag(); 256 G4double pPost = vpPost.mag(); 227 G4double pPost = vpPost.mag(); 257 #ifdef G4EVERBOSE 228 #ifdef G4EVERBOSE 258 if(iverbose >= 2) << 229 if( iverbose >= 2 ) { 259 { << 230 G4cout << "G4EP: vposPre " << vposPre << G4endl 260 G4cout << "G4EP: vposPre " << vposPre << G << 231 << "G4EP: vposPost " << vposPost << G4endl; 261 << vposPost << G4endl; << 232 G4cout << "G4EP: vpPre " << vpPre << G4endl 262 G4cout << "G4EP: vpPre " << vpPre << G4end << 233 << "G4EP: vpPost " << vpPost << G4endl; 263 << G4endl; << 264 G4cout << " err start step " << fError << 234 G4cout << " err start step " << fError << G4endl; 265 G4cout << "G4EP: stepLengthCm " << stepLen 235 G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl; 266 } 236 } 267 #endif 237 #endif 268 238 269 if(pPre == 0. || pPost == 0) << 239 if( pPre == 0. || pPost == 0 ) return 2; 270 return 2; << 240 G4double pInvPre = 1./pPre; 271 G4double pInvPre = 1. / pPre; << 241 G4double pInvPost = 1./pPost; 272 G4double pInvPost = 1. / pPost; << 273 G4double deltaPInv = pInvPost - pInvPre; 242 G4double deltaPInv = pInvPost - pInvPre; 274 if(iverbose >= 2) << 275 G4cout << "G4EP: pInvPre" << pInvPre << " << 276 << " deltaPInv:" << deltaPInv << G << 277 243 278 G4Vector3D vpPreNorm = vpPre * pInvPre; << 244 G4Vector3D vpPreNorm = vpPre * pInvPre; 279 G4Vector3D vpPostNorm = vpPost * pInvPost; 245 G4Vector3D vpPostNorm = vpPost * pInvPost; 280 if(iverbose >= 2) << 246 // if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl; 281 G4cout << "G4EP: vpPreNorm " << vpPreNorm << 247 //return if propagation along Z?? 282 << G4endl; << 248 if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4; 283 // return if propagation along Z?? << 249 G4double sinpPre = std::sin( vpPreNorm.theta() ); //cosine perpendicular to pPre = sine pPre 284 if(1. - std::fabs(vpPreNorm.z()) < kCarToler << 250 G4double sinpPost = std::sin( vpPostNorm.theta() ); //cosine perpendicular to pPost = sine pPost 285 return 4; << 251 G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() ); 286 if(1. - std::fabs(vpPostNorm.z()) < kCarTole << 287 return 4; << 288 G4double sinpPre = << 289 std::sin(vpPreNorm.theta()); // cosine pe << 290 G4double sinpPost = << 291 std::sin(vpPostNorm.theta()); // cosine p << 292 G4double sinpPostInv = 1. / std::sin(vpPostN << 293 252 294 #ifdef G4EVERBOSE 253 #ifdef G4EVERBOSE 295 if(iverbose >= 2) << 254 if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl; 296 G4cout << "G4EP: cosl " << sinpPre << " co << 297 #endif 255 #endif 298 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN 256 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR 299 //* *** NEUTRAL PARTICLE OR FIELDFREE REGION 257 //* *** NEUTRAL PARTICLE OR FIELDFREE REGION 300 G4ErrorMatrix transf(5, 5, 0); << 258 G4ErrorMatrix transf(5, 5, 0 ); 301 259 302 transf[3][2] = stepLengthCm * sinpPost; 260 transf[3][2] = stepLengthCm * sinpPost; 303 transf[4][1] = stepLengthCm; 261 transf[4][1] = stepLengthCm; 304 for(auto ii = 0; ii < 5; ++ii) << 262 for( size_t ii=0;ii < 5; ii++ ){ 305 { << 306 transf[ii][ii] = 1.; 263 transf[ii][ii] = 1.; 307 } 264 } 308 #ifdef G4EVERBOSE 265 #ifdef G4EVERBOSE 309 if(iverbose >= 2) << 266 if( iverbose >= 2 ) { 310 { << 311 G4cout << "G4EP: transf matrix neutral " < 267 G4cout << "G4EP: transf matrix neutral " << transf; 312 } 268 } 313 #endif 269 #endif 314 270 315 // charge X propagation direction 271 // charge X propagation direction 316 G4double charge = aTrack->GetDynamicParticle 272 G4double charge = aTrack->GetDynamicParticle()->GetCharge(); 317 if(G4ErrorPropagatorData::GetErrorPropagator << 273 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ) { 318 G4ErrorMode_PropBackwards) << 274 charge *= -1.; 319 { << 320 charge *= -1.; << 321 } 275 } 322 // G4cout << " charge " << charge << G4endl 276 // G4cout << " charge " << charge << G4endl; 323 // t check if particle has charge << 277 //t check if particle has charge 324 // t if( charge == 0 ) goto 45; << 278 //t if( charge == 0 ) goto 45; 325 // check if the magnetic field is = 0. 279 // check if the magnetic field is = 0. 326 280 327 // position is from geant4, it is assumed to << 281 //position is from geant4, it is assumed to be in mm (for debugging, eventually it will not be transformed) 328 // eventually it will not be transformed) it << 282 G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm; 329 // pos1[] is in mm. << 283 G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm; 330 G4double pos1[3]; << 331 pos1[0] = vposPre.x() * cm; << 332 pos1[1] = vposPre.y() * cm; << 333 pos1[2] = vposPre.z() * cm; << 334 G4double pos2[3]; << 335 pos2[0] = vposPost.x() * cm; << 336 pos2[1] = vposPost.y() * cm; << 337 pos2[2] = vposPost.z() * cm; << 338 G4double h1[3], h2[3]; 284 G4double h1[3], h2[3]; 339 285 340 const G4Field* field = G4TransportationManag << 286 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); 341 ->GetFieldManager() << 287 if( !field ) return 0; //goto 45 342 ->GetDetectorField( << 288 343 if(!field) << 344 return 0; // goto 45 << 345 << 346 // calculate transformation except it NEUTRA 289 // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION 347 if(charge != 0. && field) << 290 if( charge != 0. && field ) { 348 { << 349 field->GetFieldValue(pos1, << 350 h1); // here pos1[], << 351 field->GetFieldValue(pos2, h2); << 352 G4ThreeVector HPre = << 353 G4ThreeVector(h1[0], h1[1], h1[2]) / tes << 354 10.; // 10. is to get same dimensions a << 355 G4ThreeVector HPost = G4ThreeVector(h2[0], << 356 G4double magHPre = HPre.mag(); << 357 G4double magHPost = HPost.mag(); << 358 #ifdef G4EVERBOSE << 359 if(iverbose >= 2) << 360 { << 361 G4cout << "G4EP: h1 = " << h1[0] << ", " << 362 << G4endl; << 363 G4cout << "G4EP: pos1/mm = " << pos1[0] << 364 << pos1[2] << G4endl; << 365 G4cout << "G4EP: pos2/mm = " << pos2[0] << 366 << pos2[2] << G4endl; << 367 G4cout << "G4EP: B-filed in KGauss HPre << 368 << "G4EP: in KGauss HPost " << HP << 369 } << 370 #endif << 371 291 372 if(magHPre + magHPost != 0.) << 292 field->GetFieldValue( pos1, h1 ); 373 { << 293 field->GetFieldValue( pos2, h2 ); 374 //* *** CHECK WHETHER H*ALFA/P IS TOO DI << 294 G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; //10. is to get same dimensions as GEANT3 (kilogauss) 375 G4double gam; << 295 G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.; 376 if(magHPost != 0.) << 296 G4double magHPre = HPre.mag(); 377 { << 297 G4double magHPost = HPost.mag(); 378 gam = HPost * vpPostNorm / magHPost; << 298 #ifdef G4EVERBOSE 379 } << 299 if( iverbose >= 2 ) G4cout << "G4EP: HPre " << HPre << G4endl 380 else << 300 << "G4EP: HPost " << HPost << G4endl; 381 { << 301 #endif 382 gam = HPre * vpPreNorm / magHPre; << 302 383 } << 303 if( magHPre + magHPost != 0. ) { 384 << 304 385 // G4eMagneticLimitsProcess will limit t << 305 //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2 386 // line trajectory << 306 G4double gam; 387 G4double alphaSqr = 1. - gam * gam; << 307 if( magHPost != 0. ){ 388 G4double diffHSqr = (HPre * pInvPre - H << 308 gam = HPost * vpPostNorm / magHPost; 389 G4double delhp6Sqr = 300. * 300.; << 309 }else { 390 #ifdef G4EVERBOSE << 310 gam = HPre * vpPreNorm / magHPre; 391 if(iverbose >= 2) << 392 { << 393 G4cout << " G4EP: gam " << gam << " al << 394 << " diffHSqr " << diffHSqr << << 395 G4cout << " alpha= " << std::sqrt(alph << 396 } << 397 #endif << 398 if(diffHSqr * alphaSqr > delhp6Sqr) << 399 return 3; << 400 << 401 //* *** DEFINE AVERAGE MAGNETIC FIELD AN << 402 G4double pInvAver = 1. / (pInvPre + pInv << 403 G4double CFACT8 = 2.997925E-4; << 404 // G4double HAver << 405 G4ThreeVector vHAverNorm((HPre * pInvPre << 406 charge * CFACT8 << 407 G4double HAver = vHAverNorm.mag(); << 408 G4double invHAver = 1. / HAver; << 409 vHAverNorm *= invHAver; << 410 #ifdef G4EVERBOSE << 411 if(iverbose >= 2) << 412 G4cout << " G4EP: HaverNorm " << vHAve << 413 << " charge " << charge << G4en << 414 #endif << 415 << 416 G4double pAver = (pPre + pPost) * << 417 G4double QAver = -HAver / pAver; << 418 G4double thetaAver = QAver * stepLeng << 419 G4double sinThetaAver = std::sin(thetaAv << 420 G4double cosThetaAver = std::cos(thetaAv << 421 G4double gamma = vHAverNorm * vpP << 422 G4ThreeVector AN2 = vHAverNorm.cross << 423 << 424 #ifdef G4EVERBOSE << 425 if(iverbose >= 2) << 426 G4cout << " G4EP: AN2 " << AN2 << " g << 427 << " theta=" << thetaAver << G << 428 #endif << 429 G4double AU = 1. / vpPreNorm.perp(); << 430 // t G4ThreeVector vU( vpPreNorm.cross( << 431 G4ThreeVector vUPre(-AU * vpPreNorm.y(), << 432 G4ThreeVector vVPre(-vpPreNorm.z() * vUP << 433 vpPreNorm.x() * vUPr << 434 vpPreNorm.y() * vU << 435 << 436 // << 437 AU = 1. / vpPostNorm.perp(); << 438 // t G4ThreeVector vU( vpPostNorm.cross << 439 // ); << 440 G4ThreeVector vUPost(-AU * vpPostNorm.y( << 441 G4ThreeVector vVPost( << 442 -vpPostNorm.z() * vUPost.y(), vpPostNo << 443 vpPostNorm.x() * vUPost.y() - vpPostNo << 444 #ifdef G4EVERBOSE << 445 G4cout << " vpPostNorm " << vpPostNorm < << 446 if(iverbose >= 2) << 447 G4cout << " G4EP: AU " << AU << " vUPr << 448 << " vUPost " << vUPost << " vV << 449 #endif << 450 G4Point3D deltaPos(vposPre - vposPost); << 451 << 452 // * *** COMPLETE TRANSFORMATION MATRIX << 453 // * *** FIELD GRADIENT PERPENDICULAR TO << 454 // * *** TAKEN INTO ACCOUNT << 455 << 456 G4double QP = QAver * pAver; // = -HAve << 457 #ifdef G4EVERBOSE << 458 if(iverbose >= 2) << 459 G4cout << " G4EP: QP " << QP << " QAve << 460 << G4endl; << 461 #endif << 462 G4double ANV = << 463 -(vHAverNorm.x() * vUPost.x() + vHAver << 464 G4double ANU = << 465 (vHAverNorm.x() * vVPost.x() + vHAverN << 466 vHAverNorm.z() * vVPost.z()); << 467 G4double OMcosThetaAver = 1. - cosThetaA << 468 #ifdef G4EVERBOSE << 469 if(iverbose >= 2) << 470 G4cout << "G4EP: OMcosThetaAver " << O << 471 << cosThetaAver << " thetaAver << 472 << QAver << " stepLengthCm " << << 473 #endif << 474 G4double TMSINT = thetaAver - sinThetaAv << 475 #ifdef G4EVERBOSE << 476 if(iverbose >= 2) << 477 G4cout << " G4EP: ANV " << ANV << " AN << 478 #endif << 479 << 480 G4ThreeVector vHUPre( << 481 -vHAverNorm.z() * vUPre.y(), vHAverNor << 482 vHAverNorm.x() * vUPre.y() - vHAverNor << 483 #ifdef G4EVERBOSE << 484 // if( iverbose >= 2 ) G4cout << "G4E << 485 // << vHAverNorm.z() << " " << vUPre. << 486 #endif << 487 G4ThreeVector vHVPre( << 488 vHAverNorm.y() * vVPre.z() - vHAverNor << 489 vHAverNorm.z() * vVPre.x() - vHAverNor << 490 vHAverNorm.x() * vVPre.y() - vHAverNor << 491 #ifdef G4EVERBOSE << 492 if(iverbose >= 2) << 493 G4cout << " G4EP: HUPre " << vHUPre << << 494 #endif << 495 << 496 //------------------- COMPUTE MATRIX << 497 //---------- 1/P << 498 << 499 transf[0][0] = << 500 1. - << 501 deltaPInv * pAver * << 502 (1. + (vpPostNorm.x() * deltaPos.x() << 503 vpPostNorm.z() * deltaPos.z() << 504 stepLengthCm) + << 505 2. * deltaPInv * pAver; << 506 << 507 transf[0][1] = << 508 -deltaPInv / thetaAver * << 509 (TMSINT * gamma * << 510 (vHAverNorm.x() * vVPre.x() + vHAve << 511 vHAverNorm.z() * vVPre.z()) + << 512 sinThetaAver * << 513 (vVPre.x() * vpPostNorm.x() + vVPre << 514 vVPre.z() * vpPostNorm.z()) + << 515 OMcosThetaAver * << 516 (vHVPre.x() * vpPostNorm.x() + vHVP << 517 vHVPre.z() * vpPostNorm.z())); << 518 << 519 transf[0][2] = << 520 -sinpPre * deltaPInv / thetaAver * << 521 (TMSINT * gamma * << 522 (vHAverNorm.x() * vUPre.x() + vHAve << 523 sinThetaAver * << 524 (vUPre.x() * vpPostNorm.x() + vUPre << 525 OMcosThetaAver * << 526 (vHUPre.x() * vpPostNorm.x() + vHUP << 527 vHUPre.z() * vpPostNorm.z())); << 528 << 529 transf[0][3] = -deltaPInv / stepLengthCm << 530 (vUPre.x() * vpPostNorm.x << 531 << 532 transf[0][4] = -deltaPInv / stepLengthCm << 533 (vVPre.x() * vpPostNorm.x << 534 vVPre.z() * vpPostNorm.z << 535 << 536 // *** Lambda << 537 transf[1][0] = << 538 -QP * ANV * << 539 (vpPostNorm.x() * deltaPos.x() + vpPos << 540 vpPostNorm.z() * deltaPos.z()) * << 541 (1. + deltaPInv * pAver); << 542 #ifdef G4EVERBOSE << 543 if(iverbose >= 3) << 544 G4cout << "ctransf10= " << transf[1][0 << 545 << " " << vpPostNorm.x() << " " << 546 << vpPostNorm.y() << " " << del << 547 << " " << deltaPos.z() << " " < << 548 << G4endl; << 549 #endif << 550 << 551 transf[1][1] = << 552 cosThetaAver * (vVPre.x() * vVPost.x() << 553 vVPre.z() * vVPost.z() << 554 sinThetaAver * (vHVPre.x() * vVPost.x( << 555 vHVPre.z() * vVPost.z( << 556 OMcosThetaAver * << 557 (vHAverNorm.x() * vVPre.x() + vHAver << 558 vHAverNorm.z() * vVPre.z()) * << 559 (vHAverNorm.x() * vVPost.x() + vHAve << 560 vHAverNorm.z() * vVPost.z()) + << 561 ANV * (-sinThetaAver * << 562 (vVPre.x() * vpPostNorm.x() + << 563 vVPre.z() * vpPostNorm.z()) << 564 OMcosThetaAver * (vVPre.x() * A << 565 vVPre.z() * A << 566 TMSINT * gamma * << 567 (vHAverNorm.x() * vVPre.x() + << 568 vHAverNorm.z() * vVPre.z())) << 569 << 570 transf[1][2] = << 571 cosThetaAver * (vUPre.x() * vVPost.x() << 572 sinThetaAver * (vHUPre.x() * vVPost.x( << 573 vHUPre.z() * vVPost.z( << 574 OMcosThetaAver * << 575 (vHAverNorm.x() * vUPre.x() + vHAver << 576 (vHAverNorm.x() * vVPost.x() + vHAve << 577 vHAverNorm.z() * vVPost.z()) + << 578 ANV * (-sinThetaAver * << 579 (vUPre.x() * vpPostNorm.x() + << 580 OMcosThetaAver * (vUPre.x() * A << 581 TMSINT * gamma * << 582 (vHAverNorm.x() * vUPre.x() + << 583 transf[1][2] = sinpPre * transf[1][2]; << 584 << 585 transf[1][3] = -QAver * ANV * << 586 (vUPre.x() * vpPostNorm.x << 587 << 588 transf[1][4] = -QAver * ANV * << 589 (vVPre.x() * vpPostNorm.x << 590 vVPre.z() * vpPostNorm.z << 591 << 592 // *** Phi << 593 << 594 transf[2][0] = << 595 -QP * ANU * << 596 (vpPostNorm.x() * deltaPos.x() + vpPos << 597 vpPostNorm.z() * deltaPos.z()) * << 598 sinpPostInv * (1. + deltaPInv * pAver) << 599 #ifdef G4EVERBOSE << 600 if(iverbose >= 3) << 601 G4cout << "ctransf20= " << transf[2][0 << 602 << " " << vpPostNorm.x() << " " << 603 << vpPostNorm.y() << " " << del << 604 << " " << deltaPos.z() << " " < << 605 << " " << pAver << G4endl; << 606 #endif << 607 transf[2][1] = << 608 cosThetaAver * (vVPre.x() * vUPost.x() << 609 sinThetaAver * (vHVPre.x() * vUPost.x( << 610 OMcosThetaAver * << 611 (vHAverNorm.x() * vVPre.x() + vHAver << 612 vHAverNorm.z() * vVPre.z()) * << 613 (vHAverNorm.x() * vUPost.x() + vHAve << 614 ANU * (-sinThetaAver * << 615 (vVPre.x() * vpPostNorm.x() + << 616 vVPre.z() * vpPostNorm.z()) << 617 OMcosThetaAver * (vVPre.x() * A << 618 vVPre.z() * A << 619 TMSINT * gamma * << 620 (vHAverNorm.x() * vVPre.x() + << 621 vHAverNorm.z() * vVPre.z())) << 622 transf[2][1] = sinpPostInv * transf[2][1 << 623 << 624 transf[2][2] = << 625 cosThetaAver * (vUPre.x() * vUPost.x() << 626 sinThetaAver * (vHUPre.x() * vUPost.x( << 627 OMcosThetaAver * << 628 (vHAverNorm.x() * vUPre.x() + vHAver << 629 (vHAverNorm.x() * vUPost.x() + vHAve << 630 ANU * (-sinThetaAver * << 631 (vUPre.x() * vpPostNorm.x() + << 632 OMcosThetaAver * (vUPre.x() * A << 633 TMSINT * gamma * << 634 (vHAverNorm.x() * vUPre.x() + << 635 transf[2][2] = sinpPostInv * sinpPre * t << 636 << 637 transf[2][3] = -QAver * ANU * << 638 (vUPre.x() * vpPostNorm.x << 639 sinpPostInv; << 640 #ifdef G4EVERBOSE << 641 if(iverbose >= 3) << 642 G4cout << "ctransf23= " << transf[2][3 << 643 << " " << vUPre.x() << " " << v << 644 << " " << vpPostNorm.y() << " " << 645 #endif << 646 << 647 transf[2][4] = -QAver * ANU * << 648 (vVPre.x() * vpPostNorm.x << 649 vVPre.z() * vpPostNorm.z << 650 sinpPostInv; << 651 << 652 // *** Yt << 653 << 654 transf[3][0] = pAver * << 655 (vUPost.x() * deltaPos.x( << 656 (1. + deltaPInv * pAver); << 657 #ifdef G4EVERBOSE << 658 if(iverbose >= 3) << 659 G4cout << "ctransf30= " << transf[3][0 << 660 << vUPost.x() << " " << deltaPo << 661 << deltaPos.y() << " " << delta << 662 #endif << 663 << 664 transf[3][1] = << 665 (sinThetaAver * (vVPre.x() * vUPost.x( << 666 OMcosThetaAver * (vHVPre.x() * vUPost << 667 TMSINT * (vHAverNorm.x() * vUPost.x() << 668 (vHAverNorm.x() * vVPre.x() + vHAve << 669 vHAverNorm.z() * vVPre.z())) / << 670 QAver; << 671 << 672 transf[3][2] = << 673 (sinThetaAver * (vUPre.x() * vUPost.x( << 674 OMcosThetaAver * (vHUPre.x() * vUPost << 675 TMSINT * (vHAverNorm.x() * vUPost.x() << 676 (vHAverNorm.x() * vUPre.x() + vHAve << 677 sinpPre / QAver; << 678 #ifdef G4EVERBOSE << 679 if(iverbose >= 3) << 680 G4cout << "ctransf32= " << transf[3][2 << 681 << vUPre.x() << " " << vUPost.x << 682 << vUPost.y() << " " << OMcosTh << 683 << " " << vUPost.x() << " " << << 684 << " " << TMSINT << " " << vHAv << 685 << " " << vHAverNorm.y() << " " << 686 << vHAverNorm.x() << " " << vUP << 687 << " " << vUPre.y() << " " << s << 688 #endif << 689 << 690 transf[3][3] = (vUPre.x() * vUPost.x() + << 691 << 692 transf[3][4] = (vVPre.x() * vUPost.x() + << 693 << 694 // *** Zt << 695 transf[4][0] = pAver * << 696 (vVPost.x() * deltaPos.x( << 697 vVPost.z() * deltaPos.z( << 698 (1. + deltaPInv * pAver); << 699 << 700 transf[4][1] = << 701 (sinThetaAver * (vVPre.x() * vVPost.x( << 702 vVPre.z() * vVPost.z( << 703 OMcosThetaAver * (vHVPre.x() * vVPost << 704 vHVPre.z() * vVPost << 705 TMSINT * << 706 (vHAverNorm.x() * vVPost.x() + vHAv << 707 vHAverNorm.z() * vVPost.z()) * << 708 (vHAverNorm.x() * vVPre.x() + vHAve << 709 vHAverNorm.z() * vVPre.z())) / << 710 QAver; << 711 #ifdef G4EVERBOSE << 712 if(iverbose >= 3) << 713 G4cout << "ctransf41= " << transf[4][1 << 714 << OMcosThetaAver << " " << TMS << 715 << vVPost << " " << vHVPre << " << 716 << G4endl; << 717 #endif << 718 << 719 transf[4][2] = << 720 (sinThetaAver * (vUPre.x() * vVPost.x( << 721 OMcosThetaAver * (vHUPre.x() * vVPost << 722 vHUPre.z() * vVPost << 723 TMSINT * << 724 (vHAverNorm.x() * vVPost.x() + vHAv << 725 vHAverNorm.z() * vVPost.z()) * << 726 (vHAverNorm.x() * vUPre.x() + vHAve << 727 sinpPre / QAver; << 728 << 729 transf[4][3] = (vUPre.x() * vVPost.x() + << 730 << 731 transf[4][4] = (vVPre.x() * vVPost.x() + << 732 vVPre.z() * vVPost.z()); << 733 // if(iverbose >= 3) G4cout <<"ctransf << 734 // vVPre.x() <<" "<<vVPost.x() <<" "< << 735 // "<< vVPre.z() <<" "<< vVPost.z() << << 736 << 737 #ifdef G4EVERBOSE << 738 if(iverbose >= 1) << 739 G4cout << "G4EP: transf matrix compute << 740 #endif << 741 /* for( G4int ii=0;ii<5;ii++){ << 742 for( G4int jj=0;jj<5;jj++){ << 743 G4cout << transf[ii][jj] << " "; << 744 } << 745 G4cout << G4endl; << 746 } */ << 747 } 311 } 748 } << 312 749 // end of calculate transformation except it << 313 // G4eMagneticLimitsProcess will limit the step, but based on an straight line trajectory 750 // REGION << 314 G4double alphaSqr = 1. - gam * gam; 751 /* if( iverbose >= 1 ) G4cout << "G4EP: tra << 315 G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2(); 752 << theFirstStep << G4endl; if( theFirstStep << 316 G4double delhp6Sqr = 300.*300.; >> 317 #ifdef G4EVERBOSE >> 318 if( iverbose >= 2 ) G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr << " diffHSqr " << diffHSqr << G4endl; >> 319 #endif >> 320 if( diffHSqr * alphaSqr > delhp6Sqr ) return 3; >> 321 >> 322 >> 323 //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT >> 324 G4double pInvAver = 1./(pInvPre + pInvPost ); >> 325 G4double CFACT8 = 2.997925E-4; >> 326 //G4double HAver >> 327 G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 ); >> 328 G4double HAver = vHAverNorm.mag(); >> 329 G4double invHAver = 1./HAver; >> 330 vHAverNorm *= invHAver; >> 331 #ifdef G4EVERBOSE >> 332 if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl; >> 333 #endif >> 334 >> 335 G4double pAver = (pPre+pPost)*0.5; >> 336 G4double QAver = -HAver/pAver; >> 337 G4double thetaAver = QAver * stepLengthCm; >> 338 G4double sinThetaAver = std::sin(thetaAver); >> 339 G4double cosThetaAver = std::cos(thetaAver); >> 340 G4double gamma = vHAverNorm * vpPostNorm; >> 341 G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm ); >> 342 >> 343 #ifdef G4EVERBOSE >> 344 if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << G4endl; >> 345 #endif >> 346 G4double AU = 1./vpPreNorm.perp(); >> 347 //t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU ); >> 348 G4ThreeVector vUPre( -AU*vpPreNorm.y(), >> 349 AU*vpPreNorm.x(), >> 350 0. ); >> 351 G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(), >> 352 vpPreNorm.z()*vUPre.x(), >> 353 vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() ); >> 354 >> 355 // >> 356 AU = 1./vpPostNorm.perp(); >> 357 //t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU ); >> 358 G4ThreeVector vUPost( -AU*vpPostNorm.y(), >> 359 AU*vpPostNorm.x(), >> 360 0. ); >> 361 G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(), >> 362 vpPostNorm.z()*vUPost.x(), >> 363 vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() ); >> 364 #ifdef G4EVERBOSE >> 365 //- G4cout << " vpPostNorm " << vpPostNorm << G4endl; >> 366 if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl; >> 367 #endif >> 368 G4Point3D deltaPos( vposPre - vposPost ); >> 369 >> 370 // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2 >> 371 // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT >> 372 // * *** TAKEN INTO ACCOUNT >> 373 >> 374 G4double QP = QAver * pAver; // = -HAver >> 375 #ifdef G4EVERBOSE >> 376 if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl; >> 377 #endif >> 378 G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() ); >> 379 G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() ); >> 380 G4double OMcosThetaAver = 1. - cosThetaAver; >> 381 #ifdef G4EVERBOSE >> 382 if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl; >> 383 #endif >> 384 G4double TMSINT = thetaAver - sinThetaAver; >> 385 #ifdef G4EVERBOSE >> 386 if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl; >> 387 #endif >> 388 >> 389 G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(), >> 390 vHAverNorm.z() * vUPre.x(), >> 391 vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() ); >> 392 #ifdef G4EVERBOSE >> 393 // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " " << vHAverNorm.z() << " " << vUPre.y() << G4endl; >> 394 #endif >> 395 G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(), >> 396 vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(), >> 397 vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() ); >> 398 #ifdef G4EVERBOSE >> 399 if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl; >> 400 #endif >> 401 >> 402 //------------------- COMPUTE MATRIX >> 403 //---------- 1/P >> 404 >> 405 transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm) >> 406 +2.*deltaPInv*pAver; >> 407 >> 408 transf[0][1] = -deltaPInv/thetaAver* >> 409 ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) + >> 410 sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + >> 411 OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) ); >> 412 >> 413 transf[0][2] = -sinpPre*deltaPInv/thetaAver* >> 414 ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) + >> 415 sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + >> 416 OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) ); >> 417 >> 418 transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ); >> 419 >> 420 transf[0][4] = -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()); >> 421 >> 422 // *** Lambda >> 423 transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z()) >> 424 *(1.+deltaPInv*pAver); >> 425 #ifdef G4EVERBOSE >> 426 if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0] << " " << -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z() >> 427 << " " << deltaPInv<< " " << pAver << G4endl; >> 428 #endif >> 429 >> 430 transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + >> 431 sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + >> 432 OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* >> 433 (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + >> 434 ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + >> 435 OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - >> 436 TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) ); >> 437 >> 438 transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + >> 439 sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + >> 440 OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* >> 441 (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + >> 442 ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + >> 443 OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - >> 444 TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) ); >> 445 transf[1][2] = sinpPre*transf[1][2]; >> 446 >> 447 transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ); >> 448 >> 449 transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()); >> 450 >> 451 // *** Phi >> 452 >> 453 transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv >> 454 *(1.+deltaPInv*pAver); >> 455 #ifdef G4EVERBOSE >> 456 if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv >> 457 <<" "<<deltaPInv<<" "<<pAver<< G4endl; >> 458 #endif >> 459 transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) + >> 460 sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) + >> 461 OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* >> 462 (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) + >> 463 ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + >> 464 OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - >> 465 TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) ); >> 466 transf[2][1] = sinpPostInv*transf[2][1]; >> 467 >> 468 transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) + >> 469 sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) + >> 470 OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* >> 471 (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) + >> 472 ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + >> 473 OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - >> 474 TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) ); >> 475 transf[2][2] = sinpPostInv*sinpPre*transf[2][2]; >> 476 >> 477 transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv; >> 478 #ifdef G4EVERBOSE >> 479 if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl; >> 480 #endif >> 481 >> 482 transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv; >> 483 >> 484 // *** Yt >> 485 >> 486 transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() ) >> 487 *(1.+deltaPInv*pAver); >> 488 #ifdef G4EVERBOSE >> 489 if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y() >> 490 <<" "<<deltaPInv<<" "<<pAver<<G4endl; >> 491 #endif >> 492 >> 493 transf[3][1] = ( sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) + >> 494 OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) + >> 495 TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )* >> 496 (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver; >> 497 >> 498 transf[3][2] = ( sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) + >> 499 OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) + >> 500 TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )* >> 501 (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver; >> 502 #ifdef G4EVERBOSE >> 503 if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<< >> 504 OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<< >> 505 TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<< >> 506 vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl; >> 507 #endif >> 508 >> 509 transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ); >> 510 >> 511 transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ); >> 512 >> 513 // *** Zt >> 514 transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z()) >> 515 *(1.+deltaPInv*pAver); >> 516 >> 517 transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + >> 518 OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + >> 519 TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* >> 520 (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver; >> 521 #ifdef G4EVERBOSE >> 522 if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl; >> 523 #endif >> 524 >> 525 transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + >> 526 OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + >> 527 TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* >> 528 (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver; >> 529 >> 530 transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ); >> 531 >> 532 transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()); >> 533 // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<< vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<" "<< vVPre.z() <<" "<< vVPost.z() << G4endl; >> 534 >> 535 >> 536 #ifdef G4EVERBOSE >> 537 if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl; >> 538 #endif >> 539 /* for( G4int ii=0;ii<5;ii++){ >> 540 for( G4int jj=0;jj<5;jj++){ >> 541 G4cout << transf[ii][jj] << " "; >> 542 } >> 543 G4cout << G4endl; >> 544 } */ >> 545 } >> 546 } >> 547 // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION >> 548 /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized " << theFirstStep << G4endl; >> 549 if( theFirstStep ) { >> 550 theTransfMat = transf; 753 theFirstStep = false; 551 theFirstStep = false; 754 }else{ 552 }else{ 755 theTransfMat = theTransfMat * transf; 553 theTransfMat = theTransfMat * transf; 756 if( iverbose >= 1 ) G4cout << "G4EP: trans << 554 if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" << theTransfMat << G4endl; 757 theTransfMat << G4endl; << 555 } 758 } << 759 */ 556 */ 760 theTransfMat = std::move(transf); << 557 theTransfMat = transf; 761 #ifdef G4EVERBOSE << 762 if(iverbose >= 1) << 763 G4cout << "G4EP: error matrix before trans << 764 if(iverbose >= 2) << 765 G4cout << " tf * err " << theTransfMat * f << 766 << " transf matrix " << theTransfMa << 767 #endif << 768 << 769 fError = fError.similarity(theTransfMat).T() << 770 //- fError = transf * fError * transf.T() << 771 #ifdef G4EVERBOSE 558 #ifdef G4EVERBOSE 772 if(iverbose >= 1) << 559 if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl; 773 G4cout << "G4EP: error matrix propagated " << 560 if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl 774 #endif << 561 << " transf matrix " << theTransfMat.T() << G4endl; 775 << 562 #endif 776 //? S = B*S*BT S.similarity(B) << 563 777 //? R = S << 564 fError = fError.similarity(theTransfMat).T(); 778 // not needed * *** TRANSFORM ERROR MATRIX F << 565 //- fError = transf * fError * transf.T(); 779 // VARIABLES; << 566 #ifdef G4EVERBOSE 780 << 567 if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl; 781 PropagateErrorMSC(aTrack); << 568 #endif 782 << 569 783 PropagateErrorIoni(aTrack); << 570 //? S = B*S*BT S.similarity(B) 784 << 571 //? R = S 785 return 0; << 572 // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES; >> 573 >> 574 PropagateErrorMSC( aTrack ); >> 575 >> 576 PropagateErrorIoni( aTrack ); >> 577 >> 578 return 0; 786 } 579 } 787 580 >> 581 788 //-------------------------------------------- 582 //------------------------------------------------------------------------ 789 G4int G4ErrorFreeTrajState::PropagateErrorMSC( << 583 G4int G4ErrorFreeTrajState::PropagateErrorMSC( const G4Track* aTrack ) 790 { << 584 { 791 G4ThreeVector vpPre = aTrack->GetMomentum( << 585 G4ThreeVector vpPre = aTrack->GetMomentum()/GeV; 792 G4double pPre = vpPre.mag(); << 586 G4double pPre = vpPre.mag(); 793 G4double pBeta = pPre * pPre / (aTrac << 587 G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV); 794 G4double stepLengthCm = aTrack->GetStep()->G << 588 G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; 795 589 796 G4Material* mate = aTrack->GetVolume()->GetL 590 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial(); 797 G4double effZ, effA; 591 G4double effZ, effA; 798 CalculateEffectiveZandA(mate, effZ, effA); << 592 CalculateEffectiveZandA( mate, effZ, effA ); 799 593 800 #ifdef G4EVERBOSE 594 #ifdef G4EVERBOSE 801 if(iverbose >= 4) << 595 if( iverbose >= 4 ) G4cout << "material " << mate->GetName() 802 G4cout << "material " << 596 //<< " " << mate->GetZ() << " " << mate->GetA() 803 << mate->GetName() << 597 << " " << effZ << " " << effA 804 //<< " " << mate->GetZ() << " " << << 598 << " " << mate->GetDensity()/g*mole << " " << mate->GetRadlen()/cm << " " << mate->GetNuclearInterLength()/cm << G4endl; 805 << " effZ:" << effZ << " effA:" << << 806 << " dens(g/mole):" << mate->GetDen << 807 << " Radlen/cm:" << mate->GetRadlen << 808 << mate->GetNuclearInterLength() / << 809 #endif 599 #endif 810 600 811 G4double RI = stepLengthCm / (mate->GetRadle << 601 G4double RI = stepLengthCm / (mate->GetRadlen()/cm); 812 #ifdef G4EVERBOSE 602 #ifdef G4EVERBOSE 813 if(iverbose >= 4) << 603 if( iverbose >= 4 ) G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI " << RI << " stepLengthCm " << stepLengthCm << " radlen " << (mate->GetRadlen()/cm) << " " << RI*1.e10 << G4endl; 814 G4cout << std::setprecision(6) << std::set << 815 << " stepLengthCm " << stepLengthCm << 816 << (mate->GetRadlen() / cm) << " RI << 817 #endif 604 #endif 818 G4double charge = aTrack->GetDynamicParticle 605 G4double charge = aTrack->GetDynamicParticle()->GetCharge(); 819 G4double DD = 1.8496E-4 * RI * (charge / << 606 G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta ); 820 #ifdef G4EVERBOSE 607 #ifdef G4EVERBOSE 821 if(iverbose >= 3) << 608 if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl; 822 G4cout << "G4EP:MSC: D*1E6= " << DD * 1.E6 << 823 #endif 609 #endif 824 G4double S1 = DD * stepLengthCm * stepLength << 610 G4double S1 = DD*stepLengthCm*stepLengthCm/3.; 825 G4double S2 = DD; 611 G4double S2 = DD; 826 G4double S3 = DD * stepLengthCm / 2.; << 612 G4double S3 = DD*stepLengthCm/2.; 827 613 828 G4double CLA = << 614 G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre; 829 std::sqrt(vpPre.x() * vpPre.x() + vpPre.y( << 830 #ifdef G4EVERBOSE 615 #ifdef G4EVERBOSE 831 if(iverbose >= 2) << 616 if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 " << S2 << " S3 " << S3 << " CLA " << CLA << G4endl; 832 G4cout << std::setw(6) << "G4EP:MSC: RI " << 833 << S2 << " S3 " << S3 << " CLA " << << 834 #endif 617 #endif 835 fError[1][1] += S2; 618 fError[1][1] += S2; 836 fError[1][4] -= S3; 619 fError[1][4] -= S3; 837 fError[2][2] += S2 / CLA / CLA; << 620 fError[2][2] += S2/CLA/CLA; 838 fError[2][3] += S3 / CLA; << 621 fError[2][3] += S3/CLA; 839 fError[3][3] += S1; 622 fError[3][3] += S1; 840 fError[4][4] += S1; 623 fError[4][4] += S1; 841 624 842 #ifdef G4EVERBOSE 625 #ifdef G4EVERBOSE 843 if(iverbose >= 2) << 626 if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl; 844 G4cout << "G4EP:MSC: error matrix propagat << 845 #endif 627 #endif 846 628 847 return 0; 629 return 0; 848 } 630 } 849 631 >> 632 850 //-------------------------------------------- 633 //------------------------------------------------------------------------ 851 void G4ErrorFreeTrajState::CalculateEffectiveZ << 634 void G4ErrorFreeTrajState::CalculateEffectiveZandA( const G4Material* mate, G4double& effZ, G4double& effA ) 852 << 635 { 853 << 854 { << 855 effZ = 0.; 636 effZ = 0.; 856 effA = 0.; 637 effA = 0.; 857 auto nelem = mate->GetNumberOfElements(); << 638 G4int ii, nelem = mate->GetNumberOfElements(); 858 const G4double* fracVec = mate->GetFractionV 639 const G4double* fracVec = mate->GetFractionVector(); 859 for(G4int ii = 0; ii < (G4int)nelem; ++ii) << 640 for(ii=0; ii < nelem; ii++ ) { 860 { << 641 effZ += mate->GetElement( ii )->GetZ() * fracVec[ii]; 861 effZ += mate->GetElement(ii)->GetZ() * fra << 642 effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole; 862 effA += mate->GetElement(ii)->GetA() * fra << 863 } 643 } >> 644 864 } 645 } 865 646 >> 647 866 //-------------------------------------------- 648 //------------------------------------------------------------------------ 867 G4int G4ErrorFreeTrajState::PropagateErrorIoni << 649 G4int G4ErrorFreeTrajState::PropagateErrorIoni( const G4Track* aTrack ) 868 { << 650 { 869 G4double stepLengthCm = aTrack->GetStep()->G << 651 G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; 870 #ifdef G4EVERBOSE << 871 G4double DEDX2; 652 G4double DEDX2; 872 if(stepLengthCm < 1.E-7) << 653 if( stepLengthCm < 1.E-7 ) { 873 { << 654 DEDX2=0.; 874 DEDX2 = 0.; << 875 } 655 } 876 #endif << 877 // * Calculate xi factor (KeV). 656 // * Calculate xi factor (KeV). 878 G4Material* mate = aTrack->GetVolume()->GetL 657 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial(); 879 G4double effZ, effA; 658 G4double effZ, effA; 880 CalculateEffectiveZandA(mate, effZ, effA); << 659 CalculateEffectiveZandA( mate, effZ, effA ); 881 660 882 G4double Etot = aTrack->GetTotalEnergy() / << 661 G4double Etot = aTrack->GetTotalEnergy()/GeV; 883 G4double beta = aTrack->GetMomentum().mag() << 662 G4double beta = aTrack->GetMomentum().mag()/GeV / Etot; 884 G4double mass = aTrack->GetDynamicParticle( << 663 G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV; 885 G4double gamma = Etot / mass; 664 G4double gamma = Etot / mass; 886 << 665 887 // * Calculate xi factor (keV). << 666 // * Calculate xi factor (KeV). 888 G4double XI = 153.5 * effZ * stepLengthCm * << 667 G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) / 889 (effA * beta * beta); << 668 (effA*beta*beta); 890 669 891 #ifdef G4EVERBOSE 670 #ifdef G4EVERBOSE 892 if(iverbose >= 2) << 671 if( iverbose >= 2 ){ 893 { << 672 G4cout << "G4EP:IONI: XI " << XI << " beta " << beta << " gamma " << gamma << G4endl; 894 G4cout << "G4EP:IONI: XI/keV " << XI << " << 673 G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl; 895 << gamma << G4endl; << 896 G4cout << " density " << (mate->GetDensity << 897 << effA << " step " << stepLengthCm << 898 } 674 } 899 #endif 675 #endif 900 // * Maximum energy transfer to atomic e 676 // * Maximum energy transfer to atomic electron (KeV). 901 G4double eta = beta * gamma; << 677 G4double eta = beta*gamma; 902 G4double etasq = eta * eta; << 678 G4double etasq = eta*eta; 903 G4double eMass = 0.51099906 / GeV; << 679 G4double eMass = 0.51099906/GeV; 904 G4double massRatio = eMass / mass; 680 G4double massRatio = eMass / mass; 905 G4double F1 = 2 * eMass * etasq; << 681 G4double F1 = 2*eMass*etasq; 906 G4double F2 = 1. + 2. * massRatio * g << 682 G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio; 907 G4double Emax = 1.E+6 * F1 / F2; // no << 683 G4double Emax = 1.E+6*F1/F2; 908 684 909 // * *** and now sigma**2 in GeV 685 // * *** and now sigma**2 in GeV 910 G4double dedxSq = << 686 G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12; 911 XI * Emax * (1. - (beta * beta / 2.)) * 1. << 912 /*The above formula for var(1/p) good for d << 913 passing through a gas it leads to overesti << 914 electrons the Emax is almost equal to inci << 915 k=Xi/Emax as small as e-6 and gradually << 916 << 917 http://www2.pv.infn.it/~rotondi/kalman_1.p << 918 << 919 Since I do not have enough info at the mom << 920 sub-Landau models for k=Xi/Emax <0.01 I'll << 921 */ << 922 << 923 if(XI / Emax < 0.01) << 924 dedxSq *= << 925 XI / Emax * 100; // Quench for low Elos << 926 << 927 #ifdef G4EVERBOSE 687 #ifdef G4EVERBOSE 928 if(iverbose >= 2) << 688 if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX2 " << dedxSq << " emass " << eMass << " Emax " << Emax << G4endl; 929 G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << d << 930 << " Emax/keV: " << Emax << " k=Xi << 931 << 932 #endif 689 #endif 933 690 934 G4double pPre6 = << 691 // if( iverbose >= 2 ) G4cout << "G4EP:IONI: Etot " << Etot << " DEDX2 " << dedxSq << " emass " << eMass << G4endl; 935 (aTrack->GetStep()->GetPreStepPoint()->Get << 692 936 pPre6 = std::pow(pPre6, 6); << 693 G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag(); 937 // Apply it to error << 694 pPre6 = std::pow(pPre6, 6 ); 938 fError[0][0] += Etot * Etot * dedxSq / pPre6 << 695 //Apply it to error >> 696 fError[0][0] += Etot*Etot*dedxSq / pPre6; 939 #ifdef G4EVERBOSE 697 #ifdef G4EVERBOSE 940 if(iverbose >= 2) << 698 if( iverbose >= 2 ) G4cout << "G4:IONI getot " << Etot << " dedx2 " << dedxSq << " p " << pPre6 << G4endl; 941 G4cout << "G4:IONI Etot/GeV: " << Etot << << 699 if( iverbose >= 2 ) G4cout << "G4EP:IONI: error_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl; 942 << " p^6: " << pPre6 << G4endl; << 943 if(iverbose >= 2) << 944 G4cout << "G4EP:IONI: error2_from_ionisati << 945 << (Etot * Etot * dedxSq) / pPre6 < << 946 #endif 700 #endif 947 701 948 return 0; 702 return 0; 949 } 703 } >> 704 950 705