Compute alignment deviations, given a track as input.
43 {
44 ATH_MSG_DEBUG(
"AlignmentErrorTool::makeAlignmentDeviations()");
45
46 SG::ReadCondHandle<MuonAlignmentErrorData> readHandle{
m_readKey};
47 const MuonAlignmentErrorData* readCdo{*readHandle};
48 if (!readCdo) {
50 return;
51 }
53 std::vector<deviationSummary_t> devSumVec;
54 devSumVec.reserve(deviationVec.size());
56 for (const auto & i : deviationVec) {
57 aDevSumm.translation =
i.translation;
58 aDevSumm.rotation =
i.rotation;
59 aDevSumm.stationName =
i.stationName;
60 aDevSumm.multilayer =
i.multilayer;
61 devSumVec.emplace_back(std::move(aDevSumm));
62 }
63
65 const auto& map_struct = MuonAlignmentErrorRuleCacheVec[0];
66
68 const tsosc_t* tsosc =
track.trackStateOnSurfaces();
69
70
71 unsigned int nPrecisionHits = 0;
72 for (const auto *tsos : *tsosc) {
74 continue;
75 }
76 const Trk::MeasurementBase* meas = tsos->measurementOnTrack();
77 const auto* rot = dynamic_cast<const Trk::RIO_OnTrack*>(meas);
78
79 if (!rot) {
80 const auto* crot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(meas);
81 if (crot) {
82 unsigned int index = crot->indexOfMaxAssignProb();
83 rot = &(crot->rioOnTrack(index));
84 }
85 }
86 if (!rot) continue;
87
90
92 continue;
93 }
94
95
100 continue;
101 }
102
103
104
107 continue;
108 }
109
110 MuonCalib::MuonFixedLongId calibId =
m_idTool->idToFixedLongId(channelId);
111 if (!calibId.
isValid())
continue;
112
113 ATH_MSG_DEBUG(
"Hit is in station " << calibId <<
" (MuonFixedLongId)");
114 ++nPrecisionHits;
115
116
117 const Trk::PrepRawData* prd = rot->prepRawData();
119
122 Amg::Vector3D hitU = tsos->trackParameters()->momentum().unit();
123
124
127
128
130 if (hitP.cross(zATLAS).dot(hitV)<0.0) {
131 hitV *= -1.0;
132 }
133
134
135 bool is_matched = false;
136
137
138 int multilayer = 1;
141 }
else if (calibId.
is_mmg()) {
143 }
else if (calibId.
is_stg()) {
145 }
146
147 Identifier key_id{};
149 key_id =
m_idHelperSvc->mdtIdHelper().channelID(channelId, multilayer, 1, 1);
150 }
153 }
155 key_id =
m_idHelperSvc->mmIdHelper().channelID(channelId, multilayer, 1, 1);
156 }
158 key_id =
m_idHelperSvc->stgcIdHelper().channelID(channelId, multilayer, 1, 0, 1);
159 }
162 }
165 }
166
167
168 auto range = map_struct.id_rule_map.equal_range(key_id);
169 for (
auto i =
range.first; i !=
range.second; ++i){
171
172
173 devSumVec[rule_idx].hits.push_back(rot);
174
175 devSumVec[rule_idx].sumW2 += w2;
176 devSumVec[rule_idx].sumP += w2 * hitP;
177 devSumVec[rule_idx].sumU += w2 * hitU;
178 devSumVec[rule_idx].sumV += w2 * hitV;
179
180
181 is_matched = true;
182
183
184 }
185
186 if (!is_matched) {
187 ATH_MSG_WARNING(
"The hits in the station " << calibId <<
" (MuonFixedLongId)"
188 << " couldn't be matched to any deviation regexp in the list.");
189 }
190
191 }
192
193
194
195
196 for (auto& dev : devSumVec) {
197 if (dev.hits.size() == nPrecisionHits) {
198 dev.hits.clear();
199 }
200 }
201
202
203 std::vector<const Trk::RIO_OnTrack*> v1,
v2;
204 for (auto iti = devSumVec.begin(); iti != devSumVec.end(); ++iti) {
205
206 v1 = iti->hits;
207 if (v1.empty()) {
208 continue;
209 }
210
212
213 for (auto itj = iti+1; itj != devSumVec.end(); ++itj) {
214
215 if (iti->hits.size() != itj->hits.size()) {
216 continue;
217 }
218
221
222 if (v1 == v2) {
223 auto iDev = std::distance(devSumVec.begin(), iti);
224 auto jDev = std::distance(devSumVec.begin(), itj);
225 ATH_MSG_DEBUG(
"Found deviations " << iDev <<
" and " << jDev <<
" related to the same list of hits. Merging...");
226 ATH_MSG_DEBUG(
"old (translation, rotation) systematic uncertainties for "
227 << iDev << ": " << iti->translation << ", " << iti->rotation);
228 ATH_MSG_DEBUG(
"old (translation, rotation) systematic uncertainties for "
229 << jDev << ": " << itj->translation << ", " << itj->rotation);
230
231
232 double new_translation = std::hypot(iti->translation, itj->translation);
233 double new_rotation = std::hypot(iti->rotation, itj->rotation);
234
235
236 itj->hits.clear();
237
238
239 iti->translation = new_translation;
240 iti->rotation = new_rotation;
241 ATH_MSG_DEBUG(
"New combined (translation, rotation) systematic uncertainties: " << new_translation <<
", " << new_rotation);
242
243 }
244 }
245 }
246
247
248
249 deviations.clear();
252 for (const auto & iDev : devSumVec) {
253 if (iDev.hits.empty()) {
254 continue;
255 }
256
258 double translation = iDev.translation;
259
263 double sumW2 = iDev.sumW2;
264
265 sumP *= (1. / sumW2);
266 sumU *= (1. / sumW2);
267 sumV *= (1. / sumW2);
268
269 if (translation >= 0.001 * Gaudi::Units::mm) {
270 std::size_t hitshash = 0;
271 for (const auto *it : iDev.hits) {
272 boost::hash_combine(hitshash, (
it->identify()).get_compact());
273 }
274 deviations.push_back(
275 new AlignmentTranslationDeviation(sumU.cross(sumV), translation * Gaudi::Units::mm, iDev.hits));
276 deviations.back()->setHashOfHits(hitshash);
277
279 << sumU.x() << ", " << sumU.y() << ", " << sumU.z() << ") with sigma=" << translation * Gaudi::Units::mm
280 << " mm was applied to " << iDev.hits.size()
281 << " hits matching the station: " << iDev.stationName.str() << " and the multilayer "
282 << iDev.multilayer.str());
283 }
284 if (rotation >= 0.000001 * Gaudi::Units::rad) {
285 std::size_t hitshash = 0;
286 for (const auto *it : iDev.hits) {
287 boost::hash_combine(hitshash, (
it->identify()).get_compact());
288 }
289 deviations.push_back(new AlignmentRotationDeviation(sumP, sumV, rotation * Gaudi::Units::rad, iDev.hits));
290 deviations.back()->setHashOfHits(hitshash);
291
292 ATH_MSG_DEBUG(
"A rotation around the center = (" << sumP.x() <<
", " << sumP.y() <<
", " << sumP.z() <<
") and axis = ("
293 << sumV.x() << ", " << sumV.y() << ", " << sumV.z()
294 << ") with sigma=" << rotation / Gaudi::Units::mrad << " mrad was applied to "
295 << iDev.hits.size() << " hits matching the station "
296 << iDev.stationName.str() << " and the multilayer "
297 << iDev.multilayer.str());
298 }
299
300 }
301
304 ATH_MSG_DEBUG(
"Found " << deviations.size() <<
" nuisances after duplicates merging");
306}
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
const std::vector< MuonAlignmentErrorRuleCache > & getMuonAlignmentErrorRuleCache() const
const std::vector< MuonAlignmentErrorRule > & getAlignmentErrorRules() const
size_t MuonAlignmentErrorRuleIndex
int stgMultilayer() const
Stg specific:
int mdtMultilayer() const
Mdt specific:
int mmgMultilayer() const
Mmg specific:
bool isValid() const
check validity of the identifier.
virtual const TrkDetElementBase * detectorElement() const =0
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
Identifier identify() const
return the identifier
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
virtual const Surface & surface() const =0
Return surface associated with this detector element.
Eigen::Matrix< double, 3, 1 > Vector3D
DataVector< const Trk::TrackStateOnSurface > TrackStates
void stable_sort(DataModel_detail::iterator< DVL > beg, DataModel_detail::iterator< DVL > end)
Specialization of stable_sort for DataVector/List.
1.14.0