diff --git a/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/TrigFastTrackFinder_Config.py b/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/TrigFastTrackFinder_Config.py
index ccd1794f763d30de8da38944f0733e487fc42c59..3ffb496fd7749476b39b5940ffdedd672746065f 100755
--- a/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/TrigFastTrackFinder_Config.py
+++ b/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/TrigFastTrackFinder_Config.py
@@ -207,12 +207,6 @@ class TrigFastTrackFinderBase(TrigFastTrackFinder):
from TrigInDetConf.TrigInDetRecCommonTools import InDetTrigFastTrackSummaryTool
self.TrackSummaryTool = InDetTrigFastTrackSummaryTool
- if self.doResMon:
- from TrigInDetTrackFitter.TrigInDetTrackFitterConf import TrigL2ResidualCalculator
- resCalc = TrigL2ResidualCalculator(OfflineClusters=False)
- ToolSvc += resCalc
- self.TrigL2ResidualCalculator = resCalc
-
# Run3 monitoring
self.MonTool = TrigFastTrackFinderMonitoring(type, self.doResMon)
diff --git a/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.cxx b/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.cxx
index 38c0c5cd2b07b8082c4a54f733d0e24a217b979d..55dd30897b76ee2c3e2ba0ef878f9453e7b21f8b 100644
--- a/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.cxx
+++ b/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.cxx
@@ -54,7 +54,6 @@
#include "TrigInDetToolInterfaces/ITrigL2LayerNumberTool.h"
#include "TrigInDetToolInterfaces/ITrigSpacePointConversionTool.h"
#include "TrigInDetToolInterfaces/ITrigL2SpacePointTruthTool.h"
-#include "TrigInDetToolInterfaces/ITrigL2ResidualCalculator.h"
#include "TrigInDetToolInterfaces/TrigL2HitResidual.h"
#include "TrigInDetToolInterfaces/ITrigInDetTrackFitter.h"
@@ -86,7 +85,6 @@ TrigFastTrackFinder::TrigFastTrackFinder(const std::string& name, ISvcLocator* p
HLT::FexAlgo(name, pSvcLocator),
m_numberingTool("TrigL2LayerNumberTool"),
m_spacePointTool("TrigSpacePointConversionTool"),
- m_trigL2ResidualCalculator("TrigL2ResidualCalculator"),
m_trackMaker("InDet::SiTrackMaker_xk/InDetTrigSiTrackMaker"),
m_trigInDetTrackFitter("TrigInDetTrackFitter"),
m_trigZFinder("TrigZFinder/TrigZFinder", this ),
@@ -161,7 +159,6 @@ TrigFastTrackFinder::TrigFastTrackFinder(const std::string& name, ISvcLocator* p
declareProperty( "initialTrackMaker", m_trackMaker);
declareProperty( "trigInDetTrackFitter", m_trigInDetTrackFitter );
declareProperty( "trigZFinder", m_trigZFinder );
- declareProperty( "TrigL2ResidualCalculator", m_trigL2ResidualCalculator );
declareProperty("TrackSummaryTool", m_trackSummaryTool);
declareProperty( "doResMon", m_doResMonitoring = true);
@@ -248,16 +245,6 @@ HLT::ErrorCode TrigFastTrackFinder::hltInitialize() {
return HLT::BAD_JOB_SETUP;
}
- if (m_doResMonitoring) {
- sc = m_trigL2ResidualCalculator.retrieve();
- if ( sc.isFailure() ) {
- msg() << MSG::FATAL <<"Unable to locate Residual calculator tool " << m_trigL2ResidualCalculator << endmsg;
- return HLT::BAD_JOB_SETUP;
- }
- } else {
- m_trigL2ResidualCalculator.disable();
- }
-
//Get ID helper
if (detStore()->retrieve(m_idHelper, "AtlasID").isFailure()) {
ATH_MSG_ERROR("Could not get AtlasDetectorID helper AtlasID");
@@ -671,7 +658,7 @@ StatusCode TrigFastTrackFinder::findTracks(InDet::SiTrackMakerEventData_xk &trac
mnt_timer_TrackFitter.start(); // Run3 monitoring
- m_trigInDetTrackFitter->fit(initialTracks, outputTracks, m_particleHypothesis);
+ m_trigInDetTrackFitter->fit(initialTracks, outputTracks, ctx, m_particleHypothesis);
if( outputTracks.empty() ) {
ATH_MSG_DEBUG("REGTEST / No tracks fitted");
@@ -1018,12 +1005,12 @@ void TrigFastTrackFinder::fillMon(const TrackCollection& tracks, const TrigVerte
// tighter selection for unbiased residuals
bool goodTrack = std::fabs(pT)>1000. && (nPix + nSct/2) > 3 && nSct > 0;
if (goodTrack && m_doResMonitoring) {
- runResidualMonitoring(*track);
+ runResidualMonitoring(*track, ctx);
}
}
}
-void TrigFastTrackFinder::runResidualMonitoring(const Trk::Track& track) const {
+void TrigFastTrackFinder::runResidualMonitoring(const Trk::Track& track, const EventContext& ctx) const {
// Run3 monitoring ---------->
std::vector<float> mnt_layer_IBL;
@@ -1124,7 +1111,7 @@ void TrigFastTrackFinder::runResidualMonitoring(const Trk::Track& track) const {
std::vector<TrigL2HitResidual> vResid;
vResid.clear();
- StatusCode scRes = m_trigL2ResidualCalculator->getUnbiasedResiduals(track,vResid);
+ StatusCode scRes = m_trigInDetTrackFitter->getUnbiasedResiduals(track,vResid, ctx);
if(!scRes.isSuccess()) return;
for(std::vector<TrigL2HitResidual>::iterator it=vResid.begin();it!=vResid.end();++it) {
Identifier id = it->identify();
diff --git a/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.h b/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.h
index b4e54a9b26e0bec3e6547918436040d82be68b57..b00c6cf14acb4bea9d90bdbb860dd3db39fe88f4 100644
--- a/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.h
+++ b/Trigger/TrigAlgorithms/TrigFastTrackFinder/src/TrigFastTrackFinder.h
@@ -35,7 +35,6 @@
class ITrigL2LayerNumberTool;
class ITrigSpacePointConversionTool;
-class ITrigL2ResidualCalculator;
class ITrigInDetTrackFitter;
class ITrigZFinder;
class TrigRoiDescriptor;
@@ -100,7 +99,6 @@ protected:
ToolHandle<ITrigL2LayerNumberTool> m_numberingTool;
ToolHandle<ITrigSpacePointConversionTool> m_spacePointTool;
- ToolHandle<ITrigL2ResidualCalculator> m_trigL2ResidualCalculator;
ToolHandle<InDet::ISiTrackMaker> m_trackMaker; // Track maker
ToolHandle<ITrigInDetTrackFitter> m_trigInDetTrackFitter;
ToolHandle<ITrigZFinder> m_trigZFinder;
@@ -153,7 +151,7 @@ protected:
// Monitoring member functions
void fillMon(const TrackCollection& tracks, const TrigVertexCollection& vertices, const TrigRoiDescriptor& roi, const EventContext& ctx) const;
- void runResidualMonitoring(const Trk::Track& track) const;
+ void runResidualMonitoring(const Trk::Track& track, const EventContext&) const;
//Setup functions
void getBeamSpot(float&, float&, const EventContext&) const;
diff --git a/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigInDetTrackFitter.h b/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigInDetTrackFitter.h
index e6262dbe1b1ac9e732f063491fe9d5fde70d6895..4127c2436f8173583a98b01a7ae62ead7c6730f6 100755
--- a/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigInDetTrackFitter.h
+++ b/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigInDetTrackFitter.h
@@ -7,9 +7,14 @@
#include <vector>
#include "GaudiKernel/IAlgTool.h"
-#include "TrigInDetEvent/TrigInDetTrackCollection.h"
#include "TrkTrack/TrackCollection.h"
+namespace Trk {
+ class Track;
+}
+
+class TrigL2HitResidual;
+
static const InterfaceID IID_ITrigInDetTrackFitter("ITrigInDetTrackFitter", 1 , 0);
class ITrigInDetTrackFitter: virtual public IAlgTool
@@ -19,7 +24,8 @@ class ITrigInDetTrackFitter: virtual public IAlgTool
static const InterfaceID& interfaceID() {
return IID_ITrigInDetTrackFitter;
}
- virtual void fit(const TrackCollection&, TrackCollection&, const Trk::ParticleHypothesis& matEffects) const = 0;
+ virtual void fit(const TrackCollection&, TrackCollection&, const EventContext&, const Trk::ParticleHypothesis& matEffects) const = 0;
+ virtual StatusCode getUnbiasedResiduals(const Trk::Track&, std::vector<TrigL2HitResidual>&, const EventContext&) const = 0;
};
#endif
diff --git a/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigL2ResidualCalculator.h b/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigL2ResidualCalculator.h
deleted file mode 100644
index c6b9f0e10e7a9a386179bdab1af65be98870dfae..0000000000000000000000000000000000000000
--- a/Trigger/TrigTools/TrigInDetToolInterfaces/TrigInDetToolInterfaces/ITrigL2ResidualCalculator.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef __ITRIGL2RESIDUALCALCULATOR_H__
-#define __ITRIGL2RESIDUALCALCULATOR_H__
-
-#include "GaudiKernel/IAlgTool.h"
-#include <vector>
-
-namespace Trk {
- class Track;
-}
-class TrigL2HitResidual;
-
-static const InterfaceID IID_ITrigL2ResidualCalculator("ITrigL2ResidualCalculator",1,0);
-
-
- /** @class ITrigL2ResidualCalculator
-
- provides the abstract interface for the LVL2 ID tool that calculates smoothed residuals
- of Pixel and SCT hits
-
- @author D.Emeliyanov <http://consult.cern.ch/xwho>
- */
-
- class ITrigL2ResidualCalculator : virtual public IAlgTool {
-
- public:
- /** other standard AlgTool methods */
-
- static const InterfaceID& interfaceID () //!< the Tool's interface
- { return IID_ITrigL2ResidualCalculator; }
-
- virtual StatusCode getUnbiasedResiduals(const Trk::Track&, std::vector<TrigL2HitResidual>&) const = 0;
- };
-
-#endif
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.cxx b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.cxx
index 64b117ebe1e288fe5bf68b3d602b6ca419d0ebe7..17d64310126298f6a62f8d4664767d75c1b336fd 100644
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.cxx
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.cxx
@@ -41,6 +41,10 @@
#include "AthenaBaseComps/AthMsgStreamMacros.h"
#include "AthenaBaseComps/AthCheckMacros.h"
+#include "AtlasDetDescr/AtlasDetectorID.h"
+#include "InDetIdentifier/SCT_ID.h"
+#include "InDetIdentifier/PixelID.h"
+#include "TrigInDetToolInterfaces/TrigL2HitResidual.h"
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "TrkPrepRawData/PrepRawData.h"
@@ -50,20 +54,18 @@
#include "TrkEventPrimitives/FitQuality.h"
#include "TrkEventPrimitives/FitQualityOnSurface.h"
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
TrigInDetTrackFitter::TrigInDetTrackFitter(const std::string& t,
const std::string& n,
const IInterface* p ): AthAlgTool(t,n,p),
- m_MagFieldSvc("AtlasFieldSvc",this->name()),
m_trackMaker("TrigDkfTrackMakerTool")
{
declareInterface< ITrigInDetTrackFitter >( this );
- declareProperty( "AthenaFieldService", m_MagFieldSvc, "AlasFieldService");
declareProperty( "doMultScattering", m_doMultScatt = true);
declareProperty( "doBremmCorrection", m_doBremm=false);
declareProperty( "Chi2Cut", m_DChi2 = 1000.0);
- declareProperty( "OfflineClusters", m_offlineClusters = true);
declareProperty( "correctClusterPos", m_correctClusterPos = false);
declareProperty( "ROTcreator", m_ROTcreator, "ROTcreatorTool" );
m_nTracksTotal = 0;
@@ -74,12 +76,17 @@ TrigInDetTrackFitter::TrigInDetTrackFitter(const std::string& t,
StatusCode TrigInDetTrackFitter::initialize()
{
- ATH_CHECK(m_MagFieldSvc.retrieve());
ATH_CHECK(m_trackMaker.retrieve());
if (m_correctClusterPos) {
ATH_CHECK(m_ROTcreator.retrieve());
}
+ ATH_CHECK( m_fieldCondObjInputKey.initialize());
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "AtlasID"));
+ ATH_CHECK(detStore()->retrieve(m_pixelId, "PixelID"));
+ ATH_CHECK(detStore()->retrieve(m_sctId, "SCT_ID"));
+
return StatusCode::SUCCESS;
+
}
StatusCode TrigInDetTrackFitter::finalize()
@@ -95,10 +102,10 @@ StatusCode TrigInDetTrackFitter::finalize()
return StatusCode::SUCCESS;
}
-void TrigInDetTrackFitter::getMagneticField(double r[3],double* B) const {
+void TrigInDetTrackFitter::getMagneticField(double r[3],double* B, MagField::AtlasFieldCache& fieldCache) const {
B[0]=0.0;B[1]=0.0;B[2]=0.0;
double field[3];
- m_MagFieldSvc->getField(r,field);//field is returned in kT
+ fieldCache.getField(r,field);//field is returned in kT
for(int i=0;i<3;i++) B[i]=field[i]/Gaudi::Units::kilogauss;//convert to kG
}
@@ -126,7 +133,8 @@ void TrigInDetTrackFitter::correctScale(Trk::TrkTrackState* pTS) const {
Trk::TrkTrackState* TrigInDetTrackFitter::extrapolate(Trk::TrkTrackState* pTS,
Trk::TrkPlanarSurface* pSB,
- Trk::TrkPlanarSurface* pSE) const
+ Trk::TrkPlanarSurface* pSE,
+ MagField::AtlasFieldCache& fieldCache) const
{
const double C=0.02999975/1000.0;//using GeV internally
const double minStep=30.0;
@@ -154,7 +162,6 @@ Trk::TrkTrackState* TrigInDetTrackFitter::extrapolate(Trk::TrkTrackState* pTS,
int nStep,nStepMax;
double sl,ds,path=0.0;
- //m_numericalJacobian(pTS,pSB,pSE,J);
double sint,cost,sinf,cosf;
sint=sin(pTS->getTrackState(3));cosf=cos(pTS->getTrackState(2));
sinf=sin(pTS->getTrackState(2));cost=cos(pTS->getTrackState(3));
@@ -193,7 +200,7 @@ Trk::TrkTrackState* TrigInDetTrackFitter::extrapolate(Trk::TrkTrackState* pTS,
for(i=0;i<4;i++) D[i]=pSE->getPar(i);
for(i=0;i<3;i++) gPi[i]=gP[i];
- getMagneticField(gP,gB);
+ getMagneticField(gP,gB, fieldCache);
for(i=0;i<3;i++) gBi[i]=gB[i];
@@ -248,7 +255,7 @@ Trk::TrkTrackState* TrigInDetTrackFitter::extrapolate(Trk::TrkTrackState* pTS,
V[1]=gV[1]+Av*DVy;
V[2]=gV[2]+Av*DVz;
- getMagneticField(P,gB);
+ getMagneticField(P,gB,fieldCache);
for(i=0;i<3;i++) gBf[i]=gB[i];
for(i=0;i<3;i++)
@@ -551,18 +558,24 @@ Trk::TrkTrackState* TrigInDetTrackFitter::extrapolate(Trk::TrkTrackState* pTS,
return pTE;
}
-void TrigInDetTrackFitter::fit(const TrackCollection& inputTracks, TrackCollection& fittedTracks, const Trk::ParticleHypothesis& matEffects) const
+void TrigInDetTrackFitter::fit(const TrackCollection& inputTracks, TrackCollection& fittedTracks, const EventContext& ctx, const Trk::ParticleHypothesis& matEffects) const
{
+ MagField::AtlasFieldCache fieldCache;
+
+ SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
+ const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
+
+ fieldCondObj->getInitializedCache (fieldCache);
fittedTracks.reserve(inputTracks.size());
for(auto trIt = inputTracks.begin(); trIt != inputTracks.end(); ++trIt) {
- Trk::Track* fittedTrack = fitTrack(**trIt, matEffects);
+ Trk::Track* fittedTrack = fitTrack(**trIt, fieldCache, matEffects);
if (fittedTrack!=nullptr) {
fittedTracks.push_back(fittedTrack);
}
}
}
-Trk::Track* TrigInDetTrackFitter::fitTrack(const Trk::Track& recoTrack, const Trk::ParticleHypothesis& matEffects) const {
+Trk::Track* TrigInDetTrackFitter::fitTrack(const Trk::Track& recoTrack, MagField::AtlasFieldCache& fieldCache, const Trk::ParticleHypothesis& matEffects) const {
const Trk::TrackParameters* trackPars = recoTrack.perigeeParameters();
if(trackPars==nullptr) {
@@ -628,7 +641,7 @@ Trk::Track* TrigInDetTrackFitter::fitTrack(const Trk::Track& recoTrack, const Tr
Trk::TrkPlanarSurface* pSE=nullptr;
for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
pSE=(*pnIt)->getSurface();
- Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE);
+ Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
pSB=pSE;
if(pNS!=nullptr) {
@@ -823,3 +836,192 @@ Trk::TrackStateOnSurface* TrigInDetTrackFitter::createTrackStateOnSurface(Trk::T
pTSS = new Trk::TrackStateOnSurface(pRIO, 0, pFQ, 0, typePattern);
return pTSS;
}
+
+StatusCode TrigInDetTrackFitter::getUnbiasedResiduals(const Trk::Track& pT,
+ std::vector<TrigL2HitResidual>& vResid, const EventContext& ctx) const {
+
+ const Trk::TrackParameters* trackPars = pT.perigeeParameters();
+ if(trackPars==nullptr) {
+ ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
+ return StatusCode::FAILURE;
+ }
+
+ MagField::AtlasFieldCache fieldCache;
+
+ SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
+ const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
+
+ fieldCondObj->getInitializedCache (fieldCache);
+ std::vector<Trk::TrkBaseNode*> vpTrkNodes;
+ std::vector<Trk::TrkTrackState*> vpTrackStates;
+ vResid.clear();
+ double trk_theta = trackPars->parameters()[Trk::theta];
+ double trk_qOverP = trackPars->parameters()[Trk::qOverP];
+ double Pt = sin(trk_theta)/trk_qOverP;
+ if(fabs(Pt)<100.0)
+ {
+ ATH_MSG_DEBUG("TrigL2ResidualCalculator failed -- Estimated Pt is too low "<<Pt);
+ return StatusCode::FAILURE;
+ }
+
+ // 1. Create filtering nodes
+
+ bool trackResult = m_trackMaker->createDkfTrack(pT,vpTrkNodes, m_DChi2);
+ if(!trackResult)
+ {
+ ATH_MSG_DEBUG("TrigDkfTrackMaker failed");
+ return StatusCode::FAILURE;
+ }
+
+ int nNodeIndex=-1;
+ bool OK=true;
+ std::vector<Trk::TrkBaseNode*>::iterator pnIt,pnEnd(vpTrkNodes.end());
+
+ for(std::vector<Trk::TrkBaseNode*>::iterator pNodeIt=vpTrkNodes.begin();pNodeIt!=vpTrkNodes.end();
+ ++pNodeIt)
+ {
+ nNodeIndex++;
+ Trk::TrkBaseNode* pMaskedNode=(*pNodeIt);
+ Trk::TrkTrackState* pMaskedState=NULL;
+
+ // 2. Create initial track state:
+
+ double Rk[5];
+ Rk[0] = trackPars->parameters()[Trk::d0];
+ Rk[1] = trackPars->parameters()[Trk::z0];
+ Rk[2] = trackPars->parameters()[Trk::phi0];
+ if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
+ if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
+ trk_theta = trackPars->parameters()[Trk::theta];
+ Rk[3] = trk_theta;
+ Rk[4] = 1000.0*trackPars->parameters()[Trk::qOverP];//MeV->GeV
+ //No need to correct scale back - not returning track
+
+ // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
+
+ Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
+ double Gk[5][5] = {{100.0, 0, 0, 0, 0},
+ {0, 100.0, 0, 0, 0},
+ {0, 0, 0.01, 0, 0},
+ {0, 0, 0, 0.01, 0},
+ {0, 0, 0, 0, 0.1}};
+ pTS->setTrackCovariance(Gk);
+ if(m_doMultScatt)
+ pTS->setScatteringMode(1);
+ if(m_doBremm)
+ pTS->setScatteringMode(2);
+ vpTrackStates.push_back(pTS);
+
+ ATH_MSG_DEBUG("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
+ <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]);
+
+ OK=true;
+ Trk::TrkPlanarSurface *pSB=NULL,*pSE;
+
+ for(pnIt=vpTrkNodes.begin();pnIt!=pnEnd;++pnIt)
+ {
+ pSE=(*pnIt)->getSurface();
+ Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
+
+ pSB=pSE;
+ if(pNS!=NULL)
+ {
+ vpTrackStates.push_back(pNS);
+
+ (*pnIt)->validateMeasurement(pNS);
+ ATH_MSG_DEBUG("dChi2="<<(*pnIt)->getChi2());
+ if((*pnIt)!=pMaskedNode)
+ {
+ (*pnIt)->updateTrackState(pNS);
+ }
+ else
+ {
+ pMaskedState=pNS;
+ }
+ pTS=pNS;
+ Pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
+ if(fabs(Pt)<0.2)
+ {
+ ATH_MSG_DEBUG("Estimated Pt is too low "<<Pt<<" - skipping fit");
+ OK=false;break;
+ }
+ }
+ else
+ {
+ OK=false;break;
+ }
+ }
+ if(OK)
+ {
+ std::vector<Trk::TrkTrackState*>::reverse_iterator ptsrIt(vpTrackStates.rbegin()),
+ ptsrEnd(vpTrackStates.rend());
+
+ for(;ptsrIt!=ptsrEnd;++ptsrIt)
+ {
+ (*ptsrIt)->runSmoother();
+ }
+
+ pMaskedNode->validateMeasurement(pMaskedState);
+
+ double r[2],V[2][2];
+
+ int nSize=pMaskedNode->getResiduals(r);
+ nSize=pMaskedNode->getInverseResidualVariance(V);
+ const Trk::PrepRawData* pPRD = pMaskedNode->getPrepRawData();
+
+ Identifier id = pPRD->identify();
+
+ Region region = Region::Undefined;
+ if(m_idHelper->is_pixel(id))
+ {
+ region=(m_pixelId->is_barrel(id)) ? Region::PixBarrel: Region::PixEndcap;
+ if (m_pixelId->is_blayer(id)) {
+ region = Region::IBL;
+ }
+ }
+ if(m_idHelper->is_sct(id))
+ {
+ region=(m_sctId->is_barrel(id)) ? Region::SctBarrel : Region::SctEndcap;
+ }
+ if(nSize==1) {
+ if(V[0][0]>0.0) {
+ vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0])));
+ }
+ else {
+ OK=false;
+ break;
+ }
+ }
+ else {
+ if((V[0][0]>0.0) && (V[1][1]>0.0)) {
+ vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0]),
+ r[1],r[1]*sqrt(V[1][1])));
+ }
+ else {
+ OK=false;
+ break;
+ }
+ }
+ }
+ else
+ {
+ ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
+ vResid.clear();
+ }
+ for(std::vector<Trk::TrkTrackState*>::iterator ptsIt=vpTrackStates.begin();
+ ptsIt!=vpTrackStates.end();++ptsIt) delete (*ptsIt);
+ vpTrackStates.clear();
+ if(!OK) break;
+ }
+ pnIt=vpTrkNodes.begin();pnEnd=vpTrkNodes.end();
+ for(;pnIt!=pnEnd;++pnIt)
+ {
+ delete((*pnIt)->getSurface());
+ delete (*pnIt);
+ }
+ vpTrkNodes.clear();
+
+ if(OK) return StatusCode::SUCCESS;
+ else return StatusCode::FAILURE;
+
+}
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.h b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.h
index e24815fc4679ffe8daf251f3eeb400ba618bf8a4..31d9f0eebb625b33a39f33617f87c17453468698 100755
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.h
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFitter.h
@@ -9,14 +9,14 @@
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/ServiceHandle.h"
-#include "MagFieldInterfaces/IMagFieldSvc.h"
#include "TrigInDetToolInterfaces/ITrigInDetTrackFitter.h"
#include "TrigInDetToolInterfaces/ITrigDkfTrackMakerTool.h"
#include "TrkToolInterfaces/IRIO_OnTrackCreator.h"
+// MagField cache
+#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
+#include "MagFieldElements/AtlasFieldCache.h"
-class TrigTimer;
-
namespace Trk {
class TrkBaseNode;
class TrkTrackState;
@@ -25,9 +25,10 @@ namespace Trk {
class TrackStateOnSurface;
}
-namespace MagField {
- class IMagFieldSvc;
-}
+class AtlasDetectorID;
+class PixelID;
+class SCT_ID;
+class TrigL2HitResidual;
class TrigInDetTrackFitter: public AthAlgTool, virtual public ITrigInDetTrackFitter
{
@@ -35,14 +36,17 @@ class TrigInDetTrackFitter: public AthAlgTool, virtual public ITrigInDetTrackFit
TrigInDetTrackFitter( const std::string&, const std::string&, const IInterface* );
virtual StatusCode initialize();
virtual StatusCode finalize();
- Trk::Track* fitTrack(const Trk::Track&, const Trk::ParticleHypothesis& matEffects = Trk::pion) const;
- void fit(const TrackCollection&, TrackCollection&, const Trk::ParticleHypothesis& matEffects = Trk::pion) const;
+ Trk::Track* fitTrack(const Trk::Track&, MagField::AtlasFieldCache&, const Trk::ParticleHypothesis& matEffects = Trk::pion) const;
+ void fit(const TrackCollection&, TrackCollection&, const EventContext&, const Trk::ParticleHypothesis& matEffects = Trk::pion) const;
+ StatusCode getUnbiasedResiduals(const Trk::Track&, std::vector<TrigL2HitResidual>&, const EventContext&) const;
+
private:
Trk::TrkTrackState* extrapolate(Trk::TrkTrackState*,
Trk::TrkPlanarSurface*,
- Trk::TrkPlanarSurface*) const;
- void getMagneticField(double[3],double*) const;
+ Trk::TrkPlanarSurface*,
+ MagField::AtlasFieldCache&) const;
+ void getMagneticField(double[3],double*, MagField::AtlasFieldCache&) const;
void correctScale(Trk::TrkTrackState*) const;
@@ -57,11 +61,14 @@ private:
double m_DChi2;
bool m_doMultScatt;
bool m_doBremm;
- bool m_offlineClusters;
bool m_correctClusterPos;
- ServiceHandle<MagField::IMagFieldSvc> m_MagFieldSvc;
ToolHandle<ITrigDkfTrackMakerTool> m_trackMaker;
ToolHandle<Trk::IRIO_OnTrackCreator> m_ROTcreator;
+ SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
+ const PixelID* m_pixelId;
+ const SCT_ID* m_sctId;
+ const AtlasDetectorID* m_idHelper;
+
};
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigL2ResidualCalculator.cxx b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigL2ResidualCalculator.cxx
deleted file mode 100644
index 335c67e71b8d77f4839c4332d25f08e0cc5b9041..0000000000000000000000000000000000000000
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigL2ResidualCalculator.cxx
+++ /dev/null
@@ -1,709 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-////////////////////////////////////////////////////////////////////////////////
-// TrigL2ResidualCalculator tool
-// -------------------------------
-// ATLAS Collaboration
-//
-// 11.11.2009 Package created
-//
-// Author: Dmitry Emeliyanov, RAL
-// e-mail: D.Emeliyanov@rl.ac.uk
-//
-////////////////////////////////////////////////////////////////////////////////
-
-#include <cmath>
-#include <iostream>
-
-#include "GaudiKernel/SystemOfUnits.h"
-
-#include "AtlasDetDescr/AtlasDetectorID.h"
-#include "InDetIdentifier/SCT_ID.h"
-#include "InDetIdentifier/PixelID.h"
-
-#include "TrkTrack/Track.h"
-#include "TrkParameters/TrackParameters.h"
-
-//#include "TrkSurfaces/Surface.h"
-//#include "TrkSurfaces/TrapezoidBounds.h"
-//#include "TrkParameters/MeasuredPerigee.h"
-
-#include "TrkPrepRawData/PrepRawData.h"
-#include "TrkDistributedKalmanFilter/TrkBaseNode.h"
-#include "TrkDistributedKalmanFilter/TrkFilteringNodes.h"
-#include "TrkDistributedKalmanFilter/TrkTrackState.h"
-#include "TrkDistributedKalmanFilter/TrkPlanarSurface.h"
-#include "MagFieldInterfaces/IMagFieldSvc.h"
-
-//#include "TrigInDetToolInterfaces/ITrigL2ResidualCalculator.h"
-#include "TrigInDetToolInterfaces/TrigL2HitResidual.h"
-#include "TrigL2ResidualCalculator.h"
-
-TrigL2ResidualCalculator::TrigL2ResidualCalculator(const std::string& t,
- const std::string& n,
- const IInterface* p ) :
- AthAlgTool(t,n,p),
- m_MagFieldSvc("AtlasFieldSvc",this->name()),
- m_trackMaker("TrigDkfTrackMakerTool")
-{
- declareInterface< ITrigL2ResidualCalculator >( this );
-
- declareProperty( "AthenaFieldService", m_MagFieldSvc, "AtlasFieldSvc");
- declareProperty( "doMultScattering", m_doMultScatt = true);
- declareProperty( "doBremmCorrection", m_doBremm=false);
- declareProperty( "Chi2Cut", m_DChi2 = 1000.0);
- declareProperty( "OfflineClusters", m_offlineClusters = true);
-}
-
-StatusCode TrigL2ResidualCalculator::initialize()
-{
- StatusCode sc = AthAlgTool::initialize();
-
-
- if (detStore()->retrieve(m_idHelper, "AtlasID").isFailure()) {
- ATH_MSG_FATAL("Could not get AtlasDetectorID helper AtlasID");
- return StatusCode::FAILURE;
- }
-
- if (detStore()->retrieve(m_pixelId, "PixelID").isFailure()) {
- ATH_MSG_FATAL("Could not get Pixel ID helper");
- return StatusCode::FAILURE;
- }
- if (detStore()->retrieve(m_sctId, "SCT_ID").isFailure()) {
- ATH_MSG_FATAL("Could not get SCT ID helper");
- return StatusCode::FAILURE;
- }
- ATH_MSG_INFO("Using Athena magnetic field service");
- sc = m_MagFieldSvc.retrieve();
- if(sc.isFailure())
- {
- ATH_MSG_ERROR("Unable to retrieve Athena MagFieldService");
- return StatusCode::FAILURE;
- }
- sc=m_trackMaker.retrieve();
- if(sc.isFailure())
- {
- ATH_MSG_ERROR("Could not retrieve "<<m_trackMaker);
- return sc;
- }
- ATH_MSG_INFO("TrigL2ResidualCalculator constructed ");
- return sc;
-}
-
-
-void TrigL2ResidualCalculator::getMagneticField(double r[3],double* B) const
-{
- B[0]=0.0;B[1]=0.0;B[2]=0.0;
- double field[3];
- m_MagFieldSvc->getField(r,field);//field is returned in kT
- for(int i=0;i<3;i++) B[i]=field[i]/Gaudi::Units::kilogauss;//convert to kG
-}
-
-Trk::TrkTrackState* TrigL2ResidualCalculator::extrapolate(Trk::TrkTrackState* pTS,
- Trk::TrkPlanarSurface* pSB,
- Trk::TrkPlanarSurface* pSE) const
-{
- const double C=0.02999975/1000.0;//using GeV internally
- const double minStep=30.0;
-
- double J[5][5],Rf[5],AG[5][5],Gf[5][5],A[5][5];
- int i,j,m;
-
- bool samePlane=false;
-
- if(pSB!=NULL)
- {
- double diff=0.0;
- for(i=0;i<4;i++) diff+=fabs(pSE->getPar(i)-pSB->getPar(i));
- if(diff<1e-5) {
- samePlane=true;
- //std::cout<<"Starting plane and target plane are the same !"<<std::endl;
- }
- }
-
- if(!samePlane) {
-
- double gP[3],gPi[3],lP[3],gV[3],a,b,c,s,J0[7][5],descr,CQ,Ac,Av,Cc;
- double V[3],P[3],M[3][3],D[4],Jm[7][7],
- J1[5][7],gB[3],gBi[3],gBf[3],dBds[3],Buf[5][7],DVx,DVy,DVz;
- int nStep,nStepMax;
- double sl,ds,path=0.0;
-
- //m_numericalJacobian(pTS,pSB,pSE,J);
- double sint,cost,sinf,cosf;
- sint=sin(pTS->getTrackState(3));cosf=cos(pTS->getTrackState(2));
- sinf=sin(pTS->getTrackState(2));cost=cos(pTS->getTrackState(3));
- gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=cost;CQ=C*pTS->getTrackState(4);
-
- memset(&J0[0][0],0,sizeof(J0));
-
- if(pSB!=NULL)
- {
- double L[3][3];
- lP[0]=pTS->getTrackState(0);lP[1]=pTS->getTrackState(1);lP[2]=0.0;
- pSB->transformPointToGlobal(lP,gP);
- for(i=0;i<3;i++) for(j=0;j<3;j++) L[i][j]=pSB->getInvRotMatrix(i,j);
-
- J0[0][0]=L[0][0];J0[0][1]=L[0][1];
- J0[1][0]=L[1][0];J0[1][1]=L[1][1];
- J0[2][0]=L[2][0];J0[2][1]=L[2][1];
- J0[3][2]=-sinf*sint;J0[3][3]=cosf*cost;
- J0[4][2]= cosf*sint;J0[4][3]=sinf*cost;
- J0[5][3]=-sint;
- J0[6][4]=1.0;
- }
- else
- {
- gP[0]=-pTS->getTrackState(0)*sinf;
- gP[1]= pTS->getTrackState(0)*cosf;
- gP[2]= pTS->getTrackState(1);
- J0[0][0]=-sinf;J0[0][2]=-pTS->getTrackState(0)*cosf;
- J0[1][0]= cosf;J0[1][2]=-pTS->getTrackState(0)*sinf;
- J0[2][1]=1.0;
- J0[3][2]=-sinf*sint;J0[3][3]=cosf*cost;
- J0[4][2]= cosf*sint;J0[4][3]=sinf*cost;
- J0[5][3]=-sint;
- J0[6][4]=1.0;
- }
- for(i=0;i<4;i++) D[i]=pSE->getPar(i);
- for(i=0;i<3;i++) gPi[i]=gP[i];
-
- getMagneticField(gP,gB);
-
- for(i=0;i<3;i++) gBi[i]=gB[i];
-
- c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
- b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
- a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
- gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
- gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
-
- descr=b*b-4.0*a*c;
-
- if(descr<0.0)
- {
- // printf("D<0 - extrapolation failed\n");
- return NULL;
- }
-
- bool useExpansion=true;
- double ratio = 4*a*c/(b*b);
-
- if(fabs(ratio)>0.1)
- useExpansion = false;
-
- if(useExpansion) {
- sl=-c/b;
- sl=sl*(1-a*sl/b);
- }
- else {
- int signb = (b<0.0)?-1:1;
- sl = (-b+signb*sqrt(descr))/(2*a);
- }
-
- if(fabs(sl)<minStep) nStepMax=1;
- else
- {
- nStepMax=(int)(fabs(sl)/minStep)+1;
- }
- if((nStepMax<0)||(nStepMax>1000))
- {
- return NULL;
- }
- Av=sl*CQ;
- Ac=0.5*sl*Av;
- DVx=gV[1]*gB[2]-gV[2]*gB[1];
- DVy=gV[2]*gB[0]-gV[0]*gB[2];
- DVz=gV[0]*gB[1]-gV[1]*gB[0];
-
- P[0]=gP[0]+gV[0]*sl+Ac*DVx;
- P[1]=gP[1]+gV[1]*sl+Ac*DVy;
- P[2]=gP[2]+gV[2]*sl+Ac*DVz;
- V[0]=gV[0]+Av*DVx;
- V[1]=gV[1]+Av*DVy;
- V[2]=gV[2]+Av*DVz;
-
- getMagneticField(P,gB);
-
- for(i=0;i<3;i++) gBf[i]=gB[i];
- for(i=0;i<3;i++)
- {
- dBds[i]=(gBf[i]-gBi[i])/sl;
- gB[i]=gBi[i];
- }
- nStep=nStepMax;path=0.0;
- while(nStep>0)
- {
- c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
- b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
- a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
- gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
- gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
-
- ratio = 4*a*c/(b*b);
- if(fabs(ratio)>0.1)
- useExpansion = false;
- else useExpansion = true;
-
- if(useExpansion) {
- sl=-c/b;
- sl=sl*(1-a*sl/b);
- }
- else {
- descr=b*b-4.0*a*c;
- if(descr<0.0)
- {
- // printf("D<0 - extrapolation failed\n");
- return NULL;
- }
- int signb = (b<0.0)?-1:1;
- sl = (-b+signb*sqrt(descr))/(2*a);
- }
-
- ds=sl/nStep;path+=ds;
- Av=ds*CQ;
- Ac=0.5*ds*Av;
- DVx=gV[1]*gB[2]-gV[2]*gB[1];
- DVy=gV[2]*gB[0]-gV[0]*gB[2];
- DVz=gV[0]*gB[1]-gV[1]*gB[0];
-
- P[0]=gP[0]+gV[0]*ds+Ac*DVx;
- P[1]=gP[1]+gV[1]*ds+Ac*DVy;
- P[2]=gP[2]+gV[2]*ds+Ac*DVz;
- V[0]=gV[0]+Av*DVx;
- V[1]=gV[1]+Av*DVy;
- V[2]=gV[2]+Av*DVz;
- for(i=0;i<3;i++)
- {
- gV[i]=V[i];gP[i]=P[i];
- }
- for(i=0;i<3;i++) gB[i]+=dBds[i]*ds;
- nStep--;
- }
- pSE->transformPointToLocal(gP,lP);
- Rf[0]=lP[0];Rf[1]=lP[1];
- Rf[2]=atan2(V[1],V[0]);
-
- if(fabs(V[2])>1.0)
- {
- return NULL;
- }
-
- Rf[3]=acos(V[2]);
- Rf[4]=pTS->getTrackState(4);
-
- gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=cost;
-
- for(i=0;i<4;i++) D[i]=pSE->getPar(i);
- for(i=0;i<3;i++) gP[i]=gPi[i];
-
- for(i=0;i<3;i++)
- {
- gB[i]=0.5*(gBi[i]+gBf[i]);
- }
-
- c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
- b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
- a=CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
- gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
- gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
-
- ratio = 4*a*c/(b*b);
- if(fabs(ratio)>0.1)
- useExpansion = false;
- else useExpansion = true;
-
- if(useExpansion) {
- s=-c/b;
- s=s*(1-a*s/b);
- }
- else {
- descr=b*b-4.0*a*c;
- if(descr<0.0)
- {
- // printf("D<0 - extrapolation failed\n");
- return NULL;
- }
- int signb = (b<0.0)?-1:1;
- s = (-b+signb*sqrt(descr))/(2*a);
- }
-
- Av=s*CQ;
- Ac=0.5*s*Av;
- Cc=0.5*s*s*C;
-
- DVx=gV[1]*gB[2]-gV[2]*gB[1];
- DVy=gV[2]*gB[0]-gV[0]*gB[2];
- DVz=gV[0]*gB[1]-gV[1]*gB[0];
-
- P[0]=gP[0]+gV[0]*s+Ac*DVx;
- P[1]=gP[1]+gV[1]*s+Ac*DVy;
- P[2]=gP[2]+gV[2]*s+Ac*DVz;
-
- V[0]=gV[0]+Av*DVx;V[1]=gV[1]+Av*DVy;V[2]=gV[2]+Av*DVz;
-
- pSE->transformPointToLocal(P,lP);
-
- memset(&Jm[0][0],0,sizeof(Jm));
-
- for(i=0;i<3;i++) for(j=0;j<3;j++) M[i][j]=pSE->getRotMatrix(i,j);
-
- double coeff[3], dadVx,dadVy,dadVz,dadQ,dsdx,dsdy,dsdz,dsdVx,dsdVy,dsdVz,dsdQ;
- coeff[0]=-c*c/(b*b*b);
- coeff[1]=c*(1.0+3.0*c*a/(b*b))/(b*b);
- coeff[2]=-(1.0+2.0*c*a/(b*b))/b;
-
- dadVx=0.5*CQ*(-D[1]*gB[2]+D[2]*gB[1]);
- dadVy=0.5*CQ*( D[0]*gB[2]-D[2]*gB[0]);
- dadVz=0.5*CQ*(-D[0]*gB[1]+D[1]*gB[0]);
- dadQ=0.5*C*(D[0]*DVx+D[1]*DVy+D[2]*DVz);
-
- dsdx=coeff[2]*D[0];
- dsdy=coeff[2]*D[1];
- dsdz=coeff[2]*D[2];
- dsdVx=coeff[0]*dadVx+coeff[1]*D[0];
- dsdVy=coeff[0]*dadVy+coeff[1]*D[1];
- dsdVz=coeff[0]*dadVz+coeff[1]*D[2];
- dsdQ=coeff[0]*dadQ;
-
- Jm[0][0]=1.0+V[0]*dsdx;
- Jm[0][1]= V[0]*dsdy;
- Jm[0][2]= V[0]*dsdz;
-
- Jm[0][3]= s+V[0]*dsdVx;
- Jm[0][4]= V[0]*dsdVy+Ac*gB[2];
- Jm[0][5]= V[0]*dsdVz-Ac*gB[1];
- Jm[0][6]= V[0]*dsdQ+Cc*DVx;
-
- Jm[1][0]= V[1]*dsdx;
- Jm[1][1]=1.0+V[1]*dsdy;
- Jm[1][2]= V[1]*dsdz;
-
- Jm[1][3]= V[1]*dsdVx-Ac*gB[2];
- Jm[1][4]= s+V[1]*dsdVy;
- Jm[1][5]= V[1]*dsdVz+Ac*gB[0];
- Jm[1][6]= V[1]*dsdQ+Cc*DVy;
-
- Jm[2][0]= V[2]*dsdx;
- Jm[2][1]= V[2]*dsdy;
- Jm[2][2]=1.0+V[2]*dsdz;
- Jm[2][3]= V[2]*dsdVx+Ac*gB[1];
- Jm[2][4]= V[2]*dsdVy-Ac*gB[0];
- Jm[2][5]= s+V[2]*dsdVz;
- Jm[2][6]= V[2]*dsdQ+Cc*DVz;
-
- Jm[3][0]=dsdx*CQ*DVx;
- Jm[3][1]=dsdy*CQ*DVx;
- Jm[3][2]=dsdz*CQ*DVx;
-
- Jm[3][3]=1.0+dsdVx*CQ*DVx;
- Jm[3][4]=CQ*(dsdVy*DVx+s*gB[2]);
- Jm[3][5]=CQ*(dsdVz*DVx-s*gB[1]);
-
- Jm[3][6]=(CQ*dsdQ+C*s)*DVx;
-
- Jm[4][0]=dsdx*CQ*DVy;
- Jm[4][1]=dsdy*CQ*DVy;
- Jm[4][2]=dsdz*CQ*DVy;
-
- Jm[4][3]=CQ*(dsdVx*DVy-s*gB[2]);
- Jm[4][4]=1.0+dsdVy*CQ*DVy;
- Jm[4][5]=CQ*(dsdVz*DVy+s*gB[0]);
-
- Jm[4][6]=(CQ*dsdQ+C*s)*DVy;
-
- Jm[5][0]=dsdx*CQ*DVz;
- Jm[5][1]=dsdy*CQ*DVz;
- Jm[5][2]=dsdz*CQ*DVz;
- Jm[5][3]=CQ*(dsdVx*DVz+s*gB[1]);
- Jm[5][4]=CQ*(dsdVy*DVz-s*gB[0]);
- Jm[5][5]=1.0+dsdVz*CQ*DVz;
- Jm[5][6]=(CQ*dsdQ+C*s)*DVz;
-
- Jm[6][6]=1.0;
-
- memset(&J1[0][0],0,sizeof(J1));
-
- J1[0][0]=M[0][0];J1[0][1]=M[0][1];J1[0][2]=M[0][2];
- J1[1][0]=M[1][0];J1[1][1]=M[1][1];J1[1][2]=M[1][2];
- J1[2][3]=-V[1]/(V[0]*V[0]+V[1]*V[1]);
- J1[2][4]= V[0]/(V[0]*V[0]+V[1]*V[1]);
- J1[3][5]=-1.0/sqrt(1-V[2]*V[2]);
- J1[4][6]=1.0;
-
- for(i=0;i<7;i++)
- {
- for(j=0;j<2;j++)
- Buf[j][i]=J1[j][0]*Jm[0][i]+J1[j][1]*Jm[1][i]+J1[j][2]*Jm[2][i];
- Buf[2][i]=J1[2][3]*Jm[3][i]+J1[2][4]*Jm[4][i];
- Buf[3][i]=J1[3][5]*Jm[5][i];
- Buf[4][i]=Jm[6][i];
- }
-
- if(pSB!=NULL)
- {
- for(i=0;i<5;i++)
- {
- J[i][0]=Buf[i][0]*J0[0][0]+Buf[i][1]*J0[1][0]+Buf[i][2]*J0[2][0];
- J[i][1]=Buf[i][0]*J0[0][1]+Buf[i][1]*J0[1][1]+Buf[i][2]*J0[2][1];
- J[i][2]=Buf[i][3]*J0[3][2]+Buf[i][4]*J0[4][2];
- J[i][3]=Buf[i][3]*J0[3][3]+Buf[i][4]*J0[4][3]+Buf[i][5]*J0[5][3];
- J[i][4]=Buf[i][6];
- }
- }
- else
- {
- for(i=0;i<5;i++)
- {
- J[i][0]=Buf[i][0]*J0[0][0]+Buf[i][1]*J0[1][0];
- J[i][1]=Buf[i][2];
- J[i][2]=Buf[i][0]*J0[0][2]+Buf[i][1]*J0[1][2]+Buf[i][3]*J0[3][2]+Buf[i][4]*J0[4][2];
- J[i][3]=Buf[i][3]*J0[3][3]+Buf[i][4]*J0[4][3]+Buf[i][5]*J0[5][3];
- J[i][4]=Buf[i][6];
- }
- }
- }
- else {
- Rf[0]=pTS->getTrackState(0);
- Rf[1]=pTS->getTrackState(1);
- Rf[2]=pTS->getTrackState(2);
- Rf[3]=pTS->getTrackState(3);
- Rf[4]=pTS->getTrackState(4);
- memset(&J[0][0],0,sizeof(J));
- for(i=0;i<5;i++) J[i][i]=1.0;
- }
-
- for(i=0;i<5;i++) for(j=0;j<5;j++)
- {
- AG[i][j]=0.0;for(m=0;m<5;m++) AG[i][j]+=J[i][m]*pTS->getTrackCovariance(m,j);
- }
- for(i=0;i<5;i++) for(j=i;j<5;j++)
- {
- Gf[i][j]=0.0;
- for(m=0;m<5;m++) Gf[i][j]+=AG[i][m]*J[j][m];
- Gf[j][i]=Gf[i][j];
- }
-
- Trk::TrkTrackState* pTE=new Trk::TrkTrackState(pTS);
-
- //workaround to keep using existing TrkTrackState code
- double Rtmp[5];
-
- for(i=0;i<4;i++) Rtmp[i] = Rf[i];
- Rtmp[4] = 0.001*Rf[4];//GeV->MeV
-
- pTE->setTrackState(Rtmp);
- pTE->setTrackCovariance(Gf);
- pTE->attachToSurface(pSE);
-
- // pTE->applyMaterialEffects();
- if(m_doMultScatt)
- pTE->applyMultipleScattering();
-
- pTE->setTrackState(Rf);//restore
-
- if(m_doBremm)
- pTE->applyEnergyLoss(1);
-
- AmgSymMatrix(5) Gi;
- for(i=0;i<5;i++) for(j=i;j<5;j++)
- {
- Gi(i,j)=pTE->getTrackCovariance(i,j);
- Gi.fillSymmetric(i, j, pTE->getTrackCovariance(i,j));
- }
- Gi = Gi.inverse();
-
- for(i=0;i<5;i++) for(j=0;j<5;j++)
- {
- A[i][j]=0.0;
- for(m=0;m<5;m++) A[i][j]+=AG[m][i]*Gi(m,j);
- }
- pTE->setPreviousState(pTS);
- pTE->setSmootherGain(A);
-
- return pTE;
-}
-
-StatusCode TrigL2ResidualCalculator::getUnbiasedResiduals(const Trk::Track& pT,
- std::vector<TrigL2HitResidual>& vResid) const {
-
- const Trk::TrackParameters* trackPars = pT.perigeeParameters();
- if(trackPars==nullptr) {
- ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
- return StatusCode::FAILURE;
- }
- std::vector<Trk::TrkBaseNode*> vpTrkNodes;
- std::vector<Trk::TrkTrackState*> vpTrackStates;
- vResid.clear();
- double trk_theta = trackPars->parameters()[Trk::theta];
- double trk_qOverP = trackPars->parameters()[Trk::qOverP];
- double Pt = sin(trk_theta)/trk_qOverP;
- if(fabs(Pt)<100.0)
- {
- ATH_MSG_DEBUG("TrigL2ResidualCalculator failed -- Estimated Pt is too low "<<Pt);
- return StatusCode::FAILURE;
- }
-
- // 1. Create filtering nodes
-
- bool trackResult = m_trackMaker->createDkfTrack(pT,vpTrkNodes, m_DChi2);
- if(!trackResult)
- {
- ATH_MSG_DEBUG("TrigDkfTrackMaker failed");
- return StatusCode::FAILURE;
- }
-
- int nNodeIndex=-1;
- bool OK=true;
- std::vector<Trk::TrkBaseNode*>::iterator pnIt,pnEnd(vpTrkNodes.end());
-
- for(std::vector<Trk::TrkBaseNode*>::iterator pNodeIt=vpTrkNodes.begin();pNodeIt!=vpTrkNodes.end();
- ++pNodeIt)
- {
- nNodeIndex++;
- Trk::TrkBaseNode* pMaskedNode=(*pNodeIt);
- Trk::TrkTrackState* pMaskedState=NULL;
-
- // 2. Create initial track state:
-
- double Rk[5];
- Rk[0] = trackPars->parameters()[Trk::d0];
- Rk[1] = trackPars->parameters()[Trk::z0];
- Rk[2] = trackPars->parameters()[Trk::phi0];
- if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
- if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
- trk_theta = trackPars->parameters()[Trk::theta];
- Rk[3] = trk_theta;
- Rk[4] = 1000.0*trackPars->parameters()[Trk::qOverP];//MeV->GeV
- //No need to correct scale back - not returning track
-
- // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
-
- Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
- double Gk[5][5] = {{100.0, 0, 0, 0, 0},
- {0, 100.0, 0, 0, 0},
- {0, 0, 0.01, 0, 0},
- {0, 0, 0, 0.01, 0},
- {0, 0, 0, 0, 0.1}};
- pTS->setTrackCovariance(Gk);
- if(m_doMultScatt)
- pTS->setScatteringMode(1);
- if(m_doBremm)
- pTS->setScatteringMode(2);
- vpTrackStates.push_back(pTS);
-
- ATH_MSG_DEBUG("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
- <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]);
-
- OK=true;
- Trk::TrkPlanarSurface *pSB=NULL,*pSE;
-
- for(pnIt=vpTrkNodes.begin();pnIt!=pnEnd;++pnIt)
- {
- pSE=(*pnIt)->getSurface();
- Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE);
-
- pSB=pSE;
- if(pNS!=NULL)
- {
- vpTrackStates.push_back(pNS);
-
- (*pnIt)->validateMeasurement(pNS);
- ATH_MSG_DEBUG("dChi2="<<(*pnIt)->getChi2());
- if((*pnIt)!=pMaskedNode)
- {
- (*pnIt)->updateTrackState(pNS);
- }
- else
- {
- pMaskedState=pNS;
- }
- pTS=pNS;
- Pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
- if(fabs(Pt)<0.2)
- {
- ATH_MSG_DEBUG("Estimated Pt is too low "<<Pt<<" - skipping fit");
- OK=false;break;
- }
- }
- else
- {
- OK=false;break;
- }
- }
- if(OK)
- {
- std::vector<Trk::TrkTrackState*>::reverse_iterator ptsrIt(vpTrackStates.rbegin()),
- ptsrEnd(vpTrackStates.rend());
-
- for(;ptsrIt!=ptsrEnd;++ptsrIt)
- {
- (*ptsrIt)->runSmoother();
- }
-
- pMaskedNode->validateMeasurement(pMaskedState);
-
- double r[2],V[2][2];
-
- int nSize=pMaskedNode->getResiduals(r);
- nSize=pMaskedNode->getInverseResidualVariance(V);
- const Trk::PrepRawData* pPRD = pMaskedNode->getPrepRawData();
-
- Identifier id = pPRD->identify();
-
- Region region = Region::Undefined;
- if(m_idHelper->is_pixel(id))
- {
- region=(m_pixelId->is_barrel(id)) ? Region::PixBarrel: Region::PixEndcap;
- if (m_pixelId->is_blayer(id)) {
- region = Region::IBL;
- }
- }
- if(m_idHelper->is_sct(id))
- {
- region=(m_sctId->is_barrel(id)) ? Region::SctBarrel : Region::SctEndcap;
- }
- if(nSize==1) {
- if(V[0][0]>0.0) {
- vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0])));
- }
- else {
- OK=false;
- break;
- }
- }
- else {
- if((V[0][0]>0.0) && (V[1][1]>0.0)) {
- vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0]),
- r[1],r[1]*sqrt(V[1][1])));
- }
- else {
- OK=false;
- break;
- }
- }
- }
- else
- {
- ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
- vResid.clear();
- }
- for(std::vector<Trk::TrkTrackState*>::iterator ptsIt=vpTrackStates.begin();
- ptsIt!=vpTrackStates.end();++ptsIt) delete (*ptsIt);
- vpTrackStates.clear();
- if(!OK) break;
- }
- pnIt=vpTrkNodes.begin();pnEnd=vpTrkNodes.end();
- for(;pnIt!=pnEnd;++pnIt)
- {
- delete((*pnIt)->getSurface());
- delete (*pnIt);
- }
- vpTrkNodes.clear();
-
- if(OK) return StatusCode::SUCCESS;
- else return StatusCode::FAILURE;
-
-}
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigL2ResidualCalculator.h b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigL2ResidualCalculator.h
deleted file mode 100644
index c8d01009f81968880bb8fe230d190507d1092160..0000000000000000000000000000000000000000
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigL2ResidualCalculator.h
+++ /dev/null
@@ -1,58 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef __TRIGL2RESIDUALCALCULATOR_H__
-#define __TRIGL2RESIDUALCALCULATOR_H__
-
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "GaudiKernel/ToolHandle.h"
-#include "GaudiKernel/ServiceHandle.h"
-
-#include "TrigInDetToolInterfaces/ITrigL2ResidualCalculator.h"
-#include "TrigInDetToolInterfaces/ITrigDkfTrackMakerTool.h"
-
-class AtlasDetectorID;
-class PixelID;
-class SCT_ID;
-class TrigL2HitResidual;
-
-namespace Trk {
- class TrkBaseNode;
- class TrkTrackState;
- class TrkPlanarSurface;
- class Track;
-}
-
-namespace MagField {
- class IMagFieldSvc;
-}
-
-class TrigL2ResidualCalculator: public AthAlgTool, virtual public ITrigL2ResidualCalculator
-{
- public:
-
- TrigL2ResidualCalculator( const std::string&, const std::string&, const IInterface* );
- virtual StatusCode initialize();
-
- StatusCode getUnbiasedResiduals(const Trk::Track&, std::vector<TrigL2HitResidual>&) const;
-
-private:
-
- Trk::TrkTrackState* extrapolate(Trk::TrkTrackState*,
- Trk::TrkPlanarSurface*,
- Trk::TrkPlanarSurface*) const;
- void getMagneticField(double[3],double*) const;
-
- double m_DChi2;
- bool m_doMultScatt;
- bool m_doBremm;
- bool m_offlineClusters;
- ServiceHandle<MagField::IMagFieldSvc> m_MagFieldSvc;
- ToolHandle<ITrigDkfTrackMakerTool> m_trackMaker;
- const PixelID* m_pixelId;
- const SCT_ID* m_sctId;
- const AtlasDetectorID* m_idHelper;
-};
-
-#endif
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/components/TrigInDetTrackFitter_entries.cxx b/Trigger/TrigTools/TrigInDetTrackFitter/src/components/TrigInDetTrackFitter_entries.cxx
index eb424e52f93547b1f399cd0d03bea805114d5a3f..e4f8d28075c31fc98f2b3a78e105a3b030f4490b 100644
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/components/TrigInDetTrackFitter_entries.cxx
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/src/components/TrigInDetTrackFitter_entries.cxx
@@ -1,7 +1,5 @@
#include "../TrigInDetTrackFitter.h"
#include "../TrigDkfTrackMakerTool.h"
-#include "../TrigL2ResidualCalculator.h"
DECLARE_COMPONENT( TrigInDetTrackFitter )
DECLARE_COMPONENT( TrigDkfTrackMakerTool )
-DECLARE_COMPONENT( TrigL2ResidualCalculator )