diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc
index 2eea3be465e631fb54bede7f43334aabbd33723e..499610ea77e4290bfdd9518fe4139c181239225a 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc
@@ -33,7 +33,7 @@ template <class T> Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::hand
//EX_MSG_DEBUG(eCell.navigationStep, "Update", "char", "Update on layer with index " << m_layer->layerIndex().value() << " - corr factor = " << pathCorrection*mFactor);
const Trk::MaterialProperties* materialProperties = m_layer->layerMaterialProperties()->fullMaterial(eCell.leadParameters->position());
- if ( materialProperties ) {
+ if ( materialProperties && materialProperties->thicknessInX0()>0 ) {
m_matProp = materialProperties;
m_thicknessInX0 = mFactor*m_matProp->thicknessInX0();
m_thicknessInL0 = mFactor*m_matProp->thicknessInL0();
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/EnergyLossSamplerBetheHeitler.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/EnergyLossSamplerBetheHeitler.cxx
index ecf6a3d1dd8790b552f960a0cbd1579a2ba51ee1..48c96e36442cd53ed32ee8bdf5985a0e65153623 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/EnergyLossSamplerBetheHeitler.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/EnergyLossSamplerBetheHeitler.cxx
@@ -71,6 +71,8 @@ Trk::EnergyLoss* iFatras::EnergyLossSamplerBetheHeitler::energyLoss( const Trk::
double pathLength = pathCorrection*materialProperties.thicknessInX0();
+ if (pathLength==0.) return new Trk::EnergyLoss(0.,0.);
+
double p = pInitial;
double me = s_particleMasses.mass[Trk::electron];
double E = sqrt(p*p+me*me);
@@ -116,11 +118,12 @@ Trk::EnergyLoss* iFatras::EnergyLossSamplerBetheHeitler::energyLoss( const Trk::
double u = CLHEP::RandGamma::shoot(m_randomEngine, pathLength / log(2.), 1.);
double z = exp( -1. * u );
double deltaE(0.);
+
if ( direction == Trk::alongMomentum )
deltaE = m_scaleFactor*E * ( z - 1. );
else
deltaE = m_scaleFactor*E * ( 1. / z - 1. );
-
+
simulatedDeltaE+=fabs(deltaE);
// transform into negative energyloss
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx
index eb26e0e8b6022eb77bb093fb41b15731f32eedeb..c6f7530282cf66f07cbc8d630e668bab4dda43ef 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx
@@ -258,7 +258,7 @@ void iFatras::McMaterialEffectsEngine::multipleScatteringUpdate(const Trk::Track
// the initial values
double theta = parameters[Trk::theta];
double phi = parameters[Trk::phi];
- double sinTheta = (theta*theta > 10e-10) ? sin(theta) : 1.;
+ double sinTheta = (sin(theta)*sin(theta) > 10e-10) ? sin(theta) : 1.;
// sample them in an independent way
double deltaTheta = m_projectionFactor*simTheta;
@@ -464,7 +464,10 @@ Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer(
const Trk::TrackParameters* parm = eCell.leadParameters;
// path correction
- double pathCorrection = fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum())));
+ double pathCorrection = m_layer->surfaceRepresentation().normal().dot(parm->momentum()) !=0 ?
+ fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum()))) : 1.;
+
+ if (!pathCorrection) return Trk::ExtrapolationCode::InProgress;
// figure out if particle stopped in the layer and recalculate path limit
int iStatus = 0;
@@ -502,9 +505,9 @@ Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer(
double p = parm->momentum().mag();
double m = m_particleMasses.mass[eCell.pHypothesis];
double E = sqrt(p*p+m*m);
-
+
// radiation and ionization preceed the presampled interaction (if any)
-
+
if (m_eLoss || m_ms) {
// the updatedParameters - first a deep copy
AmgVector(5) uParameters = parm->parameters();
@@ -537,6 +540,8 @@ Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer(
uParameters[Trk::qOverP] = parm->charge()/newP;
}
+
+ EX_MSG_VERBOSE(eCell.navigationStep, "eloss", "char ", "E,deltaE:" <<E<<","<< eloss->deltaE() );
// TODO straggling
@@ -562,8 +567,8 @@ Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer(
}
- if ( m_ms ) {
-
+ if ( m_ms && m_thicknessInX0>0 ) {
+
double simTheta = m_msSampler->simTheta(*m_matProp, p, dX0/m_thicknessInX0, eCell.pHypothesis);
//do the update -> You need 2 evaluation of simTheta. The second one is used to calculate deltaphi in multipleScatteringUpdate
multipleScatteringUpdate(*(eCell.leadParameters), uParameters, simTheta,
@@ -739,7 +744,10 @@ Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer(
const Trk::NeutralParameters* parm = eCell.leadParameters;
// path correction
- double pathCorrection = fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum())));
+ double pathCorrection = m_layer->surfaceRepresentation().normal().dot(parm->momentum()) !=0 ?
+ fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum()))) : 1.;
+
+ if (!pathCorrection) return Trk::ExtrapolationCode::InProgress;
// figure out if particle stopped in the layer and recalculate path limit
int iStatus = 0;
@@ -966,6 +974,7 @@ void iFatras::McMaterialEffectsEngine::radiate(const ISF::ISFParticle* parent, A
recordBremPhotonLay(parent,eCell.time,p,deltaP,ePos,eDir,eCell.pHypothesis,dir,mFr);
p *=z ;
+ EX_MSG_VERBOSE("", "radiate", "", "brem photon emitted " << deltaP<<":updated e momentum:"<< p );
// std::cout <<"brem photon emitted, momentum update:"<< p<<","<<path<<","<<pathLim<<",mat.fraction:"<<mFr<<std::endl;
}
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx
index 5314eeb7f4a4491201e34c0e4f50e8bd028c17f5..4d1e558b8884499282e262f9447e3eaaa3dc309b 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx
@@ -431,7 +431,8 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
// recalculate if missing
m_matProp = m_matProp ? m_matProp : m_layer->fullUpdateMaterialProperties(*parm);
- double pathCorrection = fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),parm->momentum()));
+ double pathCorrection = m_layer->surfaceRepresentation().normal().dot(parm->momentum()) !=0 ?
+ fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum()))) : 1.;
// check if a bending correction has to be applied
if (m_bendingCorrection) {
@@ -477,7 +478,7 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
int iStatus = 0;
double dX0 = (1.-matFraction)*pathCorrection*m_matProp->thicknessInX0();
- if (!m_matProp || particle==Trk::geantino || particle==Trk::nonInteracting ) {
+ if (!m_matProp || particle==Trk::geantino || particle==Trk::nonInteracting || dX0==0.) {
// non-interacting - pass them back
pathLim.updateMat(dX0,m_matProp->averageZ(),dInL0);
@@ -537,7 +538,7 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
AmgVector(5) updatedParameters(parm->parameters());
const Trk::TrackParameters* updated = 0;
- if ( m_eLoss && particle==Trk::electron ) {
+ if ( m_eLoss && particle==Trk::electron && dX0>0.) {
double matTot = (1-matFraction)*pathCorrection*m_matProp->thicknessInX0();
Amg::Vector3D edir = parm->momentum().unit();
@@ -551,10 +552,10 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
matFraction += dX0/pathCorrection/m_matProp->thicknessInX0();
}
- if ( isp->charge()!=0 ) {
+ if ( isp->charge()!=0 && dX0>0.) {
if ( m_eLoss ) ionize(*parm, updatedParameters, dX0, dir, particle );
- if ( m_ms ) {
+ if ( m_ms && m_matProp->thicknessInX0()>0.) {
double sigmaMSproj = (m_use_msUpdator && m_msUpdator ) ?
sqrt(m_msUpdator->sigmaSquare(*m_matProp, p, dX0/m_matProp->thicknessInX0(), particle)) : msSigma(dX0,p,particle);
@@ -963,7 +964,7 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::update(double tim
m_deltaP = deltaP;
}
- if (m_ms){
+ if (m_ms && matprop.x0()>0 ){
// get the projected scattering angle
double sigmaMSproj = (m_use_msUpdator && m_msUpdator) ? sqrt(m_msUpdator->sigmaSquare(matprop, p, pathCorrection, particle)) :
msSigma(pathCorrection*matprop.thicknessInX0(),p,particle);
@@ -1558,7 +1559,9 @@ void iFatras::McMaterialEffectsUpdator::recordBremPhotonLay(const ISF::ISFPartic
// decide if photon leaves the layer
// amount of material to traverse
- double pathCorrection = m_layer->surfaceRepresentation().pathCorrection(bremPhoton->position(),bremPhoton->momentum());
+ double pathCorrection = m_layer->surfaceRepresentation().normal().dot(bremPhoton->momentum()) !=0 ?
+ fabs(m_layer->surfaceRepresentation().pathCorrection(bremPhoton->position(),bremPhoton->momentum())) : 1.;
+
double mTot = m_matProp->thicknessInX0()*pathCorrection*(1.-remMat);
if (mTot < pLim.x0Max) { // release