From 68b1aade510758cd0e675cd9a5101ff0f35c228f Mon Sep 17 00:00:00 2001 From: Thomas Strebler <thomas.strebler@cern.ch> Date: Sun, 9 May 2021 17:36:35 +0000 Subject: [PATCH] Revert "Updated PixelClusterOnTrackTool for ITk analogue clustering" This reverts commit 29104c2b73d098c4e3f9abf563d8ecc5219084e3. --- .../python/PixelConditionsConfig.py | 9 + .../src/ITkPixelOfflineCalibCondAlg.cxx | 210 ++++++ .../src/ITkPixelOfflineCalibCondAlg.h | 54 ++ .../PixelConditionsAlgorithms_entries.cxx | 2 + .../PixelConditionsData/CMakeLists.txt | 2 +- .../ITkPixelClusterErrorData.h | 56 ++ .../ITkPixelOfflineCalibData.h | 99 +++ .../src/ITkPixelClusterErrorData.cxx | 147 ++++ .../src/ITkPixelOfflineCalibData.cxx | 81 ++ .../ITkPixelClusterOnTrackTool.h | 171 +++++ .../src/ITkPixelClusterOnTrackTool.cxx | 692 ++++++++++++++++++ .../src/PixelClusterOnTrackTool.cxx | 14 +- .../SiClusterOnTrackTool_entries.cxx | 3 +- 13 files changed, 1528 insertions(+), 12 deletions(-) create mode 100644 InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.cxx create mode 100644 InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.h create mode 100644 InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelClusterErrorData.h create mode 100644 InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelOfflineCalibData.h create mode 100644 InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelClusterErrorData.cxx create mode 100644 InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelOfflineCalibData.cxx create mode 100755 InnerDetector/InDetRecTools/SiClusterOnTrackTool/SiClusterOnTrackTool/ITkPixelClusterOnTrackTool.h create mode 100755 InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/ITkPixelClusterOnTrackTool.cxx diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py index 8c6c42c9e4b..6b412f491d8 100644 --- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py +++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py @@ -503,3 +503,12 @@ def PixelReadoutSpeedAlgCfg(flags, name="PixelReadoutSpeedAlg", **kwargs): acc.addCondAlgo(CompFactory.PixelReadoutSpeedAlg(name, **kwargs)) return acc +def ITkPixelOfflineCalibCondAlgCfg(flags, name="ITkPixelOfflineCalibCondAlg", **kwargs): + """Return a ComponentAccumulator with configured ITkPixelOfflineCalibCondAlg""" + acc = ComponentAccumulator() + acc.merge(addFolders(flags, "/PIXEL/ITkClusterError", "PIXEL_OFL", className="CondAttrListCollection")) + kwargs.setdefault("ReadKey", "/PIXEL/ITkClusterError") + kwargs.setdefault("WriteKey", "ITkPixelOfflineCalibData") + kwargs.setdefault("InputSource", 2) + acc.addCondAlgo(CompFactory.ITkPixelOfflineCalibCondAlg(name, **kwargs)) + return acc diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.cxx new file mode 100644 index 00000000000..b9acb6e0cf5 --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.cxx @@ -0,0 +1,210 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ITkPixelOfflineCalibCondAlg.h" +#include "Identifier/Identifier.h" +#include "Identifier/IdentifierHash.h" +#include "GaudiKernel/EventIDRange.h" +#include "PathResolver/PathResolver.h" +#include <memory> +#include <sstream> + +ITkPixelOfflineCalibCondAlg::ITkPixelOfflineCalibCondAlg(const std::string& name, ISvcLocator* pSvcLocator): + ::AthReentrantAlgorithm(name, pSvcLocator) +{ +} + + +StatusCode ITkPixelOfflineCalibCondAlg::initialize() { + ATH_MSG_DEBUG("ITkPixelOfflineCalibCondAlg::initialize()"); + + ATH_CHECK(m_condSvc.retrieve()); + + ATH_CHECK(service("DetectorStore", m_detStore)); + ATH_CHECK(m_detStore->retrieve(m_pixelid, "PixelID")) ; + + if (m_inputSource==2 && m_readKey.key().empty()) { + ATH_MSG_ERROR("The database is set to be input source (2) but the ReadKey is empty."); + } + ATH_CHECK(m_readKey.initialize()); + + ATH_CHECK(m_writeKey.initialize()); + + if (m_condSvc->regHandle(this,m_writeKey).isFailure()) { + ATH_MSG_FATAL("unable to register WriteCondHandle " << m_writeKey.fullKey() << " with CondSvc"); + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + + + +StatusCode ITkPixelOfflineCalibCondAlg::execute(const EventContext& ctx) const { + ATH_MSG_DEBUG("ITkPixelOfflineCalibCondAlg::execute()"); + + SG::WriteCondHandle<ITkPixelCalib::ITkPixelOfflineCalibData> writeHandle(m_writeKey, ctx); + if (writeHandle.isValid()) { + ATH_MSG_DEBUG("CondHandle " << writeHandle.fullKey() << " is already valid.. In theory this should not be called, but may happen if multiple concurrent events are being processed out of order."); + return StatusCode::SUCCESS; + } + + // Construct the output Cond Object and fill it in + std::unique_ptr<ITkPixelCalib::ITkPixelOfflineCalibData> writeCdo(std::make_unique<ITkPixelCalib::ITkPixelOfflineCalibData>()); + + if (m_inputSource==0) { + ATH_MSG_WARNING("So far do nothing!! return StatusCode::FAILURE"); + return StatusCode::FAILURE; + } + else if (m_inputSource==1) { + ATH_MSG_WARNING("Pixel ITk constants read from text file. Only supported for local developments and debugging!"); + + ITkPixelCalib::ITkPixelOfflineCalibData* calibData = new ITkPixelCalib::ITkPixelOfflineCalibData; + + ITkPixelCalib::ITkPixelClusterErrorData* pced = calibData->getITkPixelClusterErrorData(); + + // Find and open the text file + ATH_MSG_INFO("Load ITkPixelErrorData constants from text file"); + std::string fileName = PathResolver::find_file(m_textFileName, "DATAPATH"); + if (fileName.size()==0) { ATH_MSG_WARNING("Input file " << fileName << " not found! Default (hardwired) values to be used!"); } + else { pced->load(fileName); } + + ATH_MSG_DEBUG("Get error constants"); + std::vector<float> constants = calibData->getConstants(); + if (constants.size()) { ATH_MSG_VERBOSE("constants are defined"); } + else { ATH_MSG_ERROR("constants size is NULL!!!"); } + + + const EventIDBase start{EventIDBase::UNDEFNUM, EventIDBase::UNDEFEVT, 0, 0, EventIDBase::UNDEFNUM, EventIDBase::UNDEFNUM}; + const EventIDBase stop {EventIDBase::UNDEFNUM, EventIDBase::UNDEFEVT, EventIDBase::UNDEFNUM-1, EventIDBase::UNDEFNUM-1, EventIDBase::UNDEFNUM, EventIDBase::UNDEFNUM}; + const EventIDRange rangeW{start, stop}; + + ATH_MSG_DEBUG("Range of input is " << rangeW); + + if (constants.size()) { + ATH_MSG_DEBUG("Found constants with new-style Identifier key"); + writeCdo->setConstants(constants); + } + + if (writeHandle.record(rangeW, std::move(writeCdo)).isFailure()) { + ATH_MSG_FATAL("Could not record PixelCalib::ITkPixelOfflineCalibData " << writeHandle.key() << " with EventRange " << rangeW << " into Conditions Store"); + return StatusCode::FAILURE; + } + ATH_MSG_DEBUG("recorded new CDO " << writeHandle.key() << " with range " << rangeW << " into Conditions Store"); + + if (m_dump!=0) { + ATH_MSG_DEBUG("Dump the constants to file"); + calibData->dump(); + } + delete calibData; + + } + + else if (m_inputSource==2) { + + SG::ReadCondHandle<CondAttrListCollection> readHandle{m_readKey, ctx}; + const CondAttrListCollection* readCdo{*readHandle}; + if (readCdo==nullptr) { + ATH_MSG_FATAL("Null pointer to the read conditions object"); + return StatusCode::FAILURE; + } + + // Get the validitiy range + EventIDRange rangeW; + if (not readHandle.range(rangeW)) { + ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key()); + return StatusCode::FAILURE; + } + + ATH_MSG_DEBUG("Range of input is " << rangeW); + + std::vector<float> constants; + + for(CondAttrListCollection::const_iterator attrList = readCdo->begin(); attrList != readCdo->end(); ++attrList){ + + std::ostringstream attrStr; + (*attrList).second.toOutputStream(attrStr); + ATH_MSG_DEBUG( "ChanNum " << (*attrList).first << " Attribute list " << attrStr.str() ); + + std::string stringData = (*attrList).second["data_array"].data<std::string>(); + std::string delimiter = "],"; + size_t pos = 0; + std::vector<std::string> component; + std::string buffer; + while ((pos = stringData.find(delimiter)) != std::string::npos) { + buffer = stringData.substr(0,pos); + component.push_back(buffer); + stringData.erase(0, pos + delimiter.length()); + } + component.push_back(stringData); + ATH_MSG_INFO("Last component="<<stringData); + + for (size_t i=0; i<component.size(); i++) { + std::string checkModule = component[i]; + std::vector<std::string> moduleString; + delimiter = ":["; + pos = 0; + while ((pos = checkModule.find(delimiter)) != std::string::npos) { + buffer = checkModule.substr(0,pos); + moduleString.push_back(buffer); + checkModule.erase(0, pos + delimiter.length()); + } + moduleString.push_back(checkModule); + + if (moduleString.size()!=2) { + ATH_MSG_FATAL("String size (moduleString) is not 2. " << moduleString.size() << " in " << component[i] << " channel " << (*attrList).first << " read from " << readHandle.fullKey()); + return StatusCode::FAILURE; + } + + std::stringstream checkModuleHash(moduleString[0]); + std::vector<std::string> moduleStringHash; + while (std::getline(checkModuleHash,buffer,'"')) { moduleStringHash.push_back(buffer); } + + int waferHash = std::atoi(moduleStringHash[1].c_str()); + IdentifierHash waferID_hash(waferHash); + Identifier pixelID = m_pixelid->wafer_id(waferID_hash); + constants.emplace_back( pixelID.get_compact() ); + + std::stringstream moduleConstants(moduleString[1]); + std::vector<float> moduleConstantsVec; + while (std::getline(moduleConstants,buffer,',')) { moduleConstantsVec.emplace_back(std::atof(buffer.c_str())); } + + // Format v1 with no incident angle dependance + if(moduleConstantsVec.size()==4){ + constants.emplace_back(0); // period_phi + constants.emplace_back(0); // period_sinheta + constants.emplace_back(0); // delta_x_slope + constants.emplace_back(moduleConstantsVec[0]); // delta_x_offset + constants.emplace_back(moduleConstantsVec[1]); // delta_error_x + constants.emplace_back(0); // delta_y_slope + constants.emplace_back(moduleConstantsVec[2]); // delta_y_offset + constants.emplace_back(moduleConstantsVec[3]); // delta_error_y + } + + else if(moduleConstantsVec.size()==7){ + constants.emplace_back(moduleConstantsVec[0]); // period_phi + for( auto& x : moduleConstantsVec ) constants.emplace_back(x); + } + + // Format v3 with incident angle dependance + different eta-phi periods + else if(moduleConstantsVec.size()==8){ + for( auto& x : moduleConstantsVec ) constants.emplace_back(x); + } + + } + + } + + writeCdo->setConstants(constants); + + if (writeHandle.record(rangeW, std::move(writeCdo)).isFailure()) { + ATH_MSG_FATAL("Could not record PixelCalib::ITkPixelOfflineCalibData " << writeHandle.key() << " with EventRange " << rangeW << " into Conditions Store"); + return StatusCode::FAILURE; + } + ATH_MSG_DEBUG("recorded new CDO " << writeHandle.key() << " with range " << rangeW << " into Conditions Store"); + + } + + return StatusCode::SUCCESS; +} diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.h b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.h new file mode 100644 index 00000000000..9eba9ae26ed --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/ITkPixelOfflineCalibCondAlg.h @@ -0,0 +1,54 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ITKPIXELOFFLINECALIBCONDALG +#define ITKPIXELOFFLINECALIBCONDALG + +#include "AthenaBaseComps/AthReentrantAlgorithm.h" + +#include "StoreGate/ReadCondHandleKey.h" +#include "AthenaPoolUtilities/CondAttrListCollection.h" + +#include "StoreGate/WriteCondHandleKey.h" +#include "PixelConditionsData/ITkPixelOfflineCalibData.h" + +#include "GaudiKernel/ICondSvc.h" +#include "Gaudi/Property.h" + +#include "StoreGate/StoreGateSvc.h" +#include "InDetIdentifier/PixelID.h" + +class ITkPixelOfflineCalibCondAlg : public AthReentrantAlgorithm { + public: + ITkPixelOfflineCalibCondAlg(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~ITkPixelOfflineCalibCondAlg() = default; + + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext& ctx) const override; + + private: + Gaudi::Property<int> m_inputSource + {this, "InputSource",2,"Source of data: 0 (none), 1 (text file), 2 (database)"}; + + Gaudi::Property<std::string> m_textFileName + {this, "ITkPixelClusterErrorDataFile", "ITkPixelClusterErrorData.txt","Read constants from this file"}; + + Gaudi::Property<int> m_dump + {this, "DumpConstants", 0, "Dump constants to text file"}; + + SG::ReadCondHandleKey<CondAttrListCollection> m_readKey + {this, "ReadKey", "/PIXEL/ITkClusterError", "Input key of pixreco conditions folder"}; + + SG::WriteCondHandleKey<ITkPixelCalib::ITkPixelOfflineCalibData> m_writeKey + {this, "WriteKey", "ITkPixelOfflineCalibData", "Output key of pixel module data"}; + + ServiceHandle<ICondSvc> m_condSvc{this, "CondSvc", "CondSvc"}; + + StoreGateSvc* m_detStore{nullptr}; + const PixelID* m_pixelid{nullptr}; + +}; + + +#endif diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/components/PixelConditionsAlgorithms_entries.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/components/PixelConditionsAlgorithms_entries.cxx index 512c204d515..e3b1d665e37 100644 --- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/components/PixelConditionsAlgorithms_entries.cxx +++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/components/PixelConditionsAlgorithms_entries.cxx @@ -16,6 +16,7 @@ #include "../PixelAlignCondAlg.h" #include "../PixelDetectorElementCondAlg.h" #include "../PixeldEdxAlg.h" +#include "../ITkPixelOfflineCalibCondAlg.h" DECLARE_COMPONENT( PixelDCSCondHVAlg ) DECLARE_COMPONENT( PixelDCSCondTempAlg ) @@ -35,3 +36,4 @@ DECLARE_COMPONENT( PixelHitDiscCnfgAlg ) DECLARE_COMPONENT( PixelAlignCondAlg ) DECLARE_COMPONENT( PixelDetectorElementCondAlg ) DECLARE_COMPONENT( PixeldEdxAlg ) +DECLARE_COMPONENT( ITkPixelOfflineCalibCondAlg ) diff --git a/InnerDetector/InDetConditions/PixelConditionsData/CMakeLists.txt b/InnerDetector/InDetConditions/PixelConditionsData/CMakeLists.txt index 31096f9226d..1181ac47bbf 100644 --- a/InnerDetector/InDetConditions/PixelConditionsData/CMakeLists.txt +++ b/InnerDetector/InDetConditions/PixelConditionsData/CMakeLists.txt @@ -13,5 +13,5 @@ atlas_add_library( PixelConditionsData PUBLIC_HEADERS PixelConditionsData PRIVATE_INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} PRIVATE_DEFINITIONS ${CLHEP_DEFINITIONS} - LINK_LIBRARIES AthenaKernel AthenaPoolUtilities GeoPrimitives Identifier InDetByteStreamErrors + LINK_LIBRARIES AthenaKernel AthenaPoolUtilities GeoPrimitives Identifier InDetByteStreamErrors InDetIdentifier PRIVATE_LINK_LIBRARIES ${CLHEP_LIBRARIES} ${ROOT_LIBRARIES} ) diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelClusterErrorData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelClusterErrorData.h new file mode 100644 index 00000000000..a22eca82863 --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelClusterErrorData.h @@ -0,0 +1,56 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ITKPIXELCLUSTERERRORDATA_H +#define ITKPIXELCLUSTERERRORDATA_H + +#include "AthenaKernel/CLASS_DEF.h" + +#include "Identifier/Identifier.h" +#include "InDetIdentifier/PixelID.h" +#include "StoreGate/StoreGateSvc.h" + +#include <string> +#include <vector> + +namespace ITkPixelCalib { + +class ITkPixelClusterErrorData { + + public: + ITkPixelClusterErrorData(){ Initialize(); } + ~ITkPixelClusterErrorData(){}; + + /** Methods to access the calibration data */ + + std::pair<double,double> getDelta(const Identifier* pixelId, + int sizePhi, double angle, + int sizeZ, double eta) const; + std::pair<double,double> getDeltaError(const Identifier* pixelId) const; + std::map< const Identifier, std::vector<double> > getConstMap() const {return m_constmap;} + + void setDeltaError(const Identifier* pixelId, + double period_phi, double period_sinheta, + double delta_x_slope, double delta_x_offset, double error_x, + double delta_y_slope, double delta_y_offset, double error_y ); + + void print(std::string file) const; + void load(std::string file); + + + private: + void Initialize(); + + // map to store all ITk Analogue Clustering constants and errors + std::map< const Identifier, std::vector<double> > m_constmap; + + StoreGateSvc* m_detStore{nullptr}; + const PixelID* m_pixelID{nullptr}; + +}; + +} + +#endif + diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelOfflineCalibData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelOfflineCalibData.h new file mode 100644 index 00000000000..292a094973c --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/ITkPixelOfflineCalibData.h @@ -0,0 +1,99 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +/////////////////////////////////////////////////////////////////// +// ITkPixelOfflineCalibData.h, (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// + +#ifndef ITKPIXELOFFLINECALIBDATA_H +#define ITKPIXELOFFLINECALIBDATA_H + +#include "AthenaKernel/CLASS_DEF.h" + +#include "PixelConditionsData/ITkPixelClusterErrorData.h" + +#include "AthenaKernel/CondCont.h" + +/** @class ITkPixelOfflineCalibData + + The ITkPixelOfflineCalibData is a class that designed to hold the + data used by ITk pixel offline algorithms. + +*/ + +namespace ITkPixelCalib { + +class ITkPixelOfflineCalibData{ + + public: + /** Constructor:*/ + ITkPixelOfflineCalibData(); + ITkPixelOfflineCalibData(const ITkPixelOfflineCalibData& rhs); + ITkPixelOfflineCalibData& operator=(const ITkPixelOfflineCalibData& rhs); + + /** default destructor */ + ~ITkPixelOfflineCalibData (); + + bool update(const ITkPixelClusterErrorData& idat); + + // get the pointer to pixel cluster error data + ITkPixelClusterErrorData* getITkPixelClusterErrorData(); + const ITkPixelClusterErrorData* getITkPixelClusterErrorData() const; + + std::vector<float> getConstants() const; + void setConstants(const std::vector<float> &constants); + + void dump(); + + private: + ITkPixelClusterErrorData* m_clustererrordata; + + +}; + + +inline ITkPixelOfflineCalibData::ITkPixelOfflineCalibData() { + m_clustererrordata = new ITkPixelClusterErrorData(); +} + +inline ITkPixelOfflineCalibData::ITkPixelOfflineCalibData(const ITkPixelOfflineCalibData& rhs){ + m_clustererrordata = rhs.m_clustererrordata; +} + +inline ITkPixelOfflineCalibData& ITkPixelOfflineCalibData::operator=(const ITkPixelOfflineCalibData& rhs){ + if(this != &rhs){ + m_clustererrordata = rhs.m_clustererrordata; + } + return (*this); +} + +inline bool ITkPixelOfflineCalibData::update(const ITkPixelClusterErrorData& idat){ + *m_clustererrordata = idat; + return true; +} + +inline ITkPixelClusterErrorData* ITkPixelOfflineCalibData::getITkPixelClusterErrorData() { + return m_clustererrordata; +} + +inline const ITkPixelClusterErrorData* ITkPixelOfflineCalibData::getITkPixelClusterErrorData() const { + return m_clustererrordata; +} + + + +inline ITkPixelOfflineCalibData::~ITkPixelOfflineCalibData(){ + m_clustererrordata = nullptr; //Needed to avoid segfault when destructor is called + delete m_clustererrordata; +} + +} + + +CLASS_DEF( ITkPixelCalib::ITkPixelOfflineCalibData , 222566141 , 1 ) +CLASS_DEF( CondCont<ITkPixelCalib::ITkPixelOfflineCalibData> , 115484743 , 1 ) + +#endif + + diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelClusterErrorData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelClusterErrorData.cxx new file mode 100644 index 00000000000..822c7095055 --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelClusterErrorData.cxx @@ -0,0 +1,147 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +#include "PixelConditionsData/ITkPixelClusterErrorData.h" + +#include "GaudiKernel/ISvcLocator.h" + +#include "Identifier/IdentifierHash.h" + +#include <fstream> +#include <string> +#include <stdexcept> + + +namespace ITkPixelCalib{ + +void ITkPixelClusterErrorData::Initialize(){ + + ISvcLocator* svcLoc = Gaudi::svcLocator(); + StatusCode sc = svcLoc->service("DetectorStore", m_detStore); + if(sc.isFailure()){ + throw std::runtime_error("Could not retrieve DetectorStore"); + } + sc = m_detStore->retrieve(m_pixelID, "PixelID"); + if(sc.isFailure()){ + throw std::runtime_error("Could not retrieve PixelID"); + } + +} + + +std::pair<double,double> ITkPixelClusterErrorData::getDelta(const Identifier* pixelId, + int sizePhi, double angle, + int sizeZ, double eta) const{ + + std::vector<double> value = m_constmap.at(*pixelId); + double period_phi = value[0]; + double period_sinheta = value[1]; + double delta_x_slope = value[2]; + double delta_x_offset = value[3]; + double delta_y_slope = value[5]; + double delta_y_offset = value[6]; + + double delta_x = delta_x_slope * fabs(angle - period_phi*(sizePhi-2)) + delta_x_offset; + double delta_y = delta_y_slope * fabs(sinh(fabs(eta)) - period_sinheta*(sizeZ-2)) + delta_y_offset; + + return std::make_pair(delta_x,delta_y); + +} + + +std::pair<double,double> ITkPixelClusterErrorData::getDeltaError(const Identifier* pixelId) const{ + + std::vector<double> value = m_constmap.at(*pixelId); + + double delta_x_error = value[4]; + double delta_y_error = value[7]; + + return std::make_pair(delta_x_error,delta_y_error); + +} + + +// SET METHODS + +void ITkPixelClusterErrorData::setDeltaError(const Identifier* pixelId, + double period_phi, double period_sinheta, + double delta_x_slope, double delta_x_offset, double error_x, + double delta_y_slope, double delta_y_offset, double error_y){ + + std::vector<double> linevalues = {period_phi, period_sinheta, + delta_x_slope, delta_x_offset, error_x, + delta_y_slope, delta_y_offset, error_y}; + + m_constmap[*pixelId] = linevalues; + return; + +} + + +// save all constants to file +void ITkPixelClusterErrorData::print(std::string file) const { + + std::ofstream* outfile = new std::ofstream(file.c_str()); + + for(auto& x : m_constmap){ + + std::vector<double> value = x.second; + *outfile << m_pixelID->wafer_hash(x.first) << " " << value[0] << " " << value[1] << " " << value[2] << " " << value[3] << " " << value[4] << " " << value[5] << " " << value[6] << " " << value[7] << std::endl; + + } + + outfile->close(); + delete outfile; +} + + + +// Load ITk constants from file +void ITkPixelClusterErrorData::load(std::string file){ + + std::ifstream infile( file.c_str() ); + + if(infile.is_open()){ + + // + // Data in the file is stored in the following columns: + // waferID_hash : period_phi : period_sinheta : delta_x_slope : delta_x_offset : delta_error_x : delta_y_slope : delta_y_offset : delta_error_y + // + + int waferID_hash_int; + double period_phi; + double period_sinheta; + double delta_x_slope; + double delta_x_offset; + double delta_error_x; + double delta_y_slope; + double delta_y_offset; + double delta_error_y; + + while(!infile.eof()){ + + infile >> waferID_hash_int >> period_phi >> period_sinheta >> delta_x_slope >> delta_x_offset >> delta_error_x >> delta_y_slope >> delta_y_offset >> delta_error_y; + + IdentifierHash waferID_hash(waferID_hash_int); + Identifier pixelId = m_pixelID->wafer_id(waferID_hash); + setDeltaError(&pixelId, + period_phi, period_sinheta, + delta_x_slope, delta_x_offset, delta_error_x, + delta_y_slope, delta_y_offset, delta_error_y); + + } + + infile.close(); + + }else{ + throw std::runtime_error("ITkAnalogueClusteringConstantsFile \"" + file + "\" can not be read. Unable to proceed."); + } + +} + + +} + + + diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelOfflineCalibData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelOfflineCalibData.cxx new file mode 100644 index 00000000000..5f88ee6e62c --- /dev/null +++ b/InnerDetector/InDetConditions/PixelConditionsData/src/ITkPixelOfflineCalibData.cxx @@ -0,0 +1,81 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +/////////////////////////////////////////////////////////////////// +// ITkPixelOfflineCalibData.cxx, (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// + +#include "PixelConditionsData/ITkPixelOfflineCalibData.h" +#include "Identifier/Identifier.h" + +#include <sstream> + + +namespace ITkPixelCalib{ + + std::vector<float> ITkPixelOfflineCalibData::getConstants() const { + + std::map< const Identifier, std::vector<double> > constMap = m_clustererrordata->getConstMap(); + + int entry_size = 9; // pixel Id + period_phi + period_sinheta + delta_x_slope + delta_x_offset + delta_err_x + delta_y_slope + delta_y_offset + delta_err_y + int data_size = entry_size*constMap.size(); + + std::vector<float> constants; + constants.reserve(data_size); + + for(auto& x : constMap){ + + long long pixelId(x.first.get_compact()); + std::vector<double> value = x.second; + + constants.push_back(pixelId); + for(auto& y : value) constants.push_back(y); + + } + + return constants; + + } + + void ITkPixelOfflineCalibData::dump(){ + m_clustererrordata->print("ITkPixelClusterDump.txt"); + } + + + void ITkPixelOfflineCalibData::setConstants(const std::vector<float> &constants){ + + int entry_size = 9; + int map_size = constants.size()/entry_size; + + for(int i=0;i<map_size;i++){ + + long long pixelId_long = constants[i*entry_size]; + std::ostringstream ss; + ss << "0x" << std::hex << pixelId_long; + std::string pixelId_str(ss.str()); + Identifier pixelId; + pixelId.set(pixelId_str); + + double period_phi = constants[i*entry_size + 1]; + double period_sinheta = constants[i*entry_size + 2]; + double delta_x_slope = constants[i*entry_size + 3]; + double delta_x_offset = constants[i*entry_size + 4]; + double delta_err_x = constants[i*entry_size + 5]; + double delta_y_slope = constants[i*entry_size + 6]; + double delta_y_offset = constants[i*entry_size + 7]; + double delta_err_y = constants[i*entry_size + 8]; + + m_clustererrordata->setDeltaError(&pixelId, period_phi, period_sinheta, + delta_x_slope, delta_x_offset, delta_err_x, + delta_y_slope, delta_y_offset, delta_err_y); + + } + + return; + + } + + + +} // end of namespace diff --git a/InnerDetector/InDetRecTools/SiClusterOnTrackTool/SiClusterOnTrackTool/ITkPixelClusterOnTrackTool.h b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/SiClusterOnTrackTool/ITkPixelClusterOnTrackTool.h new file mode 100755 index 00000000000..27b41ed3bee --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/SiClusterOnTrackTool/ITkPixelClusterOnTrackTool.h @@ -0,0 +1,171 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ITkPixelClusterOnTrackTool_H +#define ITkPixelClusterOnTrackTool_H + +#include "GaudiKernel/ToolHandle.h" +#include "AthenaBaseComps/AthAlgTool.h" + +#include "TrkToolInterfaces/IRIO_OnTrackCreator.h" +#include "InDetRIO_OnTrack/PixelClusterOnTrack.h" +#include "InDetRIO_OnTrack/PixelRIO_OnTrackErrorScaling.h" + +#include "InDetPrepRawData/PixelGangedClusterAmbiguities.h" +#include "TrkParameters/TrackParameters.h" +#include "GeoPrimitives/GeoPrimitives.h" + +#include "Identifier/Identifier.h" +#include "InDetIdentifier/PixelID.h" + +#include "PixelConditionsData/PixelOfflineCalibData.h" +#include "PixelConditionsData/ITkPixelOfflineCalibData.h" +#include "PixelConditionsData/PixelDistortionData.h" +#include "InDetCondTools/ISiLorentzAngleTool.h" +#include "AthenaPoolUtilities/CondAttrListCollection.h" +#include "StoreGate/ReadCondHandleKey.h" + +#include "TrkEventUtils/ClusterSplitProbabilityContainer.h" + +#include <atomic> +#include <mutex> + +namespace InDet { + + /** @brief creates 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 + + */ + + class NnClusterizationFactory; + + enum ITkPixelClusterStrategy { + ITKPIXELCLUSTER_DEFAULT=0, + ITKPIXELCLUSTER_OUTLIER=1, + ITKPIXELCLUSTER_SHARED =2, + ITKPIXELCLUSTER_SPLIT =3 + }; + + class ITkPixelClusterOnTrackTool: + public AthAlgTool, virtual public Trk::IRIO_OnTrackCreator +{ + /////////////////////////////////////////////////////////////////// + // Public methods: + /////////////////////////////////////////////////////////////////// + +public: + + //! AlgTool constructor + ITkPixelClusterOnTrackTool(const std::string&,const std::string&, + const IInterface*); + virtual ~ITkPixelClusterOnTrackTool (); + //! AlgTool initialisation + virtual StatusCode initialize() override; + //! AlgTool termination + virtual StatusCode finalize () override; + + /** @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 override; + + /////////////////////////////////////////////////////////////////// + // Private methods: + /////////////////////////////////////////////////////////////////// + +protected: + + const InDet::PixelClusterOnTrack* correctDefault(const Trk::PrepRawData&, + const Trk::TrackParameters&) const; + + const InDet::PixelClusterOnTrack* correctNN(const Trk::PrepRawData&, const Trk::TrackParameters&) const; + + bool getErrorsDefaultAmbi( const InDet::PixelCluster*, const Trk::TrackParameters&, + Amg::Vector2D&, Amg::MatrixX&) const; + + bool getErrorsTIDE_Ambi( const InDet::PixelCluster*, const Trk::TrackParameters&, + Amg::Vector2D&, Amg::MatrixX&) const; + + const InDet::PixelClusterOnTrack* correct + (const Trk::PrepRawData&, const Trk::TrackParameters&, + const InDet::ITkPixelClusterStrategy) const; + + const Trk::ClusterSplitProbabilityContainer::ProbabilityInfo &getClusterSplittingProbability(const InDet::PixelCluster*pix) const { + if (!pix || m_clusterSplitProbContainer.key().empty()) return Trk::ClusterSplitProbabilityContainer::getNoSplitProbability(); + + SG::ReadHandle<Trk::ClusterSplitProbabilityContainer> splitProbContainer(m_clusterSplitProbContainer); + if (!splitProbContainer.isValid()) { + ATH_MSG_FATAL("Failed to get cluster splitting probability container " << m_clusterSplitProbContainer); + } + return splitProbContainer->splitProbability(pix); + } + +private: + + /////////////////////////////////////////////////////////////////// + // Private data: + /////////////////////////////////////////////////////////////////// + + ToolHandle<ISiLorentzAngleTool> m_lorentzAngleTool{this, "LorentzAngleTool", "SiLorentzAngleTool", "Tool to retreive Lorentz angle"}; + + SG::ReadCondHandleKey<ITkPixelCalib::ITkPixelOfflineCalibData> m_clusterITkErrorKey{this, "ITkPixelOfflineCalibData", "ITkPixelOfflineCalibData", "Output key of ITk pixel cluster"}; + + SG::ReadCondHandleKey<RIO_OnTrackErrorScaling> m_pixelErrorScalingKey + {this,"PixelErrorScalingKey", "/Indet/TrkErrorScalingPixel", "Key for pixel error scaling conditions data."}; + + //! toolhandle for central error scaling + //! flag storing if errors need scaling or should be kept nominal + int m_positionStrategy; + mutable std::atomic_int m_errorStrategy{2}; + IntegerProperty m_errorStrategyProperty{this, "ErrorStrategy", 2, "Which calibration of cluster position errors"}; + + + /** @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) **/ + bool m_applyNNcorrection{false}; + BooleanProperty m_applyNNcorrectionProperty{this, "applyNNcorrection", false}; + + /** NN clusterizationi factory for NN based positions and errors **/ + ToolHandle<NnClusterizationFactory> m_NnClusterizationFactory; + + 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; + SG::ReadHandleKey<InDet::PixelGangedClusterAmbiguities> m_splitClusterMapKey; + + SG::ReadHandleKey<Trk::ClusterSplitProbabilityContainer> m_clusterSplitProbContainer{this, "ClusterSplitProbabilityName", "",""}; + +}; + +} // end of namespace InDet + +#endif // ITkPixelClusterOnTrackTool_H diff --git a/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/ITkPixelClusterOnTrackTool.cxx b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/ITkPixelClusterOnTrackTool.cxx new file mode 100755 index 00000000000..368ebbabddb --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/ITkPixelClusterOnTrackTool.cxx @@ -0,0 +1,692 @@ +/* + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +*/ + +///////////////////////////////////////////////////////////////// +// Implementation file for class ITkPixelClusterOnTrackTool +/////////////////////////////////////////////////////////////////// +// (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// +// AlgTool used for PixelClusterOnTrack object production +/////////////////////////////////////////////////////////////////// +// started 1.0 21/04/2004 I.Gavrilenko - see ChangeLog +/////////////////////////////////////////////////////////////////// + +#include "SiClusterOnTrackTool/ITkPixelClusterOnTrackTool.h" +#include "PixelReadoutGeometry/PixelModuleDesign.h" +#include "InDetIdentifier/PixelID.h" +#include "TrkSurfaces/PlaneSurface.h" +#include "SiClusterizationTool/NnClusterizationFactory.h" +#include "EventPrimitives/EventPrimitives.h" +#include "InDetReadoutGeometry/SiDetectorElement.h" + +#include <cmath> +#include "TrkRIO_OnTrack/check_cast.h" + +//clustermap is most likely to be removed at later date +#define __clustermap + + +/////////////////////////////////////////////////////////////////// +// 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); + } +} + +InDet::ITkPixelClusterOnTrackTool::ITkPixelClusterOnTrackTool + (const std::string &t, const std::string &n, const IInterface *p) : + ::AthAlgTool(t, n, p), + m_pixelid(nullptr), + m_NnClusterizationFactory("InDet::NnClusterizationFactory/NnClusterizationFactory", this), + m_doNotRecalibrateNN(false), + m_noNNandBroadErrors(false), + m_usingTIDE_Ambi(false), + m_splitClusterMapKey("") + { + declareInterface<IRIO_OnTrackCreator>(this); + + declareProperty("PositionStrategy", m_positionStrategy = 1, "Which calibration of cluster positions"); + declareProperty("NnClusterizationFactory", m_NnClusterizationFactory); + declareProperty("SplitClusterAmbiguityMap", m_splitClusterMapKey);//Remove Later + declareProperty("doNotRecalibrateNN", m_doNotRecalibrateNN); + declareProperty("m_noNNandBroadErrors", m_noNNandBroadErrors); + declareProperty("RunningTIDE_Ambi", m_usingTIDE_Ambi); +} + +/////////////////////////////////////////////////////////////////// +// Destructor +/////////////////////////////////////////////////////////////////// + +InDet::ITkPixelClusterOnTrackTool::~ITkPixelClusterOnTrackTool() { +} + +/////////////////////////////////////////////////////////////////// +// Initialisation +/////////////////////////////////////////////////////////////////// + +StatusCode +InDet::ITkPixelClusterOnTrackTool::initialize() { + + ATH_MSG_DEBUG(name() << " initialize()"); + + m_errorStrategy = m_errorStrategyProperty; + ATH_MSG_DEBUG("Error strategy is" << m_errorStrategy); + + m_applyNNcorrection = m_applyNNcorrectionProperty; + + ATH_CHECK(m_clusterITkErrorKey.initialize()); + + ATH_CHECK(m_clusterSplitProbContainer.initialize( !m_clusterSplitProbContainer.key().empty())); + + // get the error scaling tool + ATH_CHECK(m_pixelErrorScalingKey.initialize(!m_pixelErrorScalingKey.key().empty())); + if (!m_pixelErrorScalingKey.key().empty()) ATH_MSG_DEBUG("Detected need for scaling Pixel errors."); + + ATH_CHECK (detStore()->retrieve(m_pixelid, "PixelID")); + + m_applyNNcorrection &= !m_splitClusterMapKey.key().empty(); + ATH_CHECK(m_splitClusterMapKey.initialize(m_applyNNcorrection)); + ATH_CHECK(m_NnClusterizationFactory.retrieve( DisableTool{!m_applyNNcorrection} )); + + ATH_CHECK(m_lorentzAngleTool.retrieve()); + return StatusCode::SUCCESS; +} + + + +/////////////////////////////////////////////////////////////////// +// Finalize +/////////////////////////////////////////////////////////////////// + +StatusCode +InDet::ITkPixelClusterOnTrackTool::finalize() { +return StatusCode::SUCCESS; +} + +/////////////////////////////////////////////////////////////////// +// Trk::SiClusterOnTrack production +/////////////////////////////////////////////////////////////////// + + +const InDet::PixelClusterOnTrack * +InDet::ITkPixelClusterOnTrackTool::correct + (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar) const { + + if (not m_applyNNcorrection){ + return correctDefault(rio, trackPar); + }else { + if (m_errorStrategy == 0 || m_errorStrategy == 1) { + // version from Giacinto + if (m_noNNandBroadErrors) { + return nullptr; + } + // 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 nullptr; + } + const Trk::ClusterSplitProbabilityContainer::ProbabilityInfo &splitProb = getClusterSplittingProbability(pix); + if (splitProb.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 * +InDet::ITkPixelClusterOnTrackTool::correctDefault + (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar) const { + using CLHEP::micrometer; + + const double TOPHAT_SIGMA = 1. / std::sqrt(12.); + + const InDet::PixelCluster *pix = nullptr; + + if (!(pix = dynamic_cast<const InDet::PixelCluster *>(&rio))) { + return nullptr; + } + + 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 nullptr; + } + 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 nullptr; + } else { + const InDetDD::PixelModuleDesign *design = + dynamic_cast<const InDetDD::PixelModuleDesign *>(&element->design()); + + // get candidate track angle in module local frame + Amg::Vector3D my_track = trackPar.momentum(); + const Amg::Vector3D& my_normal = element->normal(); + const Amg::Vector3D& my_phiax = element->phiAxis(); + const 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 = std::atan2(trkphicomp, trknormcomp); + double boweta = std::atan2(trketacomp, trknormcomp); + + float tanl = m_lorentzAngleTool->getTanLorentzAngle(iH); + 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 = std::atan(std::tan(bowphi) - readoutside * tanl); + + // 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)); + + // 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(); + InDetDD::SiLocalPosition meanpos(0, 0, 0); + int rowmin = 9999; + int rowmax = -9999; + int colmin = 9999; + int colmax = -9999; + for (const auto & rId:rdos) { + const int row = m_pixelid->phi_index(rId); + const int col = m_pixelid->eta_index(rId); + rowmin = std::min(rowmin, row); + rowmax = std::max(rowmax,row); + colmin = std::min(colmin, col); + colmax = std::max(colmax, col); + meanpos += design->positionFromColumnRow(col, row); + } + meanpos = meanpos / rdos.size(); + InDetDD::SiLocalPosition pos1 = + design->positionFromColumnRow(colmin, rowmin); + InDetDD::SiLocalPosition pos2 = + design->positionFromColumnRow(colmax, rowmin); + InDetDD::SiLocalPosition pos3 = + design->positionFromColumnRow(colmin, rowmax); + InDetDD::SiLocalPosition pos4 = + design->positionFromColumnRow(colmax, rowmax); + + InDetDD::SiLocalPosition centroid = 0.25 * (pos1 + pos2 + pos3 + pos4); + double shift = m_lorentzAngleTool->getLorentzShift(iH); + int nrows = rowmax - rowmin + 1; + int ncol = colmax - colmin + 1; + + // TOT interpolation for collision data + SG::ReadCondHandle<ITkPixelCalib::ITkPixelOfflineCalibData> offlineITkCalibDataHandle(m_clusterITkErrorKey); + + if (m_positionStrategy > 0 && omegaphi > -0.5 && omegaeta > -0.5) { + localphi = centroid.xPhi() + shift; + localeta = centroid.xEta(); + + std::pair<double,double> delta = offlineITkCalibDataHandle->getITkPixelClusterErrorData()->getDelta(&element_id,nrows,angle,ncol,etaloc); + double delta_phi = nrows != 1 ? delta.first : 0.; + double delta_eta = ncol != 1 ? delta.second : 0.; + localphi += delta_phi*(omegaphi-0.5); + localeta += delta_eta*(omegaeta-0.5); + } + // digital + else { + localphi = meanpos.xPhi() + shift; + localeta = meanpos.xEta(); + } + + const 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 (std::abs(angle) > 1) { + errphi = 250 * micrometer * std::tan(std::abs(angle)) * TOPHAT_SIGMA; + erreta = width.z() > 250 * micrometer * std::tan(std::abs(boweta)) ? + width.z() * TOPHAT_SIGMA : 250 * micrometer * std::tan(std::abs(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) { + std::pair<double,double> delta_err = offlineITkCalibDataHandle->getITkPixelClusterErrorData()->getDeltaError(&element_id); + errphi = nrows != 1 ? delta_err.first : (width.phiR()/nrows)*TOPHAT_SIGMA; + erreta = ncol != 1 ? delta_err.second : (width.z()/ncol)*TOPHAT_SIGMA; + } + + Amg::Vector2D locpos = Amg::Vector2D(localphi, localeta); + 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_pixelErrorScalingKey.key().empty()) { + SG::ReadCondHandle<RIO_OnTrackErrorScaling> error_scaling( m_pixelErrorScalingKey ); + cov = check_cast<PixelRIO_OnTrackErrorScaling>(*error_scaling)->getScaledCovariance( cov, *m_pixelid, element->identify() ); + } + bool isbroad = m_errorStrategy == 0; + return new InDet::PixelClusterOnTrack(pix, locpar, cov, iH, glob, pix->gangedPixel(), isbroad); +} + + +const InDet::PixelClusterOnTrack * +InDet::ITkPixelClusterOnTrackTool::correct + (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar, + const InDet::ITkPixelClusterStrategy strategy) const { + int initial_errorStrategy; + const InDet::PixelClusterOnTrack *newROT; + + switch (strategy) { + case InDet::ITKPIXELCLUSTER_OUTLIER: // if cluster is outlier, increase errors + case InDet::ITKPIXELCLUSTER_SHARED: + initial_errorStrategy = m_errorStrategy; + 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 * +InDet::ITkPixelClusterOnTrackTool::correctNN + (const Trk::PrepRawData &rio, + const Trk::TrackParameters &trackPar) const { + const InDet::PixelCluster *pixelPrepCluster = dynamic_cast<const InDet::PixelCluster *>(&rio); + + if (pixelPrepCluster == nullptr) { + ATH_MSG_WARNING("This is not a pixel cluster, return 0."); + return nullptr; + } + + const InDetDD::SiDetectorElement *element = pixelPrepCluster->detectorElement(); + if (!element) { + ATH_MSG_WARNING("Cannot access detector element. Aborting cluster correction..."); + return nullptr; + } + + IdentifierHash iH = element->identifyHash(); + + if (m_doNotRecalibrateNN) { + Amg::Vector3D glob(pixelPrepCluster->globalPosition()); + + Amg::Vector2D locpos = pixelPrepCluster->localPosition(); + Trk::LocalParameters locpar = Trk::LocalParameters(locpos); + Amg::MatrixX cov = pixelPrepCluster->localCovariance(); + + return new InDet::PixelClusterOnTrack(pixelPrepCluster, locpar, cov, iH, glob, + pixelPrepCluster->gangedPixel(), false); + } + + + + Amg::Vector2D finalposition; + Amg::MatrixX finalerrormatrix; + + if (m_usingTIDE_Ambi) { + if (!getErrorsTIDE_Ambi(pixelPrepCluster, trackPar, finalposition, finalerrormatrix)) { + return correctDefault(rio, trackPar); + } + }else { + if (!getErrorsDefaultAmbi(pixelPrepCluster, trackPar, finalposition, finalerrormatrix)) { + return correctDefault(rio, trackPar); + } + } + + ATH_MSG_DEBUG( " Old position x: " << pixelPrepCluster->localPosition()[0] + << " +/- " << std::sqrt(pixelPrepCluster->localCovariance()(0, 0)) + << " y: " << pixelPrepCluster->localPosition()[1] + << " +/- " << std::sqrt(pixelPrepCluster->localCovariance()(1, 1)) <<"\n" + << " Final position x: " << finalposition[0] + << " +/- " << std::sqrt(finalerrormatrix(0, 0)) + << " y: " << finalposition[1] << " +/- " + <<std::sqrt(finalerrormatrix(1, 1)) ); + + Amg::MatrixX cov = finalerrormatrix; + // create new copy of error matrix + if (!m_pixelErrorScalingKey.key().empty()) { + SG::ReadCondHandle<RIO_OnTrackErrorScaling> error_scaling( m_pixelErrorScalingKey ); + cov = check_cast<PixelRIO_OnTrackErrorScaling>(*error_scaling)->getScaledCovariance( cov, *m_pixelid, element->identify() ); + } + + InDetDD::SiLocalPosition centroid = InDetDD::SiLocalPosition(finalposition[1], + finalposition[0], + 0); + Trk::LocalParameters locpar = Trk::LocalParameters(finalposition); + + const Amg::Vector3D &glob = element->globalPosition(centroid); + + + return new InDet::PixelClusterOnTrack(pixelPrepCluster, locpar, cov, iH, + glob, + pixelPrepCluster->gangedPixel(), + false); +} + +bool +InDet::ITkPixelClusterOnTrackTool::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()); + + if (m_applyNNcorrection){ + SG::ReadHandle<InDet::PixelGangedClusterAmbiguities> splitClusterMap(m_splitClusterMapKey); + InDet::PixelGangedClusterAmbiguities::const_iterator mapBegin = splitClusterMap->begin(); + InDet::PixelGangedClusterAmbiguities::const_iterator mapEnd = splitClusterMap->end(); + for (InDet::PixelGangedClusterAmbiguities::const_iterator mapIter = mapBegin; mapIter != mapEnd; ++mapIter) { + const SiCluster *first = (*mapIter).first; + const SiCluster *second = (*mapIter).second; + if (first == pixelPrepCluster && second != pixelPrepCluster) { + ATH_MSG_DEBUG("Found additional split cluster in ambiguity map (+=1)."); + numberOfSubclusters += 1; + const SiCluster *otherOne = second; + const InDet::PixelCluster *pixelAddCluster = dynamic_cast<const InDet::PixelCluster *>(otherOne); + if (pixelAddCluster == nullptr) { + ATH_MSG_WARNING("Pixel ambiguity map has empty pixel cluster. Please DEBUG!"); + continue; + } + vectorOfPositions.push_back(pixelAddCluster->localPosition()); + + ATH_MSG_DEBUG( "Found one more pixel cluster. Position x: " + << pixelAddCluster->localPosition()[0] << "y: " << pixelAddCluster->localPosition()[1]); + }// find relevant element of map + }// iterate over map + } + + // now you have numberOfSubclusters and the vectorOfPositions (Amg::Vector2D) + + if (trackPar.surfaceType() != Trk::SurfaceType::Plane || + trackPar.type() != Trk::AtaSurface) { + 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, + trackPar.associatedSurface(), + trackPar, + allErrorMatrix, + numberOfSubclusters); + + if (allLocalPositions.empty()) { + ATH_MSG_DEBUG( " Cluster cannot be treated by NN. Giving back to default clusterization " ); + + 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; + } + + + // 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); + + ATH_MSG_DEBUG( + " Old pix (1) x: " << vectorOfPositions[0][0] << " y: " << vectorOfPositions[0][1] << "\n" + << " Old pix (2) x: " << vectorOfPositions[1][0] << " y: " << vectorOfPositions[1][1] << "\n" + << " Pix (1) x: " << allLocalPositions[0][0] << " +/- " << std::sqrt(allErrorMatrix[0](0, 0)) + << " y: " << allLocalPositions[0][1] << " +/- " << std::sqrt(allErrorMatrix[0](1, 1)) <<"\n" + << " Pix (2) x: " << allLocalPositions[1][0] << " +/- " << std::sqrt(allErrorMatrix[1](0, 0)) + << " y: " << allLocalPositions[1][1] << " +/- " << std::sqrt(allErrorMatrix[1](1, 1)) << "\n" + << " Old (1) new (1) dist: " << std::sqrt(distancesq1) << " Old (1) new (2) " << std::sqrt(distancesq2) ); + + + 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 +InDet::ITkPixelClusterOnTrackTool::getErrorsTIDE_Ambi(const InDet::PixelCluster *pixelPrepCluster, + const Trk::TrackParameters &trackPar, + Amg::Vector2D &finalposition, + Amg::MatrixX &finalerrormatrix) const { + const Trk::ClusterSplitProbabilityContainer::ProbabilityInfo &splitProb = getClusterSplittingProbability(pixelPrepCluster); + std::vector<Amg::Vector2D> vectorOfPositions; + int numberOfSubclusters = 1; + if(m_applyNNcorrection){ + SG::ReadHandle<InDet::PixelGangedClusterAmbiguities> splitClusterMap(m_splitClusterMapKey); + numberOfSubclusters = 1 + splitClusterMap->count(pixelPrepCluster); + + if (splitClusterMap->count(pixelPrepCluster) == 0 && splitProb.isSplit()) { + numberOfSubclusters = 2; + } + if (splitClusterMap->count(pixelPrepCluster) != 0 && !splitProb.isSplit()) { + numberOfSubclusters = 1; + } + } + + // now you have numberOfSubclusters and the vectorOfPositions (Amg::Vector2D) + if (trackPar.surfaceType() != Trk::SurfaceType::Plane || + trackPar.type() != Trk::AtaSurface) { + ATH_MSG_WARNING("Parameters are not at a plane surface ! Aborting cluster " + "correction... "); + return false; + } + + std::vector<Amg::Vector2D> allLocalPositions; + std::vector<Amg::MatrixX> allErrorMatrix; + allLocalPositions = m_NnClusterizationFactory->estimatePositions( + *pixelPrepCluster, + trackPar.associatedSurface(), + trackPar, + allErrorMatrix, + numberOfSubclusters); + + if (allLocalPositions.empty()) { + ATH_MSG_DEBUG( + " Cluster cannot be treated by NN. Giving back to default clusterization, too big: " << + splitProb.isTooBigToBeSplit()); + 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; + } + + // 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]; + 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(std::sqrt((*trackPar.covariance())(0, 0)), std::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]; + finalerrormatrix = allErrorMatrix[index]; + return true; +} diff --git a/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/PixelClusterOnTrackTool.cxx b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/PixelClusterOnTrackTool.cxx index 3a644d60f64..c5e1b4fa139 100755 --- a/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/PixelClusterOnTrackTool.cxx +++ b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/PixelClusterOnTrackTool.cxx @@ -114,18 +114,12 @@ InDet::PixelClusterOnTrackTool::initialize() { ATH_CHECK(m_clusterSplitProbContainer.initialize( !m_clusterSplitProbContainer.key().empty())); // get the error scaling tool - if (!m_pixelErrorScalingKey.key().empty()) { - ATH_CHECK(m_pixelErrorScalingKey.initialize()); - ATH_MSG_DEBUG("Detected need for scaling Pixel errors."); - } - + ATH_CHECK(m_pixelErrorScalingKey.initialize(!m_pixelErrorScalingKey.key().empty())); + if (!m_pixelErrorScalingKey.key().empty()) ATH_MSG_DEBUG("Detected need for scaling Pixel errors."); // get the module distortions tool - if (!m_disableDistortions) { - ATH_CHECK(m_distortionKey.initialize()); - } else { - ATH_MSG_INFO("No PixelDistortions will be simulated."); - } + ATH_CHECK(m_distortionKey.initialize(!m_disableDistortions)); + if (m_disableDistortions) ATH_MSG_INFO("No PixelDistortions will be simulated."); ATH_CHECK (detStore()->retrieve(m_pixelid, "PixelID")); diff --git a/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/components/SiClusterOnTrackTool_entries.cxx b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/components/SiClusterOnTrackTool_entries.cxx index 3f10389402e..a0fd5f09720 100644 --- a/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/components/SiClusterOnTrackTool_entries.cxx +++ b/InnerDetector/InDetRecTools/SiClusterOnTrackTool/src/components/SiClusterOnTrackTool_entries.cxx @@ -1,6 +1,7 @@ #include "SiClusterOnTrackTool/PixelClusterOnTrackTool.h" #include "SiClusterOnTrackTool/SCT_ClusterOnTrackTool.h" +#include "SiClusterOnTrackTool/ITkPixelClusterOnTrackTool.h" DECLARE_COMPONENT( InDet::PixelClusterOnTrackTool ) DECLARE_COMPONENT( InDet::SCT_ClusterOnTrackTool ) - +DECLARE_COMPONENT( InDet::ITkPixelClusterOnTrackTool ) -- GitLab