6#include "GaudiKernel/SystemOfUnits.h"
7#include "GaudiKernel/PhysicalConstants.h"
15 constexpr double toRad =
M_PI/180;
16 constexpr double pitchErr = 0.425 * 0.425 / 12;
17 constexpr double reciprocalSpeedOfLight = 1. / Gaudi::Units::c_light;
20 static const std::map<std::string, float> map_transDiff {{
"ArCo2_937", 0.036},
21 {
"ArCo2_8020", 0.019}, {
"ArCo2iC4H10_9352", 0.035}};
22 static const std::map<std::string, float> map_longDiff {{
"ArCo2_937", 0.019},
23 {
"ArCo2_8020", 0.022 }, {
"ArCo2iC4H10_9352", 0.0195}};
24 static const std::map<std::string, float> map_vDrift {{
"ArCo2_937", 0.047},
25 {
"ArCo2_8020", 0.040}, {
"ArCo2iC4H10_9352", 0.045}};
29 constexpr double const MM_electronsPerfC = 6241.;
30 constexpr double const sTGC_pCPerfC = 1000.;
36 static const std::map<std::string, angleFunction> map_lorentzAngleFunctionPars {
37 {
"ArCo2_937", [](
double x) {
return 0.f + 58.87f*
x -2.983f*
x*
x -10.62f*
x*
x*
x + 2.818f*
x*
x*
x*
x;}},
38 {
"ArCo2_8020", [](
double x) {
return 0.f + 58.87f*
x -2.983f*
x*
x -10.62f*
x*
x*
x + 2.818f*
x*
x*
x*
x;}},
39 {
"ArCo2iC4H10_9352", [](
double x) {
return 0.f + 58.87f*
x -2.983f*
x*
x -10.62f*
x*
x*
x + 2.818f*
x*
x*
x*
x;}}};
42 static const std::map<std::string, float> map_interactionDensitySigma {{
"ArCo2_937", 4.04 / 5.},
43 {
"ArCo2_8020", 4.04 / 5.}, {
"ArCo2iC4H10_9352", 4.04 / 5.}};
44 static const std::map<std::string, float> map_interactionDensityMean {{
"ArCo2_937", 16.15 / 5.},
45 {
"ArCo2_8020", 16.15 / 5.}, {
"ArCo2iC4H10_9352", 16.15 / 5.}};
60 return StatusCode::SUCCESS;
66 return StatusCode::FAILURE;
75 return StatusCode::SUCCESS;
82 ATH_MSG_ERROR(
"Cannot find conditions data container for TDOs and PDOs!");
85 return readTdoPdo.
cptr();
91 if (!ctpClusterCalibDB.
isValid()) {
95 return ctpClusterCalibDB.
cptr();
101 double lorentzAngle {0.};
107 fieldCache.
getField(globalPos.data(), magneticField.data());
110 double phi = globalPos.phi();
111 double bfield = (magneticField.x()*std::sin(
phi)-magneticField.y()*std::cos(
phi))*1000.;
115 bool changeSign = ( globalPos.z() < 0. ? (gasGap==1 || gasGap==3) : (gasGap==2 || gasGap==4) );
116 if (changeSign) bfield = -bfield;
123 for (
unsigned int i = 0; i < prepData->
stripNumbers().size(); ++i){
132 calibClus.push_back(std::move(calibStrip));
134 return StatusCode::SUCCESS;
142 std::vector<NSWCalib::CalibratedStrip>& calibClus)
const {
144 double lorentzAngle {0.};
150 fieldCache.
getField(globalPos.data(), magneticField.data());
153 double phi = globalPos.phi();
154 double bfield = (magneticField.x()*std::sin(
phi)-magneticField.y()*std::cos(
phi))*1000.;
158 bool changeSign = ( globalPos.z() < 0. ? (gasGap==1 || gasGap==3) : (gasGap==2 || gasGap==4) );
159 if (changeSign) bfield = -bfield;
166 for (
unsigned int i = 0; i < prepData.
stripNumbers().size(); ++i){
175 calibClus.push_back(std::move(calibStrip));
177 return StatusCode::SUCCESS;
189 return StatusCode::FAILURE;
194 calibStrip.
time = time;
207 float max_half_drifttime = (vDrift != 0 ) ? 2.5/vDrift : 50.;
211 ATH_MSG_VERBOSE(
"Original drift time: " << time <<
" new max half drift time: " << max_half_drifttime <<
" new time: " << calibStrip.
time <<
" targett0 " <<
m_mmT0TargetValue );
215 double vDriftCorrected = vDrift * std::cos(lorentzAngle);
223 calibStrip.
dx = std::sin(lorentzAngle) * calibStrip.
time * vDrift;
225 return StatusCode::SUCCESS;
236 return StatusCode::FAILURE;
249 float time{-FLT_MAX},
charge{-FLT_MAX};
256 calibStrip.
time = time - globalPos.norm() * reciprocalSpeedOfLight -
mmPeakTime();
264 ATH_MSG_DEBUG(
"Calibrating RDO " <<
m_idHelperSvc->toString(rdoId) <<
"with pdo: " << mmRawData->
charge() <<
" tdo: "<< mmRawData->
time() <<
" relBCID "<< mmRawData->
relBcid() <<
" charge and time in counts " <<
265 mmRawData->
timeAndChargeInCounts() <<
" isData "<<
m_isData <<
" to charge: " << calibStrip.
charge <<
" electrons time after corrections " << calibStrip.
time <<
" ns time before corrections "<< time <<
"ns");
276 calibStrip.
locPos = locPos;
278 return StatusCode::SUCCESS;
296 return StatusCode::FAILURE;
300 float time{-FLT_MAX},
charge{-FLT_MAX};
315 calibStrip.
locPos = locPos;
316 return StatusCode::SUCCESS;
325 readHandle.
cptr()->getInitializedCache(fieldCache);
333 fieldCache.
getField(globalPos.data(), magneticField.data());
336 const double phi = globalPos.phi();
337 double bfield = (magneticField.x()*std::sin(
phi)-magneticField.y()*std::cos(
phi))*1000.;
341 bool changeSign = ( globalPos.z() < 0. ? (gasGap==1 || gasGap==3) : (gasGap==2 || gasGap==4) );
342 if (changeSign) bfield = -bfield;
343 double cos2_lorentzAngle = std::pow(std::cos ( (bfield>0. ? 1. : -1.)*
m_lorentzAngleFunction(std::abs(bfield)) * toRad), 2);
345 for (
const double dist : driftDistances){
346 double time = dist / (
m_vDrift*cos2_lorentzAngle);
347 driftTimes.push_back(time);
349 return StatusCode::SUCCESS;
373 pdo = c * calib.slope + calib.intercept;
394 charge = (pdo-calib.intercept)/calib.slope;
404 if (!tdoPdoData)
return false;
419 float tdoTime = -999.9;
422 if(t >= lowerBound+i_relBCID*25 && t < (lowerBound+25)+i_relBCID*25){
423 tdoTime = i_relBCID*25 - t;
428 if(tdoTime < lowerBound) {
432 tdo = tdoTime*calib.slope + calib.intercept;
445 float tdoTime = -999.9;
448 if(t >= lowerBound+i_relBCID*25 && t < (lowerBound+25)+i_relBCID*25){
449 tdoTime = i_relBCID*25 - t;
454 if(tdoTime < lowerBound) {
458 tdo = tdoTime*calib.slope + calib.intercept;
465 ATH_MSG_ERROR(
"Cannot find conditions data container for T0s!");
468 bool isGood = readT0->getT0(
id,
t0);
469 if(!isGood ||
t0==0){
470 ATH_MSG_DEBUG(
"failed to retrieve good t0 from database, skipping t0 calibration");
474 float newTime = time + (targetT0 -
t0);
475 ATH_MSG_DEBUG(
"doing T0 calibration for RDO " <<
m_idHelperSvc->toString(
id) <<
" time " << time <<
" t0 from database " <<
t0 <<
" t0 target " << targetT0 <<
" new time " << newTime);
506 time = relBCID*25. - (tdo-calib.intercept)/calib.slope + peakTime;
512 properties.driftVelocity =
m_vDrift;
513 properties.longitudinalDiffusionSigma =
m_longDiff;
Scalar phi() const
phi method
Scalar theta() const
theta method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
double charge(const T &p)
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
void getField(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field value at given position xyz[3] is in mm, bxyz[3] is in kT if deriv[9] is given,...
const Amg::Transform3D & localToGlobalTrans(const ActsTrk::GeometryContext &ctx) const
Returns the local to global transformation into the ATLAS coordinate system.
An MMReadoutElement corresponds to a single STGC module; therefore typicaly a barrel muon station con...
virtual bool stripPosition(const Identifier &id, Amg::Vector2D &pos) const override final
strip position – local or global If the strip number is outside the range of valid strips,...
const MuonChannelDesign * getDesign(const Identifier &id) const
returns the MuonChannelDesign class for the given identifier
bool stripGlobalPosition(const Identifier &id, Amg::Vector3D &gpos) const
The MuonDetectorManager stores the transient representation of the Muon Spectrometer geometry and pro...
const MMReadoutElement * getMMReadoutElement(const Identifier &id) const
access via extended identifier (requires unpacking)
const sTgcReadoutElement * getsTgcReadoutElement(const Identifier &id) const
access via extended identifier (requires unpacking)
An sTgcReadoutElement corresponds to a single STGC module; therefore typicaly a barrel muon station c...
virtual bool stripPosition(const Identifier &id, Amg::Vector2D &pos) const override final
strip position - should be renamed to channel position If the strip number is outside the range of va...
bool stripGlobalPosition(const Identifier &id, Amg::Vector3D &gpos) const
Class to represent MM measurements.
const std::vector< uint16_t > & stripNumbers() const
returns the list of strip numbers
const std::vector< short int > & stripTimes() const
returns the list of times
virtual const Amg::Vector3D & globalPosition() const override final
Returns the global position.
const std::vector< int > & stripCharges() const
returns the list of charges
Temporary class to hold the MM RDO.
static constexpr int s_BCWindow
static constexpr double s_lowerTimeBound
const Identifier & identify() const
bool timeAndChargeInCounts() const
bool timeAndChargeInCounts() const
static constexpr int s_BCWindow
static constexpr double s_lowerTimeBound
unsigned int charge() const
const Identifier identify() const
const CalibConstants * getCalibForChannel(const CalibDataType type, const Identifier &channelId) const
Retrieves the calibration constant for a particular readout channel.
const DataObjID & fullKey() const
const_pointer_type cptr()
Identifier identify() const
return the identifier
const std::vector< Identifier > & rdoList() const
return the List of rdo identifiers (pointers)
const std::vector< uint16_t > & stripNumbers() const
returns the list of strip numbers
const Identifier & identify() const
: Returns the Athena identifier of the micro mega cluster It's constructed from the measurementHash &...
IdentifierHash layerHash() const
Returns the hash of the associated layer (Needed for surface retrieval)
const std::vector< int > & stripCharges() const
returns the list of charges
const std::vector< int16_t > & stripTimes() const
returns the list of times
const MuonGMR4::MmReadoutElement * readoutElement() const
Retrieve the associated MmReadoutElement.
ConstVectorMap< N > localPosition() const
Returns the local position of the measurement.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
double channelWidth() const
calculate local channel width
std::function< double(double)> angleFunction