70 {
71
72 const Trk::Track* fullTrack = &
track;
73 const Trk::Track* refittedTrack = nullptr;
74
76
77 if (!
m_fitter.empty() &&
track.perigeeParameters() &&
track.measurementsOnTrack() &&
track.measurementsOnTrack()->size() > 3) {
80 std::vector<const Trk::MeasurementBase*>
vec;
82 for (; measIter !=
track.measurementsOnTrack()->
end(); ++measIter) {
83 const Amg::Vector3D& position = (*measIter)->globalPosition();
85 vec.push_back((*measIter));
86 }
87 ATH_MSG_DEBUG(
"Refitting ID-part of track, meast vector of size " <<
vec.size());
88 refittedTrack =
m_fitter->fit(Gaudi::Hive::currentContext(),
90 *(
track.perigeeParameters()),
91 false,
93 } else {
95 Gaudi::Hive::currentContext(), track,
false,
Trk::muon).release();
96 }
97 fullTrack = refittedTrack;
98 }
99 ATH_MSG_DEBUG(
"Detected slimmed track, " << (refittedTrack ?
"refit was successful" :
"but could not refit it!"));
100
101 if (!refittedTrack) return ScatteringAngleSignificance(0);
102 }
103
104
105
106
108 unsigned measurements = 0;
109 double neighbourRadius = 0.;
110 double neighbourSignificance = 0.;
111 unsigned neighbours = 0;
112 unsigned scatterers = 0;
113 unsigned startMeasurement = 0;
114 std::vector<bool> isMeasurement;
115 std::vector<unsigned> measurementNumber;
116 std::vector<double> radii;
117 std::vector<double> sigmas;
118 double significance = 0.;
119 double totalSigma = 0.;
120 double weightedRadius = 0.;
123
124 if ((**s).measurementOnTrack()) {
125 ++measurements;
126 } else if (!measurements)
127 continue;
128
129
130 if (!(**s).materialEffectsOnTrack()) continue;
132 const Amg::Vector3D& position = (**s).materialEffectsOnTrack()->associatedSurface().globalReferencePoint();
134
135 const Trk::MaterialEffectsOnTrack* meot = dynamic_cast<const Trk::MaterialEffectsOnTrack*>((**s).materialEffectsOnTrack());
136 if (!meot) continue;
137
138
140 if (!scattering) continue;
141
142 ++scatterers;
143 double radius = (**s).trackParameters()->position().perp();
145 if ((**s).measurementOnTrack()) {
146 isMeasurement.push_back(true);
147 } else {
148 isMeasurement.push_back(false);
149 }
150 measurementNumber.push_back(measurements);
151
152
153
154 ++neighbours;
155 significance +=
sigma;
157 if (neighbours > 1 && measurements > startMeasurement) {
158 double norm = 1. / sqrt(
static_cast<double>(neighbours));
159 if (norm * fabs(significance) > fabs(neighbourSignificance)) {
160 if (significance == 0.) {
162 } else {
163 neighbourRadius = weightedRadius / significance;
164 }
165 neighbourSignificance =
norm * significance;
166 ATH_MSG_VERBOSE(
" current maximum for neighbourSignificance " << neighbourSignificance <<
" from " << neighbours
167 << " scatterers at radius " << neighbourRadius);
168 }
169 ATH_MSG_VERBOSE(
" reset neighbour after " << significance <<
" /sqrt(" << neighbours <<
") at radius " << radius);
170 neighbours = 0;
171 significance = 0.;
172 startMeasurement = measurements;
173 weightedRadius = 0.;
174 if (isMeasurement.back()) {
175 ++neighbours;
176 significance +=
sigma;
178 }
179 }
180
181 radii.push_back(radius);
182 sigmas.push_back(sigma);
184
185 ATH_MSG_VERBOSE(scatterers <<
" " << measurements <<
" " << isMeasurement.back() <<
" radius "
186 << (**s).trackParameters()->position().perp() <<
" deltaPhi " << scattering->
deltaPhi()
187 << " significance " << sigma << " integralSignificance " << totalSigma);
188 }
189
190
191 delete refittedTrack;
193 while (scatterers && measurementNumber.back() == measurements && !isMeasurement.back()) {
194 --scatterers;
195 totalSigma -= sigmas.back();
196 isMeasurement.pop_back();
197 measurementNumber.pop_back();
198 radii.pop_back();
199 sigmas.pop_back();
200 }
201 }
202
203
204 if (!scatterers) {
206 return ScatteringAngleSignificance(0);
207 }
208
209
210
211 double curvatureSignificance = totalSigma;
212 double curvatureRadius = *radii.rbegin();
213 unsigned index = scatterers - 1;
214
215 ATH_MSG_VERBOSE(
" scatteringAngleSignificance " << totalSigma <<
" at radius " << radii[index] <<
" at scatterer " << scatterers
216 << " sigmas.size() " << sigmas.size());
217
218 for (std::vector<double>::const_reverse_iterator
r = sigmas.rbegin();
r != sigmas.rend(); ++
r, --index) {
219 totalSigma -= 2. * (*r);
220 if (index > 1)
221 ATH_MSG_VERBOSE(
" scatteringAngleSignificance " << totalSigma <<
" at radius " << radii[index - 1] <<
" at scatterer "
222 << index);
223
224 if (fabs(totalSigma) > fabs(curvatureSignificance) && index > 1) {
225 curvatureSignificance = totalSigma;
226 curvatureRadius = radii[
index - 1];
227 }
228 }
229
230
231 curvatureSignificance /= sqrt(static_cast<double>(scatterers));
232
233 ATH_MSG_DEBUG(
" scatteringAngleSignificance" << std::setiosflags(std::ios::fixed) << std::setw(7) << std::setprecision(2)
234 << curvatureSignificance << " at radius " << std::setw(6) << std::setprecision(1)
235 << curvatureRadius << " neighbourSignificance " << std::setw(7) << std::setprecision(2)
236 << neighbourSignificance << " at radius " << std::setw(6) << std::setprecision(1)
237 << neighbourRadius << " from " << scatterers << " scatterers");
238
239 return ScatteringAngleSignificance(scatterers, curvatureRadius, curvatureSignificance, neighbourRadius, neighbourSignificance);
240 }
double charge(const T &p)
std::vector< size_t > vec
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
double charge() const
Returns the charge.
double sigmaDeltaPhi() const
returns the
double deltaPhi() const
returns the
@ CaloDeposit
This TSOS contains a CaloEnergy object.
const Trk::TrackStates * trackStateOnSurfaces() const
return a pointer to a const DataVector of const TrackStateOnSurfaces.
const Perigee * perigeeParameters() const
return Perigee.
Eigen::Matrix< double, 3, 1 > Vector3D
1.14.0