Belle II Software  release-05-01-25
SVDStateVarSet.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Nils Braun *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/ckf/svd/filters/states/SVDStateVarSet.h>
12 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory3D.h>
13 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory2D.h>
14 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectorySZ.h>
15 #include <tracking/dataobjects/RecoTrack.h>
16 
17 using namespace std;
18 using namespace Belle2;
19 using namespace TrackFindingCDC;
20 
21 namespace {
23  template<class APredicate>
24  double meanOver(const std::vector<TrackFindingCDC::WithWeight<const CKFToSVDState*>>& states, const APredicate& t)
25  {
26  double sum = 0;
27  unsigned int numberOfHits = 0;
28 
29  for (const CKFToSVDState* walkState : states) {
30  const SpacePoint* spacePoint = walkState->getHit();
31  if (not spacePoint) {
32  continue;
33  }
34  for (auto& cluster : spacePoint->getRelationsTo<SVDCluster>()) {
35  numberOfHits++;
36  sum += t(cluster);
37  }
38  }
39 
40  return sum / numberOfHits;
41  }
42 
44  template<class APredicate>
45  double minOver(const std::vector<TrackFindingCDC::WithWeight<const CKFToSVDState*>>& states, const APredicate& t)
46  {
47  double minimalValue = NAN;
48 
49  for (const CKFToSVDState* walkState : states) {
50  const SpacePoint* spacePoint = walkState->getHit();
51  if (not spacePoint) {
52  continue;
53  }
54  for (auto& cluster : spacePoint->getRelationsTo<SVDCluster>()) {
55  double currentValue = t(cluster);
56  if (std::isnan(minimalValue)) {
57  minimalValue = currentValue;
58  } else {
59  minimalValue = std::min(currentValue, minimalValue);
60  }
61  }
62  }
63 
64  return minimalValue;
65  }
66 
68  template<class APredicate>
69  double stdOver(const std::vector<TrackFindingCDC::WithWeight<const CKFToSVDState*>>& states, const APredicate& t)
70  {
71  double sum = 0;
72  double sumSquared = 0;
73  unsigned int numberOfHits = 0;
74 
75  for (const CKFToSVDState* walkState : states) {
76  const SpacePoint* spacePoint = walkState->getHit();
77  if (not spacePoint) {
78  continue;
79  }
80  for (auto& cluster : spacePoint->getRelationsTo<SVDCluster>()) {
81  numberOfHits++;
82  const auto& currentValue = t(cluster);
83  sum += currentValue;
84  sumSquared += currentValue * currentValue;
85  }
86  }
87 
88  return std::sqrt((sumSquared - sum / numberOfHits) / numberOfHits - 1);
89  }
90 }
91 
92 bool SVDStateVarSet::extract(const BaseSVDStateFilter::Object* pair)
93 {
94  const std::vector<TrackFindingCDC::WithWeight<const CKFToSVDState*>>& previousStates = pair->first;
95  CKFToSVDState* state = pair->second;
96 
97  const RecoTrack* cdcTrack = previousStates.front()->getSeed();
98  B2ASSERT("Path without seed?", cdcTrack);
99 
100  const SpacePoint* spacePoint = state->getHit();
101  B2ASSERT("Path without hit?", spacePoint);
102 
103  std::vector<TrackFindingCDC::WithWeight<const CKFToSVDState*>> allStates = previousStates;
104  allStates.emplace_back(state, 0);
105 
106  const std::vector<CDCHit*>& cdcHits = cdcTrack->getSortedCDCHitList();
107 
108  var<named("seed_cdc_hits")>() = cdcHits.size();
109  var<named("seed_lowest_cdc_layer")>() = cdcHits.front()->getICLayer();
110 
111  genfit::MeasuredStateOnPlane firstMeasurement;
112  if (state->mSoPSet()) {
113  firstMeasurement = state->getMeasuredStateOnPlane();
114  } else {
115  B2ASSERT("Previous state was not fitted?", previousStates.back()->mSoPSet());
116  firstMeasurement = previousStates.back()->getMeasuredStateOnPlane();
117  }
118 
119  Vector3D position = Vector3D(firstMeasurement.getPos());
120  Vector3D momentum = Vector3D(firstMeasurement.getMom());
121 
122  const CDCTrajectory3D trajectory(position, 0, momentum, cdcTrack->getChargeSeed());
123 
124  const Vector3D& hitPosition = static_cast<Vector3D>(spacePoint->getPosition());
125 
126  const double arcLength = trajectory.calcArcLength2D(hitPosition);
127  const Vector2D& trackPositionAtHit2D = trajectory.getTrajectory2D().getPos2DAtArcLength2D(arcLength);
128  const double trackPositionAtHitZ = trajectory.getTrajectorySZ().mapSToZ(arcLength);
129 
130  Vector3D trackPositionAtHit(trackPositionAtHit2D, trackPositionAtHitZ);
131 
132  const auto calculateCharge = [](const auto & s) {
133  return s.getCharge();
134  };
135  const auto calculateSeedCharge = [](const auto & s) {
136  return s.getSeedCharge();
137  };
138  const auto calculateSize = [](const auto & s) {
139  return s.getSize();
140  };
141  const auto calculateSNR = [](const auto & s) {
142  return s.getSNR();
143  };
144  const auto calculateChargeSizeRatio = [](const auto & s) {
145  return s.getCharge() / s.getSize();
146  };
147 
148  if (spacePoint->getType() == VXD::SensorInfoBase::SensorType::SVD) {
149  const auto& clusters = spacePoint->getRelationsTo<SVDCluster>();
150 
151  B2ASSERT("Must be related to exactly 2 clusters", clusters.size() == 2);
152  const SVDCluster* firstCluster = clusters[0];
153  const SVDCluster* secondCluster = clusters[1];
154 
155  var<named("cluster_1_charge")>() = calculateCharge(*firstCluster);
156  var<named("cluster_2_charge")>() = calculateCharge(*secondCluster);
157  var<named("mean_rest_cluster_charge")>() = meanOver(allStates, calculateCharge);
158  var<named("min_rest_cluster_charge")>() = minOver(allStates, calculateCharge);
159  var<named("std_rest_cluster_charge")>() = stdOver(allStates, calculateCharge);
160 
161  var<named("cluster_1_seed_charge")>() = calculateSeedCharge(*firstCluster);
162  var<named("cluster_2_seed_charge")>() = calculateSeedCharge(*secondCluster);
163  var<named("mean_rest_cluster_seed_charge")>() = meanOver(allStates, calculateSeedCharge);
164  var<named("min_rest_cluster_seed_charge")>() = minOver(allStates, calculateSeedCharge);
165  var<named("std_rest_cluster_seed_charge")>() = stdOver(allStates, calculateSeedCharge);
166 
167  var<named("cluster_1_size")>() = calculateSize(*firstCluster);
168  var<named("cluster_2_size")>() = calculateSize(*secondCluster);
169  var<named("mean_rest_cluster_size")>() = meanOver(allStates, calculateSize);
170  var<named("min_rest_cluster_size")>() = meanOver(allStates, calculateSize);
171  var<named("std_rest_cluster_size")>() = stdOver(allStates, calculateSize);
172 
173  var<named("cluster_1_snr")>() = calculateSNR(*firstCluster);
174  var<named("cluster_2_snr")>() = calculateSNR(*secondCluster);
175  var<named("mean_rest_cluster_snr")>() = meanOver(allStates, calculateSNR);
176  var<named("min_rest_cluster_snr")>() = minOver(allStates, calculateSNR);
177  var<named("std_rest_cluster_snr")>() = stdOver(allStates, calculateSNR);
178 
179  var<named("cluster_1_charge_over_size")>() = calculateChargeSizeRatio(*firstCluster);
180  var<named("cluster_2_charge_over_size")>() = calculateChargeSizeRatio(*secondCluster);
181  var<named("mean_rest_cluster_charge_over_size")>() = meanOver(allStates, calculateChargeSizeRatio);
182  var<named("min_rest_cluster_charge_over_size")>() = minOver(allStates, calculateChargeSizeRatio);
183  var<named("std_rest_cluster_charge_over_size")>() = stdOver(allStates, calculateChargeSizeRatio);
184  }
185 
186  std::vector<const SpacePoint*> spacePoints;
187  for (const CKFToSVDState* walkState : allStates) {
188  const SpacePoint* sp = walkState->getHit();
189  if (sp) {
190  spacePoints.push_back(sp);
191  }
192  }
193 
194  if (spacePoints.size() >= 3) {
195  var<named("quality_index_triplet")>() = m_qualityTriplet.estimateQuality(spacePoints);
196  var<named("quality_index_circle")>() = m_qualityCircle.estimateQuality(spacePoints);
197  var<named("quality_index_helix")>() = m_qualityHelix.estimateQuality(spacePoints);
198  } else {
199  var<named("quality_index_triplet")>() = 0;
200  var<named("quality_index_circle")>() = 0;
201  var<named("quality_index_helix")>() = 0;
202  }
203 
204  return true;
205 }
Belle2::Vector3D
HepGeom::Vector3D< double > Vector3D
3D Vector
Definition: Cell.h:35
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
Belle2::CKFState::getMeasuredStateOnPlane
const genfit::MeasuredStateOnPlane & getMeasuredStateOnPlane() const
Get the mSoP if already set during extrapolation (or fitting)
Definition: CKFState.h:103
Belle2::CKFState::mSoPSet
bool mSoPSet() const
Is the mSoP already set? (= state was already extrapolated)
Definition: CKFState.h:116
Belle2::RelationsInterface::getRelationsTo
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
Definition: RelationsObject.h:199
Belle2::SpacePoint
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
Definition: SpacePoint.h:52
Belle2::SpacePoint::getPosition
const B2Vector3< double > & getPosition() const
return the position vector in global coordinates
Definition: SpacePoint.h:148
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::RecoTrack::getSortedCDCHitList
std::vector< Belle2::RecoTrack::UsedCDCHit * > getSortedCDCHitList() const
Return a sorted list of cdc hits. Sorted by the sortingParameter.
Definition: RecoTrack.h:461
Belle2::TrackFindingCDC::WithWeight
A mixin class to attach a weight to an object.
Definition: WithWeight.h:34
Belle2::CKFState::getHit
const Hit * getHit() const
Return the SP this state is related to. May be nullptr.
Definition: CKFState.h:76
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::RecoTrack::getChargeSeed
short int getChargeSeed() const
Return the charge seed stored in the reco track. ATTENTION: This is not the fitted charge.
Definition: RecoTrack.h:497
Belle2::CKFToSVDState
Specialized CKF State for extrapolating into the SVD.
Definition: CKFToSVDState.h:29
Belle2::SpacePoint::getType
Belle2::VXD::SensorInfoBase::SensorType getType() const
Return SensorType (PXD, SVD, ...) on which the SpacePoint lives.
Definition: SpacePoint.h:155
Belle2::TrackFindingCDC::Filter::Object
AObject Object
Type of the object to be analysed.
Definition: Filter.dcl.h:43