Belle II Software prerelease-11-00-00a
SVDClusterVariables.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <tracking/dataobjects/RecoTrack.h>
10#include <mdst/dataobjects/Track.h>
11#include <analysis/dataobjects/Particle.h>
12#include <svd/variables/SVDClusterVariables.h>
13#include <svd/dataobjects/SVDTrueHit.h>
14#include <vxd/dataobjects/VxdID.h>
15#include <vxd/geometry/GeoCache.h>
16#include <vxd/geometry/SensorInfoBase.h>
17#include <framework/dataobjects/EventT0.h>
18#include <framework/datastore/StoreObjPtr.h>
19
20namespace {
21 Belle2::RecoTrack* getRecoTrack(const Belle2::Particle* particle)
22 {
23 const Belle2::Track* track = particle->getTrack();
24 if (!track)
25 return nullptr;
26 return track->getRelatedTo<Belle2::RecoTrack>();
27 }
28
29 Belle2::SVDCluster* getSVDCluster(const Belle2::RecoTrack* recoTrack, unsigned int clusterIndex)
30 {
31 if (!recoTrack) {
32 return nullptr;
33 }
34 const std::vector<Belle2::SVDCluster*> svdClusters = recoTrack->getSVDHitList();
35 if (clusterIndex >= svdClusters.size()) {
36 return nullptr;
37 }
38 return svdClusters[clusterIndex];
39 }
40
41 Belle2::SVDCluster* getSVDCluster(const Belle2::Particle* particle, unsigned int clusterIndex)
42 {
43 const Belle2::RecoTrack* recoTrack = getRecoTrack(particle);
44 return getSVDCluster(recoTrack, clusterIndex);
45 }
46}
47
48namespace Belle2::Variable {
49
50 double SVDClusterCharge(const Particle* particle, const std::vector<double>& indices)
51 {
52 if (!particle) {
53 return Const::doubleNaN;
54 }
55 if (indices.size() != 1) {
56 B2FATAL("Exactly one parameter (cluster index) is required.");
57 }
58 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
59
60 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
61 return svdCluster ? svdCluster->getCharge() : Const::doubleNaN;
62 }
63
64 double SVDClusterSNR(const Particle* particle, const std::vector<double>& indices)
65 {
66 if (!particle) {
67 return Const::doubleNaN;
68 }
69
70 if (indices.size() != 1) {
71 B2FATAL("Exactly one parameter (cluster index) is required.");
72 }
73 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
74
75 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
76 return svdCluster ? svdCluster->getSNR() : Const::doubleNaN;
77 }
78
79 int SVDClusterSize(const Particle* particle, const std::vector<double>& indices)
80 {
81 if (!particle) {
82 return -1;
83 }
84
85 if (indices.size() != 1) {
86 B2FATAL("Exactly one parameter (cluster index) is required.");
87 }
88 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
89
90 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
91 return svdCluster ? svdCluster->getSize() : -1;
92 }
93
94 double SVDClusterTime(const Particle* particle, const std::vector<double>& indices)
95 {
96 if (!particle) {
97 return Const::doubleNaN;
98 }
99 if (indices.size() != 1) {
100 B2FATAL("Exactly one parameter (cluster index) is required.");
101 }
102 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
103
104 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
105 return svdCluster ? svdCluster->getClsTime() : Const::doubleNaN;
106 }
107
108 double SVDTrackPrime(const Particle* particle, const std::vector<double>& indices)
109 {
110 if (!particle) {
111 return Const::doubleNaN;
112 }
113 if (indices.size() != 1) {
114 B2FATAL("Exactly one parameter (cluster index) is required.");
115 }
116 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
117
118 const RecoTrack* recoTrack = getRecoTrack(particle);
119 if (!recoTrack) {
120 return Const::doubleNaN;
121 }
122 const SVDCluster* svdCluster = getSVDCluster(recoTrack, clusterIndex);
123 if (!svdCluster) {
124 return Const::doubleNaN;
125 }
126 const RecoHitInformation* recoHitInformation = recoTrack->getRecoHitInformation(svdCluster);
127 if (!recoHitInformation) {
128 return Const::doubleNaN;
129 }
130 try {
131 genfit::MeasuredStateOnPlane measuredState = recoTrack->getMeasuredStateOnPlaneFromRecoHit(recoHitInformation);
132 return svdCluster->isUCluster()
133 ? measuredState.getState()[1]
134 : measuredState.getState()[2];
135 } catch (const NoTrackFitResult&) {
136 B2WARNING("No track fit result available for this hit!");
137 return Const::doubleNaN;
138 }
139 }
140
141 double SVDTrackPositionErrorUnbiased(const Particle* particle, const std::vector<double>& indices)
142 {
143 if (!particle) {
144 return Const::doubleNaN;
145 }
146 if (indices.size() != 1) {
147 B2FATAL("Exactly one parameter (cluster index) is required.");
148 }
149 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
150
151 const RecoTrack* recoTrack = getRecoTrack(particle);
152 if (!recoTrack) {
153 return Const::doubleNaN;
154 }
155 const SVDCluster* svdCluster = getSVDCluster(recoTrack, clusterIndex);
156 if (!svdCluster) {
157 return Const::doubleNaN;
158 }
159 const RecoHitInformation* recoHitInformation = recoTrack->getRecoHitInformation(svdCluster);
160 if (!recoHitInformation) {
161 return Const::doubleNaN;
162 }
163 const genfit::TrackPoint* trackPoint = recoTrack->getCreatedTrackPoint(recoHitInformation);
164 if (!trackPoint) {
165 return Const::doubleNaN;
166 }
167 const genfit::AbsFitterInfo* fitterInfo = trackPoint->getFitterInfo();
168 if (!fitterInfo) {
169 return Const::doubleNaN;
170 }
171 try {
172 genfit::MeasuredStateOnPlane unbiasedState = fitterInfo->getFittedState(false);
173 return svdCluster->isUCluster()
174 ? sqrt(unbiasedState.getCov()[3][3])
175 : sqrt(unbiasedState.getCov()[4][4]);
176 } catch (...) {
177 B2WARNING("Could not compute SVDTrackPositionErrorUnbiased.");
178 return Const::doubleNaN;
179 }
180 }
181
182 double SVDTruePosition(const Particle* particle, const std::vector<double>& indices)
183 {
184 if (!particle) {
185 return Const::doubleNaN;
186 }
187 if (indices.size() != 1) {
188 B2FATAL("Exactly one parameter (cluster index) is required.");
189 }
190 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
191
192 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
193 if (!svdCluster) {
194 return Const::doubleNaN;
195 }
196
197 const auto trueHit = svdCluster->getRelatedTo<Belle2::SVDTrueHit>();
198
199 if (!trueHit) {
200 return Const::doubleNaN;
201 }
202 return svdCluster->isUCluster()
203 ? trueHit->getU()
204 : trueHit->getV();
205 }
206
207 double SVDResidual(const Particle* particle, const std::vector<double>& indices)
208 {
209 if (!particle) {
210 return Const::doubleNaN;
211 }
212 if (indices.size() != 1) {
213 B2FATAL("Exactly one parameter (cluster index) is required.");
214 }
215 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
216
217 const RecoTrack* recoTrack = getRecoTrack(particle);
218 if (!recoTrack) {
219 return Const::doubleNaN;
220 }
221 const SVDCluster* svdCluster = getSVDCluster(recoTrack, clusterIndex);
222 if (!svdCluster) {
223 return Const::doubleNaN;
224 }
225 const RecoHitInformation* recoHitInformation = recoTrack->getRecoHitInformation(svdCluster);
226 if (!recoHitInformation) {
227 return Const::doubleNaN;
228 }
229 const genfit::TrackPoint* trackPoint = recoTrack->getCreatedTrackPoint(recoHitInformation);
230 if (!trackPoint) {
231 return Const::doubleNaN;
232 }
233 const genfit::AbsFitterInfo* fitterInfo = trackPoint->getFitterInfo();
234 if (!fitterInfo) {
235 return Const::doubleNaN;
236 }
237 try {
238 const TVectorD residualMeasurement = fitterInfo->getResidual(0, false).getState();
239 return residualMeasurement.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
240 } catch (...) {
241 B2WARNING("Could not get track residual.");
242 return Const::doubleNaN;
243 }
244 }
245
246 int SVDLayer(const Particle* particle, const std::vector<double>& indices)
247 {
248 if (!particle) {
249 return -1;
250 }
251 if (indices.size() != 1) {
252 B2FATAL("Exactly one parameter (cluster index) is required.");
253 }
254 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
255
256 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
257 if (!svdCluster) {
258 return -1;
259 }
260 const VxdID vxdId = svdCluster->getSensorID();
261 return vxdId ? vxdId.getLayerNumber() : -1;
262 }
263
264 int SVDLadder(const Particle* particle, const std::vector<double>& indices)
265 {
266 if (!particle) {
267 return -1;
268 }
269 if (indices.size() != 1) {
270 B2FATAL("Exactly one parameter (cluster index) is required.");
271 }
272 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
273
274 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
275 if (!svdCluster) {
276 return -1;
277 }
278 const VxdID vxdId = svdCluster->getSensorID();
279 return vxdId ? vxdId.getLadderNumber() : -1;
280 }
281
282 int SVDSensor(const Particle* particle, const std::vector<double>& indices)
283 {
284 if (!particle) {
285 return -1;
286 }
287 if (indices.size() != 1) {
288 B2FATAL("Exactly one parameter (cluster index) is required.");
289 }
290 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
291
292 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
293 if (!svdCluster) {
294 return -1;
295 }
296 const VxdID vxdId = svdCluster->getSensorID();
297 return vxdId ? vxdId.getSensorNumber() : -1;
298 }
299
300 bool SVDSide(const Particle* particle, const std::vector<double>& indices)
301 {
302 if (!particle) {
303 return false;
304 }
305 if (indices.size() != 1) {
306 B2FATAL("Exactly one parameter (cluster index) is required.");
307 }
308 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
309
310 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
311 return svdCluster ? svdCluster->isUCluster() : false;
312 }
313
314 double SVDClusterTimeMinusEventT0(const Particle* particle, const std::vector<double>& indices)
315 {
316 if (!particle) {
317 return Const::doubleNaN;
318 }
319 if (indices.size() != 1) {
320 B2FATAL("Exactly one parameter (cluster index) is required.");
321 }
322 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
323
324 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
325 if (!svdCluster) {
326 return Const::doubleNaN;
327 }
328
329 StoreObjPtr<EventT0> evtT0;
330 if (!evtT0.isValid()) {
331 B2WARNING("StoreObjPtr<EventT0> does not exist, are you running over cDST data?");
332 return Const::doubleNaN;
333 }
334 if (!evtT0->hasEventT0()) {
335 return Const::doubleNaN;
336 }
337
338 return svdCluster->getClsTime() - evtT0->getEventT0();
339 }
340
341 double SVDClusterChargeNormTrkLength(const Particle* particle, const std::vector<double>& indices)
342 {
343 if (!particle) {
344 return Const::doubleNaN;
345 }
346 if (indices.size() != 1) {
347 B2FATAL("Exactly one parameter (cluster index) is required.");
348 }
349 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
350
351 const RecoTrack* recoTrack = getRecoTrack(particle);
352 if (!recoTrack) {
353 return Const::doubleNaN;
354 }
355 const SVDCluster* svdCluster = getSVDCluster(recoTrack, clusterIndex);
356 if (!svdCluster) {
357 return Const::doubleNaN;
358 }
359 const RecoHitInformation* recoHitInformation = recoTrack->getRecoHitInformation(svdCluster);
360 if (!recoHitInformation) {
361 return Const::doubleNaN;
362 }
363 if (!recoTrack->wasFitSuccessful()) {
364 return Const::doubleNaN;
365 }
366 if (!recoHitInformation->useInFit()) {
367 return Const::doubleNaN;
368 }
369
370 try {
371 // Get measured state on plane to extract track angles
372 genfit::MeasuredStateOnPlane measuredState = recoTrack->getMeasuredStateOnPlaneFromRecoHit(recoHitInformation);
373 const TVectorD& state = measuredState.getState();
374
375 // Extract track prime (du/dw) and track prime OS (dv/dw)
376 // State vector elements: [0]=q/p, [1]=du/dw, [2]=dv/dw, [3]=u, [4]=v
377 const double trackPrime = state[1];
378 const double trackPrimeOS = state[2];
379
380 const double normTrkLength = sqrt(1.0 + trackPrime * trackPrime + trackPrimeOS * trackPrimeOS);
381
382 // Get cluster charge and normalize by track length
383 const double clusterCharge = svdCluster->getCharge();
384
385 return clusterCharge / normTrkLength;
386 } catch (const NoTrackFitResult&) {
387 B2WARNING("No track fit result available for this hit!");
388 return Const::doubleNaN;
389 } catch (...) {
390 B2WARNING("Could not compute SVDClusterChargeNormTrkLength.");
391 return Const::doubleNaN;
392 }
393 }
394
395 double SVDClusterTimeMinusCDCEventT0(const Particle* particle, const std::vector<double>& indices)
396 {
397 if (!particle) {
398 return Const::doubleNaN;
399 }
400 if (indices.size() != 1) {
401 B2FATAL("Exactly one parameter (cluster index) is required.");
402 }
403 const auto clusterIndex = static_cast<unsigned int>(indices[0]);
404
405 SVDCluster* svdCluster = getSVDCluster(particle, clusterIndex);
406 if (!svdCluster) {
407 return Const::doubleNaN;
408 }
409
410 StoreObjPtr<EventT0> evtT0;
411 if (!evtT0.isValid()) {
412 B2WARNING("StoreObjPtr<EventT0> does not exist, are you running over cDST data?");
413 return Const::doubleNaN;
414 }
415 auto cdcT0 = evtT0->getBestCDCTemporaryEventT0();
416 if (!cdcT0) {
417 return Const::doubleNaN;
418 }
419
420 return svdCluster->getClsTime() - cdcT0->eventT0;
421 }
422
423 double CDCEventT0(const Particle* /*particle*/, const std::vector<double>& /*indices*/)
424 {
425 StoreObjPtr<EventT0> evtT0;
426 if (!evtT0.isValid()) {
427 B2WARNING("StoreObjPtr<EventT0> does not exist, are you running over cDST data?");
428 return Const::doubleNaN;
429 }
430 auto cdcT0 = evtT0->getBestCDCTemporaryEventT0();
431 if (!cdcT0) {
432 return Const::doubleNaN;
433 }
434 return cdcT0->eventT0;
435 }
436
437 double SVDEventT0(const Particle* /*particle*/, const std::vector<double>& /*indices*/)
438 {
439 StoreObjPtr<EventT0> evtT0;
440 if (!evtT0.isValid()) {
441 B2WARNING("StoreObjPtr<EventT0> does not exist, are you running over cDST data?");
442 return Const::doubleNaN;
443 }
444 auto svdT0 = evtT0->getBestSVDTemporaryEventT0();
445 if (!svdT0) {
446 return Const::doubleNaN;
447 }
448 return svdT0->eventT0;
449 }
450
451 VARIABLE_GROUP("SVD Validation");
452
453 REGISTER_VARIABLE("SVDClusterCharge(i)", SVDClusterCharge,
454 "Returns the charge of the i-th SVD cluster related to the Particle.");
455 REGISTER_VARIABLE("SVDClusterSNR(i)", SVDClusterSNR,
456 "Returns the SNR of the i-th SVD cluster related to the Particle.");
457 REGISTER_VARIABLE("SVDClusterSize(i)", SVDClusterSize,
458 "Returns the size of the i-th SVD cluster related to the Particle.");
459 REGISTER_VARIABLE("SVDClusterTime(i)", SVDClusterTime,
460 "Returns the time of the i-th SVD cluster related to the Particle.");
461 REGISTER_VARIABLE("SVDTrackPrime(i)", SVDTrackPrime,
462 "Returns the tan of the incident angle projected on U/V of the i-th SVD cluster related to the Particle.");
463 REGISTER_VARIABLE("SVDResidual(i)", SVDResidual,
464 "Returns the track residual of the i-th SVD cluster related to the Particle.");
465 REGISTER_VARIABLE("SVDTrackPositionErrorUnbiased(i)", SVDTrackPositionErrorUnbiased,
466 "Returns the unbiased track position error of the i-th SVD cluster related to the Particle.");
467 REGISTER_VARIABLE("SVDTruePosition(i)", SVDTruePosition,
468 "Returns the true position of the i-th SVD cluster related to the Particle.");
469 REGISTER_VARIABLE("SVDClusterChargeNormTrkLength(i)", SVDClusterChargeNormTrkLength,
470 "Returns the cluster charge normalized to track length in SVD for the i-th SVD cluster related to the Particle.");
471 REGISTER_VARIABLE("SVDClusterTimeMinusEventT0(i)", SVDClusterTimeMinusEventT0,
472 "Returns the cluster time minus the event T0 for the i-th SVD cluster related to the Particle.");
473 REGISTER_VARIABLE("SVDLayer(i)", SVDLayer,
474 "Returns the layer number of the i-th SVD cluster related to the Particle. If no SVD cluster is found, returns -1.");
475 REGISTER_VARIABLE("SVDLadder(i)", SVDLadder,
476 "Returns the ladder number of the i-th SVD cluster related to the Particle. If no SVD cluster is found, returns -1.");
477 REGISTER_VARIABLE("SVDSensor(i)", SVDSensor,
478 "Returns the sensor number of the i-th SVD cluster related to the Particle. If no SVD cluster is found, returns -1.");
479 REGISTER_VARIABLE("SVDSide(i)", SVDSide,
480 "Returns true if the i-th SVD cluster related to the Particle is a U cluster.");
481 REGISTER_VARIABLE("SVDClusterTimeMinusCDCEventT0(i)", SVDClusterTimeMinusCDCEventT0,
482 "Returns the cluster time minus the best CDC-based event T0 for the i-th SVD cluster related to the Particle.");
483 REGISTER_VARIABLE("CDCEventT0", CDCEventT0,
484 "Returns the best CDC-based event T0.");
485 REGISTER_VARIABLE("SVDEventT0", SVDEventT0,
486 "Returns the best SVD-based event T0.");
487}
Belle2::Const::doubleNaN
static const double doubleNaN
quiet_NaN
Definition Const.h:703
Belle2::Particle
Class to store reconstructed particles.
Definition Particle.h:76
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition RecoTrack.h:79
Belle2::RecoTrack::getSVDHitList
std::vector< Belle2::RecoTrack::UsedSVDHit * > getSVDHitList() const
Return an unsorted list of svd hits.
Definition RecoTrack.h:452
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition SVDCluster.h:29
Belle2::Track
Class that bundles various TrackFitResults.
Definition Track.h:25
Belle2::Unit::convertValueToUnit
static double convertValueToUnit(double value, const std::string &unitString)
Converts a floating point value from the standard framework unit to the given unit.
Definition UnitConst.cc:138
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28