ATLAS Offline Software
Loading...
Searching...
No Matches
MuonPhiHitSelector.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "MuonSegmentCleaner/MuonPhiHitSelector.h"
6
7#include <ctime>
8#include <iostream>
9#include <list>
10#include <sstream>
11#include <vector>
12
13#include "CxxUtils/sincos.h"
14#include "MuonCompetingRIOsOnTrack/CompetingMuonClustersOnTrack.h"
15#include "MuonPrepRawData/MuonCluster.h"
16#include "MuonRIO_OnTrack/CscClusterOnTrack.h"
17#include "MuonRIO_OnTrack/MdtDriftCircleOnTrack.h"
18#include "MuonRIO_OnTrack/MuonClusterOnTrack.h"
19#include "MuonRIO_OnTrack/RpcClusterOnTrack.h"
20#include "MuonRIO_OnTrack/TgcClusterOnTrack.h"
21#include "MuonReadoutGeometry/CscReadoutElement.h"
22#include "MuonReadoutGeometry/MdtReadoutElement.h"
23#include "MuonReadoutGeometry/RpcReadoutElement.h"
24#include "MuonReadoutGeometry/TgcReadoutElement.h"
25#include "TrkPrepRawData/PrepRawData.h"
26#include "TrkRIO_OnTrack/RIO_OnTrack.h"
27#include "TrkSurfaces/PlaneSurface.h"
28
29MuonPhiHitSelector::MuonPhiHitSelector(const std::string& type, const std::string& name, const IInterface* parent)
30 : AthAlgTool(type, name, parent)
31{
32 declareInterface<IMuonHitSelector>(this);
33
34
35}
36
37StatusCode
38MuonPhiHitSelector::initialize()
39{
40 ATH_MSG_VERBOSE("MuonPhiHitSelector::Initializing");
41 ATH_CHECK(m_competingRIOsOnTrackTool.retrieve());
42 ATH_CHECK(m_clusterCreator.retrieve());
43 ATH_CHECK(m_idHelperSvc.retrieve());
44 ATH_CHECK(m_cscRotCreator.retrieve(EnableTool{!m_cscRotCreator.empty()}));
45 ATH_MSG_VERBOSE("End of Initializing");
46 return StatusCode::SUCCESS;
47}
48
49std::vector<std::unique_ptr<const Trk::MeasurementBase>>
50MuonPhiHitSelector::select_rio(const double pmom, const std::vector<const Trk::RIO_OnTrack*>& associatedHits,
51 const std::vector<const Trk::PrepRawData*>& unassociatedHits) const
52{
53 // Make Rios on Track
54
55 int time_start = std::clock() / 1000;
56
57 std::vector<std::unique_ptr<const Trk::MeasurementBase>> selectedHits{}, selectedClusters{};
58
59 ATH_MSG_VERBOSE("Executing MuonPhiHitSelectorTool select_rio ");
60
61 int nhits = associatedHits.size() + unassociatedHits.size();
62
63 ATH_MSG_DEBUG("Executing MuonPhiHitSelectorTool nhits select_rio " << nhits);
64
65 std::vector<double> phiHitx(nhits);
66 std::vector<double> phiHity(nhits);
67 std::vector<double> phiHitz(nhits);
68 std::vector<double> phiError(nhits);
69 std::vector<Identifier> phiId(nhits);
70 std::vector<double> phiPull(nhits);
71 std::vector<int> phiSelect(nhits);
72 std::vector<int> phiMult(nhits);
73 std::vector<int> quality(nhits);
74 std::vector<const Trk::PrepRawData*> phiPrep(nhits);
75
76 std::map<Identifier, int> phiMapId;
77 int nphi = 0;
78
79 for (const Trk::RIO_OnTrack* rot : associatedHits) {
80 Identifier id = rot->identify();
81 phiId[nphi] = id;
82 Amg::Vector3D gHitPos = rot->globalPosition();
83 if (m_idHelperSvc->isRpc(id)) {
84 phiSelect[nphi] = 1;
85 } else if (m_idHelperSvc->isTgc(id)) {
86 phiSelect[nphi] = 2;
87 } else if (m_idHelperSvc->isCsc(id)) {
88 phiSelect[nphi] = 3;
89 }
90 phiHitx[nphi] = gHitPos.x();
91 phiHity[nphi] = gHitPos.y();
92 phiHitz[nphi] = gHitPos.z();
93
94 const Amg::MatrixX& cov = rot->localCovariance();
95 double error = cov(0, 0);
96
97 // for the TGCs diagonalize the error matrix and use the smallest component
98 if (cov.cols() > 1) {
99 AmgSymMatrix(2) Er{AmgSymMatrix(2)::Zero()};
100 Er(0, 0) = cov(0, 0);
101 Er(0, 1) = cov(0, 1);
102 Er(1, 1) = cov(1, 1);
103 Er(1, 0) = Er(0, 1);
104
105 double chi = Er(0, 0) != Er(1, 1) ? std::atan(-2 * Er(0, 1) / (Er(0, 0) - Er(1, 1))) / 2. : 0.;
106
107 CxxUtils::sincos scchi(chi);
108
109 AmgSymMatrix(2) Rot{AmgSymMatrix(2)::Zero()};
110 Rot(0, 0) = scchi.cs;
111 Rot(1, 1) = Rot(0, 0);
112 Rot(0, 1) = scchi.sn;
113 Rot(1, 0) = -Rot(0, 1);
114 AmgMatrix(2, 2) D = Rot.transpose() * Er * Rot;
115 ATH_MSG_DEBUG("Diagonalized error matrix " << D);
116 error = std::min(D(0, 0),D(1, 1));
117 }
118 phiError[nphi] = std::sqrt(error);
119 quality[nphi] = 1000;
120 phiMapId[id] = 1;
121 phiPrep[nphi] = rot->prepRawData();
122 double phipos = std::atan2(phiHity[nphi], phiHitx[nphi]);
123 ATH_MSG_DEBUG("phi Segment Hit " << nphi << " det " << phiSelect[nphi] << " phi " << phipos);
124 nphi++;
125 }
126 int nphiseg = nphi;
127
128 for (const Trk::PrepRawData* prd : unassociatedHits) {
129 Identifier id = prd->identify();
130 phiId[nphi] = id;
131 // Skip phi hits already on segment
132 if (phiMapId.count(id)) continue;
133 const Muon::MuonCluster* clus = dynamic_cast<const Muon::MuonCluster*>(prd);
134 if (!clus) continue;
135 if (m_idHelperSvc->isRpc(id))
136 phiSelect[nphi] = 1;
137 else if (m_idHelperSvc->isTgc(id))
138 phiSelect[nphi] = 2;
139 else if (m_idHelperSvc->isCsc(id))
140 phiSelect[nphi] = 3;
141 Amg::Vector3D gHitPos = clus->globalPosition();
142 phiHitx[nphi] = gHitPos.x();
143 phiHity[nphi] = gHitPos.y();
144 phiHitz[nphi] = gHitPos.z();
145 phiError[nphi] = prd->localCovariance()(Trk::locX);
146 quality[nphi] = 10;
147 phiPrep[nphi] = prd;
148 double phipos = std::atan2(phiHity[nphi], phiHitx[nphi]);
149 ATH_MSG_DEBUG("phi Pattern Hit " << nphi << " phi " << phipos);
150 nphi++;
151 }
152
153 double chi2(0);
154 double r0(0);
155 int nfit;
156 std::vector<double> errorM(4);
157 double phi(DBL_MAX);
158 fitRecPhi(pmom, phiId, phiHitx, phiHity, phiHitz, phiError, quality, nphi, phiPull, phiMult, phiSelect, chi2, r0,
159 phi, errorM, nfit);
160
161 // Define global track parameters (not used 27-8 JS)
162
163 for (int i = 0; i < nphi; ++i) {
164 if (phiSelect[i] > 0) {
165 if (phiSelect[i] == 1) {
166 const Muon::RpcPrepData* prd = dynamic_cast<const Muon::RpcPrepData*>(phiPrep[i]);
167 const Amg::Vector3D globalpos(phiHitx[i], phiHity[i], phiHitz[i]);
168 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
169 if (rio) selectedHits.push_back(std::move(rio));
170 } else if (phiSelect[i] == 2) {
171 const Muon::TgcPrepData* prd = dynamic_cast<const Muon::TgcPrepData*>(phiPrep[i]);
172 const Amg::Vector3D globalpos(phiHitx[i], phiHity[i], phiHitz[i]);
173 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
174 if (rio) selectedHits.push_back(std::move(rio));
175 } else if (phiSelect[i] == 3) {
176 const Muon::CscPrepData* prd = dynamic_cast<const Muon::CscPrepData*>(phiPrep[i]);
177 const Amg::Vector3D globalpos(phiHitx[i], phiHity[i], phiHitz[i]);
178 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_cscRotCreator->createRIO_OnTrack(*prd, globalpos)};
179 if (rio) selectedHits.push_back(std::move(rio));
180 }
181
182 ATH_MSG_DEBUG("Make ONE rio per PrepData");
183 }
184 }
185 ATH_MSG_DEBUG("Fit hit results phi " << phi << " chi2 " << chi2 << " segment hits " << nphiseg
186 << " pattern hits " << nphi - nphiseg << " nfit " << nfit << " rio size "
187 << selectedHits.size());
188
189
190 std::vector<double> clusterX(nphi);
191 std::vector<double> clusterY(nphi);
192 std::vector<double> clusterZ(nphi);
193 std::vector<double> clusterError(nphi);
194 std::vector<Identifier> clusterId(nphi);
195 std::vector<int> clusterHits(nphi);
196 std::vector<double> clusterPull(nphi);
197 std::vector<int> clusterSelect(nphi);
198 // Link from hit to cluster
199 std::vector<int> clusterInt(nphi);
200
201 int ncl, imax;
202 double chi2cl, r0cl, phicl;
203 std::vector<double> errorMcl(4);
204 clusterPhi(phiId, phiHitx, phiHity, phiHitz, phiError, phiPull, phiSelect, nphi, clusterX, clusterY, clusterZ,
205 clusterError, clusterId, clusterHits, clusterSelect, clusterInt, ncl);
206
207
208 for (int ic = 0; ic < ncl; ++ic) {
209 std::list<const Trk::PrepRawData*> prdList;
210 int iic = -1;
211 double avError = 0.;
212 int ip = -1;
213 int np = 0;
214 for (int i = 0; i < nphi; ++i) {
215 if (clusterInt[i] == ic) {
216 ip = i;
217 prdList.push_back(phiPrep[i]);
218 avError += 1. / (phiError[i] * phiError[i]);
219 if (clusterId[ic] == phiId[i]) iic = i;
220 np++;
221 }
222 }
223 if (iic > -1) {
224 ATH_MSG_DEBUG("Phi cluster found np " << np << " ip " << ip);
225 if (np == 1) {
226 // Only one PrepData: create RIO on Track
227 const Amg::Vector3D globalpos(clusterX[ic], clusterY[ic], clusterZ[ic]);
228 if (phiSelect[ip] == 1) {
229 ATH_MSG_DEBUG("Phi RPC rio");
230 const Muon::RpcPrepData* prd = dynamic_cast<const Muon::RpcPrepData*>(phiPrep[ip]);
231 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
232 if (rio) selectedClusters.push_back(std::move(rio));
233 } else if (phiSelect[ip] == 2) {
234 ATH_MSG_DEBUG("Phi TGC rio");
235 const Muon::TgcPrepData* prd = dynamic_cast<const Muon::TgcPrepData*>(phiPrep[ip]);
236 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
237 if (rio) selectedClusters.push_back(std::move(rio));
238 } else if (phiSelect[ip] == 3) {
239 ATH_MSG_DEBUG("Phi CSC rio");
240 const Muon::CscPrepData* prd = dynamic_cast<const Muon::CscPrepData*>(phiPrep[ip]);
241 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_cscRotCreator->createRIO_OnTrack(*prd, globalpos)};
242 if (rio) selectedClusters.push_back(std::move(rio));
243 }
244 } else {
245
246 if (m_competingRios) {
247 // More PrepData's: create Competing RIOs on Track
248 avError = std::sqrt(1. / avError);
249 double scaleFactor = clusterError[ic] / avError;
250 std::unique_ptr<const Trk::CompetingRIOsOnTrack> rio =
251 m_competingRIOsOnTrackTool->createBroadCluster(prdList, scaleFactor);
252 selectedClusters.push_back(std::move(rio));
253 if (msgLvl(MSG::DEBUG))
254 ATH_MSG_DEBUG("Make competing rio/cluster "
255 << " scale factor " << scaleFactor << " number of rios " << prdList.size());
256 } else {
257 // Make one Rio for central cluster
258 ip = iic;
259 const Amg::Vector3D globalpos(clusterX[ic], clusterY[ic], clusterZ[ic]);
260 if (phiSelect[ip] == 1) {
261 ATH_MSG_DEBUG("Phi RPC rio central cluster");
262 const Muon::RpcPrepData* prd = dynamic_cast<const Muon::RpcPrepData*>(phiPrep[ip]);
263 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
264 if (rio) selectedClusters.push_back(std::move(rio));
265 } else if (phiSelect[ip] == 2) {
266 ATH_MSG_DEBUG("Phi TGC rio central cluster");
267 const Muon::TgcPrepData* prd = dynamic_cast<const Muon::TgcPrepData*>(phiPrep[ip]);
268 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{ m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
269 if (rio) selectedClusters.push_back(std::move(rio));
270 } else if (phiSelect[ip] == 3) {
271 ATH_MSG_DEBUG("Phi CSC rio central cluster");
272 const Muon::CscPrepData* prd = dynamic_cast<const Muon::CscPrepData*>(phiPrep[ip]);
273 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
274 if (rio) selectedClusters.push_back(std::move(rio));
275 }
276 }
277 }
278 } else {
279 ATH_MSG_DEBUG("Phi cluster NOT found ");
280 }
281 }
282
283 fitPhiSL(pmom, clusterId, clusterX, clusterY, clusterZ, clusterError, clusterSelect, ncl, clusterPull, imax, chi2cl,
284 r0cl, phicl, errorMcl, false);
285
286
287 if (msgLvl(MSG::DEBUG) || m_summary) {
288 ATH_MSG_DEBUG("PhiHitSelector Time spent " << std::clock() / 1000 - time_start << " nhits " << nhits
289 << " segment hits " << associatedHits.size() << " nfit " << nfit
290 << " nclusters " << ncl);
291 ATH_MSG_DEBUG("Fit cluster results phi " << phicl << " chi2 " << chi2cl << " number of clusters " << ncl
292 << " size cluster Hits " << selectedClusters.size());
293 }
294 if (m_makeClusters) {
295 return selectedClusters;
296 }
297 return selectedHits;
298}
299
300
301void
302MuonPhiHitSelector::clusterPhi(const std::vector<Identifier>& id, const std::vector<double>& hitx,
303 const std::vector<double>& hity, const std::vector<double>& hitz,
304 const std::vector<double>& error, const std::vector<double>& pull,
305 std::vector<int>& select, const int n, std::vector<double>& clusterX,
306 std::vector<double>& clusterY, std::vector<double>& clusterZ,
307 std::vector<double>& clusterError, std::vector<Identifier>& clusterId,
308 std::vector<int>& clusterHits, std::vector<int>& clusterSelect,
309 std::vector<int>& clusterInt, int& ncl) const
310{
311
312 //
313 // Use hits (select > 0) and pulls to make clusters
314 //
315 //
316 // Inputs
317 // id = identifiers hits
318 // hitx hity hitz = position in space
319 // error = associated error (in x-y plane)
320 // pull (from fit)= residual (hit -fit) /error
321 // select = > 0 for selected hits
322 // n = total number of hits
323
324 // Outputs
325 // clusterX Y Z = cluster positions
326 // clusterError = errors
327 // clusterId = cluster identifier (smallest pull)
328 // clusterHits = number of hits per cluster
329 // ncl = number of clusters
330 // chi2 = total chi2
331 // r0 = perigee parameter of fit (0,0)
332 // phi = azimuthal angle of fit at perigee
333
334
335 ATH_MSG_DEBUG("Start phi clustering");
336
337 ncl = 0;
338 if (n == 0) return;
339
340 std::vector<int> scode(n);
341
342 for (int i = 0; i < n; ++i) {
343 Identifier idi = id[i];
344 int code = 0;
345 if (m_idHelperSvc->isRpc(idi)) {
346 int doubZ = m_idHelperSvc->rpcIdHelper().doubletZ(idi);
347 int doubPhi = m_idHelperSvc->rpcIdHelper().doubletPhi(idi);
348 code = 100000000 * (m_idHelperSvc->rpcIdHelper().stationName(idi))
349 + 1000000 * (m_idHelperSvc->rpcIdHelper().stationPhi(idi))
350 + 10000 * ((m_idHelperSvc->rpcIdHelper().stationEta(idi)) + 1000);
351 code += 1000 * (doubZ - 1) + 100 * (doubPhi - 1);
352 code += 2 * ((m_idHelperSvc->rpcIdHelper().doubletR(idi)) - 1)
353 + 16 * ((m_idHelperSvc->rpcIdHelper().gasGap(idi)) - 1);
354 } else if (m_idHelperSvc->isTgc(idi)) {
355 code = 1000000 * (m_idHelperSvc->tgcIdHelper().stationName(idi))
356 + 10000 * (m_idHelperSvc->tgcIdHelper().stationPhi(idi))
357 + 100 * ((m_idHelperSvc->tgcIdHelper().stationEta(idi)) + 10);
358 code = code + m_idHelperSvc->tgcIdHelper().gasGap(idi);
359 } else if (m_idHelperSvc->isCsc(idi)) {
360 code = 1000000 * (m_idHelperSvc->cscIdHelper().stationName(idi))
361 + 10000 * (m_idHelperSvc->cscIdHelper().stationPhi(idi))
362 + 100 * ((m_idHelperSvc->cscIdHelper().stationEta(idi)) + 10);
363 code = code + m_idHelperSvc->cscIdHelper().wireLayer(idi);
364 }
365 scode[i] = code;
366 }
367
368 // std::vector<int> clusterInt(n);
369
370 for (int i = 0; i < n; ++i) {
371 clusterInt[i] = -1;
372 }
373
374 int icl = -1;
375 for (int i = 0; i < n; ++i) {
376 if (error[i] != 0 && select[i] > 0 && clusterInt[i] == -1) {
377 icl++;
378 clusterInt[i] = icl;
379 for (int j = i + 1; j < n; ++j) {
380 if (clusterInt[j] == -1) {
381 if (error[j] != 0 && select[j] > 0) {
382 if (scode[i] == scode[j]) clusterInt[j] = icl;
383 }
384 }
385 }
386 }
387 }
388
389 std::vector<double> clusterCommon2Error(icl + 1);
390 std::vector<int> clusterCode(icl + 1);
391 ncl = icl + 1;
392 for (int ic = 0; ic < icl + 1; ++ic) {
393 clusterX[ic] = 0.;
394 clusterY[ic] = 0.;
395 clusterZ[ic] = 0.;
396 clusterError[ic] = 0.;
397 clusterCommon2Error[ic] = 0.;
398 clusterHits[ic] = 0;
399 clusterCode[ic] = 0;
400 clusterSelect[ic] = 0;
401 double pullMax = 10.;
402 for (int i = 0; i < n; ++i) {
403 if (select[i] > 0) {
404 if (ic == clusterInt[i]) {
405 clusterSelect[ic] = select[i];
406 double w = 1. / (error[i] * error[i]);
407 clusterX[ic] += hitx[i] * w;
408 clusterY[ic] += hity[i] * w;
409 clusterZ[ic] += hitz[i] * w;
410 clusterError[ic] += w;
411 if (std::abs(pull[i]) < std::abs(pullMax)) {
412 pullMax = pull[i];
413 clusterId[ic] = id[i];
414 clusterCode[ic] = scode[i];
415 clusterSelect[ic] = select[i];
416 }
417 clusterHits[ic]++;
418 if (clusterHits[ic] == 1) clusterCommon2Error[ic] = 0.;
419 }
420 }
421 }
422 clusterX[ic] = clusterX[ic] / clusterError[ic];
423 clusterY[ic] = clusterY[ic] / clusterError[ic];
424 clusterZ[ic] = clusterZ[ic] / clusterError[ic];
425 // Don't assume improvement on errors due to clustering
426 clusterError[ic] = std::sqrt(clusterHits[ic]) / std::sqrt(clusterError[ic]);
427 {
428 ATH_MSG_DEBUG("cluster phi " << ic << " x " << clusterX[ic] << " y " << clusterY[ic] << " z "
429 << clusterZ[ic] << " error " << clusterError[ic] << " hits " << clusterHits[ic]
430 << " select " << clusterSelect[ic] << " Code " << clusterCode[ic]);
431 }
432 }
433}
434void
435MuonPhiHitSelector::fitRecPhi(const double pmom, const std::vector<Identifier>& phiId,
436 const std::vector<double>& phiHitx, const std::vector<double>& phiHity,
437 const std::vector<double>& phiHitz, const std::vector<double>& phiError,
438 std::vector<int>& quality, const int nphi, std::vector<double>& phiPull,
439 std::vector<int>& phiMult, std::vector<int>& phiSelect, double& chi2, double& r0,
440 double& phi, std::vector<double>& errorM, int& nfit) const
441{
442
443 //
444 // Use reconstructed hits to perform fit for phi
445 //
446
447 ATH_MSG_DEBUG("Start phi fit reconstruction");
448
449 chi2 = 0.;
450 r0 = 0.;
451 nfit = 0;
452 phi = 0.;
453
454 int ncsc = 0;
455 int ntgc = 0;
456 int nrpc = 0;
457
458
459 if (nphi == 0) return;
460
461 std::vector<double> error0(nphi);
462 std::vector<double> error(nphi);
463 std::vector<double> errorf(nphi);
464 std::vector<int> scode(nphi);
465 std::vector<int> srcode(nphi);
466 std::vector<int> phiSelectKeep(nphi);
467 std::map<int, int> clusters;
468 std::map<int, int> clustersr;
469 std::map<int, int> clusterspat;
470
471 for (int i = 0; i < nphi; ++i) {
472
473 Identifier idi = phiId[i];
474 int code = 0;
475 int rcode = 0;
476 if (m_idHelperSvc->isRpc(idi)) {
477 code = 1000000 * (m_idHelperSvc->rpcIdHelper().stationName(idi))
478 + 10000 * (m_idHelperSvc->rpcIdHelper().stationPhi(idi))
479 + 100 * ((m_idHelperSvc->rpcIdHelper().stationEta(idi)) + 10);
480 code = code + 2 * ((m_idHelperSvc->rpcIdHelper().doubletR(idi)) - 1)
481 + 16 * ((m_idHelperSvc->rpcIdHelper().gasGap(idi)) - 1);
482 rcode = 1000000 * (m_idHelperSvc->rpcIdHelper().stationName(idi))
483 + 10000 * (m_idHelperSvc->rpcIdHelper().stationPhi(idi))
484 + 0 * ((m_idHelperSvc->rpcIdHelper().stationEta(idi)) + 10);
485 rcode = rcode + 2 * ((m_idHelperSvc->rpcIdHelper().doubletR(idi)) - 1)
486 + 16 * ((m_idHelperSvc->rpcIdHelper().gasGap(idi)) - 1);
487 } else if (m_idHelperSvc->isTgc(idi)) {
488 code = 1000000 * (m_idHelperSvc->tgcIdHelper().stationName(idi))
489 + 10000 * (m_idHelperSvc->tgcIdHelper().stationPhi(idi))
490 + 100 * ((m_idHelperSvc->tgcIdHelper().stationEta(idi)) + 10);
491 code = code + m_idHelperSvc->tgcIdHelper().gasGap(idi);
492 rcode = 1000000 * (m_idHelperSvc->tgcIdHelper().stationName(idi))
493 + 10000 * (m_idHelperSvc->tgcIdHelper().stationPhi(idi))
494 + 0 * ((m_idHelperSvc->tgcIdHelper().stationEta(idi)) + 10);
495 rcode = rcode + m_idHelperSvc->tgcIdHelper().gasGap(idi);
496 } else if (m_idHelperSvc->isCsc(idi)) {
497 code = 1000000 * (m_idHelperSvc->cscIdHelper().stationName(idi))
498 + 10000 * (m_idHelperSvc->cscIdHelper().stationPhi(idi))
499 + 100 * ((m_idHelperSvc->cscIdHelper().stationEta(idi)) + 10);
500 code = code + m_idHelperSvc->cscIdHelper().wireLayer(idi);
501 rcode = 1000000 * (m_idHelperSvc->cscIdHelper().stationName(idi))
502 + 10000 * (m_idHelperSvc->cscIdHelper().stationPhi(idi))
503 + 0 * ((m_idHelperSvc->cscIdHelper().stationEta(idi)) + 10);
504 rcode = rcode + m_idHelperSvc->cscIdHelper().wireLayer(idi);
505 }
506
507 scode[i] = code;
508 srcode[i] = rcode;
509 int idet = 0;
510 if (m_idHelperSvc->isRpc(idi))
511 idet = 1;
512 else if (m_idHelperSvc->isTgc(idi))
513 idet = 2;
514 else if (m_idHelperSvc->isCsc(idi))
515 idet = 3;
516 phiSelect[i] = idet;
517 phiSelectKeep[i] = idet;
518 }
519 // Hits on segments
520 for (int i = 0; i < nphi; ++i) {
521 if (phiError[i] != 0 && quality[i] > 100) {
522 clusters[scode[i]]++;
523 clustersr[srcode[i]]++;
524 }
525 }
526 // Drop hits on patterns that are in same station and layer as segment hit
527 // Avoid adding again (solved) ambiguous hits
528
529 for (int i = 0; i < nphi; ++i) {
530 if (phiError[i] != 0 && quality[i] > 0 && quality[i] < 100) {
531 if (clustersr.count(srcode[i]) > 0) {
532 quality[i] = 0;
533 } else {
534 clusterspat[scode[i]]++;
535 }
536 }
537 }
538
539 // Assign errors according to multiplicities
540 if (msgLvl(MSG::DEBUG))
541 ATH_MSG_DEBUG("phi hits " << nphi << " segment clusters " << clusters.size() << " pattern clusters "
542 << clusterspat.size());
543
544 for (int i = 0; i < nphi; ++i) {
545 error0[i] = 0;
546 Identifier id = phiId[i];
547 phiMult[i] = 0;
548 if (phiError[i] != 0 && quality[i] > 0) {
549 int n = 0;
550 if (quality[i] > 100) {
551 n = clusters[scode[i]];
552 // Treat phi hits from segment with high multiplicity > 10 as lying on patterm
553 if (n > 10) quality[i] = 10;
554 } else if (quality[i] < 100) {
555 n = clusterspat[scode[i]];
556 // Drop phi hits patterns with high multiplicity
557 if (clusters.count(scode[i]) == 1 || n > 10) {
558 n = 0;
559 // drop phi hits from pattern if already segment hits in same layer
560 quality[i] = 0;
561 phiSelect[i] = 0;
562 phiSelectKeep[i] = 0;
563 continue;
564 }
565 }
566 phiMult[i] = n;
567 double fact = 1.;
568 if (m_idHelperSvc->isRpc(id))
569 fact = 1.2;
570 else if (m_idHelperSvc->isTgc(id))
571 n = 1;
572 else if (m_idHelperSvc->isCsc(id))
573 n = 1;
574
575 error0[i] = phiError[i] * std::sqrt(n) * fact;
576 error[i] = phiError[i] * std::sqrt(n) * fact;
577 double phiHit = std::atan2(phiHity[i], phiHitx[i]);
578 {
579 ATH_MSG_DEBUG(i << " Station " << int(scode[i] / 1000000) << " Hit x " << phiHitx[i] << " Hit y "
580 << phiHity[i] << " Hit z " << phiHitz[i] << " error " << phiError[i] << " phi Hit "
581 << phiHit);
582 ATH_MSG_DEBUG("station " << phiSelect[i]);
583 ATH_MSG_DEBUG("code " << scode[i] << " multiplicity " << n << " error " << error0[i] << " quality "
584 << quality[i]);
585 if (error0[i] < 1.) ATH_MSG_DEBUG("TOO small error ");
586 }
587 }
588 }
589
590 // Count layers hit
591
592 std::map<int, int> layersHit;
593 for (int i = 0; i < nphi; ++i) {
594 if (phiError[i] != 0 && quality[i] > 0) {
595 layersHit[srcode[i]]++;
596 }
597 }
598 int allLayerHits = layersHit.size();
599 int allLayerRecoHits = 0;
600 double pfit = 20000.;
601
602 for (int iqua = 0; iqua < 3; ++iqua) {
603
604 double quacut = 10;
605 if (iqua == 1)
606 quacut = 0;
607 else if (iqua == 2) {
608 quacut = 10;
609 pfit = pmom;
610 }
611
612 ATH_MSG_DEBUG("Quality loop " << iqua << " quality cut " << quacut);
613 int nsel = 0;
614 int nselseg = 0;
615 for (int i = 0; i < nphi; ++i) {
616
617 if (iqua == 1) phiSelect[i] = phiSelectKeep[i];
618
619 if (phiError[i] != 0 && quality[i] > quacut) {
620 nsel++;
621 if (quality[i] > 100) nselseg++;
622 if (quality[i] == 10 && iqua == 1) quality[i] = 11;
623 } else {
624 phiSelect[i] = 0;
625 }
626 if (msgLvl(MSG::DEBUG))
627 ATH_MSG_DEBUG("index i " << i << " phiSelect " << phiSelect[i] << " Quality " << quality[i] << " error "
628 << error[i]);
629 }
630
631 int imax = -1;
632 if (iqua == 1 && nselseg > 0) {
633 // Test and drop pattern Hits if far off
634 double errorScaleFactor = 25.;
635 std::vector<int> phiPatSelect(nphi, 0);
636 for (int i = 0; i < nphi; ++i) {
637 phiPatSelect[i] = 0;
638 if (phiSelect[i] > 0 && quality[i] > 0 && quality[i] < 100) {
639 phiPatSelect[i] = 1;
640 error[i] = errorScaleFactor * error[i];
641 }
642 if (msgLvl(MSG::DEBUG))
643 ATH_MSG_DEBUG("select " << phiSelect[i] << " quality " << quality[i] << " error " << error[i]);
644 }
645 ATH_MSG_DEBUG("performing outlier removal for pattern hits ");
646 fitPhiSL(pfit, phiId, phiHitx, phiHity, phiHitz, error, phiSelect, nphi, phiPull, imax, chi2, r0, phi,
647 errorM, false);
648 for (int i = 0; i < nphi; ++i) {
649 if (phiPatSelect[i] == 1) {
650 error[i] = error[i] / errorScaleFactor;
651 double rescaledPull = phiPull[i] * errorScaleFactor;
652 // 3 sigma cut
653 if (std::abs(rescaledPull) < 3.) {
654 phiSelect[i] = phiSelectKeep[i];
655 } else {
656 phiSelect[i] = 0;
657 phiSelectKeep[i] = 0;
658 if (msgLvl(MSG::DEBUG))
659 ATH_MSG_DEBUG("Drop Pattern Hits with Quality == 1 "
660 << i << " quality " << quality[i] << " Pull " << rescaledPull << " phiSelect "
661 << phiSelect[i]);
662 }
663 }
664 }
665 }
666
667 const double pfitc = pfit;
668 imax = -1;
669
670 if (iqua == 2) {
671 // low momentum fit with scaled error (factor 10) for dropped segment hits
672 std::vector<int> phiReSelect(nphi);
673 for (int i = 0; i < nphi; ++i) {
674 ATH_MSG_DEBUG("select " << phiSelect[i] << " quality " << quality[i]);
675 phiReSelect[i] = 0;
676 if (phiSelect[i] == 0 && quality[i] > 99) {
677 phiReSelect[i] = 1;
678 phiSelect[i] = phiSelectKeep[i];
679 error[i] = 10. * error[i];
680 }
681 }
682 fitPhiSL(pfitc, phiId, phiHitx, phiHity, phiHitz, error, phiSelect, nphi, phiPull, imax, chi2, r0, phi,
683 errorM, false);
684 for (int i = 0; i < nphi; ++i) {
685 if (phiReSelect[i] == 1) {
686 error[i] = error[i] / 10.;
687 // 10 sigma cut (error rescale = 10)
688 if (std::abs(phiPull[i]) < 1) {
689 phiSelect[i] = phiSelectKeep[i];
690 } else {
691 phiSelect[i] = 0;
692 }
693 if (msgLvl(MSG::DEBUG))
694 ATH_MSG_DEBUG("Low momentum Quality == 2 add hit nr "
695 << i << " quality " << quality[i] << " Pull " << phiPull[i] << " phiSelect "
696 << phiSelect[i]);
697 }
698 }
699 }
700 if (iqua == 1 && msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG("Quality loop ");
701 nsel = 0;
702 for (int i = 0; i < nphi; ++i) {
703 errorf[i] = error[i];
704 if (iqua == 1) phiSelect[i] = phiSelectKeep[i];
705
706 if (phiError[i] != 0 && quality[i] > quacut) {
707 nsel++;
708 if (quality[i] == 10 && iqua == 1) quality[i] = 11;
709 } else {
710 phiSelect[i] = 0;
711 }
712 }
713
714 ATH_MSG_DEBUG("Selected PHI hits in fit " << nsel << " iqua " << iqua);
715 if (nsel == 0) continue;
716
717 int niter = -1;
718 // do hit dropping in maximal 10 iterations by putting quality to 0
719
720 for (int iter = 0; iter < 100; ++iter) {
721 niter++;
722 double power = (iter - 10) / 20.;
723 if (power < 0.) power = 0.;
724 chi2 = 0.;
725 nfit = 0;
726 if (iter > 10) {
727 // Shower treatment inflate errors with multiplicity
728 for (int i = 0; i < nphi; ++i) {
729 errorf[i] = error[i] * std::pow(phiMult[i], power);
730 }
731 }
732 fitPhiSL(pfitc, phiId, phiHitx, phiHity, phiHitz, errorf, phiSelect, nphi, phiPull, imax, chi2, r0, phi,
733 errorM, false);
734
735 ncsc = 0;
736 ntgc = 0;
737 nrpc = 0;
738
739 // Count layers hit in Reconstruction
740
741 std::map<int, int> layersRecoHit;
742
743 for (int i = 0; i < nphi; ++i) {
744 Identifier id = phiId[i];
745 if (error[i] == 0 || quality[i] < quacut) phiSelect[i] = 0;
746 if (error[i] != 0 && quality[i] > quacut) {
747 layersRecoHit[srcode[i]]++;
748 {
749 if (m_idHelperSvc->isRpc(id))
750 nrpc++;
751 else if (m_idHelperSvc->isTgc(id))
752 ntgc++;
753 else if (m_idHelperSvc->isCsc(id))
754 ncsc++;
755 }
756 nfit++;
757 }
758 }
759 allLayerRecoHits = layersRecoHit.size();
760 double frac = allLayerRecoHits / (allLayerHits + 0.0001);
761
762 if (nfit == 1) break;
763
764 if (imax < 0 || imax > nphi) {
765 ATH_MSG_DEBUG("Fitphi imax " << imax);
766 break;
767 }
768
769 if (chi2 < 5 * (nfit + 1) || std::abs(phiPull[imax]) < 3.0) {
770
771 ATH_MSG_DEBUG("Final phi " << phi << " frac " << frac << " chi2 " << chi2);
772 break;
773 }
774
775 phiSelect[imax] = 0;
776
777 {
778 ATH_MSG_DEBUG("Start hit dropping " << imax << " pullmax " << phiPull[imax] << " phi " << phi
779 << " chi2 " << chi2);
780 }
781 }
782
783 {
784 ATH_MSG_DEBUG("Fit results phi " << phi << " chi2 " << chi2 << " ndof " << nfit);
785 ATH_MSG_DEBUG("Reco RPC " << nrpc << " TGC " << ntgc << " CSC " << ncsc);
786 }
787
788
789 int nacc = 0;
790 int nshowerdrop = 0;
791 for (int i = 0; i < nphi; ++i) {
792 double power = (niter - 10) / 20.;
793 if (power < 0.) power = 0.;
794 double pull = phiPull[i] * std::pow(phiMult[i], power);
795 if (niter > 10 && std::abs(pull) > 3.0 && phiSelect[i] > 0) {
796 phiSelect[i] = 0;
797 quality[i] = 0;
798 nshowerdrop++;
799 if (msgLvl(MSG::DEBUG))
800 ATH_MSG_DEBUG("Drop shower hit i " << i << " with pull " << pull << " iterations " << niter
801 << " power " << power);
802 }
803 if (phiSelect[i] != 0) nacc++;
804 }
805 if (msgLvl(MSG::DEBUG))
806 ATH_MSG_DEBUG("phi hits " << nphi << " selected for fit " << nfit << " iqua " << iqua << " iterations "
807 << niter << " accepted hits " << nacc << " nshower drop " << nshowerdrop);
808 }
809}
810
811void
812MuonPhiHitSelector::fitPhiSL(const double pmom, const std::vector<Identifier>& /*id*/, const std::vector<double>& hitx,
813 const std::vector<double>& hity, const std::vector<double>& hitz,
814 const std::vector<double>& error, std::vector<int>& select, const int n,
815 std::vector<double>& pull, int& imax, double& chi2, double& r0, double& phi,
816 std::vector<double>& errorM, bool fast) const
817{
818
819 // Perform straight line fit to hits: good hits have select > 0
820 // in the fit scattering centres are added for nfit-1 angles
821 // WITH beamspot constraint
822 // degrees of freedom = 2*nfit
823
824 // Fit is based on matrix inversions formulae
825
826 // Inputs pmom = estimate of momentum
827 // id = identifiers hits
828 // hitx hity hitz = position in space
829 // error = associated error (in x-y plane)
830 // select = > 0 for selected hits
831 // n = total number of hits
832
833
834 // Outputs pull = residual (hit -fit) /error
835 // imax = index for hit with maximum pull
836 // chi2 = total chi2
837 // r0 = perigee parameter of fit (0,0)
838 // phi = azimuthal angle of fit at perigee
839
840 double pest = pmom;
841 if (pest > 20000.) pest = 20000.;
842
843 r0 = 0.;
844 phi = 0.;
845 chi2 = 0.;
846 imax = 0;
847
848 // Calculate mean position
849 double xm = 0.;
850 double ym = 0.;
851 double dtot = 0.;
852 double em = 0.;
853 for (int i = 0; i < n; ++i) {
854 if (error[i] != 0 && select[i] > 0) {
855 double inver2 = 1. / (error[i] * error[i]);
856 xm += hitx[i] * inver2;
857 ym += hity[i] * inver2;
858 dtot += std::sqrt(hitx[i] * hitx[i] + hity[i] * hity[i] + hitz[i] * hitz[i]) * inver2;
859 em += inver2;
860 }
861 }
862
863 if (em == 0) return;
864
865 dtot = dtot / em;
866
867 // Beamspot error 10 mm for cosmics 10000
868
869 double ebs = 0.1;
870 if (m_cosmics) ebs = 10000.;
871
872 ATH_MSG_DEBUG("pmom " << pmom << " error beam " << ebs);
873 double ebs2 = ebs * ebs;
874 double invebs2 = 1. / ebs2;
875 double xmc = xm / (em + invebs2);
876 double ymc = ym / (em + invebs2);
877 xm = xm / em;
878 ym = ym / em;
879
880 // Constraint on beam spot
881
882 double len2 = xmc * xmc + ymc * ymc;
883 double xcc = len2 * xmc * invebs2;
884 double ycc = len2 * ymc * invebs2;
885
886 for (int i = 0; i < n; ++i) {
887 if (error[i] != 0 && select[i] > 0) {
888 double inver2 = 1. / (error[i] * error[i]);
889 double xdiff = hitx[i] - xmc;
890 double ydiff = hity[i] - ymc;
891 double xdiff2 = xdiff * xdiff;
892 double ydiff2 = ydiff * ydiff;
893 len2 = xdiff2 + ydiff2;
894 double sign = 1.;
895 // Non Cosmics assume IP at 0 0
896 if (xdiff * hitx[i] + ydiff * hity[i] < 0 && !m_cosmics) sign = -1;
897 // Cosmics assume down going
898 if (ydiff < 0 && m_cosmics) sign = -1;
899 xcc += len2 * sign * xdiff * inver2;
900 ycc += len2 * sign * ydiff * inver2;
901 }
902 }
903
904 if (em > 0) phi = std::atan2(ycc, xcc);
905 CxxUtils::sincos scphi(phi);
906
907 r0 = xmc * scphi.sn - ymc * scphi.cs;
908 double x0 = r0 * scphi.sn;
909 double y0 = -r0 * scphi.cs;
910
911 if (msgLvl(MSG::DEBUG))
912 ATH_MSG_DEBUG("Constraint r0 " << r0 << " xpos " << xmc << " ypos " << ymc << " phi " << phi);
913 // assume 0,0
914 std::vector<double> d(n);
915 std::vector<double> dist(n);
916 std::map<double, int> distanceSort;
917 double pullmax = 0.;
918 for (int i = 0; i < n; ++i) {
919 if (error[i] != 0 && select[i] > 0) {
920 double xdiff = hitx[i] - x0;
921 double ydiff = hity[i] - y0;
922 double xdiff2 = xdiff * xdiff;
923 double ydiff2 = ydiff * ydiff;
924 d[i] = std::sqrt(xdiff2 + ydiff2);
925 dist[i] = std::sqrt(xdiff2 + ydiff2 + hitz[i] * hitz[i]);
926 distanceSort[dist[i]] = i;
927 pull[i] = hitx[i] * scphi.sn - hity[i] * scphi.cs - r0;
928 if (std::abs(pull[i]) > std::abs(pullmax)) {
929 pullmax = pull[i];
930 imax = i;
931 }
932 }
933 }
934
935 if (fast) return;
936
937 std::map<double, int>::iterator it = distanceSort.begin();
938 std::map<double, int>::iterator it_end = distanceSort.end();
939
940 int nfit = 0;
941 std::vector<double> xf(2 * n);
942 std::vector<double> lf(2 * n);
943 std::vector<double> yf(2 * n);
944 std::vector<double> ef(2 * n);
945 std::vector<int> indexf(n);
946 //
947 // measurements yf error ef at distance xf (0:nfit)
948 // scattering centra angle zero yf error ef at distance xf(nfit+1..2nfit-1)
949 // beamspot at yf(2 nfit) = 0 error ebs2 at distance xf(2 nfit)
950
951 for (; it != it_end; ++it) {
952 int index = it->second;
953 xf[nfit] = d[index];
954 lf[nfit] = dist[index];
955 yf[nfit] = (hitx[index] - xmc) * scphi.sn - (hity[index] - ymc) * scphi.cs;
956 ef[nfit] = error[index];
957 indexf[nfit] = index;
958 nfit++;
959 }
960
961 // NB start at 1 to add scattering centra
962
963 double erang = 0.030 * 5000. / (pest + 1000.);
964 for (int i = 1; i < nfit; ++i) {
965 xf[nfit + i - 1] = xf[i - 1];
966 yf[nfit + i - 1] = 0.;
967 double scale = 1.;
968 if (select[i] == 1)
969 scale = 1.;
970 else if (select[i] == 3)
971 scale = 0.5;
972 else if (select[i] == 2)
973 scale = 2.5;
974 ef[nfit + i - 1] = scale * erang;
975 }
976 // Beamspot
977 yf[2 * nfit - 1] = 0.;
978 xf[2 * nfit - 1] = 0.;
979 ef[2 * nfit - 1] = ebs;
980
981 Amg::MatrixX v(nfit + 1, 1);
982 v.setIdentity();
983
984 if (msgLvl(MSG::DEBUG))
985 ATH_MSG_DEBUG("fitPhiSL "
986 << " nfit " << nfit);
987
988 for (int i = 0; i < nfit + 1; ++i) {
989 v(i, 0) = 0.;
990 for (int j = 0; j < nfit; ++j) {
991 double inver2 = 1. / (ef[j] * ef[j]);
992 if (i == 0)
993 v(i, 0) += yf[j] * inver2;
994 else if (i == 1)
995 v(i, 0) += yf[j] * xf[j] * inver2;
996 else if (i > 1 && j > i - 2) {
997 v(i, 0) += yf[j] * (lf[j] - lf[i - 2]) * inver2;
998 }
999 }
1000 }
1001
1002 // Track Model Matrix
1003
1004 Amg::MatrixX model(nfit + 1, 2 * nfit);
1005 model.setIdentity();
1006 // Measurements related to position and slope
1007
1008 for (int i = 0; i < nfit + 1; ++i) {
1009 for (int j = 0; j < nfit; ++j) {
1010 model(i, j) = 0.;
1011 if (i == 0)
1012 model(i, j) = 1.;
1013 else if (i == 1)
1014 model(i, j) = xf[j];
1015 // scattering angle
1016 else if (i > 1 && j > i - 2)
1017 model(i, j) = lf[j] - lf[i - 2];
1018 }
1019 }
1020
1021 // Constraints on scattering angles and beamspot
1022
1023 for (int i = 0; i < nfit + 1; ++i) {
1024 for (int j = nfit; j < 2 * nfit; ++j) {
1025 model(i, j) = 0.;
1026 // scattering angle
1027 if (i == j - nfit + 2) model(i, j) = 1.;
1028 // Beam spot
1029 if (i == 0 && j == 2 * nfit - 1) model(i, j) = 1.;
1030 }
1031 }
1032
1033 // Covariance Inverse of Track parameters
1034
1035 Amg::MatrixX covT(nfit + 1, nfit + 1);
1036 for (int i = 0; i < nfit + 1; ++i) {
1037 for (int j = 0; j < nfit + 1; ++j) {
1038 covT(i, j) = 0.;
1039 for (int k = 0; k < 2 * nfit; ++k) {
1040 double er2 = ef[k] * ef[k];
1041 covT(i, j) += model(i, k) * model(j, k) / er2;
1042 }
1043 }
1044 }
1045
1046 // Invert covariance matrix and replace it (should be CovT)
1047 Amg::MatrixX covTI = covT.inverse();
1048
1049 Amg::MatrixX t(nfit + 1, 1);
1050 // Solution for Track parameters
1051 t = covTI * v;
1052
1053 if (msgLvl(MSG::DEBUG) && std::abs(t(1, 0)) > 0.2) {
1054 ATH_MSG_DEBUG("Don't trust fit result " << t(1, 0) << " Keep Old result");
1055 }
1056 if (std::abs(t(1, 0)) > 0.2) return;
1057
1058 // calculate residuals and chi2
1059 std::vector<double> resi(2 * nfit);
1060 std::vector<double> pulli(2 * nfit);
1061 std::vector<double> errf(2 * nfit);
1062 std::vector<double> pullf(2 * nfit);
1063 std::vector<double> resiOut(2 * nfit);
1064 std::vector<double> pullOut(2 * nfit);
1065 pullmax = 0.;
1066 int jmax = 0;
1067
1068 errorM[0] = covTI(0, 0);
1069 errorM[1] = covTI(0, 1);
1070 errorM[2] = covTI(1, 1);
1071 errorM[3] = 0.;
1072 if (nfit > 2) {
1073 double invlt = 1. / (lf[nfit - 1] - lf[1]);
1074 for (int i = 1; i < nfit - 1; ++i) {
1075 double w = (lf[nfit - 1] - lf[i]) * invlt;
1076 errorM[3] += covTI(i + 1, i + 1) * w * w;
1077 }
1078 }
1079
1080 {
1081 if (nfit >= 3) {
1082 ATH_MSG_DEBUG("Error angle " << covTI(3, 3));
1083 } // covTI has dim nfit+1
1084 ATH_MSG_DEBUG("errorM[3] " << errorM[3]);
1085 }
1086
1087 for (int i = 0; i < 2 * nfit; ++i) {
1088
1089 // Calculate prediction at each measurement i
1090 // propagate error of track parameters to measurement i
1091 double error2 = 0.;
1092 double ypred = 0.;
1093 for (int j = 0; j < nfit + 1; ++j) {
1094 if (msgLvl(MSG::DEBUG) && i == 0) ATH_MSG_DEBUG("Parameter j " << j << " t(j,0) " << t(j, 0));
1095 if (msgLvl(MSG::DEBUG) && model(j, i) != 0) ATH_MSG_DEBUG("i " << i << " model ij " << model(j, i));
1096 ypred += model(j, i) * t(j, 0);
1097 for (int k = 0; k < nfit + 1; ++k) {
1098 error2 += model(j, i) * covTI(j, k) * model(k, i);
1099 }
1100 }
1101 double ef_i2 = ef[i] * ef[i];
1102 double inv_ef_i2 = 1. / ef_i2;
1103
1104 resi[i] = ypred - yf[i];
1105 pulli[i] = resi[i] / ef[i];
1106
1107 // errf propagated error and pullf
1108 errf[i] = std::sqrt(error2);
1109 pullf[i] = resi[i] / errf[i];
1110
1111 // calculate residual without hit and error without hit
1112 // Think of Kalmanm method to exclude hit and error
1113 double err2invOut = 1. / error2 - inv_ef_i2;
1114 if (err2invOut > 0) {
1115 resiOut[i] = (ypred / error2 - yf[i] * inv_ef_i2) / err2invOut - yf[i];
1116 pullOut[i] = resiOut[i] / std::sqrt(1. / err2invOut + ef_i2);
1117 }
1118
1119 if ((i < nfit) and (std::abs(pullOut[i]) > std::abs(pullmax)) ) {
1120 imax = indexf[i];
1121 jmax = i;
1122 pullmax = pullOut[i];
1123 }
1124 chi2 += resi[i] * resi[i] * inv_ef_i2;
1125
1126 if (i < nfit) {
1127 pull[indexf[i]] = pullOut[i];
1128 }
1129 if (msgLvl(MSG::DEBUG) && i < nfit)
1130 ATH_MSG_DEBUG("i " << i << " index " << indexf[i] << " det " << select[indexf[i]] << " ypred " << ypred
1131 << " mst " << yf[i] << " residual " << resi[i] << " error " << ef[i] << " dist "
1132 << dist[i] << " hitz " << hitz[i] << " Pull " << pulli[i] << " Pullf " << pullf[i]
1133 << " resi out " << resiOut[i] << " pull out " << pullOut[i]);
1134 if (msgLvl(MSG::DEBUG) && i > nfit)
1135 ATH_MSG_DEBUG("i " << i << " ypred " << ypred << " mst " << yf[i] << " residual " << resi[i] << " error "
1136 << ef[i]);
1137 }
1138 r0 = r0 + t(0, 0);
1139 phi = phi + t(1, 0);
1140
1141 ATH_MSG_DEBUG("delta phi " << t(1, 0));
1142 if (msgLvl(MSG::DEBUG) && std::abs(t(1, 0)) > 0.1) ATH_MSG_DEBUG("ALARM delta phi " << t(1, 0));
1143
1144 if (msgLvl(MSG::DEBUG))
1145 ATH_MSG_DEBUG("Track parameters r0 " << r0 << " phi " << phi << " chi2 " << chi2 << " jmax " << jmax << " imax "
1146 << imax << " pullmax " << pullmax);
1147}
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition EventPrimitives.h:50
AmgMatrix
#define AmgMatrix(rows, cols)
Definition EventPrimitives.h:49
layersHit
static int layersHit(FPGATrackSimRoad &r)
Definition FPGATrackSimHough1DShiftTool.cxx:716
RIO_OnTrack.h
sign
int sign(int a)
Definition TRT_StrawNeighbourSvc.h:108
imax
int imax(int i, int j)
Definition TileLaserTimingTool.cxx:33
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition AthCommonMsg.h:30
MuonPhiHitSelector::m_competingRIOsOnTrackTool
ToolHandle< Muon::IMuonCompetingClustersOnTrackCreator > m_competingRIOsOnTrackTool
Toolhandle to CompetingRIOsOnTrackTool creator.
Definition MuonPhiHitSelector.h:44
MuonPhiHitSelector::m_makeClusters
Gaudi::Property< bool > m_makeClusters
flag that performs a clusterization and return clusters (default: false)
Definition MuonPhiHitSelector.h:59
MuonPhiHitSelector::m_cscRotCreator
ToolHandle< Muon::IMuonClusterOnTrackCreator > m_cscRotCreator
Definition MuonPhiHitSelector.h:52
MuonPhiHitSelector::MuonPhiHitSelector
MuonPhiHitSelector(const std::string &, const std::string &, const IInterface *)
Definition MuonPhiHitSelector.cxx:29
MuonPhiHitSelector::m_clusterCreator
ToolHandle< Muon::IMuonClusterOnTrackCreator > m_clusterCreator
Toolhandle to ClusterOnTrackTool creator.
Definition MuonPhiHitSelector.h:50
MuonPhiHitSelector::fitPhiSL
void fitPhiSL(const double pmom, const std::vector< Identifier > &id, const std::vector< double > &hitx, const std::vector< double > &hity, const std::vector< double > &hitz, const std::vector< double > &error, std::vector< int > &select, const int n, std::vector< double > &pull, int &imax, double &chi2, double &r0, double &phi, std::vector< double > &errorM, bool fast) const
fit method straight line model
Definition MuonPhiHitSelector.cxx:812
MuonPhiHitSelector::m_summary
Gaudi::Property< bool > m_summary
flag to print out a summary of what comes in and what comes out
Definition MuonPhiHitSelector.h:55
MuonPhiHitSelector::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition MuonPhiHitSelector.h:37
MuonPhiHitSelector::m_competingRios
Gaudi::Property< bool > m_competingRios
flag that build competing rios on track for amibguous trigger hits (default: false)
Definition MuonPhiHitSelector.h:61
MuonPhiHitSelector::initialize
StatusCode initialize() override
Definition MuonPhiHitSelector.cxx:38
MuonPhiHitSelector::clusterPhi
void clusterPhi(const std::vector< Identifier > &id, const std::vector< double > &hitx, const std::vector< double > &hity, const std::vector< double > &hitz, const std::vector< double > &error, const std::vector< double > &pull, std::vector< int > &select, const int n, std::vector< double > &clusterX, std::vector< double > &clusterY, std::vector< double > &clusterZ, std::vector< double > &clusterError, std::vector< Identifier > &clusterId, std::vector< int > &clusterHits, std::vector< int > &clusterSelect, std::vector< int > &clusterInt, int &ncl) const
clusterization method
Definition MuonPhiHitSelector.cxx:302
MuonPhiHitSelector::fitRecPhi
void fitRecPhi(const double pmom, const std::vector< Identifier > &phiId, const std::vector< double > &phiHitx, const std::vector< double > &phiHity, const std::vector< double > &phiHitz, const std::vector< double > &phiError, std::vector< int > &quality, const int nphi, std::vector< double > &phiPull, std::vector< int > &phiMult, std::vector< int > &phiSelect, double &chi2, double &r0, double &phi, std::vector< double > &errorM, int &nfit) const
fit method curved track model
Definition MuonPhiHitSelector.cxx:435
MuonPhiHitSelector::m_cosmics
Gaudi::Property< bool > m_cosmics
flag for use of cosmics, straight line model will be used, no interaction point constraint
Definition MuonPhiHitSelector.h:57
MuonPhiHitSelector::select_rio
std::vector< std::unique_ptr< const Trk::MeasurementBase > > select_rio(const double pmom, const std::vector< const Trk::RIO_OnTrack * > &associatedHits, const std::vector< const Trk::PrepRawData * > &unassociatedHits) const override
Selects and builds a cleaned vector of RIO fits the associatedHits and build new RIOs,...
Definition MuonPhiHitSelector.cxx:50
Muon::CscPrepData
Class representing clusters from the CSC.
Definition CscPrepData.h:39
Muon::MuonCluster::globalPosition
virtual const Amg::Vector3D & globalPosition() const =0
Returns the global position of the measurement (calculated on the fly)
Muon::RpcPrepData
Class to represent RPC measurements.
Definition RpcPrepData.h:35
Muon::TgcPrepData
Class to represent TGC measurements.
Definition TgcPrepData.h:32
Trk::RIO_OnTrack
Class to handle RIO On Tracks ROT) for InDet and Muons, it inherits from the common MeasurementBase.
Definition RIO_OnTrack.h:70
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition EventPrimitives.h:27
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
Trk::locX
@ locX
Definition ParamDefs.h:37
index
Definition index.py:1
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition sincos.h:39
type
Generated on for ATLAS Offline Software by doxygen 1.14.0