diff --git a/Tools/PROCTools/data/master_q431_AOD_digest.ref b/Tools/PROCTools/data/master_q431_AOD_digest.ref
index 409bddb5d77767f87168ed0676328515813e2329..997e935bd5782d51fb237592f6d83a4b0745d99b 100644
--- a/Tools/PROCTools/data/master_q431_AOD_digest.ref
+++ b/Tools/PROCTools/data/master_q431_AOD_digest.ref
@@ -11,10 +11,10 @@
330470 1183738949 368 444 9 1
330470 1183742489 152 127 2 1
330470 1183743040 285 326 5 0
- 330470 1183746343 492 495 12 1
+ 330470 1183746343 492 495 12 0
330470 1183746710 6 0 0 0
330470 1183751782 239 246 2 0
- 330470 1183752624 347 366 7 3
+ 330470 1183752624 347 366 7 2
330470 1183753006 357 398 11 3
330470 1183754806 470 423 14 0
330470 1183769295 342 336 7 2
diff --git a/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/TrkExSTEP_Propagator/STEP_Propagator.h b/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/TrkExSTEP_Propagator/STEP_Propagator.h
index 15545780fc84580c192b9ac6d98698834cfb7baf..34f894acdbe4b246fb149dc032e8b2807ee4fcfb 100755
--- a/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/TrkExSTEP_Propagator/STEP_Propagator.h
+++ b/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/TrkExSTEP_Propagator/STEP_Propagator.h
@@ -442,7 +442,7 @@ namespace Trk {
double& h,
double* P,
double* dDir,
- float* BG1,
+ double* BG1,
bool& firstStep,
double& distanceStepped) const;
@@ -456,7 +456,7 @@ namespace Trk {
getMagneticField(Cache& cache,
const Amg::Vector3D& position,
bool getGradients,
- float* BG) const;
+ double* BG) const;
/////////////////////////////////////////////////////////////////////////////////
diff --git a/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/src/STEP_Propagator.cxx b/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/src/STEP_Propagator.cxx
index e6bcfae97a26a8583d2e0add7aa34908217a5e9e..c7e7b2ed361339ad9561ac18bf2dbb400b747dc6 100755
--- a/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/src/STEP_Propagator.cxx
+++ b/Tracking/TrkExtrapolation/TrkExSTEP_Propagator/src/STEP_Propagator.cxx
@@ -69,7 +69,7 @@ Trk::STEP_Propagator::STEP_Propagator
m_maxPath(100000.), //Maximum propagation length in mm.
m_maxSteps(10000), //Maximum number of allowed steps (to avoid infinite loops).
m_layXmax(1.), // maximal layer thickness for multiple scattering calculations
- m_simulation(false), //flag for simulation mode
+ m_simulation(false), //flag for simulation mode
m_rndGenSvc("AtDSFMTGenSvc", n),
m_randomEngine(nullptr),
m_randomEngineName("FatrasRnd")
@@ -111,7 +111,7 @@ StatusCode Trk::STEP_Propagator::initialize()
// Read handle for AtlasFieldCacheCondObj
ATH_CHECK( m_fieldCacheCondObjInputKey.initialize() );
-
+
if (!m_materialEffects) { //override all material interactions
m_multipleScattering = false;
m_energyLoss = false;
@@ -137,8 +137,8 @@ StatusCode Trk::STEP_Propagator::initialize()
m_randomEngine = m_rndGenSvc->GetEngine(m_randomEngineName);
if (!m_randomEngine) {
ATH_MSG_WARNING( "Could not get random engine '" << m_randomEngineName << "', no smearing done." );
- }
- }
+ }
+ }
}
return StatusCode::SUCCESS;
@@ -160,7 +160,7 @@ Trk::STEP_Propagator::propagate (const Trk::NeutralParameters&,
const Trk::BoundaryCheck&,
bool) const
{
- ATH_MSG_WARNING( "[STEP_Propagator] STEP_Propagator does not handle neutral track parameters."
+ ATH_MSG_WARNING( "[STEP_Propagator] STEP_Propagator does not handle neutral track parameters."
<<"Use the StraightLinePropagator instead." );
return nullptr;
}
@@ -191,7 +191,7 @@ Trk::STEP_Propagator::propagate (const EventContext& ctx,
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
//Check for tracking volume (materialproperties)
cache.m_trackingVolume = tVol;
@@ -214,7 +214,7 @@ Trk::STEP_Propagator::propagate (const EventContext& ctx,
/////////////////////////////////////////////////////////////////////////////////
// Main function for track parameters and covariance matrix propagation
-// with search of closest surface (ST)
+// with search of closest surface (ST)
/////////////////////////////////////////////////////////////////////////////////
Trk::TrackParameters*
@@ -239,7 +239,7 @@ Trk::STEP_Propagator::propagate (const EventContext& ctx,
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
//Check for tracking volume (materialproperties)
cache.m_trackingVolume = tVol;
@@ -264,19 +264,19 @@ Trk::STEP_Propagator::propagate (const EventContext& ctx,
} else {
cache.m_propagateWithPathLimit = 0;
cache.m_pathLimit = -1.;
- path = 0.;
- }
+ path = 0.;
+ }
if ( particle==Trk::neutron || particle==Trk::photon || particle==Trk::pi0 || particle==Trk::k0 )
return propagateNeutral(trackParameters,targetSurfaces,propagationDirection,solutions,path,usePathLimit,returnCurv);
return propagateRungeKutta( cache,true, trackParameters, targetSurfaces, propagationDirection, magneticFieldProperties,
- particle, solutions, path, returnCurv);
+ particle, solutions, path, returnCurv);
}
/////////////////////////////////////////////////////////////////////////////////
// Main function for track parameters and covariance matrix propagation
-// with search of closest surface and time info (ST)
+// with search of closest surface and time info (ST)
/////////////////////////////////////////////////////////////////////////////////
Trk::TrackParameters*
@@ -302,7 +302,7 @@ Trk::STEP_Propagator::propagateT (const EventContext& ctx,
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
// cache particle mass
cache.m_particleMass = s_particleMasses.mass[particle]; //Get particle mass from ParticleHypothesis
@@ -329,16 +329,16 @@ Trk::STEP_Propagator::propagateT (const EventContext& ctx,
double dMat = pathLim.x0Max-pathLim.x0Collected;
double path = dMat>0 && cache.m_matPropOK && cache.m_material->x0()>0. ? dMat * cache.m_material->x0() : -1.;
- double dTim = timeLim.tMax - timeLim.time;
+ double dTim = timeLim.tMax - timeLim.time;
double beta = 1.;
if (dTim>0.) {
- double mom = trackParameters.momentum().mag();
+ double mom = trackParameters.momentum().mag();
beta = mom/std::sqrt(mom*mom+cache.m_particleMass*cache.m_particleMass);
}
double timMax = dTim>0 ? dTim*beta*Gaudi::Units::c_light : -1.;
if ( timMax>0. && timMax < path ) path = timMax;
- bool usePathLimit = (path > 0.);
+ bool usePathLimit = (path > 0.);
// resolve path limit input
if (path>0.) {
@@ -348,32 +348,32 @@ Trk::STEP_Propagator::propagateT (const EventContext& ctx,
} else {
cache.m_propagateWithPathLimit = 0;
cache.m_pathLimit = -1.;
- path = 0.;
- }
+ path = 0.;
+ }
Trk::TrackParameters* nextPar = nullptr;
if ( particle==Trk::neutron || particle==Trk::photon || particle==Trk::pi0 || particle==Trk::k0 ){
nextPar = propagateNeutral(trackParameters,targetSurfaces,propagationDirection,
solutions,path,usePathLimit,returnCurv);
- } else{
- nextPar = propagateRungeKutta( cache,true, trackParameters, targetSurfaces, propagationDirection,
+ } else{
+ nextPar = propagateRungeKutta( cache,true, trackParameters, targetSurfaces, propagationDirection,
magneticFieldProperties,
- particle, solutions, path,returnCurv);
+ particle, solutions, path,returnCurv);
}
// update material path
if (cache.m_matPropOK && cache.m_material->x0()>0. && path>0.) {
- pathLim.updateMat( path/cache.m_material->x0(),cache.m_material->averageZ(),0.);
+ pathLim.updateMat( path/cache.m_material->x0(),cache.m_material->averageZ(),0.);
}
// return value
timeLim.time +=cache.m_timeOfFlight;
- return nextPar;
+ return nextPar;
}
/////////////////////////////////////////////////////////////////////////////////
// Main function for track parameters and covariance matrix propagation
-// with search of closest surface and material collection (ST)
+// with search of closest surface and material collection (ST)
/////////////////////////////////////////////////////////////////////////////////
Trk::TrackParameters*
@@ -401,7 +401,7 @@ Trk::STEP_Propagator::propagateM (const EventContext&
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
//Check for tracking volume (materialproperties)
cache.m_trackingVolume = tVol;
@@ -430,15 +430,15 @@ Trk::STEP_Propagator::propagateM (const EventContext&
} else {
cache.m_propagateWithPathLimit = 0;
cache.m_pathLimit = -1.;
- path = 0.;
- }
+ path = 0.;
+ }
if ( particle==Trk::neutron || particle==Trk::photon || particle==Trk::pi0 || particle==Trk::k0 ){
return propagateNeutral(trackParameters,targetSurfaces,propagationDirection,solutions,path,
usePathLimit,returnCurv);
}
- return propagateRungeKutta( cache,true, trackParameters, targetSurfaces,
+ return propagateRungeKutta( cache,true, trackParameters, targetSurfaces,
propagationDirection, magneticFieldProperties,
- particle, solutions, path,returnCurv);
+ particle, solutions, path,returnCurv);
}
/////////////////////////////////////////////////////////////////////////////////
@@ -468,7 +468,7 @@ Trk::STEP_Propagator::propagate (const EventContext& ctx,
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
//Check for tracking volume (materialproperties)
cache.m_trackingVolume = tVol;
@@ -485,8 +485,8 @@ Trk::STEP_Propagator::propagate (const EventContext& ctx,
cache.m_matupd_lastpath = 0.;
cache.m_matdump_lastpath = 0.;
- Trk::TrackParameters* parameters = propagateRungeKutta(cache, true, trackParameters, targetSurface,
- propagationDirection,magneticFieldProperties,
+ Trk::TrackParameters* parameters = propagateRungeKutta(cache, true, trackParameters, targetSurface,
+ propagationDirection,magneticFieldProperties,
particle, boundaryCheck, Jacobian, returnCurv);
if (parameters) {
@@ -525,7 +525,7 @@ Trk::STEP_Propagator::propagateParameters (const EventContext& c
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
//Check for tracking volume (materialproperties)
@@ -538,9 +538,9 @@ Trk::STEP_Propagator::propagateParameters (const EventContext& c
cache.m_matstates = nullptr;
cache.m_hitVector = nullptr;
- return propagateRungeKutta( cache,false, trackParameters,
+ return propagateRungeKutta( cache,false, trackParameters,
targetSurface, propagationDirection,
- magneticFieldProperties, particle, boundaryCheck,
+ magneticFieldProperties, particle, boundaryCheck,
Jacobian, returnCurv);
}
@@ -563,7 +563,7 @@ Trk::STEP_Propagator::propagateParameters (const EventContext& c
{
// ATH_MSG_WARNING( "[STEP_Propagator] enter 7");
-
+
double Jacobian[25];
Cache cache{};
@@ -572,7 +572,7 @@ Trk::STEP_Propagator::propagateParameters (const EventContext& c
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
//Check for tracking volume (materialproperties)
cache.m_trackingVolume = tVol;
@@ -585,11 +585,11 @@ Trk::STEP_Propagator::propagateParameters (const EventContext& c
cache.m_extrapolationCache = nullptr;
cache.m_hitVector = nullptr;
- Trk::TrackParameters* parameters = propagateRungeKutta( cache,true, trackParameters, targetSurface,
+ Trk::TrackParameters* parameters = propagateRungeKutta( cache,true, trackParameters, targetSurface,
propagationDirection,magneticFieldProperties,
particle, boundaryCheck, Jacobian, returnCurv);
- if (parameters) {
+ if (parameters) {
Jacobian[24]=Jacobian[20]; Jacobian[23]=0.; Jacobian[22]=0.; Jacobian[21]=0.; Jacobian[20]=0.;
jacobian = new Trk::TransportJacobian(Jacobian);
} else {
@@ -619,7 +619,7 @@ Trk::STEP_Propagator::intersect (const EventContext& ctx,
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
cache.m_particle = particle; //Store for later use
@@ -636,7 +636,7 @@ Trk::STEP_Propagator::intersect (const EventContext& ctx,
cache.m_hitVector = nullptr;
// Bfield mode
- mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
+ mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
//Check inputvalues
if (m_tolerance <= 0.) return nullptr;
@@ -656,51 +656,51 @@ Trk::STEP_Propagator::intersect (const EventContext& ctx,
if (!Trk::RungeKuttaUtils::transformLocalToGlobal( false, trackParameters, cache.m_P)) return nullptr;
double path = 0.;
- const Amg::Transform3D& T = targetSurface.transform();
- Trk::Surface::SurfaceType ty = targetSurface.type();
+ const Amg::Transform3D& T = targetSurface.transform();
+ Trk::Surface::SurfaceType ty = targetSurface.type();
- if (ty == Trk::Surface::Plane || ty == Trk::Surface::Disc ) {
+ if (ty == Trk::Surface::Plane || ty == Trk::Surface::Disc ) {
double s[4];
- double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+ double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
- if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
- else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
if (!propagateWithJacobian(cache, false, ty, s, cache.m_P, path)) return nullptr;
}
- else if (ty == Trk::Surface::Line ) {
+ else if (ty == Trk::Surface::Line ) {
- double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
if (!propagateWithJacobian( cache,false, ty, s, cache.m_P, path)) return nullptr;
}
- else if (ty == Trk::Surface::Cylinder ) {
+ else if (ty == Trk::Surface::Cylinder ) {
- const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(&targetSurface);
- double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),Trk::alongMomentum,0.};
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(&targetSurface);
+ double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),Trk::alongMomentum,0.};
if (!propagateWithJacobian(cache, false, ty, s, cache.m_P, path)) return nullptr;
}
- else if (ty == Trk::Surface::Cone ) {
+ else if (ty == Trk::Surface::Cone ) {
- double k = static_cast<const Trk::ConeSurface*>(&targetSurface)->bounds().tanAlpha(); k = k*k+1.;
- double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,Trk::alongMomentum,0.};
+ double k = static_cast<const Trk::ConeSurface*>(&targetSurface)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,Trk::alongMomentum,0.};
if (!propagateWithJacobian(cache, false, ty, s, cache.m_P, path)) return nullptr;
}
- else if (ty == Trk::Surface::Perigee ) {
+ else if (ty == Trk::Surface::Perigee ) {
- double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
if (!propagateWithJacobian( cache,false, ty, s, cache.m_P, path)) return nullptr;
}
else { // presumably curvilinear
double s[4];
- double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+ double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
- if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
- else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
if (!propagateWithJacobian( cache,false, ty, s, cache.m_P, path)) return nullptr;
}
@@ -715,11 +715,11 @@ Trk::STEP_Propagator::intersect (const EventContext& ctx,
const Trk::TrackSurfaceIntersection*
Trk::STEP_Propagator::intersectSurface(const EventContext& ctx,
const Trk::Surface& surface,
- const Trk::TrackSurfaceIntersection*
+ const Trk::TrackSurfaceIntersection*
trackIntersection,
const double qOverP,
const Trk::MagneticFieldProperties& mft,
- ParticleHypothesis particle) const
+ ParticleHypothesis particle) const
{
const Amg::Vector3D& origin = trackIntersection->position();
@@ -743,10 +743,10 @@ Trk::STEP_Propagator::intersectSurface(const EventContext& ctx,
Trk::IntersectionSolutionIter output_iter = solution->begin();
if(*output_iter) {
Trk::TrackSurfaceIntersection* result = new Trk::TrackSurfaceIntersection(*(*output_iter));
- delete solution;
- return result;
- }
- delete solution;
+ delete solution;
+ return result;
+ }
+ delete solution;
return nullptr;
}
@@ -771,7 +771,7 @@ Trk::STEP_Propagator::globalPositions ( const EventContext& ctx,
getFieldCacheObject(cache, ctx);
cache.m_detailedElossFlag=m_detailedEloss;
- clearCache(cache);
+ clearCache(cache);
cache.m_particle = particle; //Store for later use
@@ -787,7 +787,7 @@ Trk::STEP_Propagator::globalPositions ( const EventContext& ctx,
cache.m_matPropOK = false;
}
- mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
+ mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
//Check inputvalues
if (m_tolerance <= 0.) return;
@@ -798,7 +798,7 @@ Trk::STEP_Propagator::globalPositions ( const EventContext& ctx,
double maxPath = m_maxPath; // Max path allowed
double dDir[3] = { 0., 0., 0.}; // Start directions derivs. Zero in case of no RK steps
double distanceStepped = 0.;
- float BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
+ double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
// dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz
bool firstStep = true; // Poll B1, else recycle B4
double path = 0.; // path of the trajectory
@@ -874,7 +874,7 @@ Trk::STEP_Propagator::globalPositions ( const EventContext& ctx,
__attribute__ ((optimize(2)))
#endif
Trk::TrackParameters*
-Trk::STEP_Propagator::propagateRungeKutta (Cache& cache,
+Trk::STEP_Propagator::propagateRungeKutta (Cache& cache,
bool errorPropagation,
const Trk::TrackParameters& inputTrackParameters,
const Trk::Surface& targetSurface,
@@ -892,7 +892,7 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
const Trk::TrackParameters* trackParameters = nullptr;
// Bfield mode
- mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
+ mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
//Check inputvalues
if (m_tolerance <= 0.) return nullptr;
@@ -935,70 +935,70 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
}
double path = 0.;
- const Amg::Transform3D& T = targetSurface.transform();
- Trk::Surface::SurfaceType ty = targetSurface.type();
+ const Amg::Transform3D& T = targetSurface.transform();
+ Trk::Surface::SurfaceType ty = targetSurface.type();
- if (ty == Trk::Surface::Plane || ty == Trk::Surface::Disc ) {
+ if (ty == Trk::Surface::Plane || ty == Trk::Surface::Disc ) {
double s[4];
- double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+ double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
- if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
- else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
if (!propagateWithJacobian(cache,errorPropagation, ty, s, cache.m_P, path)) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
return nullptr;
}
}
- else if (ty == Trk::Surface::Line ) {
+ else if (ty == Trk::Surface::Line ) {
- double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
if (!propagateWithJacobian( cache,errorPropagation, ty, s, cache.m_P, path)) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
return nullptr;
}
}
- else if (ty == Trk::Surface::Cylinder ) {
+ else if (ty == Trk::Surface::Cylinder ) {
- const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(&targetSurface);
- double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),(double)propagationDirection,0.};
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(&targetSurface);
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),(double)propagationDirection,0.};
if (!propagateWithJacobian( cache,errorPropagation, ty, s, cache.m_P, path)) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
return nullptr;
- }
+ }
}
- else if (ty == Trk::Surface::Cone ) {
+ else if (ty == Trk::Surface::Cone ) {
- double k = static_cast<const Trk::ConeSurface*>(&targetSurface)->bounds().tanAlpha(); k = k*k+1.;
- double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,(double)propagationDirection,0.};
+ double k = static_cast<const Trk::ConeSurface*>(&targetSurface)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,(double)propagationDirection,0.};
if (!propagateWithJacobian(cache, errorPropagation, ty, s, cache.m_P, path)) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
return nullptr;
- }
+ }
}
- else if (ty == Trk::Surface::Perigee ) {
+ else if (ty == Trk::Surface::Perigee ) {
- double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
if (!propagateWithJacobian(cache, errorPropagation, ty, s, cache.m_P, path)) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
return nullptr;
- }
+ }
}
else { // presumably curvilinear
double s[4];
- double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+ double d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
- if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
- else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
if (!propagateWithJacobian(cache, errorPropagation, ty, s, cache.m_P, path)) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
return nullptr;
- }
+ }
}
if (propagationDirection * path < 0.) {
@@ -1007,7 +1007,7 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
}
double localp[5];
- // output in curvilinear parameters
+ // output in curvilinear parameters
if (returnCurv || ty==Trk::Surface::Cone) {
Trk::RungeKuttaUtils::transformGlobalToLocal(cache.m_P,localp);
@@ -1015,7 +1015,7 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
if ( boundaryCheck && !targetSurface.isOnSurface(gp) ) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
- return nullptr;
+ return nullptr;
}
if ( !errorPropagation || !trackParameters->covariance() ) {
@@ -1028,14 +1028,14 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
AmgSymMatrix(5)* measurementCovariance = Trk::RungeKuttaUtils::newCovarianceMatrix(
Jacobian, *trackParameters->covariance());
- if (cache.m_matPropOK && (m_multipleScattering || m_straggling) && fabs(path)>0. )
+ if (cache.m_matPropOK && (m_multipleScattering || m_straggling) && fabs(path)>0. )
covarianceContribution( cache,trackParameters, path, fabs( 1./cache.m_P[6]), measurementCovariance);
if (trackParameters != &inputTrackParameters) delete trackParameters;
- return new Trk::CurvilinearParameters(gp,localp[2],localp[3],localp[4],measurementCovariance);
+ return new Trk::CurvilinearParameters(gp,localp[2],localp[3],localp[4],measurementCovariance);
}
- // Common transformation for all surfaces
+ // Common transformation for all surfaces
Trk::RungeKuttaUtils::transformGlobalToLocal(&targetSurface,errorPropagation,cache.m_P,localp,Jacobian);
if (boundaryCheck) {
@@ -1047,11 +1047,11 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
}
Trk::TrackParameters* onTargetSurf = targetSurface.createTrackParameters(localp[0],localp[1],localp[2],
- localp[3],localp[4],nullptr);
+ localp[3],localp[4],nullptr);
if ( !errorPropagation || !trackParameters->covariance() ) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
- return onTargetSurf;
+ return onTargetSurf;
}
//Errormatrix is included. Use Jacobian to calculate new covariance
@@ -1059,12 +1059,12 @@ Trk::STEP_Propagator::propagateRungeKutta (Cache& c
Jacobian, *trackParameters->covariance());
//Calculate multiple scattering and straggling covariance contribution.
- if (cache.m_matPropOK && (m_multipleScattering || m_straggling) && fabs(path)>0. )
+ if (cache.m_matPropOK && (m_multipleScattering || m_straggling) && fabs(path)>0. )
covarianceContribution( cache,trackParameters, path, onTargetSurf, measurementCovariance);
delete onTargetSurf; // the covariance matrix can be just added instead of recreating ?
if (trackParameters != &inputTrackParameters) delete trackParameters;
- return targetSurface.createTrackParameters(localp[0],localp[1],localp[2],localp[3],localp[4],measurementCovariance);
+ return targetSurface.createTrackParameters(localp[0],localp[1],localp[2],localp[3],localp[4],measurementCovariance);
}
/////////////////////////////////////////////////////////////////////////////////
@@ -1085,14 +1085,14 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
bool returnCurv) const
{
//Store for later use
- cache.m_particle = particle;
+ cache.m_particle = particle;
cache.m_charge = inputTrackParameters.charge();
cache.m_inputThetaVariance = 0.;
const Trk::TrackParameters* trackParameters = nullptr;
// Bfield mode
- mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
+ mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
//Check inputvalues
if (m_tolerance <= 0.) return nullptr;
@@ -1152,9 +1152,9 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
// Common transformation for all surfaces (angles and momentum)
double localp[5];
double Jacobian[21];
- while ( validStep ) {
+ while ( validStep ) {
// propagation to next surface
- validStep = propagateWithJacobian( cache,errorPropagation, targetSurfaces, cache.m_P,
+ validStep = propagateWithJacobian( cache,errorPropagation, targetSurfaces, cache.m_P,
propagationDirection, solutions, path, totalPath);
if (!validStep) {
if (trackParameters != &inputTrackParameters) delete trackParameters;
@@ -1165,14 +1165,14 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
return nullptr;
}
totalPath += path; cache.m_timeOfFlight += cache.m_timeStep;
- if (cache.m_propagateWithPathLimit>1 || cache.m_binMat ) {
+ if (cache.m_propagateWithPathLimit>1 || cache.m_binMat ) {
// make sure that for sliding surfaces the result does not get distorted
// return curvilinear parameters
Trk::CurvilinearParameters* cPar = nullptr;
Trk::RungeKuttaUtils::transformGlobalToLocal(cache.m_P, localp);
- if (!errorPropagation) {
+ if (!errorPropagation) {
cPar = new Trk::CurvilinearParameters(Amg::Vector3D(cache.m_P[0],cache.m_P[1],cache.m_P[2]),
- localp[2],localp[3],localp[4]);
+ localp[2],localp[3],localp[4]);
} else {
double useless[2];
Trk::RungeKuttaUtils::transformGlobalToCurvilinear( true, cache.m_P, useless, Jacobian);
@@ -1189,8 +1189,8 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
}
// material collection : first iteration, bin material averaged
// collect material
- if ( cache.m_binMat && (cache.m_matstates ||
- (errorPropagation && cache.m_extrapolationCache)) &&
+ if ( cache.m_binMat && (cache.m_matstates ||
+ (errorPropagation && cache.m_extrapolationCache)) &&
fabs(totalPath-cache.m_matdump_lastpath)>1.) {
dumpMaterialEffects( cache,cPar, totalPath);
}
@@ -1199,23 +1199,23 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
if (cache.m_propagateWithPathLimit>0) cache.m_pathLimit -= path;
// boundary check
// take into account that there may be many identical surfaces with different boundaries
- Amg::Vector3D gp(cache.m_P[0],cache.m_P[1],cache.m_P[2]);
+ Amg::Vector3D gp(cache.m_P[0],cache.m_P[1],cache.m_P[2]);
bool solution = false;
std::vector<unsigned int> valid_solutions;
valid_solutions.reserve(solutions.size());
std::vector<unsigned int>::iterator iSol= solutions.begin();
- while ( iSol != solutions.end() ) {
+ while ( iSol != solutions.end() ) {
if ( targetSurfaces[*iSol].first->isOnSurface(gp,targetSurfaces[*iSol].second ,0.001,0.001) ) {
if (!solution) {
Trk::RungeKuttaUtils::transformGlobalToLocal(cache.m_P, localp);
if (returnCurv || targetSurfaces[*iSol].first->type()==Trk::Surface::Cone) {
- Trk::RungeKuttaUtils::transformGlobalToCurvilinear(errorPropagation,cache.m_P,localp,Jacobian);
+ Trk::RungeKuttaUtils::transformGlobalToCurvilinear(errorPropagation,cache.m_P,localp,Jacobian);
} else Trk::RungeKuttaUtils::transformGlobalToLocal(targetSurfaces[*iSol].first,errorPropagation,cache.m_P,localp,Jacobian);
solution = true;
}
valid_solutions.push_back( *iSol );
- }
+ }
iSol++;
}
solutions = valid_solutions;
@@ -1233,7 +1233,7 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
smear(cache,localp[2],localp[3],trackParameters,radDist);
}
- Trk::TrackParameters* onTargetSurf = (returnCurv || targetSurfaces[solutions[0]].first->type()==Trk::Surface::Cone) ?
+ Trk::TrackParameters* onTargetSurf = (returnCurv || targetSurfaces[solutions[0]].first->type()==Trk::Surface::Cone) ?
nullptr : targetSurfaces[solutions[0]].first->createTrackParameters(localp[0],localp[1],localp[2],localp[3],localp[4],nullptr);
if (!errorPropagation) {
@@ -1266,7 +1266,7 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
delete onTargetSurf; // the covariance matrix can be just added instead of recreating ?
return targetSurfaces[solutions[0]].first->createTrackParameters(localp[0],localp[1],localp[2],localp[3],localp[4],
- measurementCovariance);
+ measurementCovariance);
}
@@ -1276,11 +1276,11 @@ Trk::STEP_Propagator::propagateRungeKutta ( Cache&
/////////////////////////////////////////////////////////////////////////////////
bool
-Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
+Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
bool errorPropagation,
Trk::Surface::SurfaceType surfaceType,
double* targetSurface,
- double* P,
+ double* P,
double& path) const
{
double maxPath = m_maxPath; // Max path allowed
@@ -1291,7 +1291,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
double previousDistance = 0.;
double distanceStepped = 0.;
bool distanceEstimationSuccessful = false; // Was the distance estimation successful?
- float BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
+ double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
// dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz at first point
bool firstStep = true; // Poll BG1, else recycle BG4
double absolutePath = 0.; // Absolute path to register oscillating behaviour
@@ -1322,15 +1322,15 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if(fabs(distanceStepped)>0.001) {
cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL*log(fabs(distanceStepped));// the non-linear term
}
- // update straggling covariance
+ // update straggling covariance
if (errorPropagation && m_straggling) {
double sigTot2 = cache.m_sigmaIoni*cache.m_sigmaIoni + cache.m_sigmaRad*cache.m_sigmaRad;
// /(beta*beta*p*p*p*p) transforms Var(E) to Var(q/p)
mom = fabs(1./P[6]); beta = mom/std::sqrt(mom*mom+cache.m_particleMass*cache.m_particleMass);
double bp2 = beta*mom*mom;
- cache.m_stragglingVariance += sigTot2/(bp2*bp2)*distanceStepped*distanceStepped;
+ cache.m_stragglingVariance += sigTot2/(bp2*bp2)*distanceStepped*distanceStepped;
}
- if (cache.m_matstates || errorPropagation)
+ if (cache.m_matstates || errorPropagation)
cache.m_combinedEloss.update(cache.m_delIoni*distanceStepped,cache.m_sigmaIoni*fabs(distanceStepped),
cache.m_delRad *distanceStepped,cache.m_sigmaRad *fabs(distanceStepped),m_MPV);
//Calculate new distance to target
@@ -1352,7 +1352,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if (steps++ > m_maxSteps) return false; //Too many steps, something is wrong
}
- if (cache.m_material && cache.m_material->x0()!=0.) cache.m_combinedThickness += fabs(path)/cache.m_material->x0();
+ if (cache.m_material && cache.m_material->x0()!=0.) cache.m_combinedThickness += fabs(path)/cache.m_material->x0();
//Use Taylor expansions to step the remaining distance (typically microns).
path = path + distanceToTarget;
@@ -1385,7 +1385,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
/////////////////////////////////////////////////////////////////////////////////
bool
-Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
+Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
bool errorPropagation,
std::vector<DestSurf >& sfs,
double* P,
@@ -1401,7 +1401,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
//int targetPassed = 0; // how many times have we passed the target?
double previousDistance = 0.;
double distanceStepped = 0.;
- float BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
+ double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
// dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz at first point
bool firstStep = true; // Poll BG1, else recycle BG4
int steps = 0;
@@ -1416,8 +1416,8 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
// limit number of recovery attempts
int restartLimit =10;
- Amg::Vector3D position(P[0],P[1],P[2]);
- Amg::Vector3D direction0(P[3],P[4],P[5]);
+ Amg::Vector3D position(P[0],P[1],P[2]);
+ Amg::Vector3D direction0(P[3],P[4],P[5]);
// binned material ?
cache.m_binMat = nullptr;
@@ -1429,7 +1429,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
// closest distance estimate
// maintain count of oscilations and previous distance for each surface;
// skip initial trivial solutions (input parameters at surface) - should be treated before call to the propagator !
- double tol = 0.001;
+ double tol = 0.001;
solutions.clear();
double distanceToTarget = propDir*maxPath;
cache.m_currentDist.resize(sfs.size()); // keep size through the call
@@ -1438,7 +1438,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
int nextSfCand=nextSf;
std::vector<DestSurf >::iterator sIter = sfs.begin();
std::vector<DestSurf >::iterator sBeg = sfs.begin();
- unsigned int numSf=0;
+ unsigned int numSf=0;
unsigned int iCurr=0; // index for m_currentDist
int startSf = -99;
for (; sIter!=sfs.end(); sIter++) {
@@ -1448,18 +1448,18 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if (distSol.numberOfSolutions()>0 ) {
distEst = distSol.first();
dist1Est = distSol.first();
- if ( distSol.numberOfSolutions()>1 && ( fabs(distEst) < tol ||
- (propDir*distEst<-tol && propDir*distSol.second()>tol)) )
+ if ( distSol.numberOfSolutions()>1 && ( fabs(distEst) < tol ||
+ (propDir*distEst<-tol && propDir*distSol.second()>tol)) )
distEst = distSol.second();
}
// select input surfaces;
// do not accept trivial solutions (at the surface)
- // but include them into further distance estimate (aca return to the same surface)
+ // but include them into further distance estimate (aca return to the same surface)
if ( distEst*propDir > -tol ) {
if (distSol.currentDistance()>500.) {
cache.m_currentDist[iCurr]=std::pair<int,std::pair<double,double> >
(0,std::pair<double,double>(distEst,distSol.currentDistance(true)));
- }else{
+ }else{
cache.m_currentDist[iCurr]=std::pair<int,std::pair<double,double> >
(1,std::pair<double,double>(distEst,distSol.currentDistance(true)));
}
@@ -1475,9 +1475,9 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
}
if(fabs(dist1Est)<tol) startSf = (int) iCurr;
iCurr++;
- }
+ }
- if (distanceToTarget == maxPath || numSf == 0 ) return false;
+ if (distanceToTarget == maxPath || numSf == 0 ) return false;
// these do not change
std::vector< std::pair<int,std::pair<double,double> > >::iterator vsBeg = cache.m_currentDist.begin();
@@ -1497,7 +1497,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
binIDMat = cache.m_binMat->material(position);
std::pair<size_t,float> dist2next = lbu->distanceToNext(position,propDir*direction0);
if (dist2next.first < lbu->bins() && fabs(dist2next.second)>1. && fabs(dist2next.second)< fabs(h) ){
- h = dist2next.second*propDir;
+ h = dist2next.second*propDir;
}
if (binIDMat) cache.m_material = binIDMat->first;
}
@@ -1511,22 +1511,22 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
// bremstrahlung : sample if activated
if (cache.m_brem) {
- mom = fabs(1./P[6]);
+ mom = fabs(1./P[6]);
sampleBrem(cache,mom);
}
- while ( numSf > 0 && ( fabs( distanceToTarget) > distanceTolerance ||
+ while ( numSf > 0 && ( fabs( distanceToTarget) > distanceTolerance ||
fabs(path+distanceStepped)<tol) ) { // Step until within tolerance
//Do the step. Stop the propagation if the energy goes below m_momentumCutOff
if (!rungeKuttaStep(cache, errorPropagation, h, P, dDir, BG1, firstStep, distanceStepped)) {
// emit brem photon before stopped ?
if (cache.m_brem) {
if ( m_momentumCutOff < cache.m_bremEmitThreshold && m_simMatUpdator ) {
- Amg::Vector3D position(P[0],P[1],P[2]);
- Amg::Vector3D direction(P[3],P[4],P[5]);
- m_simMatUpdator->recordBremPhoton(cache.m_timeIn+cache.m_timeOfFlight+cache.m_timeStep,
+ Amg::Vector3D position(P[0],P[1],P[2]);
+ Amg::Vector3D direction(P[3],P[4],P[5]);
+ m_simMatUpdator->recordBremPhoton(cache.m_timeIn+cache.m_timeOfFlight+cache.m_timeStep,
m_momentumCutOff, cache.m_bremMom, position, direction, cache.m_particle);
- // the recoil can be skipped here
+ // the recoil can be skipped here
for (int i=0; i<3; i++) P[3+i] = direction[i];
// end recoil ( momentum not adjusted here ! continuous energy loss maintained for the moment)
}
@@ -1540,13 +1540,13 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if(fabs(distanceStepped)>0.001) {
cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL*log(fabs(distanceStepped));
}
- // update straggling covariance
+ // update straggling covariance
if (errorPropagation && m_straggling) {
// 15% of the Radition moves the MOP value thus only 85% is accounted for by the Mean-MOP shift
double sigTot2 = cache.m_sigmaIoni*cache.m_sigmaIoni + cache.m_sigmaRad*cache.m_sigmaRad;
// /(beta*beta*p*p*p*p) transforms Var(E) to Var(q/p)
double bp2 = beta*mom*mom;
- cache.m_stragglingVariance += sigTot2/(bp2*bp2)*distanceStepped*distanceStepped;
+ cache.m_stragglingVariance += sigTot2/(bp2*bp2)*distanceStepped*distanceStepped;
}
if (cache.m_matstates||errorPropagation){
cache.m_combinedEloss.update(cache.m_delIoni*distanceStepped,
@@ -1554,7 +1554,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
cache.m_delRad *distanceStepped,cache.m_sigmaRad *fabs(distanceStepped),m_MPV);
}
if (cache.m_material && cache.m_material->x0()!=0.) {
- cache.m_combinedThickness += propDir*distanceStepped/cache.m_material->x0();
+ cache.m_combinedThickness += propDir*distanceStepped/cache.m_material->x0();
}
return false;
@@ -1567,43 +1567,43 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
cache.m_timeStep += distanceStepped/beta/Gaudi::Units::c_light;
if(fabs(distanceStepped)>0.001) cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL*log(fabs(distanceStepped));
- // update straggling covariance
+ // update straggling covariance
if (errorPropagation && m_straggling) {
// 15% of the Radition moves the MOP value thus only 85% is accounted for by the Mean-MOP shift
double sigTot2 = cache.m_sigmaIoni*cache.m_sigmaIoni + cache.m_sigmaRad*cache.m_sigmaRad;
// /(beta*beta*p*p*p*p) transforms Var(E) to Var(q/p)
double bp2 = beta*mom*mom;
- cache.m_stragglingVariance += sigTot2/(bp2*bp2)*distanceStepped*distanceStepped;
+ cache.m_stragglingVariance += sigTot2/(bp2*bp2)*distanceStepped*distanceStepped;
}
if (cache.m_matstates||errorPropagation){
cache.m_combinedEloss.update(cache.m_delIoni*distanceStepped,cache.m_sigmaIoni*fabs(distanceStepped),
cache.m_delRad *distanceStepped,cache.m_sigmaRad *fabs(distanceStepped),m_MPV);
}
if (cache.m_material && cache.m_material->x0()!=0.) {
- cache.m_combinedThickness += propDir*distanceStepped/cache.m_material->x0();
+ cache.m_combinedThickness += propDir*distanceStepped/cache.m_material->x0();
}
if (absPath > maxPath) return false;
// path limit implemented
- if (cache.m_propagateWithPathLimit>0 && cache.m_pathLimit<= path) { cache.m_propagateWithPathLimit++; return true; }
+ if (cache.m_propagateWithPathLimit>0 && cache.m_pathLimit<= path) { cache.m_propagateWithPathLimit++; return true; }
- bool restart = false;
+ bool restart = false;
// in case of problems, make shorter steps
if ( propDir*path < -tol || absPath-fabs(path)>10.) {
helpSoft = fabs(path)/absPath > 0.5 ? fabs(path)/absPath : 0.5;
- }
+ }
- Amg::Vector3D position(P[0],P[1],P[2]);
- Amg::Vector3D direction(P[3],P[4],P[5]);
+ Amg::Vector3D position(P[0],P[1],P[2]);
+ Amg::Vector3D direction(P[3],P[4],P[5]);
//Adapt step size to the material binning : change of bin layer triggers dump of material effects
- float distanceToNextBin = h; // default
+ float distanceToNextBin = h; // default
if (cache.m_binMat) {
const Trk::BinUtility* lbu = cache.m_binMat->layerBinUtility(position);
if (lbu) {
size_t layerBin = cache.m_binMat->layerBin(position);
- const Trk::IdentifiedMaterial* iMat = cache.m_binMat->material(position);
+ const Trk::IdentifiedMaterial* iMat = cache.m_binMat->material(position);
std::pair<size_t,float> dist2next = lbu->distanceToNext(position,propDir*direction);
distanceToNextBin = dist2next.second;
if (layerBin != cache.m_currentLayerBin ) { // step into new bin
@@ -1613,8 +1613,8 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
//std::cout <<" STEP overshoots bin boundary by:"<< stepOver<<" :w.r.t. bin:" << dist2previous.first<< std::endl;
double localp[5];
Trk::RungeKuttaUtils::transformGlobalToLocal(P, localp);
- const Trk::CurvilinearParameters* cPar =
- new Trk::CurvilinearParameters(Amg::Vector3D(P[0],P[1],P[2]),localp[2],localp[3],localp[4]);
+ const Trk::CurvilinearParameters* cPar =
+ new Trk::CurvilinearParameters(Amg::Vector3D(P[0],P[1],P[2]),localp[2],localp[3],localp[4]);
if (cache.m_identifiedParameters) {
if (binIDMat && binIDMat->second>0 && !iMat ) { // exit from active layer
cache.m_identifiedParameters->push_back(std::pair<const Trk::TrackParameters*,int> (cPar->clone(),-binIDMat->second));
@@ -1623,7 +1623,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
} else if (iMat && iMat->second>0) { // entry active layer
cache.m_identifiedParameters->push_back(std::pair<const Trk::TrackParameters*,int> (cPar->clone(),iMat->second));
}
- }
+ }
if (cache.m_hitVector) {
double hitTiming = cache.m_timeIn+cache.m_timeOfFlight+cache.m_timeStep;
if (binIDMat && binIDMat->second>0 && !iMat ) { // exit from active layer
@@ -1633,21 +1633,21 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
} else if (iMat && iMat->second>0) { // entry active layer
cache.m_hitVector->push_back(Trk::HitInfo(cPar->clone(), hitTiming, iMat->second,0.));
}
- }
+ }
delete cPar;
cache.m_currentLayerBin = layerBin;
binIDMat = iMat;
- if (binIDMat) {
- // change of material triggers update of the cache
- // @TODO Coverity complains about a possible NULL pointer dereferencing here
+ if (binIDMat) {
+ // change of material triggers update of the cache
+ // @TODO Coverity complains about a possible NULL pointer dereferencing here
// because the code above does not explicitly forbid m_material to be NULL and m_material is used
// unchecked inside updateMaterialEffects.
// Can m_material be NULL at this point ?
if (cache.m_material) {
updateMaterialEffects(cache,mom, sin(direction.theta()), sumPath+path-stepOver);
}
- cache.m_material = binIDMat->first;
+ cache.m_material = binIDMat->first;
}
// recalculate distance to next bin
if (distanceToNextBin<h) {
@@ -1658,16 +1658,16 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
} else if ( dist2next.first < lbu->bins() && fabs(distanceToNextBin) < 0.01 && h>0.01 ) { // tolerance 10 microns ?
double localp[5];
Trk::RungeKuttaUtils::transformGlobalToLocal(P, localp);
- const Trk::CurvilinearParameters* cPar =
- new Trk::CurvilinearParameters(Amg::Vector3D(P[0],P[1],P[2]),localp[2],localp[3],localp[4]);
+ const Trk::CurvilinearParameters* cPar =
+ new Trk::CurvilinearParameters(Amg::Vector3D(P[0],P[1],P[2]),localp[2],localp[3],localp[4]);
const Trk::IdentifiedMaterial* nextMat = binIDMat;
// need to know what comes next
Amg::Vector3D probe = position + (distanceToNextBin+0.01)*propDir*direction.normalized();
- nextMat = cache.m_binMat->material(probe);
+ nextMat = cache.m_binMat->material(probe);
//if (m_identifiedParameters || m_hitVector) {
- // std::cout <<"dump of identified parameters? step before:"<<distanceToNextBin<<": current:"<<
+ // std::cout <<"dump of identified parameters? step before:"<<distanceToNextBin<<": current:"<<
// m_currentLayerBin<<","<<dist2next.first<< ":"<<position.perp()<<","<<position.z()<<std::endl;
// if (binIDMat) std::cout <<"layer identification current:"<<binIDMat->second <<std::endl;
// if (nextMat) std::cout <<"layer identification next:"<<nextMat->second <<std::endl;
@@ -1675,7 +1675,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if (cache.m_identifiedParameters ) {
if (binIDMat && binIDMat->second>0 && !nextMat ) { // exit from active layer
cache.m_identifiedParameters->push_back(std::pair<const Trk::TrackParameters*,int> (cPar->clone(),-binIDMat->second));
- } else if (binIDMat && binIDMat->second>0 && (nextMat->second==0 || nextMat->second==binIDMat->second) ) {
+ } else if (binIDMat && binIDMat->second>0 && (nextMat->second==0 || nextMat->second==binIDMat->second) ) {
// exit from active layer
cache.m_identifiedParameters->push_back(std::pair<const Trk::TrackParameters*,int> (cPar->clone(),-binIDMat->second));
} else if (nextMat && nextMat->second>0) { // entry active layer
@@ -1695,12 +1695,12 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
delete cPar;
cache.m_currentLayerBin = dist2next.first;
- if (binIDMat!=nextMat) { // change of material triggers update of the cache
+ if (binIDMat!=nextMat) { // change of material triggers update of the cache
binIDMat = nextMat;
if (binIDMat) {
assert( cache.m_material);
updateMaterialEffects(cache,mom, sin(direction.theta()), sumPath+path);
- cache.m_material = binIDMat->first;
+ cache.m_material = binIDMat->first;
}
}
// recalculate distance to next bin
@@ -1725,7 +1725,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if ( mom < cache.m_bremEmitThreshold && m_simMatUpdator) {
// ST : strictly speaking, the emission point should be shifted backwards a bit (mom-m_bremEmitThreshold)
// this seems to be a minor point
- m_simMatUpdator->recordBremPhoton(cache.m_timeIn+cache.m_timeOfFlight+cache.m_timeStep,
+ m_simMatUpdator->recordBremPhoton(cache.m_timeIn+cache.m_timeOfFlight+cache.m_timeStep,
mom, cache.m_bremMom, position, direction, cache.m_particle);
cache.m_bremEmitThreshold = 0.;
}
@@ -1755,14 +1755,14 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
distanceEst = distSol.second();
}
// Peter Kluit: avoid jumping into other (distSol.first->second) distance solution for start surface with negative distance solution
- // negative distanceEst will trigger flipDirection = true and will iterate to the start surface
- // this will lead to very similar positions for multiple propagator calls and many tiny X0 scatterers
+ // negative distanceEst will trigger flipDirection = true and will iterate to the start surface
+ // this will lead to very similar positions for multiple propagator calls and many tiny X0 scatterers
if(ic==startSf&&distanceEst<0&&distSol.first()>0) distanceEst = distSol.first();
}
// eliminate close surface if path too small
if (ic==nextSf && fabs(distanceEst)<tol && fabs(path)<tol) {
(*vsIter).first=-1; vsIter=vsBeg; restart=true; distanceToTarget=maxPath; nextSf=-1;
- continue;
+ continue;
}
//If h and distance are in opposite directions, target is passed. Flip propagation direction
@@ -1779,12 +1779,12 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
helpSoft = fmax(0.05,1.-0.05*(*vsIter).first);
if ((*vsIter).first>20) helpSoft=1./(*vsIter).first;
}
- // take care of eliminating when number of flips even - otherwise it may end up at the start !
+ // take care of eliminating when number of flips even - otherwise it may end up at the start !
if ((*vsIter).first>50 && h*propDir>0 ) {
// fabs(distanceEst) >= fabs(previousDistance) ) {
(*vsIter).first = -1; vsIter = vsBeg; restart = true;
- continue;
- }
+ continue;
+ }
if ((*vsIter).first!=-1) flipDirection = true;
} else if ( fabs((*vsIter).second.second)>tol && fabs(distSol.currentDistance(true))>tol ) {
// here we need to compare with distance from current closest
@@ -1806,7 +1806,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
}
} else if (ic==nextSf) {
vsIter = vsBeg; restart = true;
- continue;
+ continue;
}
}
@@ -1816,7 +1816,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
//std::cout <<"iterating over surfaces:"<<vsIter-vsBeg<<","<<(*vsIter).second.first<<","<<(*vsIter).second.second<< std::endl;
// find closest surface: the step may have been beyond several surfaces
- // from all surfaces with 'negative' distance, consider only the one currently designed as 'closest'
+ // from all surfaces with 'negative' distance, consider only the one currently designed as 'closest'
// mw if ((*vsIter).first!=-1 && ( distanceEst>-tol || ic==nextSf ) ) {
if ((*vsIter).first!=-1 && ( distanceEst > 0. || ic==nextSf ) ) {
numSf++;
@@ -1850,10 +1850,10 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
}
// additional protection - return to the same surface
// eliminate the surface and restart the search
- // 04/10/10 ST:infinite loop due to distanceTolerance>tol fixed;
+ // 04/10/10 ST:infinite loop due to distanceTolerance>tol fixed;
if ( fabs(distanceToTarget)<=distanceTolerance && path*propDir<distanceTolerance ) {
(*vsIter).first = -1; vsIter = vsBeg; restart = true;
- continue;
+ continue;
}
sIter++; ic++;
}
@@ -1874,7 +1874,7 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
//don't step beyond bin boundary - adjust step
if (cache.m_binMat && fabs( h) > fabs(distanceToNextBin)+0.001 ) {
- if ( distanceToNextBin>0 ) { // TODO : investigate source of negative distance in BinningData
+ if ( distanceToNextBin>0 ) { // TODO : investigate source of negative distance in BinningData
//std::cout <<"adjusting step because of bin boundary:"<< h<<"->"<< distanceToNextBin*propDir<< std::endl;
h = distanceToNextBin*propDir;
}
@@ -1883,11 +1883,11 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
if ( helpSoft<1.) h*=helpSoft;
// don't step much beyond path limit
- if (cache.m_propagateWithPathLimit>0 && h > cache.m_pathLimit ) h = cache.m_pathLimit+tol;
+ if (cache.m_propagateWithPathLimit>0 && h > cache.m_pathLimit ) h = cache.m_pathLimit+tol;
//std::cout <<"current closest estimate: distanceToTarget: step size :"<< nextSf<<":"<< distanceToTarget <<":"<<h<< std::endl;
- //Abort if maxPath is reached
+ //Abort if maxPath is reached
if (fabs( path) > maxPath) return false;
if (steps++ > m_maxSteps) return false; //Too many steps, something is wrong
@@ -1925,12 +1925,12 @@ Trk::STEP_Propagator::propagateWithJacobian (Cache& cache,
// collect all surfaces with distance below tolerance
std::vector< std::pair<int,std::pair<double,double> > >::iterator vsIter = vsBeg;
- int index = 0;
- for ( ; vsIter!= vsEnd; vsIter++) {
- if ( (*vsIter).first != -1 && propDir*(*vsIter).second.first>=propDir*distanceToTarget-tol &&
+ int index = 0;
+ for ( ; vsIter!= vsEnd; vsIter++) {
+ if ( (*vsIter).first != -1 && propDir*(*vsIter).second.first>=propDir*distanceToTarget-tol &&
propDir*(*vsIter).second.first < 0.01 && index!= nextSf) {
solutions.push_back(index);
- }
+ }
if (index==nextSf) solutions.push_back(index);
index++;
}
@@ -1958,7 +1958,7 @@ Trk::STEP_Propagator::rungeKuttaStep( Cache& cache,
double& h,
double* P,
double* dDir,
- float* BG1,
+ double* BG1,
bool& firstStep,
double& distanceStepped) const
{
@@ -1980,28 +1980,34 @@ Trk::STEP_Propagator::rungeKuttaStep( Cache& cache,
double initialMomentum = fabs( 1./P[6]); // Set initial momentum
Amg::Vector3D initialPos( P[0], P[1], P[2]); // Set initial values for position
Amg::Vector3D initialDir( P[3], P[4], P[5]); // Set initial values for direction.
+// Directions at the different points. Used by the error propagation
Amg::Vector3D dir1;
Amg::Vector3D dir2;
Amg::Vector3D dir3;
- Amg::Vector3D dir4; // Directions at the different points. Used by the error propagation
- float BG23[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz, dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz
- float BG4[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // The gradients are used by the error propagation
- double g = 0.; // Energyloss in Mev/mm.
+ Amg::Vector3D dir4;
+ // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz, dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy,dBz/dz
+ double BG23[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
+ // The gradients are used by the error propagation
+ double BG4[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
+ double g = 0.; // Energyloss in Mev/mm.
double dgdl = 0.; // dg/dlambda in Mev/mm.
double particleMass = s_particleMasses.mass[cache.m_particle]; //Get particle mass from ParticleHypothesis
int steps = 0;
- //POINT 1. Only calculate this once per step, even if step is rejected. This point is independant of the step length, h
- double momentum = initialMomentum; // Current momentum
+ // POINT 1. Only calculate this once per step, even if step is rejected. This
+ // point is independant of the step length, h
+ double momentum = initialMomentum; // Current momentum
if (m_energyLoss && cache.m_matPropOK) {
- g = dEds( cache,momentum); //Use the same energy loss throughout the step.
- double E = std::sqrt( momentum*momentum + particleMass*particleMass);
- dP1 = g*E/momentum;
+ g = dEds(cache, momentum); // Use the same energy loss throughout the step.
+ double E = std::sqrt(momentum * momentum + particleMass * particleMass);
+ dP1 = g * E / momentum;
if (errorPropagation) {
if (m_includeGgradient) {
- dgdl = dgdlambda(cache, lambda1); //Use this value throughout the step.
+ dgdl = dgdlambda(cache, lambda1); // Use this value throughout the step.
}
- dL1 = -lambda1*lambda1*g*E*(3.-(momentum*momentum)/(E*E)) - lambda1*lambda1*lambda1*E*dgdl;
+ dL1 =
+ -lambda1 * lambda1 * g * E * (3. - (momentum * momentum) / (E * E)) -
+ lambda1 * lambda1 * lambda1 * E * dgdl;
}
}
@@ -2010,7 +2016,7 @@ Trk::STEP_Propagator::rungeKuttaStep( Cache& cache,
dir1 = dir;
if (firstStep) { // Poll BG1 if this is the first step, else use recycled BG4
firstStep = false;
- // MT Field cache is stored in cache
+ //Get the gradients needed for the error propagation if errorPropagation=true
getMagneticField( cache, pos, errorPropagation, BG1); //Get the gradients needed for the error propagation if errorPropagation=true
}
// Lorentz force, d2r/ds2 = lambda * (dir x B)
@@ -2034,7 +2040,6 @@ Trk::STEP_Propagator::rungeKuttaStep( Cache& cache,
dir = initialDir + (h/2.)*k1;
dir2 = dir;
- // MT Field cache is stored in cache
getMagneticField( cache, pos, errorPropagation, BG23);
Amg::Vector3D k2( dir.y()*BG23[2] - dir.z()*BG23[1], dir.z()*BG23[0] - dir.x()*BG23[2],
dir.x()*BG23[1] - dir.y()*BG23[0]);
@@ -2270,7 +2275,7 @@ void
Trk::STEP_Propagator::getMagneticField( Cache& cache,
const Amg::Vector3D& position,
bool getGradients,
- float* BG) const
+ double* BG) const
{
double pos[3] = {0.}; //Auxiliary variable, needed for fieldGradient_XYZ_in_mm.
pos[0] = position.x();
@@ -2283,31 +2288,27 @@ Trk::STEP_Propagator::getMagneticField( Cache& cache,
if (getGradients && m_includeBgradients) { // field gradients needed and available
- // MT Field cache is stored in cache
getFieldGradient(cache, R,H,dH);
-
- BG[0]=H[0]*magScale;
- BG[1]=H[1]*magScale;
- BG[2]=H[2]*magScale;
- BG[3] =dH[0]*magScale;
- BG[4] =dH[1]*magScale;
- BG[5] =dH[2]*magScale;
- BG[6] =dH[3]*magScale;
- BG[7] =dH[4]*magScale;
- BG[8] =dH[5]*magScale;
- BG[9] =dH[6]*magScale;
- BG[10]=dH[7]*magScale;
- BG[11]=dH[8]*magScale;
+ BG[0]=H[0]*magScale;
+ BG[1]=H[1]*magScale;
+ BG[2]=H[2]*magScale;
+ BG[3] =dH[0]*magScale;
+ BG[4] =dH[1]*magScale;
+ BG[5] =dH[2]*magScale;
+ BG[6] =dH[3]*magScale;
+ BG[7] =dH[4]*magScale;
+ BG[8] =dH[5]*magScale;
+ BG[9] =dH[6]*magScale;
+ BG[10]=dH[7]*magScale;
+ BG[11]=dH[8]*magScale;
}
else { //Homogenous field or no gradients needed, only retrieve the field strength.
- // MT Field cache is stored in cache
- getField(cache, R,H);
-
+ getField(cache, R,H);
BG[0]=H[0]*magScale;
- BG[1]=H[1]*magScale;
- BG[2]=H[2]*magScale;
+ BG[1]=H[1]*magScale;
+ BG[2]=H[2]*magScale;
for (int i=3; i<12; i++) { //Set gradients to zero
BG[i] = 0.;
@@ -2332,11 +2333,11 @@ Trk::STEP_Propagator::distance (Surface::SurfaceType surfaceType,
if (surfaceType == Trk::Surface::Cylinder)
return Trk::RungeKuttaUtils::stepEstimatorToCylinder(
targetSurface, P, distanceEstimationSuccessful);
-
+
if (surfaceType == Trk::Surface::Line || surfaceType == Trk::Surface::Perigee)
return Trk::RungeKuttaUtils::stepEstimatorToStraightLine(
targetSurface, P, distanceEstimationSuccessful);
-
+
if (surfaceType == Trk::Surface::Cone)
return Trk::RungeKuttaUtils::stepEstimatorToCone(targetSurface, P,
distanceEstimationSuccessful);
@@ -2353,8 +2354,8 @@ Trk::STEP_Propagator::distance (Surface::SurfaceType surfaceType,
double Trk::STEP_Propagator::dgdlambda( Cache& cache,double l) const
{
if (cache.m_particle == Trk::geantino || cache.m_particle == Trk::nonInteractingMuon) return 0.;
- if (cache.m_material->x0()==0 || cache.m_material->averageZ()==0) return 0.;
- if (cache.m_material->zOverAtimesRho()==0) return 0.;
+ if (cache.m_material->x0()==0 || cache.m_material->averageZ()==0) return 0.;
+ if (cache.m_material->zOverAtimesRho()==0) return 0.;
double p = fabs( 1./l);
double m = s_particleMasses.mass[cache.m_particle];
@@ -2406,7 +2407,7 @@ double Trk::STEP_Propagator::dgdlambda( Cache& cache,double l) const
// Multiple scattering and straggling contribution to the covariance matrix
/////////////////////////////////////////////////////////////////////////////////
-void Trk::STEP_Propagator::covarianceContribution( Cache& cache,
+void Trk::STEP_Propagator::covarianceContribution( Cache& cache,
const Trk::TrackParameters* trackParameters,
double path,
const Trk::TrackParameters* targetParms,
@@ -2427,7 +2428,7 @@ void Trk::STEP_Propagator::covarianceContribution( Cache& cache,
*measurementCovariance += *localMSCov;
- delete localMSCov;
+ delete localMSCov;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -2476,17 +2477,17 @@ void Trk::STEP_Propagator::dumpMaterialEffects( Cache& cache,
cache.m_combinedEloss.meanIoni(), cache.m_combinedEloss.sigmaIoni(),
cache.m_combinedEloss.meanRad(), cache.m_combinedEloss.sigmaRad(), path ) ;
- Trk::ScatteringAngles* sa = new Trk::ScatteringAngles(0.,0.,std::sqrt(cache.m_covariance(2,2)),
- std::sqrt(cache.m_covariance(3,3)));
+ Trk::ScatteringAngles* sa = new Trk::ScatteringAngles(0.,0.,std::sqrt(cache.m_covariance(2,2)),
+ std::sqrt(cache.m_covariance(3,3)));
Trk::CurvilinearParameters* cvlTP = parms->clone();
Trk::MaterialEffectsOnTrack* mefot = new Trk::MaterialEffectsOnTrack(cache.m_combinedThickness,sa,eloss,
- cvlTP->associatedSurface());
+ cvlTP->associatedSurface());
cache.m_matstates->push_back(new TrackStateOnSurface(nullptr,cvlTP,nullptr,mefot));
}
- cache.m_matdump_lastpath = path;
+ cache.m_matdump_lastpath = path;
// clean-up
cache.m_combinedCovariance += cache.m_covariance;
@@ -2499,14 +2500,14 @@ void Trk::STEP_Propagator::dumpMaterialEffects( Cache& cache,
// Multiple scattering and straggling contribution to the covariance matrix in global coordinates
/////////////////////////////////////////////////////////////////////////////////////////////////////////
-void Trk::STEP_Propagator::updateMaterialEffects( Cache& cache,
- double mom,
+void Trk::STEP_Propagator::updateMaterialEffects( Cache& cache,
+ double mom,
double sinTheta,
double path ) const
{
// collects material effects in between material dumps ( m_covariance )
- // collect straggling
+ // collect straggling
if (m_straggling) {
cache.m_covariance(4,4) += cache.m_stragglingVariance;
cache.m_stragglingVariance = 0.;
@@ -2515,7 +2516,7 @@ void Trk::STEP_Propagator::updateMaterialEffects( Cache& cache,
if (!m_multipleScattering) return;
double totalMomentumLoss = mom - cache.m_matupd_lastmom;
- double pathSinceLastUpdate = path - cache.m_matupd_lastpath;
+ double pathSinceLastUpdate = path - cache.m_matupd_lastpath;
double pathAbs = fabs(pathSinceLastUpdate);
@@ -2523,8 +2524,8 @@ void Trk::STEP_Propagator::updateMaterialEffects( Cache& cache,
double matX0 = cache.m_material->x0();
- // calculate number of layers : TODO : thickness to be adjusted
- int msLayers = int(pathAbs/m_layXmax/matX0)+1 ;
+ // calculate number of layers : TODO : thickness to be adjusted
+ int msLayers = int(pathAbs/m_layXmax/matX0)+1 ;
double massSquared = s_particleMasses.mass[cache.m_particle]*s_particleMasses.mass[cache.m_particle];
double layerThickness = pathAbs/msLayers;
@@ -2577,7 +2578,7 @@ void Trk::STEP_Propagator::updateMaterialEffects( Cache& cache,
if(!useCache) {
cumulatedX0 = cumulatedVariance/C0;
if (cumulatedX0>0.001) {
- double lX0 = log(cumulatedX0);
+ double lX0 = log(cumulatedX0);
cumulatedX0 = cumulatedX0/(1+2*0.038*lX0+0.038*0.038*lX0*lX0);
}
}
@@ -2586,8 +2587,8 @@ void Trk::STEP_Propagator::updateMaterialEffects( Cache& cache,
thetaVariance = C0*MS2 + C1*MS2s*log(MS2s) + C2*MS2s*log(MS2s)*log(MS2s);
if (cumulatedX0>0.001) {
- double lX0 = log(cumulatedX0);
- thetaVariance += -C1*cumulatedX0*lX0 - C2*cumulatedX0*lX0*lX0;
+ double lX0 = log(cumulatedX0);
+ thetaVariance += -C1*cumulatedX0*lX0 - C2*cumulatedX0*lX0*lX0;
}
//Calculate ms covariance contributions and scale
@@ -2620,12 +2621,12 @@ Trk::TrackParameters* Trk::STEP_Propagator::createStraightLine( const Trk::Track
// ATH_MSG_VERBOSE("STEP propagator detects invalid input parameters (q/p=0 ), resetting momentum to 1.e10");
if (inputTrackParameters->type()==Trk::Curvilinear) {
- return new Trk::CurvilinearParameters(inputTrackParameters->position(), lp[2], lp[3], lp[4],
- (inputTrackParameters->covariance() ?
+ return new Trk::CurvilinearParameters(inputTrackParameters->position(), lp[2], lp[3], lp[4],
+ (inputTrackParameters->covariance() ?
new AmgSymMatrix(5)(*inputTrackParameters->covariance()) : nullptr ) );
- }
+ }
return inputTrackParameters->associatedSurface().createTrackParameters(lp[0], lp[1], lp[2], lp[3], lp[4],
- (inputTrackParameters->covariance() ?
+ (inputTrackParameters->covariance() ?
new AmgSymMatrix(5)(*inputTrackParameters->covariance()) : nullptr ) );
}
@@ -2639,10 +2640,10 @@ double Trk::STEP_Propagator::dEds( Cache& cache, double p) const
cache.m_delIoni = 0.; cache.m_delRad = 0.; cache.m_kazL = 0.;
if (cache.m_particle == Trk::geantino || cache.m_particle == Trk::nonInteractingMuon) return 0.;
- if (cache.m_material->x0()==0 || cache.m_material->averageZ()==0) return 0.;
+ if (cache.m_material->x0()==0 || cache.m_material->averageZ()==0) return 0.;
- cache.m_delIoni = cache.m_matInt.dEdl_ionization(p, cache.m_material, cache.m_particle,
- cache.m_sigmaIoni, cache.m_kazL);
+ cache.m_delIoni = cache.m_matInt.dEdl_ionization(p, cache.m_material, cache.m_particle,
+ cache.m_sigmaIoni, cache.m_kazL);
cache.m_delRad = cache.m_matInt.dEdl_radiation(p, cache.m_material, cache.m_particle, cache.m_sigmaRad);
@@ -2654,10 +2655,10 @@ double Trk::STEP_Propagator::dEds( Cache& cache, double p) const
// Smear momentum ( simulation mode )
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-void Trk::STEP_Propagator::smear(Cache& cache,
- double& phi,
- double& theta,
- const Trk::TrackParameters* parms,
+void Trk::STEP_Propagator::smear(Cache& cache,
+ double& phi,
+ double& theta,
+ const Trk::TrackParameters* parms,
double radDist) const
{
if (cache.m_particle == Trk::geantino) return ;
@@ -2694,11 +2695,11 @@ void Trk::STEP_Propagator::sampleBrem(Cache& cache, double mom) const
double rndx = CLHEP::RandFlat::shoot(m_randomEngine);
double rnde = CLHEP::RandFlat::shoot(m_randomEngine);
- // sample visible fraction of the mother momentum taken according to 1/f
+ // sample visible fraction of the mother momentum taken according to 1/f
double eps = m_momentumCutOff/mom;
- cache.m_bremMom = pow(eps,pow(rndx,exp(1.)))*mom; // adjustment here ?
+ cache.m_bremMom = pow(eps,pow(rndx,exp(1.)))*mom; // adjustment here ?
cache.m_bremSampleThreshold = mom - cache.m_bremMom;
- cache.m_bremEmitThreshold = mom - rnde*cache.m_bremMom;
+ cache.m_bremEmitThreshold = mom - rnde*cache.m_bremMom;
}
Trk::TrackParameters*
@@ -2711,16 +2712,16 @@ Trk::STEP_Propagator::propagateNeutral(const Trk::TrackParameters& parm,
bool returnCurv) const {
// find nearest valid intersection
- double tol = 0.001;
+ double tol = 0.001;
solutions.clear();
std::vector<std::pair<unsigned int,double> > currentDist;
- currentDist.clear();
+ currentDist.clear();
std::vector<std::pair<const Trk::Surface*,Trk::BoundaryCheck> >::iterator sIter = targetSurfaces.begin();
std::vector<std::pair<unsigned int,double> >::iterator oIter = currentDist.begin();
std::vector<DestSurf >::iterator sBeg = targetSurfaces.begin();
- const Amg::Vector3D& position(parm.position());
- Amg::Vector3D direction(parm.momentum().normalized());
+ const Amg::Vector3D& position(parm.position());
+ Amg::Vector3D direction(parm.momentum().normalized());
for (; sIter!=targetSurfaces.end(); sIter++) {
Trk::DistanceSolution distSol = (*sIter).first->straightLineDistanceEstimate(position,direction);
@@ -2753,16 +2754,16 @@ Trk::STEP_Propagator::propagateNeutral(const Trk::TrackParameters& parm,
if (targetSurfaces[(*oIter).first].first->isOnSurface(xsct, (*oIter).second, 0.001, 0.001) ) {
if ( usePathLimit && path>0. && path<((*oIter).second) ) {
Amg::Vector3D endpoint = position+propDir*direction*path;
- return new Trk::CurvilinearParameters(endpoint,parm.momentum(),parm.charge());
+ return new Trk::CurvilinearParameters(endpoint,parm.momentum(),parm.charge());
}
path = (*oIter).second;
solutions.push_back((*oIter).first);
const Trk::Surface* sf = targetSurfaces[(*oIter).first].first;
- if( returnCurv || sf->type()==Trk::Surface::Cone) return new Trk::CurvilinearParameters(xsct,parm.momentum(),parm.charge());
+ if( returnCurv || sf->type()==Trk::Surface::Cone) return new Trk::CurvilinearParameters(xsct,parm.momentum(),parm.charge());
return sf->createTrackParameters(xsct,parm.momentum(),parm.charge(),nullptr);
- }
+ }
}
return nullptr;
diff --git a/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_v1Dev_build.ref b/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_v1Dev_build.ref
index 1ce4d677498b729b9de1380e7242bcd68e298bd9..8c12610dcf4e52e6329ba2dc6c465815a91fddee 100644
--- a/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_v1Dev_build.ref
+++ b/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_v1Dev_build.ref
@@ -2171,25 +2171,25 @@ HLT_mu50_L1MU20:
stepCounts:
0: 8
1: 6
- 2: 2
+ 2: 1
3: 1
stepFeatures:
0: 10
1: 7
- 2: 2
+ 2: 1
3: 1
HLT_mu50_RPCPEBSecondaryReadout_L1MU20:
eventCount: 1
stepCounts:
0: 8
1: 6
- 2: 2
+ 2: 1
3: 1
4: 1
stepFeatures:
0: 10
1: 7
- 2: 2
+ 2: 1
3: 1
4: 1
HLT_mu60_0eta105_msonly_L1MU20:
@@ -2295,14 +2295,14 @@ HLT_mu6_mu6noL1_L1MU6:
1: 10
2: 10
3: 10
- 4: 6
+ 4: 5
5: 3
stepFeatures:
0: 23
1: 22
2: 23
3: 23
- 4: 21
+ 4: 18
5: 11
HLT_mu6_xe30_mht_L1XE10:
eventCount: 9