Project
Loading...
Searching...
No Matches
Smoother.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11//
12// \author matteo.concas@cern.ch
13
15
16namespace o2
17{
18namespace its
19{
20
21constexpr std::array<double, 3> getInverseSymm2D(const std::array<double, 3>& mat)
22{
23 const double det = mat[0] * mat[2] - mat[1] * mat[1];
24 return std::array<double, 3>{mat[2] / det, -mat[1] / det, mat[0] / det};
25}
26
27// Smoother
28// template <unsigned int D>
29// Smoother<D>::Smoother(TrackITSExt& track, size_t smoothingLayer, const ROframe& event, float bZ, o2::base::PropagatorF::MatCorrType corr) : mLayerToSmooth{smoothingLayer},
30// mBz(bZ),
31// mCorr(corr)
32// {
33//
34// auto propInstance = o2::base::Propagator::Instance();
35// const TrackingFrameInfo& originalTf = event.getTrackingFrameInfoOnLayer(mLayerToSmooth).at(track.getClusterIndex(mLayerToSmooth));
36//
37// mOutwardsTrack = track; // This track will be propagated outwards inside the smoother! (as last step of fitting did inward propagation)
38// mInwardsTrack = {track.getParamOut(), // This track will be propagated inwards inside the smoother!
39// static_cast<short>(mOutwardsTrack.getNumberOfClusters()), -999, static_cast<std::uint32_t>(event.getROFrameId()),
40// mOutwardsTrack.getParamOut(), mOutwardsTrack.getClusterIndexes()};
41//
42// mOutwardsTrack.resetCovariance();
43// mOutwardsTrack.setChi2(0);
44// mInwardsTrack.resetCovariance();
45// mInwardsTrack.setChi2(0);
46//
47// bool statusOutw{false};
48// bool statusInw{false};
49//
50// //////////////////////
51// // Outward propagation
52// for (size_t iLayer{0}; iLayer < mLayerToSmooth; ++iLayer) {
53// if (mOutwardsTrack.getClusterIndex(iLayer) == constants::UnusedIndex) { // Shorter tracks
54// continue;
55// }
56// const TrackingFrameInfo& tF = event.getTrackingFrameInfoOnLayer(iLayer).at(mOutwardsTrack.getClusterIndex(iLayer));
57// statusOutw = mOutwardsTrack.rotate(tF.alphaTrackingFrame);
58// statusOutw &= propInstance->propagateToX(mOutwardsTrack,
59// tF.xTrackingFrame,
60// mBz,
61// o2::base::PropagatorImpl<float>::MAX_SIN_PHI,
62// o2::base::PropagatorImpl<float>::MAX_STEP,
63// mCorr);
64// mOutwardsTrack.setChi2(mOutwardsTrack.getChi2() + mOutwardsTrack.getPredictedChi2(tF.positionTrackingFrame, tF.covarianceTrackingFrame));
65// statusOutw &= mOutwardsTrack.o2::track::TrackParCov::update(tF.positionTrackingFrame, tF.covarianceTrackingFrame);
66// // LOG(info) << "Outwards loop on inwards track, layer: " << iLayer << " x: " << mOutwardsTrack.getX();
67// }
68//
69// // Prediction on the previously outwards-propagated track is done on a copy, as the process seems to be not reversible
70// auto outwardsClone = mOutwardsTrack;
71// statusOutw = outwardsClone.rotate(originalTf.alphaTrackingFrame);
72// statusOutw &= propInstance->propagateToX(outwardsClone,
73// originalTf.xTrackingFrame,
74// mBz,
75// o2::base::PropagatorImpl<float>::MAX_SIN_PHI,
76// o2::base::PropagatorImpl<float>::MAX_STEP,
77// mCorr);
78// /////////////////////
79// // Inward propagation
80// for (size_t iLayer{D - 1}; iLayer > mLayerToSmooth; --iLayer) {
81// if (mInwardsTrack.getClusterIndex(iLayer) == constants::UnusedIndex) { // Shorter tracks
82// continue;
83// }
84// const TrackingFrameInfo& tF = event.getTrackingFrameInfoOnLayer(iLayer).at(mInwardsTrack.getClusterIndex(iLayer));
85// statusInw = mInwardsTrack.rotate(tF.alphaTrackingFrame);
86// statusInw &= propInstance->propagateToX(mInwardsTrack,
87// tF.xTrackingFrame,
88// mBz,
89// o2::base::PropagatorImpl<float>::MAX_SIN_PHI,
90// o2::base::PropagatorImpl<float>::MAX_STEP,
91// mCorr);
92// mInwardsTrack.setChi2(mInwardsTrack.getChi2() + mInwardsTrack.getPredictedChi2(tF.positionTrackingFrame, tF.covarianceTrackingFrame));
93// statusInw &= mInwardsTrack.o2::track::TrackParCov::update(tF.positionTrackingFrame, tF.covarianceTrackingFrame);
94// // LOG(info) << "Inwards loop on outwards track, layer: " << iLayer << " x: " << mInwardsTrack.getX();
95// }
96//
97// // Prediction on the previously inwards-propagated track is done on a copy, as the process seems to be not revesible
98// auto inwardsClone = mInwardsTrack;
99// statusInw = inwardsClone.rotate(originalTf.alphaTrackingFrame);
100// statusInw &= propInstance->propagateToX(inwardsClone,
101// originalTf.xTrackingFrame,
102// mBz,
103// o2::base::PropagatorImpl<float>::MAX_SIN_PHI,
104// o2::base::PropagatorImpl<float>::MAX_STEP,
105// mCorr);
106// // Compute weighted local chi2
107// mInitStatus = statusInw && statusOutw;
108// if (mInitStatus) {
109// mBestChi2 = computeSmoothedPredictedChi2(inwardsClone, outwardsClone, originalTf.positionTrackingFrame, originalTf.covarianceTrackingFrame);
110// mLastChi2 = mBestChi2;
111// LOG(info) << "Smoothed chi2 on original cluster: " << mBestChi2;
112// }
113// }
114
115template <unsigned int D>
116Smoother<D>::~Smoother() = default;
117
118template <unsigned int D>
119float Smoother<D>::computeSmoothedPredictedChi2(const o2::track::TrackParCov& firstTrack, // outwards track: from innermost cluster to outermost
120 const o2::track::TrackParCov& secondTrack, // inwards track: from outermost cluster to innermost
121 const std::array<float, 2>& cls,
122 const std::array<float, 3>& clCov)
123{
124 // Tracks need to be already propagated, compute only chi2
125 // Symmetric covariances assumed
126
127 if (firstTrack.getX() != secondTrack.getX()) {
128 LOG(fatal) << "Tracks need to be propagated to the same point! secondTrack.X=" << secondTrack.getX() << " firstTrack.X=" << firstTrack.getX();
129 }
130
131 std::array<double, 2> pp1 = {static_cast<double>(firstTrack.getY()), static_cast<double>(firstTrack.getZ())}; // P1: predicted Y,Z points
132 std::array<double, 2> pp2 = {static_cast<double>(secondTrack.getY()), static_cast<double>(secondTrack.getZ())}; // P2: predicted Y,Z points
133
134 std::array<double, 3> c1 = {static_cast<double>(firstTrack.getSigmaY2()),
135 static_cast<double>(firstTrack.getSigmaZY()),
136 static_cast<double>(firstTrack.getSigmaZ2())}; // Cov. track 1
137
138 std::array<double, 3> c2 = {static_cast<double>(secondTrack.getSigmaY2()),
139 static_cast<double>(secondTrack.getSigmaZY()),
140 static_cast<double>(secondTrack.getSigmaZ2())}; // Cov. track 2
141
142 std::array<double, 3> w1 = getInverseSymm2D(c1); // weight matrices
143 std::array<double, 3> w2 = getInverseSymm2D(c2);
144
145 std::array<double, 3> w1w2 = {w1[0] + w2[0], w1[1] + w2[1], w1[2] + w2[2]}; // (W1 + W2)
146 std::array<double, 3> C = getInverseSymm2D(w1w2); // C = (W1+W2)^-1
147
148 std::array<double, 2> w1pp1 = {w1[0] * pp1[0] + w1[1] * pp1[1], w1[1] * pp1[0] + w1[2] * pp1[1]}; // W1 * P1
149 std::array<double, 2> w2pp2 = {w2[0] * pp2[0] + w2[1] * pp2[1], w2[1] * pp2[0] + w2[2] * pp2[1]}; // W2 * P2
150
151 double Y = C[0] * (w1pp1[0] + w2pp2[0]) + C[1] * (w1pp1[1] + w2pp2[1]); // Pp: weighted normalized combination of the predictions:
152 double Z = C[1] * (w1pp1[0] + w2pp2[0]) + C[2] * (w1pp1[1] + w2pp2[1]); // Pp = [(W1 * P1) + (W2 * P2)] / (W1 + W2)
153
154 std::array<double, 2> delta = {Y - cls[0], Z - cls[1]}; // Δ = Pp - X, X: space point of cluster (Y,Z)
155 std::array<double, 3> CCp = {C[0] + static_cast<double>(clCov[0]), C[1] + static_cast<double>(clCov[1]), C[2] + static_cast<double>(clCov[2])}; // Transformation of cluster covmat: CCp = C + Cov
156 std::array<double, 3> Wp = getInverseSymm2D(CCp); // Get weight matrix: Wp = CCp^-1
157
158 float chi2 = static_cast<float>(delta[0] * (Wp[0] * delta[0] + Wp[1] * delta[1]) + delta[1] * (Wp[1] * delta[0] + Wp[2] * delta[1])); // chi2 = tΔ * (Wp * Δ)
159
160 // #ifdef CA_DEBUG
161 LOG(info) << "Cluster_y: " << cls[0] << " Cluster_z: " << cls[1];
162 LOG(info) << "\t\t- Covariance cluster: Y2: " << clCov[0] << " YZ: " << clCov[1] << " Z2: " << clCov[2];
163 LOG(info) << "\t\t- Propagated t1_y: " << pp1[0] << " t1_z: " << pp1[1];
164 LOG(info) << "\t\t- Propagated t2_y: " << pp2[0] << " t2_z: " << pp2[1];
165 LOG(info) << "\t\t- Covariance t1: sY2: " << c1[0] << " sYZ: " << c1[1] << " sZ2: " << c1[2];
166 LOG(info) << "\t\t- Covariance t2: sY2: " << c2[0] << " sYZ: " << c2[1] << " sZ2: " << c2[2];
167 LOG(info) << "Smoother prediction Y: " << Y << " Z: " << Z;
168 LOG(info) << "\t\t- Delta_y: " << delta[0] << " Delta_z: " << delta[1];
169 LOG(info) << "\t\t- Covariance Pr: Y2: " << C[0] << " YZ: " << C[1] << " Z2: " << C[2];
170 LOG(info) << "\t\t- predicted chi2 t1: " << firstTrack.getPredictedChi2(cls, clCov);
171 LOG(info) << "\t\t- predicted chi2 t2: " << secondTrack.getPredictedChi2(cls, clCov);
172 // #endif
173 return chi2;
174}
175
176// template <unsigned int D>
177// bool Smoother<D>::testCluster(const int clusterId, const ROframe& event)
178// {
179// if (!mInitStatus) {
180// return false;
181// }
182// auto propInstance = o2::base::Propagator::Instance();
183// const TrackingFrameInfo& testTf = event.getTrackingFrameInfoOnLayer(mLayerToSmooth).at(clusterId);
184//
185// bool statusOutw{false};
186// bool statusInw{false};
187//
188// // Prediction on the previously outwards-propagated track is done on a copy, as the process seems to be not reversible
189// auto outwardsClone = mOutwardsTrack;
190// statusOutw = outwardsClone.rotate(testTf.alphaTrackingFrame);
191// statusOutw &= propInstance->propagateToX(outwardsClone,
192// testTf.xTrackingFrame,
193// mBz,
194// o2::base::PropagatorImpl<float>::MAX_SIN_PHI,
195// o2::base::PropagatorImpl<float>::MAX_STEP,
196// mCorr);
197//
198// // Prediction on the previously inwards-propagated track is done on a copy, as the process seems to be not reversible
199// auto inwardsClone = mInwardsTrack;
200// statusInw = inwardsClone.rotate(testTf.alphaTrackingFrame);
201// statusInw &= propInstance->propagateToX(inwardsClone,
202// testTf.xTrackingFrame,
203// mBz,
204// o2::base::PropagatorImpl<float>::MAX_SIN_PHI,
205// o2::base::PropagatorImpl<float>::MAX_STEP,
206// mCorr);
207// if (!(statusOutw && statusInw)) {
208// LOG(warning) << "Failed propagation in smoother!";
209// return false;
210// }
211//
212// // Compute weighted local chi2
213// mLastChi2 = computeSmoothedPredictedChi2(inwardsClone, outwardsClone, testTf.positionTrackingFrame, testTf.covarianceTrackingFrame);
214// LOG(info) << "Smoothed chi2 on tested cluster: " << mLastChi2;
215//
216// return true;
217// }
218
219template class Smoother<7>;
220
221} // namespace its
222} // namespace o2
bool const GPUTPCGMMerger::trackCluster * c1
bool const GPUTPCGMMerger::trackCluster const clcomparestruct * c2
Class to handle Kalman smoothing for ITS tracking. Its instance stores the state of the track to the ...
constexpr std::array< double, 3 > getInverseSymm2D(const std::array< double, 3 > &mat)
Definition Smoother.cxx:21
@ C
Definition Defs.h:36
TrackParCovF TrackParCov
Definition Track.h:33
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"