12 #include "CLHEP/Random/RandGaussZiggurat.h"
13 #include "CoralBase/Attribute.h"
14 #include "CoralBase/AttributeListSpecification.h"
15 #include "GaudiKernel/PhysicalConstants.h"
71 ATH_MSG_DEBUG(
" initialize Random Number Service: running with t0 shift "
76 ATH_MSG_ERROR(
"Could not get random number engine from AthRNGSvc. Abort.");
77 return StatusCode::FAILURE;
82 ATH_MSG_INFO(
"************************************" << std::endl
83 <<
" Running with Calibration Deformations! " << std::endl
84 <<
" For performance studies only!" << std::endl
85 <<
" **************************************");
96 return StatusCode::SUCCESS;
103 if (
key.empty())
continue;
105 if (!readHandle.isValid()) {
106 ATH_MSG_FATAL(
"Failed to retrieve conditions object "<<readHandle.fullKey());
107 return StatusCode::FAILURE;
110 ATH_MSG_INFO(
"Size of CondAttrListCollection " << readHandle.fullKey() <<
" readCdoRt->size()= " << readHandle->size());
111 ATH_MSG_INFO(
"Range of input is " << readHandle.getRange());
115 if (!readHandle.isValid()) {
117 return StatusCode::FAILURE;
120 return StatusCode::SUCCESS;
125 if (writeHandle.isValid()) {
126 ATH_MSG_DEBUG(
"CondHandle " << writeHandle.fullKey() <<
" is already valid."
127 <<
". In theory this should not be called, but may happen"
128 <<
" if multiple concurrent events are being processed out of order.");
129 return StatusCode::SUCCESS;
136 gran = RegionGranularity::OneRt;
138 ATH_MSG_DEBUG(
"Save one set of calibration constants per multi layer");
139 gran = RegionGranularity::OnePerMultiLayer;
140 }
else ATH_MSG_DEBUG(
"Save one set of calibration constants per chamber");
141 std::unique_ptr<MuonCalib::MdtCalibDataContainer> writeCdo = std::make_unique<MuonCalib::MdtCalibDataContainer>(
m_idHelperSvc.get(), gran);
145 ATH_CHECK(writeHandle.record(std::move(writeCdo)));
146 return StatusCode::SUCCESS;
155 return StatusCode::FAILURE;
168 for (
unsigned int n = 0;
n <
rts.nRts(); ++
n) {
169 std::unique_ptr<MuonCalib::RtDataFromFile::RtRelation> rt(
rts.getRt(
n));
176 if (
times.size() < 2) {
178 return StatusCode::FAILURE;
181 if (
times.size() != radii.size() ||
times.size() != reso.size()) {
183 return StatusCode::FAILURE;
186 double t_min =
times[0];
187 double bin_size =
times[1] - t_min;
190 return StatusCode::FAILURE;
196 rtPars.push_back(t_min);
197 rtPars.push_back(bin_size);
200 rtPars.insert(rtPars.end(), radii.begin(), radii.end());
206 resoPars.insert(resoPars.end(), reso.begin(), reso.end());
218 if (!resoRel || !rtRel) {
231 RtRelationPtr MdtRt = std::make_unique<MuonCalib::MdtRtRelation>(std::move(rtRelRegion), std::move(resoRelRegion), 0.);
233 for(
auto itr = idHelper.detectorElement_begin();
234 itr!= idHelper.detectorElement_end();++itr){
235 const Identifier detElId{*itr};
245 if (idHelper.multilayer(detElId) == 2) {
246 if (writeCdo.
granularity() != RegionGranularity::OnePerMultiLayer)
continue;
247 const Identifier firstML = idHelper.multilayerID(detElId, 1);
257 if (!writeCdo.
storeData(detElId, storeMe, msgStream())) {
259 return StatusCode::FAILURE;
263 loadedRts[detElId] = MdtRt;
269 int npoints = rtRel->nPar() - 2;
271 for (
int ipt = 0; ipt <
npoints; ++ipt) {
272 double t = t_min + ipt * bin_size;
273 ATH_MSG_VERBOSE(
" " << ipt <<
" " <<
t <<
" " << rtRel->radius(
t) <<
" " << resoRel->resolution(
t));
280 return StatusCode::SUCCESS;
283 static std::atomic<bool> rtWarningPrinted =
false;
287 if (detEl) {
return std::make_optional<double>(detEl->
innerTubeRadius()); }
291 if (detEl) {
return std::make_optional<double>(detEl->
innerTubeRadius()); }
293 if (!rtWarningPrinted) {
294 ATH_MSG_WARNING(
"getInnerTubeRadius() - no Muon station known under the name "
296 rtWarningPrinted =
true;
302 std::string
data{}, delim{};
307 return StatusCode::FAILURE;
311 data = *(
static_cast<const std::string *
>((attr[
"data"]).addressOfData()));
317 return StatusCode::FAILURE;
323 unsigned int numPoints{0};
326 channel[
"appliedRT"] = rt_ts_applied;
331 if(tokensHeader.size()< 2){
332 ATH_MSG_FATAL(
"Failed to deduce extract number of points & calib Identifier from "<<
header);
333 return StatusCode::FAILURE;
335 unsigned int calibId = tokensHeader[0];
336 numPoints = tokensHeader[1];
340 return StatusCode::FAILURE;
344 ATH_MSG_WARNING(
"The translation from the calibration ID with station: "
345 <<
id.stationNameString()<<
"("<<
id.
stationName()<<
") "
346 <<
" eta:"<<
id.
eta()<<
" phi: "<<
id.
phi());
351 channel[
"ml"] = idHelper.multilayer(athenaId);
352 channel[
"layer"] = idHelper.tubeLayer(athenaId);
353 channel[
"tube"] = idHelper.tube(athenaId);
357 std::vector<double> radii{},
times{}, resos{};
358 radii.reserve(numPoints);
359 times.reserve(numPoints);
360 resos.reserve(numPoints);
364 for (
unsigned int k = 0 ;
k < dataPoints.size(); ++
k) {
365 const double value = dataPoints[
k];
368 radii.push_back(
value);
374 resos.push_back(
value);
381 if (radii.size() != numPoints ||
382 times.size() != numPoints ||
383 resos.size() != numPoints) {
384 ATH_MSG_FATAL(
"Payload "<<
payload<<
" does not lead to the expected number of points "<<numPoints<<
" vs. "<<dataPoints.size());
385 return StatusCode::FAILURE;
387 channel[
"radii"] = std::move(radii);
389 channel[
"resolutions"] = std::move(resos);
391 return StatusCode::SUCCESS;
399 if (!readHandleRt.isValid()) {
401 return StatusCode::FAILURE;
409 itr != readHandleRt->end(); ++itr) {
416 return StatusCode::FAILURE;
420 data = *(
static_cast<const std::string *
>((atr[
"data"]).addressOfData()));
424 for (
auto &
it :
yy.items()) {
426 rtCalibJson.push_back(yx);
433 itr != readHandleRt->end(); ++itr) {
440 for (
const auto&
payload : rtCalibJson) {
441 const bool rt_ts_applied =
payload[
"appliedRT"];
443 const std::string stName =
payload[
"station"];
444 const Identifier athenaId = idHelper.channelID(stName,
payload[
"eta"],
payload[
"phi"],
447 std::optional<double> innerTubeRadius =
getInnerTubeRadius(idHelper.multilayerID(athenaId, 1));
448 if (!innerTubeRadius)
continue;
451 const std::vector<double> radii =
payload[
"radii"];
453 const std::vector<double> resolutions =
payload[
"resolutions"];
466 std::vector<MuonCalib::SamplePoint> tr_points{}, ts_points{};
471 for (
unsigned int k = 0;
k < radii.size(); ++
k) {
477 radius = oldradius + rshift;
478 ATH_MSG_DEBUG(
"DEFORM RT: old radius " << oldradius <<
" new radius " <<
radius <<
" shift " << rshift
493 float sigma = resolutions[
k];
497 if (tr_point.
x2() < -99) {
498 multilayer_tmax_diff = tr_point.
x1();
499 }
else if (
k == 0 || (tr_points[
k - 1].
x1() < tr_point.
x1() && tr_points[
k - 1].x2() < tr_point.
x2())) {
500 tr_points.push_back(tr_point);
501 ts_points.push_back(ts_point);
506 if (ts_points.size() < 3) {
508 return StatusCode::FAILURE;
512 float sign(rt_ts_applied ? -1.0 : 1.0);
514 for (
auto & tr_point : tr_points) {
515 int slice_number =
static_cast<int>(std::floor(tr_point.
x2() / slice_width));
516 if (slice_number < 0) slice_number = 0;
528 ATH_MSG_VERBOSE(point.x1() <<
"|" << point.x2() <<
"|" << point.error());
536 if (!reso || !rt) {
continue; }
538 if (rt->
par(1) == 0.) {
542 return StatusCode::FAILURE;
545 if (multilayer_tmax_diff > -8e8) { rt->
SetTmaxDiff(multilayer_tmax_diff); }
547 RtRelationPtr rt_rel = std::make_unique<MuonCalib::MdtRtRelation>(std::move(rt), std::move(reso), 0.);
549 if (!writeCdo.
storeData(athenaId ,rt_rel, msgStream()))
return StatusCode::FAILURE;
551 loadedRtRel[athenaId] = rt_rel;
556 if (loadedRtRel.empty()) {
557 return StatusCode::SUCCESS;
560 ATH_MSG_DEBUG(
"Initializing " << loadedRtRel.size()<<
" b-field functions");
565 if (readCondHandleDb->hasDCS()) {
566 condDbData = readCondHandleDb.cptr();
568 ATH_MSG_INFO(
"Do not retrieve the HV from DCS. Fall back to 2730 & 3080");
573 for (
const auto& [athenaId, rtRelation] : loadedRtRel) {
574 CorrectionPtr corrFuncSet = std::make_unique<MuonCalib::MdtCorFuncSet>();
577 std::vector<double> corr_params(2);
578 bool loadDefault{
false};
581 corr_params[0] = dcs.readyVolt;
583 if (corr_params[0] < std::numeric_limits<float>::epsilon()) {
587 }
else loadDefault =
true;
590 corr_params[0] = 2730.0;
592 corr_params[0] = 3080.0;
595 corr_params[1] = 0.11;
596 corrFuncSet->
setBField(std::make_unique<MuonCalib::BFieldCorFunc>(
"medium", corr_params, rtRelation->rt()));
602 if (!writeCdo.
storeData(athenaId, corrFuncSet, msgStream()))
return StatusCode::FAILURE;
605 return StatusCode::SUCCESS;
618 for (;
it != it_end; ++
it) {
629 if (!writeCdo.
storeData(*
it, tubes, msgStream()))
return StatusCode::FAILURE;
635 unsigned int nlayers = tubes->
numLayers();
636 unsigned int ntubes = tubes->
numTubes();
637 int size = nml * nlayers * ntubes;
640 <<
" size " <<
size <<
" ml " << nml <<
" l " << nlayers <<
" t " << ntubes);
641 for (
unsigned int ml = 1; ml <= nml; ++ml) {
642 for (
unsigned int l = 1;
l <= nlayers; ++
l) {
643 for (
unsigned int t = 1;
t <= ntubes; ++
t) {
645 const Identifier tubeId = id_helper.channelID(*
it, ml,
l,
t);
648 data.inversePropSpeed = inversePropSpeed;
654 return StatusCode::SUCCESS;
664 return StatusCode::FAILURE;
668 data = *(
static_cast<const std::string *
>((attr[
"data"]).addressOfData()));
673 return StatusCode::FAILURE;
682 const std::string stName =
header.substr(2,3);
683 int eta{0},
phi{0}, nTubes{0};
686 const std::vector<std::string> headerTokens =
tokenize(
header,
",");
689 nTubes =
atoi(headerTokens[5]);
693 const Identifier chamID = idHelper.elementID(stName,
eta,
phi,
isValid);
695 static std::atomic<bool> idWarningPrinted =
false;
696 if (!idWarningPrinted) {
698 <<
" is invalid. Skipping");
699 idWarningPrinted.store(
true, std::memory_order_relaxed);
701 return StatusCode::SUCCESS;
705 channel[
"appliedT0"] = t0_ts_applied;
711 std::vector<double> tzeros{}, meanAdcs{};
712 std::vector<int> statusCodes{};
715 for (
unsigned int k = 0;
k < payLoadData.size(); ++
k){
716 const double value = payLoadData[
k];
719 tzeros.push_back(
value);
722 statusCodes.push_back(
value);
725 meanAdcs.push_back(
value);
731 if (statusCodes.size() != tzeros.size() ||
732 statusCodes.size() != meanAdcs.size() ||
733 statusCodes.empty()) {
735 return StatusCode::FAILURE;
740 const int numMl = idHelper.numberOfMultilayers(chamID);
741 const Identifier secondMlID = idHelper.multilayerID(chamID, numMl);
742 const int tubesPerLay =
std::max(idHelper.tubeMax(chamID), idHelper.tubeMax(secondMlID));
743 const int numLayers =
std::max(idHelper.tubeLayerMax(chamID), idHelper.tubeLayerMax(secondMlID));
745 ATH_MSG_FATAL(
"Calibration database differs in terms of number of tubes for chamber "
747 <<
" vs. observed "<<nTubes);
748 return StatusCode::FAILURE;
751 for (
unsigned int k = 0;
k < tzeros.size(); ++
k) {
753 channelData[
"ml"] = ml;
754 channelData[
"layer"] =
layer;
755 channelData[
"tube"] =
tube;
756 channelData[
"t0"] = tzeros[
k];
757 channelData[
"meanAdc"] = meanAdcs[
k];
758 channelData[
"status"] = statusCodes[
k];
760 if (
tube > tubesPerLay){
768 calibData.push_back(std::move(channelData));
770 channel[
"calibConstants"] = std::move(calibData);
772 return StatusCode::SUCCESS;
784 itr != readHandleTube->end(); ++itr) {
791 return StatusCode::FAILURE;
795 data = *(
static_cast<const std::string *
>((atr[
"data"]).addressOfData()));
799 for (
auto &
it :
yy.items()) {
801 t0CalibJson.push_back(yx);
808 itr != readHandleTube->end(); ++itr) {
818 for (
const auto& chambChannel : t0CalibJson) {
819 const std::string stName = chambChannel[
"station"];
820 const int ieta = chambChannel[
"eta"];
821 const int iphi = chambChannel[
"phi"];
822 const bool t0_ts_applied = chambChannel[
"appliedT0"];
826 const Identifier chId = idHelper.elementID(stName, ieta, iphi,
isValid);
828 static std::atomic<bool> idWarningPrinted =
false;
829 if (!idWarningPrinted) {
830 ATH_MSG_WARNING(
"Element Identifier " << chId.get_compact() <<
" retrieved for station name " << stName
831 <<
" is not valid, skipping");
832 idWarningPrinted.store(
true, std::memory_order_relaxed);
846 if (!writeCdo.
storeData(chId, tubes, msgStream())) {
848 <<
" ID fields: "<<stName<<
","<<ieta<<
","<<iphi);
849 return StatusCode::FAILURE;
852 const nlohmann::json& tubeConstants = chambChannel[
"calibConstants"];
853 for (
const auto& tubeChannel : tubeConstants) {
854 const int ml = tubeChannel[
"ml"];
855 const int l = tubeChannel[
"layer"];
856 const int t = tubeChannel[
"tube"];
857 double tzero = tubeChannel[
"t0"];
864 double sh = CLHEP::RandGaussZiggurat::shoot(engine, 0.,
m_t0Spread);
866 ATH_MSG_VERBOSE(
"T0 spread " <<
sh <<
" t0 " << tzero <<
" id " << ml <<
" " <<
l <<
" " <<
t);
872 const double meanAdc = tubeChannel[
"meanAdc"];
877 datatube.
adcCal = meanAdc;
878 const Identifier tubeId = idHelper.channelID(chId, ml,
l,
t);
879 tubes->
setCalib(std::move(datatube), tubeId, msgStream());
885 return StatusCode::SUCCESS;
892 std::vector<Double_t>
x(sample_points.size(),0);
893 std::vector<Double_t>
y(sample_points.size(),0);
895 for (
unsigned int i = 0;
i < sample_points.size();
i++) {
896 x[
i] = sample_points[
i].x1();
897 y[
i] = sample_points[
i].x2();
899 TSpline3 sp(
"Rt Res Tmp",
x.data(),
y.data(), sample_points.size());
903 unsigned int nb_points(100);
904 std::vector<double> res_param(nb_points + 2);
905 Double_t
bin_width = (
x[sample_points.size() - 1] -
x[0]) /
static_cast<Double_t
>(nb_points);
909 for (
unsigned int k = 0;
k < nb_points;
k++) {
911 res_param[
k + 2] = sp.Eval(xx);
912 if (std::isnan(res_param[
k + 2])) {
913 TFile outf(
"kacke.root",
"RECREATE");
915 throw std::runtime_error(
"MdtCalibDbAlg::getRtResolutionInterpolation "
916 "encountered nan element");
919 return std::make_unique<MuonCalib::RtResolutionLookUp>(std::move(res_param));
924 std::vector<double> corr_params(0);
925 funcSet.
wireSag(std::make_unique<MuonCalib::WireSagCorFunc>(corr_params));