ATLAS Offline Software
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 
29 MuonPhiHitSelector::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 
37 StatusCode
38 MuonPhiHitSelector::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 
49 std::vector<std::unique_ptr<const Trk::MeasurementBase>>
50 MuonPhiHitSelector::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 
301 void
302 MuonPhiHitSelector::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 }
434 void
435 MuonPhiHitSelector::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 
811 void
812 MuonPhiHitSelector::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 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
MuonCluster.h
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
MdtReadoutElement.h
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
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:29
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
Trk::locX
@ locX
Definition: ParamDefs.h:43
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
index
Definition: index.py:1
hist_file_dump.d
d
Definition: hist_file_dump.py:137
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
skel.it
it
Definition: skel.GENtoEVGEN.py:423
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
CompetingMuonClustersOnTrack.h
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
Trk::RIO_OnTrack
Definition: RIO_OnTrack.h:70
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
Muon::MuonCluster::globalPosition
virtual const Amg::Vector3D & globalPosition() const =0
Returns the global position of the measurement (calculated on the fly)
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
MdtDriftCircleOnTrack.h
PlotPulseshapeFromCool.np
np
Definition: PlotPulseshapeFromCool.py:64
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:52
PrepRawData.h
AmgMatrix
#define AmgMatrix(rows, cols)
Definition: EventPrimitives.h:51
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
TgcClusterOnTrack.h
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
RpcClusterOnTrack.h
MuonPhiHitSelector::m_cscRotCreator
ToolHandle< Muon::IMuonClusterOnTrackCreator > m_cscRotCreator
Definition: MuonPhiHitSelector.h:52
CxxUtils::sincos::cs
double cs
Definition: sincos.h:95
MuonPhiHitSelector::initialize
StatusCode initialize() override
Definition: MuonPhiHitSelector.cxx:38
MuonPhiHitSelector.h
lumiFormat.i
int i
Definition: lumiFormat.py:92
Muon::RpcPrepData
Class to represent RPC measurements.
Definition: RpcPrepData.h:35
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
beamspotman.n
n
Definition: beamspotman.py:731
Muon::CscPrepData
Class representing clusters from the CSC.
Definition: CscPrepData.h:39
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
checkxAOD.frac
frac
Definition: Tools/PyUtils/bin/checkxAOD.py:256
CscClusterOnTrack.h
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
CscReadoutElement.h
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition: TRT_PAI_physicsConstants.h:20
test_pyathena.parent
parent
Definition: test_pyathena.py:15
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:127
find_tgc_unfilled_channelids.ip
ip
Definition: find_tgc_unfilled_channelids.py:3
python.StandardJetMods.pull
pull
Definition: StandardJetMods.py:264
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
MuonPhiHitSelector::m_clusterCreator
ToolHandle< Muon::IMuonClusterOnTrackCreator > m_clusterCreator
Toolhandle to ClusterOnTrackTool creator.
Definition: MuonPhiHitSelector.h:50
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
MuonPhiHitSelector::m_competingRIOsOnTrackTool
ToolHandle< Muon::IMuonCompetingClustersOnTrackCreator > m_competingRIOsOnTrackTool
Toolhandle to CompetingRIOsOnTrackTool creator.
Definition: MuonPhiHitSelector.h:44
min
#define min(a, b)
Definition: cfImp.cxx:40
grepfile.ic
int ic
Definition: grepfile.py:33
Trk::PrepRawData
Definition: PrepRawData.h:62
pmontree.code
code
Definition: pmontree.py:443
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
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:191
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
RIO_OnTrack.h
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
CxxUtils::sincos::sn
double sn
Definition: sincos.h:92
fast
bool fast
Definition: TrigGlobEffCorrValidation.cxx:190
python.PyAthena.v
v
Definition: PyAthena.py:157
DeMoScan.index
string index
Definition: DeMoScan.py:362
correlationModel::model
model
Definition: AsgElectronEfficiencyCorrectionTool.cxx:46
PlaneSurface.h
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
RunTileMonitoring.clusters
clusters
Definition: RunTileMonitoring.py:133
Muon::TgcPrepData
Class to represent TGC measurements.
Definition: TgcPrepData.h:32
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:76
MuonPhiHitSelector::MuonPhiHitSelector
MuonPhiHitSelector(const std::string &, const std::string &, const IInterface *)
Definition: MuonPhiHitSelector.cxx:29
MuonPhiHitSelector::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: MuonPhiHitSelector.h:37
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
Muon::MuonCluster
Class representing clusters in the muon system.
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonPrepRawData/MuonPrepRawData/MuonCluster.h:37
TgcReadoutElement.h
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
MuonClusterOnTrack.h
AthAlgTool
Definition: AthAlgTool.h:26
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
get_generator_info.error
error
Definition: get_generator_info.py:40
DerivationFramework::ClustersInCone::select
void select(const xAOD::IParticle *particle, const float coneSize, const xAOD::CaloClusterContainer *clusters, std::vector< bool > &mask)
Definition: ClustersInCone.cxx:14
error
Definition: IImpactPoint3dEstimator.h:70
LheEventFiller_Common.ef
ef
Definition: SFGen_i/share/common/LheEventFiller_Common.py:7
MuonPhiHitSelector::m_makeClusters
Gaudi::Property< bool > m_makeClusters
flag that performs a clusterization and return clusters (default: false)
Definition: MuonPhiHitSelector.h:59
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
RpcReadoutElement.h
Generated on Sun Jun 30 2024 21:22:55 for ATLAS Offline Software by   doxygen 1.8.18