Commit 8d75f32a authored by Sarka Todorova's avatar Sarka Todorova Committed by Graeme Stewart
Browse files

attempt to fix invalid read (TrkExTools-03-01-11)

2014-12-18 Sharka Todorova <sarka.todorova@cern.ch>
	* attempt to fix invalid read (ATLRECTS-1433) by cleanup at the entry
	* Tagged as TrkExTools-03-01-11

2014-12-11  Massimiliano Bellomo  <massimiliano.bellomo@cern.ch>
	* EnergyLossUpdator fix for creating extended EnergyLoss with correct uncertainties
	* tagged as TrkExTools-03-01-10

2014-12-09 Andreas Salzburger < Andreas.Salzburger -at- cern.ch >
    * add boundary surfaces to extrapolateM()
    * tagged as TrkExTools-03-01-09

2014-12-08 Andreas Salzburger < Andreas.Salzburger -at- cern.ch >
    * plug memory leak, part 2
    * tagges as TrkExTools-03-01-08

2014-12-05 Andreas Salzburger < Andreas.Salzburger -at- cern.ch >
    * plug memory leak
    * Tagged as TrkExTools-03-01-07

...
(Long ChangeLog diff - truncated)
parent 1d8c30fc
......@@ -140,6 +140,14 @@ namespace Trk {
bool mpv = false) const;
/** Method to recalculate Eloss values for the fit setting an elossFlag using as an input
the detailed Eloss information Calorimeter energy, error momentum and momentum error */
EnergyLoss* updateEnergyLoss(EnergyLoss* eLoss, double caloEnergy, double caloEnergyError,
double pCaloEntry, double momentumError, int & elossFlag) const;
/** Routine to calculate X0 and Eloss scale factors for the Calorimeter and Muon System */
void getX0ElossScales(int icalo, double eta, double phi, double & X0Scale, double & ElossScale ) const;
/** Method to return the variance of the change in q/p for the Bethe-Heitler parameterisation */
double varianceDeltaQoverP(const MaterialProperties&,
double p,
......@@ -147,13 +155,6 @@ namespace Trk {
PropDirection direction = alongMomentum,
ParticleHypothesis particleHypothesis = electron ) const;
/** Method to recalculate Eloss values for the fit using as an input
the detailed Eloss information Calorimeter energy, error and momentum error*/
EnergyLoss* updateEnergyLoss(
EnergyLoss* eLoss, double caloEnergy, double caloEnergyError, double momentumError) const;
private:
......
......@@ -821,7 +821,6 @@ namespace Trk {
// ------------------------------- static members --------------------------------------------------------------------
static double s_distIncreaseTolerance; //!< distance increatse tolerance to account for straight line approx.
static double s_distEntryLayerMax; //!< maximal allowed distance to the entry layer
unsigned int m_maxNavigSurf;
unsigned int m_maxNavigVol;
......
......@@ -40,6 +40,7 @@ namespace Trk {
class Layer;
class Volume;
class DetachedTrackingVolume;
class AlignableTrackingVolume;
class TrackingGeometry;
class TrackParticleBase;
class IPropagator;
......@@ -63,6 +64,30 @@ namespace Trk {
{}
};
struct BoundaryTrackParameters
{
const TrackParameters* trPar;
const TrackingVolume* exitVol;
const TrackingVolume* entryVol;
//
BoundaryTrackParameters( const TrackParameters* parms,
const TrackingVolume* exitTV,
const TrackingVolume* entryTV ) :
trPar(parms), exitVol(exitTV), entryVol(entryTV)
{}
};
struct IdentifiedIntersection
{
float distance;
int identifier;
const Trk::Material* material;
//
IdentifiedIntersection( float dist, int id, const Trk::Material* mat) :
distance(dist), identifier(id), material(mat)
{}
};
/** @struct ParametersAtBoundarySurface
has only three member
- BoundarySurface
......@@ -179,6 +204,14 @@ namespace Trk {
Trk::GeometrySignature& nextGeoID,
const Trk::TrackingVolume* destVol) const;
BoundaryTrackParameters extrapolateInAlignableTV(
const Trk::TrackParameters& parm,
Trk::TimeLimit& time,
Trk::PropDirection dir,
Trk::ParticleHypothesis particle,
Trk::GeometrySignature& nextGeoId,
const Trk::AlignableTrackingVolume* aliTV) const;
const Trk::TrackParameters* transportToVolumeWithPathLimit(
const Trk::TrackParameters& parm,
......@@ -188,6 +221,14 @@ namespace Trk {
Trk::GeometrySignature& nextGeoId,
const Trk::TrackingVolume* boundaryVol) const;
BoundaryTrackParameters transportInAlignableTV(
const Trk::TrackParameters& parm,
Trk::TimeLimit& time,
Trk::PropDirection dir,
Trk::ParticleHypothesis particle,
Trk::GeometrySignature& nextGeoId,
const Trk::AlignableTrackingVolume* aliTV) const;
/** Access the subPropagator to the given volume*/
const IPropagator* subPropagator(const TrackingVolume& tvol) const;
......@@ -311,6 +352,7 @@ namespace Trk {
mutable bool m_robustSampling;
mutable PathLimit m_path;
mutable double m_time;
mutable size_t m_currentLayerBin;
//------------------------- NAVIGATION -------- ----------------------------------------------//
mutable int m_methodSequence;
......
......@@ -85,7 +85,7 @@ const Trk::TrackParameters* Trk::DummyMaterialEffectsUpdator::update(const Trac
if (tvol && updateProperties){
double correctionFactor = lay.pathCorrection(*parm);
double correctionFactor = fabs(lay.surfaceRepresentation().pathCorrection(parm->position(),parm->momentum()));
// material properties
double pathInX0 = updateProperties->thicknessInX0();
double x0 = updateProperties->x0();
......@@ -131,7 +131,7 @@ const Trk::TrackParameters* Trk::DummyMaterialEffectsUpdator::preUpdate(const T
// retrun if the preFactor is too small
if (preFactor < 0.1) return (parm);
double correctionFactor = lay.pathCorrection(*parm);
double correctionFactor = fabs(lay.surfaceRepresentation().pathCorrection(parm->position(),parm->momentum()));
bool surfaceValidaiton = false;
// get the material properties
......@@ -188,7 +188,7 @@ const Trk::TrackParameters* Trk::DummyMaterialEffectsUpdator::postUpdate(const
if (postFactor < 0.1 ) return 0;
// the correction Factor to the layer
double correctionFactor = lay.pathCorrection(parm);
double correctionFactor = fabs(lay.surfaceRepresentation().pathCorrection(parm.position(),parm.momentum()));
// get the material properties
const Trk::MaterialProperties* updateProperties = lay.fullUpdateMaterialProperties(parm);
correctionFactor *= postFactor;
......@@ -241,6 +241,6 @@ const Trk::TrackParameters* Trk::DummyMaterialEffectsUpdator::update(const Trac
void Trk::DummyMaterialEffectsUpdator::validationAction() const
{
// first record the values
if (m_validationMode && m_materialMapper)
m_materialMapper->finalizeEvent(m_etaFinalize, m_phiFinalize);
// if (m_validationMode && m_materialMapper)
// m_materialMapper->finalizeEvent(m_etaFinalize, m_phiFinalize);
}
......@@ -41,7 +41,7 @@ Trk::EnergyLossUpdator::EnergyLossUpdator(const std::string& t, const std::strin
m_stragglingErrorScale(1.),
m_mpvScale(0.98),
m_mpvSigmaParametric(false),
m_detailedEloss(false)
m_detailedEloss(true)
{
declareInterface<Trk::IEnergyLossUpdator>(this);
// scale from outside
......@@ -128,7 +128,7 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::energyLoss(const MaterialProperties& ma
m_transTmax = 0.;
// preparation
double sign = (dir==Trk::oppositeMomentum) ? 1. : -1.;
double sign = (dir==Trk::oppositeMomentum) ? -1. : 1.;
double pathLength = pathcorrection * mat.thicknessInX0()*mat.x0();
......@@ -138,8 +138,8 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::energyLoss(const MaterialProperties& ma
double meanIoni = m_matInt.dEdl_ionization(p, &(mat.material()), particle, sigIoni, kazL);
double meanRad = m_matInt.dEdl_radiation(p, &(mat.material()), particle, sigRad);
meanIoni = pathLength*meanIoni;
meanRad = pathLength*meanRad;
meanIoni = sign*pathLength*meanIoni;
meanRad = sign*pathLength*meanRad;
sigIoni = pathLength*sigIoni;
sigRad = pathLength*sigRad;
kazL = pathLength*kazL;
......@@ -150,11 +150,10 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::energyLoss(const MaterialProperties& ma
deltaE = meanIoni + meanRad;
sigmaDeltaE = sqrt(sigIoni*sigIoni+sigRad*sigRad);
ATH_MSG_DEBUG( " Energy loss updator deltaE " << deltaE << " meanIoni " << meanIoni << " meanRad " << meanRad << " sigIoni " << sigIoni << " sigRad " << sigRad );
ATH_MSG_DEBUG( " Energy loss updator deltaE " << deltaE << " meanIoni " << meanIoni << " meanRad " << meanRad << " sigIoni " << sigIoni << " sigRad " << sigRad << " sign " << sign << " pathLength " << pathLength );
Trk::EnergyLoss* eloss = !m_detailedEloss ? new Trk::EnergyLoss(deltaE,sigmaDeltaE):
new Trk::EnergyLoss(deltaE, sigmaDeltaE, meanIoni, sigIoni, meanRad, sigRad );
new Trk::EnergyLoss(deltaE, sigmaDeltaE, sigmaDeltaE, sigmaDeltaE, meanIoni, sigIoni, meanRad, sigRad, pathLength );
if(m_useTrkUtils) return eloss;
// Code below will not be used if the parameterization of TrkUtils is used
......@@ -260,9 +259,9 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::energyLoss(const MaterialProperties& ma
sigmaDeltaE = sqrt( sigmaDeltaE_rad*sigmaDeltaE_rad + sigmaDeltaE_ioni*sigmaDeltaE_ioni);
return ( m_detailedEloss ? new EnergyLoss(deltaE, sigmaDeltaE,
return ( m_detailedEloss ? new EnergyLoss(deltaE, sigmaDeltaE, sigmaDeltaE, sigmaDeltaE,
(mpvSwitch? deltaE_ioni :0.9*deltaE_ioni), sigmaDeltaE_ioni,
deltaE_rad, sigmaDeltaE_rad ):
deltaE_rad, sigmaDeltaE_rad, pathLength ):
new EnergyLoss(deltaE, sigmaDeltaE) );
}
......@@ -343,20 +342,29 @@ double Trk::EnergyLossUpdator::dEdXBetheHeitler(const MaterialProperties& mat,
return initialE/mat.x0()*mfrac;
}
Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss, double caloEnergy, double caloEnergyError, double momentumError) const {
// public interface method
Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss, double caloEnergy, double caloEnergyError, double pCaloEntry, double momentumError, int & elossFlag) const {
//
// Input the detailed EnergyLoss object that contains the different Eloss terms and their uncertainties
// and the momentumError (all in MeV)
// Input: the detailed EnergyLoss object in the Calorimeter that contains the different Eloss terms
// and their uncertainties; caloEnergy and error; and the muon momentumError (all in MeV)
//
// For use in the MuonSystem
// Input: caloEnergy = 0. caloEnergyError = 0. and pCaloEntry = pMuonEntry momentum at MuonEntry
//
// Output: an updated Energy loss values deltaE()
// that can be used in the track fit and corresponds to the Most Probable EnergyLoss value
// taking into account the ionization, radiation and
// smearing due to the errors including the momentumError (in MeV)
//
// elossFlag = false if Calorimeter Energy is NOT stored (and later fitted) on the Eloss object
// = true Calorimeter Energy is stored and will be fitted
//
// deltaE is used in the final fit
//
elossFlag = 0;
int isign = 1;
if(eLoss->deltaE()<0) isign = -1;
......@@ -364,22 +372,26 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss
double sigmaDeltaE = eLoss->sigmaDeltaE();
// Detailed Eloss
double deltaE_ioni = eLoss->meanIoni();
double sigmaDeltaE_ioni = eLoss->sigmaIoni(); // sigma Landau
// double deltaE_rad = eLoss->meanRad();
double sigmaDeltaE_ioni = 0.45*eLoss->sigmaIoni(); // sigma Landau
double deltaE_rad = eLoss->meanRad();
double sigmaDeltaE_rad = eLoss->sigmaRad(); // rms and mean of steep exponential
double depth = eLoss->length();
double sigmaPlusDeltaE = eLoss->sigmaPlusDeltaE();
double sigmaMinusDeltaE = eLoss->sigmaMinusDeltaE();
double MOP = deltaE_ioni - isign*3.59524*sigmaDeltaE_ioni;
// std::cout << " update Energyloss old deltaE " << deltaE << " sigmaPlusDeltaE " << sigmaPlusDeltaE << " sigmaMinusDeltaE " << sigmaMinusDeltaE << std::endl;
//
// MOP shift due to ionization and radiation
//
double MOPshift = isign*0.05*sqrt(sigmaDeltaE_ioni*sigmaDeltaE_rad);
double MOPshift = isign*50*10000./pCaloEntry + isign*0.75*sqrt(sigmaDeltaE_ioni*sigmaDeltaE_rad);
//
// define sigmas for Landau convoluted with exponential
//
double sigmaL = sigmaDeltaE_ioni + sigmaDeltaE_rad/3.59524;
double fracErad = sigmaDeltaE_rad + fabs(deltaE_rad)*pCaloEntry/(800000.+pCaloEntry);
double sigmaL = sigmaDeltaE_ioni + fracErad/3.59524;
double sigmaMinus = 1.02*sigmaDeltaE_ioni + 0.08*sigmaDeltaE_rad;
double sigmaPlus = 4.65*sigmaDeltaE_ioni + 1.16*sigmaDeltaE_rad;
......@@ -391,7 +403,8 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss
//
// Shift of MOP due to momentum resolution
//
double xc = 0.87388*momentumError/(3.59524*sigmaL);
double scale_xc = 2.3;
double xc = scale_xc*0.87388*momentumError/(3.59524*sigmaL);
double correction = (1.747*xc*xc + 0.97*0.938*xc*xc*xc)/(1+4.346*xc+5.371*xc*xc+0.938*xc*xc*xc); // correction ranges from 0 to 0.97
double MOPreso = isign*3.59524*sigmaL*correction;
......@@ -412,10 +425,11 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss
sigmaMinusDeltaE = caloEnergyError;
sigmaPlusDeltaE = caloEnergyError;
sigmaDeltaE = caloEnergyError;
elossFlag = 1;
} else {
// Use MOp after corrections
// Use MOP after corrections
//
// MOPCalo is correction to MOP for Calorimeter energy cut
......@@ -441,8 +455,89 @@ Trk::EnergyLoss* Trk::EnergyLossUpdator::updateEnergyLoss(Trk::EnergyLoss* eLoss
}
// std::cout << " update Energyloss new deltaE " << deltaE << " sigmaPlusDeltaE " << sigmaPlusDeltaE << " sigmaMinusDeltaE " << sigmaMinusDeltaE << std::endl;
return new EnergyLoss(deltaE, sigmaDeltaE, sigmaMinusDeltaE, sigmaPlusDeltaE, deltaE_ioni, sigmaDeltaE_ioni, deltaE_rad, sigmaDeltaE_rad, depth );
}
// public interface method
void Trk::EnergyLossUpdator::getX0ElossScales(int icalo, double eta, double phi, double & X0Scale, double & ElossScale ) const {
//
// for Calorimeter icalo = 1
// Muon System icalo = 0
// convention eta, phi is at Calorimeter Exit (or Muon Entry)
// eta and phi are from the position (not direction)
//
// input X0 and ElossScale = 1
// output updated X0Scale and ElossScale
//
double X0CaloGirder[4] = {-1.02877e-01, -2.74322e-02, 8.12989e-02,9.73551e-01};
double X0CaloScale[60] = {1.03178 ,1.03178 ,1.03178 ,1.03178 ,1.01512 ,1.01456 ,1.01266 ,1.02299 ,1.01437 ,1.01561 ,1.01731 ,1.02502 ,1.01171 ,0.987469 ,0.979134 ,1.01018 ,1.01813 ,0.986147 ,0.973686 ,0.983771 ,0.960334 ,0.967629 ,0.982935 ,0.993162 ,1.00241 ,1.00012 ,0.9814 ,0.937632 ,0.927918 ,0.93491 ,0.934184 ,0.921438 ,0.933788 ,0.954208 ,0.973684 ,0.977765 ,1.00697 ,1.01779 ,1.02017 ,1.01604 ,1.02121 ,1.02344 ,1.02888 ,1.02278 ,1.03096 ,1.02409 ,1.03149 ,1.03797 ,1.04254 ,1.04943 ,1.02433 ,1.02379 ,1.02177 ,1.02271 ,1.02011 ,1.02414 ,1.03649 ,1.03739 ,1.03739 ,1.03739};
double ElossCaloScale[30] = { 1.03504 ,1.03504 ,1.03504 ,1.02631 ,1.04258 ,1.03506 ,1.05307 ,1.02399 ,1.04237 ,1.04368 ,1.04043 ,1.05899 ,1.07933 ,1.08604 ,1.08984 ,1.03564 ,1.04158 ,1.05983 ,1.06291 ,1.06853 ,1.0674 ,1.05427 ,1.06466 ,1.06274 ,1.06141 ,1.06314 ,1.06868 ,1.07242 ,1.07242 ,1.07242};
double X0MuonScale[60] = {-0.0320612 ,-0.0320612 ,-0.0320612 ,-0.0320612 ,-0.0693796 ,-0.0389677 ,-0.0860891 ,-0.124606 ,-0.0882329 ,-0.100014 ,-0.0790912 ,-0.0745538 ,-0.099088 ,-0.0933711 ,-0.0618782 ,-0.0619762 ,-0.0658361 ,-0.109704 ,-0.129547 ,-0.143364 ,-0.0774768 ,-0.0739859 ,-0.0417835 ,-0.022119 ,0.00308797 ,0.0197657 ,-0.0137871 ,-0.036848 ,-0.0643794 ,-0.0514949 ,-0.0317105 ,0.016539 ,0.0308435 ,-0.00056883 ,-0.00756813 ,-0.00760612 ,-0.0234571 ,-0.0980915 ,-0.101175 ,-0.102354 ,-0.0920337 ,-0.100337 ,-0.0887628 ,0.0660931 ,0.228999 ,0.260675 ,0.266301 ,0.267907 ,0.281668 ,0.194433 ,0.132954 ,0.20707 ,0.220466 ,0.20936 ,0.191441 ,0.191441 ,0.191441 ,0.191441 ,0.191441 ,0.191441};
int i60 = fabs(eta)*20.;
if(i60<0) i60 = 0;
if(i60>59) i60 = 59;
if(icalo==1) {
//
// Girder parametrization
//
double x = phi+ 3.1416-3.1416/32.*int((3.1416+phi)/(3.1416/32.));
double scale = 0.;
double pi = acos(-1.);
if(x>pi/64.) x = pi/32.-x;
if(x<0.005) {
scale = X0CaloGirder[0]*(1-x/0.005)+X0CaloGirder[1] + X0CaloGirder[3];
} else if(x<0.017) {
scale = X0CaloGirder[1] + X0CaloGirder[3];
} else if(x<0.028) {
scale = X0CaloGirder[2] + X0CaloGirder[3];
} else {
scale = X0CaloGirder[3];
}
if(fabs(eta)>1.3) scale = 1.;
//
// eta dependence of X0
//
scale *= X0CaloScale[i60];
X0Scale = scale;
//
// eta dependence of Eloss
//
int i30 = fabs(eta)*10.;
if(i30<0) i30 = 0;
if(i30>29) i30 = 29;
return new EnergyLoss(deltaE, sigmaDeltaE, sigmaPlusDeltaE, sigmaMinusDeltaE);
double nfactor = 0.987363/1.05471;
ElossScale = nfactor*ElossCaloScale[i30]*scale;
} else {
//
// Muon system
//
// eta dependence of X0
//
double scale = 1. + X0MuonScale[i60];
//
// Muon scale is now 1 with MuonTrackingGeometry and TrkDetDescrGeoModelCnv fixes
//
scale = 1.0;
X0Scale = scale;
ElossScale = 0.93*scale;
}
}
......@@ -19,7 +19,6 @@
#include "TrkSurfaces/StraightLineSurface.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkTrack/Track.h"
#include "TrkGeometry/EntryLayerProvider.h"
#include "TrkGeometry/DetachedTrackingVolume.h"
#include "TrkGeometry/AlignableTrackingVolume.h"
#include "TrkGeometry/Layer.h"
......@@ -49,6 +48,8 @@
// Trk
#include "TrkSurfaces/PlaneSurface.h"
#include <memory>
// screen output measures
// "[+] Text describing layer - with " << layerRZoutput()
// "[+] Text describing position - at " << positionOutput()
......@@ -57,7 +58,6 @@
// reference surface for Blind extrapolation
//Trk::PlaneSurface Trk::Extrapolator::s_referenceSurface(new Amg::Transform3D(Trk::s_idTransform), 0.,0.);
double Trk::Extrapolator::s_distIncreaseTolerance = 100. * Gaudi::Units::millimeter;
double Trk::Extrapolator::s_distEntryLayerMax = 100. * Gaudi::Units::millimeter;
// constructor
Trk::Extrapolator::Extrapolator(const std::string& t, const std::string& n, const IInterface* p) :
AthAlgTool(t,n,p),
......@@ -628,7 +628,11 @@ const Trk::TrackParameters* Trk::Extrapolator::extrapolate(const IPropagator& p
}
// start from the nextParameter (which are at volume boundary)
if (nextParameters) {
resultParameters = extrapolateWithinDetachedVolumes(*currentPropagator,
if (!m_stepPropagator) {
ATH_MSG_ERROR("extrapolation in Calo/MS called without configured STEP propagator, aborting");
return 0;
}
resultParameters = extrapolateWithinDetachedVolumes(*m_stepPropagator,
*nextParameters,
sf,
*nextVolume,
......@@ -1353,14 +1357,10 @@ const Trk::TrackParameters* Trk::Extrapolator::extrapolateToNextMaterialLayer(co
// reset remaining counters
m_currentDense = m_dense ? m_currentStatic : m_highestVolume;
m_navigBoundaries.clear();
if (m_denseVols.size()>m_denseResolved.first) {
m_denseVols.resize(m_denseResolved.first);
m_denseBoundaries.resize(m_denseResolved.second);
}
if (m_layers.size()>m_layerResolved) {
m_layers.resize(m_layerResolved);
m_navigLays.resize(m_layerResolved);
}
if (m_denseVols.size()>m_denseResolved.first) m_denseVols.resize(m_denseResolved.first);
while (m_denseBoundaries.size()>m_denseResolved.second) m_denseBoundaries.pop_back();
if (m_layers.size()>m_layerResolved) m_navigLays.resize(m_layerResolved);
while (m_layers.size()>m_layerResolved) m_layers.pop_back();
// current detached volumes
// collect : subvolume boundaries, ordered/unordered layers, confined dense volumes
......@@ -1613,6 +1613,78 @@ const Trk::TrackParameters* Trk::Extrapolator::extrapolateToNextMaterialLayer(co
unsigned int iSol = 0;
while ( iSol < solutions.size() ) {
if ( solutions[iSol] < iDest + m_staticBoundaries.size() ) {
// material attached ?
const Trk::Layer* mb = m_navigSurfs[solutions[iSol]].first->materialLayer();
if (mb) {
if (mb->layerMaterialProperties() && mb->layerMaterialProperties()->fullMaterial(nextPar->position()) ) {
const IMaterialEffectsUpdator* currentUpdator = subMaterialEffectsUpdator(*m_currentStatic);
if ( currentUpdator) {
const Trk::TrackParameters* upNext = currentUpdator->update(nextPar, *mb, dir, particle, matupmode);
if (upNext && upNext!= nextPar ) throwIntoGarbageBin(upNext);
nextPar = upNext;
}
if (!nextPar) {
ATH_MSG_VERBOSE( " [+] Update may have killed track - return." );
m_parametersAtBoundary.resetBoundaryInformation();
return returnParameters;
}
// collect material
const Trk::MaterialProperties* lmat = mb->fullUpdateMaterialProperties(*nextPar);
double lx0 = lmat->x0();
double layThick = mb->thickness();
double thick = 0.;
double costr = fabs(nextPar->momentum().normalized().dot(mb->surfaceRepresentation().normal())) ;
if ( mb->surfaceRepresentation().isOnSurface(mb->surfaceRepresentation().center(),false,0.,0.) )
thick = fmin(mb->surfaceRepresentation().bounds().r(),
layThick/fabs(nextPar->momentum().normalized().dot(mb->surfaceRepresentation().normal())) );
else {
//const Trk::CylinderBounds* cyl = dynamic_cast<const Trk::CylinderBounds*> (&(nextLayer->surfaceRepresentation().bounds()));
//double hmax = cyl ? cyl->halflengthZ() : nextLayer->surfaceRepresentation().bounds().r();
thick = fmin(2*mb->thickness(), layThick/(1-costr));
}
if (!m_matstates && m_extrapolationCache) {
if(checkCache(" extrapolateToNextMaterialLayer thin")) {
double dInX0 = thick/lx0;
if(m_dumpCache) dumpCache(" extrapolateToNextMaterialLayer thin ");
m_extrapolationCache->updateX0(dInX0);
double currentqoverp=nextPar->parameters()[Trk::qOverP];
EnergyLoss* eloss = m_elossupdators[0]->energyLoss(*lmat,fabs(1./currentqoverp),1./costr,dir,particle);
m_extrapolationCache->updateEloss(eloss->meanIoni(),eloss->sigmaIoni(),eloss->meanRad(),eloss->sigmaRad());
if(m_dumpCache) dumpCache(" After");
delete eloss;
}
}
if (m_matstates) {
double dInX0 = thick/lx0;
double scatsigma=sqrt(m_msupdators[0]->sigmaSquare(*lmat,1./fabs(nextPar->parameters()[qOverP]),1.,particle));
Trk::ScatteringAngles *newsa=new Trk::ScatteringAngles(0,0,scatsigma/sin(nextPar->parameters()[Trk::theta]),scatsigma);
//energy loss
double currentqoverp=nextPar->parameters()[Trk::qOverP];
EnergyLoss* eloss = m_elossupdators[0]->energyLoss(*lmat,fabs(1./currentqoverp),1./costr,dir,particle);
// use curvilinear TPs to simplify retrieval by fitters
Trk::CurvilinearParameters* cvlTP = new Trk::CurvilinearParameters(nextPar->position(),nextPar->momentum(),nextPar->charge());
Trk::MaterialEffectsOnTrack* mefot = new Trk::MaterialEffectsOnTrack(dInX0,newsa,eloss,cvlTP->associatedSurface());
if(m_extrapolationCache) {
if(checkCache(" mat states extrapolateToNextMaterialLayer thin" )) {
if(m_dumpCache) dumpCache(" extrapolateToNextMaterialLayer thin");
m_extrapolationCache->updateX0(dInX0);
m_extrapolationCache->updateEloss(eloss->meanIoni(),eloss->sigmaIoni(),eloss->meanRad(),eloss->sigmaRad());
if(m_dumpCache) dumpCache(" After");
}
}
m_matstates->push_back(new TrackStateOnSurface(0,cvlTP,0,mefot));
}
}
} // end material update at massive (static volume) boundary
// static volume boundary; return to the main loop
unsigned int index = solutions[iSol]-iDest;
// use global coordinates to retrieve attached volume (just for static!)
......@@ -1646,16 +1718,11 @@ const Trk::TrackParameters* Trk::Extrapolator::extrapolateToNextMaterialLayer(co
// material update HERE and NOW (pre/post udpdate ? )
// don't repeat if identical to last update && input parameters on the layer
bool collect = true;
if ( nextLayer == m_lastMaterialLayer ) {
// cylinder layers !
const Trk::CylinderLayer* cyl = dynamic_cast<const Trk::CylinderLayer*> (nextLayer);
const Trk::SubtractedCylinderLayer* scyl = dynamic_cast<const Trk::SubtractedCylinderLayer*> (nextLayer);
if (!cyl && !scyl) {
if ( nextLayer == m_lastMaterialLayer && nextLayer->surfaceRepresentation().type()!=Trk::Surface::Cylinder ) {
ATH_MSG_DEBUG( " [!] This layer is identical to the one with last material update, return layer without repeating the update" );
collect = false;
if (!destSurf && (nextLayer->layerType()>0 || m_returnPassiveLayers) ) return nextPar->clone();
}
}
double layThick = nextLayer->thickness();
if (collect && layThick>0.) { // collect material
// get the right updator
......@@ -1975,6 +2042,7 @@ const Trk::TrackParameters* Trk::Extrapolator::extrapolateInAlignableTV(const IP
unsigned int iSol = 0;
while ( iSol < solutions.size() ) {
if ( solutions[iSol] < iDest + m_staticBoundaries.size() ) {
// TODO if massive boundary coded, add the material effects here
// static volume boundary; return to the main loop : TODO move from misaligned to static
unsigned int index = solutions[iSol]-iDest;
// use global coordinates to retrieve attached volume (just for static!)
......@@ -2090,7 +2158,7 @@ std::pair<const Trk::TrackParameters*,const Trk::Layer*> Trk::Extrapolator::extr
for (; dIter != detVols->end(); dIter++) {
const Trk::Layer* lay = (*dIter)->layerRepresentation();
if ( lay ) {
bool checkBounds = lay->layerType() > 0 ? bcheck : true;
Trk::BoundaryCheck checkBounds = lay->layerType() > 0 ? bcheck : Trk::BoundaryCheck(true);
std::pair<const Trk::Surface*, Trk::BoundaryCheck> newSurf(&(lay->surfaceRepresentation()),checkBounds);
m_navigSurfs.push_back(newSurf);
m_navigVols.push_back(*dIter);
......@@ -2810,12 +2878,13 @@ const Trk::TrackParameters* Trk::Extrapolator::insideVolumeStaticLayers(
// ============================ RESOLVE STARTPOINT =============================
// only if you do not have an input associated Layer
// - this means that a volume step has been done
if (!associatedLayer){
ATH_MSG_VERBOSE( " [+] Volume switch has happened, searching for entry layers." );
// reset the exitFace
exitFace = m_parametersAtBoundary.exitFace;
// Step [1] Check for entry layers -------------------------------------------------
associatedLayer = tvol.entryLayer(navParameters->position(),dir*navParameters->momentum().normalized());
associatedLayer = tvol.associatedLayer(navParameters->position());
if ( associatedLayer && associatedLayer->layerMaterialProperties() ){
//-------------------------------------------------------------------------------
ATH_MSG_VERBOSE( " [+] Entry layer to volume found with " << layerRZoutput(*associatedLayer) );
......@@ -2926,48 +2995,13 @@ const Trk::TrackParameters* Trk::Extrapolator::insideVolumeStaticLayers(
return nextParameters;
}
// ----------------------------------------------------------------------------------------------------------
// Case II: Volume EXIT
// same game as before: find the best parameters for the navigation search
// - if any action has been taken in this volume (i.e. nextParameters != pars), take nextParameters
// - if no action has been taken, use the navigation parameters since they give a better estimate
navParameters = (nextParameters == (&parm) && m_parametersAtBoundary.navParameters) ?
m_parametersAtBoundary.navParameters : nextParameters;
if (navParameters && tvol.entryLayerProvider()){
// call for the exit layer search