diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibInterfaces/MdtCalibInterfaces/IShiftMapTools.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibInterfaces/MdtCalibInterfaces/IShiftMapTools.h new file mode 100644 index 0000000000000000000000000000000000000000..b0bb26ddfc13f8292cfdf51f8cf84625bacc40a9 --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibInterfaces/MdtCalibInterfaces/IShiftMapTools.h @@ -0,0 +1,46 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +// IShiftMapTools.h +// Header file for class IShiftMapTools +/////////////////////////////////////////////////////////////////// +// (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// +// der.andi@cern.ch +/////////////////////////////////////////////////////////////////// + +#ifndef ISHIFTMAPTOOLS_H +#define ISHIFTMAPTOOLS_H + +// Include files +// from Gaudi +#include "GaudiKernel/IService.h" + +// forward declarations +class Identifier; + +namespace MuonCalib { + +/** + @class IShiftMapTools + Interface for the ShiftMapTools inheriting from MdtCalibrationShiftMapBase + + @author der.andi@cern.ch +*/ + +class IShiftMapTools : virtual public IService { + public: + IShiftMapTools() { ; } + virtual ~IShiftMapTools() { ; } + + /* Creates the InterfaceID and interfaceID() method */ + DeclareInterfaceID(IShiftMapTools, 1, 0); + + /* get shift value */ + virtual float getValue(const Identifier& id) = 0; +}; + +} // namespace MuonCalib + +#endif diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/CMakeLists.txt b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/CMakeLists.txt index 86a93a21015e44206a02663502586ffee8fca63a..82f43ac6457c142f7a9a16a563ea83ce27481400 100644 --- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/CMakeLists.txt +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/CMakeLists.txt @@ -12,6 +12,8 @@ atlas_depends_on_subdirs( PUBLIC DetectorDescription/GeoPrimitives DetectorDescription/Identifier GaudiKernel + MuonSpectrometer/MuonCablings/MuonMDT_Cabling + MuonSpectrometer/MuonCablings/MuonCablingData MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibData MuonSpectrometer/MuonCalib/MuonCalibIdentifier MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonPrepRawData @@ -21,7 +23,8 @@ atlas_depends_on_subdirs( PUBLIC MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibInterfaces MuonSpectrometer/MuonCalib/MuonCalibEvent MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry - MuonSpectrometer/MuonIdHelpers ) + MuonSpectrometer/MuonIdHelpers + Tools/PathResolver ) # External dependencies: find_package( Eigen ) @@ -31,14 +34,15 @@ atlas_add_library( MdtCalibSvcLib src/*.cxx PUBLIC_HEADERS MdtCalibSvc INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS} - LINK_LIBRARIES ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier GaudiKernel MdtCalibData MuonCalibIdentifier MuonPrepRawData StoreGateLib SGtests MuonIdHelpersLib + LINK_LIBRARIES ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier GaudiKernel MuonMDT_CablingLib MdtCalibData MuonCalibIdentifier MuonPrepRawData StoreGateLib SGtests MuonIdHelpersLib PathResolver PRIVATE_LINK_LIBRARIES MagFieldInterfaces MuonCalibEvent MuonReadoutGeometry ) atlas_add_component( MdtCalibSvc src/components/*.cxx INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS} - LINK_LIBRARIES ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier GaudiKernel MdtCalibData MuonCalibIdentifier MuonPrepRawData StoreGateLib SGtests MagFieldInterfaces MuonCalibEvent MuonReadoutGeometry MuonIdHelpersLib MdtCalibSvcLib ) + LINK_LIBRARIES MdtCalibSvcLib ) # Install files from the package: atlas_install_runtime( share/DC2_t0.dat share/DC2_rt.dat ) +atlas_install_joboptions( share/*.py ) diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationShiftMapBase.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationShiftMapBase.h new file mode 100644 index 0000000000000000000000000000000000000000..a8bde36667d50b4486297eb5aa2c202bef2fa775 --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationShiftMapBase.h @@ -0,0 +1,73 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MDTCALIBSVC_MDTCALIBRATIONSHIFTMAPBASE_H +#define MDTCALIBSVC_MDTCALIBRATIONSHIFTMAPBASE_H + +/* + * + * Author Andreas Hoenle (der.andi@cern.ch) + * + */ + +#include <map> + +// Framework includes +#include "AthenaBaseComps/AthService.h" + +// MDT interfaces +#include "MdtCalibInterfaces/IShiftMapTools.h" + +// MDT includes +#include "MuonMDT_Cabling/MuonMDT_CablingSvc.h" + +class Identifier; +class TTree; + +/** + @class MdtCalibrationShiftMapBase + Provides the base class for the per-tube shifting tools, like + MdtCalibT0ShiftTool & MdtCalibTMaxShiftTool. + @author Andreas Hoenle +*/ +class MdtCalibrationShiftMapBase + : public extends<AthService, MuonCalib::IShiftMapTools> { + public: + /* constructor */ + MdtCalibrationShiftMapBase(const std::string& name, ISvcLocator* pSvcLocator); + + /* destructor */ + ~MdtCalibrationShiftMapBase(); + + /* get shift value, override from IShiftMapTools */ + float getValue(const Identifier& id) override; + + /* + * initialization of map cannot happen before first event + * special function required + * we need the cabling service to be ready first + */ + virtual StatusCode initializeMap(); + + /* dump the map in binary file, given a path */ + StatusCode dumpMapAsFile(); + + /* load the map from a binary file, given a path */ + StatusCode loadMapFromFile(); + + protected: + ServiceHandle<MuonMDT_CablingSvc> m_cablingSvc; + std::map<Identifier, float> m_shiftValues; + bool m_mapIsInitialized; + + /* Make map configurable in job options */ + std::string m_mapFileName; + float m_centralValue; + float m_sigma; + bool m_forceMapRecreate; + + private: +}; + +#endif diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationT0ShiftSvc.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationT0ShiftSvc.h new file mode 100644 index 0000000000000000000000000000000000000000..f076a3e65e3ed045f16baa89de8b54e00e47691f --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationT0ShiftSvc.h @@ -0,0 +1,45 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MDTCALIBSVC_MDTCALIBRATIONT0SHIFTSVC_H +#define MDTCALIBSVC_MDTCALIBRATIONT0SHIFTSVC_H + +/* + * + * Author Andreas Hoenle (der.andi@cern.ch) + * + */ + +#include <map> + +#include "AthenaBaseComps/AthService.h" + +#include "MdtCalibSvc/MdtCalibrationShiftMapBase.h" + +class Identifier; +class TTree; + +/** + @class MdtCalibrationT0ShiftSvc + Provides a per-tube smearing of the T0 value. + @author Andreas Hoenle +*/ +class MdtCalibrationT0ShiftSvc : virtual public MdtCalibrationShiftMapBase { +public: + /** constructor */ + MdtCalibrationT0ShiftSvc(const std::string &name, ISvcLocator* pSvcLocator); + + /** destructor */ + ~MdtCalibrationT0ShiftSvc(); + + /* + * initalization of map cannot happen before first event + * special function required + */ + StatusCode initializeMap() override; + +private: +}; + +#endif diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationTMaxShiftSvc.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationTMaxShiftSvc.h new file mode 100644 index 0000000000000000000000000000000000000000..8550aa9ddf9e6d4455342ec36441f13e84afec7d --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/MdtCalibSvc/MdtCalibrationTMaxShiftSvc.h @@ -0,0 +1,48 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef MDTCALIBSVC_MDTCALIBRATIONTMAXSHIFTSVC_H +#define MDTCALIBSVC_MDTCALIBRATIONTMAXSHIFTSVC_H + +/* + * + * Author Andreas Hoenle (der.andi@cern.ch) + * + */ + +#include <map> + +#include "AthenaBaseComps/AthService.h" + +#include "MdtCalibSvc/MdtCalibrationShiftMapBase.h" + +class Identifier; + +/* + @class MdtCalibrationTMaxShiftSvc + Provides a per-tube shifting of the TMax value. + @author Andreas Hoenle +*/ +class MdtCalibrationTMaxShiftSvc : virtual public MdtCalibrationShiftMapBase { +public: + /* constructor */ + MdtCalibrationTMaxShiftSvc(const std::string &name, ISvcLocator* pSvcLocator); + + /* destructor */ + ~MdtCalibrationTMaxShiftSvc(); + + /* + * initalization of map cannot happen before first event + * special function required + */ + StatusCode initializeMap() override; + + StatusCode setTUpper(float tUpper); + float getTUpper() { return m_tUpper; } + +private: + float m_tUpper; +}; + +#endif diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/share/MdtCalibrationT0ShiftExample.py b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/share/MdtCalibrationT0ShiftExample.py new file mode 100644 index 0000000000000000000000000000000000000000..1ca4ef26f70604fbb46fcb604aea03ae494c9a59 --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/share/MdtCalibrationT0ShiftExample.py @@ -0,0 +1,20 @@ +# example configuration of MdtCalibT0ShiftSvc and MdtCalibTMaxShiftSvc +# in this example the 0ns shift maps are loaded, which means that the +# modification code runs, but it does not change anything + + +""" +NEVER recreate the shift maps when you're running multiple jobs. +If you want to recreate, run a single job to create the map and configure the +jobs to load the new map when submitting them. +It's all just random numbers and you cannot be sure that the tubes will get the +same shifts in all jobs. +(Yes, the RNGs are seeded, but what if someone changes this later?) +""" + +from MdtCalibT0ShiftSvc.MdtCalibT0ShiftSvcConf import MdtCalibrationT0ShiftSvc +ServiceMgr += MdtCalibrationT0ShiftSvc() +ServiceMgr.MdtCalibrationT0ShiftSvc.MapFile = 'shift_t0_0ns.dat' # which shift map to load / how to save it if not existing (in the package's share folder) +ServiceMgr.MdtCalibrationT0ShiftSvc.CentralValue = 0 # central value of Gaussian RNG - only relevant if you recreate the map +ServiceMgr.MdtCalibrationT0ShiftSvc.Sigma = 0 # sigma value of Gaussian RNG - only relevant if you recreate the map +ServiceMgr.MdtCalibrationT0ShiftSvc.ForceMapRecreate = False # force recreation of map - never do this when submitting multiple jobs diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationShiftMapBase.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationShiftMapBase.cxx new file mode 100644 index 0000000000000000000000000000000000000000..921757b230177c79d6abfe5f7e212c255dc42446 --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationShiftMapBase.cxx @@ -0,0 +1,175 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +// std +#include <sys/stat.h> +#include <cstring> +#include <fstream> +#include <iostream> +#include <string> + +// other packages +#include "GaudiKernel/GaudiException.h" +#include "GaudiKernel/MsgStream.h" +#include "Identifier/Identifier.h" +#include "MuonCablingData/MuonMDT_CablingMap.h" +#include "MuonIdHelpers/MdtIdHelper.h" +#include "PathResolver/PathResolver.h" +#include "StoreGate/DataHandle.h" +#include "StoreGate/StoreGateSvc.h" + +// this package +#include "MdtCalibSvc/MdtCalibrationShiftMapBase.h" + +#include "TTree.h" +// TRandom3 for random shifting +// should later be replaced with std +#include "TRandom3.h" + +// +// private helper functions +// + +MdtCalibrationShiftMapBase::MdtCalibrationShiftMapBase(const std::string &name, + ISvcLocator *pSvcLocator) + : base_class(name, pSvcLocator), + m_cablingSvc("MuonMDT_CablingSvc", name), + m_mapIsInitialized(false), + m_mapFileName(""), + m_centralValue(0), + m_sigma(10), + m_forceMapRecreate(false) { + declareProperty("MapFile", m_mapFileName); + declareProperty("CentralValue", m_centralValue = 0); + declareProperty("Sigma", m_sigma = 10); + declareProperty("ForceMapRecreate", m_forceMapRecreate = false); +} + +MdtCalibrationShiftMapBase::~MdtCalibrationShiftMapBase() { ; } + +// return failure if not overloaded +StatusCode MdtCalibrationShiftMapBase::initializeMap() { + return StatusCode::FAILURE; +} + +StatusCode MdtCalibrationShiftMapBase::dumpMapAsFile() { + /* initialize map if it's not there */ + if (!m_mapIsInitialized) { + initializeMap(); + } + + /* write the map to a file */ + { + std::ofstream file(m_mapFileName.c_str(), std::ios::out | std::ios::trunc); + + /* see if opening the file was successful */ + if (!file.is_open()) { + ATH_MSG_FATAL( + "Cannot open map output file for writing. Tried accessing file at " + "\"./" + << m_mapFileName.c_str() << "\""); + return StatusCode::FAILURE; + } + + /* write header, comments can start with '#' or '//' */ + file << "# This file contains shift values for MDT tubes\n"; + file << "# Each Identifier is mapped to a float\n"; + file << "# Below are comma separated values with formatting\n"; + file << "# Identifier.get_compact(),float\n"; + file << "# ------------------------------------------------\n"; + + /* dump map contents */ + for (auto shift : m_shiftValues) { + Identifier::value_type identifierCompact = shift.first.get_compact(); + float shiftValue = shift.second; + file << identifierCompact; + file << ","; + file << shiftValue; + file << "\n"; + } + } // '}' flushes file + return StatusCode::SUCCESS; +} + +StatusCode MdtCalibrationShiftMapBase::loadMapFromFile() { + /* check if map was already initialized */ + if (m_mapIsInitialized) { + ATH_MSG_WARNING("Map already initialized, won't overwrite it."); + return StatusCode::SUCCESS; + } + + /* resolve path */ + std::string fileWithPath = PathResolver::find_file(m_mapFileName, "DATAPATH"); + if (fileWithPath == "") { + ATH_MSG_ERROR("Cannot find file " << m_mapFileName << " in $DATAPATH"); + return StatusCode::FAILURE; + } + + /* check if the file exists */ + struct stat buffer; + if (stat(fileWithPath.c_str(), &buffer) != 0) { + ATH_MSG_ERROR("Cannot stat the file \"" + << fileWithPath.c_str() + << "\" -> map can not be initialized from this file."); + return StatusCode::FAILURE; + } + + /* check if the map is already stored */ + std::ifstream fin(fileWithPath.c_str(), std::ios::in); + std::string line; + bool initializedWithWarnings = false; + // get all lines in file + while (std::getline(fin, line)) { + // check if file is empty, begins with '#', or '//' + if (line.empty() || line.compare(0, 1, "#") == 0 || + line.compare(0, 2, "//")) { + continue; + } + // need a stringstream for readline + std::stringstream lineStream(line); + // get all csv tokens in line + std::string token; + std::vector<std::string> tokenVector; + while (std::getline(lineStream, token, ',')) { + tokenVector.push_back(token); + } + // expect exactly two tokens: identifier and value + if (tokenVector.size() == 2) { + // we are careful about the type here + // watch out for compiler warnings warning about casting ull to value_type + // they might occur of Identifier::value_type is changed at some point + Identifier::value_type identifierCompact = std::stoull(tokenVector[0]); + Identifier id(identifierCompact); + float shift = std::stof(tokenVector[1]); + m_shiftValues[id] = shift; + } else { + ATH_MSG_WARNING("Unexpected input format in shift map file " + << fileWithPath.c_str()); + ATH_MSG_WARNING("Expected 2 tokens, got " << tokenVector.size()); + ATH_MSG_WARNING("Broken line: \"" << line.c_str() << "\""); + initializedWithWarnings = true; + } + } + m_mapIsInitialized = true; + if (initializedWithWarnings) { + ATH_MSG_WARNING("Initialized shift map WITH WARNINGS from file \"" + << fileWithPath.c_str() << "\""); + } else { + ATH_MSG_INFO("Successfully initialized shift map from file \"" + << fileWithPath.c_str() << "\""); + } + return StatusCode::SUCCESS; +} + +float MdtCalibrationShiftMapBase::getValue(const Identifier &id) { + if (!m_mapIsInitialized) { + StatusCode sc = initializeMap(); + if (sc.isFailure()) { + ATH_MSG_FATAL("Could not initialize the shift map."); + throw GaudiException("Cannot run without shift map.", "", sc); + } + } + + return m_shiftValues[id]; +} diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationSvc.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationSvc.cxx index 8db26ab5922e1b87a258295efee174e5baf3c7b8..1f729502aac13b3a3480838bdb9fecced9250a43 100644 --- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationSvc.cxx +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationSvc.cxx @@ -38,6 +38,8 @@ #include "MdtCalibData/TrRelation.h" #include "MdtCalibData/RtScaleFunction.h" +#include "MdtCalibInterfaces/IShiftMapTools.h" + #include "GaudiKernel/ServiceHandle.h" // @@ -46,7 +48,7 @@ class MdtCalibrationSvc::Imp { public: - + Imp( std::string name ); const MuonGM::MdtReadoutElement* getGeometry( const Identifier &id ) { @@ -64,10 +66,10 @@ public: const MuonGM::MuonDetectorManager *m_muonGeoManager; const MdtIdHelper *m_mdtIdHelper; MdtCalibrationSvcSettings settings; - + double m_inverseSpeedOfLight; // in ns/mm double m_inversePropagationSpeed; // in ns/mm - + MsgStream *m_log; //!< Pointer to MsgStream - MdtCalibrationSvc now owns this bool m_verbose; //!< True if we should print MSG::VERBOSE messages bool m_debug; //!< True if we should print MSG::DEBUG messages @@ -75,7 +77,16 @@ public: //handle of b-field service - job option ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc; ServiceHandle<MdtCalibrationDbSvc> m_dbSvc; - + + /* T0 Shift Service -- Per-tube offsets of t0 value */ + ServiceHandle<MuonCalib::IShiftMapTools> m_t0ShiftSvc; + /* TMax Shift Service -- Per-tube offsets of Tmax */ + ServiceHandle<MuonCalib::IShiftMapTools> m_tMaxShiftSvc; + + // tools should only be retrieved if they are used + bool m_doT0Shift; + bool m_doTMaxShift; + double m_unphysicalHitRadiusUpperBound; double m_unphysicalHitRadiusLowerBound; double m_resTwin; @@ -93,6 +104,10 @@ MdtCalibrationSvc::Imp::Imp(std::string name) : m_debug(false), m_magFieldSvc("AtlasFieldSvc", name), m_dbSvc("MdtCalibrationDbSvc", name), + m_t0ShiftSvc("MdtCalibrationT0ShiftSvc", name), + m_tMaxShiftSvc("MdtCalibrationTMaxShiftSvc", name), + m_doT0Shift(false), + m_doTMaxShift(false), m_unphysicalHitRadiusUpperBound(-1.), m_unphysicalHitRadiusLowerBound(-1.), m_resTwin(-1.), @@ -118,35 +133,37 @@ MdtCalibrationSvc::MdtCalibrationSvc(const std::string &name,ISvcLocator *sl) declareProperty("MagFieldSvc", m_imp->m_magFieldSvc); declareProperty("UpperBoundHitRadius", m_imp->m_unphysicalHitRadiusUpperBound = 20. ); declareProperty("LowerBoundHitRadius", m_imp->m_unphysicalHitRadiusLowerBound = 0. ); + declareProperty("DoT0Shift", m_imp->m_doT0Shift = false ); + declareProperty("DoTMaxShift", m_imp->m_doTMaxShift = false ); } MdtCalibrationSvc::~MdtCalibrationSvc() { delete m_imp; } -// queryInterface -StatusCode MdtCalibrationSvc::queryInterface(const InterfaceID &riid, void **ppvIF) { - if ( interfaceID().versionMatch(riid) ) { - *ppvIF = dynamic_cast<MdtCalibrationSvc*>(this); - } else { - return AthService::queryInterface(riid, ppvIF); +// queryInterface +StatusCode MdtCalibrationSvc::queryInterface(const InterfaceID &riid, void **ppvIF) { + if ( interfaceID().versionMatch(riid) ) { + *ppvIF = dynamic_cast<MdtCalibrationSvc*>(this); + } else { + return AthService::queryInterface(riid, ppvIF); } addRef(); return StatusCode::SUCCESS; } -StatusCode MdtCalibrationSvc::initialize() { - ATH_MSG_DEBUG( "Initializing" ); +StatusCode MdtCalibrationSvc::initialize() { + ATH_MSG_DEBUG( "Initializing" ); if( AthService::initialize().isFailure() ) return StatusCode::FAILURE; m_imp->settings.initialize(); // initialize the MdtIdHelper ServiceHandle<StoreGateSvc> detStore("StoreGateSvc/DetectorStore", name()); if ( detStore.retrieve().isFailure() ) { - ATH_MSG_FATAL( "Can't locate the DetectorStore" ); + ATH_MSG_FATAL( "Can't locate the DetectorStore" ); return StatusCode::FAILURE; } - + if ( detStore->retrieve(m_imp->m_mdtIdHelper, "MDTIDHELPER" ).isFailure() ) { ATH_MSG_FATAL( "Can't retrieve MdtIdHelper" ); return StatusCode::FAILURE; @@ -162,7 +179,7 @@ StatusCode MdtCalibrationSvc::initialize() { ATH_MSG_FATAL( "Can't retrieve MuonDetectorManager" ); return StatusCode::FAILURE; } - + if( m_imp->m_dbSvc.retrieve().isFailure() ) { ATH_MSG_FATAL( "Could not get MdtCalibrationDbSvc" ); return StatusCode::FAILURE; @@ -171,14 +188,18 @@ StatusCode MdtCalibrationSvc::initialize() { if ( m_imp->m_magFieldSvc.retrieve().isFailure() ) { ATH_MSG_FATAL( "Could not find MagFieldSvc" ); return StatusCode::FAILURE; - } + } + + if (m_imp->m_doT0Shift) ATH_CHECK( m_imp->m_t0ShiftSvc.retrieve() ); + + if (m_imp->m_doTMaxShift) ATH_CHECK ( m_imp->m_tMaxShiftSvc.retrieve() ); m_imp->m_log = &msg(); m_imp->m_verbose = msgLvl(MSG::VERBOSE); m_imp->m_debug = msgLvl(MSG::DEBUG); ATH_MSG_DEBUG("Settings:"); - ATH_MSG_DEBUG(" Time window: lower bound " << m_imp->settings.timeWindowLowerBound() + ATH_MSG_DEBUG(" Time window: lower bound " << m_imp->settings.timeWindowLowerBound() << " upper bound " << m_imp->settings.timeWindowUpperBound() ); if( m_imp->settings.doTofCorrection() ) ATH_MSG_DEBUG(" Using TOF "); if( m_imp->settings.doPropCorrection() ) ATH_MSG_DEBUG(" Using Prop "); @@ -192,7 +213,7 @@ StatusCode MdtCalibrationSvc::initialize() { return StatusCode::SUCCESS; } //end MdtCalibrationSvc::initialize -StatusCode MdtCalibrationSvc::finalize() { +StatusCode MdtCalibrationSvc::finalize() { ATH_MSG_DEBUG( "Finalizing " ); return AthService::finalize(); } @@ -206,7 +227,7 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, inputData.triggerOffset = triggerTime; return driftRadiusFromTime( hit, inputData, m_imp->settings, resolFromRtrack ); } - + bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, const MdtCalibrationSvcInput &inputData, bool resolFromRtrack ) { @@ -217,13 +238,13 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, const MdtCalibrationSvcInput &inputData, const MdtCalibrationSvcSettings &settings, bool resolFromRtrack ) { - + if( settings.timeWindowUpperBound() < 0. || settings.timeWindowLowerBound() < 0. ) { // Should be an ERROR, but can't return StatusCode::FAILURE ATH_MSG_WARNING( "Uninitialized settings object at " << &settings ); return false; } - + // Easiest MessageSvc implementation, not necessarily the best if (msgLvl(MSG::VERBOSE)) { msg(MSG::VERBOSE) << "Input: tof " << inputData.tof << " trigger offset " << inputData.triggerOffset; @@ -256,8 +277,8 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, // get calibration constants from DbSvc MuonCalib::MdtFullCalibData data = m_imp->m_dbSvc->getCalibration( geo->collectionHash(), - geo->detectorElementHash() ); - + geo->detectorElementHash() ); + // require at least the MdtRtRelation to be available const MuonCalib::MdtRtRelation *rtRelation = data.rtRelation; // Hardcoded MDT tube radius 14.6mm here - not correct for sMDT @@ -277,14 +298,14 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, const int ml = m_imp->m_mdtIdHelper->multilayer(id)-1; const int layer = m_imp->m_mdtIdHelper->tubeLayer(id)-1; const int tube = m_imp->m_mdtIdHelper->tube(id)-1; - + if ( ml<0 || layer<0 || tube<0 ){ ATH_MSG_WARNING( "Oops negative index....." ); return false; } - + // extract calibration constants for single tube - const MuonCalib::MdtTubeCalibContainer::SingleTubeCalib *singleTubeData = + const MuonCalib::MdtTubeCalibContainer::SingleTubeCalib *singleTubeData = data.tubeCalib->getCalib( ml, layer, tube ); if( singleTubeData ){ t0 = singleTubeData->t0; @@ -294,12 +315,15 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, msg(MSG::WARNING) << "failed to access tubedata for " << ml << " " << layer << " " << tube << " using defaults.. " << endmsg; - if ( geo ) + if ( geo ) msg(MSG::WARNING) << "detel " << geo->getMultilayer() - << " lay " << geo->getNLayers() + << " lay " << geo->getNLayers() << " tubes " << geo->getNtubesperlayer() << endmsg; t0 = 800.; } + + // get t0 shift from tool (default: no shift, value is zero) + if (m_imp->m_doT0Shift) t0 += m_imp->m_t0ShiftSvc->getValue(id); } else { msg(MSG::WARNING) << "MdtTubeCalibContainer not found for " << m_imp->m_mdtIdHelper->print_to_string( id ) << endmsg; @@ -312,7 +336,7 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, hit.setTubeT0(t0); hit.setTubeAdcCal(adcCal); - // in order to have clean info in the ntuple set to 0 the + // in order to have clean info in the ntuple set to 0 the // corrections which are not used /* hit.setSlewingTime( 0 ); hit.setLorentzTime( 0 ); @@ -331,13 +355,13 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, // set time-of-flight double triggerTime = inputData.tof + inputData.triggerOffset; - hit.setTimeOfFlight( settings.doTof ? triggerTime : 0. ); + hit.setTimeOfFlight( settings.doTof ? triggerTime : 0. ); // calculate drift time double driftTime = hit.tdcCount() * tdcBinSize(id) - hit.timeOfFlight() - hit.tubeT0() - hit.propagationTime(); hit.setDriftTime( driftTime ); - + // apply corrections double corTime(0.); if ( settings.doCorrections() ) { @@ -353,35 +377,42 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, bool calibOk = true; Muon::MdtDriftCircleStatus timeStatus = driftTimeStatus(t, rtRelation, settings); if( rtRelation->rt() ){ - + + r = rtRelation->rt()->radius(t); + + // apply tUpper gshift + if (m_imp->m_doTMaxShift) { + float tShift = m_imp->m_tMaxShiftSvc->getValue(id); + r = rtRelation->rt()->radius( t * (1 + tShift) ); + } + // check whether drift times are within range, if not fix them to the min/max range - r = rtRelation->rt()->radius( t ); if ( t < rtRelation->rt()->tLower() ) { t_inrange = rtRelation->rt()->tLower(); double rmin = rtRelation->rt()->radius( t_inrange ); double drdt = (rtRelation->rt()->radius( t_inrange + 30. ) - rmin)/30.; - + // now check whether we are outside the time window if ( timeStatus == Muon::MdtStatusBeforeSpectrum ) { t = rtRelation->rt()->tLower() - settings.timeWindowLowerBound(); calibOk = false; } - // if we get here we are outside the rt range but inside the window. + // if we get here we are outside the rt range but inside the window. r = rmin + drdt*(t-t_inrange); if( r < m_imp->m_unphysicalHitRadiusLowerBound ) r = m_imp->m_unphysicalHitRadiusLowerBound; } else if( t > rtRelation->rt()->tUpper() ) { t_inrange = rtRelation->rt()->tUpper(); double rmax = rtRelation->rt()->radius( t_inrange ); double drdt = (rmax - rtRelation->rt()->radius( t_inrange - 30. ))/30.; - + // now check whether we are outside the time window if ( timeStatus == Muon::MdtStatusAfterSpectrum ) { t = rtRelation->rt()->tUpper() + settings.timeWindowUpperBound(); calibOk = false; } - // if we get here we are outside the rt range but inside the window. + // if we get here we are outside the rt range but inside the window. r = rmax + drdt*(t-t_inrange); } } else { @@ -402,19 +433,19 @@ bool MdtCalibrationSvc::driftRadiusFromTime( MdtCalibHit &hit, hit.setDriftRadius( r, reso ); // summary - ATH_MSG_VERBOSE( "driftRadiusFromTime for tube " << m_imp->m_mdtIdHelper->print_to_string(id) + ATH_MSG_VERBOSE( "driftRadiusFromTime for tube " << m_imp->m_mdtIdHelper->print_to_string(id) << (calibOk ? " OK" : "FAILED") ); ATH_MSG_VERBOSE( " raw drift time " << hit.tdcCount() * tdcBinSize(id) << " TriggerOffset " << inputData.triggerOffset << endmsg << "Tof " << inputData.tof << " Propagation Delay " - << hit.propagationTime() << " T0 " << hit.tubeT0() + << hit.propagationTime() << " T0 " << hit.tubeT0() << " invProp " << inversePropSpeed << endmsg << "Drift time after corrections " << driftTime - << " time cor " << corTime + << " time cor " << corTime << " Drift radius " << hit.driftRadius() << " sigma " << hit.sigmaDriftRadius() ); - + return calibOk; } //end MdtCalibrationSvc::driftRadiusFromTime @@ -428,7 +459,7 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, MdtCalibrationSvcInput inputData; inputData.tof = signedTrackLength*m_imp->m_inverseSpeedOfLight; inputData.triggerOffset = triggerTime; - + MdtCalibrationSvcInput secondInputData; secondInputData.tof = secondSignedTrackLength*m_imp->m_inverseSpeedOfLight; secondInputData.triggerOffset = triggerTime; @@ -443,12 +474,12 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, const MdtCalibrationSvcSettings &settings, bool &secondDigitIsPrompt ) { - // 13/02/2009 A.Koutsman: after discussion with P.Kluit rebuilt this function to use the standard way - // of calibrating a MdtCalibHit with driftRadiusFromTime(...) - + // 13/02/2009 A.Koutsman: after discussion with P.Kluit rebuilt this function to use the standard way + // of calibrating a MdtCalibHit with driftRadiusFromTime(...) + driftRadiusFromTime( hit, inputData, settings ); driftRadiusFromTime( secondHit, secondInputData, settings ); - + // get Identifier and MdtReadOutElement for twin tubes const Identifier &id = hit.identify(); const Identifier &idSecond = secondHit.identify(); @@ -458,7 +489,7 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, // get 'raw' drifttimes of twin pair; we don't use timeofFlight or propagationTime cause they are irrelevant for twin coordinate double driftTime = hit.tdcCount()*tdcBinSize(id) - hit.tubeT0(); double driftTimeSecond = secondHit.tdcCount()*tdcBinSize(idSecond) - secondHit.tubeT0(); - + if(!geo) { ATH_MSG_WARNING( "Geometry not set for first hit" ); return false; @@ -468,31 +499,31 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, return false; } // get calibration constants from DbSvc - MuonCalib::MdtFullCalibData data = m_imp->m_dbSvc->getCalibration( geo->collectionHash(), geo->detectorElementHash() ); - MuonCalib::MdtFullCalibData dataSecond = m_imp->m_dbSvc->getCalibration( geoSecond->collectionHash(), geoSecond->detectorElementHash() ); - + MuonCalib::MdtFullCalibData data = m_imp->m_dbSvc->getCalibration( geo->collectionHash(), geo->detectorElementHash() ); + MuonCalib::MdtFullCalibData dataSecond = m_imp->m_dbSvc->getCalibration( geoSecond->collectionHash(), geoSecond->detectorElementHash() ); + //double t0(0.); double inversePropSpeed = m_imp->m_inversePropagationSpeed; //double adcCal(1.); - + //double t0Second(0.); double inversePropSpeedSecond = m_imp->m_inversePropagationSpeed; //double adcCalSecond(1.); - + // access t0 for the give tube if( data.tubeCalib ){ const int ml = m_imp->m_mdtIdHelper->multilayer(id)-1; const int layer = m_imp->m_mdtIdHelper->tubeLayer(id)-1; const int tube = m_imp->m_mdtIdHelper->tube(id)-1; - + if( ml<0 || layer<0 || tube<0 ){ ATH_MSG_WARNING( "Oops negative index....." ); return false; } - + // extract calibration constants for single tube - const MuonCalib::MdtTubeCalibContainer::SingleTubeCalib *singleTubeData = + const MuonCalib::MdtTubeCalibContainer::SingleTubeCalib *singleTubeData = data.tubeCalib->getCalib( ml, layer, tube ); if ( singleTubeData ){ //t0 = singleTubeData->t0; @@ -500,11 +531,11 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, //adcCal = singleTubeData->adcCal; } else { ATH_MSG_WARNING( "failed to access tubedata for " - << ml << " " << layer << " " << tube + << ml << " " << layer << " " << tube << " using defaults.. " ); - if ( geo ) + if ( geo ) ATH_MSG_WARNING( "detel " - << geo->getMultilayer() << " lay " << geo->getNLayers() + << geo->getMultilayer() << " lay " << geo->getNLayers() << " tubes " << geo->getNtubesperlayer() ); } } else { @@ -513,20 +544,20 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, ATH_MSG_WARNING( "Tube cannot be calibrated!!!" ); return false; } - + // access t0 for the given second tube if ( dataSecond.tubeCalib ){ const int mlSecond = m_imp->m_mdtIdHelper->multilayer(idSecond)-1; const int layerSecond = m_imp->m_mdtIdHelper->tubeLayer(idSecond)-1; const int tubeSecond = m_imp->m_mdtIdHelper->tube(idSecond)-1; - + if( mlSecond<0 || layerSecond<0 || tubeSecond<0 ){ ATH_MSG_WARNING( "Oops negative index for second tube....." ); return false; } - + // extract calibration constants for single tube - const MuonCalib::MdtTubeCalibContainer::SingleTubeCalib *singleTubeDataSecond = + const MuonCalib::MdtTubeCalibContainer::SingleTubeCalib *singleTubeDataSecond = dataSecond.tubeCalib->getCalib( mlSecond, layerSecond, tubeSecond ); if ( singleTubeDataSecond ){ //t0Second = singleTubeDataSecond->t0; @@ -534,11 +565,11 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, //adcCalSecond = singleTubeDataSecond->adcCal; } else { ATH_MSG_WARNING( "failed to access (second) tubedata for " << mlSecond - << " " << layerSecond << " " << tubeSecond + << " " << layerSecond << " " << tubeSecond << " using defaults.. " ); - if ( geoSecond ) + if ( geoSecond ) ATH_MSG_WARNING( "detel " << geoSecond->getMultilayer() - << " lay " << geoSecond->getNLayers() + << " lay " << geoSecond->getNLayers() << " tubes " << geoSecond->getNtubesperlayer() ); //t0Second = 800.; } @@ -563,20 +594,20 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, const MuonGM::MdtReadoutElement *prompthit_geo = geo; int prompthit_tdc = 0; //double prompthit_adc = 0; int twinhit_tdc = 0; //double twinhit_adc = 0; - if ( driftTime < driftTimeSecond) { - twin_timedif = driftTimeSecond - driftTime; - //prompthit_tdc = hit.tdcCount(); prompthit_adc = adcCal; - //twinhit_tdc = secondHit.tdcCount(); twinhit_adc = adcCalSecond; + if ( driftTime < driftTimeSecond) { + twin_timedif = driftTimeSecond - driftTime; + //prompthit_tdc = hit.tdcCount(); prompthit_adc = adcCal; + //twinhit_tdc = secondHit.tdcCount(); twinhit_adc = adcCalSecond; secondDigitIsPrompt = false; - } else if (driftTimeSecond <= driftTime){ - prompthit_id = idSecond; - prompthit_geo = geoSecond; - twin_timedif = driftTime - driftTimeSecond; - //prompthit_tdc = secondHit.tdcCount(); prompthit_adc = adcCalSecond; - //twinhit_tdc = hit.tdcCount(); twinhit_adc = adcCal; + } else if (driftTimeSecond <= driftTime){ + prompthit_id = idSecond; + prompthit_geo = geoSecond; + twin_timedif = driftTime - driftTimeSecond; + //prompthit_tdc = secondHit.tdcCount(); prompthit_adc = adcCalSecond; + //twinhit_tdc = hit.tdcCount(); twinhit_adc = adcCal; secondDigitIsPrompt = true; } - + // get tubelength and set HV-delay (~6ns) double tubelength = prompthit_geo->tubeLength(prompthit_id); @@ -585,22 +616,22 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, // twin_timedif must be between min and max of possible time-difference // between prompt and twin signals // accounting for 3 std.dev. of twin time resolution - if ( twin_timedif < (HVdelay - 5.*m_imp->m_resTwin) + if ( twin_timedif < (HVdelay - 5.*m_imp->m_resTwin) || twin_timedif > (tubelength*inversePropSpeed + tubelength*inversePropSpeedSecond + HVdelay + 5.*m_imp->m_resTwin)){ if ( m_imp->m_debug ) { - ATH_MSG_DEBUG( " TIME DIFFERENCE OF TWIN PAIR OUT OF RANGE(" - << (HVdelay - 5.*m_imp->m_resTwin) + ATH_MSG_DEBUG( " TIME DIFFERENCE OF TWIN PAIR OUT OF RANGE(" + << (HVdelay - 5.*m_imp->m_resTwin) << "-" << (tubelength*inversePropSpeed + tubelength*inversePropSpeedSecond - + HVdelay + 5.*m_imp->m_resTwin) + + HVdelay + 5.*m_imp->m_resTwin) << ") time difference = " << twin_timedif ); - } + } } - + // Make ONLY a twin PrepData if twin time difference is physical (within tubelength) if (twin_timedif < (tubelength*inversePropSpeed + tubelength*inversePropSpeedSecond @@ -617,22 +648,22 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, ATH_MSG_DEBUG( " TWIN HIT outside acceptance with time difference " << twin_timedif << " Z local hit " << z_hit_sign_from_twin - << " Z local minimum " << -tubelength/2 ); + << " Z local minimum " << -tubelength/2 ); } - z_hit_sign_from_twin = - tubelength/2.; + z_hit_sign_from_twin = - tubelength/2.; } // do sign management just like in MdtDigitizationTool.cxx - double distRO = prompthit_geo->tubeFrame_localROPos(prompthit_id).z(); - double sign(-1.); + double distRO = prompthit_geo->tubeFrame_localROPos(prompthit_id).z(); + double sign(-1.); if (distRO < 0.) sign = 1.; double z_hit_geo_from_twin = sign*z_hit_sign_from_twin; - + zTwin = z_hit_geo_from_twin; errZTwin = m_imp->m_resTwin*inversePropSpeed; - + if ( m_imp->m_verbose ){ - ATH_MSG_VERBOSE( " TWIN TUBE " + ATH_MSG_VERBOSE( " TWIN TUBE " << " tube: " << m_imp->m_mdtIdHelper->print_to_string(id) << " twintube: " << m_imp->m_mdtIdHelper->print_to_string(idSecond) ); @@ -648,16 +679,16 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, << " distRO = " << distRO ); } - + } // end if(twin_timedif < (tubelength*inversePropSpeed + tubelength*inversePropSpeedSecond + HVdelay + 10.*m_imp->m_resTwin)){ else { - + if ( m_imp->m_verbose ){ - ATH_MSG_VERBOSE( " TIME DIFFERENCE OF TWIN PAIR UNPHYSICAL OUT OF RANGE(" + ATH_MSG_VERBOSE( " TIME DIFFERENCE OF TWIN PAIR UNPHYSICAL OUT OF RANGE(" << (HVdelay - 5*m_imp->m_resTwin) << "-" - << (2*tubelength*inversePropSpeed + HVdelay + 5*m_imp->m_resTwin) + << (2*tubelength*inversePropSpeed + HVdelay + 5*m_imp->m_resTwin) << ") time difference = " - << twin_timedif ); + << twin_timedif ); } zTwin = 0.; errZTwin = tubelength/2.; @@ -674,17 +705,17 @@ bool MdtCalibrationSvc::twinPositionFromTwinHits( MdtCalibHit &hit, double MdtCalibrationSvc::tdcBinSize(const Identifier &id) { //BMG which uses HPTDC instead of AMT, and has 0.2ns TDC ticksize //if( m_imp->m_mdtIdHelper->stationName(id) == 54 ) //BMG - //return 0.2; + //return 0.2; // Alternative method if you don't like hardcoding BMG stationName (54) // don't do string comparisons! //if( m_imp->m_mdtIdHelper->stationNameString( m_imp->m_mdtIdHelper->stationName(id) ) == "BMG" ) if( m_imp->m_mdtIdHelper->stationName(id) == m_imp->m_BMGid && m_imp->m_BMGpresent) return 0.1953125; // 25/128 - return 0.78125; //25/32; exact number: (1000.0/40.079)/32.0 + return 0.78125; //25/32; exact number: (1000.0/40.079)/32.0 } double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, - const MdtCalibrationSvcInput &inputData, + const MdtCalibrationSvcInput &inputData, const MdtCalibrationSvcSettings &settings, double /*adcCal*/, const MuonCalib::MdtCorFuncSet *corrections, const MuonCalib::IRtRelation *rt) const { @@ -692,25 +723,25 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, double corTime(0.); // apply corrections if ( corrections ){ - + if (m_verbose) *m_log << MSG::VERBOSE << "There are correction functions." << endmsg; - + // slewing corrections if ( settings.doSlew && corrections->slewing() ){ double slew_time=corrections->slewing()->correction( hit.driftTime(), hit.adcCount()); corTime -=slew_time; hit.setSlewingTime(slew_time); } - + if( settings.doField && corrections->bField()){ if( inputData.trackDirection && inputData.pointOfClosestApproach ) { Amg::Transform3D gToStation = hit.geometry()->GlobalToAmdbLRSTransform(); double BGXYZ[3]; - double xyz[3] = { (*inputData.pointOfClosestApproach)(0), - (*inputData.pointOfClosestApproach)(1), + double xyz[3] = { (*inputData.pointOfClosestApproach)(0), + (*inputData.pointOfClosestApproach)(1), (*inputData.pointOfClosestApproach)(2) }; m_magFieldSvc->getField(xyz, BGXYZ); - for (int i=0; i<3; i++) BGXYZ[i] *= 1000.; // convert kT to Tesla + for (int i=0; i<3; i++) BGXYZ[i] *= 1000.; // convert kT to Tesla Amg::Vector3D B_global(BGXYZ[0], BGXYZ[1], BGXYZ[2]); Amg::Vector3D B_loc(gToStation.linear()*B_global); double Bpar = B_loc.x(); @@ -718,7 +749,7 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, Amg::Vector3D dir(0.0, loc_dir.y(), loc_dir.z()); double Bper = B_loc.dot(dir.unit()); hit.setBFieldPerp(Bper); - hit.setBFieldPara(Bpar); + hit.setBFieldPara(Bpar); if( m_verbose ) *m_log << MSG::VERBOSE << "doing b-field correction" << endmsg; if(hit.bFieldPara() != MdtCalibHit::kNoValue && hit.bFieldPerp() != MdtCalibHit::kNoValue) { hit.setLorentzTime(corrections->bField()->correction( hit.driftTime(), hit.bFieldPara(), hit.bFieldPerp() )); @@ -726,11 +757,11 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, hit.setLorentzTime(0); } corTime -= hit.lorentzTime(); - } - } - + } + } + // temperature corrections - // NOTE: Use this temporarily for ML dependent scaling. + // NOTE: Use this temporarily for ML dependent scaling. /* if ( settings.doTemp && corrections->temperature() ){ if(hit.temperature() != MdtCalibHit::kNoValue) { @@ -742,7 +773,7 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, } corTime += hit.TemperatureTime(); }*/ - + // Scale RT function from Tmax difference if( settings.doTemp && rt && rt->HasTmaxDiff()) { float scle_time=MuonCalib::RtScaleFunction(hit.driftTime(), m_mdtIdHelper->multilayer(hit.identify())==2, *rt); @@ -751,13 +782,13 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, } else { hit.setTemperatureTime(0); } - + // background corrections if ( settings.doBkg && corrections->background() ){ double bgLevel = 0.; corTime += corrections->background()->correction( hit.driftTime(), bgLevel ); } - + // wire sag corrections // First some debug output if (m_verbose) { @@ -774,10 +805,10 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, } // Wire sag corrections if ( settings.doWireSag && corrections->wireSag() ){ - + if (m_verbose) { *m_log << MSG::VERBOSE << "Performing Rt Corrections due to Wire Sag." << endmsg; - + if (inputData.pointOfClosestApproach) { *m_log << MSG::VERBOSE << "Have a point of closest approach." << endmsg; } else { @@ -799,12 +830,12 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, *m_log << MSG::VERBOSE << "No sagged wire surface!" << endmsg; } } - + const Amg::Vector3D *pointOfClosestApproach = inputData.pointOfClosestApproach; const Amg::Vector3D *trackDirection = inputData.trackDirection; const Trk::StraightLineSurface *nominalWireSurface = inputData.nominalWireSurface; const Trk::StraightLineSurface *wireSurface = inputData.wireSurface; - + // check whether all input data is there if ( !pointOfClosestApproach || !trackDirection || !nominalWireSurface ){ *m_log << MSG::WARNING << " cannot perform wire sag correction: "; @@ -825,7 +856,7 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, // store pointer to sagged surface as we get ownership if we recalculate it const Trk::StraightLineSurface *tempSaggedSurface = 0; - // if wire surface is missing, calculate sagged wire position + // if wire surface is missing, calculate sagged wire position if ( !wireSurface ){ const Trk::SaggedLineSurface *nominalSurf = dynamic_cast<const Trk::SaggedLineSurface*>(nominalWireSurface); if ( nominalSurf ){ @@ -856,7 +887,7 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, // get to transformation matrix from global into the tube frame Amg::Transform3D globalToTube = nominalWireSurface->transform().inverse(); - // calculate the local position of the nominal and sagged wire + // calculate the local position of the nominal and sagged wire Amg::Vector3D locNominalPos(0.,0.,0.); Amg::Vector3D locSaggedPos = globalToTube*wireSurface->center(); locSaggedPos[2] = 0.; @@ -888,7 +919,7 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, // above/below wire double signedDriftRadius = fabs(hit.driftRadius())*deltaY/fabs(deltaY); - + // calculate the magnitude of the wire sag double effectiveSag = nominalWireSurface->center().y() @@ -912,11 +943,11 @@ double MdtCalibrationSvc::Imp::applyCorrections(MdtCalibHit &hit, Muon::MdtDriftCircleStatus MdtCalibrationSvc::driftTimeStatus( double driftTime, const MuonCalib::MdtRtRelation *rtRelation, const MdtCalibrationSvcSettings &settings ) { if( settings.timeWindowUpperBound() < 0. || settings.timeWindowLowerBound() < 0. ) { - ATH_MSG_WARNING( " Unphysical time window detected, both ranges should be positive: " + ATH_MSG_WARNING( " Unphysical time window detected, both ranges should be positive: " << " lower bound " << settings.timeWindowLowerBound() << " upper bound " << settings.timeWindowUpperBound() ); } - + if( rtRelation && rtRelation->rt() ) { if( driftTime < rtRelation->rt()->tLower() - settings.timeWindowLowerBound() ) { // Will this produce too much output from ROT creators? @@ -934,7 +965,7 @@ Muon::MdtDriftCircleStatus MdtCalibrationSvc::driftTimeStatus( double driftTime, } } else { ATH_MSG_WARNING( "No valid rt relation supplied for driftTimeStatus method" ); - return Muon::MdtStatusUnDefined; + return Muon::MdtStatusUnDefined; } return Muon::MdtStatusDriftTime; } //end MdtCalibrationSvc::driftTimeStatus diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationT0ShiftSvc.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationT0ShiftSvc.cxx new file mode 100644 index 0000000000000000000000000000000000000000..63f48bc47728b8f94f5a875e9a57d6d0ec3c52c9 --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationT0ShiftSvc.cxx @@ -0,0 +1,175 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +// std +#include <fstream> +#include <string> +#include <iostream> +#include <cstring> +#include <sys/stat.h> + +// other packages +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/GaudiException.h" +#include "Identifier/Identifier.h" +#include "StoreGate/DataHandle.h" +#include "StoreGate/StoreGateSvc.h" +#include "MuonCablingData/MuonMDT_CablingMap.h" +#include "MuonIdHelpers/MdtIdHelper.h" + +// this package +#include "TTree.h" +#include "MdtCalibSvc/MdtCalibrationShiftMapBase.h" +#include "MdtCalibSvc/MdtCalibrationT0ShiftSvc.h" + +// TRandom3 for random smearing +// should later be replaced with std +#include "TRandom3.h" + +// +// private helper functions +// + +MdtCalibrationT0ShiftSvc::MdtCalibrationT0ShiftSvc(const std::string &name, + ISvcLocator* pSvcLocator) + : MdtCalibrationShiftMapBase(name, pSvcLocator) {} + +MdtCalibrationT0ShiftSvc::~MdtCalibrationT0ShiftSvc() {} + +StatusCode MdtCalibrationT0ShiftSvc::initializeMap() +{ + if (m_mapIsInitialized) { + ATH_MSG_WARNING("Map already initalized. Multiple calls of initalizeMap should not happen."); + return StatusCode::SUCCESS; + } + + /* First try loading the map from its default location */ + if (!m_forceMapRecreate) { + ATH_CHECK(loadMapFromFile()); + return StatusCode::SUCCESS; + } + + /* map was not found as a file or needs to be recreated */ + // TODO: What if multiple jobs started in parallel? + // (DataRace on MapFile because of parallel jobs) + + /* initialize random number generator that creates the shift values */ + TRandom3 rng(/*seed*/ 20120704); + + /* idHelper to retrieve channel Identifiers */ + const MdtIdHelper *idhelper = nullptr; + const ServiceHandle<StoreGateSvc> detStore("StoreGateSvc/DetectorStore", "detStore"); + ATH_CHECK( detStore->retrieve(idhelper, "MDTIDHELPER") ); + + + // Cabling Map (part of cablingSvc) must be retrieved AFTER the first event + // i.e. this can NOT happen in initialize + /* initialize cabling from which IDs are read */ + ATH_CHECK( m_cablingSvc.retrieve() ); + + /* Get ROBs */ + std::vector<uint32_t> robVector = m_cablingSvc->getAllROBId(); + + std::map<uint8_t, MdtSubdetectorMap*> *listOfSubdet; + std::map<uint8_t, MdtSubdetectorMap*>::const_iterator it_sub; + + std::map<uint8_t, MdtRODMap*> *listOfROD; + std::map<uint8_t, MdtRODMap*>::const_iterator it_rod; + + std::map<uint8_t, MdtCsmMap*> *listOfCsm; + std::map<uint8_t, MdtCsmMap*>::const_iterator it_csm; + + std::map<uint8_t, MdtAmtMap*> *listOfAmt; + std::map<uint8_t, MdtAmtMap*>::const_iterator it_amt; + + DataHandle<MuonMDT_CablingMap> cablingMap = m_cablingSvc->getCablingMap(); + listOfSubdet = cablingMap->getListOfElements(); + + int subdetectorId, rodId, csmId, amtId; + int stationName, stationEta, stationPhi, multiLayer, layer, tube; + + for (it_sub = listOfSubdet->begin(); it_sub != listOfSubdet->end(); ++it_sub) { + subdetectorId = it_sub->first; + listOfROD = it_sub->second->getListOfElements(); + + for (it_rod = listOfROD->begin(); it_rod != listOfROD->end(); ++it_rod) { + rodId = it_rod->first; + listOfCsm = it_rod->second->getListOfElements(); + + for (it_csm = listOfCsm->begin(); it_csm != listOfCsm->end(); ++it_csm) { + csmId = it_csm->first; + listOfAmt = it_csm->second->getListOfElements(); + + for (it_amt = listOfAmt->begin(); it_amt != listOfAmt->end(); ++it_amt) { + amtId = it_amt->first; + // amtId = it_amt->second->moduleId(); + + for (int tubeId = 0; tubeId < 24; ++tubeId) { + + /* Get the offline ID, given the current detector element */ + if (!m_cablingSvc->getOfflineId( + subdetectorId, rodId, csmId, amtId, tubeId, stationName, + stationEta, stationPhi, multiLayer, layer, tube)) { + std::ostringstream ss; + ss << "Trying to initialize non-existing tube with\n"; + ss << " subdetectorId " << subdetectorId << "\n"; + ss << " rodId " << rodId << "\n"; + ss << " csmId " << csmId << "\n"; + ss << " amtId " << amtId << "\n"; + ss << " tubeId " << tubeId; + ATH_MSG_WARNING(ss.str()); + continue; + } + + bool isValid = false; + const Identifier channelIdentifier = + idhelper->channelID(stationName, stationEta, stationPhi, + multiLayer, layer, tube, true, &isValid); + // this debug msg can be removed eventually + std::ostringstream ss; + ss << "Trying to set from online IDs\n"; + ss << " subdetectorId " << subdetectorId << "\n"; + ss << " rodId " << rodId << "\n"; + ss << " csmId " << csmId << "\n"; + ss << " amtId " << amtId << "\n"; + ss << " tubeId " << tubeId << "\n"; + ss << " .. that are converted to offline IDs\n"; + ss << " stationName " << stationName << "\n"; + ss << " stationEta " << stationEta << "\n"; + ss << " stationPhi " << stationPhi << "\n"; + ss << " multiLayer " << multiLayer << "\n"; + ss << " layer " << layer << "\n"; + ss << " tube " << tube << "\n"; + ss << " .. with channelID\n"; + ss << " Identifier " << channelIdentifier << "\n"; + ss << " Identifier32 " << channelIdentifier.get_identifier32(); + ATH_MSG_DEBUG( ss.str() ); + + if (!isValid) { + ATH_MSG_FATAL( ss.str() ); + ATH_MSG_FATAL("Conversion to Identifier is NOT valid."); + return StatusCode::FAILURE; + } + + float rn = rng.Gaus(m_centralValue, m_sigma); + // Fatal if entry exists + if (m_shiftValues.find(channelIdentifier) != m_shiftValues.end()) { + ATH_MSG_FATAL("Double counting in initialization of map"); + return StatusCode::FAILURE; + } + m_shiftValues[channelIdentifier] = rn; + } + } + } + } + } + + /* initalization was successful */ + m_mapIsInitialized = true; + + /* write map to the default location */ + ATH_CHECK( dumpMapAsFile() ); + + return StatusCode::SUCCESS; +} diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationTMaxShiftSvc.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationTMaxShiftSvc.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c9fc712f3a8ee5ac965b2d5401c093ba632067cc --- /dev/null +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/MdtCalibrationTMaxShiftSvc.cxx @@ -0,0 +1,191 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +// std +#include <sys/stat.h> +#include <cstring> +#include <fstream> +#include <iostream> +#include <string> + +// other packages +#include "GaudiKernel/GaudiException.h" +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/ServiceHandle.h" +#include "Identifier/Identifier.h" +#include "MuonCablingData/MuonMDT_CablingMap.h" +#include "MuonIdHelpers/MdtIdHelper.h" +#include "StoreGate/DataHandle.h" +#include "StoreGate/StoreGateSvc.h" + +// this package +#include "MdtCalibSvc/MdtCalibrationShiftMapBase.h" +#include "MdtCalibSvc/MdtCalibrationTMaxShiftSvc.h" + +// TRandom2 for random smearing +// should later be replaced with std +#include "TRandom2.h" +// temporary +#include "TCanvas.h" +#include "TH1D.h" + +// +// private helper functions +// + +MdtCalibrationTMaxShiftSvc::MdtCalibrationTMaxShiftSvc(const std::string &name, + ISvcLocator *pSvcLocator) + : MdtCalibrationShiftMapBase(name, pSvcLocator), m_tUpper(688.1818) {} + +MdtCalibrationTMaxShiftSvc::~MdtCalibrationTMaxShiftSvc() { ; } + +StatusCode MdtCalibrationTMaxShiftSvc::initializeMap() { + if (m_mapIsInitialized) { + ATH_MSG_WARNING( + "Map already initalized. Multiple calls of initalizeMap " + "should not happen."); + return StatusCode::SUCCESS; + } + + /* First try loading the map from its default location */ + if (!m_forceMapRecreate) { + ATH_CHECK(loadMapFromFile()); + return StatusCode::SUCCESS; + } + + /* map was not found as a file */ + // TODO: What if multiple jobs started in parallel? + // (DataRace on MapFile because of parallel jobs) + + /* initialize random number generator that creates the shift values */ + TRandom2 rng(/*seed*/ 20160211); + + /* idHelper to retrieve channel Identifiers */ + const MdtIdHelper *idhelper = nullptr; + const ServiceHandle<StoreGateSvc> detStore("StoreGateSvc/DetectorStore", + "detStore"); + ATH_CHECK(detStore->retrieve(idhelper, "MDTIDHELPER")); + + // Cabling Map (part of cablingSvc) must be retrieved AFTER the first event + // i.e. this can NOT happen in initialize + /* initialize cabling from which IDs are read */ + ATH_CHECK(m_cablingSvc.retrieve()); + + /* Get ROBs */ + std::vector<uint32_t> robVector = m_cablingSvc->getAllROBId(); + + std::map<uint8_t, MdtSubdetectorMap *> *listOfSubdet; + std::map<uint8_t, MdtSubdetectorMap *>::const_iterator it_sub; + + std::map<uint8_t, MdtRODMap *> *listOfROD; + std::map<uint8_t, MdtRODMap *>::const_iterator it_rod; + + std::map<uint8_t, MdtCsmMap *> *listOfCsm; + std::map<uint8_t, MdtCsmMap *>::const_iterator it_csm; + + std::map<uint8_t, MdtAmtMap *> *listOfAmt; + std::map<uint8_t, MdtAmtMap *>::const_iterator it_amt; + + DataHandle<MuonMDT_CablingMap> cablingMap = m_cablingSvc->getCablingMap(); + listOfSubdet = cablingMap->getListOfElements(); + + int subdetectorId, rodId, csmId, amtId; + int stationName, stationEta, stationPhi, multiLayer, layer, tube; + + for (it_sub = listOfSubdet->begin(); it_sub != listOfSubdet->end(); + ++it_sub) { + subdetectorId = it_sub->first; + listOfROD = it_sub->second->getListOfElements(); + + for (it_rod = listOfROD->begin(); it_rod != listOfROD->end(); ++it_rod) { + rodId = it_rod->first; + listOfCsm = it_rod->second->getListOfElements(); + + for (it_csm = listOfCsm->begin(); it_csm != listOfCsm->end(); ++it_csm) { + csmId = it_csm->first; + listOfAmt = it_csm->second->getListOfElements(); + + for (it_amt = listOfAmt->begin(); it_amt != listOfAmt->end(); + ++it_amt) { + amtId = it_amt->first; + amtId = it_amt->second->moduleId(); + + for (int tubeId = 0; tubeId < 24; ++tubeId) { + /* Get the offline ID, given the current detector element */ + if (!m_cablingSvc->getOfflineId( + subdetectorId, rodId, csmId, amtId, tubeId, stationName, + stationEta, stationPhi, multiLayer, layer, tube)) { + std::ostringstream ss; + ss << "Trying to initialize non-existing tube with\n"; + ss << " subdetectorId " << subdetectorId << "\n"; + ss << " rodId " << rodId << "\n"; + ss << " csmId " << csmId << "\n"; + ss << " amtId " << amtId << "\n"; + ss << " tubeId " << tubeId; + ATH_MSG_WARNING(ss.str()); + continue; + } + + bool isValid = false; + const Identifier channelIdentifier = + idhelper->channelID(stationName, stationEta, stationPhi, + multiLayer, layer, tube, true, &isValid); + + // this debug msg can be removed eventually + std::ostringstream ss; + ss << "Trying to set from online IDs\n"; + ss << " subdetectorId " << subdetectorId << "\n"; + ss << " rodId " << rodId << "\n"; + ss << " csmId " << csmId << "\n"; + ss << " amtId " << amtId << "\n"; + ss << " tubelId " << tubeId << "\n"; + ss << " .. that are converted to offline IDs\n"; + ss << " stationName " << stationName << "\n"; + ss << " stationEta " << stationEta << "\n"; + ss << " stationPhi " << stationPhi << "\n"; + ss << " multiLayer " << multiLayer << "\n"; + ss << " layer " << layer << "\n"; + ss << " tube " << tube << "\n"; + ss << " .. with channelID\n"; + ss << " Identifier " << channelIdentifier << "\n"; + ss << " Identifier32 " << channelIdentifier.get_identifier32(); + ATH_MSG_DEBUG(ss.str()); + // Throw fatal msg if entry exists + if (!isValid) { + ATH_MSG_FATAL(ss.str()); + ATH_MSG_FATAL("Conversion to Identifier is NOT valid."); + return StatusCode::FAILURE; + } + if (m_shiftValues.find(channelIdentifier) != m_shiftValues.end()) { + ATH_MSG_FATAL("Double counting in initilization of map"); + return StatusCode::FAILURE; + } + float rn = rng.Gaus(m_centralValue, m_sigma); + float shift = + rn / m_tUpper; // store relative variation: Delta_t / t + m_shiftValues[channelIdentifier] = shift; + } + } + } + } + } + + /* initalization was successful */ + m_mapIsInitialized = true; + + /* write map to the default location */ + ATH_CHECK(dumpMapAsFile()); + + return StatusCode::SUCCESS; +} + +StatusCode MdtCalibrationTMaxShiftSvc::setTUpper(float tUpper) { + if (m_mapIsInitialized) { + ATH_MSG_FATAL("You cannot change m_tUpper once the map is initialized."); + return StatusCode::FAILURE; + } + + m_tUpper = tUpper; + return StatusCode::SUCCESS; +} diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/components/MdtCalibSvc_entries.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/components/MdtCalibSvc_entries.cxx index d2a114fde8d388d8ac4d00f13744d9aa95356738..17090fa3c4c0c1dcd84fb44836840ea92269c5e6 100644 --- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/components/MdtCalibSvc_entries.cxx +++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc/src/components/MdtCalibSvc_entries.cxx @@ -3,13 +3,19 @@ #include "MdtCalibSvc/MdtCalibrationSvc.h" #include "MdtCalibSvc/MdtCalibrationRegionSvc.h" #include "MdtCalibSvc/MdtCalibrationDbSvc.h" +#include "MdtCalibSvc/MdtCalibrationT0ShiftSvc.h" +#include "MdtCalibSvc/MdtCalibrationTMaxShiftSvc.h" DECLARE_SERVICE_FACTORY ( MdtCalibrationSvc ) DECLARE_SERVICE_FACTORY ( MdtCalibrationRegionSvc ) DECLARE_SERVICE_FACTORY ( MdtCalibrationDbSvc ) +DECLARE_SERVICE_FACTORY ( MdtCalibrationT0ShiftSvc ) +DECLARE_SERVICE_FACTORY ( MdtCalibrationTMaxShiftSvc ) DECLARE_FACTORY_ENTRIES(MdtCalibSvc) { DECLARE_SERVICE(MdtCalibrationSvc) DECLARE_SERVICE(MdtCalibrationRegionSvc) DECLARE_SERVICE(MdtCalibrationDbSvc) + DECLARE_SERVICE(MdtCalibrationT0ShiftSvc) + DECLARE_SERVICE(MdtCalibrationTMaxShiftSvc) }