00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include <iomanip>
00034
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4Field.hh"
00038 #include "G4FieldManager.hh"
00039 #include "G4TransportationManager.hh"
00040 #include "G4GeometryTolerance.hh"
00041 #include "G4Material.hh"
00042 #include "G4ErrorPropagatorData.hh"
00043 #include "G4ErrorFreeTrajState.hh"
00044 #include "G4ErrorFreeTrajParam.hh"
00045 #include "G4ErrorSurfaceTrajState.hh"
00046 #include "G4ErrorMatrix.hh"
00047
00048
00049 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partName, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partName, pos, mom, errmat )
00050 {
00051 fTrajParam = G4ErrorFreeTrajParam( pos, mom );
00052 Init();
00053 }
00054
00055
00056
00057 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() )
00058 {
00059
00060
00061
00062
00063
00064 fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum );
00065 Init();
00066
00067
00068 G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
00069 G4double mom = fMomentum.mag();
00070 G4double mom2 = fMomentum.mag2();
00071 G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPW()*tpSDparam.GetPW()) );
00072 G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 );
00073 G4Vector3D vectorU = tpSDparam.GetVectorV().cross( tpSDparam.GetVectorW() );
00074 G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW();
00075
00076 #ifdef G4EVERBOSE
00077 if( iverbose >= 5){
00078 G4double pc2 = std::asin( vTN.z() );
00079 G4double pc3 = std::atan (vTN.y()/vTN.x());
00080
00081 G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() << " diff " << pc2-GetParameters().GetLambda() << G4endl;
00082 G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() << " diff " << pc3-GetParameters().GetPhi() << G4endl;
00083 }
00084 #endif
00085
00086
00087 G4double cosl = std::cos( GetParameters().GetLambda() );
00088 if (cosl < 1.E-30) cosl = 1.E-30;
00089 G4double cosl1 = 1./cosl;
00090 G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. );
00091 G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl );
00092
00093 G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.);
00094 G4Vector3D vVperp = vUperp.cross( fMomentum );
00095 vUperp *= 1./vUperp.mag();
00096 vVperp *= 1./vVperp.mag();
00097
00098 #ifdef G4EVERBOSE
00099 if( iverbose >= 5 ){
00100 G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff " << (vUN-vUperp).mag() << G4endl;
00101 G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff " << (vVN-vVperp).mag() << G4endl;
00102 }
00103 #endif
00104
00105
00106 G4double dUU = vUperp * tpSD.GetVectorV();
00107 G4double dUV = vUperp * tpSD.GetVectorW();
00108 G4double dVU = vVperp * tpSD.GetVectorV();
00109 G4double dVV = vVperp * tpSD.GetVectorW();
00110
00111
00112 G4ErrorMatrix transfM(5, 5, 1 );
00113
00114 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
00115 G4ThreeVector dir = fTrajParam.GetDirection();
00116 G4double invCosTheta = 1./std::cos( dir.theta() );
00117 G4cout << " dir="<<dir<<" invCosTheta "<<invCosTheta << G4endl;
00118
00119 if( fCharge != 0 && field ) {
00120 G4double pos1[3]; pos1[0] = fPosition.x()*cm; pos1[1] = fPosition.y()*cm; pos1[2] = fPosition.z()*cm;
00121 G4double h1[3];
00122 field->GetFieldValue( pos1, h1 );
00123 G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.;
00124 G4double magHPre = HPre.mag();
00125 G4double invP = 1./fMomentum.mag();
00126 G4double magHPreM = magHPre * invP;
00127 if( magHPre != 0. ) {
00128 G4double magHPreM2 = fCharge / magHPre;
00129
00130 G4double Q = -magHPreM * c_light;
00131 G4double sinz = -HPre*vUperp * magHPreM2;
00132 G4double cosz = HPre*vVperp * magHPreM2;
00133
00134 transfM[1][3] = -Q*dir.y()*sinz;
00135 transfM[1][4] = -Q*dir.z()*sinz;
00136 transfM[2][3] = -Q*dir.y()*cosz*invCosTheta;
00137 transfM[2][4] = -Q*dir.z()*cosz*invCosTheta;
00138 }
00139 }
00140
00141 transfM[0][0] = 1.;
00142 transfM[1][1] = dir.x()*dVU;
00143 transfM[1][2] = dir.x()*dVV;
00144 transfM[2][1] = dir.x()*dUU*invCosTheta;
00145 transfM[2][2] = dir.x()*dUV*invCosTheta;
00146 transfM[3][3] = dUU;
00147 transfM[3][4] = dUV;
00148 transfM[4][3] = dVU;
00149 transfM[4][4] = dVV;
00150
00151 fError = G4ErrorTrajErr( tpSD.GetError().similarity( transfM ) );
00152
00153 #ifdef G4EVERBOSE
00154 if( iverbose >= 1) G4cout << "error matrix SD2SC " << fError << G4endl;
00155 if( iverbose >= 4) G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl;
00156 #endif
00157 }
00158
00159
00160
00161 void G4ErrorFreeTrajState::Init()
00162 {
00163 theTSType = G4eTS_FREE;
00164 BuildCharge();
00165 theTransfMat = G4ErrorMatrix(5,5,0);
00166 theFirstStep = true;
00167 }
00168
00169
00170 void G4ErrorFreeTrajState::Dump( std::ostream& out ) const
00171 {
00172 out << *this;
00173 }
00174
00175
00176 G4int G4ErrorFreeTrajState::Update( const G4Track* aTrack )
00177 {
00178 G4int ierr = 0;
00179 fTrajParam.Update( aTrack );
00180 UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() );
00181 return ierr;
00182 }
00183
00184
00185
00186 std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts)
00187 {
00188 std::ios::fmtflags orig_flags = out.flags();
00189
00190 out.setf(std::ios::fixed,std::ios::floatfield);
00191
00192 ts.DumpPosMomError( out );
00193
00194 out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
00195
00196 out.flags(orig_flags);
00197
00198 return out;
00199 }
00200
00201
00202
00203 G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack )
00204 {
00205 G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
00206 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() == G4ErrorStage_Deflation ) stepLengthCm *= -1.;
00207
00208 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00209
00210 if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0;
00211
00212 #ifdef G4EVERBOSE
00213 if( iverbose >= 2 )G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl;
00214 G4cout << "G4EP: iverbose="<< iverbose << G4endl;
00215 #endif
00216
00217
00218 G4Point3D vposPost = aTrack->GetPosition()/cm;
00219 G4Vector3D vpPost = aTrack->GetMomentum()/GeV;
00220
00221
00222 G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm;
00223 G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV;
00224
00225 if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV );
00226 if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV );
00227
00228 G4double pPre = vpPre.mag();
00229 G4double pPost = vpPost.mag();
00230 #ifdef G4EVERBOSE
00231 if( iverbose >= 2 ) {
00232 G4cout << "G4EP: vposPre " << vposPre << G4endl
00233 << "G4EP: vposPost " << vposPost << G4endl;
00234 G4cout << "G4EP: vpPre " << vpPre << G4endl
00235 << "G4EP: vpPost " << vpPost << G4endl;
00236 G4cout << " err start step " << fError << G4endl;
00237 G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
00238 }
00239 #endif
00240
00241 if( pPre == 0. || pPost == 0 ) return 2;
00242 G4double pInvPre = 1./pPre;
00243 G4double pInvPost = 1./pPost;
00244 G4double deltaPInv = pInvPost - pInvPre;
00245 if( iverbose >= 2 ) G4cout << "G4EP: pInvPre" << pInvPre<< " pInvPost:" << pInvPost<<" deltaPInv:" << deltaPInv<< G4endl;
00246
00247
00248 G4Vector3D vpPreNorm = vpPre * pInvPre;
00249 G4Vector3D vpPostNorm = vpPost * pInvPost;
00250 if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl;
00251
00252 if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4;
00253 G4double sinpPre = std::sin( vpPreNorm.theta() );
00254 G4double sinpPost = std::sin( vpPostNorm.theta() );
00255 G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() );
00256
00257 #ifdef G4EVERBOSE
00258 if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
00259 #endif
00260
00261
00262 G4ErrorMatrix transf(5, 5, 0 );
00263
00264 transf[3][2] = stepLengthCm * sinpPost;
00265 transf[4][1] = stepLengthCm;
00266 for( size_t ii=0;ii < 5; ii++ ){
00267 transf[ii][ii] = 1.;
00268 }
00269 #ifdef G4EVERBOSE
00270 if( iverbose >= 2 ) {
00271 G4cout << "G4EP: transf matrix neutral " << transf;
00272 }
00273 #endif
00274
00275
00276 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
00277 if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ) {
00278 charge *= -1.;
00279 }
00280
00281
00282
00283
00284
00285
00286
00287 G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm;
00288 G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm;
00289 G4double h1[3], h2[3];
00290
00291 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
00292 if( !field ) return 0;
00293
00294
00295
00296
00297 if( charge != 0. && field ) {
00298
00299 field->GetFieldValue( pos1, h1 );
00300 field->GetFieldValue( pos2, h2 );
00301 G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.;
00302 G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.;
00303 G4double magHPre = HPre.mag();
00304 G4double magHPost = HPost.mag();
00305 #ifdef G4EVERBOSE
00306 if( iverbose >= 2 ) {
00307 G4cout << "G4EP: h1 = "
00308 << h1[0] << ", " << h1[1] << ", " << h1[2] << G4endl;
00309 G4cout << "G4EP: pos1/mm = "
00310 << pos1[0] << ", " << pos1[1] << ", " << pos1[2] << G4endl;
00311 G4cout << "G4EP: pos2/mm = "
00312 << pos2[0] << ", " << pos2[1] << ", " << pos2[2] << G4endl;
00313 G4cout << "G4EP: B-filed in KGauss HPre " << HPre << G4endl
00314 << "G4EP: in KGauss HPost " << HPost << G4endl;
00315 }
00316 #endif
00317
00318 if( magHPre + magHPost != 0. ) {
00319
00320
00321 G4double gam;
00322 if( magHPost != 0. ){
00323 gam = HPost * vpPostNorm / magHPost;
00324 }else {
00325 gam = HPre * vpPreNorm / magHPre;
00326 }
00327
00328
00329 G4double alphaSqr = 1. - gam * gam;
00330 G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2();
00331 G4double delhp6Sqr = 300.*300.;
00332 #ifdef G4EVERBOSE
00333 if( iverbose >= 2 ) {
00334 G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr
00335 << " diffHSqr " << diffHSqr << G4endl;
00336 G4cout << " alpha= " << sqrt(alphaSqr) << G4endl;
00337 }
00338 #endif
00339 if( diffHSqr * alphaSqr > delhp6Sqr ) return 3;
00340
00341
00342
00343 G4double pInvAver = 1./(pInvPre + pInvPost );
00344 G4double CFACT8 = 2.997925E-4;
00345
00346 G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 );
00347 G4double HAver = vHAverNorm.mag();
00348 G4double invHAver = 1./HAver;
00349 vHAverNorm *= invHAver;
00350 #ifdef G4EVERBOSE
00351 if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl;
00352 #endif
00353
00354 G4double pAver = (pPre+pPost)*0.5;
00355 G4double QAver = -HAver/pAver;
00356 G4double thetaAver = QAver * stepLengthCm;
00357 G4double sinThetaAver = std::sin(thetaAver);
00358 G4double cosThetaAver = std::cos(thetaAver);
00359 G4double gamma = vHAverNorm * vpPostNorm;
00360 G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm );
00361
00362 #ifdef G4EVERBOSE
00363 if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << " gamma:"<<gamma<< " theta="<< thetaAver<<G4endl;
00364 #endif
00365 G4double AU = 1./vpPreNorm.perp();
00366
00367 G4ThreeVector vUPre( -AU*vpPreNorm.y(),
00368 AU*vpPreNorm.x(),
00369 0. );
00370 G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(),
00371 vpPreNorm.z()*vUPre.x(),
00372 vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() );
00373
00374
00375 AU = 1./vpPostNorm.perp();
00376
00377 G4ThreeVector vUPost( -AU*vpPostNorm.y(),
00378 AU*vpPostNorm.x(),
00379 0. );
00380 G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(),
00381 vpPostNorm.z()*vUPost.x(),
00382 vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() );
00383 #ifdef G4EVERBOSE
00384 G4cout << " vpPostNorm " << vpPostNorm << G4endl;
00385 if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
00386 #endif
00387 G4Point3D deltaPos( vposPre - vposPost );
00388
00389
00390
00391
00392
00393 G4double QP = QAver * pAver;
00394 #ifdef G4EVERBOSE
00395 if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl;
00396 #endif
00397 G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() );
00398 G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() );
00399 G4double OMcosThetaAver = 1. - cosThetaAver;
00400 #ifdef G4EVERBOSE
00401 if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl;
00402 #endif
00403 G4double TMSINT = thetaAver - sinThetaAver;
00404 #ifdef G4EVERBOSE
00405 if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
00406 #endif
00407
00408 G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(),
00409 vHAverNorm.z() * vUPre.x(),
00410 vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() );
00411 #ifdef G4EVERBOSE
00412
00413 #endif
00414 G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
00415 vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
00416 vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() );
00417 #ifdef G4EVERBOSE
00418 if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
00419 #endif
00420
00421
00422
00423
00424 transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm)
00425 +2.*deltaPInv*pAver;
00426
00427 transf[0][1] = -deltaPInv/thetaAver*
00428 ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) +
00429 sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
00430 OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) );
00431
00432 transf[0][2] = -sinpPre*deltaPInv/thetaAver*
00433 ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) +
00434 sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
00435 OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) );
00436
00437 transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
00438
00439 transf[0][4] = -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
00440
00441
00442 transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())
00443 *(1.+deltaPInv*pAver);
00444 #ifdef G4EVERBOSE
00445 if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0] << " " << -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z()
00446 << " " << deltaPInv<< " " << pAver << G4endl;
00447 #endif
00448
00449 transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
00450 sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
00451 OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
00452 (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
00453 ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
00454 OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
00455 TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
00456
00457 transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) +
00458 sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
00459 OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )*
00460 (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
00461 ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
00462 OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) -
00463 TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
00464 transf[1][2] = sinpPre*transf[1][2];
00465
00466 transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
00467
00468 transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
00469
00470
00471
00472 transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv
00473 *(1.+deltaPInv*pAver);
00474 #ifdef G4EVERBOSE
00475 if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv
00476 <<" "<<deltaPInv<<" "<<pAver<< G4endl;
00477 #endif
00478 transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) +
00479 sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) +
00480 OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
00481 (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) +
00482 ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
00483 OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
00484 TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
00485 transf[2][1] = sinpPostInv*transf[2][1];
00486
00487 transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) +
00488 sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) +
00489 OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )*
00490 (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) +
00491 ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
00492 OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) -
00493 TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
00494 transf[2][2] = sinpPostInv*sinpPre*transf[2][2];
00495
00496 transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv;
00497 #ifdef G4EVERBOSE
00498 if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl;
00499 #endif
00500
00501 transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv;
00502
00503
00504
00505 transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() )
00506 *(1.+deltaPInv*pAver);
00507 #ifdef G4EVERBOSE
00508 if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y()
00509 <<" "<<deltaPInv<<" "<<pAver<<G4endl;
00510 #endif
00511
00512 transf[3][1] = ( sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) +
00513 OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) +
00514 TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )*
00515 (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
00516
00517 transf[3][2] = ( sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) +
00518 OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) +
00519 TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )*
00520 (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
00521 #ifdef G4EVERBOSE
00522 if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<<
00523 OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<<
00524 TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<<
00525 vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl;
00526 #endif
00527
00528 transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() );
00529
00530 transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() );
00531
00532
00533 transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z())
00534 *(1.+deltaPInv*pAver);
00535
00536 transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
00537 OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
00538 TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
00539 (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
00540 #ifdef G4EVERBOSE
00541 if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl;
00542 #endif
00543
00544 transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) +
00545 OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
00546 TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
00547 (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
00548
00549 transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() );
00550
00551 transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z());
00552
00553
00554
00555 #ifdef G4EVERBOSE
00556 if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl;
00557 #endif
00558
00559
00560
00561
00562
00563
00564 }
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576 theTransfMat = transf;
00577 #ifdef G4EVERBOSE
00578 if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
00579 if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl
00580 << " transf matrix " << theTransfMat.T() << G4endl;
00581 #endif
00582
00583 fError = fError.similarity(theTransfMat).T();
00584
00585 #ifdef G4EVERBOSE
00586 if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl;
00587 #endif
00588
00589
00590
00591
00592
00593 PropagateErrorMSC( aTrack );
00594
00595 PropagateErrorIoni( aTrack );
00596
00597 return 0;
00598 }
00599
00600
00601
00602 G4int G4ErrorFreeTrajState::PropagateErrorMSC( const G4Track* aTrack )
00603 {
00604 G4ThreeVector vpPre = aTrack->GetMomentum()/GeV;
00605 G4double pPre = vpPre.mag();
00606 G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV);
00607 G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
00608
00609 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
00610 G4double effZ, effA;
00611 CalculateEffectiveZandA( mate, effZ, effA );
00612
00613 #ifdef G4EVERBOSE
00614 if( iverbose >= 4 ) G4cout << "material " << mate->GetName()
00615
00616 << " effZ:" << effZ << " effA:" << effA
00617 << " dens(g/mole):" << mate->GetDensity()/g*mole << " Radlen/cm:" << mate->GetRadlen()/cm << " nuclLen/cm" << mate->GetNuclearInterLength()/cm << G4endl;
00618 #endif
00619
00620 G4double RI = stepLengthCm / (mate->GetRadlen()/cm);
00621 #ifdef G4EVERBOSE
00622 if( iverbose >= 4 ) G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI=X/X0 " << RI << " stepLengthCm " << stepLengthCm << " radlen/cm " << (mate->GetRadlen()/cm) << " RI*1.e10:" << RI*1.e10 << G4endl;
00623 #endif
00624 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
00625 G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta );
00626 #ifdef G4EVERBOSE
00627 if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl;
00628 #endif
00629 G4double S1 = DD*stepLengthCm*stepLengthCm/3.;
00630 G4double S2 = DD;
00631 G4double S3 = DD*stepLengthCm/2.;
00632
00633 G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre;
00634 #ifdef G4EVERBOSE
00635 if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 " << S2 << " S3 " << S3 << " CLA " << CLA << G4endl;
00636 #endif
00637 fError[1][1] += S2;
00638 fError[1][4] -= S3;
00639 fError[2][2] += S2/CLA/CLA;
00640 fError[2][3] += S3/CLA;
00641 fError[3][3] += S1;
00642 fError[4][4] += S1;
00643
00644 #ifdef G4EVERBOSE
00645 if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl;
00646 #endif
00647
00648 return 0;
00649 }
00650
00651
00652
00653 void G4ErrorFreeTrajState::CalculateEffectiveZandA( const G4Material* mate, G4double& effZ, G4double& effA )
00654 {
00655 effZ = 0.;
00656 effA = 0.;
00657 G4int ii, nelem = mate->GetNumberOfElements();
00658 const G4double* fracVec = mate->GetFractionVector();
00659 for(ii=0; ii < nelem; ii++ ) {
00660 effZ += mate->GetElement( ii )->GetZ() * fracVec[ii];
00661 effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole;
00662 }
00663
00664 }
00665
00666
00667
00668 G4int G4ErrorFreeTrajState::PropagateErrorIoni( const G4Track* aTrack )
00669 {
00670 G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
00671 #ifdef G4EVERBOSE
00672 G4double DEDX2;
00673 if( stepLengthCm < 1.E-7 ) {
00674 DEDX2=0.;
00675 }
00676 #endif
00677
00678 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
00679 G4double effZ, effA;
00680 CalculateEffectiveZandA( mate, effZ, effA );
00681
00682 G4double Etot = aTrack->GetTotalEnergy()/GeV;
00683 G4double beta = aTrack->GetMomentum().mag()/GeV / Etot;
00684 G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV;
00685 G4double gamma = Etot / mass;
00686
00687
00688 G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) /
00689 (effA*beta*beta);
00690
00691 #ifdef G4EVERBOSE
00692 if( iverbose >= 2 ){
00693 G4cout << "G4EP:IONI: XI/keV " << XI << " beta " << beta << " gamma " << gamma << G4endl;
00694 G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl;
00695 }
00696 #endif
00697
00698 G4double eta = beta*gamma;
00699 G4double etasq = eta*eta;
00700 G4double eMass = 0.51099906/GeV;
00701 G4double massRatio = eMass / mass;
00702 G4double F1 = 2*eMass*etasq;
00703 G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio;
00704 G4double Emax = 1.E+6*F1/F2;
00705
00706
00707 G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12;
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718 if (XI/Emax<0.01) dedxSq *=XI/Emax*100 ;
00719
00720 #ifdef G4EVERBOSE
00721 if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << dedxSq << " emass/GeV: " << eMass << " Emax/keV: " << Emax
00722 <<" k=Xi/Emax="<< XI/Emax<< G4endl;
00723
00724 #endif
00725
00726 G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag();
00727 pPre6 = std::pow(pPre6, 6 );
00728
00729 fError[0][0] += Etot*Etot*dedxSq / pPre6;
00730 #ifdef G4EVERBOSE
00731 if( iverbose >= 2 ) G4cout << "G4:IONI Etot/GeV: " << Etot << " err_dedx^2/GeV^2: " << dedxSq << " p^6: " << pPre6 << G4endl;
00732 if( iverbose >= 2 ) G4cout << "G4EP:IONI: error2_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl;
00733 #endif
00734
00735 return 0;
00736 }