6 #include "GaudiKernel/SystemOfUnits.h"
7 #include "GaudiKernel/PhysicalConstants.h"
12 constexpr
double toRad =
M_PI/180;
13 constexpr
double pitchErr = 0.425 * 0.425 / 12;
17 static const std::map<std::string, float> map_transDiff {{
"ArCo2_937", 0.036},
18 {
"ArCo2_8020", 0.019}, {
"ArCo2iC4H10_9352", 0.035}};
19 static const std::map<std::string, float> map_longDiff {{
"ArCo2_937", 0.019},
20 {
"ArCo2_8020", 0.022 }, {
"ArCo2iC4H10_9352", 0.0195}};
21 static const std::map<std::string, float> map_vDrift {{
"ArCo2_937", 0.047},
22 {
"ArCo2_8020", 0.040}, {
"ArCo2iC4H10_9352", 0.045}};
26 constexpr
double const MM_electronsPerfC = 6241.;
27 constexpr
double const sTGC_pCPerfC = 1000.;
33 static const std::map<std::string, angleFunction> map_lorentzAngleFunctionPars {
34 {
"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;}},
35 {
"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;}},
36 {
"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;}}};
39 static const std::map<std::string, float> map_interactionDensitySigma {{
"ArCo2_937", 4.04 / 5.},
40 {
"ArCo2_8020", 4.04 / 5.}, {
"ArCo2iC4H10_9352", 4.04 / 5.}};
41 static const std::map<std::string, float> map_interactionDensityMean {{
"ArCo2_937", 16.15 / 5.},
42 {
"ArCo2_8020", 16.15 / 5.}, {
"ArCo2iC4H10_9352", 16.15 / 5.}};
48 declareInterface<INSWCalibTool>(
this);
57 ATH_CHECK(m_condT0Key.initialize(m_applyMmT0Calib || m_applysTgcT0Calib));
58 ATH_CHECK(m_fieldCondObjInputKey.initialize());
61 ATH_CHECK(m_ctpClusterCalibKey.initialize(m_CalibDriftVelocityFromData));
63 return StatusCode::SUCCESS;
67 if (!map_vDrift.count(m_gasMixture)) {
68 ATH_MSG_FATAL(
"Configured Micromegas with unkown gas mixture: " << m_gasMixture);
69 return StatusCode::FAILURE;
72 m_vDrift = map_vDrift.find(m_gasMixture)->second;
73 m_transDiff = map_transDiff.find(m_gasMixture)->second;
74 m_longDiff = map_longDiff.find(m_gasMixture)->second;
75 m_interactionDensitySigma = map_interactionDensitySigma.find(m_gasMixture)->second;
76 m_interactionDensityMean = map_interactionDensityMean.find(m_gasMixture)->second;
77 m_lorentzAngleFunction = map_lorentzAngleFunctionPars.find(m_gasMixture)->second;
78 return StatusCode::SUCCESS;
84 if(!readTdoPdo.isValid()){
85 ATH_MSG_ERROR(
"Cannot find conditions data container for TDOs and PDOs!");
88 return readTdoPdo.cptr();
94 if (!ctpClusterCalibDB.isValid()) {
95 ATH_MSG_FATAL(
"Failed to retrieve the parameterized errors "<<ctpClusterCalibDB.fullKey());
98 return ctpClusterCalibDB.cptr();
104 double lorentzAngle {0.};
105 if(m_applyMmBFieldCalib){
108 if (!loadMagneticField(ctx, fieldCache))
return StatusCode::FAILURE;
110 fieldCache.
getField(globalPos.data(), magneticField.data());
113 double phi = globalPos.phi();
114 double bfield = (magneticField.x()*
std::sin(phi)-magneticField.y()*
std::cos(phi))*1000.;
117 int gasGap = m_idHelperSvc->mmIdHelper().gasGap(prepData->
identify());
119 if (changeSign) bfield = -bfield;
122 lorentzAngle = (bfield>0. ? 1. : -1.)*m_lorentzAngleFunction(std::abs(bfield)) * toRad;
133 ATH_CHECK(calibrateStrip(ctx,
id, time,
charge, (globPos.theta() / toRad) , lorentzAngle, calibStrip));
135 calibClus.push_back(std::move(calibStrip));
137 return StatusCode::SUCCESS;
145 if(!localStripPosition(
id,locPos)) {
146 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
" Failed to retrieve local strip position "<<m_idHelperSvc->toString(
id));
147 return StatusCode::FAILURE;
152 calibStrip.
time = time;
156 Identifier gasGapId = m_idHelperSvc->mmIdHelper().channelID(
id, m_idHelperSvc->mmIdHelper().multilayer(
id),
157 m_idHelperSvc->mmIdHelper().gasGap(
id), 1);
159 double vDrift = m_vDrift;
161 if(m_CalibDriftVelocityFromData){
162 vDrift = getCTPClusterCalibData(ctx)->getCTPCorrectedDriftVelocity(gasGapId, theta);
165 float max_half_drifttime = (vDrift != 0 ) ? 2.5/vDrift : 50.;
168 calibStrip.
time = time + (max_half_drifttime - m_mmT0TargetValue);
169 ATH_MSG_VERBOSE(
"Original drift time: " << time <<
" new max half drift time: " << max_half_drifttime <<
" new time: " << calibStrip.
time <<
" targett0 " << m_mmT0TargetValue );
173 double vDriftCorrected = vDrift *
std::cos(lorentzAngle);
181 calibStrip.
dx =
std::sin(lorentzAngle) * calibStrip.
time * vDrift;
183 return StatusCode::SUCCESS;
192 if(!localStripPosition(rdoId,locPos)) {
193 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
" Failed to retrieve local strip position "<<m_idHelperSvc->toString(rdoId));
194 return StatusCode::FAILURE;
207 float time{-FLT_MAX},
charge{-FLT_MAX};
214 calibStrip.
time = time - globalPos.norm() * reciprocalSpeedOfLight - mmPeakTime();
216 if(m_applyMmT0Calib){
217 calibStrip.
time = applyT0Calibration(ctx, rdoId, calibStrip.
time);
222 ATH_MSG_DEBUG(
"Calibrating RDO " << m_idHelperSvc->toString(rdoId) <<
"with pdo: " << mmRawData->
charge() <<
" tdo: "<< mmRawData->
time() <<
" relBCID "<< mmRawData->
relBcid() <<
" charge and time in counts " <<
223 mmRawData->
timeAndChargeInCounts() <<
" isData "<< m_isData <<
" to charge: " << calibStrip.
charge <<
" electrons time after corrections " << calibStrip.
time <<
" ns time before corrections "<< time <<
"ns");
234 calibStrip.
locPos = locPos;
236 return StatusCode::SUCCESS;
252 if(!localStripPosition(rdoId,locPos)) {
253 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
" Failed to retrieve local strip position "<<m_idHelperSvc->toString(rdoId));
254 return StatusCode::FAILURE;
258 float time{-FLT_MAX},
charge{-FLT_MAX};
266 calibStrip.
time = time - stgcPeakTime();
268 if(m_applysTgcT0Calib){
269 calibStrip.
time = applyT0Calibration(ctx, rdoId, calibStrip.
time);
273 calibStrip.
locPos = locPos;
274 return StatusCode::SUCCESS;
279 if (!readHandle.isValid()) {
280 ATH_MSG_ERROR(
"doDigitization: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
283 readHandle.cptr()->getInitializedCache(fieldCache);
289 if (!loadMagneticField(ctx, fieldCache))
return StatusCode::FAILURE;
291 fieldCache.
getField(globalPos.data(), magneticField.data());
294 const double phi = globalPos.phi();
295 double bfield = (magneticField.x()*
std::sin(phi)-magneticField.y()*
std::cos(phi))*1000.;
298 int gasGap = m_idHelperSvc->mmIdHelper().gasGap(prepData->
identify());
300 if (changeSign) bfield = -bfield;
301 double cos2_lorentzAngle =
std::pow(
std::cos ( (bfield>0. ? 1. : -1.)*m_lorentzAngleFunction(std::abs(bfield)) * toRad), 2);
303 for (
const double dist : driftDistances){
304 double time = dist / (m_vDrift*cos2_lorentzAngle);
305 driftTimes.push_back(time);
307 return StatusCode::SUCCESS;
325 if (m_idHelperSvc->isMM (chnlId))
c /= MM_electronsPerfC;
326 else if(m_idHelperSvc->issTgc(chnlId))
c *= sTGC_pCPerfC;
353 if (m_idHelperSvc->isMM (chnlId))
charge *= MM_electronsPerfC;
354 else if(m_idHelperSvc->issTgc(chnlId))
charge /= sTGC_pCPerfC;
362 if (!tdoPdoData)
return false;
363 if (m_idHelperSvc->isMM (chnlId))
return timeToTdoMM (tdoPdoData, time, chnlId, tdo, relBCID);
364 else if(m_idHelperSvc->issTgc(chnlId))
return timeToTdoSTGC(tdoPdoData, time, chnlId, tdo, relBCID);
370 const float t = time - m_mmPeakTime - m_mmLatencyMC;
377 float tdoTime = -999.9;
380 if(
t >= lowerBound+i_relBCID*25 &&
t < (lowerBound+25)+i_relBCID*25){
381 tdoTime = i_relBCID*25 -
t;
386 if(tdoTime < lowerBound) {
396 const float t = time - m_stgcPeakTime - m_stgcLatencyMC;
403 float tdoTime = -999.9;
406 if(
t >= lowerBound+i_relBCID*25 &&
t < (lowerBound+25)+i_relBCID*25){
407 tdoTime = i_relBCID*25 -
t;
412 if(tdoTime < lowerBound) {
422 if(!readT0.isValid()){
423 ATH_MSG_ERROR(
"Cannot find conditions data container for T0s!");
426 bool isGood = readT0->getT0(
id, t0);
427 if(!isGood || t0==0){
428 ATH_MSG_DEBUG(
"failed to retrieve good t0 from database, skipping t0 calibration");
431 float targetT0 = (m_idHelperSvc->isMM(
id) ? m_mmT0TargetValue : m_stgcT0TargetValue);
432 float newTime = time + (targetT0 - t0);
433 ATH_MSG_DEBUG(
"doing T0 calibration for RDO " << m_idHelperSvc->toString(
id) <<
" time " << time <<
" t0 from database " << t0 <<
" t0 target " << targetT0 <<
" new time " << newTime);
460 float mmLatency = (m_isData? m_mmLatencyData : m_mmLatencyMC );
461 float stgcLatency = (m_isData? m_stgcLatencyData : m_stgcLatencyMC);
463 const float peakTime = m_idHelperSvc->isMM(chnlId) ? mmPeakTime() + mmLatency : stgcPeakTime() + stgcLatency;
464 time = relBCID*25. - (tdo-
calib.intercept)/
calib.slope + peakTime;
471 properties.longitudinalDiffusionSigma = m_longDiff;
472 properties.transverseDiffusionSigma = m_transDiff;
473 properties.interactionDensityMean = m_interactionDensityMean;
474 properties.interactionDensitySigma = m_interactionDensitySigma;
475 properties.lorentzAngleFunction = m_lorentzAngleFunction;
484 if(m_idHelperSvc->isMM(
id)){
488 }
else if(m_idHelperSvc->issTgc(
id)){