12 #include "CLHEP/Random/RandGaussZiggurat.h"
13 #include "CoralBase/Attribute.h"
14 #include "CoralBase/AttributeListSpecification.h"
15 #include "GaudiKernel/PhysicalConstants.h"
44 #include "GeoModelHelpers/throwExcept.h"
72 ATH_MSG_DEBUG(
" initialize Random Number Service: running with t0 shift "
77 ATH_MSG_ERROR(
"Could not get random number engine from AthRNGSvc. Abort.");
78 return StatusCode::FAILURE;
83 ATH_MSG_INFO(
"************************************" << std::endl
84 <<
" Running with Calibration Deformations! " << std::endl
85 <<
" For performance studies only!" << std::endl
86 <<
" **************************************");
97 return StatusCode::SUCCESS;
104 if (
key.empty())
continue;
106 if (!readHandle.isValid()) {
107 ATH_MSG_FATAL(
"Failed to retrieve conditions object "<<readHandle.fullKey());
108 return StatusCode::FAILURE;
111 ATH_MSG_INFO(
"Size of CondAttrListCollection " << readHandle.fullKey() <<
" readCdoRt->size()= " << readHandle->size());
112 ATH_MSG_INFO(
"Range of input is " << readHandle.getRange());
116 if (!readHandle.isValid()) {
118 return StatusCode::FAILURE;
121 return StatusCode::SUCCESS;
126 if (writeHandle.isValid()) {
127 ATH_MSG_DEBUG(
"CondHandle " << writeHandle.fullKey() <<
" is already valid."
128 <<
". In theory this should not be called, but may happen"
129 <<
" if multiple concurrent events are being processed out of order.");
130 return StatusCode::SUCCESS;
137 gran = RegionGranularity::OneRt;
139 ATH_MSG_DEBUG(
"Save one set of calibration constants per multi layer");
140 gran = RegionGranularity::OnePerMultiLayer;
141 }
else ATH_MSG_DEBUG(
"Save one set of calibration constants per chamber");
142 std::unique_ptr<MuonCalib::MdtCalibDataContainer> writeCdo = std::make_unique<MuonCalib::MdtCalibDataContainer>(
m_idHelperSvc.get(), gran);
146 ATH_CHECK(writeHandle.record(std::move(writeCdo)));
147 return StatusCode::SUCCESS;
156 return StatusCode::FAILURE;
169 for (
unsigned int n = 0;
n <
rts.nRts(); ++
n) {
170 std::unique_ptr<MuonCalib::RtDataFromFile::RtRelation> rt(
rts.getRt(
n));
177 if (
times.size() < 2) {
179 return StatusCode::FAILURE;
182 if (
times.size() != radii.size() ||
times.size() != reso.size()) {
184 return StatusCode::FAILURE;
187 double t_min =
times[0];
188 double bin_size =
times[1] - t_min;
191 return StatusCode::FAILURE;
197 rtPars.push_back(t_min);
198 rtPars.push_back(bin_size);
201 rtPars.insert(rtPars.end(), radii.begin(), radii.end());
207 resoPars.insert(resoPars.end(), reso.begin(), reso.end());
219 if (!resoRel || !rtRel) {
232 RtRelationPtr MdtRt = std::make_unique<MuonCalib::MdtRtRelation>(std::move(rtRelRegion), std::move(resoRelRegion), 0.);
234 for(
auto itr = idHelper.detectorElement_begin();
235 itr!= idHelper.detectorElement_end();++itr){
246 if (idHelper.multilayer(detElId) == 2) {
247 if (writeCdo.
granularity() != RegionGranularity::OnePerMultiLayer)
continue;
248 const Identifier firstML = idHelper.multilayerID(detElId, 1);
258 if (!writeCdo.
storeData(detElId, storeMe, msgStream())) {
260 return StatusCode::FAILURE;
264 loadedRts[detElId] = MdtRt;
270 int npoints = rtRel->nPar() - 2;
272 for (
int ipt = 0; ipt <
npoints; ++ipt) {
273 double t = t_min + ipt * bin_size;
274 ATH_MSG_VERBOSE(
" " << ipt <<
" " <<
t <<
" " << rtRel->radius(
t) <<
" " << resoRel->resolution(
t));
281 return StatusCode::SUCCESS;
284 static std::atomic<bool> rtWarningPrinted =
false;
288 if (detEl) {
return std::make_optional<double>(detEl->
innerTubeRadius()); }
292 if (detEl) {
return std::make_optional<double>(detEl->
innerTubeRadius()); }
294 if (!rtWarningPrinted) {
295 ATH_MSG_WARNING(
"getInnerTubeRadius() - no Muon station known under the name "
297 rtWarningPrinted =
true;
303 std::string
data{}, delim{};
308 return StatusCode::FAILURE;
312 data = *(
static_cast<const std::string *
>((attr[
"data"]).addressOfData()));
318 return StatusCode::FAILURE;
324 unsigned int numPoints{0};
327 channel[
"appliedRT"] = rt_ts_applied;
332 if(tokensHeader.size()< 2){
333 ATH_MSG_FATAL(
"Failed to deduce extract number of points & calib Identifier from "<<
header);
334 return StatusCode::FAILURE;
336 unsigned int calibId = tokensHeader[0];
337 numPoints = tokensHeader[1];
341 return StatusCode::FAILURE;
345 ATH_MSG_WARNING(
"The translation from the calibration ID with station: "
346 <<
id.stationNameString()<<
"("<<
id.
stationName()<<
") "
347 <<
" eta:"<<
id.
eta()<<
" phi: "<<
id.
phi());
352 channel[
"ml"] = idHelper.multilayer(athenaId);
353 channel[
"layer"] = idHelper.tubeLayer(athenaId);
354 channel[
"tube"] = idHelper.tube(athenaId);
358 std::vector<double> radii{},
times{}, resos{};
359 radii.reserve(numPoints);
360 times.reserve(numPoints);
361 resos.reserve(numPoints);
365 for (
unsigned int k = 0 ;
k < dataPoints.size(); ++
k) {
366 const double value = dataPoints[
k];
369 radii.push_back(
value);
375 resos.push_back(
value);
382 if (radii.size() != numPoints ||
383 times.size() != numPoints ||
384 resos.size() != numPoints) {
385 ATH_MSG_FATAL(
"Payload "<<
payload<<
" does not lead to the expected number of points "<<numPoints<<
" vs. "<<dataPoints.size());
386 return StatusCode::FAILURE;
388 channel[
"radii"] = std::move(radii);
390 channel[
"resolutions"] = std::move(resos);
392 return StatusCode::SUCCESS;
400 if (!readHandleRt.isValid()) {
402 return StatusCode::FAILURE;
410 itr != readHandleRt->end(); ++itr) {
417 return StatusCode::FAILURE;
421 data = *(
static_cast<const std::string *
>((atr[
"data"]).addressOfData()));
425 for (
auto &
it :
yy.items()) {
427 rtCalibJson.push_back(yx);
434 itr != readHandleRt->end(); ++itr) {
441 for (
const auto&
payload : rtCalibJson) {
442 const bool rt_ts_applied =
payload[
"appliedRT"];
444 const std::string stName =
payload[
"station"];
448 std::optional<double> innerTubeRadius =
getInnerTubeRadius(idHelper.multilayerID(athenaId, 1));
449 if (!innerTubeRadius)
continue;
452 const std::vector<double> radii =
payload[
"radii"];
454 const std::vector<double> resolutions =
payload[
"resolutions"];
467 std::vector<MuonCalib::SamplePoint> tr_points{}, ts_points{};
472 for (
unsigned int k = 0;
k < radii.size(); ++
k) {
478 radius = oldradius + rshift;
479 ATH_MSG_DEBUG(
"DEFORM RT: old radius " << oldradius <<
" new radius " <<
radius <<
" shift " << rshift
494 float sigma = resolutions[
k];
498 if (tr_point.
x2() < -99) {
499 multilayer_tmax_diff = tr_point.
x1();
500 }
else if (
k == 0 || (tr_points[
k - 1].
x1() < tr_point.
x1() && tr_points[
k - 1].x2() < tr_point.
x2())) {
501 tr_points.push_back(tr_point);
502 ts_points.push_back(ts_point);
507 if (ts_points.size() < 3) {
509 return StatusCode::FAILURE;
513 float sign(rt_ts_applied ? -1.0 : 1.0);
515 for (
auto & tr_point : tr_points) {
516 int slice_number =
static_cast<int>(std::floor(tr_point.
x2() / slice_width));
517 if (slice_number < 0) slice_number = 0;
529 ATH_MSG_VERBOSE(point.x1() <<
"|" << point.x2() <<
"|" << point.error());
537 if (!reso || !rt) {
continue; }
539 if (rt->
par(1) == 0.) {
543 return StatusCode::FAILURE;
546 if (multilayer_tmax_diff > -8e8) { rt->
SetTmaxDiff(multilayer_tmax_diff); }
548 RtRelationPtr rt_rel = std::make_unique<MuonCalib::MdtRtRelation>(std::move(rt), std::move(reso), 0.);
550 if (!writeCdo.
storeData(athenaId ,rt_rel, msgStream()))
return StatusCode::FAILURE;
552 loadedRtRel[athenaId] = rt_rel;
557 if (loadedRtRel.empty()) {
558 return StatusCode::SUCCESS;
561 ATH_MSG_DEBUG(
"Initializing " << loadedRtRel.size()<<
" b-field functions");
566 if (readCondHandleDb->hasDCS()) {
567 condDbData = readCondHandleDb.cptr();
569 ATH_MSG_INFO(
"Do not retrieve the HV from DCS. Fall back to 2730 & 3080");
574 for (
const auto& [athenaId, rtRelation] : loadedRtRel) {
575 CorrectionPtr corrFuncSet = std::make_unique<MuonCalib::MdtCorFuncSet>();
578 std::vector<double> corr_params(2);
579 bool loadDefault{
false};
582 corr_params[0] = dcs.readyVolt;
584 if (corr_params[0] < std::numeric_limits<float>::epsilon()) {
588 }
else loadDefault =
true;
591 corr_params[0] = 2730.0;
593 corr_params[0] = 3080.0;
596 corr_params[1] = 0.11;
597 corrFuncSet->
setBField(std::make_unique<MuonCalib::BFieldCorFunc>(
"medium", corr_params, rtRelation->rt()));
603 if (!writeCdo.
storeData(athenaId, corrFuncSet, msgStream()))
return StatusCode::FAILURE;
606 return StatusCode::SUCCESS;
619 for (;
it != it_end; ++
it) {
630 if (!writeCdo.
storeData(*
it, tubes, msgStream()))
return StatusCode::FAILURE;
636 unsigned int nlayers = tubes->
numLayers();
637 unsigned int ntubes = tubes->
numTubes();
638 int size = nml * nlayers * ntubes;
641 <<
" size " <<
size <<
" ml " << nml <<
" l " << nlayers <<
" t " << ntubes);
642 for (
unsigned int ml = 1; ml <= nml; ++ml) {
643 for (
unsigned int l = 1;
l <= nlayers; ++
l) {
644 for (
unsigned int t = 1;
t <= ntubes; ++
t) {
649 data.inversePropSpeed = inversePropSpeed;
655 return StatusCode::SUCCESS;
665 return StatusCode::FAILURE;
669 data = *(
static_cast<const std::string *
>((attr[
"data"]).addressOfData()));
674 return StatusCode::FAILURE;
683 const std::string stName =
header.substr(2,3);
684 int eta{0},
phi{0}, nTubes{0};
687 const std::vector<std::string> headerTokens =
tokenize(
header,
",");
690 nTubes =
atoi(headerTokens[5]);
696 static std::atomic<bool> idWarningPrinted =
false;
697 if (!idWarningPrinted) {
699 <<
" is invalid. Skipping");
700 idWarningPrinted.store(
true, std::memory_order_relaxed);
702 return StatusCode::SUCCESS;
706 channel[
"appliedT0"] = t0_ts_applied;
712 std::vector<double> tzeros{}, meanAdcs{};
713 std::vector<int> statusCodes{};
716 for (
unsigned int k = 0;
k < payLoadData.size(); ++
k){
717 const double value = payLoadData[
k];
720 tzeros.push_back(
value);
723 statusCodes.push_back(
value);
726 meanAdcs.push_back(
value);
732 if (statusCodes.size() != tzeros.size() ||
733 statusCodes.size() != meanAdcs.size() ||
734 statusCodes.empty()) {
736 return StatusCode::FAILURE;
741 const int numMl = idHelper.numberOfMultilayers(chamID);
742 const Identifier secondMlID = idHelper.multilayerID(chamID, numMl);
743 const int tubesPerLay =
std::max(idHelper.tubeMax(chamID), idHelper.tubeMax(secondMlID));
744 const int numLayers =
std::max(idHelper.tubeLayerMax(chamID), idHelper.tubeLayerMax(secondMlID));
746 ATH_MSG_FATAL(
"Calibration database differs in terms of number of tubes for chamber "
748 <<
" vs. observed "<<nTubes);
749 return StatusCode::FAILURE;
752 for (
unsigned int k = 0;
k < tzeros.size(); ++
k) {
754 channelData[
"ml"] = ml;
755 channelData[
"layer"] =
layer;
756 channelData[
"tube"] =
tube;
757 channelData[
"t0"] = tzeros[
k];
758 channelData[
"meanAdc"] = meanAdcs[
k];
759 channelData[
"status"] = statusCodes[
k];
761 if (
tube > tubesPerLay){
769 calibData.push_back(std::move(channelData));
771 channel[
"calibConstants"] = std::move(calibData);
773 return StatusCode::SUCCESS;
785 itr != readHandleTube->end(); ++itr) {
792 return StatusCode::FAILURE;
796 data = *(
static_cast<const std::string *
>((atr[
"data"]).addressOfData()));
800 for (
auto &
it :
yy.items()) {
802 t0CalibJson.push_back(yx);
809 itr != readHandleTube->end(); ++itr) {
819 for (
const auto& chambChannel : t0CalibJson) {
820 const std::string stName = chambChannel[
"station"];
821 const int ieta = chambChannel[
"eta"];
822 const int iphi = chambChannel[
"phi"];
823 const bool t0_ts_applied = chambChannel[
"appliedT0"];
829 static std::atomic<bool> idWarningPrinted =
false;
830 if (!idWarningPrinted) {
831 ATH_MSG_WARNING(
"Element Identifier " << chId.get_compact() <<
" retrieved for station name " << stName
832 <<
" is not valid, skipping");
833 idWarningPrinted.store(
true, std::memory_order_relaxed);
847 if (!writeCdo.
storeData(chId, tubes, msgStream())) {
849 <<
" ID fields: "<<stName<<
","<<ieta<<
","<<iphi);
850 return StatusCode::FAILURE;
853 const nlohmann::json& tubeConstants = chambChannel[
"calibConstants"];
854 for (
const auto& tubeChannel : tubeConstants) {
855 const int ml = tubeChannel[
"ml"];
856 const int l = tubeChannel[
"layer"];
857 const int t = tubeChannel[
"tube"];
858 double tzero = tubeChannel[
"t0"];
865 double sh = CLHEP::RandGaussZiggurat::shoot(engine, 0.,
m_t0Spread);
867 ATH_MSG_VERBOSE(
"T0 spread " <<
sh <<
" t0 " << tzero <<
" id " << ml <<
" " <<
l <<
" " <<
t);
873 const double meanAdc = tubeChannel[
"meanAdc"];
878 datatube.
adcCal = meanAdc;
879 const Identifier tubeId = idHelper.channelID(chId, ml,
l,
t);
880 tubes->
setCalib(std::move(datatube), tubeId, msgStream());
886 return StatusCode::SUCCESS;
893 std::vector<Double_t>
x(sample_points.size(),0);
894 std::vector<Double_t>
y(sample_points.size(),0);
896 for (
unsigned int i = 0;
i < sample_points.size();
i++) {
897 x[
i] = sample_points[
i].x1();
898 y[
i] = sample_points[
i].x2();
900 TSpline3 sp(
"Rt Res Tmp",
x.data(),
y.data(), sample_points.size());
904 unsigned int nb_points(100);
905 std::vector<double> res_param(nb_points + 2);
906 Double_t
bin_width = (
x[sample_points.size() - 1] -
x[0]) /
static_cast<Double_t
>(nb_points);
910 for (
unsigned int k = 0;
k < nb_points;
k++) {
912 res_param[
k + 2] = sp.Eval(xx);
913 if (std::isnan(res_param[
k + 2])) {
914 TFile outf(
"kacke.root",
"RECREATE");
917 "encountered nan element");
920 return std::make_unique<MuonCalib::RtResolutionLookUp>(std::move(res_param));
925 std::vector<double> corr_params(0);
926 funcSet.
wireSag(std::make_unique<MuonCalib::WireSagCorFunc>(corr_params));