12 #include <analysis/variables/TrackVariables.h>
13 #include <analysis/VariableManager/Manager.h>
16 #include <framework/datastore/StoreObjPtr.h>
17 #include <framework/dataobjects/Helix.h>
20 #include <mdst/dataobjects/Track.h>
21 #include <mdst/dataobjects/MCParticle.h>
22 #include <mdst/dataobjects/TrackFitResult.h>
23 #include <mdst/dataobjects/EventLevelTrackingInfo.h>
24 #include <mdst/dataobjects/HitPatternCDC.h>
25 #include <mdst/dataobjects/HitPatternVXD.h>
26 #include <mdst/dataobjects/ECLCluster.h>
29 #include <framework/logging/Logger.h>
46 TrackFitResult
const* getTrackFitResultFromParticle(Particle
const* particle)
48 const Track* track = particle->getTrack();
53 const TrackFitResult* trackFit = track->getTrackFitResultWithClosestMass(Const::ChargedStable(abs(
54 particle->getPDGCode())));
58 double trackNHits(
const Particle* part,
const Const::EDetector& det)
60 auto trackFit = getTrackFitResultFromParticle(part);
62 return std::numeric_limits<double>::quiet_NaN();
65 if (det == Const::EDetector::CDC) {
66 return trackFit->getHitPatternCDC().getNHits();
67 }
else if (det == Const::EDetector::SVD) {
68 return trackFit->getHitPatternVXD().getNSVDHits();
69 }
else if (det == Const::EDetector::PXD) {
70 return trackFit->getHitPatternVXD().getNPXDHits();
72 return std::numeric_limits<double>::quiet_NaN();
76 double trackNCDCHits(
const Particle* part)
78 return trackNHits(part, Const::EDetector::CDC);
81 double trackNSVDHits(
const Particle* part)
83 return trackNHits(part, Const::EDetector::SVD);
86 double trackNPXDHits(
const Particle* part)
88 return trackNHits(part, Const::EDetector::PXD);
91 double trackNVXDHits(
const Particle* part)
93 return trackNPXDHits(part) + trackNSVDHits(part);
96 double trackNDF(
const Particle* part)
98 auto trackFit = getTrackFitResultFromParticle(part);
100 return std::numeric_limits<double>::quiet_NaN();
102 return trackFit->getNDF();
105 double trackChi2(
const Particle* part)
107 auto trackFit = getTrackFitResultFromParticle(part);
109 return std::numeric_limits<double>::quiet_NaN();
111 return trackFit->getChi2();
114 double trackFirstSVDLayer(
const Particle* part)
116 auto trackFit = getTrackFitResultFromParticle(part);
118 return std::numeric_limits<double>::quiet_NaN();
120 return trackFit->getHitPatternVXD().getFirstSVDLayer();
123 double trackFirstPXDLayer(
const Particle* part)
125 auto trackFit = getTrackFitResultFromParticle(part);
127 return std::numeric_limits<double>::quiet_NaN();
129 return trackFit->getHitPatternVXD().getFirstPXDLayer(HitPatternVXD::PXDMode::normal);
132 double trackFirstCDCLayer(
const Particle* part)
134 auto trackFit = getTrackFitResultFromParticle(part);
136 return std::numeric_limits<double>::quiet_NaN();
138 return trackFit->getHitPatternCDC().getFirstLayer();
141 double trackLastCDCLayer(
const Particle* part)
143 auto trackFit = getTrackFitResultFromParticle(part);
145 return std::numeric_limits<double>::quiet_NaN();
147 return trackFit->getHitPatternCDC().getLastLayer();
150 double trackD0(
const Particle* part)
152 auto trackFit = getTrackFitResultFromParticle(part);
154 return std::numeric_limits<double>::quiet_NaN();
156 return trackFit->getD0();
159 double trackPhi0(
const Particle* part)
161 auto trackFit = getTrackFitResultFromParticle(part);
163 return std::numeric_limits<double>::quiet_NaN();
165 return trackFit->getPhi0();
168 double trackOmega(
const Particle* part)
170 auto trackFit = getTrackFitResultFromParticle(part);
172 return std::numeric_limits<double>::quiet_NaN();
174 return trackFit->getOmega();
177 double trackZ0(
const Particle* part)
179 auto trackFit = getTrackFitResultFromParticle(part);
181 return std::numeric_limits<double>::quiet_NaN();
183 return trackFit->getZ0();
186 double trackTanLambda(
const Particle* part)
188 auto trackFit = getTrackFitResultFromParticle(part);
190 return std::numeric_limits<double>::quiet_NaN();
192 return trackFit->getTanLambda();
195 double trackD0Error(
const Particle* part)
197 auto trackFit = getTrackFitResultFromParticle(part);
199 return std::numeric_limits<double>::quiet_NaN();
202 double errorSquared = trackFit->getCovariance5()[0][0];
203 if (errorSquared > 0.0)
204 return sqrt(errorSquared);
206 return std::numeric_limits<double>::quiet_NaN();
209 double trackPhi0Error(
const Particle* part)
211 auto trackFit = getTrackFitResultFromParticle(part);
213 return std::numeric_limits<double>::quiet_NaN();
216 double errorSquared = trackFit->getCovariance5()[1][1];
217 if (errorSquared > 0.0)
218 return sqrt(errorSquared);
220 return std::numeric_limits<double>::quiet_NaN();
223 double trackOmegaError(
const Particle* part)
225 auto trackFit = getTrackFitResultFromParticle(part);
227 return std::numeric_limits<double>::quiet_NaN();
230 double errorSquared = trackFit->getCovariance5()[2][2];
231 if (errorSquared > 0.0)
232 return sqrt(errorSquared);
234 return std::numeric_limits<double>::quiet_NaN();
237 double trackZ0Error(
const Particle* part)
239 auto trackFit = getTrackFitResultFromParticle(part);
241 return std::numeric_limits<double>::quiet_NaN();
244 double errorSquared = trackFit->getCovariance5()[3][3];
245 if (errorSquared > 0.0)
246 return sqrt(errorSquared);
248 return std::numeric_limits<double>::quiet_NaN();
251 double trackTanLambdaError(
const Particle* part)
253 auto trackFit = getTrackFitResultFromParticle(part);
255 return std::numeric_limits<double>::quiet_NaN();
258 double errorSquared = trackFit->getCovariance5()[4][4];
259 if (errorSquared > 0.0)
260 return sqrt(errorSquared);
262 return std::numeric_limits<double>::quiet_NaN();
265 double trackPValue(
const Particle* part)
267 auto trackFit = getTrackFitResultFromParticle(part);
272 return trackFit->getPValue();
275 double trackFitHypothesisPDG(
const Particle* part)
277 auto trackFit = getTrackFitResultFromParticle(part);
282 return trackFit->getParticleType().getPDGCode();
285 double trackNECLClusters(
const Particle* part)
287 const Track* track = part->getTrack();
289 return std::numeric_limits<double>::quiet_NaN();
293 for (
const ECLCluster& cluster : track->getRelationsTo<ECLCluster>())
294 if (cluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons))
296 return double(count);
300 double trackHelixExtTheta(
const Particle* part,
const std::vector<double>& pars)
302 if (pars.size() != 3) {
303 B2WARNING(
"Exactly three parameters (r, zfwd, zbwd) required.");
304 return std::numeric_limits<double>::quiet_NaN();
307 const double r = pars[0];
308 const double zfwd = pars[1];
309 const double zbwd = pars[2];
312 auto trackFit = getTrackFitResultFromParticle(part);
313 if (!trackFit)
return std::numeric_limits<double>::quiet_NaN();
316 const double z0 = trackFit->getZ0();
317 const double tanlambda = trackFit->getTanLambda();
318 const Helix h = trackFit->getHelix();
321 const double arcLength = h.getArcLength2DAtCylindricalR(r);
322 const double lHelixRadius = arcLength > 0 ? arcLength : std::numeric_limits<double>::max();
325 const double lFWD = (zfwd - z0) / tanlambda > 0 ? (zfwd - z0) / tanlambda : std::numeric_limits<double>::max();
328 const double lBWD = (zbwd - z0) / tanlambda > 0 ? (zbwd - z0) / tanlambda : std::numeric_limits<double>::max();
331 const double l = std::min(std::min(lHelixRadius, lFWD), lBWD);
333 return atan2(h.getPositionAtArcLength2D(l).Perp(), h.getPositionAtArcLength2D(l).Z());
337 double trackHelixExtPhi(
const Particle* part,
const std::vector<double>& pars)
339 if (pars.size() != 3) {
340 B2WARNING(
"Exactly three parameters (r, zfwd, zbwd) required.");
341 return std::numeric_limits<double>::quiet_NaN();
344 const double r = pars[0];
345 const double zfwd = pars[1];
346 const double zbwd = pars[2];
349 auto trackFit = getTrackFitResultFromParticle(part);
350 if (!trackFit)
return std::numeric_limits<double>::quiet_NaN();
353 const double z0 = trackFit->getZ0();
354 const double tanlambda = trackFit->getTanLambda();
355 const Helix h = trackFit->getHelix();
358 const double arcLength = h.getArcLength2DAtCylindricalR(r);
359 const double lHelixRadius = arcLength > 0 ? arcLength : std::numeric_limits<double>::max();
362 const double lFWD = (zfwd - z0) / tanlambda > 0 ? (zfwd - z0) / tanlambda : std::numeric_limits<double>::max();
365 const double lBWD = (zbwd - z0) / tanlambda > 0 ? (zbwd - z0) / tanlambda : std::numeric_limits<double>::max();
368 const double l = std::min(std::min(lHelixRadius, lFWD), lBWD);
370 return atan2(h.getPositionAtArcLength2D(l).Y(), h.getPositionAtArcLength2D(l).X());
379 double nExtraCDCHits(
const Particle*)
381 StoreObjPtr<EventLevelTrackingInfo> elti;
382 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
383 return elti->getNCDCHitsNotAssigned();
388 double nExtraCDCHitsPostCleaning(
const Particle*)
390 StoreObjPtr<EventLevelTrackingInfo> elti;
391 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
392 return elti->getNCDCHitsNotAssignedPostCleaning();
396 double hasExtraCDCHitsInLayer(
const Particle*,
const std::vector<double>& layer)
398 StoreObjPtr<EventLevelTrackingInfo> elti;
399 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
400 int ilayer = int(std::lround(layer[0]));
401 return elti->hasCDCLayer(ilayer);
405 double hasExtraCDCHitsInSuperLayer(
const Particle*,
const std::vector<double>& layer)
407 StoreObjPtr<EventLevelTrackingInfo> elti;
408 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
409 int ilayer = int(std::lround(layer[0]));
410 return elti->hasCDCSLayer(ilayer);
414 double nExtraCDCSegments(
const Particle*)
416 StoreObjPtr<EventLevelTrackingInfo> elti;
417 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
418 return elti->getNCDCSegments();
422 double nExtraVXDHitsInLayer(
const Particle*,
const std::vector<double>& layer)
424 StoreObjPtr<EventLevelTrackingInfo> elti;
425 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
426 int ilayer = int(std::lround(layer[0]));
427 return elti->getNVXDClustersInLayer(ilayer);
431 double nExtraVXDHits(
const Particle*)
433 StoreObjPtr<EventLevelTrackingInfo> elti;
434 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
436 for (uint16_t ilayer = 1; ilayer < 7; ilayer++)
437 out += elti->getNVXDClustersInLayer(ilayer);
442 double svdFirstSampleTime(
const Particle*)
444 StoreObjPtr<EventLevelTrackingInfo> elti;
445 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
446 return elti->getSVDFirstSampleTime();
453 double trackFindingFailureFlag(
const Particle*)
455 StoreObjPtr<EventLevelTrackingInfo> elti;
456 if (!elti)
return std::numeric_limits<double>::quiet_NaN();
457 return elti->hasAnErrorFlag();
460 double getHelixParameterPullAtIndex(
const Particle* particle,
const int index)
463 if (!particle) {
return std::numeric_limits<double>::quiet_NaN(); }
465 const MCParticle* mcparticle = particle->getRelatedTo<MCParticle>();
466 if (!mcparticle) {
return std::numeric_limits<double>::quiet_NaN(); }
469 if (!track) {
return std::numeric_limits<double>::quiet_NaN(); }
472 particle->getPDGCode())));
473 if (!trackfit) {
return std::numeric_limits<double>::quiet_NaN(); }
476 const TMatrixDSym measCovariance = measHelix.
getCovariance();
477 const TVector3 mcProdVertex = mcparticle->getVertex();
478 const TVector3 mcMomentum = mcparticle->getMomentum();
482 const double mcParticleCharge = mcparticle->getCharge();
485 const std::vector<double> mcHelixPars = {mcHelix.getD0(), mcHelix.getPhi0(), mcHelix.getOmega(), mcHelix.getZ0(), mcHelix.getTanLambda()};
486 const std::vector<double> measHelixPars = {measHelix.getD0(), measHelix.getPhi0(), measHelix.getOmega(), measHelix.getZ0(), measHelix.getTanLambda()};
487 const std::vector<double> measErrSquare = {measCovariance[0][0], measCovariance[1][1], measCovariance[2][2], measCovariance[3][3], measCovariance[4][4]};
489 return (mcHelixPars.at(index) - measHelixPars.at(index)) / std::sqrt(measErrSquare.at(index));
492 double getHelixD0Pull(
const Particle* part)
494 return getHelixParameterPullAtIndex(part, 0);
497 double getHelixPhi0Pull(
const Particle* part)
499 return getHelixParameterPullAtIndex(part, 1);
502 double getHelixOmegaPull(
const Particle* part)
504 return getHelixParameterPullAtIndex(part, 2);
507 double getHelixZ0Pull(
const Particle* part)
509 return getHelixParameterPullAtIndex(part, 3);
511 double getHelixTanLambdaPull(
const Particle* part)
513 return getHelixParameterPullAtIndex(part, 4);
517 VARIABLE_GROUP(
"Tracking");
518 REGISTER_VARIABLE(
"d0Pull", getHelixD0Pull,
"mc-meas/err_meas for d0");
519 REGISTER_VARIABLE(
"phi0Pull", getHelixPhi0Pull,
"mc-meas/err_meas for phi0");
520 REGISTER_VARIABLE(
"omegaPull", getHelixOmegaPull,
"mc-meas/err_meas for omega");
521 REGISTER_VARIABLE(
"z0Pull", getHelixZ0Pull,
"mc-meas/err_meas for z0");
522 REGISTER_VARIABLE(
"tanLambdaPull", getHelixTanLambdaPull,
"mc-meas/err_meas for tanLambda");
524 REGISTER_VARIABLE(
"nCDCHits", trackNCDCHits,
"Number of CDC hits associated to the track");
525 REGISTER_VARIABLE(
"nSVDHits", trackNSVDHits,
"Number of SVD hits associated to the track");
526 REGISTER_VARIABLE(
"nPXDHits", trackNPXDHits,
"Number of PXD hits associated to the track");
527 REGISTER_VARIABLE(
"nVXDHits", trackNVXDHits,
"Number of PXD and SVD hits associated to the track");
528 REGISTER_VARIABLE(
"ndf", trackNDF,
529 R
"DOC(Number of degrees of freedom of the track fit. Note that it is not NHIT-5 due to outlier hit rejection.
530 For mdst versions < 5.1, returns quiet_NaN().)DOC"
532 REGISTER_VARIABLE("chi2", trackChi2,
533 R
"DOC(Chi2 of the track fit.
534 Computed based on pValue and ndf. Note that for pValue exactly equal to 0 it returns infinity().
535 For mdst versions < 5.1, returns quiet_NaN().)DOC");
536 REGISTER_VARIABLE("firstSVDLayer", trackFirstSVDLayer,
"First activated SVD layer associated to the track");
537 REGISTER_VARIABLE(
"firstPXDLayer", trackFirstPXDLayer,
"First activated PXD layer associated to the track");
538 REGISTER_VARIABLE(
"firstCDCLayer", trackFirstCDCLayer,
"First activated CDC layer associated to the track");
539 REGISTER_VARIABLE(
"lastCDCLayer", trackLastCDCLayer,
"Last CDC layer associated to the track");
540 REGISTER_VARIABLE(
"d0", trackD0,
"Signed distance to the POCA in the r-phi plane");
541 REGISTER_VARIABLE(
"phi0", trackPhi0,
"Angle of the transverse momentum in the r-phi plane");
542 REGISTER_VARIABLE(
"omega", trackOmega,
"Curvature of the track");
543 REGISTER_VARIABLE(
"z0", trackZ0,
"z coordinate of the POCA");
544 REGISTER_VARIABLE(
"tanlambda", trackTanLambda,
"Slope of the track in the r-z plane");
545 REGISTER_VARIABLE(
"d0Err", trackD0Error,
"Error of signed distance to the POCA in the r-phi plane");
546 REGISTER_VARIABLE(
"phi0Err", trackPhi0Error,
"Error of angle of the transverse momentum in the r-phi plane");
547 REGISTER_VARIABLE(
"omegaErr", trackOmegaError,
"Error of curvature of the track");
548 REGISTER_VARIABLE(
"z0Err", trackZ0Error,
"Error of z coordinate of the POCA");
549 REGISTER_VARIABLE(
"tanlambdaErr", trackTanLambdaError,
"Error of slope of the track in the r-z plane");
550 REGISTER_VARIABLE(
"pValue", trackPValue,
"chi2 probalility of the track fit");
551 REGISTER_VARIABLE(
"trackFitHypothesisPDG", trackFitHypothesisPDG,
"PDG code of the track hypothesis actually used for the fit");
552 REGISTER_VARIABLE(
"trackNECLClusters", trackNECLClusters,
553 "Number ecl clusters matched to the track. This is always 0 or 1 with newer versions of ECL reconstruction.");
554 REGISTER_VARIABLE(
"helixExtTheta", trackHelixExtTheta,
555 "Returns theta of extrapolated helix parameters (parameters (in cm): radius, z fwd, z bwd)");
556 REGISTER_VARIABLE(
"helixExtPhi", trackHelixExtPhi,
557 "Returns phi of extrapolated helix parameters (parameters (in cm): radius, z fwd, z bwd)");
559 REGISTER_VARIABLE(
"nExtraCDCHits", nExtraCDCHits,
"[Eventbased] The number of CDC hits in the event not assigned to any track");
560 REGISTER_VARIABLE(
"nExtraCDCHitsPostCleaning", nExtraCDCHitsPostCleaning,
561 "[Eventbased] The number of CDC hits in the event not assigned to any track nor very likely beam background (i.e. hits that survive a cleanup selection)");
562 REGISTER_VARIABLE(
"hasExtraCDCHitsInLayer(i)", hasExtraCDCHitsInLayer,
563 "[Eventbased] Returns 1 if a non-assigned hit exists in the specified CDC layer");
564 REGISTER_VARIABLE(
"hasExtraCDCHitsInSuperLayer(i)", hasExtraCDCHitsInSuperLayer,
565 "[Eventbased] Returns 1 if a non-assigned hit exists in the specified CDC SuperLayer");
566 REGISTER_VARIABLE(
"nExtraCDCSegments", nExtraCDCSegments,
"[Eventbased] The number of CDC segments not assigned to any track");
573 REGISTER_VARIABLE(
"trackFindingFailureFlag", trackFindingFailureFlag,
574 "[Eventbased] A flag set by the tracking if there is reason to assume there was a track in the event missed by the tracking, "
575 "or the track finding was (partly) aborted for this event.");