[815] | 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 | // $Id: G4ErrorFreeTrajState.cc,v 1.7 2007/09/24 16:25:57 arce Exp $ |
---|
[850] | 27 | // GEANT4 tag $Name: HEAD $ |
---|
[815] | 28 | // |
---|
| 29 | // ------------------------------------------------------------ |
---|
| 30 | // GEANT 4 class implementation file |
---|
| 31 | // ------------------------------------------------------------ |
---|
| 32 | // |
---|
| 33 | #include "G4ErrorFreeTrajState.hh" |
---|
| 34 | #include "G4ErrorFreeTrajParam.hh" |
---|
| 35 | #include "G4ErrorSurfaceTrajState.hh" |
---|
| 36 | |
---|
| 37 | #include "G4ErrorMatrix.hh" |
---|
| 38 | #include <iomanip> |
---|
| 39 | |
---|
| 40 | #include "G4Field.hh" |
---|
| 41 | #include "G4FieldManager.hh" |
---|
| 42 | #include "G4TransportationManager.hh" |
---|
| 43 | #include "G4GeometryTolerance.hh" |
---|
| 44 | #include "G4Material.hh" |
---|
| 45 | #include "G4ErrorPropagatorData.hh" |
---|
| 46 | |
---|
| 47 | //------------------------------------------------------------------------ |
---|
| 48 | G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partType, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partType, pos, mom, errmat ) |
---|
| 49 | { |
---|
| 50 | fTrajParam = G4ErrorFreeTrajParam( pos, mom ); |
---|
| 51 | Init(); |
---|
| 52 | } |
---|
| 53 | |
---|
| 54 | |
---|
| 55 | //------------------------------------------------------------------------ |
---|
| 56 | G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() ) |
---|
| 57 | { |
---|
| 58 | // G4ThreeVector planeNormal = tpSD.GetPlaneNormal(); |
---|
| 59 | // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to plane |
---|
| 60 | // G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters(); |
---|
| 61 | // G4ThreeVector Psc = fPt * planeNormal + tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW(); |
---|
| 62 | |
---|
| 63 | fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum ); |
---|
| 64 | Init(); |
---|
| 65 | |
---|
| 66 | //----- Get the error matrix in SC coordinates |
---|
| 67 | G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters(); |
---|
| 68 | G4double mom = fMomentum.mag(); |
---|
| 69 | G4double mom2 = fMomentum.mag2(); |
---|
| 70 | G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPV()*tpSDparam.GetPV()) ); |
---|
| 71 | G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 ); |
---|
| 72 | G4Vector3D vectorU = tpSDparam.GetVectorV().cross( tpSDparam.GetVectorW() ); |
---|
| 73 | G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW(); |
---|
| 74 | |
---|
| 75 | #ifdef G4EVERBOSE |
---|
| 76 | if( iverbose >= 5){ |
---|
| 77 | G4double pc2 = std::asin( vTN.z() ); |
---|
| 78 | G4double pc3 = std::atan (vTN.y()/vTN.x()); |
---|
| 79 | |
---|
| 80 | G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() << " diff " << pc2-GetParameters().GetLambda() << G4endl; |
---|
| 81 | G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() << " diff " << pc3-GetParameters().GetPhi() << G4endl; |
---|
| 82 | } |
---|
| 83 | #endif |
---|
| 84 | |
---|
| 85 | //--- Get the unit vectors perp to P |
---|
| 86 | G4double cosl = std::cos( GetParameters().GetLambda() ); |
---|
| 87 | if (cosl < 1.E-30) cosl = 1.E-30; |
---|
| 88 | G4double cosl1 = 1./cosl; |
---|
| 89 | G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. ); |
---|
| 90 | G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl ); |
---|
| 91 | |
---|
| 92 | G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.); |
---|
| 93 | G4Vector3D vVperp = vUperp.cross( fMomentum ); |
---|
| 94 | vUperp *= 1./vUperp.mag(); |
---|
| 95 | vVperp *= 1./vVperp.mag(); |
---|
| 96 | |
---|
| 97 | #ifdef G4EVERBOSE |
---|
| 98 | if( iverbose >= 5){ |
---|
| 99 | G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff " << (vUN-vUperp).mag() << G4endl; |
---|
| 100 | G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff " << (vVN-vVperp).mag() << G4endl; |
---|
| 101 | } |
---|
| 102 | #endif |
---|
| 103 | |
---|
| 104 | //get the dot products of vectors perpendicular to direction and vector defining SD plane |
---|
| 105 | G4double dUU = vUperp * tpSD.GetVectorV(); |
---|
| 106 | G4double dUV = vUperp * tpSD.GetVectorW(); |
---|
| 107 | G4double dVU = vVperp * tpSD.GetVectorV(); |
---|
| 108 | G4double dVV = vVperp * tpSD.GetVectorW(); |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | //--- Get transformation first |
---|
| 112 | G4ErrorMatrix transfM(5, 5, 1 ); |
---|
| 113 | //--- Get magnetic field |
---|
| 114 | const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); |
---|
| 115 | G4ThreeVector dir = fTrajParam.GetDirection(); |
---|
| 116 | G4double invCosTheta = 1./std::cos( dir.theta() ); |
---|
| 117 | |
---|
| 118 | if( fCharge != 0 |
---|
| 119 | && field ) { |
---|
| 120 | G4double pos1[3]; pos1[0] = fPosition.x()*cm; pos1[1] = fPosition.y()*cm; pos1[2] = fPosition.z()*cm; |
---|
| 121 | G4double h1[3]; |
---|
| 122 | field->GetFieldValue( pos1, h1 ); |
---|
| 123 | G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; |
---|
| 124 | G4double magHPre = HPre.mag(); |
---|
| 125 | G4double invP = 1./fMomentum.mag(); |
---|
| 126 | G4double magHPreM = magHPre * invP; |
---|
| 127 | if( magHPre != 0. ) { |
---|
| 128 | G4double magHPreM2 = fCharge / magHPre; |
---|
| 129 | |
---|
| 130 | G4double Q = -magHPreM * c_light; |
---|
| 131 | G4double sinz = -HPre*vUperp * magHPreM2; |
---|
| 132 | G4double cosz = HPre*vVperp * magHPreM2; |
---|
| 133 | |
---|
| 134 | transfM[1][3] = -Q*dir.y()*sinz; |
---|
| 135 | transfM[1][4] = -Q*dir.z()*sinz; |
---|
| 136 | transfM[2][3] = -Q*dir.y()*cosz*invCosTheta; |
---|
| 137 | transfM[2][4] = -Q*dir.z()*cosz*invCosTheta; |
---|
| 138 | } |
---|
| 139 | } |
---|
| 140 | |
---|
| 141 | transfM[0][0] = 1.; |
---|
| 142 | transfM[1][1] = dir.x()*dVU; |
---|
| 143 | transfM[1][2] = dir.x()*dVV; |
---|
| 144 | transfM[2][1] = dir.x()*dUU*invCosTheta; |
---|
| 145 | transfM[2][2] = dir.x()*dUV*invCosTheta; |
---|
| 146 | transfM[3][3] = dUU; |
---|
| 147 | transfM[3][4] = dUV; |
---|
| 148 | transfM[4][3] = dVU; |
---|
| 149 | transfM[4][4] = dVV; |
---|
| 150 | |
---|
| 151 | fError = G4ErrorTrajErr( tpSD.GetError().similarity( transfM ) ); |
---|
| 152 | |
---|
| 153 | #ifdef G4EVERBOSE |
---|
| 154 | if( iverbose >= 1) G4cout << "error matrix SD2SC " << fError << G4endl; |
---|
| 155 | if( iverbose >= 4) G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl; |
---|
| 156 | #endif |
---|
| 157 | } |
---|
| 158 | |
---|
| 159 | |
---|
| 160 | //------------------------------------------------------------------------ |
---|
| 161 | void G4ErrorFreeTrajState::Init() |
---|
| 162 | { |
---|
| 163 | theTSType = G4eTS_FREE; |
---|
| 164 | BuildCharge(); |
---|
| 165 | theTransfMat = G4ErrorMatrix(5,5,0); |
---|
| 166 | //- theFirstStep = true; |
---|
| 167 | } |
---|
| 168 | |
---|
| 169 | //------------------------------------------------------------------------ |
---|
| 170 | void G4ErrorFreeTrajState::Dump( std::ostream& out ) const |
---|
| 171 | { |
---|
| 172 | out << *this; |
---|
| 173 | } |
---|
| 174 | |
---|
| 175 | //------------------------------------------------------------------------ |
---|
| 176 | G4int G4ErrorFreeTrajState::Update( const G4Track* aTrack ) |
---|
| 177 | { |
---|
| 178 | G4int ierr = 0; |
---|
| 179 | fTrajParam.Update( aTrack ); |
---|
| 180 | UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() ); |
---|
| 181 | return ierr; |
---|
| 182 | |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | |
---|
| 186 | //------------------------------------------------------------------------ |
---|
| 187 | std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts) |
---|
| 188 | { |
---|
| 189 | out.setf(std::ios::fixed,std::ios::floatfield); |
---|
| 190 | |
---|
| 191 | |
---|
| 192 | ts.DumpPosMomError( out ); |
---|
| 193 | |
---|
| 194 | out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl; |
---|
| 195 | |
---|
| 196 | return out; |
---|
| 197 | |
---|
| 198 | } |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | //------------------------------------------------------------------------ |
---|
| 202 | G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack ) |
---|
| 203 | { |
---|
| 204 | G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; |
---|
| 205 | G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
| 206 | |
---|
| 207 | if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0; |
---|
| 208 | |
---|
| 209 | #ifdef G4EVERBOSE |
---|
| 210 | if( iverbose >= 2 )G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl; |
---|
| 211 | #endif |
---|
| 212 | |
---|
| 213 | // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES |
---|
| 214 | G4Point3D vposPost = aTrack->GetPosition()/cm; |
---|
| 215 | G4Vector3D vpPost = aTrack->GetMomentum()/GeV; |
---|
| 216 | // G4Point3D vposPre = fPosition/cm; |
---|
| 217 | // G4Vector3D vpPre = fMomentum/GeV; |
---|
| 218 | G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm; |
---|
| 219 | G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV; |
---|
| 220 | //correct to avoid propagation along Z |
---|
| 221 | if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV ); |
---|
| 222 | if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV ); |
---|
| 223 | |
---|
| 224 | G4double pPre = vpPre.mag(); |
---|
| 225 | G4double pPost = vpPost.mag(); |
---|
| 226 | #ifdef G4EVERBOSE |
---|
| 227 | if( iverbose >= 2 ) { |
---|
| 228 | G4cout << "G4EP: vposPre " << vposPre << G4endl |
---|
| 229 | << "G4EP: vposPost " << vposPost << G4endl; |
---|
| 230 | G4cout << "G4EP: vpPre " << vpPre << G4endl |
---|
| 231 | << "G4EP: vpPost " << vpPost << G4endl; |
---|
| 232 | G4cout << " err start step " << fError << G4endl; |
---|
| 233 | G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl; |
---|
| 234 | } |
---|
| 235 | #endif |
---|
| 236 | |
---|
| 237 | if( pPre == 0. || pPost == 0 ) return 2; |
---|
| 238 | G4double pInvPre = 1./pPre; |
---|
| 239 | G4double pInvPost = 1./pPost; |
---|
| 240 | G4double deltaPInv = pInvPost - pInvPre; |
---|
| 241 | |
---|
| 242 | G4Vector3D vpPreNorm = vpPre * pInvPre; |
---|
| 243 | G4Vector3D vpPostNorm = vpPost * pInvPost; |
---|
| 244 | // if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl; |
---|
| 245 | //return if propagation along Z?? |
---|
| 246 | if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4; |
---|
| 247 | G4double sinpPre = std::sin( vpPreNorm.theta() ); //cosine perpendicular to pPre = sine pPre |
---|
| 248 | G4double sinpPost = std::sin( vpPostNorm.theta() ); //cosine perpendicular to pPost = sine pPost |
---|
| 249 | G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() ); |
---|
| 250 | |
---|
| 251 | #ifdef G4EVERBOSE |
---|
| 252 | if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl; |
---|
| 253 | #endif |
---|
| 254 | //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR |
---|
| 255 | //* *** NEUTRAL PARTICLE OR FIELDFREE REGION |
---|
| 256 | G4ErrorMatrix transf(5, 5, 0 ); |
---|
| 257 | |
---|
| 258 | transf[3][2] = stepLengthCm * sinpPost; |
---|
| 259 | transf[4][1] = stepLengthCm; |
---|
| 260 | for( size_t ii=0;ii < 5; ii++ ){ |
---|
| 261 | transf[ii][ii] = 1.; |
---|
| 262 | } |
---|
| 263 | #ifdef G4EVERBOSE |
---|
| 264 | if( iverbose >= 2 ) { |
---|
| 265 | G4cout << "G4EP: transf matrix neutral " << transf; |
---|
| 266 | } |
---|
| 267 | #endif |
---|
| 268 | |
---|
| 269 | // charge X propagation direction |
---|
| 270 | G4double charge = aTrack->GetDynamicParticle()->GetCharge(); |
---|
| 271 | if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ) { |
---|
| 272 | charge *= -1.; |
---|
| 273 | } |
---|
| 274 | // G4cout << " charge " << charge << G4endl; |
---|
| 275 | //t check if particle has charge |
---|
| 276 | //t if( charge == 0 ) goto 45; |
---|
| 277 | // check if the magnetic field is = 0. |
---|
| 278 | |
---|
| 279 | //position is from geant4, it is assumed to be in mm (for debugging, eventually it will not be transformed) |
---|
| 280 | G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm; |
---|
| 281 | G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm; |
---|
| 282 | G4double h1[3], h2[3]; |
---|
| 283 | |
---|
| 284 | const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); |
---|
| 285 | if( !field ) return 0; //goto 45 |
---|
| 286 | |
---|
| 287 | // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION |
---|
| 288 | if( charge != 0. && field ) { |
---|
| 289 | |
---|
| 290 | field->GetFieldValue( pos1, h1 ); |
---|
| 291 | field->GetFieldValue( pos2, h2 ); |
---|
| 292 | G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; //10. is to get same dimensions as GEANT3 (kilogauss) |
---|
| 293 | G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.; |
---|
| 294 | G4double magHPre = HPre.mag(); |
---|
| 295 | G4double magHPost = HPost.mag(); |
---|
| 296 | #ifdef G4EVERBOSE |
---|
| 297 | if( iverbose >= 2 ) G4cout << "G4EP: HPre " << HPre << G4endl |
---|
| 298 | << "G4EP: HPost " << HPost << G4endl; |
---|
| 299 | #endif |
---|
| 300 | |
---|
| 301 | if( magHPre + magHPost != 0. ) { |
---|
| 302 | |
---|
| 303 | //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2 |
---|
| 304 | G4double gam; |
---|
| 305 | if( magHPost != 0. ){ |
---|
| 306 | gam = HPost * vpPostNorm / magHPost; |
---|
| 307 | }else { |
---|
| 308 | gam = HPre * vpPreNorm / magHPre; |
---|
| 309 | } |
---|
| 310 | |
---|
| 311 | // G4eMagneticLimitsProcess will limit the step, but based on an straight line trajectory |
---|
| 312 | G4double alphaSqr = 1. - gam * gam; |
---|
| 313 | G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2(); |
---|
| 314 | G4double delhp6Sqr = 300.*300.; |
---|
| 315 | #ifdef G4EVERBOSE |
---|
| 316 | if( iverbose >= 2 ) G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr << " diffHSqr " << diffHSqr << G4endl; |
---|
| 317 | #endif |
---|
| 318 | if( diffHSqr * alphaSqr > delhp6Sqr ) return 3; |
---|
| 319 | |
---|
| 320 | |
---|
| 321 | //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT |
---|
| 322 | G4double pInvAver = 1./(pInvPre + pInvPost ); |
---|
| 323 | G4double CFACT8 = 2.997925E-4; |
---|
| 324 | //G4double HAver |
---|
| 325 | G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 ); |
---|
| 326 | G4double HAver = vHAverNorm.mag(); |
---|
| 327 | G4double invHAver = 1./HAver; |
---|
| 328 | vHAverNorm *= invHAver; |
---|
| 329 | #ifdef G4EVERBOSE |
---|
| 330 | if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl; |
---|
| 331 | #endif |
---|
| 332 | |
---|
| 333 | G4double pAver = (pPre+pPost)*0.5; |
---|
| 334 | G4double QAver = -HAver/pAver; |
---|
| 335 | G4double thetaAver = QAver * stepLengthCm; |
---|
| 336 | G4double sinThetaAver = std::sin(thetaAver); |
---|
| 337 | G4double cosThetaAver = std::cos(thetaAver); |
---|
| 338 | G4double gamma = vHAverNorm * vpPostNorm; |
---|
| 339 | G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm ); |
---|
| 340 | |
---|
| 341 | #ifdef G4EVERBOSE |
---|
| 342 | if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << G4endl; |
---|
| 343 | #endif |
---|
| 344 | G4double AU = 1./vpPreNorm.perp(); |
---|
| 345 | //t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU ); |
---|
| 346 | G4ThreeVector vUPre( -AU*vpPreNorm.y(), |
---|
| 347 | AU*vpPreNorm.x(), |
---|
| 348 | 0. ); |
---|
| 349 | G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(), |
---|
| 350 | vpPreNorm.z()*vUPre.x(), |
---|
| 351 | vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() ); |
---|
| 352 | |
---|
| 353 | // |
---|
| 354 | AU = 1./vpPostNorm.perp(); |
---|
| 355 | //t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU ); |
---|
| 356 | G4ThreeVector vUPost( -AU*vpPostNorm.y(), |
---|
| 357 | AU*vpPostNorm.x(), |
---|
| 358 | 0. ); |
---|
| 359 | G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(), |
---|
| 360 | vpPostNorm.z()*vUPost.x(), |
---|
| 361 | vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() ); |
---|
| 362 | #ifdef G4EVERBOSE |
---|
| 363 | //- G4cout << " vpPostNorm " << vpPostNorm << G4endl; |
---|
| 364 | if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl; |
---|
| 365 | #endif |
---|
| 366 | G4Point3D deltaPos( vposPre - vposPost ); |
---|
| 367 | |
---|
| 368 | // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2 |
---|
| 369 | // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT |
---|
| 370 | // * *** TAKEN INTO ACCOUNT |
---|
| 371 | |
---|
| 372 | G4double QP = QAver * pAver; // = -HAver |
---|
| 373 | #ifdef G4EVERBOSE |
---|
| 374 | if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl; |
---|
| 375 | #endif |
---|
| 376 | G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() ); |
---|
| 377 | G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() ); |
---|
| 378 | G4double OMcosThetaAver = 1. - cosThetaAver; |
---|
| 379 | #ifdef G4EVERBOSE |
---|
| 380 | if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl; |
---|
| 381 | #endif |
---|
| 382 | G4double TMSINT = thetaAver - sinThetaAver; |
---|
| 383 | #ifdef G4EVERBOSE |
---|
| 384 | if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl; |
---|
| 385 | #endif |
---|
| 386 | |
---|
| 387 | G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(), |
---|
| 388 | vHAverNorm.z() * vUPre.x(), |
---|
| 389 | vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() ); |
---|
| 390 | #ifdef G4EVERBOSE |
---|
| 391 | // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " " << vHAverNorm.z() << " " << vUPre.y() << G4endl; |
---|
| 392 | #endif |
---|
| 393 | G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(), |
---|
| 394 | vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(), |
---|
| 395 | vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() ); |
---|
| 396 | #ifdef G4EVERBOSE |
---|
| 397 | if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl; |
---|
| 398 | #endif |
---|
| 399 | |
---|
| 400 | //------------------- COMPUTE MATRIX |
---|
| 401 | //---------- 1/P |
---|
| 402 | |
---|
| 403 | transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm) |
---|
| 404 | +2.*deltaPInv*pAver; |
---|
| 405 | |
---|
| 406 | transf[0][1] = -deltaPInv/thetaAver* |
---|
| 407 | ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) + |
---|
| 408 | sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + |
---|
| 409 | OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) ); |
---|
| 410 | |
---|
| 411 | transf[0][2] = -sinpPre*deltaPInv/thetaAver* |
---|
| 412 | ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) + |
---|
| 413 | sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + |
---|
| 414 | OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) ); |
---|
| 415 | |
---|
| 416 | transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ); |
---|
| 417 | |
---|
| 418 | transf[0][4] = -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()); |
---|
| 419 | |
---|
| 420 | // *** Lambda |
---|
| 421 | transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z()) |
---|
| 422 | *(1.+deltaPInv*pAver); |
---|
| 423 | #ifdef G4EVERBOSE |
---|
| 424 | if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0] << " " << -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z() |
---|
| 425 | << " " << deltaPInv<< " " << pAver << G4endl; |
---|
| 426 | #endif |
---|
| 427 | |
---|
| 428 | transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + |
---|
| 429 | sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + |
---|
| 430 | OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* |
---|
| 431 | (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + |
---|
| 432 | ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + |
---|
| 433 | OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - |
---|
| 434 | TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) ); |
---|
| 435 | |
---|
| 436 | transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + |
---|
| 437 | sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + |
---|
| 438 | OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* |
---|
| 439 | (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + |
---|
| 440 | ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + |
---|
| 441 | OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - |
---|
| 442 | TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) ); |
---|
| 443 | transf[1][2] = sinpPre*transf[1][2]; |
---|
| 444 | |
---|
| 445 | transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ); |
---|
| 446 | |
---|
| 447 | transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()); |
---|
| 448 | |
---|
| 449 | // *** Phi |
---|
| 450 | |
---|
| 451 | transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv |
---|
| 452 | *(1.+deltaPInv*pAver); |
---|
| 453 | #ifdef G4EVERBOSE |
---|
| 454 | if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv |
---|
| 455 | <<" "<<deltaPInv<<" "<<pAver<< G4endl; |
---|
| 456 | #endif |
---|
| 457 | transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) + |
---|
| 458 | sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) + |
---|
| 459 | OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* |
---|
| 460 | (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) + |
---|
| 461 | ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + |
---|
| 462 | OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - |
---|
| 463 | TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) ); |
---|
| 464 | transf[2][1] = sinpPostInv*transf[2][1]; |
---|
| 465 | |
---|
| 466 | transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) + |
---|
| 467 | sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) + |
---|
| 468 | OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* |
---|
| 469 | (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) + |
---|
| 470 | ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + |
---|
| 471 | OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - |
---|
| 472 | TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) ); |
---|
| 473 | transf[2][2] = sinpPostInv*sinpPre*transf[2][2]; |
---|
| 474 | |
---|
| 475 | transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv; |
---|
| 476 | #ifdef G4EVERBOSE |
---|
| 477 | if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl; |
---|
| 478 | #endif |
---|
| 479 | |
---|
| 480 | transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv; |
---|
| 481 | |
---|
| 482 | // *** Yt |
---|
| 483 | |
---|
| 484 | transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() ) |
---|
| 485 | *(1.+deltaPInv*pAver); |
---|
| 486 | #ifdef G4EVERBOSE |
---|
| 487 | if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y() |
---|
| 488 | <<" "<<deltaPInv<<" "<<pAver<<G4endl; |
---|
| 489 | #endif |
---|
| 490 | |
---|
| 491 | transf[3][1] = ( sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) + |
---|
| 492 | OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) + |
---|
| 493 | TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )* |
---|
| 494 | (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver; |
---|
| 495 | |
---|
| 496 | transf[3][2] = ( sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) + |
---|
| 497 | OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) + |
---|
| 498 | TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )* |
---|
| 499 | (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver; |
---|
| 500 | #ifdef G4EVERBOSE |
---|
| 501 | if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<< |
---|
| 502 | OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<< |
---|
| 503 | TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<< |
---|
| 504 | vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl; |
---|
| 505 | #endif |
---|
| 506 | |
---|
| 507 | transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ); |
---|
| 508 | |
---|
| 509 | transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ); |
---|
| 510 | |
---|
| 511 | // *** Zt |
---|
| 512 | transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z()) |
---|
| 513 | *(1.+deltaPInv*pAver); |
---|
| 514 | |
---|
| 515 | transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + |
---|
| 516 | OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + |
---|
| 517 | TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* |
---|
| 518 | (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver; |
---|
| 519 | #ifdef G4EVERBOSE |
---|
| 520 | if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl; |
---|
| 521 | #endif |
---|
| 522 | |
---|
| 523 | transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + |
---|
| 524 | OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + |
---|
| 525 | TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* |
---|
| 526 | (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver; |
---|
| 527 | |
---|
| 528 | transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ); |
---|
| 529 | |
---|
| 530 | transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()); |
---|
| 531 | // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<< vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<" "<< vVPre.z() <<" "<< vVPost.z() << G4endl; |
---|
| 532 | |
---|
| 533 | |
---|
| 534 | #ifdef G4EVERBOSE |
---|
| 535 | if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl; |
---|
| 536 | #endif |
---|
| 537 | /* for( G4int ii=0;ii<5;ii++){ |
---|
| 538 | for( G4int jj=0;jj<5;jj++){ |
---|
| 539 | G4cout << transf[ii][jj] << " "; |
---|
| 540 | } |
---|
| 541 | G4cout << G4endl; |
---|
| 542 | } */ |
---|
| 543 | } |
---|
| 544 | } |
---|
| 545 | // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION |
---|
| 546 | /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized " << theFirstStep << G4endl; |
---|
| 547 | if( theFirstStep ) { |
---|
| 548 | theTransfMat = transf; |
---|
| 549 | theFirstStep = false; |
---|
| 550 | }else{ |
---|
| 551 | theTransfMat = theTransfMat * transf; |
---|
| 552 | if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" << theTransfMat << G4endl; |
---|
| 553 | } |
---|
| 554 | */ |
---|
| 555 | theTransfMat = transf; |
---|
| 556 | #ifdef G4EVERBOSE |
---|
| 557 | if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl; |
---|
| 558 | if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl |
---|
| 559 | << " transf matrix " << theTransfMat.T() << G4endl; |
---|
| 560 | #endif |
---|
| 561 | |
---|
| 562 | fError = fError.similarity(theTransfMat).T(); |
---|
| 563 | //- fError = transf * fError * transf.T(); |
---|
| 564 | #ifdef G4EVERBOSE |
---|
| 565 | if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl; |
---|
| 566 | #endif |
---|
| 567 | |
---|
| 568 | //? S = B*S*BT S.similarity(B) |
---|
| 569 | //? R = S |
---|
| 570 | // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES; |
---|
| 571 | |
---|
| 572 | PropagateErrorMSC( aTrack ); |
---|
| 573 | |
---|
| 574 | PropagateErrorIoni( aTrack ); |
---|
| 575 | |
---|
| 576 | return 0; |
---|
| 577 | } |
---|
| 578 | |
---|
| 579 | |
---|
| 580 | //------------------------------------------------------------------------ |
---|
| 581 | G4int G4ErrorFreeTrajState::PropagateErrorMSC( const G4Track* aTrack ) |
---|
| 582 | { |
---|
| 583 | G4ThreeVector vpPre = aTrack->GetMomentum()/GeV; |
---|
| 584 | G4double pPre = vpPre.mag(); |
---|
| 585 | G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV); |
---|
| 586 | G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; |
---|
| 587 | |
---|
| 588 | G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial(); |
---|
| 589 | G4double effZ, effA; |
---|
| 590 | CalculateEffectiveZandA( mate, effZ, effA ); |
---|
| 591 | |
---|
| 592 | #ifdef G4EVERBOSE |
---|
| 593 | if( iverbose >= 4 ) G4cout << "material " << mate->GetName() |
---|
| 594 | //<< " " << mate->GetZ() << " " << mate->GetA() |
---|
| 595 | << " " << effZ << " " << effA |
---|
| 596 | << " " << mate->GetDensity()/g*mole << " " << mate->GetRadlen()/cm << " " << mate->GetNuclearInterLength()/cm << G4endl; |
---|
| 597 | #endif |
---|
| 598 | |
---|
| 599 | G4double RI = stepLengthCm / (mate->GetRadlen()/cm); |
---|
| 600 | #ifdef G4EVERBOSE |
---|
| 601 | if( iverbose >= 4 ) G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI " << RI << " stepLengthCm " << stepLengthCm << " radlen " << (mate->GetRadlen()/cm) << " " << RI*1.e10 << G4endl; |
---|
| 602 | #endif |
---|
| 603 | G4double charge = aTrack->GetDynamicParticle()->GetCharge(); |
---|
| 604 | G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta ); |
---|
| 605 | #ifdef G4EVERBOSE |
---|
| 606 | if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl; |
---|
| 607 | #endif |
---|
| 608 | G4double S1 = DD*stepLengthCm*stepLengthCm/3.; |
---|
| 609 | G4double S2 = DD; |
---|
| 610 | G4double S3 = DD*stepLengthCm/2.; |
---|
| 611 | |
---|
| 612 | G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre; |
---|
| 613 | #ifdef G4EVERBOSE |
---|
| 614 | if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 " << S2 << " S3 " << S3 << " CLA " << CLA << G4endl; |
---|
| 615 | #endif |
---|
| 616 | fError[1][1] += S2; |
---|
| 617 | fError[1][4] -= S3; |
---|
| 618 | fError[2][2] += S2/CLA/CLA; |
---|
| 619 | fError[2][3] += S3/CLA; |
---|
| 620 | fError[3][3] += S1; |
---|
| 621 | fError[4][4] += S1; |
---|
| 622 | |
---|
| 623 | #ifdef G4EVERBOSE |
---|
| 624 | if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl; |
---|
| 625 | #endif |
---|
| 626 | |
---|
| 627 | return 0; |
---|
| 628 | } |
---|
| 629 | |
---|
| 630 | |
---|
| 631 | //------------------------------------------------------------------------ |
---|
| 632 | void G4ErrorFreeTrajState::CalculateEffectiveZandA( const G4Material* mate, G4double& effZ, G4double& effA ) |
---|
| 633 | { |
---|
| 634 | effZ = 0.; |
---|
| 635 | effA = 0.; |
---|
| 636 | G4int ii, nelem = mate->GetNumberOfElements(); |
---|
| 637 | const G4double* fracVec = mate->GetFractionVector(); |
---|
| 638 | for(ii=0; ii < nelem; ii++ ) { |
---|
| 639 | effZ += mate->GetElement( ii )->GetZ() * fracVec[ii]; |
---|
| 640 | effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole; |
---|
| 641 | } |
---|
| 642 | |
---|
| 643 | } |
---|
| 644 | |
---|
| 645 | |
---|
| 646 | //------------------------------------------------------------------------ |
---|
| 647 | G4int G4ErrorFreeTrajState::PropagateErrorIoni( const G4Track* aTrack ) |
---|
| 648 | { |
---|
| 649 | G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; |
---|
| 650 | G4double DEDX2; |
---|
| 651 | if( stepLengthCm < 1.E-7 ) { |
---|
| 652 | DEDX2=0.; |
---|
| 653 | } |
---|
| 654 | // * Calculate xi factor (KeV). |
---|
| 655 | G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial(); |
---|
| 656 | G4double effZ, effA; |
---|
| 657 | CalculateEffectiveZandA( mate, effZ, effA ); |
---|
| 658 | |
---|
| 659 | G4double Etot = aTrack->GetTotalEnergy()/GeV; |
---|
| 660 | G4double beta = aTrack->GetMomentum().mag()/GeV / Etot; |
---|
| 661 | G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV; |
---|
| 662 | G4double gamma = Etot / mass; |
---|
| 663 | |
---|
| 664 | // * Calculate xi factor (KeV). |
---|
| 665 | G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) / |
---|
| 666 | (effA*beta*beta); |
---|
| 667 | |
---|
| 668 | #ifdef G4EVERBOSE |
---|
| 669 | if( iverbose >= 2 ){ |
---|
| 670 | G4cout << "G4EP:IONI: XI " << XI << " beta " << beta << " gamma " << gamma << G4endl; |
---|
| 671 | G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl; |
---|
| 672 | } |
---|
| 673 | #endif |
---|
| 674 | // * Maximum energy transfer to atomic electron (KeV). |
---|
| 675 | G4double eta = beta*gamma; |
---|
| 676 | G4double etasq = eta*eta; |
---|
| 677 | G4double eMass = 0.51099906/GeV; |
---|
| 678 | G4double massRatio = eMass / mass; |
---|
| 679 | G4double F1 = 2*eMass*etasq; |
---|
| 680 | G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio; |
---|
| 681 | G4double Emax = 1.E+6*F1/F2; |
---|
| 682 | |
---|
| 683 | // * *** and now sigma**2 in GeV |
---|
| 684 | G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12; |
---|
| 685 | #ifdef G4EVERBOSE |
---|
| 686 | if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX2 " << dedxSq << " emass " << eMass << " Emax " << Emax << G4endl; |
---|
| 687 | #endif |
---|
| 688 | |
---|
| 689 | // if( iverbose >= 2 ) G4cout << "G4EP:IONI: Etot " << Etot << " DEDX2 " << dedxSq << " emass " << eMass << G4endl; |
---|
| 690 | |
---|
| 691 | G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag(); |
---|
| 692 | pPre6 = std::pow(pPre6, 6 ); |
---|
| 693 | //Apply it to error |
---|
| 694 | fError[0][0] += Etot*Etot*dedxSq / pPre6; |
---|
| 695 | #ifdef G4EVERBOSE |
---|
| 696 | if( iverbose >= 2 ) G4cout << "G4:IONI getot " << Etot << " dedx2 " << dedxSq << " p " << pPre6 << G4endl; |
---|
| 697 | if( iverbose >= 2 ) G4cout << "G4EP:IONI: error_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl; |
---|
| 698 | #endif |
---|
| 699 | |
---|
| 700 | return 0; |
---|
| 701 | } |
---|
| 702 | |
---|