diff --git a/Trigger/TrigFTK/FTK_RecTools/CMakeLists.txt b/Trigger/TrigFTK/FTK_RecTools/CMakeLists.txt
index b6f47d680a707522ef9ee448358e5e5f9222ad8f..6e68377a4954dee3fa433e7c838cdb7e9fe0af37 100644
--- a/Trigger/TrigFTK/FTK_RecTools/CMakeLists.txt
+++ b/Trigger/TrigFTK/FTK_RecTools/CMakeLists.txt
@@ -7,18 +7,34 @@ atlas_subdir( FTK_RecTools )
# Declare the package's dependencies:
atlas_depends_on_subdirs( PUBLIC
+ AtlasPolicy
Control/AthenaBaseComps
+ Event/xAOD/xAODTracking
GaudiKernel
Tools/PyJobTransforms
+ DetectorDescription/GeoPrimitives
Tracking/TrkEvent/TrkTrack
Tracking/TrkEvent/VxVertex
+ Tracking/TrkVertexFitter/TrkVxEdmCnv
Trigger/TrigFTK/FTK_DataProviderInterfaces
Trigger/TrigFTK/FTK_RecToolInterfaces
Trigger/TrigFTK/TrigFTK_RawData
- Event/xAOD/xAODTracking
- Tracking/TrkVertexFitter/TrkVxEdmCnv
- PRIVATE
- Tracking/TrkEvent/TrkParameters)
+ InnerDetector/InDetDetDescr/SCT_ModuleDistortions
+ InnerDetector/InDetRecEvent/InDetPrepRawData
+ InnerDetector/InDetRecEvent/InDetRIO_OnTrack
+ Tracking/TrkEvent/TrkParameters
+ Tracking/TrkTools/TrkToolInterfaces
+ Tracking/TrkTools/TrkAmbiguityProcessor
+ PRIVATE
+ Event/EventPrimitives
+ InnerDetector/InDetConditions/PixelConditionsServices
+ InnerDetector/InDetConditions/PixelConditionsTools
+ InnerDetector/InDetDetDescr/InDetIdentifier
+ InnerDetector/InDetDetDescr/InDetReadoutGeometry
+ InnerDetector/InDetDetDescr/PixelGeoModel
+ InnerDetector/InDetRecTools/SiClusterizationTool
+ Tracking/TrkDetDescr/TrkSurfaces )
+
# External dependencies:
find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread )
@@ -28,10 +44,10 @@ atlas_add_library( FTK_RecToolsLib
src/*.cxx
PUBLIC_HEADERS FTK_RecTools
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel TrkTrack TrkParameters VxVertex FTK_DataProviderInterfaces TrigFTK_RawData xAODTracking TrkVxEdmCnvLib)
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GeoPrimitives xAODTracking TrkTrack VxVertex FTK_DataProviderInterfaces TrigFTK_RawData TrkVxEdmCnvLib InDetPrepRawData InDetRIO_OnTrack TrkParameters TrkToolInterfaces TrkAmbiguityProcessorLib StoreGateLib SGtests EventPrimitives InDetIdentifier InDetReadoutGeometry SiClusterizationToolLib TrkSurfaces)
atlas_add_component( FTK_RecTools
- src/components/*.cxx
+ src/*.cxx src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel TrkTrack TrkParameters VxVertex FTK_DataProviderInterfaces TrigFTK_RawData xAODTracking TrkVxEdmCnvLib FTK_RecToolsLib )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GeoPrimitives xAODTracking TrkTrack VxVertex TrkVxEdmCnvLib FTK_DataProviderInterfaces TrigFTK_RawData TrkParameters FTK_RecToolsLib InDetPrepRawData InDetRIO_OnTrack TrkParameters TrkToolInterfaces TrkAmbiguityProcessorLib StoreGateLib SGtests EventPrimitives InDetIdentifier InDetReadoutGeometry SiClusterizationToolLib TrkSurfaces)
diff --git a/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_PixelClusterOnTrackTool.h b/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_PixelClusterOnTrackTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..8533a1e0c21dc0867eef4fb84dd22bca97a57017
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_PixelClusterOnTrackTool.h
@@ -0,0 +1,207 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// Header file for class FTK_PixelClusterOnTrackTool
+///////////////////////////////////////////////////////////////////
+
+#ifndef FTK_PixelClusterOnTrackTool_H
+#define FTK_PixelClusterOnTrackTool_H
+
+#include "GaudiKernel/ToolHandle.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+
+
+#include "TrkToolInterfaces/IRIO_OnTrackCreator.h"
+#include "TrkToolInterfaces/IRIO_OnTrackErrorScalingTool.h"
+#include "InDetRIO_OnTrack/PixelClusterOnTrack.h"
+
+#include "InDetPrepRawData/PixelGangedClusterAmbiguities.h"
+#include "TrkParameters/TrackParameters.h"
+//#include "InDetIdentifier/PixelID.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "TrkAmbiguityProcessor/dRMap.h"
+#include "SiClusterizationTool/NnClusterizationFactory.h"
+
+//#include "PixelConditionsServices/IPixelOfflineCalibSvc.h"
+//#include "PixelConditionsTools/IModuleDistortionsTool.h"
+
+#include "AthenaPoolUtilities/CondAttrListCollection.h"
+class PixelID;
+class IPixelOfflineCalibSvc;
+class IModuleDistortionsTool;
+
+class StoreGateSvc;
+class IBLParameterSvc;
+
+//namespace InDet {
+
+ /** @brief creates FTK_PixelClusterOnTrack objects allowing to
+ calibrate cluster position and error using a given track hypothesis.
+
+ See doxygen of Trk::RIO_OnTrackCreator for details.
+ Different strategies to calibrate the cluster error can be chosen
+ by job Option. Also the handle to the general hit-error scaling
+ is implemented.
+
+ Special strategies for correction can be invoked by calling the
+ correct method with an additional argument from the
+ PixelClusterStrategy enumeration
+
+ */
+
+
+namespace InDet {
+ enum PixelClusterStrategy {
+ PIXELCLUSTER_DEFAULT=0,
+ PIXELCLUSTER_OUTLIER=1,
+ PIXELCLUSTER_SHARED =2,
+ PIXELCLUSTER_SPLIT =3
+ };
+}
+
+ class FTK_PixelClusterOnTrackTool:
+ public AthAlgTool, virtual public Trk::IRIO_OnTrackCreator
+{
+ ///////////////////////////////////////////////////////////////////
+ // Public methods:
+ ///////////////////////////////////////////////////////////////////
+
+public:
+
+ //! AlgTool constructor
+ FTK_PixelClusterOnTrackTool(const std::string&,const std::string&,
+ const IInterface*);
+ virtual ~FTK_PixelClusterOnTrackTool ();
+ //! AlgTool initialisation
+ virtual StatusCode initialize() override;
+ //! AlgTool termination
+ virtual StatusCode finalize () override;
+
+
+
+ void correctBow(const Identifier&, Amg::Vector2D& locpos, const double tanphi, const double taneta) const;
+
+ double splineIBLPullX(float x, int layer) const;
+
+ /** @brief produces a PixelClusterOnTrack (object factory!).
+
+ Depending on job options it changes the pixel cluster position
+ and error according to the parameters (in particular, the angle)
+ of the intersecting track.
+ */
+ virtual const InDet::PixelClusterOnTrack* correct(const Trk::PrepRawData&,
+ const Trk::TrackParameters&) const;
+
+ virtual const InDet::PixelClusterOnTrack* correctDefault(const Trk::PrepRawData&,
+ const Trk::TrackParameters&) const;
+
+ virtual const InDet::PixelClusterOnTrack* correctNN(const Trk::PrepRawData&, const Trk::TrackParameters&) const;
+
+ virtual bool getErrorsDefaultAmbi( const InDet::PixelCluster*, const Trk::TrackParameters&,
+ Amg::Vector2D&, Amg::MatrixX&) const;
+
+ virtual bool getErrorsTIDE_Ambi( const InDet::PixelCluster*, const Trk::TrackParameters&,
+ Amg::Vector2D&, Amg::MatrixX&) const;
+
+ virtual const InDet::PixelClusterOnTrack* correct
+ (const Trk::PrepRawData&, const Trk::TrackParameters&,
+ const InDet::PixelClusterStrategy) const;
+
+
+ ///////////////////////////////////////////////////////////////////
+ // Private methods:
+ ///////////////////////////////////////////////////////////////////
+
+ private:
+
+ /** @brief parametrizes the pixel cluster position error as a function of
+ the track angle alpha and the cluster width (number of rows) deltax */
+ // double getBarrelPhiError(double& alpha, int& deltax) const;
+ // double getBarrelEtaError(double eta, int deltax, int deltay) const;
+ // double getEndcapPhiError(int etasize, int phisize) const;
+ // double getEndcapEtaError(int etasize, int phisize) const;
+
+ void FillFromDataBase() const;
+
+ ///////////////////////////////////////////////////////////////////
+ // Private data:
+ ///////////////////////////////////////////////////////////////////
+
+ ToolHandle<IModuleDistortionsTool> m_pixDistoTool ;
+ ToolHandle<Trk::IRIO_OnTrackErrorScalingTool> m_errorScalingTool;
+ ServiceHandle<IPixelOfflineCalibSvc> m_calibSvc ;
+ StoreGateSvc* m_detStore ;
+ /* ME: Test histos have nothing to do with production code, use a flag
+ IHistogram1D* m_h_Resx;
+ IHistogram1D* m_h_Resy;
+ IHistogram1D* m_h_Locx;
+ IHistogram1D* m_h_Locy;
+ IHistogram1D* m_h_PhiTrack;
+ IHistogram1D* m_h_ThetaTrack;
+ IHistogram1D* m_h_Rad;
+ IHistogram1D* m_h_Slope;
+ */
+
+ //! toolhandle for central error scaling
+ //! flag storing if errors need scaling or should be kept nominal
+ bool m_scalePixelCov ;
+ bool m_disableDistortions;
+ bool m_rel13like ;
+ int m_positionStrategy ;
+ mutable int m_errorStrategy ;
+
+
+ /** @brief Flag controlling how module distortions are taken into account:
+
+ case 0 -----> No distorsions implemented;
+
+ case 1 -----> Set curvature (in 1/meter) and twist (in radiant) equal for all modules;
+
+ case 2 -----> Read curvatures and twists from textfile containing Survey data;
+
+ case 3 -----> Set curvature and twist from Gaussian random generator with mean and RMS coming from Survey data;
+
+ case 4 -----> Read curvatures and twists from database (not ready yet);
+ */
+ //! identifier-helper
+ const PixelID* m_pixelid;
+
+ /** Enable NN based calibration (do only if NN calibration is applied) **/
+ mutable bool m_applyNNcorrection;
+ mutable bool m_applydRcorrection;
+ bool m_NNIBLcorrection;
+ bool m_IBLAbsent;
+
+ /** NN clusterizationi factory for NN based positions and errors **/
+ ToolHandle<InDet::NnClusterizationFactory> m_NnClusterizationFactory;
+ ServiceHandle<StoreGateSvc> m_storeGate; //!< Event store
+ ServiceHandle<IBLParameterSvc> m_IBLParameterSvc;
+
+
+ SG::ReadHandleKey<InDet::DRMap> m_dRMap; //!< the actual dR map
+ std::string m_dRMapName;
+
+ bool m_doNotRecalibrateNN;
+ bool m_noNNandBroadErrors;
+
+ /** Enable different treatment of cluster errors based on NN information (do only if TIDE ambi is run) **/
+ bool m_usingTIDE_Ambi;
+ std::string m_splitClusterMapName; //no longer used
+ mutable std::vector< std::vector<float> > m_fX, m_fY, m_fB, m_fC, m_fD;
+
+ //moved from static to member variable
+ static constexpr int nbinphi=9;
+ static constexpr int nbineta=6;
+ double calphi[nbinphi];
+ double caleta[nbineta][3];
+ double calerrphi[nbinphi][3];
+ double calerreta[nbineta][3];
+ double phix[nbinphi+1];
+ double etax[nbineta+1];
+};
+
+//} // end of namespace FTK
+
+#endif // FTK_PixelClusterOnTrackTool_H
diff --git a/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_SCTClusterOnTrackTool.h b/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_SCTClusterOnTrackTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..c6f76ff562169df433d7cc2e78bfa5787b743000
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_SCTClusterOnTrackTool.h
@@ -0,0 +1,107 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// Header file for class FTK_SCTClusterOnTrackTool
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+// Interface for SCT_ClusterOnTrack production
+///////////////////////////////////////////////////////////////////
+// started 1/05/2004 I.Gavrilenko - see ChangeLog for details
+///////////////////////////////////////////////////////////////////
+
+#ifndef FTK_SCTClusterOnTrackTool_H
+#define FTK_SCTClusterOnTrackTool_H
+
+#include "GaudiKernel/ToolHandle.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "TrkToolInterfaces/IRIO_OnTrackCreator.h"
+#include "TrkToolInterfaces/IRIO_OnTrackErrorScalingTool.h"
+#include "TrkParameters/TrackParameters.h"
+#include "InDetRIO_OnTrack/SCT_ClusterOnTrack.h"
+#include "SCT_ModuleDistortions/ISCT_ModuleDistortionsTool.h"
+
+//namespace InDet {
+
+
+/** @brief creates SCT_ClusterOnTrack objects allowing to
+ calibrate cluster position and error using a given track hypothesis.
+
+ See doxygen of Trk::RIO_OnTrackCreator for details.
+ Different strategies to calibrate the cluster error can be chosen
+ by job Option. Also the handle to the general hit-error scaling
+ is implemented.
+*/
+ class FTK_SCTClusterOnTrackTool:
+ public AthAlgTool,virtual public Trk::IRIO_OnTrackCreator
+{
+ ///////////////////////////////////////////////////////////////////
+ // Public methods:
+ ///////////////////////////////////////////////////////////////////
+
+public:
+
+ //! AlgTool constructor
+ FTK_SCTClusterOnTrackTool(const std::string&,const std::string&,const IInterface*);
+ virtual ~FTK_SCTClusterOnTrackTool ();
+ //! AlgTool initialisation
+ virtual StatusCode initialize();
+ //! AlgTool termination
+ virtual StatusCode finalize ();
+
+
+/** @brief produces an SCT_ClusterOnTrack using the measured
+ SCT_Cluster and the track prediction.
+
+ This method is a factory, so the client has to take care
+ of management/deletion of the SCT_ClusterOnTrack.
+ */
+ virtual const InDet::SCT_ClusterOnTrack* correct
+ (const Trk::PrepRawData&, const Trk::TrackParameters&) const;
+
+
+/** @brief Returns a correction to be applied to the SCT cluster local x position
+ in simulated events to remove a position bias introduced by the SCT digitisation.
+
+ @param [in] phi angle of track relative to Lorentz drift direction, in transverse plane
+ @param [in] nstrip SCT cluster size (number of strips)
+ */
+ double getCorrection(double phi, int nstrip) const;
+
+
+/** @brief Returns the resolution on the reconstructed position (local x) of SCT clusters
+ in simulated events.
+
+ @param [in] phi angle of track relative to Lorentz drift direction, in transverse plane
+ @param [in] nstrip SCT cluster size (number of strips)
+
+ The parameterisation of the resolution contained in getError() was derived from SCT
+ barrel clusters (80 micron pitch). It can be applied also to endcap clusters, after
+ rescaling to the appropriate pitch.
+ */
+ double getError(double phi, int nstrip) const;
+
+ private:
+
+ ///////////////////////////////////////////////////////////////////
+ // Private data:
+ ///////////////////////////////////////////////////////////////////
+
+ //! toolhandle for central error scaling
+ ToolHandle< Trk::IRIO_OnTrackErrorScalingTool > m_errorScalingTool;
+ ToolHandle<ISCT_ModuleDistortionsTool> m_distortionsTool;
+ //! flag storing if errors need scaling or should be kept nominal
+ bool m_scaleSctCov;
+
+ //! job options
+ bool m_option_make2dimBarrelClusters;
+ bool m_doDistortions ;//!< Flag to set Distortions
+ int m_option_errorStrategy;
+ int m_option_correctionStrategy;
+};
+
+//} // end of namespace FTK
+
+#endif // FTK_SCTClusterOnTrackTool_H
diff --git a/Trigger/TrigFTK/FTK_RecTools/src/FTK_PixelClusterOnTrackTool.cxx b/Trigger/TrigFTK/FTK_RecTools/src/FTK_PixelClusterOnTrackTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..aa756c64222fa1a075e9ddbb72d652e1218d9b81
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/src/FTK_PixelClusterOnTrackTool.cxx
@@ -0,0 +1,1058 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/////////////////////////////////////////////////////////////////
+// Implementation file for class FTK_PixelClusterOnTrackTool
+///////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////
+// AlgTool used for FTK_PixelClusterOnTrack object production
+///////////////////////////////////////////////////////////////////
+// Created 9/11/2016 J. Baines based on version by I.Gavrilenko
+///////////////////////////////////////////////////////////////////
+
+#include "FTK_RecTools/FTK_PixelClusterOnTrackTool.h"
+#include "InDetReadoutGeometry/SiDetectorManager.h"
+#include "InDetReadoutGeometry/PixelModuleDesign.h"
+#include "InDetIdentifier/PixelID.h"
+#include "PixelConditionsServices/IPixelOfflineCalibSvc.h"
+#include "PixelConditionsTools/IModuleDistortionsTool.h"
+#include "TrkSurfaces/PlaneSurface.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "GaudiKernel/IIncidentSvc.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include "PixelGeoModel/IBLParameterSvc.h"
+#include "InDetReadoutGeometry/SiDetectorElement.h"
+
+#include "CoralBase/AttributeListSpecification.h"
+#include "CoralBase/Attribute.h"
+#include "AthenaPoolUtilities/CondAttrListCollection.h"
+#include "AthenaPoolUtilities/AthenaAttributeList.h"
+using CLHEP::mm;
+using CLHEP::micrometer;
+
+///////////////////////////////////////////////////////////////////
+// Constructor
+///////////////////////////////////////////////////////////////////
+
+namespace
+{
+ // using x*x might be quicker than pow(x,2), depends on compiler optimisation
+ inline double
+ square(const double x) {
+ return x * x;
+ }
+
+ double
+ distance(const std::vector<Amg::Vector2D> &vectorOfPositions,
+ const std::vector<Amg::Vector2D> &allLocalPositions,
+ const std::vector<Amg::MatrixX> &allErrorMatrix,
+ const int i, const int j, const int k) {
+ return
+ square(vectorOfPositions[i][0] - allLocalPositions[0][0]) / allErrorMatrix[0](0, 0) +
+ square(vectorOfPositions[j][0] - allLocalPositions[1][0]) / allErrorMatrix[1](0, 0) +
+ square(vectorOfPositions[k][0] - allLocalPositions[2][0]) / allErrorMatrix[2](0, 0) +
+ square(vectorOfPositions[i][1] - allLocalPositions[0][1]) / allErrorMatrix[0](1, 1) +
+ square(vectorOfPositions[j][1] - allLocalPositions[1][1]) / allErrorMatrix[1](1, 1) +
+ square(vectorOfPositions[k][1] - allLocalPositions[2][1]) / allErrorMatrix[2](1, 1);
+ }
+
+ // avoid a lot of '/sqrt(12)' in loops, declare the constant here
+ // constexpr sqrt is a gcc extension, unsupported by clang.
+ //constexpr double TOPHAT_SIGMA = 1. / sqrt(12.);
+ static double TOPHAT_SIGMA = 1. / sqrt(12.);
+}
+
+
+FTK_PixelClusterOnTrackTool::FTK_PixelClusterOnTrackTool
+ (const std::string &t, const std::string &n, const IInterface *p) :
+ ::AthAlgTool(t, n, p),
+ m_pixDistoTool("PixelDistortionsTool", this),
+ m_errorScalingTool("Trk::RIO_OnTrackErrorScalingTool/RIO_OnTrackErrorScalingTool", this),
+ m_calibSvc("PixelOfflineCalibSvc", n),
+ m_detStore(nullptr),
+ m_scalePixelCov(false),
+ m_disableDistortions(false),
+ m_rel13like(false),
+ m_pixelid(nullptr),
+ m_applyNNcorrection(false),
+ m_applydRcorrection(false),
+ m_NNIBLcorrection(false),
+ m_IBLAbsent(true),
+ m_NnClusterizationFactory("InDet::NnClusterizationFactory/NnClusterizationFactory", this),
+ m_storeGate("StoreGateSvc", n),
+ m_IBLParameterSvc("IBLParameterSvc", n),
+ m_dRMap(""),
+ m_dRMapName("dRMap"),
+ m_doNotRecalibrateNN(false),
+ m_noNNandBroadErrors(false),
+ m_usingTIDE_Ambi(false),
+ m_splitClusterMapName("SplitClusterAmbiguityMap") {
+ declareInterface<IRIO_OnTrackCreator>(this);
+
+ declareProperty("ErrorScalingTool", m_errorScalingTool, "The error scaling tool");
+ declareProperty("PixelDistortionsTool", m_pixDistoTool, "Tool to retrieve pixel distortions");
+ declareProperty("PositionStrategy", m_positionStrategy = 1, "Which calibration of cluster positions");
+ declareProperty("ErrorStrategy", m_errorStrategy = 2, "Which calibration of cluster position errors");
+ declareProperty("DisableDistortions", m_disableDistortions, "Disable simulation of module distortions");
+ declareProperty("Release13like", m_rel13like, "Activate release-13 like settigs");
+ declareProperty("PixelOfflineCalibSvc", m_calibSvc, "Offline calibration svc");
+ declareProperty("applyNNcorrection", m_applyNNcorrection);
+ declareProperty("applydRcorrection", m_applydRcorrection);
+ declareProperty("NNIBLcorrection", m_NNIBLcorrection);
+ declareProperty("EventStore", m_storeGate);
+ declareProperty("NnClusterizationFactory", m_NnClusterizationFactory);
+ declareProperty("SplitClusterAmbiguityMap", m_splitClusterMapName);//Remove Later
+ declareProperty("dRMapName", m_dRMapName); //This is a string to prevent the scheduler seeing trkAmbSolver as both creating and requiring the map
+ declareProperty("doNotRecalibrateNN", m_doNotRecalibrateNN);
+ declareProperty("m_noNNandBroadErrors", m_noNNandBroadErrors);
+ declareProperty("RunningTIDE_Ambi", m_usingTIDE_Ambi);
+}
+
+///////////////////////////////////////////////////////////////////
+// Destructor
+///////////////////////////////////////////////////////////////////
+
+FTK_PixelClusterOnTrackTool::~FTK_PixelClusterOnTrackTool() {
+}
+
+///////////////////////////////////////////////////////////////////
+// Initialisation
+///////////////////////////////////////////////////////////////////
+
+StatusCode
+FTK_PixelClusterOnTrackTool::initialize() {
+ StatusCode sc = AthAlgTool::initialize();
+
+ ATH_MSG_INFO(name() << " initialize()");
+ ATH_MSG_DEBUG("Error strategy is" << m_errorStrategy);
+ ATH_MSG_DEBUG("Release 13 flag is" << m_rel13like);
+
+ if (m_IBLParameterSvc.retrieve().isFailure()) {
+ ATH_MSG_WARNING("Could not retrieve IBLParameterSvc");
+ } else {
+ m_IBLParameterSvc->setBoolParameters(m_applyNNcorrection, "doPixelClusterSplitting");
+ m_IBLParameterSvc->setBoolParameters(m_IBLAbsent, "IBLAbsent");
+ }
+
+ // get the offline calibration service
+ sc = m_calibSvc.retrieve();
+ if (sc.isFailure() || !m_calibSvc) {
+ ATH_MSG_ERROR(m_calibSvc.type() << " not found! ");
+ } else {
+ ATH_MSG_INFO("Retrieved tool " << m_calibSvc.type());
+ }
+
+ // get the error scaling tool
+ if (m_errorScalingTool.retrieve().isFailure()) {
+ msg(MSG::FATAL) << "Failed to retrieve tool " << m_errorScalingTool << endmsg;
+ return StatusCode::FAILURE;
+ } else {
+ ATH_MSG_INFO("Retrieved tool " << m_errorScalingTool);
+ m_scalePixelCov = m_errorScalingTool->needToScalePixel();
+ if (m_scalePixelCov) {
+ ATH_MSG_DEBUG("Detected need for scaling Pixel errors.");
+ }
+ }
+
+ // the NN corrections
+ if (m_applyNNcorrection) {
+ msg(MSG::WARNING) << " NNcorrection not possible for FTK clusters - using default correction" << endmsg;
+ m_applyNNcorrection=false;
+ }
+
+ // get the module distortions tool
+ if (!m_disableDistortions) {
+ if (!m_pixDistoTool.empty()) {
+ sc = m_pixDistoTool.retrieve();
+ if (sc != StatusCode::SUCCESS) {
+ msg(MSG::ERROR) << "Can't get pixel distortions tool " << endmsg;
+ } else {
+ ATH_MSG_INFO("Pixel distortions tool retrieved");
+ }
+ } else {
+ ATH_MSG_INFO("No PixelDistortionsTool selected.");
+ }
+ } else {
+ ATH_MSG_INFO("No PixelDistortions will be simulated.");
+ }
+
+ if (detStore()->retrieve(m_pixelid, "PixelID").isFailure()) {
+ ATH_MSG_FATAL("Could not get Pixel ID helper");
+ return StatusCode::FAILURE;
+ }
+
+ sc = service("DetectorStore", m_detStore);
+ if (!sc.isSuccess() || 0 == m_detStore) {
+ return msg(MSG:: ERROR) << "Could not find DetStore" << endmsg, StatusCode::FAILURE;
+ }
+ if (m_storeGate.retrieve().isFailure()) {
+ ATH_MSG_WARNING("Can not retrieve " << m_storeGate << ". Exiting.");
+ return StatusCode::FAILURE;
+ }
+
+ /* ME : booking of test histos is something to be hidden from production code !
+ sc = service("HistogramDataSvc", m_HistSvc, true);
+ // Define here the histograms;
+ m_h_Resx = m_HistSvc->book(m_foldername,"Resx","X shift (um)",400,-400.,400.);
+ m_h_Resy = m_HistSvc->book(m_foldername,"Resy","Y shift (um)",400,-400.,400.);
+ m_h_Locx = m_HistSvc->book(m_foldername,"Locx","X position (mm)",400,-20.,20.);
+ m_h_Locy = m_HistSvc->book(m_foldername,"Locy","Y position (mm)",800,-40.,40.);
+ m_h_ThetaTrack = m_HistSvc->book(m_foldername,"ThetaTrack","Theta Track (rad)",140,-7.,7);
+ m_h_PhiTrack = m_HistSvc->book(m_foldername,"PhiTrack","Phi Track (rad)",140,-7.,7);
+ */
+
+
+ m_dRMap = SG::ReadHandleKey<InDet::DRMap>(m_dRMapName);
+ ATH_CHECK( m_dRMap.initialize() );
+ // Include IBL calibration constants
+ //Moved to initialize to remove statics and prevent repitition
+
+ constexpr double phimin=-0.27, phimax=0.27;
+ for (int i=0; i<=nbinphi; i++) phix[i]=phimin+i*(phimax-phimin)/nbinphi;
+ constexpr double etacen[nbineta]={-0.,1.,1.55,1.9,2.15,2.35};
+ etax[0]=0.; etax[nbineta]=2.7;
+ for (int i=0; i<nbineta-1; i++) etax[i+1]=(etacen[i]+etacen[i+1])/2.;
+
+ ///UGLY!
+#include "IBL_calibration.h"
+ ///
+
+
+ return sc;
+}
+
+
+
+///////////////////////////////////////////////////////////////////
+// Finalize
+///////////////////////////////////////////////////////////////////
+
+StatusCode
+FTK_PixelClusterOnTrackTool::finalize() {
+ return AlgTool::finalize();
+}
+
+void
+FTK_PixelClusterOnTrackTool::FillFromDataBase() const{
+ if (m_fX.empty()) {
+ const CondAttrListCollection *atrlistcol = nullptr;
+ if (StatusCode::SUCCESS == m_detStore->retrieve(atrlistcol, "/PIXEL/PixelClustering/PixelCovCorr")) {
+ // loop over objects in collection
+ for (CondAttrListCollection::const_iterator citr = atrlistcol->begin(); citr != atrlistcol->end(); ++citr) {
+ std::vector<float> fx, fy, fb, fc, fd;
+ const coral::AttributeList &atrlist = citr->second;
+
+ fx.push_back(atrlist["fX1"].data<float>());
+ fx.push_back(atrlist["fX2"].data<float>());
+ fx.push_back(atrlist["fX3"].data<float>());
+ fx.push_back(atrlist["fX4"].data<float>());
+ fx.push_back(atrlist["fX5"].data<float>());
+ fx.push_back(atrlist["fX6"].data<float>());
+ fx.push_back(atrlist["fX7"].data<float>());
+ m_fX.emplace_back(std::move(fx));
+
+ fy.push_back(atrlist["fY1"].data<float>());
+ fy.push_back(atrlist["fY2"].data<float>());
+ fy.push_back(atrlist["fY3"].data<float>());
+ fy.push_back(atrlist["fY4"].data<float>());
+ fy.push_back(atrlist["fY5"].data<float>());
+ fy.push_back(atrlist["fY6"].data<float>());
+ fy.push_back(atrlist["fY7"].data<float>());
+ m_fY.emplace_back(std::move(fy));
+
+ fb.push_back(atrlist["fB1"].data<float>());
+ fb.push_back(atrlist["fB2"].data<float>());
+ fb.push_back(atrlist["fB3"].data<float>());
+ fb.push_back(atrlist["fB4"].data<float>());
+ fb.push_back(atrlist["fB5"].data<float>());
+ fb.push_back(atrlist["fB6"].data<float>());
+ fb.push_back(atrlist["fB7"].data<float>());
+ m_fB.emplace_back(std::move(fb));
+
+ fc.push_back(atrlist["fC1"].data<float>());
+ fc.push_back(atrlist["fC2"].data<float>());
+ fc.push_back(atrlist["fC3"].data<float>());
+ fc.push_back(atrlist["fC4"].data<float>());
+ fc.push_back(atrlist["fC5"].data<float>());
+ fc.push_back(atrlist["fC6"].data<float>());
+ fc.push_back(atrlist["fC7"].data<float>());
+ m_fC.emplace_back(std::move(fc));
+
+ fd.push_back(atrlist["fD1"].data<float>());
+ fd.push_back(atrlist["fD2"].data<float>());
+ fd.push_back(atrlist["fD3"].data<float>());
+ fd.push_back(atrlist["fD4"].data<float>());
+ fd.push_back(atrlist["fD5"].data<float>());
+ fd.push_back(atrlist["fD6"].data<float>());
+ fd.push_back(atrlist["fD7"].data<float>());
+ m_fD.emplace_back(std::move(fd));
+ }
+ }else {
+ msg(MSG::ERROR) << "Cannot find covariance corrections for key "
+ << "/PIXEL/PixelClustering/PixelCovCorr" << " - no correction " << endmsg;
+
+ }
+ }
+}
+
+
+///////////////////////////////////////////////////////////////////
+// Trk::SiClusterOnTrack production
+///////////////////////////////////////////////////////////////////
+
+
+const InDet::PixelClusterOnTrack *
+FTK_PixelClusterOnTrackTool::correct
+ (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar) const {
+
+
+ if (!m_applyNNcorrection ||
+ ((dynamic_cast<const InDetDD::SiDetectorElement *>(rio.detectorElement()))->isBlayer() && !m_NNIBLcorrection &&
+ !m_IBLAbsent)) {
+ return correctDefault(rio, trackPar);
+ }else {
+ if (m_errorStrategy == 0 || m_errorStrategy == 1) {
+ // version from Giacinto
+ if (m_noNNandBroadErrors) {
+ return 0;
+ }
+ // if we try broad errors, get Pixel Cluster to test if it is split
+ const InDet::PixelCluster *pix = dynamic_cast<const InDet::PixelCluster *>(&rio);
+ if (!pix) {
+ return 0;
+ }
+ if (pix->isSplit()) {
+ return correctNN(rio, trackPar);
+ } else {
+ return correctDefault(rio, trackPar);
+ }
+ } else {
+ return correctNN(rio, trackPar);
+ }
+ }
+}
+
+/** The correct method produces a PixelClusterOnTrack using the
+ * measured PixelCluster and the track prediction.
+ */
+const InDet::PixelClusterOnTrack *
+FTK_PixelClusterOnTrackTool::correctDefault
+ (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar) const {
+ // const InDet::SiCluster* SC = dynamic_cast<const InDet::SiCluster*> (&rio);
+ const InDet::PixelCluster *pix = nullptr;
+
+ if (!(pix = dynamic_cast<const InDet::PixelCluster *>(&rio))) {
+ return 0;
+ }
+
+ ATH_MSG_VERBOSE("Correct called with Error strategy " << m_errorStrategy);
+
+ // PixelClusterOnTrack production
+ //
+ Trk::LocalParameters locpar;
+ Amg::Vector3D glob(pix->globalPosition());
+
+
+ // Get pointer to detector element
+ const InDetDD::SiDetectorElement *element = pix->detectorElement();
+ if (!element) {
+ return 0;
+ }
+ bool blayer = element->isBlayer();
+ IdentifierHash iH = element->identifyHash();
+
+ double errphi = -1;
+ double erreta = -1;
+
+ if (pix->rdoList().empty()) {
+ ATH_MSG_WARNING("Pixel RDO-list size is 0, check integrity of pixel clusters! stop ROT creation.");
+ return NULL;
+ } else {
+
+ // get candidate track angle in module local frame
+ Amg::Vector3D my_track = trackPar.momentum();
+ Amg::Vector3D my_normal = element->normal();
+ Amg::Vector3D my_phiax = element->phiAxis();
+ Amg::Vector3D my_etaax = element->etaAxis();
+ float trkphicomp = my_track.dot(my_phiax);
+ float trketacomp = my_track.dot(my_etaax);
+ float trknormcomp = my_track.dot(my_normal);
+ double bowphi = atan2(trkphicomp, trknormcomp);
+ double boweta = atan2(trketacomp, trknormcomp);
+ float etatrack = trackPar.eta();
+
+ float tanl = element->getTanLorentzAnglePhi();
+ int readoutside = element->design().readoutSide();
+
+ // map the angles of inward-going tracks onto [-PI/2, PI/2]
+ if (bowphi > M_PI *0.5) {
+ bowphi -= M_PI;
+ }
+ if (bowphi < -M_PI *0.5) {
+ bowphi += M_PI;
+ }
+ // finally, subtract the Lorentz angle effect
+ // the readoutside term is needed because of a bug in old
+ // geometry versions (CSC-01-* and CSC-02-*)
+ double angle = atan(tan(bowphi) - readoutside * tanl);
+
+ // try to understand...
+ const Identifier element_id = element->identify();
+ int PixEtaModule = m_pixelid->eta_module(element_id);
+ int PixPhiModule = m_pixelid->phi_module(element_id);
+ double PixTrkPt = trackPar.pT();
+ double PixTrkEta = trackPar.eta();
+ ATH_MSG_VERBOSE("tanl = " << tanl << " readout side is " << readoutside <<
+ " module " << PixEtaModule << " " << PixPhiModule <<
+ " track pt, eta = " << PixTrkPt << " " << PixTrkEta <<
+ " track momentum phi, norm = " << trkphicomp << " " <<
+ trknormcomp << " bowphi = " << bowphi << " angle = " << angle);
+
+ float omegaphi = pix->omegax();
+ float omegaeta = pix->omegay();
+ double localphi = -9999.;
+ double localeta = -9999.;
+
+ const std::vector<Identifier> & rdos = pix->rdoList();
+ std::vector<Identifier>::const_iterator oneRDO = rdos.begin();
+ Identifier rId = *oneRDO;
+ const int row = m_pixelid->phi_index(rId);
+ const int col = m_pixelid->eta_index(rId);
+
+ InDetDD::SiLocalPosition centroid = pix->localPosition();
+ double shift = element->getLorentzCorrection();
+ int nrows = pix->width().colRow()[Trk::locX];
+ int ncol = pix->width().colRow()[Trk::locX];
+ double ang = 999.;
+ double eta = 999.;
+ ATH_MSG_VERBOSE(" row " << row << " col " << col << " nrows " << nrows <<
+ " ncol " << ncol << " locX " << centroid.xPhi() << " locY " << centroid.xEta());
+
+ // ATH_MSG_VERBOSE ( << "Position strategy = "
+ // << m_positionStrategy << "omegaphi = " << omegaphi )
+
+ // TOT interpolation for collision data
+ // Force IBL to use digital clustering and broad errors.
+ if (m_positionStrategy > 0 && omegaphi > -0.5 && omegaeta > -0.5) {
+ localphi = centroid.xPhi() + shift;
+ localeta = centroid.xEta();
+ // barrel
+ if (element->isBarrel()) {
+ ang = 180 * angle * M_1_PI; //M_1_PI in cmath, = 1/pi
+ double delta = 0.;
+ if (m_IBLAbsent || !blayer) {
+ delta = m_calibSvc->getBarrelDeltaX(nrows, ang);
+ } else { // special calibration for IBL
+ if (angle < phix[0] || angle > phix[nbinphi] || nrows != 2) {
+ delta = 0.;
+ }else {
+ int bin = -1;
+ while (angle > phix[bin + 1]) {
+ bin++;
+ }
+ if ((bin >= 0)and(bin < nbinphi)) {
+ delta = calphi[bin];
+ } else {
+ ATH_MSG_ERROR("bin out of range in line " << __LINE__ << " of FTK_PixelClusterOnTrackTool.cxx.");
+ }
+ }
+ if (m_calibSvc->includesIBLParams()) {
+ delta = m_calibSvc->getIBLDeltaX(nrows, ang);
+ }
+ }
+ localphi += delta * (omegaphi - 0.5);
+ ATH_MSG_VERBOSE(" shifting localphi by delta * (omegaphi - 0.5) = " << delta << " * " << (omegaphi - 0.5));
+ // settle the sign/pi periodicity issues
+ double thetaloc = -999.;
+ if (boweta > -0.5 * M_PI && boweta < M_PI / 2.) { //M_PI_2 in cmath
+ thetaloc = M_PI_2 - boweta;
+ }else if (boweta > M_PI_2 && boweta < M_PI) {
+ thetaloc = 1.5 * M_PI - boweta;
+ } else { // 3rd quadrant
+ thetaloc = -M_PI_2 - boweta;
+ }
+ double etaloc = -1 * log(tan(thetaloc * 0.5));
+ if (m_IBLAbsent || !blayer) {
+ delta = m_calibSvc->getBarrelDeltaY(ncol, etaloc);
+ ATH_MSG_VERBOSE(" Barrel Pixel: shifting localeta by delta * (omegaeta - 0.5) = " << delta << " * " << (omegaeta - 0.5));
+ } else { // special calibration for IBL
+ etaloc = fabs(etaloc);
+ if (etaloc < etax[0] || etaloc > etax[nbineta]) {
+ delta = 0.;
+ } else {
+ int bin = -1;
+ while (etaloc > etax[bin + 1]) {
+ bin++;
+ }
+ if ((bin >= 0)and(bin < nbineta)) {
+ if (ncol == bin) {
+ delta = caleta[bin][0];
+ } else if (ncol == bin + 1) {
+ delta = caleta[bin][1];
+ } else if (ncol == bin + 2) {
+ delta = caleta[bin][2];
+ } else {
+ delta = 0.;
+ }
+ } else {// bin out of range of array indices
+ ATH_MSG_ERROR("bin out of range in line " << __LINE__ << " of FTK_PixelClusterOnTrackTool.cxx.");
+ }
+
+ }
+ if (m_calibSvc->includesIBLParams()) {
+ delta = m_calibSvc->getIBLDeltaY(ncol, etaloc);
+ }
+ ATH_MSG_VERBOSE(" IBL: shifting localeta by delta * (omegaeta - 0.5) = " << delta << " * " << (omegaeta - 0.5));
+ }
+ localeta += delta * (omegaeta - 0.5);
+ }else {
+ // collision endcap data
+ if (m_positionStrategy == 1) {
+ double deltax = m_calibSvc->getEndcapDeltaX();
+ double deltay = m_calibSvc->getEndcapDeltaY();
+ localphi += deltax * (omegaphi - 0.5);
+ localeta += deltay * (omegaeta - 0.5);
+ ATH_MSG_VERBOSE(" endcap pixel shifting localphi by deltax * (omegaphi - 0.5) = " << deltax << " * " << (omegaphi - 0.5));
+ ATH_MSG_VERBOSE(" endcap pixel shifting localeta by deltay * (omegaeta - 0.5) = " << deltay << " * " << (omegaeta - 0.5));
+ }
+ // SR1 cosmics endcap data
+ // some parametrization dependent on track angle
+ // would be better here
+ else if (m_positionStrategy == 20) {
+ double deltax = 35 * micrometer;
+ double deltay = 35 * micrometer;
+ localphi += deltax * (omegaphi - 0.5);
+ localeta += deltay * (omegaeta - 0.5);
+ }
+ }
+ }
+// digital
+ else {
+ localphi = centroid.xPhi() + shift;
+ localeta = centroid.xEta();
+ }
+
+ InDet::SiWidth width = pix->width();
+
+ // Error strategies
+
+ // For very shallow tracks the cluster can easily break as
+ // the average charge per pixel is of the order of the threshold
+ // In this case, an error equal to the geometrical projection
+ // of the track path in silicon onto the module surface seems
+ // appropriate
+ if (fabs(angle) > 1) {
+ errphi = 250 * micrometer * tan(fabs(angle)) * TOPHAT_SIGMA;
+ erreta = width.z() > 250 * micrometer * tan(fabs(boweta)) ?
+ width.z() * TOPHAT_SIGMA : 250 * micrometer * tan(fabs(boweta)) * TOPHAT_SIGMA;
+ ATH_MSG_VERBOSE("Shallow track with tanl = " << tanl << " bowphi = " <<
+ bowphi << " angle = " << angle << " width.z = " << width.z() <<
+ " errphi = " << errphi << " erreta = " << erreta);
+ }else if (m_errorStrategy == 0) {
+ errphi = width.phiR() * TOPHAT_SIGMA;
+ erreta = width.z() * TOPHAT_SIGMA;
+ }else if (m_errorStrategy == 1) {
+ errphi = (width.phiR() / nrows) * TOPHAT_SIGMA;
+ erreta = (width.z() / ncol) * TOPHAT_SIGMA;
+ }else if (m_errorStrategy == 2) {
+ if (element->isBarrel()) {
+ if (m_IBLAbsent || !blayer) {
+ errphi = m_calibSvc->getBarrelNewErrorPhi(ang, nrows);
+ } else { // special calibration for IBL
+ if (angle < phix[0] || angle > phix[nbinphi]) {
+ errphi = width.phiR() * TOPHAT_SIGMA;
+ } else {
+ int bin = -1;// cannot be used as array index, which will happen if angle<phix[bin+1]
+ while (angle > phix[bin + 1]) {
+ bin++;
+ }
+ if ((bin >= 0)and(bin < nbinphi)) {
+ if (nrows == 1) {
+ errphi = calerrphi[bin][0];
+ } else if (nrows == 2) {
+ errphi = calerrphi[bin][1];
+ } else {
+ errphi = calerrphi[bin][2];
+ }
+ } else {
+ ATH_MSG_ERROR("bin out of range in line " << __LINE__ << " of FTK_PixelClusterOnTrackTool.cxx.");
+ }
+ }
+ }
+
+ if (m_rel13like) {
+ erreta = m_calibSvc->getBarrelErrorEta(eta, ncol, nrows);
+ }else if (m_IBLAbsent || !blayer) {
+ erreta = m_calibSvc->getBarrelNewErrorEta(fabs(etatrack), nrows, ncol);
+ } else { // special calibration for IBL
+ double etaloc = fabs(etatrack);
+ if (etaloc < etax[0] || etaloc > etax[nbineta]) {
+ erreta = width.z() * TOPHAT_SIGMA;
+ } else {
+ int bin = 0;
+ while (etaloc > etax[bin + 1]) {
+ ++bin;
+ }
+ if (bin >= nbineta) {
+ ATH_MSG_ERROR("bin out of range in line " << __LINE__ << " of FTK_PixelClusterOnTrackTool.cxx.");
+ } else {
+ if (ncol == bin) {
+ erreta = calerreta[bin][0];
+ } else if (ncol == bin + 1) {
+ erreta = calerreta[bin][1];
+ } else if (ncol == bin + 2) {
+ erreta = calerreta[bin][2];
+ } else {
+ erreta = width.z() * TOPHAT_SIGMA;
+ }
+ }
+ }
+ }
+ }else {
+ errphi = m_calibSvc->getEndCapErrorPhi(ncol, nrows);
+ erreta = m_calibSvc->getEndCapErrorEta(ncol, nrows);
+ }
+ if (errphi > erreta) {
+ erreta = width.z() * TOPHAT_SIGMA;
+ }
+ }
+
+ Amg::Vector2D locpos = Amg::Vector2D(localphi, localeta);
+ if (element->isBarrel() && !m_disableDistortions) {
+ correctBow(element->identify(), locpos, bowphi, boweta);
+ }
+
+
+ locpar = Trk::LocalParameters(locpos);
+ centroid = InDetDD::SiLocalPosition(localeta, localphi, 0.);
+ glob = element->globalPosition(centroid);
+ }
+
+ // Error matrix production
+
+ Amg::MatrixX cov = pix->localCovariance();
+
+ // corrected phi error
+ if (errphi > 0) {
+ cov(0, 0) = errphi * errphi;
+ }
+ if (erreta > 0) {
+ cov(1, 1) = erreta * erreta;
+ }
+
+ ATH_MSG_VERBOSE(" errphi = " << errphi << " erreta = " << erreta);
+
+ // create new copy of error matrix
+ if (m_scalePixelCov) {
+ Amg::MatrixX *newCov = m_errorScalingTool->createScaledPixelCovariance(cov, element->identify());
+ if (!newCov) {
+ ATH_MSG_WARNING("Failed to create scaled error for Pixel");
+ return 0;
+ }
+ cov = *newCov;
+ delete newCov;
+ }
+ bool isbroad = (m_errorStrategy == 0) ? true : false;
+
+ ATH_MSG_VERBOSE("returning corrected cluster at locX " << locpar.get(Trk::locX) << " locY " << locpar.get(Trk::locY));
+
+ return new InDet::PixelClusterOnTrack(pix, locpar, cov, iH, glob, pix->gangedPixel(), isbroad);
+}
+
+void
+FTK_PixelClusterOnTrackTool::correctBow(const Identifier &id, Amg::Vector2D &localpos, const double phi,
+ const double theta) const {
+ Amg::Vector3D dir(tan(phi), tan(theta), 1.);
+ Amg::Vector2D newpos =
+ m_pixDistoTool->correctReconstruction(id, localpos, dir);
+
+ localpos = newpos;
+ return;
+}
+
+const InDet::PixelClusterOnTrack *
+FTK_PixelClusterOnTrackTool::correct
+ (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar,
+ const InDet::PixelClusterStrategy strategy) const {
+ int initial_errorStrategy;
+ const InDet::PixelClusterOnTrack *newROT;
+
+ switch (strategy) {
+ case InDet::PIXELCLUSTER_OUTLIER: // if cluster is outlier, increase errors
+ case InDet::PIXELCLUSTER_SHARED:
+ initial_errorStrategy = m_errorStrategy;
+ if (!m_rel13like) {
+ m_errorStrategy = 0; // error as size of cluster /sqrt(12)
+ }
+ newROT = correct(rio, trackPar);
+ m_errorStrategy = initial_errorStrategy;
+ return newROT;
+
+ default:
+ return correct(rio, trackPar);
+ }
+}
+
+// GP: NEW correct() method in case of NN based calibration */
+const InDet::PixelClusterOnTrack *
+FTK_PixelClusterOnTrackTool::correctNN
+ (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar) const {
+
+ return correct(rio, trackPar);
+}
+
+bool
+FTK_PixelClusterOnTrackTool::getErrorsDefaultAmbi(const InDet::PixelCluster *pixelPrepCluster,
+ const Trk::TrackParameters &trackPar,
+ Amg::Vector2D &finalposition,
+ Amg::MatrixX &finalerrormatrix) const {
+ std::vector<Amg::Vector2D> vectorOfPositions;
+ int numberOfSubclusters = 1;
+
+ vectorOfPositions.push_back(pixelPrepCluster->localPosition());
+
+
+
+ // A.S. hack for the moment: isSplit is also set for modified (non-split) 1-pixel-clusters ... is this needed anyway ?
+ //
+ // if (pixelPrepCluster->isSplit() && numberOfSubclusters<2)
+ // {
+ // msg(MSG::WARNING) << " Cluster is split but no further clusters found in split cluster map." << endmsg;
+ // }
+
+ // now you have numberOfSubclusters and the vectorOfPositions (Amg::Vector2D)
+
+ const Trk::AtaPlane *parameters = dynamic_cast<const Trk::AtaPlane *>(&trackPar);
+ if (parameters == 0) {
+ msg(MSG::WARNING) << "Parameters are not at a plane ! Aborting cluster correction... " << endmsg;
+ return false;
+ }
+
+ std::vector<Amg::Vector2D> allLocalPositions;
+ std::vector<Amg::MatrixX> allErrorMatrix;
+
+ allLocalPositions =
+ m_NnClusterizationFactory->estimatePositions(*pixelPrepCluster,
+ parameters->associatedSurface(),
+ *parameters,
+ allErrorMatrix,
+ numberOfSubclusters);
+
+ if (allLocalPositions.size() == 0) {
+ if (msgLvl(MSG::DEBUG)) {
+ msg(MSG::DEBUG) << " Cluster cannot be treated by NN. Giving back to default clusterization " << endmsg;
+ }
+ return false;
+ }
+
+ if (allLocalPositions.size() != size_t(numberOfSubclusters)) {
+ msg(MSG::WARNING) << "Returned position vector size " << allLocalPositions.size() <<
+ " not according to expected number of subclusters: " << numberOfSubclusters << ". Abort cluster correction..." <<
+ endmsg;
+ return false;
+ }
+
+
+ // GP: now the not so nice part of matching the new result with the old one...
+ // Takes the error into account to improve the matching
+
+ if (numberOfSubclusters == 1) {
+ finalposition = allLocalPositions[0];
+ finalerrormatrix = allErrorMatrix[0];
+ }
+
+ else if (numberOfSubclusters == 2) {
+ double distancesq1 =
+ square(vectorOfPositions[0][0] - allLocalPositions[0][0]) / allErrorMatrix[0](0, 0) +
+ square(vectorOfPositions[1][0] - allLocalPositions[1][0]) / allErrorMatrix[1](0, 0) +
+ square(vectorOfPositions[0][1] - allLocalPositions[0][1]) / allErrorMatrix[0](1, 1) +
+ square(vectorOfPositions[1][1] - allLocalPositions[1][1]) / allErrorMatrix[1](1, 1);
+
+ double distancesq2 =
+ square(vectorOfPositions[1][0] - allLocalPositions[0][0]) / allErrorMatrix[0](0, 0) +
+ square(vectorOfPositions[0][0] - allLocalPositions[1][0]) / allErrorMatrix[1](0, 0) +
+ square(vectorOfPositions[1][1] - allLocalPositions[0][1]) / allErrorMatrix[0](1, 1) +
+ square(vectorOfPositions[0][1] - allLocalPositions[1][1]) / allErrorMatrix[1](1, 1);
+
+ if (msgLvl(MSG::DEBUG)) {
+ msg(MSG::DEBUG) << " Old pix (1) x: " << vectorOfPositions[0][0] << " y: " << vectorOfPositions[0][1] << endmsg;
+ msg(MSG::DEBUG) << " Old pix (2) x: " << vectorOfPositions[1][0] << " y: " << vectorOfPositions[1][1] << endmsg;
+ msg(MSG::DEBUG) << " Pix (1) x: " << allLocalPositions[0][0] << " +/- " << sqrt(allErrorMatrix[0](0, 0))
+ << " y: " << allLocalPositions[0][1] << " +/- " << sqrt(allErrorMatrix[0](1, 1)) << endmsg;
+ msg(MSG::DEBUG) << " Pix (2) x: " << allLocalPositions[1][0] << " +/- " << sqrt(allErrorMatrix[1](0, 0))
+ << " y: " << allLocalPositions[1][1] << " +/- " << sqrt(allErrorMatrix[1](1, 1)) << endmsg;
+ msg(MSG::DEBUG) << " Old (1) new (1) dist: " << sqrt(distancesq1) << " Old (1) new (2) " << sqrt(distancesq2) <<
+ endmsg;
+ }
+
+ if (distancesq1 < distancesq2) {
+ finalposition = allLocalPositions[0];
+ finalerrormatrix = allErrorMatrix[0];
+ }else {
+ finalposition = allLocalPositions[1];
+ finalerrormatrix = allErrorMatrix[1];
+ }
+ }
+
+
+ else if (numberOfSubclusters == 3) {
+ double distances[6];
+
+ distances[0] = distance(vectorOfPositions, allLocalPositions, allErrorMatrix, 0, 1, 2);
+ distances[1] = distance(vectorOfPositions, allLocalPositions, allErrorMatrix, 0, 2, 1);
+ distances[2] = distance(vectorOfPositions, allLocalPositions, allErrorMatrix, 1, 0, 2);
+ distances[3] = distance(vectorOfPositions, allLocalPositions, allErrorMatrix, 1, 2, 0);
+ distances[4] = distance(vectorOfPositions, allLocalPositions, allErrorMatrix, 2, 0, 1);
+ distances[5] = distance(vectorOfPositions, allLocalPositions, allErrorMatrix, 2, 1, 0);
+
+ int smallestDistanceIndex = -10;
+ double minDistance = 1e10;
+
+ for (int i = 0; i < 6; i++) {
+ ATH_MSG_VERBOSE(" distance n.: " << i << " distance is: " << distances[i]);
+
+ if (distances[i] < minDistance) {
+ minDistance = distances[i];
+ smallestDistanceIndex = i;
+ }
+ }
+
+ ATH_MSG_DEBUG(" The minimum distance is : " << minDistance << " for index: " << smallestDistanceIndex);
+
+ if (smallestDistanceIndex == 0 || smallestDistanceIndex == 1) {
+ finalposition = allLocalPositions[0];
+ finalerrormatrix = allErrorMatrix[0];
+ }
+ if (smallestDistanceIndex == 2 || smallestDistanceIndex == 4) {
+ finalposition = allLocalPositions[1];
+ finalerrormatrix = allErrorMatrix[1];
+ }
+ if (smallestDistanceIndex == 3 || smallestDistanceIndex == 5) {
+ finalposition = allLocalPositions[2];
+ finalerrormatrix = allErrorMatrix[2];
+ }
+ }
+ return true;
+}
+
+bool
+FTK_PixelClusterOnTrackTool::getErrorsTIDE_Ambi(const InDet::PixelCluster *pixelPrepCluster,
+ const Trk::TrackParameters &trackPar,
+ Amg::Vector2D &finalposition,
+ Amg::MatrixX &finalerrormatrix) const {
+ std::vector<Amg::Vector2D> vectorOfPositions;
+
+ int numberOfSubclusters = 1;
+ // now you have numberOfSubclusters and the vectorOfPositions (Amg::Vector2D)
+
+ const Trk::AtaPlane *parameters = dynamic_cast<const Trk::AtaPlane *>(&trackPar);
+ if (parameters == 0) {
+ ATH_MSG_WARNING("Parameters are not at a plane ! Aborting cluster correction... ");
+ return false;
+ }
+
+ std::vector<Amg::Vector2D> allLocalPositions;
+ std::vector<Amg::MatrixX> allErrorMatrix;
+ allLocalPositions =
+ m_NnClusterizationFactory->estimatePositions(*pixelPrepCluster,
+ parameters->associatedSurface(),
+ *parameters,
+ allErrorMatrix,
+ numberOfSubclusters);
+
+ if (allLocalPositions.empty()) {
+ ATH_MSG_DEBUG(
+ " Cluster cannot be treated by NN. Giving back to default clusterization, too big: " <<
+ pixelPrepCluster->tooBigToBeSplit());
+ return false;
+ }
+
+ if (allLocalPositions.size() != size_t(numberOfSubclusters)) {
+ ATH_MSG_WARNING(
+ "Returned position vector size " << allLocalPositions.size() << " not according to expected number of subclusters: " << numberOfSubclusters <<
+ ". Abort cluster correction...");
+ return false;
+ }
+
+ bool correctdR = false;
+
+ std::array<float, 3> correctiondR = {0., 0., 0.};
+
+ if(m_applydRcorrection){
+ SG::ReadHandle<InDet::DRMap> mapHandle(m_dRMap);
+ if(!mapHandle.isValid()){
+ ATH_MSG_ERROR("Error retreving " << m_dRMap.key());
+ return false;
+ }
+
+ if(!mapHandle->empty() && pixelPrepCluster->isSplit() ){
+
+ const InDetDD::SiDetectorElement *det = pixelPrepCluster->detectorElement();
+
+ if (det->isBarrel()) {
+ int pixelLayer = m_pixelid->layer_disk(det->identify());
+
+ float correctiondRX = 1.0; // correction factor to be applied to covariance(!)
+ float correctiondRZ = 1.0; // correction factor to be applied to covariance(!)
+
+ const InDetDD::PixelModuleDesign *design(dynamic_cast<const InDetDD::PixelModuleDesign *>(&det->design()));
+ if (! design){
+ ATH_MSG_WARNING("Design cannot be cast to a pixel design");
+ return false;
+ }
+ const float pitchX = design->phiPitch();
+ const float pitchZ = design->etaPitch();
+
+ ATH_MSG_DEBUG("FTK_PixelClusterOnTrackTool: Trying to apply dR correction for pixelLayer " << pixelLayer);
+
+ InDet::DRMap::const_iterator it = mapHandle->find(pixelPrepCluster);
+ if (it != mapHandle->end()) {
+ correctdR = true;
+ float mindX = it->second.first;
+ float mindZ = it->second.second;
+ float mindXOverPitch = mindX / pitchX;
+ float mindZOverPitch = mindZ / pitchZ;
+
+ float dR = sqrt(mindX * mindX + mindZ * mindZ);
+
+ ATH_MSG_DEBUG(" ---- Min dX ---- " << mindX << " mindZ " << mindZ);
+
+ if (dR > 0. && dR < 2.0) {
+ correctiondRX = splineIBLPullX(mindXOverPitch, pixelLayer);
+ correctiondRZ = splineIBLPullX(mindZOverPitch, pixelLayer);
+
+ if (correctiondRX > 3.75) {
+ correctiondRX = 3.75;
+ }
+ if (correctiondRZ > 3.75) {
+ correctiondRZ = 3.75;
+ }
+ }else {
+ correctiondRX = 2.35;
+ correctiondRZ = 1.70;
+ }
+
+ ATH_MSG_DEBUG(
+ "FTK_PixelClusterOnTrackTool: Correction factor calculation for distX/pitch= " << mindXOverPitch << " -> " << correctiondRX << " distZ/pitch= " << mindZOverPitch << " -> " <<
+ correctiondRZ);
+ }else {
+ ATH_MSG_WARNING("Split Pixel cluster not found in dRmap! " << pixelPrepCluster);
+ }
+
+ ATH_MSG_DEBUG(
+ " ++++ Hit Error ++++ Layer " << pixelLayer << " Correction Flag " << m_applydRcorrection << " Split " << pixelPrepCluster->isSplit() << " Scaling " << correctiondRX << " " <<
+ correctiondRZ);
+
+ correctiondRX *= correctiondRX;
+ correctiondRZ *= correctiondRZ;
+
+ correctiondR[0] = correctiondRX * correctiondRX;
+ correctiondR[1] = correctiondRZ * correctiondRZ;
+ correctiondR[2] = correctiondRX * correctiondRZ;
+ }
+ }
+ }
+ // AKM: now the not so nice part find the best match position option
+ // Takes the error into account to scale the importance of the measurement
+
+ if (numberOfSubclusters == 1) {
+ finalposition = allLocalPositions[0];
+
+ if (correctdR) {
+ allErrorMatrix[0](0, 0) *= correctiondR[0];
+ allErrorMatrix[0](1, 1) *= correctiondR[1];
+ allErrorMatrix[0](0, 1) *= correctiondR[2];
+ allErrorMatrix[0](1, 0) *= correctiondR[2];
+ ATH_MSG_DEBUG(
+ " ++++ Hit Error ++++ " << numberOfSubclusters << " Split " << pixelPrepCluster->isSplit() << " Scaling " <<
+ sqrt(correctiondR[0]) << " " << sqrt(correctiondR[1]) << " Rescaled ErrX " << sqrt(allErrorMatrix[0](0,
+ 0)) << " Rescaled ErrZ " <<
+ sqrt(allErrorMatrix[0](1, 1)));
+ }
+
+ finalerrormatrix = allErrorMatrix[0];
+ return true;
+ }
+
+ // Get the track parameters local position
+ const Amg::Vector2D localpos = trackPar.localPosition();
+ // Use the track parameters cov to weight distcances
+ Amg::Vector2D localerr(0.01, 0.05);
+ if (trackPar.covariance()) {
+ localerr = Amg::Vector2D(sqrt((*trackPar.covariance())(0, 0)), sqrt((*trackPar.covariance())(1, 1)));
+ }
+
+ double minDistance(1e300);
+ int index(0);
+
+ for (unsigned int i(0); i < allLocalPositions.size(); ++i) {
+ double distance =
+ square(localpos[0] - allLocalPositions[i][0]) / localerr[0]
+ + square(localpos[1] - allLocalPositions[i][1]) / localerr[1];
+
+ if (distance < minDistance) {
+ index = i;
+ minDistance = distance;
+ }
+ }
+
+ finalposition = allLocalPositions[index];
+
+ if (correctdR) {
+ allErrorMatrix[index](0, 0) *= correctiondR[0];
+ allErrorMatrix[index](1, 1) *= correctiondR[1];
+ allErrorMatrix[index](0, 1) *= correctiondR[2];
+ allErrorMatrix[index](1, 0) *= correctiondR[2];
+ ATH_MSG_DEBUG(
+ " ++++ Hit Error ++++ " << numberOfSubclusters << " Split " << pixelPrepCluster->isSplit() << " Scaling " <<
+ sqrt(correctiondR[0]) << " " << sqrt(correctiondR[1]) << " Rescaled ErrX " << sqrt(allErrorMatrix[index](0,
+ 0)) << " Rescaled ErrZ " <<
+ sqrt(allErrorMatrix[index](1, 1)));
+ }
+
+ finalerrormatrix = allErrorMatrix[index];
+ return true;
+}
+
+double
+FTK_PixelClusterOnTrackTool::splineIBLPullX(float x, int layer) const {
+ FillFromDataBase();
+ const int fNp = m_fY[layer].size();
+ const int fKstep = 0;
+ const double fDelta = -1;
+ const auto result = std::minmax_element(std::begin(m_fX[layer]), std::end(m_fX[layer]));
+ const double fXmin = *(result.first);
+ const double fXmax = *(result.second);
+ //const double fXmin = *std::min_element(std::begin(m_fX[layer]), std::end(m_fX[layer]));
+ //const double fXmax = *std::max_element(std::begin(m_fX[layer]), std::end(m_fX[layer]));
+
+ int klow = 0;
+
+ if (x <= fXmin) {
+ klow = 0;
+ } else if (x >= fXmax) {
+ klow = fNp - 1;
+ } else {
+ if (fKstep) {
+ // Equidistant knots, use histogramming
+ klow = int((x - fXmin) / fDelta);
+ if (klow < fNp - 1) {
+ klow = fNp - 1;
+ }
+ } else {
+ int khig = fNp - 1, khalf;
+ // Non equidistant knots, binary search
+ while (khig - klow > 1) {
+ if (x > m_fX[layer][khalf = (klow + khig) / 2]) {
+ klow = khalf;
+ } else {
+ khig = khalf;
+ }
+ }
+ }
+ }
+ // Evaluate now
+ double dx = x - m_fX[layer][klow];
+ return(m_fY[layer][klow] + dx * (m_fB[layer][klow] + dx * (m_fC[layer][klow] + dx * m_fD[layer][klow])));
+}
diff --git a/Trigger/TrigFTK/FTK_RecTools/src/FTK_SCTClusterOnTrackTool.cxx b/Trigger/TrigFTK/FTK_RecTools/src/FTK_SCTClusterOnTrackTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f98a727a8359e64d766ee63c992d8aa059e5aeec
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/src/FTK_SCTClusterOnTrackTool.cxx
@@ -0,0 +1,411 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// Implementation file for class FTK_SCTClusterOnTrackTool
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+// AlgTool used for SCT_ClusterOnTrack object production
+///////////////////////////////////////////////////////////////////
+// started 21/04/2004 I.Gavrilenko -- see ChangeLog
+///////////////////////////////////////////////////////////////////
+
+#include "FTK_RecTools/FTK_SCTClusterOnTrackTool.h"
+
+#include "TrkSurfaces/RectangleBounds.h"
+#include "TrkSurfaces/TrapezoidBounds.h"
+
+using CLHEP::micrometer;
+using CLHEP::deg;
+
+///////////////////////////////////////////////////////////////////
+// Constructor
+///////////////////////////////////////////////////////////////////
+
+FTK_SCTClusterOnTrackTool::FTK_SCTClusterOnTrackTool
+ (const std::string &t, const std::string &n, const IInterface *p) :
+ AthAlgTool(t, n, p),
+ m_errorScalingTool("Trk::RIO_OnTrackErrorScalingTool/RIO_OnTrackErrorScalingTool", this),
+ m_distortionsTool("SCT_DistortionsTool", this),
+ m_scaleSctCov(false),
+ m_option_make2dimBarrelClusters(false),
+ m_doDistortions(false),
+ m_option_errorStrategy(-1),
+ m_option_correctionStrategy(-1) {
+ // declareInterface<FTK_SCTClusterOnTrackTool>(this);
+ declareInterface<IRIO_OnTrackCreator>(this);
+
+ declareProperty("ErrorScalingTool", m_errorScalingTool,
+ "The toolhandle for central error scaling");
+ declareProperty("MakeTwoDimBarrelClusters", m_option_make2dimBarrelClusters,
+ "flag if strip length should be part of the measurement");
+ declareProperty("ErrorStrategy", m_option_errorStrategy,
+ "if ErrorStrategy < 0, keep previous errors else recompute");
+ declareProperty("CorrectionStrategy", m_option_correctionStrategy,
+ "if CorrectionStrategy >= 0, apply a correction to the cluster position");
+ declareProperty("doDistortions", m_doDistortions,
+ "Simulation of module distortions");
+ declareProperty("SCTDistortionsTool", m_distortionsTool,
+ "Tool to retrieve SCT distortions");
+}
+
+///////////////////////////////////////////////////////////////////
+// Destructor
+///////////////////////////////////////////////////////////////////
+
+FTK_SCTClusterOnTrackTool::~FTK_SCTClusterOnTrackTool() {
+}
+
+///////////////////////////////////////////////////////////////////
+// Initialisation
+///////////////////////////////////////////////////////////////////
+
+StatusCode
+FTK_SCTClusterOnTrackTool::initialize() {
+ StatusCode sc = AlgTool::initialize();
+
+ msg(MSG::INFO) << "A strategy to ";
+ switch (m_option_errorStrategy) {
+ case -1: msg(MSG::INFO) << "keep the PRD errors";
+ break;
+
+ case 0: msg(MSG::INFO) << "apply simple pitch errors";
+ break;
+
+ case 1: msg(MSG::INFO) << "assign tuned SCT errors";
+ break;
+
+ case 2: msg(MSG::INFO) << "assign tuned, angle-dependent SCT errors";
+ break;
+
+ default: msg(MSG::INFO) << " -- NO, UNKNOWN. Pls check jobOptions!";
+ break;
+ }
+ msg(MSG::INFO) << " will be applied during SCT_ClusterOnTrack making" << endmsg;
+ if (m_option_correctionStrategy == 0) {
+ msg(MSG::INFO) << "SCT cluster positions will be corrected" << endmsg;
+ }
+
+ if (m_errorScalingTool.retrieve().isFailure()) {
+ msg(MSG::FATAL) << "Failed to retrieve tool " << m_errorScalingTool << endmsg;
+ return StatusCode::FAILURE;
+ } else {
+ msg(MSG::INFO) << "Retrieved tool " << m_errorScalingTool << endmsg;
+ m_scaleSctCov = m_errorScalingTool->needToScaleSct();
+ if (m_scaleSctCov) {
+ msg(MSG::DEBUG) << "Detected need for scaling SCT errors." << endmsg;
+ }
+ }
+
+ // Get ISCT_ModuleDistortionsTool
+ if (m_distortionsTool.retrieve().isFailure()) {
+ msg(MSG::FATAL) << "Could not retrieve distortions tool: " << m_distortionsTool.name() << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ return sc;
+}
+
+///////////////////////////////////////////////////////////////////
+// Finalize
+///////////////////////////////////////////////////////////////////
+
+StatusCode
+FTK_SCTClusterOnTrackTool::finalize() {
+ StatusCode sc = AlgTool::finalize();
+
+ return sc;
+}
+
+///////////////////////////////////////////////////////////////////
+// Trk::SCT_ClusterOnTrack production
+///////////////////////////////////////////////////////////////////
+
+const InDet::SCT_ClusterOnTrack *
+FTK_SCTClusterOnTrackTool::correct
+ (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar) const {
+ const InDet::SCT_Cluster *SC = 0;
+
+ if (!(SC = dynamic_cast<const InDet::SCT_Cluster *> (&rio))) {
+ msg(MSG::VERBOSE) << "not a InDet::SCT_Cluster, return null" <<endmsg;
+ return 0;
+ }
+
+ const InDet::SiWidth width = SC->width();
+ const Amg::Vector2D &colRow = width.colRow();
+
+ // Get pointer to detector element
+ //
+ const InDetDD::SiDetectorElement *EL = SC->detectorElement();
+ if (!EL) {
+ msg(MSG::VERBOSE) << "can't get SiDetectorElement, return null" <<endmsg;
+ return 0;
+ }
+ IdentifierHash iH = EL->identifyHash();
+
+ // Get local position of track
+ //
+ Amg::Vector2D loct = trackPar.localPosition();
+
+ // Find phi angle of track relative to Lorentz drift direction, if needed
+ //
+ double dphi(0.);
+ double sinAlpha = EL->sinStereoLocal(SC->localPosition());
+ double cosAlpha = sqrt(1 - sinAlpha * sinAlpha);
+ if (m_option_errorStrategy == 2 || m_option_correctionStrategy == 0) {
+ double pNormal = trackPar.momentum().dot(EL->normal());
+ double pPhi = trackPar.momentum().dot(Amg::AngleAxis3D(asin(-sinAlpha), Amg::Vector3D::UnitZ()) * EL->phiAxis());
+ dphi = atan(pPhi / pNormal) - atan(EL->getTanLorentzAnglePhi());
+ }
+
+ // SCT_ClusterOnTrack production
+ Amg::MatrixX oldcov = SC->localCovariance();
+ // Local position and error matrix production
+ //
+ // let's start to re-compute cluster error if errorStrategy >=0
+ if (m_option_errorStrategy > -1) {
+ Amg::MatrixX mat(2, 2);
+ mat.setZero();
+ switch (m_option_errorStrategy) {
+ case 0:
+ mat(0, 0) = pow(width.phiR(), 2) / 12;
+ mat(1, 1) = pow(width.z(), 2) / 12;
+ break;
+
+ case 1:
+ if (colRow.x() == 1) {
+ mat(0, 0) = pow(1.05 * width.phiR(), 2) / 12;
+ } else if (colRow.x() == 2) {
+ mat(0, 0) = pow(0.27 * width.phiR(), 2) / 12;
+ } else {
+ mat(0, 0) = pow(width.phiR(), 2) / 12;
+ }
+ mat(1, 1) = pow(width.z() / colRow.y(), 2) / 12;
+ break;
+
+ case 2:
+ mat(0, 0) = pow(getError(dphi, int(colRow.x())) * (EL->phiPitch() / 0.080), 2);
+ mat(1, 1) = pow(width.z() / colRow.y(), 2) / 12;
+ break;
+
+ default:
+ // don't do anything....
+ break;
+ }
+ // rotation for endcap SCT
+ if (EL->design().shape() == InDetDD::Trapezoid) {
+ double sn = EL->sinStereoLocal(SC->localPosition());
+ double sn2 = sn * sn;
+ double cs2 = 1. - sn2;
+ double w = EL->phiPitch(SC->localPosition()) / EL->phiPitch();
+ double v0 = mat(0, 0) * w * w;
+ double v1 = mat(1, 1);
+ mat(0, 0) = (cs2 * v0 + sn2 * v1);
+ mat(1, 0) = (sn * sqrt(cs2) * (v0 - v1));
+ mat(0, 1) = mat(1, 0);
+ mat(1, 1) = (sn2 * v0 + cs2 * v1);
+ }
+ oldcov = mat;
+ }
+
+ Amg::MatrixX cov(oldcov);
+ Trk::LocalParameters locpar;
+
+ if (EL->design().shape() != InDetDD::Trapezoid) { // barrel
+ Trk::DefinedParameter lpos1dim(SC->localPosition().x(), Trk::locX);
+ locpar = (m_option_make2dimBarrelClusters) ?
+ Trk::LocalParameters(SC->localPosition()) :// PRDformation does 2-dim
+ Trk::LocalParameters(lpos1dim);
+ if (!m_option_make2dimBarrelClusters) {
+ cov = Amg::MatrixX(1, 1);
+ cov(0, 0) = oldcov(0, 0);
+ }
+
+ if (m_scaleSctCov) {
+ Amg::MatrixX *newCov = m_errorScalingTool->createScaledSctCovariance(cov, false, 0.0);
+ if (!newCov) {
+ ATH_MSG_WARNING("Failed to create scaled error for SCT");
+ return 0;
+ }
+ cov = *newCov;
+ delete newCov;
+ }
+ }else { // endcap
+ locpar = Trk::LocalParameters(SC->localPosition());
+ if (m_scaleSctCov) {
+ Amg::MatrixX *newCov = m_errorScalingTool->createScaledSctCovariance(cov, true,
+ EL->sinStereoLocal(SC->localPosition()));
+ if (!newCov) {
+ ATH_MSG_WARNING("Failed to create scaled error for SCT");
+ return 0;
+ }
+ cov = *newCov;
+ delete newCov;
+ }
+ double Sn = EL->sinStereoLocal(SC->localPosition());
+ double Sn2 = Sn * Sn;
+ double Cs2 = (1. - Sn) * (1. + Sn);
+ double SC = Sn * sqrt(Cs2);
+ double W = EL->phiPitch(loct) / EL->phiPitch();
+ double dV0 = (Cs2 * cov(0, 0) + Sn2 * cov(1, 1) +
+ 2. * SC * cov(1, 0)) * (W * W - 1.);
+ cov(0, 0) += (Cs2 * dV0);
+ cov(1, 0) += (SC * dV0);
+ cov(0, 1) = cov(1, 0);
+ cov(1, 1) += (Sn2 * dV0);
+ }
+
+ if (m_doDistortions) {
+ if (EL->isBarrel() == 1) {// Only apply disortions to barrel modules
+ locpar[Trk::locX] -= m_distortionsTool->correctReconstruction(trackPar, *EL, locpar, loct);
+ }
+ }
+
+
+ // Apply correction for cluster position bias
+ //
+ if (m_option_correctionStrategy == 0) {
+ double correction = getCorrection(dphi, int(colRow.x())) * EL->hitDepthDirection();
+ locpar[Trk::locX] += correction;
+ }
+ bool isbroad = (m_option_errorStrategy == 0) ? true : false;
+
+ msg(MSG::VERBOSE) << "returning corrected SCT cluster-on-track with position locX= " << locpar[Trk::locX] <<endmsg;
+
+ InDet::SCT_ClusterOnTrack* clus=nullptr;
+
+ Amg::Vector3D localstripdir(-sinAlpha, cosAlpha, 0.);
+ Amg::Vector3D globalstripdir = trackPar.associatedSurface().transform().linear() * localstripdir;
+ double distance = (trackPar.position() - SC->globalPosition()).mag();
+ const Trk::TrapezoidBounds *tbounds =
+ dynamic_cast<const Trk::TrapezoidBounds *>(&trackPar.associatedSurface().bounds());
+ const Trk::RectangleBounds *rbounds =
+ dynamic_cast<const Trk::RectangleBounds *>(&trackPar.associatedSurface().bounds());
+
+ if (tbounds || rbounds) {
+
+ double boundsy = rbounds ? rbounds->halflengthY() : tbounds->halflengthY();
+ if (distance * cosAlpha > boundsy) {
+ distance = boundsy / cosAlpha - 1;
+ }
+
+
+ if (loct.y() < 0) {
+ distance = -distance;
+ }
+ Amg::Vector3D glob (SC->globalPosition() + distance * globalstripdir);
+ clus = new InDet::SCT_ClusterOnTrack(SC, locpar, cov, iH, glob, isbroad);
+ } else {
+ clus = new InDet::SCT_ClusterOnTrack(SC, locpar, cov, iH, SC->globalPosition(), isbroad);
+ }
+
+ return clus;
+}
+
+double
+FTK_SCTClusterOnTrackTool::getCorrection(double phi, int nstrip) const {
+ float corr1[30] = {
+ 0.3, 0.8, 1.1, 1.5, 1.9, 1.9, 2.1, 2.4, 2.3, 2.6,
+ 2.6, 2.7, 2.8, 2.7, 2.5, 2.6, 2.8, 2.6, 2.6, 2.7,
+ 2.2, 1.8, 1.8, 1.6, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0
+ };
+ float corr2[30] = {
+ 0.0, 0.0, 0.0, 1.0, 1.5, 1.7, 1.7, 2.3, 2.1, 2.5,
+ 2.5, 2.7, 2.7, 2.9, 3.0, 3.0, 3.0, 3.0, 3.4, 3.4,
+ 3.0, 3.2, 2.6, 2.6, 3.0, 2.7, 2.5, 2.4, 1.7, 1.3
+ };
+
+ // Phi bins have 1 degree width, and cover 0-30 degrees
+ int phiBin = int(fabs(phi) / deg);
+
+ float correction(0.);
+
+ if (phiBin < 30) {
+ if (nstrip == 1) {
+ correction = corr1[phiBin];
+ }
+ if (nstrip == 2) {
+ correction = corr2[phiBin];
+ }
+ }
+
+ if (phi > 0.) {
+ correction *= -1.;
+ }
+
+ return correction * micrometer;
+}
+
+double
+FTK_SCTClusterOnTrackTool::getError(double phi, int nstrip) const {
+ float sigma1[60] = {
+ 22.1, 21.8, 21.4, 21.0, 20.5, 20.0, 19.6, 19.1, 18.5, 18.0,
+ 17.4, 17.0, 16.4, 15.8, 15.4, 14.9, 14.4, 14.1, 13.3, 13.1,
+ 12.9, 12.4, 12.6, 12.2, 12.3, 12.6, 13.4, 14.2, 15.6, 19.3,
+ 22.8, 29.5, 33.2, 41.8, 44.3, 48.4, 49.9, 54.0, 53.0, 56.3,
+ 57.5, 56.3, 64.5, 65.7, 66.1, 69.4, 74.8, 78.3, 78.8, 79.8,
+ 73.5, 73.8, 75.8, 84.3, 87.0, 99.9, 86.3, 0.0, 0.0, 0.0
+ };
+ float sigma2[60] = {
+ 22.2, 20.3, 18.8, 16.0, 14.6, 13.8, 12.9, 12.9, 12.7, 12.3,
+ 12.7, 12.6, 13.0, 13.3, 14.0, 14.6, 15.3, 15.9, 16.6, 17.6,
+ 18.4, 19.3, 19.9, 20.5, 21.0, 21.2, 21.5, 21.4, 21.3, 21.3,
+ 20.9, 20.8, 20.6, 20.7, 20.3, 20.7, 21.7, 24.4, 26.5, 29.5,
+ 34.6, 41.6, 48.5, 52.3, 54.5, 58.4, 61.8, 66.7, 69.9, 72.1,
+ 78.9, 79.2, 81.8, 80.9, 87.5, 99.2, 0.0, 0.0, 0.0, 0.0
+ };
+ float sigma3[60] = {
+ 70.1, 73.6, 71.7, 66.9, 68.3, 66.8, 66.2, 64.8, 66.6, 63.3,
+ 63.3, 60.4, 59.0, 57.1, 56.4, 54.4, 54.2, 54.4, 50.3, 48.9,
+ 48.1, 41.9, 38.0, 31.8, 28.3, 23.1, 23.0, 20.3, 18.5, 17.6,
+ 17.7, 16.8, 18.3, 19.3, 19.0, 20.0, 20.9, 21.6, 22.0, 22.2,
+ 22.7, 22.4, 24.3, 24.8, 24.6, 27.0, 29.8, 37.0, 47.7, 49.3,
+ 58.2, 60.2, 66.8, 70.8, 77.3, 80.6, 0.0, 0.0, 0.0, 0.0
+ };
+ float sigma4[60] = {
+ 103.2, 100.4, 100.7, 101.2, 107.4, 100.6, 100.9, 100.4, 96.3, 98.2,
+ 96.7, 94.5, 96.9, 91.7, 90.5, 89.5, 86.3, 90.6, 82.4, 89.3,
+ 87.3, 77.6, 75.7, 77.2, 77.3, 84.1, 80.1, 66.9, 73.7, 72.3,
+ 58.1, 65.6, 64.2, 54.7, 47.2, 44.4, 34.6, 36.4, 29.1, 25.8,
+ 18.8, 21.6, 18.6, 20.3, 22.7, 23.3, 24.1, 22.4, 24.7, 24.7,
+ 27.3, 30.4, 37.0, 46.4, 59.4, 62.6, 65.3, 0.0, 0.0, 0.0
+ };
+ float sigma5[60] = {
+ 150.9, 139.7, 133.9, 139.8, 141.4, 134.9, 138.4, 129.3, 137.9, 128.7,
+ 132.4, 130.1, 124.2, 115.8, 131.4, 115.2, 128.7, 112.8, 130.7, 129.0,
+ 115.8, 101.3, 115.9, 116.1, 121.7, 109.9, 110.0, 97.2, 96.4, 107.3,
+ 98.2, 80.0, 73.2, 87.0, 97.0, 88.5, 72.2, 73.9, 80.8, 75.7,
+ 69.5, 67.1, 54.1, 58.9, 47.3, 50.6, 29.5, 26.6, 25.8, 20.9,
+ 20.6, 21.9, 22.1, 21.1, 27.9, 41.6, 0.0, 0.0, 0.0, 0.0
+ };
+
+ // Phi bins have 1 degree width, and cover 0-60 degrees
+ int phiBin = int(fabs(phi) / deg);
+
+ float sigma(0.);
+
+ if (phiBin < 60) {
+ if (nstrip == 1) {
+ sigma = sigma1[phiBin];
+ }
+ if (nstrip == 2) {
+ sigma = sigma2[phiBin];
+ }
+ if (nstrip == 3) {
+ sigma = sigma3[phiBin];
+ }
+ if (nstrip == 4) {
+ sigma = sigma4[phiBin];
+ }
+ if (nstrip == 5) {
+ sigma = sigma5[phiBin];
+ }
+ }
+ if (sigma < 0.1) {
+ sigma = std::max(100., float(nstrip) * 80. / sqrt(12));
+ }
+
+ return sigma * micrometer;
+}
diff --git a/Trigger/TrigFTK/FTK_RecTools/src/IBL_calibration.h b/Trigger/TrigFTK/FTK_RecTools/src/IBL_calibration.h
new file mode 100644
index 0000000000000000000000000000000000000000..36a61072052d8329dc463ab646143d5061f97f54
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/src/IBL_calibration.h
@@ -0,0 +1,76 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+calphi[0] = 0.0463822;
+calerrphi[0][0] = 0.00458704;
+calerrphi[0][1] = 0.00240123;
+calerrphi[0][2] = 0.0381553;
+calphi[1] = 0.0395409;
+calerrphi[1][0] = 0.00571953;
+calerrphi[1][1] = 0.00217718;
+calerrphi[1][2] = 0.0443544;
+calphi[2] = 0.0291725;
+calerrphi[2][0] = 0.00753389;
+calerrphi[2][1] = 0.00192153;
+calerrphi[2][2] = 0.0502262;
+calphi[3] = 0.0221493;
+calerrphi[3][0] = 0.00908161;
+calerrphi[3][1] = 0.00234422;
+calerrphi[3][2] = 0.052859;
+calphi[4] = 0.0183585;
+calerrphi[4][0] = 0.009916;
+calerrphi[4][1] = 0.00307481;
+calerrphi[4][2] = 0.0560763;
+calphi[5] = 0.0228795;
+calerrphi[5][0] = 0.00886146;
+calerrphi[5][1] = 0.00260481;
+calerrphi[5][2] = 0.0516711;
+calphi[6] = 0.0296685;
+calerrphi[6][0] = 0.00746781;
+calerrphi[6][1] = 0.00197842;
+calerrphi[6][2] = 0.0497699;
+calphi[7] = 0.0395992;
+calerrphi[7][0] = 0.00557493;
+calerrphi[7][1] = 0.00213298;
+calerrphi[7][2] = 0.0454933;
+calphi[8] = 0.0497221;
+calerrphi[8][0] = 0.00461163;
+calerrphi[8][1] = 0.00251781;
+calerrphi[8][2] = 0.0309285;
+caleta[0][0] = 0;
+caleta[0][1] = 0;
+caleta[0][2] = 0.0669391;
+calerreta[0][0] = 0;
+calerreta[0][1] = 0.0614102;
+calerreta[0][2] = 0.00891571;
+caleta[1][0] = 0.0669391;
+caleta[1][1] = 0.199526;
+caleta[1][2] = 0.0642114;
+calerreta[1][0] = 0.0445204;
+calerreta[1][1] = 0.05385;
+calerreta[1][2] = 0.0104638;
+caleta[2][0] = 0.322464;
+caleta[2][1] = 0.207504;
+caleta[2][2] = 0.0699703;
+calerreta[2][0] = 0.045776;
+calerreta[2][1] = 0.0222315;
+calerreta[2][2] = 0.0106168;
+caleta[3][0] = 0.299711;
+caleta[3][1] = 0.215319;
+caleta[3][2] = 0.0653161;
+calerreta[3][0] = 0.0277008;
+calerreta[3][1] = 0.0218932;
+calerreta[3][2] = 0.0108854;
+caleta[4][0] = 0.295881;
+caleta[4][1] = 0.215595;
+caleta[4][2] = 0.0284735;
+calerreta[4][0] = 0.0281094;
+calerreta[4][1] = 0.0215104;
+calerreta[4][2] = 0.0159464;
+caleta[5][0] = 0.321041;
+caleta[5][1] = 0.216542;
+caleta[5][2] = 0.172694;
+calerreta[5][0] = 0.0312848;
+calerreta[5][1] = 0.0236938;
+calerreta[5][2] = 0.0202331;
diff --git a/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx b/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx
index 5c4c30317cd5afafce2cf4cafc893f2f8a88d0e5..f376046dab2dd39d49f034fe50f97b6bb725c999 100644
--- a/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx
+++ b/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx
@@ -1,9 +1,15 @@
#include "GaudiKernel/DeclareFactoryEntries.h"
#include "FTK_RecTools/FTK_VertexFinderTool.h"
+#include "FTK_RecTools/FTK_PixelClusterOnTrackTool.h"
+#include "FTK_RecTools/FTK_SCTClusterOnTrackTool.h"
DECLARE_TOOL_FACTORY(FTK_VertexFinderTool)
+DECLARE_TOOL_FACTORY(FTK_PixelClusterOnTrackTool)
+DECLARE_TOOL_FACTORY(FTK_SCTClusterOnTrackTool)
DECLARE_FACTORY_ENTRIES(FTK_RecTools) {
DECLARE_TOOL(FTK_VertexFinderTool)
+ DECLARE_TOOL(FTK_PixelClusterOnTrackTool)
+ DECLARE_TOOL(FTK_SCTClusterOnTrackTool)
}