User will overwrite this function. Histogram booking is no longer done in C++. This function is called in execute once the filters are all passed.
103 {
104 using namespace Monitored;
105
106 const auto & numberGroup =
getGroup(
"PixelNNNumber");
107 const auto & posSummary =
getGroup(
"PixelNNPosSummary");
108 const auto & splitGroup =
getGroup(
"PixelNNSplitFrac");
109
110 SG::ReadHandle<InDet::PixelClusterContainer> pixelClusters(
m_pixelClusterKey, ctx);
111 if (!pixelClusters.isValid()) {
113 return StatusCode::SUCCESS;
114 }
115
116
117 std::vector<std::vector<const SiHit*>> siHitsByHash;
119 siHitsByHash.resize(
m_pixelID->wafer_hash_max());
120 SG::ReadHandle<SiHitCollection> siHits(
m_siHitKey, ctx);
121 if (siHits.isValid()) {
122 for (const SiHit& h : *siHits) {
123 if (!
h.isPixel()) {
continue; }
124 Identifier wid =
m_pixelID->wafer_id(
h.getBarrelEndcap(),
h.getLayerDisk(),
125 h.getPhiModule(),
h.getEtaModule());
127 if (wh < m_pixelID->wafer_hash_max()) siHitsByHash[
wh].push_back(&h);
128 }
129 } else {
130 ATH_MSG_DEBUG(
"SiHit collection not available; truth plots disabled this event");
131 }
132 }
133
134
135
136
138 std::unordered_map<Identifier::value_type, std::vector<OnTrackInfo>> onTrackClusters;
140 if (tracks.isValid()) {
141 for (const Trk::Track* track : *tracks) {
142 if (!track) { continue; }
143 for (
const auto* tsos : *
track->trackStateOnSurfaces()) {
145 const auto* rio = dynamic_cast<const Trk::RIO_OnTrack*>(tsos->measurementOnTrack());
146 if (!rio || !rio->prepRawData()) { continue; }
147 const auto* pixClus = dynamic_cast<const InDet::PixelCluster*>(rio->prepRawData());
148 if (!pixClus || !tsos->trackParameters()) { continue; }
149 onTrackClusters[pixClus->identify().get_compact()].push_back(
150 {tsos->trackParameters(), &rio->associatedSurface()});
151 }
152 }
153 }
154
155
156
157 const Trk::ClusterSplitProbabilityContainer* splitProbs = nullptr;
159 SG::ReadHandle<Trk::ClusterSplitProbabilityContainer>
h(
m_splitProbKey, ctx);
160 if (
h.isValid()) splitProbs =
h.cptr();
161 }
162
163 int nClusters = 0;
164
165
166
167 double evtMinErrX = std::numeric_limits<double>::max();
168 double evtMaxErrX = 0.;
169 double evtMinErrY = std::numeric_limits<double>::max();
170 double evtMaxErrY = 0.;
171 double evtMaxAbsDeltaX = 0.;
172 double evtMaxAbsDeltaY = 0.;
173 double evtMaxProb2 = 0.;
174 bool evtHasPos = false;
175
176 for (const auto* coll : *pixelClusters) {
177 if (!coll) { continue; }
178 for (const auto* cluster : *coll) {
179 if (!cluster) { continue; }
180 const InDetDD::SiDetectorElement* element = cluster->detectorElement();
181 if (!element) { continue; }
182 ++nClusters;
183
184 const double eta = cluster->globalPosition().eta();
185 const int nCell =
static_cast<int>(cluster->rdoList().
size());
186
187
188 std::vector<Amg::Vector2D> truths;
189 if (
m_doTruth && !siHitsByHash.empty()) {
191 }
192 const int trueN = static_cast<int>(truths.size());
193 const bool haveTruth = (trueN >= 1 && trueN <= 3);
194
195
196
197
198
199 auto it = onTrackClusters.find(cluster->identify().get_compact());
200 if (it == onTrackClusters.end() ||
it->second.empty()) {
continue; }
201
202
203
204
205
206 const Trk::Surface* surf0 =
it->second.front().surface;
208 std::vector<double> probs =
209 m_nnFactory->estimateNumberOfParticles(*cluster, *surf0, *tp0);
210 if (probs.size() < 3) { continue; }
211
212
213
214 const int predN = splitDecision(std::span<const double, 3>(probs.data(), 3));
215
216
217
218
219
220 constexpr double refPathLength = 0.250;
223 const double cosTheta = std::cos(localIntersection.theta());
224
225 if (std::abs(cosTheta) < 1e-6) { continue; }
226 localIntersection *= refPathLength / cosTheta;
227 const double trkTheta = std::atan2(localIntersection.y(), refPathLength);
228 double trkPhi = std::atan2(localIntersection.x(), refPathLength);
230 trkPhi = std::atan(std::tan(trkPhi) - tanl);
231
232 evtMaxProb2 = std::max(evtMaxProb2, probs[1]);
233 auto monPredN = Scalar<int> ("predN", predN);
234 auto monProb1 = Scalar<float>("prob1", probs[0]);
235 auto monProb2 = Scalar<float>("prob2", probs[1]);
236 auto monProb3 = Scalar<float>("prob3", probs[2]);
237 fill(numberGroup, monPredN, monProb1, monProb2, monProb3);
238
239
240
241
242 double leadPt = 0.;
243 for (
const OnTrackInfo& t :
it->second)
244 if (
t.params) leadPt = std::max(leadPt,
t.params->momentum().perp());
245 const int nnArgmax = 1 + static_cast<int>(
246 std::max_element(probs.begin(), probs.begin() + 3) - probs.begin());
247
248 auto monPt = Scalar<float>("trackPt", leadPt / Gaudi::Units::GeV);
249 auto monPhi = Scalar<float>("trkPhi", trkPhi);
250 auto monTheta = Scalar<float>("trkTheta", trkTheta);
251 auto monClusEta = Scalar<float>(
"clusEta",
eta);
252 auto monNn = Scalar<int> ("nnSplit", nnArgmax >= 2 ? 1 : 0);
253 fill(splitGroup, monPt, monPhi, monTheta, monClusEta, monNn);
254 if (splitProbs) {
256 auto monReco = Scalar<int>(
"recoSplit",
sp.isSplit() ? 1 : 0);
257 fill(splitGroup, monPt, monPhi, monTheta, monClusEta, monReco);
258 }
259 if (haveTruth) {
260 auto monTruth = Scalar<int>("truthSplit", trueN >= 2 ? 1 : 0);
261 fill(splitGroup, monPt, monPhi, monTheta, monClusEta, monTruth);
262 }
263
264 if (haveTruth) {
265 auto monTrueN = Scalar<int> ("trueN", trueN);
266 auto monPredNc = Scalar<int> ("predNconf", predN);
267 auto monIsCorrect = Scalar<int> ("isCorrect", trueN == predN ? 1 : 0);
268 auto monEta = Scalar<float>(
"eta",
eta);
269 auto monNCell = Scalar<int> ("nCell", nCell);
270 auto monProbMulti = Scalar<float>("probMulti", probs[1] + probs[2]);
271 fill(numberGroup, monTrueN, monPredNc, monIsCorrect, monEta, monNCell, monProbMulti);
272 }
273
274
275
276
277
278
279
280 const int numberOfSubclusters = predN;
281
282 const auto & posDQ =
getGroup(
"PixelNNPosDQ");
283 const auto & posGroup =
getGroup(
"PixelNNPosN" + std::to_string(numberOfSubclusters));
284 for (
const OnTrackInfo& trk :
it->second) {
285 if (!trk.params || !trk.surface) { continue; }
286 std::vector<Amg::MatrixX>
errors;
287 std::vector<Amg::Vector2D> positions =
m_nnFactory->estimatePositions(
288 *cluster, *trk.surface, *trk.params, errors, numberOfSubclusters);
289 if (static_cast<int>(positions.size()) != numberOfSubclusters ||
290 static_cast<int>(
errors.size()) != numberOfSubclusters) {
continue; }
291
292
293
294 for (
int i = 0;
i < numberOfSubclusters; ++
i) {
295 if (errors[i].
rows() < 2) {
continue; }
296 const double sigX = std::sqrt(errors[i](0, 0));
297 const double sigY = std::sqrt(errors[i](1, 1));
298 if (sigX <= 0 || sigY <= 0) { continue; }
301 auto monDeltaX = Scalar<float>("posDeltaXWide", dX);
302 auto monDeltaY = Scalar<float>("posDeltaYWide", dY);
303 auto monErrXWide = Scalar<float>("posErrXWide", sigX);
304 auto monErrYWide = Scalar<float>("posErrYWide", sigY);
305 auto monPrecX = Scalar<float>("posPrecX", 1.0 / (sigX * sigX));
306 auto monPrecY = Scalar<float>("posPrecY", 1.0 / (sigY * sigY));
307 fill(posDQ, monDeltaX, monDeltaY, monErrXWide, monErrYWide, monPrecX, monPrecY);
308 evtMinErrX = std::min(evtMinErrX, sigX);
309 evtMaxErrX = std::max(evtMaxErrX, sigX);
310 evtMinErrY = std::min(evtMinErrY, sigY);
311 evtMaxErrY = std::max(evtMaxErrY, sigY);
312 evtMaxAbsDeltaX = std::max(evtMaxAbsDeltaX, std::abs(dX));
313 evtMaxAbsDeltaY = std::max(evtMaxAbsDeltaY, std::abs(dY));
314 evtHasPos = true;
315 }
316
317 if (!haveTruth || predN != trueN) { continue; }
318
319
320
322
323 constexpr double fallbackErrX = 0.01;
324 constexpr double fallbackErrY = 0.05;
326 if (trk.params->covariance()) {
327 trkErr =
Amg::Vector2D(std::sqrt((*trk.params->covariance())(0, 0)),
328 std::sqrt((*trk.params->covariance())(1, 1)));
329 }
330 int sub = 0;
331 double best = std::numeric_limits<double>::max();
332 for (
int i = 0;
i < numberOfSubclusters; ++
i) {
333 double d = std::pow(trkPos[0] - positions[i][0], 2) / trkErr[0]
334 + std::pow(trkPos[1] - positions[i][1], 2) / trkErr[1];
335 if (d < best) {
best =
d; sub =
i; }
336 }
337 if (errors[sub].
rows() < 2) {
continue; }
338 const double errX = std::sqrt(errors[sub](0, 0));
339 const double errY = std::sqrt(errors[sub](1, 1));
340 if (errX <= 0 || errY <= 0) { continue; }
341
342
343
344
345
346 std::vector<int> perm(numberOfSubclusters);
347 std::iota(perm.begin(), perm.end(), 0);
348 std::vector<int> bestPerm = perm;
349 double bestDist = std::numeric_limits<double>::max();
350 do {
351 double dsum = 0.;
352 for (
int i = 0;
i < numberOfSubclusters; ++
i)
353 dsum += (positions[i] - truths[perm[i]]).squaredNorm();
354 if (dsum < bestDist) { bestDist = dsum; bestPerm = perm; }
355 } while (std::next_permutation(perm.begin(), perm.end()));
356 const int tr = bestPerm[sub];
357
360
361 auto monResX = Scalar<float>("resX", resX / Gaudi::Units::micrometer);
362 auto monResY = Scalar<float>("resY", resY / Gaudi::Units::micrometer);
363 auto monPullX = Scalar<float>("pullX", resX / errX);
364 auto monPullY = Scalar<float>("pullY", resY / errY);
365 auto monErrX = Scalar<float>("errX", errX / Gaudi::Units::micrometer);
366 auto monErrY = Scalar<float>("errY", errY / Gaudi::Units::micrometer);
367 auto monEta = Scalar<float>(
"eta",
eta);
368 auto monNCell = Scalar<int> ("nCell", nCell);
369 fill(posGroup, monResX, monResY, monPullX, monPullY, monErrX, monErrY, monEta, monNCell);
370
371 auto monPosN = Scalar<int>("posN", numberOfSubclusters);
372 fill(posSummary, monPosN, monPullX, monPullY, monErrX, monErrY);
373 }
374 }
375 }
376
377 auto monNClusters = Scalar<int>("nClusters", nClusters);
378 fill(numberGroup, monNClusters);
379
380 if (evtHasPos) {
381 const auto & extremes =
getGroup(
"PixelNNExtremes");
382 auto monMinErrX = Scalar<float>("evtMinErrX", evtMinErrX);
383 auto monMaxErrX = Scalar<float>("evtMaxErrX", evtMaxErrX);
384 auto monMinErrY = Scalar<float>("evtMinErrY", evtMinErrY);
385 auto monMaxErrY = Scalar<float>("evtMaxErrY", evtMaxErrY);
386 auto monMaxDX = Scalar<float>("evtMaxAbsDeltaX", evtMaxAbsDeltaX);
387 auto monMaxDY = Scalar<float>("evtMaxAbsDeltaY", evtMaxAbsDeltaY);
388 auto monMaxP2 = Scalar<float>("evtMaxProb2", evtMaxProb2);
389 fill(extremes, monMinErrX, monMaxErrX, monMinErrY, monMaxErrY, monMaxDX, monMaxDY, monMaxP2);
390 }
391 return StatusCode::SUCCESS;
392}
Scalar eta() const
pseudorapidity method
size_t size() const
Number of registered mappings.
const ToolHandle< GenericMonitoringTool > & getGroup(const std::string &name) const
Get a specific monitoring tool from the tool handle array.
virtual IdentifierHash identifyHash() const override final
identifier hash (inline)
std::vector< Amg::Vector2D > truthPositions(const InDet::PixelCluster &cluster, const InDetDD::SiDetectorElement &element, const std::vector< std::vector< const SiHit * > > &siHitsByHash) const
True particle positions (local mm) in the cluster, one per particle, from the matching Geant4 SiHits ...
ToolHandle< NnClusterizationFactory > m_nnFactory
ToolHandle< ISiLorentzAngleTool > m_lorentzTool
Gaudi::Property< bool > m_doTruth
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > m_splitProbKey
SG::ReadHandleKey< InDet::PixelClusterContainer > m_pixelClusterKey
SG::ReadHandleKey< SiHitCollection > m_siHitKey
const PixelID * m_pixelID
SG::ReadHandleKey< TrackCollection > m_trackCollectionKey
const ProbabilityInfo & splitProbability(const PrepRawData *cluster) const
const Amg::Vector3D & momentum() const
Access method for the momentum.
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
void fill(const ToolHandle< GenericMonitoringTool > &groupHandle, std::vector< std::reference_wrapper< Monitored::IMonitoredVariable > > &&variables) const
Fills a vector of variables to a group by reference.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
ParametersBase< TrackParametersDim, Charged > TrackParameters
best(iterable, priorities=[3, 2, 1, -1, 0])
1.16.1