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;
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;
64 if (!map_vDrift.count(m_gasMixture)) {
65 ATH_MSG_FATAL(
"Configured Micromegas with unkown gas mixture: " << m_gasMixture);
66 return StatusCode::FAILURE;
69 m_vDrift = map_vDrift.find(m_gasMixture)->second;
70 m_transDiff = map_transDiff.find(m_gasMixture)->second;
71 m_longDiff = map_longDiff.find(m_gasMixture)->second;
72 m_interactionDensitySigma = map_interactionDensitySigma.find(m_gasMixture)->second;
73 m_interactionDensityMean = map_interactionDensityMean.find(m_gasMixture)->second;
74 m_lorentzAngleFunction = map_lorentzAngleFunctionPars.find(m_gasMixture)->second;
75 return StatusCode::SUCCESS;
81 if(!readTdoPdo.isValid()){
82 ATH_MSG_ERROR(
"Cannot find conditions data container for TDOs and PDOs!");
85 return readTdoPdo.cptr();
91 if (!ctpClusterCalibDB.isValid()) {
92 ATH_MSG_FATAL(
"Failed to retrieve the parameterized errors "<<ctpClusterCalibDB.fullKey());
95 return ctpClusterCalibDB.cptr();
101 double lorentzAngle {0.};
102 if(m_applyMmBFieldCalib){
105 if (!loadMagneticField(ctx, fieldCache))
return StatusCode::FAILURE;
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.;
114 int gasGap = m_idHelperSvc->mmIdHelper().gasGap(prepData->
identify());
116 if (changeSign) bfield = -bfield;
119 lorentzAngle = (bfield>0. ? 1. : -1.)*m_lorentzAngleFunction(std::abs(bfield)) * toRad;
130 ATH_CHECK(calibrateStrip(ctx,
id, time,
charge, (globPos.theta() / toRad) , lorentzAngle, calibStrip));
132 calibClus.push_back(std::move(calibStrip));
134 return StatusCode::SUCCESS;
140 double lorentzAngle {0.};
141 if(m_applyMmBFieldCalib){
144 if (!loadMagneticField(ctx, fieldCache))
return StatusCode::FAILURE;
146 fieldCache.
getField(globalPos.data(), magneticField.data());
149 double phi = globalPos.phi();
150 double bfield = (magneticField.x()*
std::sin(phi)-magneticField.y()*
std::cos(phi))*1000.;
153 int gasGap = m_idHelperSvc->mmIdHelper().gasGap(prepData->
identify());
155 if (changeSign) bfield = -bfield;
158 lorentzAngle = (bfield>0. ? 1. : -1.)*m_lorentzAngleFunction(std::abs(bfield)) * toRad;
169 ATH_CHECK(calibrateStrip(ctx,
id, time,
charge, (globPos.theta() / toRad) , lorentzAngle, calibStrip));
171 calibClus.push_back(std::move(calibStrip));
173 return StatusCode::SUCCESS;
183 if(!localStripPosition(
id,locPos)) {
184 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
" Failed to retrieve local strip position "<<m_idHelperSvc->toString(
id));
185 return StatusCode::FAILURE;
190 calibStrip.
time = time;
194 Identifier gasGapId = m_idHelperSvc->mmIdHelper().channelID(
id, m_idHelperSvc->mmIdHelper().multilayer(
id),
195 m_idHelperSvc->mmIdHelper().gasGap(
id), 1);
197 double vDrift = m_vDrift;
199 if(m_CalibDriftVelocityFromData){
200 vDrift = getCTPClusterCalibData(ctx)->getCTPCorrectedDriftVelocity(gasGapId, theta);
203 float max_half_drifttime = (vDrift != 0 ) ? 2.5/vDrift : 50.;
206 calibStrip.
time = time + (max_half_drifttime - m_mmT0TargetValue);
207 ATH_MSG_VERBOSE(
"Original drift time: " << time <<
" new max half drift time: " << max_half_drifttime <<
" new time: " << calibStrip.
time <<
" targett0 " << m_mmT0TargetValue );
211 double vDriftCorrected = vDrift *
std::cos(lorentzAngle);
219 calibStrip.
dx =
std::sin(lorentzAngle) * calibStrip.
time * vDrift;
221 return StatusCode::SUCCESS;
230 if(!localStripPosition(rdoId,locPos)) {
231 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
" Failed to retrieve local strip position "<<m_idHelperSvc->toString(rdoId));
232 return StatusCode::FAILURE;
245 float time{-FLT_MAX},
charge{-FLT_MAX};
252 calibStrip.
time = time - globalPos.norm() * reciprocalSpeedOfLight - mmPeakTime();
254 if(m_applyMmT0Calib){
255 calibStrip.
time = applyT0Calibration(ctx, rdoId, calibStrip.
time);
260 ATH_MSG_DEBUG(
"Calibrating RDO " << m_idHelperSvc->toString(rdoId) <<
"with pdo: " << mmRawData->
charge() <<
" tdo: "<< mmRawData->
time() <<
" relBCID "<< mmRawData->
relBcid() <<
" charge and time in counts " <<
261 mmRawData->
timeAndChargeInCounts() <<
" isData "<< m_isData <<
" to charge: " << calibStrip.
charge <<
" electrons time after corrections " << calibStrip.
time <<
" ns time before corrections "<< time <<
"ns");
272 calibStrip.
locPos = locPos;
274 return StatusCode::SUCCESS;
290 if(!localStripPosition(rdoId,locPos)) {
291 ATH_MSG_WARNING(__FILE__<<
":"<<__LINE__<<
" Failed to retrieve local strip position "<<m_idHelperSvc->toString(rdoId));
292 return StatusCode::FAILURE;
296 float time{-FLT_MAX},
charge{-FLT_MAX};
304 calibStrip.
time = time - stgcPeakTime();
306 if(m_applysTgcT0Calib){
307 calibStrip.
time = applyT0Calibration(ctx, rdoId, calibStrip.
time);
311 calibStrip.
locPos = locPos;
312 return StatusCode::SUCCESS;
317 if (!readHandle.isValid()) {
318 ATH_MSG_ERROR(
"doDigitization: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
321 readHandle.cptr()->getInitializedCache(fieldCache);
327 if (!loadMagneticField(ctx, fieldCache))
return StatusCode::FAILURE;
329 fieldCache.
getField(globalPos.data(), magneticField.data());
332 const double phi = globalPos.phi();
333 double bfield = (magneticField.x()*
std::sin(phi)-magneticField.y()*
std::cos(phi))*1000.;
336 int gasGap = m_idHelperSvc->mmIdHelper().gasGap(prepData->
identify());
338 if (changeSign) bfield = -bfield;
339 double cos2_lorentzAngle =
std::pow(
std::cos ( (bfield>0. ? 1. : -1.)*m_lorentzAngleFunction(std::abs(bfield)) * toRad), 2);
341 for (
const double dist : driftDistances){
342 double time = dist / (m_vDrift*cos2_lorentzAngle);
343 driftTimes.push_back(time);
345 return StatusCode::SUCCESS;
363 if (m_idHelperSvc->isMM (chnlId))
c /= MM_electronsPerfC;
364 else if(m_idHelperSvc->issTgc(chnlId))
c *= sTGC_pCPerfC;
391 if (m_idHelperSvc->isMM (chnlId))
charge *= MM_electronsPerfC;
392 else if(m_idHelperSvc->issTgc(chnlId))
charge /= sTGC_pCPerfC;
400 if (!tdoPdoData)
return false;
401 if (m_idHelperSvc->isMM (chnlId))
return timeToTdoMM (tdoPdoData, time, chnlId, tdo, relBCID);
402 else if(m_idHelperSvc->issTgc(chnlId))
return timeToTdoSTGC(tdoPdoData, time, chnlId, tdo, relBCID);
408 const float t = time - m_mmPeakTime - m_mmLatencyMC;
415 float tdoTime = -999.9;
418 if(
t >= lowerBound+i_relBCID*25 &&
t < (lowerBound+25)+i_relBCID*25){
419 tdoTime = i_relBCID*25 -
t;
424 if(tdoTime < lowerBound) {
434 const float t = time - m_stgcPeakTime - m_stgcLatencyMC;
441 float tdoTime = -999.9;
444 if(
t >= lowerBound+i_relBCID*25 &&
t < (lowerBound+25)+i_relBCID*25){
445 tdoTime = i_relBCID*25 -
t;
450 if(tdoTime < lowerBound) {
460 if(!readT0.isValid()){
461 ATH_MSG_ERROR(
"Cannot find conditions data container for T0s!");
464 bool isGood = readT0->getT0(
id, t0);
465 if(!isGood || t0==0){
466 ATH_MSG_DEBUG(
"failed to retrieve good t0 from database, skipping t0 calibration");
469 float targetT0 = (m_idHelperSvc->isMM(
id) ? m_mmT0TargetValue : m_stgcT0TargetValue);
470 float newTime = time + (targetT0 - t0);
471 ATH_MSG_DEBUG(
"doing T0 calibration for RDO " << m_idHelperSvc->toString(
id) <<
" time " << time <<
" t0 from database " << t0 <<
" t0 target " << targetT0 <<
" new time " << newTime);
498 float mmLatency = (m_isData? m_mmLatencyData : m_mmLatencyMC );
499 float stgcLatency = (m_isData? m_stgcLatencyData : m_stgcLatencyMC);
501 const float peakTime = m_idHelperSvc->isMM(chnlId) ? mmPeakTime() + mmLatency : stgcPeakTime() + stgcLatency;
502 time = relBCID*25. - (tdo-
calib.intercept)/
calib.slope + peakTime;
509 properties.longitudinalDiffusionSigma = m_longDiff;
510 properties.transverseDiffusionSigma = m_transDiff;
511 properties.interactionDensityMean = m_interactionDensityMean;
512 properties.interactionDensitySigma = m_interactionDensitySigma;
513 properties.lorentzAngleFunction = m_lorentzAngleFunction;
522 if(m_idHelperSvc->isMM(
id)){
526 }
else if(m_idHelperSvc->issTgc(
id)){