ATLAS Offline Software
CscSegmentUtilTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "CscSegmentUtilTool.h"
6 
7 #include <cmath>
8 #include <iostream>
9 
10 #include "CscClusterization/ICscClusterFitter.h"
11 #include "CscClusterization/ICscStripFitter.h"
12 #include "CscSegmentMakers/ICscSegmentFinder.h"
13 #include "EventPrimitives/EventPrimitivesHelpers.h"
14 #include "MuonPrepRawData/CscPrepData.h"
15 #include "MuonPrepRawData/CscPrepDataContainer.h"
16 #include "MuonPrepRawData/CscStripPrepData.h"
17 #include "MuonRIO_OnTrack/CscClusterOnTrack.h"
18 #include "MuonReadoutGeometry/CscReadoutElement.h"
19 #include "MuonSegment/MuonSegment.h"
20 #include "MuonSegment/MuonSegmentCombination.h"
21 #include "MuonSegment/MuonSegmentCombinationCollection.h"
22 #include "TMath.h" // for TMath::Landau()
23 #include "TrkEventPrimitives/FitQuality.h"
24 #include "TrkEventPrimitives/LocalDirection.h"
25 #include "TrkEventPrimitives/ParamDefs.h"
26 #include "TrkSurfaces/PlaneSurface.h"
27 #include "TrkSurfaces/RotatedTrapezoidBounds.h"
28 #include "TrkSurfaces/TrapezoidBounds.h"
29 
30 using Muon::CscClusterOnTrack;
31 using Muon::CscPrepData;
32 using Muon::MuonSegment;
33 using Muon::MuonSegmentCombination;
34 using MuonGM::CscReadoutElement;
35 using Trk::PlaneSurface;
36 using Trk::RIO_OnTrack;
37 
38 namespace {
39  std::string station_name(const int station) {
40  if (station == 1)
41  return "CSS";
42  else if (station == 2)
43  return "CSL";
44  else
45  return "UNKNOWN_STATION";
46  }
47 
48  std::string measphi_name(const bool measphi) {
49  if (measphi)
50  return "phi";
51  else
52  return "eta";
53  }
54 
55  double alignConst(const bool measphi, const int wlay) {
56  if (measphi) return 0.;
57  const double aConst[4] = {0, -0.2289, -0.620181, -0.6534445};
58  return aConst[wlay - 1];
59  }
60 
61  bool enoughHitLayers(const ICscSegmentFinder::ChamberTrkClusters& eta_clus, const ICscSegmentFinder::ChamberTrkClusters& phi_clus) {
62  int nHitLayer_eta = 0;
63  int nHitLayer_phi = 0;
64  for (int i = 0; i < 4; ++i) {
65  if (!eta_clus[i].empty()) ++nHitLayer_eta;
66  if (!phi_clus[i].empty()) ++nHitLayer_phi;
67  }
68  return nHitLayer_eta >= 2 || nHitLayer_phi >= 2;
69  }
70 
71  bool IsUnspoiled(const Muon::CscClusterStatus status) {
72  return status == Muon::CscStatusUnspoiled || status == Muon::CscStatusSplitUnspoiled;
73  }
74 } // namespace
75 
76 //******************************************************************************
77 // Constructor.
78 CscSegmentUtilTool::CscSegmentUtilTool(const std::string& type, const std::string& name, const IInterface* parent) :
79  AthAlgTool(type, name, parent) {
80  declareInterface<ICscSegmentUtilTool>(this);
81 }
82 
83 /*********************************/
84 StatusCode CscSegmentUtilTool::initialize() {
85  ATH_MSG_DEBUG("Initializing " << name());
86  ATH_MSG_DEBUG(" Max chi-square: " << m_max_chisquare);
87  ATH_MSG_DEBUG(" chi-square tight: " << m_max_chisquare_tight);
88  ATH_MSG_DEBUG(" chi-square loose: " << m_max_chisquare_loose);
89  ATH_MSG_DEBUG(" Max r:phi slope: " << m_max_slope_r << " : " << m_max_slope_phi);
90  ATH_MSG_DEBUG(" Max segments/chamber: " << m_max_seg_per_chamber);
91  ATH_MSG_DEBUG(" ROT tan(theta) tolerance: " << m_fitsegment_tantheta_tolerance);
92  ATH_MSG_DEBUG(" cluster_error_scaler " << m_cluster_error_scaler);
93  ATH_MSG_DEBUG(" ROT creator: " << m_rotCreator.typeAndName());
94 
95  if (m_TightenChi2) m_IPerror = 2.;
96  if (m_TightenChi2) ATH_MSG_DEBUG(" Chi2 cuts are tightened and m_IPerror is: " << m_IPerror);
97 
98  if (m_x5data) ATH_MSG_DEBUG(" Things for X5Data analysis is applied such as alignment ");
99 
100  ATH_CHECK(m_rotCreator.retrieve());
101 
102  ATH_CHECK(m_readKey.initialize());
103 
104  ATH_CHECK(m_DetectorManagerKey.initialize());
105 
106  ATH_CHECK(m_idHelperSvc.retrieve());
107 
108  ATH_CHECK(m_eventInfo.initialize());
109 
110  return StatusCode::SUCCESS;
111 }
112 
114 // Note: this code was required for old MuGirl but is not used by the current reconstruction
115 // Leaving it in place for backwards compatibility
116 std::unique_ptr<std::vector<std::unique_ptr<MuonSegment> > > CscSegmentUtilTool::getMuonSegments(
117  Identifier eta_id, Identifier phi_id, ICscSegmentFinder::ChamberTrkClusters& eta_clus, ICscSegmentFinder::ChamberTrkClusters& phi_clus,
118  const Amg::Vector3D& lpos000, const EventContext& ctx) const {
119  if (!enoughHitLayers(eta_clus, phi_clus)) {
120  ATH_MSG_DEBUG(" Could not find at least two individual layer hits! ");
121  std::unique_ptr<std::vector<std::unique_ptr<MuonSegment> > > segments;
122  return segments;
123  }
124 
125  for (int i = 0; i < 4; ++i)
126  ATH_MSG_DEBUG("getMuonSegments2: No of clusters in layer " << i << " " << eta_clus[i].size() << " " << phi_clus[i].size());
127 
128  ATH_MSG_DEBUG("getMuonSegments called get2dMuonSegmentCombination");
129  std::unique_ptr<MuonSegmentCombination> Muon2dSegComb(get2dMuonSegmentCombination(eta_id, phi_id, eta_clus, phi_clus, lpos000, ctx));
130 
131  ATH_MSG_DEBUG("getMuonSegments called get4dMuonSegmentCombination");
132  std::unique_ptr<MuonSegmentCombination> Muon4dSegComb(get4dMuonSegmentCombination(Muon2dSegComb.get(), ctx));
133 
134  std::unique_ptr<std::vector<std::unique_ptr<MuonSegment> > > segments_clone(new std::vector<std::unique_ptr<MuonSegment> >);
135 
136  if (Muon4dSegComb) {
137  std::vector<std::unique_ptr<MuonSegment> >* segments = Muon4dSegComb->stationSegments(0);
138 
139  if (!segments->empty()) {
140  for (unsigned int i = 0; i < segments->size(); ++i) { segments_clone->push_back(std::move(segments->at(i))); }
141  }
142  }
143 
144  return segments_clone;
145 }
146 
148 MuonSegmentCombination* CscSegmentUtilTool::get2dMuonSegmentCombination(Identifier eta_id, Identifier phi_id,
149  ICscSegmentFinder::ChamberTrkClusters& eta_clus,
150  ICscSegmentFinder::ChamberTrkClusters& phi_clus,
151  const Amg::Vector3D& lpos000, const EventContext& ctx, int etaStat,
152  int phiStat) const {
153  int nGoodEta = 0, nGoodPhi = 0;
154  for (int i = 0; i < 4; ++i) {
155  ATH_MSG_DEBUG("get2dMuonSegmentCombination2: No of clusters in layer " << i << " " << eta_clus[i].size() << " "
156  << phi_clus[i].size());
157  if ((etaStat % (int)std::pow(10, i + 1)) / (int)std::pow(10, i) == 0) nGoodEta++;
158  if ((phiStat % (int)std::pow(10, i + 1)) / (int)std::pow(10, i) == 0) nGoodPhi++;
159  }
160 
161  if (nGoodEta < 2 && nGoodPhi < 2) {
162  ATH_MSG_DEBUG(" Not enough hits in either eta or phi to form a segment");
163  MuonSegmentCombination* pcol = nullptr;
164  return pcol;
165  }
166 
167  MuonSegmentCombination* pcol = new MuonSegmentCombination;
168  pcol->setNGoodCscLayers(nGoodEta, nGoodPhi);
169  // Find 2D segments.
170  ICscSegmentFinder::Segments eta_segs;
171  ICscSegmentFinder::Segments phi_segs;
172 
173  // get2dSegments does : find_2dsegments -> find_2dseg3hit -> add_2dsegments
174  get2dSegments(eta_id, phi_id, eta_clus, phi_clus, eta_segs, phi_segs, lpos000, ctx, etaStat, phiStat);
175  std::unique_ptr<MuonSegmentCombination::SegmentVec> psegs(new MuonSegmentCombination::SegmentVec);
176  for (ICscSegmentFinder::Segments::const_iterator iseg = eta_segs.begin(); iseg != eta_segs.end(); ++iseg) {
177  std::unique_ptr<MuonSegment> pseg(build_segment(*iseg, false, eta_id, nGoodEta == 2, ctx)); // build_segment does getRios
178  if (pseg) {
179  ATH_MSG_DEBUG(" =============================> get2dMuonSegmentCombination:: MuonSegment time (eta) from build_segment is "
180  << pseg->time());
181  psegs->push_back(std::move(pseg));
182  }
183  }
184  ATH_MSG_DEBUG("added " << psegs->size() << " eta segments");
185  pcol->addSegments(std::move(psegs));
186 
187  // Insert phi-segments.
188  std::unique_ptr<MuonSegmentCombination::SegmentVec> phisegs(new MuonSegmentCombination::SegmentVec);
189  for (ICscSegmentFinder::Segments::const_iterator iseg = phi_segs.begin(); iseg != phi_segs.end(); ++iseg) {
190  std::unique_ptr<MuonSegment> pseg(build_segment(*iseg, true, phi_id, nGoodPhi == 2, ctx));
191  if (pseg) {
192  ATH_MSG_DEBUG(" get2dMuonSegmentCombination:: MuonSegment time (phi) from build_segment is " << pseg->time());
193  phisegs->push_back(std::move(pseg));
194  }
195  }
196  ATH_MSG_DEBUG("added " << phisegs->size() << " phi segments");
197  pcol->addSegments(std::move(phisegs));
198 
199  // Add to SG container.
200  ATH_MSG_DEBUG("Added " << eta_segs.size() << " r-segments and " << phi_segs.size() << " phi-segments.");
201 
202  return pcol;
203 }
204 
205 /**********************************************/
206 // Fit a segment using a list of clusters.
207 // Filling is least squares with the usual 1/d**2 weighting.
208 // local z = -38.51 -12.82 12.87 38.56
209 void CscSegmentUtilTool::fit_segment(const ICscSegmentFinder::TrkClusters& clus, const Amg::Vector3D& lpos000, double& s0, double& s1,
210  double& d0, double& d1, double& d01, double& chsq, double& time, double& dtime, double& zshift,
211  const EventContext& ctx, int outlierHitLayer) const {
212  ATH_MSG_DEBUG(" fit_segment called ");
213 
214  bool measphi = false;
215  fit_detailCalcPart1(clus, lpos000, s0, s1, d0, d1, d01, chsq, measphi, time, dtime, zshift, false, outlierHitLayer,
216  ctx); // IsSlopeGiven should be false
217 
218  if (measphi) return; // No need to do the second try
219  fit_detailCalcPart1(clus, lpos000, s0, s1, d0, d1, d01, chsq, measphi, time, dtime, zshift, true, outlierHitLayer,
220  ctx); }
221 
223 // if outlierid == iclu... error should be blown out in case of precision measurements..
224 // if IsSlopeGiven true, then.... error is recalculated.....
225 
226 void CscSegmentUtilTool::fit_detailCalcPart1(const ICscSegmentFinder::TrkClusters& clus, const Amg::Vector3D& lpos000, double& s0,
227  double& s1, double& d0, double& d1, double& d01, double& chsq, bool& measphi, double& time,
228  double& dtime, double& zshift, bool IsSlopeGiven, int outlierHitLayer,
229  const EventContext& ctx) const {
230  // if (IsSlopeGiven)
231  // measure zshift
232 
233  double q0 = 0.0;
234  double q1 = 0.0;
235  for (ICscSegmentFinder::TrkClusters::const_iterator iclu = clus.begin(); iclu != clus.end(); ++iclu) {
236  const ICscSegmentFinder::Cluster& cl = *iclu;
237  const CscClusterOnTrack* clu = cl.cl;
238  const CscPrepData* prd = clu->prepRawData();
239  Identifier id = clu->identify();
240 
241  measphi = m_idHelperSvc->cscIdHelper().measuresPhi(id);
242  // Cluster position.
243  double y = cl.locY();
244  if (m_x5data) { y = y - alignConst(measphi, m_idHelperSvc->cscIdHelper().wireLayer(id)); }
245  double x = cl.locX();
246 
247  bool isunspoiled = IsUnspoiled(prd->status());
248 
249  // Error in cluster position.
250  // Slope is given and Unspoiled Status..... then error should be recalculated depending on slope...
251  // Outlier treatment is only working in case of unspoiled...
252  double d = Amg::error(clu->localCovariance(), Trk::locX);
253  if (isunspoiled) {
254  if (IsSlopeGiven && outlierHitLayer != (iclu - clus.begin()))
255  d = m_rotCreator->GetICscClusterFitter()->getCorrectedError(prd, s1);
256  }
257 
258  d *= m_cluster_error_scaler; // This is for error scaler for cosmic!!!
259  double w = 1.0 / (d * d);
260  q0 += w;
261  q1 += w * x;
262  }
263  zshift = 0.0;
264  if (q0 > 0.) zshift = q1 / q0;
265  if (!m_zshift) zshift = 0.0;
266 
267  ATH_MSG_VERBOSE(" fit_detailCalcPart1 zshift " << zshift);
268 
269  q0 = 0.0;
270  q1 = 0.0;
271  double q2 = 0.0;
272  double q01 = 0.0;
273  double q11 = 0.0;
274  double q02 = 0.0;
275 
276  if (m_IPconstraint) {
277  double x = lpos000.z();
278  double y = measphi ? lpos000.x() : lpos000.y();
279 
280  ATH_MSG_DEBUG(" constraint x " << x << " y " << y << " error " << m_IPerror);
281 
282  double d = m_IPerror;
283  double w = 1.0 / (d * d);
284  q0 += w;
285  q1 += w * (x - zshift);
286  q2 += w * (x - zshift) * (x - zshift);
287  q01 += w * y;
288  q11 += w * (x - zshift) * y;
289  q02 += w * y * y;
290  }
291 
292  int cntSuccessHit = 0;
293  double sumHitTimes = 0.;
294  double sumHitTimeSquares = 0.;
295 
296  int cntLateHit = 0;
297  int cntEarlyHit = 0;
298  float latestEarlyTime = -9999;
299  float earliestLateTime = 9999;
300  for (ICscSegmentFinder::TrkClusters::const_iterator iclu = clus.begin(); iclu != clus.end(); ++iclu) {
301  const ICscSegmentFinder::Cluster& cl = *iclu;
302  const CscClusterOnTrack* clu = cl.cl;
303  const CscPrepData* prd = clu->prepRawData();
304  Identifier id = clu->identify();
305 
306  measphi = m_idHelperSvc->cscIdHelper().measuresPhi(id);
307  // Cluster position.
308  double y = cl.locY();
309  if (m_x5data) { y = y - alignConst(measphi, m_idHelperSvc->cscIdHelper().wireLayer(id)); }
310  double x = cl.locX();
311 
312  bool isunspoiled = IsUnspoiled(prd->status());
313 
314  // Error in cluster position.
315  // Slope is given and Unspoiled Status..... then error should be recalculated depending on slope...
316  // Outlier treatment is only working in case of unspoiled...
317  double d = Amg::error(clu->localCovariance(), Trk::locX);
318  if (isunspoiled) {
319  if (IsSlopeGiven) d = m_rotCreator->GetICscClusterFitter()->getCorrectedError(prd, s1);
320  if (outlierHitLayer == (iclu - clus.begin())) d = getDefaultError(id, measphi, prd, ctx);
321  }
322  d *= m_cluster_error_scaler; // This is for error scaler for cosmic!!!
323  ATH_MSG_VERBOSE(" +++fit_segment() x/y/d = " << x << " " << y << " " << d);
324 
325  double w = 1.0 / (d * d);
326  q0 += w;
327  q1 += w * (x - zshift);
328  q2 += w * (x - zshift) * (x - zshift);
329  q01 += w * y;
330  q11 += w * (x - zshift) * y;
331  q02 += w * y * y;
332 
333  if (prd->timeStatus() == Muon::CscTimeSuccess) {
334  ++cntSuccessHit;
335  sumHitTimes += prd->time();
336  sumHitTimeSquares += std::pow(prd->time(), 2);
337  }
338 
339  if (prd->timeStatus() == Muon::CscTimeEarly) {
340  ++cntEarlyHit;
341  if (latestEarlyTime < prd->time()) latestEarlyTime = prd->time();
342  }
343 
344  if (prd->timeStatus() == Muon::CscTimeLate) {
345  ++cntLateHit;
346  if (earliestLateTime > prd->time()) earliestLateTime = prd->time();
347  }
348  }
349 
350  ATH_MSG_DEBUG(" Time Calculation!!! " << cntSuccessHit << " " << cntEarlyHit << " " << cntLateHit << " " << sumHitTimes << " "
351  << latestEarlyTime << " " << earliestLateTime);
352 
353  if (cntSuccessHit > 0) {
354  time = sumHitTimes / float(cntSuccessHit);
355  float timeSquare = sumHitTimeSquares / float(cntSuccessHit);
356  dtime = timeSquare - time * time;
357  if (dtime < 0.0)
358  dtime = 0.0;
359  else
360  dtime = std::sqrt(dtime);
361  } else {
362  if (cntEarlyHit > 0 && cntLateHit > 0) {
363  time = 99999;
364  dtime = 99999;
365  ATH_MSG_DEBUG("Segment has nonzero earlyHits and nonzero lateHits. This should be backgrounds!!");
366  } else if (cntEarlyHit > 0) {
367  time = latestEarlyTime;
368  dtime = std::abs(latestEarlyTime);
369  } else if (cntLateHit > 0) {
370  time = earliestLateTime;
371  dtime = earliestLateTime;
372  } else {
373  time = 99999;
374  dtime = 99999;
375  ATH_MSG_DEBUG("Every Hit is CscTimeUnavailable or CscTimeStatusUndefined!!! Should be investigated");
376  }
377  }
378 
379  ATH_MSG_DEBUG(" calc result " << time << " " << dtime);
380 
381  // IsSlopeGiven (true) means 2nd try and somehow denominator is zero
382  // then return only available calc..
383  if (IsSlopeGiven && q2 * q0 - q1 * q1 == 0.0) return;
384 
385  fit_detailCalcPart2(q0, q1, q2, q01, q11, q02, s0, s1, d0, d1, d01, chsq);
386 }
387 
388 //******************************************************************************
389 //
390 int CscSegmentUtilTool::find_outlier_cluster(const ICscSegmentFinder::TrkClusters& clus, const Amg::Vector3D& lpos000,
391  double& returned_chsq, const EventContext& ctx) const {
392  int nunspoil = 0;
393 
394  double chsq_reference = 10000;
395  int remclu_ind = -1;
396  ATH_MSG_VERBOSE(" find_outlier_cluster clus size " << clus.size());
397 
398  // check all four cases of outlier ... found the best chsq!!!
399  for (unsigned int ire = 0; ire < clus.size(); ire++) {
400  ICscSegmentFinder::TrkClusters fitclus;
401  bool isunspoiled = IsUnspoiled(clus[ire].cl->status());
402  if (isunspoiled) ++nunspoil;
403 
404  for (unsigned int itr = 0; itr < clus.size(); itr++)
405  if (ire != itr) fitclus.push_back(clus[itr]);
406 
407  ATH_MSG_VERBOSE(" find_outlier_cluster drop ire " << ire << " fitclus size " << fitclus.size());
408 
409  double s0, s1, d0, d1, d01, chsq, time, dtime, zshift;
410  fit_segment(fitclus, lpos000, s0, s1, d0, d1, d01, chsq, time, dtime, zshift, ctx);
411  if (chsq < chsq_reference && isunspoiled) { // require outlier should be precision measurement
412  chsq_reference = chsq;
413  remclu_ind = ire;
414  }
415 
416  if (remclu_ind >= 0) ATH_MSG_VERBOSE(" new chsq " << ire << " " << chsq << " @" << remclu_ind);
417  }
418 
419  // If nunspoil clusters are only 2 or less, we don't want to treat this segment outlier removal process!!
420  // So, set remclu_ind = -1 meaning no outlier is found...
421  if (nunspoil < m_nunspoil) {
422  ATH_MSG_VERBOSE(" Number of unspoiled cluster is " << nunspoil << " which is not enough!! At least 3 unspoiled hits required!! ");
423  remclu_ind = -1;
424  }
425 
426  returned_chsq = 99999;
427  if (remclu_ind >= 0) {
428  ICscSegmentFinder::TrkClusters fitclus; // Muon::CscClusterOnTrack!!
429  // loading clusters into fitclus with outlier error blown up ....
430  for (unsigned int ire = 0; ire < clus.size(); ire++) fitclus.push_back(clus[ire]);
431 
432  ATH_MSG_VERBOSE(" final fit outlier_cluster drop cluster index " << remclu_ind << " fitclus.size " << fitclus.size());
433 
434  double s0, s1, d0, d1, d01, chsq, time, dtime, zshift;
435  fit_segment(fitclus, lpos000, s0, s1, d0, d1, d01, chsq, time, dtime, zshift, ctx, remclu_ind);
436  returned_chsq = chsq;
437  }
438 
439  return remclu_ind;
440 }
441 
442 //******************************************************************************
443 
444 // Fit a segment residual with a list of clusters.
445 // Segment is fit excluding cluster irclu and then the residual is calculated
446 // as the difference between the cluster position and that predicted by
447 // the segment. The error in the residual includes those of the segment and the
448 // excluded cluster.
449 
450 void CscSegmentUtilTool::fit_residual(const ICscSegmentFinder::TrkClusters& clus, const Amg::Vector3D& lpos000, unsigned int irclu,
451  double& res, double& dres, const EventContext& ctx) const {
452  ATH_MSG_DEBUG("CscSegmentUtilTool::fit_residual called ");
453 
454  ICscSegmentFinder::TrkClusters fitclus;
455  for (unsigned int iclu = 0; iclu < clus.size(); ++iclu) {
456  if (iclu != irclu) fitclus.push_back(clus[iclu]);
457  }
458  // Fit cluster.
459  double s0, s1, d0, d1, d01, chsq, time, dtime, zshift;
460  fit_segment(fitclus, lpos000, s0, s1, d0, d1, d01, chsq, time, dtime, zshift, ctx);
461 
462  // Extract excluded cluster paramters.
463  const CscClusterOnTrack* cot = clus[irclu].cl;
464 
465  Trk::ParamDefs ierr = Trk::loc1;
466  double y = clus[irclu].locY();
467  if (m_x5data) {
468  const Identifier& id = cot->identify();
469  y = y - alignConst(m_idHelperSvc->cscIdHelper().measuresPhi(id), m_idHelperSvc->cscIdHelper().wireLayer(id));
470  }
471  double x = clus[irclu].locX();
472  // Error in cluster position.
473  double d = Amg::error(cot->localCovariance(), ierr) * m_cluster_error_scaler;
474 
475  // This is for CscClusterCollection
476  if (d0 < 0 || d1 < 0) {
477  res = -99.;
478  dres = -9.;
479  } else {
480  // Calculate predicted position and error.
481  double seg_y = s0 + s1 * x;
482  double seg_dsquare = d0 * d0 + 2.0 * x * d01 + d1 * d1 * x * x;
483  // Calculate residual and error.
484  res = y - seg_y;
485  dres = d * d + seg_dsquare;
486 
487  ATH_MSG_VERBOSE(" fit_residual d0:d01:d1 :: seg_dsquare:dres " << d0 << " " << d01 << " " << d1 << " :: " << seg_dsquare << " "
488  << dres);
489 
490  if (dres < 0)
491  dres = -1.0 * std::sqrt(-1.0 * dres);
492  else if (dres >= 0)
493  dres = std::sqrt(dres);
494  else // in case of nan
495  dres = -9.;
496  }
497 }
498 
499 //******************************************************************************
500 
501 // Fit a segment using a list of RIO's.
502 // Filling is least squares with the usual 1/d**2 weighting.
503 
504 void CscSegmentUtilTool::fit_rio_segment(const Trk::PlaneSurface& ssrf, bool /*dump*/, const ICscSegmentFinder::RioList& rios, double& s0,
505  double& s1, double& d0, double& d1, double& d01, double& chsq, double& zshift) const {
506  ATH_MSG_DEBUG("CscSegmentUtilTool::fit_rio_segment called: ");
507 
508  // measure zshift
509 
510  double q0 = 0.0;
511  double q1 = 0.0;
512  for (ICscSegmentFinder::RioList::const_iterator irio = rios.begin(); irio != rios.end(); ++irio) {
513  const RIO_OnTrack& rio = **irio;
514  // Fetch the measurement.
515  Trk::ParamDefs iloc = Trk::loc1;
516  const Trk::LocalParameters& msmt = rio.localParameters();
517  double y = msmt[iloc];
518  if (m_x5data) {
519  y = y - alignConst(m_idHelperSvc->cscIdHelper().measuresPhi(rio.identify()),
520  m_idHelperSvc->cscIdHelper().wireLayer(rio.identify()));
521  }
522  // Fetch the measurement error.
523  const Amg::MatrixX& cov = rio.localCovariance();
524  int dim = cov.rows();
525  Trk::ParamDefs ierr = dim == 1 ? Trk::loc1 : iloc;
526  double d = Amg::error(cov, ierr);
527  // Fetch the position of the measurement coordinate.
528  const Amg::Vector3D& gpos = rio.globalPosition();
529  Amg::Vector3D lpos = ssrf.transform().inverse() * gpos;
530  double x = lpos.z();
531  if (msgLvl(MSG::VERBOSE)) {
532  ATH_MSG_VERBOSE(" RIO global pos: " << gpos.x() << " " << gpos.y() << " " << gpos.z());
533  ATH_MSG_VERBOSE(" RIO local pos: " << lpos.x() << " " << lpos.y() << " " << lpos.z());
534  }
535  ATH_MSG_DEBUG(" RIO: " << x << " " << y << " " << d);
536 
537  // Update least-square sums.
538  double w = 1.0 / (d * d);
539  q0 += w;
540  q1 += w * x;
541  }
542 
543  zshift = 0.0;
544  if (q0 > 0) zshift = q1 / q0;
545  if (!m_zshift) zshift = 0.0;
546  ATH_MSG_VERBOSE(" fit_rio_segment: zshift " << zshift);
547  // Initialize least-square sums.
548 
549  q0 = 0.;
550  q1 = 0.;
551  double q2 = 0.0;
552  double q01 = 0.0;
553  double q11 = 0.0;
554  double q02 = 0.0;
555  for (ICscSegmentFinder::RioList::const_iterator irio = rios.begin(); irio != rios.end(); ++irio) {
556  const RIO_OnTrack& rio = **irio;
557  // Fetch the measurement.
558  Trk::ParamDefs iloc = Trk::loc1;
559  const Trk::LocalParameters& msmt = rio.localParameters();
560  double y = msmt[iloc];
561  if (m_x5data) {
562  y = y - alignConst(m_idHelperSvc->cscIdHelper().measuresPhi(rio.identify()),
563  m_idHelperSvc->cscIdHelper().wireLayer(rio.identify()));
564  }
565  // Fetch the measurement error.
566  // Fetch the measurement error.
567  const Amg::MatrixX& cov = rio.localCovariance();
568  int dim = cov.rows();
569  Trk::ParamDefs ierr = dim == 1 ? Trk::loc1 : iloc;
570  double d = Amg::error(cov, ierr);
571  // Fetch the position of the measurement coordinate.
572  const Amg::Vector3D& gpos = rio.globalPosition();
573  Amg::Vector3D lpos = ssrf.transform().inverse() * gpos;
574  double x = lpos.z();
575  if (msgLvl(MSG::VERBOSE)) {
576  ATH_MSG_VERBOSE(" RIO global pos: " << gpos.x() << " " << gpos.y() << " " << gpos.z());
577  ATH_MSG_VERBOSE(" RIO local pos: " << lpos.x() << " " << lpos.y() << " " << lpos.z());
578  }
579  ATH_MSG_DEBUG(" RIO: " << x << " " << y << " " << d);
580 
581  // Update least-square sums.
582  double w = 1.0 / (d * d);
583  q0 += w;
584  q1 += w * (x - zshift);
585  q2 += w * (x - zshift) * (x - zshift);
586  q01 += w * y;
587  q11 += w * (x - zshift) * y;
588  q02 += w * y * y;
589  }
590 
591  // detailed calc part2
592  fit_detailCalcPart2(q0, q1, q2, q01, q11, q02, s0, s1, d0, d1, d01, chsq);
593 }
594 
595 //******************************************************************************
596 
597 // Fit a RIO segment residual with a list of clusters.
598 // Segment is fit excluding RIO irio and then the residual is calculated
599 // as the difference between the cluster position and that predicted by
600 // the segment. The error in the residual includes those of the segment and the
601 // excluded cluster.
602 
603 void CscSegmentUtilTool::fit_rio_residual(const Trk::PlaneSurface& ssrf, bool dump, const ICscSegmentFinder::RioList& rios,
604  unsigned int irrio, double& res, double& dres, double& rs, double& drs) const {
605  ATH_MSG_DEBUG("CscSegmentUtilTool::fit_rio_segment called ");
606 
607  ICscSegmentFinder::RioList fitrios;
608  ICscSegmentFinder::RioList fitrios3pt;
609  for (unsigned int irio = 0; irio < rios.size(); ++irio) {
610  if (irio != irrio) fitrios.push_back(rios[irio]);
611  // For three point method.... 0 x 2; 1 x 3;
612  if ((irrio == 1) && (irio == 0 || irio == 2)) fitrios3pt.push_back(rios[irio]);
613  if ((irrio == 2) && (irio == 1 || irio == 3)) fitrios3pt.push_back(rios[irio]);
614  }
615  // Fit cluster.
616  double s0, s1, d0, d1, d01, chsq, zshift;
617  fit_rio_segment(ssrf, dump, fitrios, s0, s1, d0, d1, d01, chsq, zshift);
618 
619  // Extract excluded cluster paramters (see fit_rio_segment).
620  const RIO_OnTrack& rio = *rios[irrio];
621  Trk::ParamDefs iloc = Trk::loc1;
622  const Trk::LocalParameters& msmt = rio.localParameters();
623  double y = msmt[iloc];
624  if (m_x5data) {
625  y = y -
626  alignConst(m_idHelperSvc->cscIdHelper().measuresPhi(rio.identify()), m_idHelperSvc->cscIdHelper().wireLayer(rio.identify()));
627  }
628  const Amg::MatrixX& cov = rio.localCovariance();
629  int dim = cov.rows();
630  Trk::ParamDefs ierr = dim == 1 ? Trk::loc1 : iloc;
631  double d = Amg::error(cov, ierr);
632  const Amg::Vector3D& gpos = rio.globalPosition();
633  Amg::Vector3D lpos = ssrf.transform().inverse() * gpos;
634  double x = lpos.z();
635  ATH_MSG_VERBOSE(" excluded RIO: " << x << " " << y << " " << d);
636 
637  // Calculate predicted position and error.
638  double seg_y = s0 + s1 * x;
639  double seg_dsquare = d0 * d0 + 2.0 * x * d01 + d1 * d1 * x * x;
640  // Calculate residual and error.
641  res = y - seg_y;
642  dres = d * d + seg_dsquare;
643 
644  ATH_MSG_VERBOSE(" fit_residual d0:d01:d1 :: seg_dsquare:dres " << d0 << " " << d01 << " " << d1 << " :: " << seg_dsquare << " "
645  << dres);
646 
647  if (dres < 0)
648  dres = -1.0 * std::sqrt(-1.0 * dres);
649  else if (dres >= 0)
650  dres = std::sqrt(dres);
651  else // in case of nan
652  dres = -9.;
653 
654  // 3pt method
655  // Calculate predicted position and error.
656  if (irrio == 1 || irrio == 2) {
657  // For 3 point method
658  double s0_3pt, s1_3pt, d0_3pt, d1_3pt, d01_3pt, chsq_3pt, zshift_3pt;
659  fit_rio_segment(ssrf, dump, fitrios3pt, s0_3pt, s1_3pt, d0_3pt, d1_3pt, d01_3pt, chsq_3pt, zshift_3pt);
660 
661  if (d0_3pt < 0 || d1_3pt < 0) {
662  rs = -99.;
663  drs = -9.;
664  } else {
665  seg_y = s0_3pt + s1_3pt * x;
666  seg_dsquare = d0_3pt * d0_3pt + 2.0 * x * d01_3pt + d1_3pt * d1_3pt * x * x;
667  // Calculate residual and error.
668  rs = y - seg_y;
669  drs = d * d + seg_dsquare;
670 
671  ATH_MSG_VERBOSE(" fit_rio_residual d0:d01:d1 :: seg_dsquare:dres " << d0_3pt << " " << d01_3pt << " " << d1_3pt
672  << " :: " << seg_dsquare << " " << drs);
673 
674  if (drs < 0)
675  drs = -1.0 * std::sqrt(-1.0 * drs);
676  else if (drs >= 0)
677  drs = std::sqrt(drs);
678  else // in case of nan
679  drs = -999.;
680  }
681  } else {
682  rs = -99.;
683  drs = -9.;
684  }
685 
686  ATH_MSG_DEBUG("Residual " << irrio << ": " << res << " " << dres);
687 }
688 
689 //******************************************************************************
690 
691 // Extract then numbers of spoiled and unspoiled measurements from a list
692 // of RIO's.
693 
694 void CscSegmentUtilTool::spoiled_count(const ICscSegmentFinder::RioList& rios, double threshold, int& nspoil, int& nunspoil) {
695  nspoil = 0;
696  nunspoil = 0;
697  for (ICscSegmentFinder::RioList::const_iterator irio = rios.begin(); irio != rios.end(); ++irio) {
698  const RIO_OnTrack& rio = **irio;
699  // Fetch the measurement.
700  Trk::ParamDefs iloc = Trk::loc1;
701  const Amg::MatrixX& cov = rio.localCovariance();
702  int dim = cov.rows();
703  Trk::ParamDefs ierr = dim == 1 ? Trk::loc1 : iloc;
704  double d = Amg::error(cov, ierr);
705  if (d < threshold)
706  ++nunspoil;
707  else
708  ++nspoil;
709  }
710 }
711 
712 void CscSegmentUtilTool::spoiled_count(const ICscSegmentFinder::RioList& rios, int& nspoil, int& nunspoil) {
713  nspoil = 0;
714  nunspoil = 0;
715  for (ICscSegmentFinder::RioList::const_iterator irio = rios.begin(); irio != rios.end(); ++irio) {
716  const RIO_OnTrack& rio = **irio;
717  const CscPrepData* pclu = dynamic_cast<const CscPrepData*>(rio.prepRawData());
718 
719  if (pclu) {
720  if (IsUnspoiled(pclu->status()))
721  ++nunspoil;
722  else
723  ++nspoil;
724  } else {
725  ATH_MSG_WARNING(" Could not find CscPrepData from RIO_OnTrack ");
726  }
727  }
728 }
729 
730 void CscSegmentUtilTool::spoiled_count(const ICscSegmentFinder::RioList& rios, int& nspoil, int& nunspoil, int& spoilmap) {
731  nspoil = 0;
732  nunspoil = 0;
733  spoilmap = 0;
734 
735  for (ICscSegmentFinder::RioList::const_iterator irio = rios.begin(); irio != rios.end(); ++irio) {
736  const RIO_OnTrack& rio = **irio;
737  const CscPrepData* pclu = dynamic_cast<const CscPrepData*>(rio.prepRawData());
738  if (pclu) {
739  if (IsUnspoiled(pclu->status()))
740  ++nunspoil;
741  else {
742  int wlay = m_idHelperSvc->cscIdHelper().wireLayer(pclu->identify());
743  ++nspoil;
744  spoilmap += int(std::pow(2, wlay - 1));
745  }
746  } else {
747  ATH_MSG_WARNING(" Could not find CscPrepData from RIO_OnTrack ");
748  }
749  }
750 }
751 
752 //******************************************************************************
753 // Build an ATLAS segment.
754 // Using the muon convention that x is normal to the plane of measurement:
755 // My r segment: (z, dz/dx)
756 // My phi segment: (y, dy/dx)
757 // ATLAS segment: (u, v, atan(dx/du), atan(dx/dv)) = (y, z, axu, axv)
758 // d(axy)/d(dx/dy) = -1/(1+(dx/dy)**2)
759 // where u is the measured coordinate.
760 // Use 0+/-1000 for the missing position and pi/2+/-1 for the missing direction.
761 
762 MuonSegment* CscSegmentUtilTool::build_segment(const ICscSegmentFinder::Segment& seg, bool measphi, Identifier chid, bool use2Lay,
763  const EventContext& ctx) const {
764  // chid from any last cluster in given chamber
765 
766  ATH_MSG_DEBUG("Building csc segment.");
767 
768  SG::ReadCondHandle<MuonGM::MuonDetectorManager> DetectorManagerHandle{m_DetectorManagerKey, ctx};
769  const MuonGM::MuonDetectorManager* MuonDetMgr{*DetectorManagerHandle};
770  if (MuonDetMgr == nullptr) {
771  ATH_MSG_ERROR("Null pointer to the read MuonDetectorManager conditions object");
772  return nullptr;
773  }
774  const CscReadoutElement* pro = MuonDetMgr->getCscReadoutElement(chid);
775  Amg::Transform3D gToLocal = pro->GlobalToAmdbLRSTransform();
776  Amg::Vector3D lpos000 = gToLocal * Amg::Vector3D(0.0, 0.0, 0.0);
777  Amg::Transform3D lToGlobal = gToLocal.inverse();
778 
779  // Surface.
780  // Position is the local origin transformed to the global coordinate system.
781  // We cannot use chamber position because it is at the first plane, not the
782  // origin of the local coordinate system!
783  ATH_MSG_VERBOSE(" build_segment: zshift " << seg.zshift);
784  Amg::Vector3D locp(0.0, 0.0, seg.zshift);
785  Amg::Vector3D glop = lToGlobal * locp;
786  // Use chamber rotation.
787  Amg::Transform3D xf;
788  xf = pro->transform(chid).rotation();
789  xf.pretranslate(glop);
790  // Use chamber bounds.
791  Trk::SurfaceBounds* pbnd = pro->bounds(chid).clone();
792  Trk::TrapezoidBounds* pbnd_trap = dynamic_cast<Trk::TrapezoidBounds*>(pbnd);
793  Trk::RotatedTrapezoidBounds* pbnd_rtrap = dynamic_cast<Trk::RotatedTrapezoidBounds*>(pbnd);
794  Trk::PlaneSurface* psrf = nullptr;
795 
796  ATH_MSG_VERBOSE(" Surface: ");
797  if (pbnd_trap) {
798  psrf = new Trk::PlaneSurface(xf, pbnd_trap);
799  ATH_MSG_VERBOSE("trapezoid");
800  } else if (pbnd_rtrap) {
801  psrf = new Trk::PlaneSurface(xf, pbnd_rtrap);
802  ATH_MSG_VERBOSE("rotated trapezoid");
803  } else {
804  ATH_MSG_FATAL(" Invalid boundary: " << *pbnd);
805  abort();
806  }
807  ATH_MSG_VERBOSE(" Segment surface: " << *psrf);
808 
809  Amg::Transform3D gToSurf = psrf->transform().inverse();
810  double lx = measphi ? seg.s0 : 0.;
811  double ly = measphi ? 0. : seg.s0;
812  double lz = seg.zshift;
813  double dlx = measphi ? seg.s1 : 0.;
814  double dly = measphi ? 0. : seg.s1;
815  double dlz = 1.;
816  Amg::Vector3D lposAMDB(lx, ly, lz);
817  Amg::Vector3D gpos = lToGlobal * lposAMDB;
818  Amg::Vector3D lposRef = gToSurf * gpos;
819  Amg::Vector3D ldirAMDB(dlx, dly, dlz);
820  Amg::Vector3D gdir = lToGlobal.linear() * ldirAMDB;
821  Amg::Vector3D ldirRef = gToSurf.linear() * gdir;
822  Amg::Vector3D lposRefShift = lposRef - ldirRef * lposRef.z();
823 
824  ATH_MSG_VERBOSE(" extrapolation to lposRef.z() " << lposRef.z());
825  double s0 = lposRefShift.x();
826  ATH_MSG_VERBOSE(" Input position, slope: " << s0 << " " << seg.s1);
827  ATH_MSG_VERBOSE(" Error: " << seg.d0 << " " << seg.d1 << " " << seg.d01);
828 
829  // Build list of RIO on track objects.
830  auto prios = ICscSegmentFinder::MbaseList{};
831  getRios(seg, &prios, measphi, ctx); // if hit is in outlier, error is estimated in width/sqrt(12)
832 
833  // Fit quality.
834  int ndof = int(prios.size()) - 2;
835  if (m_IPconstraint) ndof = ndof + 1;
836  if (use2Lay) ndof = 1;
837  Trk::FitQuality* pfq = new Trk::FitQuality(seg.chsq, ndof);
838  // Build segment.
839  // Build position vector.
840  Amg::Vector2D pos(s0, 0.0);
841  // Build direction vector.
842  double ameas = std::atan2(1.0, seg.s1);
843  Trk::LocalDirection pdir(ameas, M_PI_2);
844  // Error matrix.
845  double dfac = -1.0 / (1 + seg.s1 * seg.s1);
846  double e_pos_pos = seg.d0 * seg.d0;
847  double e_pos_dir = seg.d01 * dfac;
848  const double e_dir_dir = seg.d1 * seg.d1 * dfac * dfac;
849  Amg::MatrixX cov(4, 4);
850  cov.setIdentity();
851  cov(0, 0) = e_pos_pos;
852  cov(0, 2) = e_pos_dir;
853  cov(2, 0) = cov(0, 2);
854  cov(2, 2) = e_dir_dir;
855  cov(1, 1) = 1000000.0;
856  cov(1, 3) = 0.0;
857  cov(3, 1) = cov(1, 3);
858  cov(3, 3) = 1.0;
859  MuonSegment* pseg_ref = new MuonSegment(pos,
860  pdir,
861  std::move(cov),
862  psrf,
863  std::move(prios),
864  pfq,
865  Trk::Segment::Csc2dSegmentMaker);
866  pseg_ref->setT0Error(float(seg.time), float(seg.dtime));
867 
868  ATH_MSG_DEBUG(" build_segment:: right after ctor* " << pseg_ref->time());
869 
870  Amg::Transform3D globalToLocal = pro->transform(chid).inverse();
871  Amg::Vector3D d(globalToLocal.linear() * pseg_ref->globalDirection());
872 
873  double tantheta = 0;
874  if (d.z() == 0) {
875  ATH_MSG_WARNING("build_segment() - segment z is 0, set tantheta=0");
876  } else
877  tantheta = d.x() / d.z(); // is equal to seg.s1
878 
879  ATH_MSG_VERBOSE(" Position: " << pos[Trk::loc1] << " " << pos[Trk::loc2]);
880  ATH_MSG_VERBOSE("seg.s1 : ameas " << seg.s1 << " " << ameas);
881  ATH_MSG_VERBOSE(" Direction: " << pdir.angleXZ() << " " << pdir.angleYZ());
882  ATH_MSG_VERBOSE(" d.x : d.z : tantheta " << d.x() << " " << d.z() << " " << tantheta);
883  ATH_MSG_VERBOSE(resetiosflags(std::ios_base::floatfield));
884 
885  MuonSegment* pseg = pseg_ref->clone();
886  ATH_MSG_DEBUG(" build_segment:: right after copying * " << pseg->time());
887 
888  const double initialDiffTanTheta = 99;
889  double diff_tantheta = initialDiffTanTheta;
890  const unsigned int nTrials = 10;
891  unsigned int n_update = 0;
892  while (diff_tantheta > m_fitsegment_tantheta_tolerance && n_update < nTrials) {
893  // Loop over collections in the container.
894  auto prios_new = ICscSegmentFinder::MbaseList{};
895  ICscSegmentFinder::TrkClusters fitclus;
896  ICscSegmentFinder::RioList oldrios;
897  for (unsigned int irot = 0; irot < pseg_ref->numberOfContainedROTs(); irot++) oldrios.push_back(pseg_ref->rioOnTrack(irot));
898 
899  int cnt = 0;
900  for (ICscSegmentFinder::RioList::size_type irio = 0; irio < pseg_ref->numberOfContainedROTs(); ++irio) {
901  const Trk::RIO_OnTrack* pold = oldrios[irio];
902  const Trk::RIO_OnTrack* cot =
903  (seg.outlierid == cnt)
904  ? pold->clone() // No update for outlier owing to error blown up
905  : m_rotCreator->createRIO_OnTrack(*pold->prepRawData(), pold->globalPosition(), pseg->globalDirection());
906  if (!cot) {
907  ATH_MSG_WARNING("Failed to create CscClusterOnTrack");
908  continue;
909  }
910  const CscClusterOnTrack* pcl = dynamic_cast<const CscClusterOnTrack*>(cot);
911  if (!pcl) {
912  ATH_MSG_WARNING("Failed to cast to CscClusterOnTrack");
913  delete cot;
914  continue;
915  }
916  prios_new.push_back(cot);
917 
918  // Create new calibrated hit to put into fitclus
919  const CscPrepData* prd = pcl->prepRawData();
920  Amg::Vector3D lpos = gToLocal * pcl->globalPosition();
921 
922  ATH_MSG_DEBUG(" ---+++----> build_segment each rios time " << pcl->time() << " " << prd->time());
923 
924  ATH_MSG_VERBOSE(cnt << " " << n_update << " error ");
925  if (seg.outlierid == cnt)
926  ATH_MSG_VERBOSE(" !! outlier !! ");
927  else
928  ATH_MSG_VERBOSE(" !! !! ");
929  ATH_MSG_VERBOSE(Amg::error(pcl->localCovariance(), Trk::loc1));
930 
931  fitclus.emplace_back(lpos, pcl, measphi);
932  cnt++;
933  }
934 
935  ICscSegmentFinder::Segment seg_new; // Reconstruct segment in 2nd try with re-calibrated cluster position error
936  seg_new.outlierid = seg.outlierid;
937 
938  fit_segment(fitclus, lpos000, seg_new.s0, seg_new.s1, seg_new.d0, seg_new.d1, seg_new.d01, seg_new.chsq, seg_new.time,
939  seg_new.dtime, seg_new.zshift, ctx, seg.outlierid);
940  ATH_MSG_DEBUG("build_segments:: " << seg_new.time << " " << seg_new.dtime << " fitclus size " << fitclus.size());
941 
942  if (seg.outlierid >= 0) ATH_MSG_VERBOSE(n_update << " outlierid =" << seg.outlierid << " new chsq ==> " << seg_new.chsq);
943 
944  Trk::FitQuality* pfq_new = new Trk::FitQuality(seg_new.chsq, ndof);
945 
946  double lx = measphi ? seg_new.s0 : 0.;
947  double ly = measphi ? 0. : seg_new.s0;
948  double lz = seg_new.zshift;
949  double dlx = measphi ? seg_new.s1 : 0.;
950  double dly = measphi ? 0. : seg_new.s1;
951  double dlz = 1.;
952  Amg::Vector3D lposAMDB(lx, ly, lz);
953  Amg::Vector3D gpos = lToGlobal * lposAMDB;
954  Amg::Vector3D lposRef = gToSurf * gpos;
955  Amg::Vector3D ldirAMDB(dlx, dly, dlz);
956  Amg::Vector3D gdir = lToGlobal.linear() * ldirAMDB;
957  Amg::Vector3D ldirRef = gToSurf.linear() * gdir;
958  Amg::Vector3D lposRefShift = lposRef - ldirRef * lposRef.z();
959  double s0 = lposRefShift.x();
960 
961  Amg::Vector2D pos_new(s0, 0.0);
962  // Build direction vector.
963  double ameas_new = atan2(1.0, seg_new.s1);
964  Trk::LocalDirection pdir_new(ameas_new, M_PI_2);
965  // Error matrix.
966  double dfac_new = -1.0 / (1 + seg_new.s1 * seg_new.s1);
967  double e_pos_pos_new = seg_new.d0 * seg_new.d0;
968  double e_pos_dir_new = seg_new.d01 * dfac_new;
969  const double e_dir_dir_new = seg_new.d1 * seg_new.d1 * dfac_new * dfac_new;
970  Amg::MatrixX cov(4, 4);
971  cov.setIdentity();
972  cov(0, 0) = e_pos_pos_new;
973  cov(0, 2) = e_pos_dir_new;
974  cov(2, 2) = e_dir_dir_new;
975  cov(1, 1) = 1000000.0;
976  cov(1, 3) = 0.0;
977  cov(3, 3) = 1.0;
978 
979  MuonSegment* pseg_new =
980  new MuonSegment(pos_new,
981  pdir_new,
982  std::move(cov),
983  pseg->associatedSurface().clone(),
984  std::move(prios_new),
985  pfq_new,
986  Trk::Segment::Csc2dSegmentMaker);
987  pseg_new->setT0Error(float(seg_new.time), float(seg_new.dtime));
988  ATH_MSG_DEBUG(" build_segment:: right after recreating * " << pseg_new->time());
989 
990  Amg::Vector3D dnew(globalToLocal.linear() * pseg->globalDirection());
991  double tanthetanew = 0;
992  if (dnew.z() == 0) {
993  ATH_MSG_WARNING("build_segment() - new segment z is 0, staying with old segment and continue loop");
994  // we need to delete the new segment we just created
995  delete pseg_new;
996  // we set diff_tantheta again to its initial value and continue the loop
997  diff_tantheta = initialDiffTanTheta;
998  // if nTrials is reached anyway in the next iteration, break
999  if (n_update == (nTrials - 1)) break;
1000  continue;
1001  } else
1002  tanthetanew = dnew.x() / dnew.z();
1003 
1004  // this gets only called if the new segment did not have tanthetanew=0
1005  delete pseg;
1006  pseg = pseg_new;
1007  ATH_MSG_DEBUG(" build_segment:: right after new assigning and resetting * " << pseg->time());
1008 
1009  if (tantheta == 0) {
1010  ATH_MSG_WARNING("build_segment() - tantheta=0 but tanthetanew=" << tanthetanew << ", set diff_tantheta=tanthetanew");
1011  diff_tantheta = tanthetanew;
1012  } else {
1013  // both tanthetanew and tantheta are not 0
1014  diff_tantheta = std::abs(tanthetanew - tantheta) / tantheta;
1015  }
1016 
1017  ATH_MSG_VERBOSE(" tantheta change in segment: " << tantheta << " ==> " << tanthetanew << " -ratio- " << diff_tantheta);
1018 
1019  tantheta = tanthetanew;
1020  n_update++;
1021  }
1022  ATH_MSG_VERBOSE(" chisq " << (pseg->fitQuality())->chiSquared());
1023 
1024  if ((pseg->fitQuality())->chiSquared() > m_max_chisquare) {
1025  ATH_MSG_DEBUG(" Segment dropped chisq " << (pseg->fitQuality())->chiSquared() << " cut value " << m_max_chisquare);
1026  delete pseg;
1027  delete pseg_ref;
1028  return nullptr;
1029  }
1030 
1031  const Trk::LocalParameters& l = pseg->localParameters();
1032  ATH_MSG_VERBOSE(" Seg local key: " << l.parameterKey());
1033  ATH_MSG_VERBOSE(" Seg local params: " << l[Trk::loc1] << " " << l[Trk::loc2] << " " << l[Trk::phi] << " " << l[Trk::theta]);
1034  const Trk::LocalDirection& dr = pseg->localDirection();
1035  ATH_MSG_VERBOSE(" Seg local dir: " << dr.angleXZ() << " " << dr.angleYZ());
1036  const Amg::Vector3D& g = pseg->globalPosition();
1037  ATH_MSG_VERBOSE(" Seg global pos: " << g.x() << " " << g.y() << " " << g.z());
1038 
1039  delete pseg_ref;
1040 
1041  ATH_MSG_DEBUG(" build_segment:: right just before returning * " << pseg->time());
1042  return pseg;
1043 }
1044 
1045 /**************************************************************/
1046 void CscSegmentUtilTool::find_2dsegments(bool measphi, int station, int eta, int phi, const ICscSegmentFinder::ChamberTrkClusters& chclus,
1047  const Amg::Vector3D& lpos000, ICscSegmentFinder::Segments& segs, double lpos, double lslope,
1048  const EventContext& ctx) const {
1049  if (msgLvl(MSG::DEBUG)) {
1050  ATH_MSG_DEBUG("find_2dsegments called!! ID: " << measphi_name(measphi) << " " << std::showpos << eta << " "
1051  << station_name(station) << " " << phi << " ");
1052  ATH_MSG_DEBUG(" Counts: " << chclus[0].size() << " " << chclus[1].size() << " " << chclus[2].size() << " " << chclus[3].size());
1053  }
1054 
1055  const ICscSegmentFinder::TrkClusters& clus1 = chclus[0];
1056  const ICscSegmentFinder::TrkClusters& clus2 = chclus[1];
1057  const ICscSegmentFinder::TrkClusters& clus3 = chclus[2];
1058  const ICscSegmentFinder::TrkClusters& clus4 = chclus[3];
1059  ICscSegmentFinder::TrkClusters fitclus;
1060  fitclus.reserve(4);
1061  ICscSegmentFinder::Segment seg;
1062 
1063  // Loop over clusters in each layer and find all combinatoric segments.
1064  for (ICscSegmentFinder::TrkClusters::const_iterator icl1 = clus1.begin(); icl1 != clus1.end(); ++icl1) {
1065  for (ICscSegmentFinder::TrkClusters::const_iterator icl2 = clus2.begin(); icl2 != clus2.end(); ++icl2) {
1066  for (ICscSegmentFinder::TrkClusters::const_iterator icl3 = clus3.begin(); icl3 != clus3.end(); ++icl3) {
1067  for (ICscSegmentFinder::TrkClusters::const_iterator icl4 = clus4.begin(); icl4 != clus4.end(); ++icl4) {
1068  fitclus = {*icl1, *icl2, *icl3, *icl4};
1069 
1070  Trk::ParamDefs ierr = Trk::loc1;
1071  ATH_MSG_VERBOSE(" +++++++ find_2dsegments ++ Errors " << Amg::error((*icl1).cl->localCovariance(), ierr) << " "
1072  << Amg::error((*icl2).cl->localCovariance(), ierr) << " "
1073  << Amg::error((*icl3).cl->localCovariance(), ierr) << " "
1074  << Amg::error((*icl4).cl->localCovariance(), ierr) << " ");
1075 
1076  seg.s0 = lpos;
1077  seg.s1 = lslope;
1078  fit_segment(fitclus, lpos000, seg.s0, seg.s1, seg.d0, seg.d1, seg.d01, seg.chsq, seg.time, seg.dtime, seg.zshift, ctx);
1079  seg.clus[0] = *icl1;
1080  seg.clus[1] = *icl2;
1081  seg.clus[2] = *icl3;
1082  seg.clus[3] = *icl4;
1083 
1084  int nunspoil = 0;
1085  for (int i = 0; i < 4; ++i) {
1086  if (IsUnspoiled(seg.clus[i].cl->status())) ++nunspoil;
1087  ATH_MSG_DEBUG("Status of HIT #" << i << " " << seg.clus[i].cl->status());
1088  }
1089 
1090  seg.nunspoil = nunspoil;
1091  seg.outlierid = -1;
1092  seg.nclus = 4;
1093 
1094  double local_max_chi = 0.;
1095  if (nunspoil > 2)
1096  local_max_chi = m_max_chisquare_loose;
1097  else
1098  local_max_chi = m_max_chisquare;
1099 
1100  bool keep = true;
1101  if (seg.chsq > m_max_chisquare_tight) keep = false;
1102 
1103  ATH_MSG_VERBOSE(" find_2dsegments:: nunspoil and local_max_chi " << nunspoil << " " << local_max_chi << " " << keep
1104  << " " << m_max_chisquare_tight << " " << seg.chsq);
1105 
1106  if (!keep) {
1107  // chi2 after outlier removal on 4 hit segments.
1108  double outlierRemoved_chsq;
1109  int ioutlierid = find_outlier_cluster(seg.clus, lpos000, outlierRemoved_chsq, ctx);
1110  if (ioutlierid > 0 && outlierRemoved_chsq < m_max_chisquare_tight) {
1111  keep = true;
1112  seg.chsq = outlierRemoved_chsq;
1113  seg.outlierid = ioutlierid;
1114  ATH_MSG_DEBUG("find_2dsegments: keep outlier ");
1115  } else {
1116  ATH_MSG_DEBUG("find_2dsegments: reject segment bad chi2 ");
1117  }
1118  }
1119  ATH_MSG_DEBUG("find_2dsegments:: " << seg.time << " " << seg.dtime);
1120 
1121  // reject segments with slope of 0 (should only apply to NCB segments)
1122  if (seg.s1 == 0 && !m_IPconstraint) {
1123  ATH_MSG_DEBUG("slope too small, rejecting");
1124  keep = false;
1125  }
1126 
1127  if (keep) segs.push_back(seg);
1128 
1129  ATH_MSG_VERBOSE(" nunspoil and local_max_chi " << nunspoil << " " << local_max_chi << " " << keep << " " << seg.s1
1130  << " >< " << m_max_slope_r << " + " << seg.d1);
1131 
1132  ATH_MSG_VERBOSE(" Segment measphi? " << measphi << " chsq:" << seg.chsq << " abs(seg.s1) " << std::abs(seg.s1)
1133  << " req eta:phi " << m_max_slope_r + seg.d1 << " : " << m_max_slope_phi + seg.d1
1134  << " keep? " << keep);
1135  }
1136  }
1137  }
1138  }
1139 
1140  // Sort segments by nunspoil then chi-square.
1141  std::sort(segs.begin(), segs.end(), ICscSegmentFinder::sortByNunspoilAndChsq());
1142 }
1143 
1145 
1146 void CscSegmentUtilTool::add_2dsegments(ICscSegmentFinder::Segments& segs4, ICscSegmentFinder::Segments& segs3) const {
1147  // Limit number of segments
1148 
1149  if (segs4.empty() and segs3.empty()) return;
1150 
1151  ATH_MSG_DEBUG(" Total Input seg4 segment size " << segs4.size());
1152  ATH_MSG_DEBUG(" Total Input seg3 segment size " << segs3.size());
1153 
1154  std::vector<int> isegs4OK(segs4.size(), 1);
1155  std::vector<int> isegs3OK(segs3.size(), 1);
1156  ICscSegmentFinder::Segments segs4All{segs4};
1157 
1158  ICscSegmentFinder::Segments::const_iterator iseg;
1159  ICscSegmentFinder::Segments::const_iterator iseg2;
1160 
1161  segs4.clear();
1162 
1163 
1164 
1165  int iiseg = -1;
1166  for (iseg = segs4All.begin(); iseg != segs4All.end(); ++iseg) {
1167  iiseg++;
1168  int iiseg2 = iiseg;
1169  for (iseg2 = iseg + 1; iseg2 != segs4All.end(); ++iseg2) {
1170  int nhits_common = 0;
1171  iiseg2++;
1172  for (int iclus = 0; iclus < iseg->nclus; iclus++) {
1173  const Muon::CscClusterOnTrack* cot = iseg->clus[iclus].cl;
1174  for (int iclus2 = 0; iclus2 < iseg2->nclus; iclus2++) {
1175  const Muon::CscClusterOnTrack* cot2 = iseg2->clus[iclus2].cl;
1176  if (cot->identify() == cot2->identify()) nhits_common++;
1177  }
1178  }
1179  if (nhits_common > 0 && m_remove4Overlap) {
1180  isegs4OK[iiseg2] = 0;
1181  ATH_MSG_DEBUG(" seg4 segment nr " << iiseg2 << " dropped with nhits_common " << nhits_common);
1182  }
1183  }
1184  }
1185 
1186  iiseg = -1;
1187  for (iseg = segs4All.begin(); iseg != segs4All.end(); ++iseg) {
1188  iiseg++;
1189  const Muon::CscClusterOnTrack* cot = iseg->clus[0].cl;
1190  Identifier id = cot->identify();
1191  if (!m_idHelperSvc->cscIdHelper().measuresPhi(id) && iseg->nunspoil < m_nunspoil) {
1192  ATH_MSG_DEBUG(" seg4 eta segment rejected with nclusters too few unspoiled hits " << iseg->nclus << " chi2 " << iseg->chsq
1193  << " unspoiled " << iseg->nunspoil);
1194  isegs4OK[iiseg] = 0;
1195  }
1196  if (isegs4OK[iiseg] == 1 && segs4.size() < m_max_seg_per_chamber) {
1197  segs4.push_back(*iseg);
1198  ATH_MSG_DEBUG(" seg4 segment accepted with nclusters " << iseg->nclus << " chi2 " << iseg->chsq << " unspoiled "
1199  << iseg->nunspoil << " measuresPhi "
1200  << m_idHelperSvc->cscIdHelper().measuresPhi(id));
1201  } else {
1202  ATH_MSG_DEBUG(" seg4 segment rejected with nclusters " << iseg->nclus << " chi2 " << iseg->chsq << " unspoiled "
1203  << iseg->nunspoil << " measuresPhi "
1204  << m_idHelperSvc->cscIdHelper().measuresPhi(id));
1205  }
1206  }
1207  int segs4Size = segs4.size();
1208  ATH_MSG_DEBUG(" Total seg4 segment accepted size " << segs4Size);
1209 
1210  if (segs4.size() == m_max_seg_per_chamber) { return; }
1211 
1212  for (iseg = segs4.begin(); iseg != segs4.end(); ++iseg) {
1213  int iiseg2 = -1;
1214  for (iseg2 = segs3.begin(); iseg2 != segs3.end(); ++iseg2) {
1215  iiseg2++;
1216  int nhits_common = 0;
1217  for (int iclus = 0; iclus < iseg->nclus; iclus++) {
1218  const Muon::CscClusterOnTrack* cot = iseg->clus[iclus].cl;
1219  for (int iclus2 = 0; iclus2 < iseg2->nclus; iclus2++) {
1220  const Muon::CscClusterOnTrack* cot2 = iseg2->clus[iclus2].cl;
1221  if (cot->identify() == cot2->identify()) nhits_common++;
1222  }
1223  }
1224  if (nhits_common > 0 && m_remove3Overlap) {
1225  isegs3OK[iiseg2] = 0;
1226  ATH_MSG_DEBUG(" seg3 segment nr " << iiseg2 << " dropped with nhits_common " << nhits_common);
1227  }
1228  }
1229  }
1230 
1231  iiseg = -1;
1232  for (iseg = segs3.begin(); iseg != segs3.end(); ++iseg) {
1233  iiseg++;
1234  int iiseg2 = iiseg;
1235  for (iseg2 = iseg + 1; iseg2 != segs3.end(); ++iseg2) {
1236  iiseg2++;
1237  int nhits_common = 0;
1238  for (int iclus = 0; iclus < iseg->nclus; iclus++) {
1239  const Muon::CscClusterOnTrack* cot = iseg->clus[iclus].cl;
1240  for (int iclus2 = 0; iclus2 < iseg2->nclus; iclus2++) {
1241  const Muon::CscClusterOnTrack* cot2 = iseg2->clus[iclus2].cl;
1242  if (cot->identify() == cot2->identify()) nhits_common++;
1243  }
1244  }
1245  if (nhits_common > 0 && m_remove3Overlap) {
1246  isegs3OK[iiseg2] = 0;
1247  ATH_MSG_DEBUG(" seg3 segment nr " << iiseg2 << " dropped with nhits_common " << nhits_common);
1248  }
1249  }
1250  }
1251 
1252  // As long as we don't exceed max size, add segments from segs3 to end of segs4
1253  iiseg = -1;
1254  for (iseg = segs3.begin(); iseg != segs3.end(); ++iseg) {
1255  iiseg++;
1256  const Muon::CscClusterOnTrack* cot = iseg->clus[0].cl;
1257  Identifier id = cot->identify();
1258  if (!m_idHelperSvc->cscIdHelper().measuresPhi(id) && iseg->nunspoil < m_nunspoil) {
1259  ATH_MSG_DEBUG(" seg3 eta segment rejected with nclusters too few unspoiled hits " << iseg->nclus << " chi2 " << iseg->chsq
1260  << " unspoiled " << iseg->nunspoil);
1261  isegs3OK[iiseg] = 0;
1262  }
1263  if (isegs3OK[iiseg] == 1 && segs4.size() < m_max_seg_per_chamber) {
1264  segs4.push_back(*iseg);
1265  ATH_MSG_DEBUG(" seg3 segment accepted with nclusters " << iseg->nclus << " chi2 " << iseg->chsq << " unspoiled "
1266  << iseg->nunspoil << " measuresPhi "
1267  << m_idHelperSvc->cscIdHelper().measuresPhi(id));
1268  } else {
1269  ATH_MSG_DEBUG(" seg3 segment rejected with nclusters " << iseg->nclus << " chi2 " << iseg->chsq << " unspoiled "
1270  << iseg->nunspoil << " measuresPhi "
1271  << m_idHelperSvc->cscIdHelper().measuresPhi(id));
1272  }
1273  }
1274  ATH_MSG_DEBUG(" Total seg3 segment size " << segs4.size() - segs4Size);
1275 }
1276 
1278 // stores 2-hit segments
1279 void CscSegmentUtilTool::add_2dseg2hits(ICscSegmentFinder::Segments& segs, ICscSegmentFinder::Segments& segs2, int layStat) const {
1280  if (segs2.empty()) return;
1281  ATH_MSG_DEBUG(" Total Input 2-layer segment size " << segs2.size());
1282 
1283  int lay0 = -1, lay1 = -1;
1284  for (int i = 0; i < 4; i++) {
1285  if ((layStat % (int)std::pow(10, i + 1)) / (int)std::pow(10, i) == 0) {
1286  if (lay0 == -1)
1287  lay0 = i;
1288  else if (lay1 == -1)
1289  lay1 = i;
1290  }
1291  }
1292  bool checkCrossTalk = false;
1293  if (std::abs(lay0 - lay1) == 1 && lay0 + lay1 != 3)
1294  checkCrossTalk = true; // if we have layers 0 and 1 or 2 and 3 there could be cross-talk creating fake 2-layer segments
1295 
1296  std::vector<int> isegs2OK(segs2.size(), 1);
1297  ICscSegmentFinder::Segments::const_iterator iseg;
1298  ICscSegmentFinder::Segments::const_iterator iseg2;
1299  ICscSegmentFinder::Segments segsAll{segs};
1300 
1301  segs.clear();
1302 
1303  int iiseg = -1;
1304  for (iseg = segs2.begin(); iseg != segs2.end(); ++iseg) {
1305  iiseg++;
1306  if (!isegs2OK[iiseg]) continue;
1307  int iiseg2 = iiseg;
1308  for (iseg2 = iseg + 1; iseg2 != segs2.end(); ++iseg2) {
1309  int nhits_common = 0;
1310  iiseg2++;
1311  if (!isegs2OK[iiseg2]) continue;
1312  double charges[2] = {0, 0};
1313  for (int iclus = 0; iclus < iseg->nclus; iclus++) {
1314  const Muon::CscClusterOnTrack* cot = iseg->clus[iclus].cl;
1315  if (checkCrossTalk) {
1316  const Muon::CscPrepData* prep = cot->prepRawData();
1317  std::vector<const Muon::CscStripPrepData*> strips = m_rotCreator->GetICscClusterUtilTool()->getStrips(prep);
1318  std::vector<double> stripCharges;
1319  for (unsigned int s = 0; s < strips.size(); ++s) {
1320  ICscClusterFitter::StripFit sfit;
1321  sfit = m_rotCreator->GetICscStripFitter()->fit(*strips[s]);
1322  stripCharges.push_back(sfit.charge);
1323  }
1324  double maxCharge = 0, centCharge = 0;
1325  for (unsigned int s = 0; s < stripCharges.size(); s++) {
1326  if (stripCharges[s] > maxCharge) {
1327  maxCharge = stripCharges[s];
1328  centCharge = stripCharges[s];
1329  if (s > 0) centCharge += stripCharges[s - 1];
1330  if (s < stripCharges.size() - 1) centCharge += stripCharges[s + 1];
1331  }
1332  }
1333  charges[iclus] = centCharge;
1334  if (iclus == 1) {
1335  float chargeRatio = charges[0] / charges[1];
1336  if (charges[0] > charges[1]) chargeRatio = charges[1] / charges[0];
1337  if (chargeRatio < .01) {
1338  nhits_common = -1;
1339  break;
1340  } // ratio this small means crosstalk, kill this segment
1341  }
1342  }
1343  for (int iclus2 = 0; iclus2 < iseg2->nclus; iclus2++) {
1344  const Muon::CscClusterOnTrack* cot2 = iseg2->clus[iclus2].cl;
1345  if (cot->identify() == cot2->identify()) nhits_common++;
1346  }
1347  }
1348  if (nhits_common != 0) { //>0, overlap; <0, cross-talk
1349  isegs2OK[iiseg2] = 0;
1350  ATH_MSG_DEBUG(" seg2 segment nr " << iiseg2 << " dropped with nhits_common " << nhits_common);
1351  }
1352  }
1353  }
1354  iiseg = -1;
1355  for (iseg = segsAll.begin(); iseg != segsAll.end(); ++iseg) {
1356  iiseg++;
1357  int iiseg2 = -1;
1358  for (iseg2 = segs2.begin(); iseg2 != segs2.end(); ++iseg2) {
1359  int nhits_common = 0;
1360  iiseg2++;
1361  if (isegs2OK[iiseg2] == 0) continue; // already rejected this segment
1362  for (int iclus = 0; iclus < iseg->nclus; iclus++) {
1363  const Muon::CscClusterOnTrack* cot = iseg->clus[iclus].cl;
1364  int wlay = m_idHelperSvc->cscIdHelper().wireLayer(cot->identify());
1365  if ((layStat % (int)std::pow(10, wlay + 1)) / (int)std::pow(10, wlay) ==
1366  1) { // this 3-layer segment has a hit in a bad layer: dump it
1367  nhits_common = -1;
1368  break;
1369  }
1370  for (int iclus2 = 0; iclus2 < iseg2->nclus; iclus2++) {
1371  const Muon::CscClusterOnTrack* cot2 = iseg2->clus[iclus2].cl;
1372  if (cot->identify() == cot2->identify()) nhits_common++;
1373  }
1374  }
1375  if (nhits_common > 0) {
1376  isegs2OK[iiseg2] = 0;
1377  ATH_MSG_DEBUG(" seg2 segment nr " << iiseg2 << " dropped with nhits_common " << nhits_common);
1378  if (iseg2 == segs2.begin()) segs.push_back(*iseg); // no hits in bad layers, add this segment (but only once)
1379  } else if (nhits_common == 0) {
1380  if (iseg2 == segs2.begin()) segs.push_back(*iseg); // no hits in bad layers, add this segment (but only once)
1381  } else if (nhits_common == -1)
1382  break; // this segment has a hit in a bad layer, skip to the next one
1383  }
1384  }
1385  iiseg = -1;
1386  for (iseg = segs2.begin(); iseg != segs2.end(); ++iseg) {
1387  iiseg++;
1388  const Muon::CscClusterOnTrack* cot = iseg->clus[0].cl;
1389  Identifier id = cot->identify();
1390  if (isegs2OK[iiseg] == 1 && segs.size() < m_max_seg_per_chamber) {
1391  segs.push_back(*iseg);
1392  ATH_MSG_DEBUG(" seg2 accepted, nclusters " << iseg->nclus << " chi2 " << iseg->chsq << " unspoiled " << iseg->nunspoil
1393  << " mPhi " << m_idHelperSvc->cscIdHelper().measuresPhi(id));
1394  } else {
1395  ATH_MSG_DEBUG(" seg2 rejected, nclusters " << iseg->nclus << " chi2 " << iseg->chsq << " unspoiled " << iseg->nunspoil
1396  << " mPhi " << m_idHelperSvc->cscIdHelper().measuresPhi(id));
1397  }
1398  }
1399  ATH_MSG_DEBUG(" Total seg2 accepted: " << segs.size());
1400 }
1401 
1403 void CscSegmentUtilTool::find_2dseg3hit(bool measphi, int station, int eta, int phi, const ICscSegmentFinder::ChamberTrkClusters& chclus,
1404  const Amg::Vector3D& lpos000, ICscSegmentFinder::Segments& segs,
1405  ICscSegmentFinder::Segments& segs4hit, double lpos, double lslope, const EventContext& ctx) const {
1406  ATH_MSG_DEBUG("find_2dseg3hit called");
1407 
1408  // List of possible combinations for three hits.
1409  constexpr int maxcomb = 4;
1410  constexpr std::array<int, maxcomb> layAcomb{1, 2, 3, 0};
1411  constexpr std::array<int, maxcomb> layBcomb{2, 3, 0, 1};
1412  constexpr std::array<int, maxcomb> layCcomb{3, 0, 1, 2};
1413 
1414  ATH_MSG_VERBOSE("station " << station << " eta " << eta << " phi " << phi);
1415 
1416  // Maximum number of hits per segment
1417  constexpr int maxhits = 3;
1418  ICscSegmentFinder::TrkClusters fitclus;
1419  fitclus.reserve(3);
1420  ICscSegmentFinder::Segment seg;
1421 
1422  // Loop over four possible combinations of clusters for three segments.
1423  for (int icomb = 0; icomb < maxcomb; icomb++) {
1424  // iterator over clusters
1425  ICscSegmentFinder::TrkClusters::const_iterator icl[maxhits];
1426  auto & pCluster1 = icl[0];
1427  auto & pCluster2 = icl[1];
1428  auto & pCluster3 = icl[2];
1429  // Select appropriate layer for first layer
1430  const ICscSegmentFinder::TrkClusters& clus1 = chclus[layAcomb[icomb]];
1431  for (pCluster1 = clus1.begin(); pCluster1 != clus1.end(); ++pCluster1) {
1432  // Select appropriate layer for second layer
1433  const ICscSegmentFinder::TrkClusters& clus2 = chclus[layBcomb[icomb]];
1434  for (pCluster2 = clus2.begin(); pCluster2 != clus2.end(); ++pCluster2) {
1435  // Select appropriate layer for third layer
1436  const ICscSegmentFinder::TrkClusters& clus3 = chclus[layCcomb[icomb]];
1437  for (pCluster3 = clus3.begin(); pCluster3 != clus3.end(); ++pCluster3) {
1438  // Use these three clusters as a segment.
1439  fitclus = {*pCluster1, *pCluster2, *pCluster3};
1440 
1441  // Check if these hits are used by any other segments
1442  if (!unique_hits(fitclus, segs4hit)) {
1443  ATH_MSG_VERBOSE(" hits already used by other segments ");
1444  continue;
1445  }
1446  // Calculate chi2 for this segment.
1447  seg.s0 = lpos;
1448  seg.s1 = lslope;
1449  fit_segment(fitclus, lpos000, seg.s0, seg.s1, seg.d0, seg.d1, seg.d01, seg.chsq, seg.time, seg.dtime, seg.zshift, ctx);
1450 
1451  // Count number of unspoiled clusters
1452  int nunspoil = 0;
1453  for (int i = 0; i < maxhits; i++) {
1454  seg.clus[i] = *(icl[i]);
1455  if (IsUnspoiled(seg.clus[i].cl->status())) ++nunspoil;
1456  }
1457 
1458  seg.nunspoil = nunspoil;
1459  seg.outlierid = -1;
1460  seg.nclus = 3;
1461 
1462  // Cut on chi2
1463  double local_max_chi = 0.;
1464  if (nunspoil > 2)
1465  local_max_chi = m_max_chisquare_loose;
1466  else
1467  local_max_chi = m_max_chisquare;
1468 
1469  // tighten chi2 cut
1470  if (m_TightenChi2) local_max_chi = (2./3.) * m_max_chisquare ;
1471 
1472  bool keep = true;
1473  if (seg.chsq > local_max_chi) keep = false;
1474 
1475  ATH_MSG_VERBOSE(" nunspoil and local_max_chi " << nunspoil << " " << local_max_chi << " " << keep);
1476  ATH_MSG_DEBUG("find_2dseg3hit:: " << seg.time << " " << seg.dtime);
1477 
1478  if (!keep) {
1479  // chi2 after outlier removal on 3 hit segments.
1480  double outlierRemoved_chsq;
1481  int ioutlierid = find_outlier_cluster(fitclus, lpos000, outlierRemoved_chsq, ctx);
1482  if (ioutlierid > -1 && outlierRemoved_chsq < m_max_chisquare_tight) {
1483  keep = true;
1484  seg.chsq = outlierRemoved_chsq;
1485  seg.outlierid = ioutlierid;
1486  ATH_MSG_DEBUG("find_2dseg3hit: keep outlier ");
1487  } else {
1488  ATH_MSG_DEBUG("find_2dseg3hit: reject segment bad chi2");
1489  }
1490  }
1491  // Add to segment list.
1492  if (keep) segs.push_back(seg);
1493 
1494  ATH_MSG_VERBOSE(" Segment measphi? " << measphi << " chsq:" << seg.chsq << " abs(seg.s1) " << std::abs(seg.s1)
1495  << " keep? " << keep);
1496  }
1497  }
1498  }
1499  }
1500 
1501  // Sort segments by nunspoil then chi-square.
1502  std::sort(segs.begin(), segs.end(), ICscSegmentFinder::sortByNunspoilAndChsq());
1503 }
1504 void CscSegmentUtilTool::find_2dseg2hit(bool measphi, int station, int eta, int phi, int layStat,
1505  const ICscSegmentFinder::ChamberTrkClusters& chclus, const Amg::Vector3D& lpos000,
1506  ICscSegmentFinder::Segments& segs, double lpos, double lslope, const EventContext& ctx) const {
1507  ATH_MSG_DEBUG("find_2dseg2hit called");
1508  // List of possible combinations for three hits.
1509 
1510  ATH_MSG_VERBOSE("station " << station << " eta " << eta << " phi " << phi);
1511 
1512  int lay0 = -1, lay1 = -1;
1513  for (int i = 0; i < 4; i++) {
1514  if ((layStat % (int)std::pow(10, i + 1)) / (int)std::pow(10, i) == 0) {
1515  if (lay0 == -1)
1516  lay0 = i;
1517  else if (lay1 == -1)
1518  lay1 = i;
1519  else {
1520  ATH_MSG_WARNING("can't do 2-layer segment finding, more than 2 layers marked as good");
1521  return;
1522  }
1523  }
1524  }
1525  if (lay0 == -1 || lay1 == -1) {
1526  ATH_MSG_WARNING("can't do 2-layer segment finding, fewer than 2 layers marked as good");
1527  return;
1528  }
1529  // Maximum number of hits per segment
1530  const int maxhits = 2;
1531 
1532  // iterator over clusters
1533  ICscSegmentFinder::TrkClusters::const_iterator icl[maxhits] {};
1534 
1535  // Select appropriate layer for first layer
1536  const ICscSegmentFinder::TrkClusters& clus1 = chclus[lay0];
1537  for (icl[0] = clus1.begin(); icl[0] != clus1.end(); ++icl[0]) {
1538  // Select appropriate layer for second layer
1539  const ICscSegmentFinder::TrkClusters& clus2 = chclus[lay1];
1540  for (icl[1] = clus2.begin(); icl[1] != clus2.end(); ++icl[1]) {
1541  // Use these two clusters as a segment.
1542  ATH_MSG_DEBUG("got 2 clusters for segment");
1543 
1544  ICscSegmentFinder::TrkClusters fitclus;
1545  for (int i = 0; i < maxhits; i++) { fitclus.push_back(*(icl[i])); }
1546 
1547  // Calculate chi2 for this segment.
1548  ICscSegmentFinder::Segment seg;
1549  seg.s0 = lpos;
1550  seg.s1 = lslope;
1551  fit_segment(fitclus, lpos000, seg.s0, seg.s1, seg.d0, seg.d1, seg.d01, seg.chsq, seg.time, seg.dtime, seg.zshift, ctx);
1552 
1553  // Count number of unspoiled clusters
1554  int nunspoil = 0;
1555  for (int i = 0; i < maxhits; i++) {
1556  seg.clus[i] = *(icl[i]);
1557  if (IsUnspoiled(seg.clus[i].cl->status())) ++nunspoil;
1558  }
1559  seg.nunspoil = nunspoil;
1560  seg.outlierid = -1;
1561  seg.nclus = 2;
1562 
1563  // Cut on chi2
1564  double local_max_chi = 0.;
1565  if (nunspoil > 2)
1566  local_max_chi = m_max_chisquare_loose;
1567  else
1568  local_max_chi = m_max_chisquare;
1569  // tighten chi2 cut
1570  if (m_TightenChi2) local_max_chi = 1. * m_max_chisquare / 3.;
1571 
1572  bool keep = true;
1573  if (seg.chsq > local_max_chi) keep = false;
1574  ATH_MSG_VERBOSE(" nunspoil, chi2, and local_max_chi " << nunspoil << " " << seg.chsq << " " << local_max_chi << " " << keep);
1575  ATH_MSG_DEBUG("find_2dseg2hit:: " << seg.time << " " << seg.dtime);
1576  // No outlier done on 3 hit segments.
1577 
1578  ATH_MSG_VERBOSE(seg.s1 << " >< " << m_max_slope_r << " + " << seg.d1);
1579 
1580  // Add to segment list.
1581  if (keep) {
1582  ATH_MSG_DEBUG("good segment found");
1583  segs.push_back(seg);
1584  } else
1585  ATH_MSG_DEBUG("bad segment, not keeping");
1586 
1587  ATH_MSG_VERBOSE(" Segment measphi? " << measphi << " chsq:" << seg.chsq << " abs(seg.s1) " << std::abs(seg.s1)
1588  << " req eta:phi " << m_max_slope_r + seg.d1 << " : " << m_max_slope_phi + seg.d1
1589  << " keep? " << keep);
1590  }
1591  }
1592 
1593  // Sort segments by nunspoil then chi-square.
1594  std::sort(segs.begin(), segs.end(), ICscSegmentFinder::sortByNunspoilAndChsq());
1595 }
1596 
1597 void CscSegmentUtilTool::fit_detailCalcPart2(double q0, double q1, double q2, double q01, double q11, double q02, double& s0, double& s1,
1598  double& d0, double& d1, double& d01, double& chsq) const {
1599  ATH_MSG_VERBOSE(" q1 " << q1);
1600 
1601  if (q2 * q0 - q1 * q1 == 0.0) {
1602  s0 = -999;
1603  s1 = -999;
1604  d0 = -99;
1605  d1 = -99;
1606  d01 = -99;
1607  chsq = 999;
1608  return;
1609  }
1610 
1611  double r00 = q2 / (q2 * q0 - q1 * q1);
1612  double r01 = q1 / (q1 * q1 - q2 * q0);
1613  double r10 = q1 / (q1 * q1 - q0 * q2);
1614  double r11 = q0 / (q0 * q2 - q1 * q1);
1615  s0 = r01 * q11 + r00 * q01;
1616  s1 = r11 * q11 + r10 * q01;
1617  d0 = r01 * r01 * q2 + 2.0 * r00 * r01 * q1 + r00 * r00 * q0;
1618  if (d0 < 0)
1619  d0 = -1.0 * std::sqrt(-1.0 * d0);
1620  else
1621  d0 = std::sqrt(d0);
1622 
1623  d1 = r11 * r11 * q2 + 2.0 * r10 * r11 * q1 + r10 * r10 * q0;
1624  if (d1 < 0)
1625  d1 = -1.0 * std::sqrt(-1.0 * d1);
1626  else
1627  d1 = std::sqrt(d1);
1628  d01 = r01 * r11 * q2 + (r01 * r10 + r00 * r11) * q1 + r00 * r10 * q0;
1629  chsq = q02 + s1 * s1 * q2 + 2 * s0 * s1 * q1 + s0 * s0 * q0 - 2 * s0 * q01 - 2 * s1 * q11;
1630 
1631  ATH_MSG_VERBOSE(" details s0 = " << s0 << " " << d0 << " => " << r00);
1632  ATH_MSG_VERBOSE(" details s1 = " << s1 << " " << d1 << " " << d01 << " => " << r11 << " " << r10 << "chsq = " << chsq);
1633 }
1634 
1635 MuonSegmentCombination* CscSegmentUtilTool::get4dMuonSegmentCombination(const MuonSegmentCombination* insegs,
1636  const EventContext& ctx) const {
1637  ATH_MSG_DEBUG("get4dMuonSegmentCombination called");
1638  const ICscSegmentFinder::SegmentVec& rsegs = *insegs->stationSegments(0);
1639  const ICscSegmentFinder::SegmentVec& phisegs = *insegs->stationSegments(1);
1640 
1641  MuonSegmentCombination* pcol = nullptr;
1642 
1643  ATH_MSG_DEBUG(" Combination has r/phi segments " << rsegs.size() << " / " << phisegs.size());
1644  ATH_MSG_DEBUG("chamber has " << insegs->getNGoodCscLayers(1) << " good eta layers and " << insegs->getNGoodCscLayers(0)
1645  << " good phi layers");
1646 
1647  if ((rsegs.empty() && insegs->useStripsInSegment(1)) || (phisegs.empty() && insegs->useStripsInSegment(0))) {
1648  ATH_MSG_DEBUG(" Station does not have a r/phi segment, cannot make 4d segment ");
1649  return pcol;
1650  }
1651 
1652  SG::ReadHandle<xAOD::EventInfo> eventInfo(m_eventInfo, ctx);
1653 
1654  std::unique_ptr<ICscSegmentFinder::SegmentVec> pnewsegs(new ICscSegmentFinder::SegmentVec);
1655  if (insegs->useStripsInSegment(1) && insegs->useStripsInSegment(0)) {
1656  for (ICscSegmentFinder::SegmentVec::const_iterator irsg = rsegs.begin(); irsg != rsegs.end(); ++irsg) {
1657  for (ICscSegmentFinder::SegmentVec::const_iterator ipsg = phisegs.begin(); ipsg != phisegs.end(); ++ipsg) {
1658  // If there are too many segments, we break from the inner loop.
1659  // There will be one additional message for each remaining entry in
1660  // the outer loop.
1661  if (pnewsegs->size() >= m_max_seg_per_chamber) {
1662  ATH_MSG_DEBUG("Too many segments: Discarding segments m_max_seg_per_chamber=" << m_max_seg_per_chamber);
1663  break;
1664  }
1665  const MuonSegment& rsg = **irsg;
1666  const MuonSegment& psg = **ipsg;
1667  const Trk::FitQuality& rfq = *rsg.fitQuality();
1668  ATH_MSG_DEBUG("got fit quality for r segment");
1669  const Trk::FitQuality& phifq = *psg.fitQuality();
1670  ATH_MSG_DEBUG("and phi segment");
1671  // Fit quality.
1672  double chsq = rfq.chiSquared() + phifq.chiSquared();
1673  ATH_MSG_DEBUG("total chi2=" << chsq);
1674  if (chsq > m_max_chisquare) {
1675  ATH_MSG_DEBUG("Segment rejected too large chsq: " << chsq);
1676  continue;
1677  }
1678  // ECC - require good x-y matching.
1679  // - a better way would be to calculate this value for all segments,
1680  // sort the segments (best to worst) and keep only the top N best ones.
1681  double xylike = matchLikelihood(rsg, psg);
1682  ATH_MSG_DEBUG("xy likelihood: " << xylike << " min: " << m_min_xylike);
1683  if (xylike < m_min_xylike) {
1684  ATH_MSG_DEBUG("Segment rejected too low likelihood: " << xylike);
1685  continue;
1686  }
1687  // if don't use phi, make segments from eta only; if don't use eta, make segments from phi only; else make segments from
1688  // both
1689  std::unique_ptr<MuonSegment> pseg(make_4dMuonSegment(rsg, psg, insegs->use2LayerSegments(1), insegs->use2LayerSegments(0)));
1690  if (pseg) {
1691  ATH_MSG_DEBUG("created new 4d segment");
1692  pnewsegs->push_back(std::move(pseg));
1693  }
1694  } // for phisegs
1695  } // for rsegs
1696  } else if (!insegs->useStripsInSegment(0)) {
1697  for (ICscSegmentFinder::SegmentVec::const_iterator irsg = rsegs.begin(); irsg != rsegs.end(); ++irsg) {
1698  const MuonSegment& rsg = **irsg;
1699  unsigned int nMinRIOs = 3;
1700  if (insegs->use2LayerSegments(1)) nMinRIOs = 2;
1701  if (rsg.numberOfContainedROTs() < nMinRIOs) {
1702  ATH_MSG_DEBUG("only " << rsg.numberOfContainedROTs()
1703  << ", RIO's, insufficient to build the 4d segment from a single eta segment");
1704  continue;
1705  }
1706  std::unique_ptr<MuonSegment> pseg(new MuonSegment(rsg));
1707  ATH_MSG_DEBUG("created new 4d segment from eta hits only");
1708  pnewsegs->push_back(std::move(pseg));
1709  }
1710  } else if (!insegs->useStripsInSegment(1)) {
1711  for (ICscSegmentFinder::SegmentVec::const_iterator ipsg = phisegs.begin(); ipsg != phisegs.end(); ++ipsg) {
1712  const MuonSegment& psg = **ipsg;
1713  unsigned int nMinRIOs = 3;
1714  if (insegs->use2LayerSegments(0)) nMinRIOs = 2;
1715  if (psg.numberOfContainedROTs() < nMinRIOs) {
1716  ATH_MSG_DEBUG("only " << psg.numberOfContainedROTs()
1717  << ", RIO's, insufficient to build the 4d segment from a single phi segment");
1718  continue;
1719  }
1720  std::unique_ptr<MuonSegment> pseg(new MuonSegment(psg));
1721  ATH_MSG_DEBUG("created new 4d segment from phi hits only");
1722  pnewsegs->push_back(std::move(pseg));
1723  }
1724  }
1725 
1726  if (pnewsegs->empty()) {
1727  return pcol;
1728  } else
1729  ATH_MSG_DEBUG("created " << pnewsegs->size() << " new 4d segments");
1730 
1731  pcol = new MuonSegmentCombination;
1732  pcol->addSegments(std::move(pnewsegs));
1733 
1734  return pcol;
1735 }
1736 
1738 MuonSegmentCombination* CscSegmentUtilTool::get4dMuonSegmentCombination(Identifier eta_id, Identifier phi_id,
1739  ICscSegmentFinder::ChamberTrkClusters& eta_clus,
1740  ICscSegmentFinder::ChamberTrkClusters& phi_clus,
1741  const Amg::Vector3D& lpos000, const EventContext& ctx) const {
1742  ATH_MSG_DEBUG(" get4dMuonSegmentCombination called ");
1743 
1744  MuonSegmentCombination* Muon2dSegComb = get2dMuonSegmentCombination(eta_id, phi_id, eta_clus, phi_clus, lpos000, ctx);
1745 
1746  MuonSegmentCombination* Muon4dSegComb = get4dMuonSegmentCombination(Muon2dSegComb, ctx);
1747 
1748  delete Muon2dSegComb;
1749 
1750  return Muon4dSegComb;
1751 }
1752 
1753 /*************** PRIVATE FUNCTION *******************/
1754 MuonSegment* CscSegmentUtilTool::make_4dMuonSegment(const MuonSegment& rsg, const MuonSegment& psg, bool use2LaySegsEta,
1755  bool use2LaySegsPhi) const {
1756  ATH_MSG_DEBUG("make_4dMuonSegment called");
1757 
1758  double rpos = rsg.localParameters()[Trk::locX];
1759  double rdir = rsg.localDirection().angleXZ();
1760  double rtime = rsg.time();
1761  double rerrorTime = rsg.errorTime();
1762 
1763  const Amg::MatrixX& rcov = rsg.localCovariance();
1764  const Trk::PlaneSurface& etasrf = rsg.associatedSurface();
1765  ICscSegmentFinder::RioList etarios;
1766  for (unsigned int irot = 0; irot < rsg.numberOfContainedROTs(); irot++) etarios.push_back(rsg.rioOnTrack(irot));
1767  const Trk::FitQuality& rfq = *rsg.fitQuality();
1768 
1769  const Trk::FitQuality& phifq = *psg.fitQuality();
1770  double phipos = psg.localParameters()[Trk::locX];
1771  double phidir = psg.localDirection().angleXZ();
1772  const Amg::MatrixX& phicov = psg.localCovariance();
1773  const Trk::PlaneSurface& phisrf = psg.associatedSurface();
1774  ICscSegmentFinder::RioList phirios;
1775  for (unsigned int irot = 0; irot < psg.numberOfContainedROTs(); irot++) phirios.push_back(psg.rioOnTrack(irot));
1776 
1777  // Fit quality.
1778  double chsq = rfq.chiSquared() + phifq.chiSquared();
1779 
1780  Amg::MatrixX cov(4, 4);
1781  cov.setIdentity();
1782  // Local position and direction.
1783  Trk::LocalDirection pdir(phidir, rdir);
1784  // Error matrix.
1785  cov(0, 0) = phicov(0, 0);
1786  cov(2, 2) = phicov(2, 2);
1787  cov(0, 2) = phicov(0, 2);
1788  cov(2, 0) = cov(0, 2);
1789  cov(1, 1) = rcov(0, 0);
1790  cov(3, 3) = rcov(2, 2);
1791  cov(1, 3) = rcov(0, 2);
1792  cov(3, 1) = rcov(1, 3);
1793  ATH_MSG_VERBOSE("Error matrix is " << cov.cols() << "D");
1794  ATH_MSG_VERBOSE("( " << cov(0, 0) << " " << cov(0, 1) << " " << cov(0, 2) << " " << cov(0, 3));
1795  ATH_MSG_VERBOSE("( " << cov(1, 0) << " " << cov(1, 1) << " " << cov(1, 2) << " " << cov(1, 3));
1796  ATH_MSG_VERBOSE("( " << cov(2, 0) << " " << cov(2, 1) << " " << cov(2, 2) << " " << cov(2, 3));
1797  ATH_MSG_VERBOSE("( " << cov(3, 0) << " " << cov(3, 1) << " " << cov(3, 2) << " " << cov(3, 3));
1798 
1799  // Surface -- for rotations use phi surface for positions use eta surface
1800 
1801  Trk::SurfaceBounds* bnd = phisrf.bounds().clone();
1802  Trk::TrapezoidBounds* bnd_trap = dynamic_cast<Trk::TrapezoidBounds*>(bnd);
1803  Trk::RotatedTrapezoidBounds* bnd_rtrap = dynamic_cast<Trk::RotatedTrapezoidBounds*>(bnd);
1804  Amg::Vector3D glop = etasrf.transform().translation();
1805  Amg::Transform3D xf(phisrf.transform().rotation());
1806  xf.pretranslate(glop);
1807  Trk::PlaneSurface* psrf = nullptr;
1808  if (bnd_trap) {
1809  psrf = new Trk::PlaneSurface(xf, bnd_trap);
1810  } else if (bnd_rtrap) {
1811  psrf = new Trk::PlaneSurface(xf, bnd_rtrap);
1812  } else {
1813  ATH_MSG_WARNING(" no SurfaceBounds bounds for 4D segment found keep old bound ");
1814  psrf = new Trk::PlaneSurface(phisrf);
1815  }
1816 
1817  Amg::Vector3D philpos(phipos, 0., 0.);
1818  Amg::Vector3D phigpos = phisrf.transform() * philpos;
1819 
1820  Amg::Vector3D phietalpos = psrf->transform().inverse() * phigpos;
1821 
1822  ATH_MSG_VERBOSE(" positions in NEW Eta frame for phi measurement x " << phietalpos.x() << " y " << phietalpos.y() << " z shift "
1823  << phietalpos.z() << " angleXZ " << phidir);
1824  double phiposNew = phietalpos.x() - phidir * phietalpos.z();
1825 
1826  ATH_MSG_VERBOSE(" positions old z " << phipos << " New frame " << phietalpos.x() << " corrected " << phiposNew);
1827 
1828  Amg::Vector2D pos(phiposNew, rpos);
1829 
1830  // List of RIO's.
1831  auto rios = ICscSegmentFinder::MbaseList();
1832  // If each orientation has four RIOS's (as expected) then order the
1833  // RIO's eta1 phi1 eta2 .... to make fitting easier.
1834 
1835  // ECC - allow for 3-hit segments.
1836  unsigned int nMinRIOsEta = 3, nMinRIOsPhi = 3;
1837  if (use2LaySegsEta) nMinRIOsEta = 2;
1838  if (use2LaySegsPhi) nMinRIOsPhi = 2;
1839  if (etarios.size() >= nMinRIOsEta && phirios.size() >= nMinRIOsPhi) {
1840  ATH_MSG_DEBUG("Using new RIO order.");
1841  ATH_MSG_DEBUG(" eta/phi segment sizes: " << etarios.size() << " " << phirios.size());
1842  // ECC - try to match eta and phi layers
1843  // for ( RioList::size_type irio=0; irio<4; ++irio ) {
1844  int maxeta = etarios.size();
1845  int maxphi = phirios.size();
1846 
1847  // This is used to determine whether any hits are left unmatched.
1848  int eta_matchcode = 6;
1849  int phi_matchcode = 6;
1850  if (maxeta == 3) eta_matchcode = 3;
1851  if (maxphi == 3) phi_matchcode = 3;
1852  int eta_match = 0, phi_match = 0;
1853 
1854  ATH_MSG_DEBUG("make_4dMuonSegment:: r/phi seg time " << rsg.time() << " " << psg.time());
1855  // Loop over all hits in eta and then in phi
1856  for (int ieta = 0; ieta < maxeta; ieta++) {
1857  for (int iphi = 0; iphi < maxphi; iphi++) {
1858  // Get RIO_OnTrack
1859  const Trk::RIO_OnTrack* etapold = etarios[ieta];
1860  const Trk::RIO_OnTrack* phipold = phirios[iphi];
1861  // ECC figure out wire layer.
1862  Identifier id_eta = etapold->identify();
1863  Identifier id_phi = phipold->identify();
1864  int iw_eta = m_idHelperSvc->cscIdHelper().wireLayer(id_eta);
1865  int iw_phi = m_idHelperSvc->cscIdHelper().wireLayer(id_phi);
1866 
1867  // Check to see if these are the same layers.
1868  if (iw_eta != iw_phi) { continue; }
1869  ATH_MSG_DEBUG(" id_eta: " << m_idHelperSvc->toString(id_eta) << " "
1870  << " id_phi: " << m_idHelperSvc->toString(id_phi));
1871 
1872  // get the reference surface of the eta hit
1873  const Trk::Surface& surf = etapold->associatedSurface();
1874 
1875  // transform the global position of the phi hit into the reference frame of the eta hit
1876  std::optional<Amg::Vector2D> lpos = surf.globalToLocal(phipold->globalPosition());
1877 
1878  // calculate the 2D space point
1879  if (!lpos) {
1880  Amg::Vector3D locPos = surf.transform().inverse() * phipold->globalPosition();
1881  Amg::Vector3D locPosS = surf.transform().inverse() * phipold->associatedSurface().center();
1882  Amg::Vector3D locNorm = surf.transform().inverse().linear() * phipold->associatedSurface().normal();
1883  Amg::Vector2D lpn(locPos.x(), locPos.y());
1884  ATH_MSG_WARNING("phipold->globalPosition() not on surface!"
1885  << std::endl
1886  << std::setprecision(9) << " etapos: r " << etapold->globalPosition().perp() << " z "
1887  << etapold->globalPosition().z() << " surfz " << surf.center().z() << std::endl
1888  << surf << std::endl
1889  << " phipos: r " << phipold->globalPosition().perp() << " z " << phipold->globalPosition().z()
1890  << " surfz " << phipold->associatedSurface().center().z() << std::endl
1891  << phipold->associatedSurface() << std::endl
1892  << " locpos " << locPos.x() << " " << locPos.y() << " " << locPos.z() << " locposS " << locPosS.x()
1893  << " " << locPosS.y() << " " << locPosS.z() << " inbounds " << surf.insideBounds(lpn) << " normals "
1894  << std::setprecision(9) << surf.normal().dot(phipold->associatedSurface().normal()) << " locN "
1895  << locNorm.x() << " " << locNorm.y() << " " << locNorm.z());
1896  continue;
1897  }
1898 
1899  // Keep track of matched hits.
1900  phi_match++;
1901  eta_match++;
1902  eta_matchcode -= ieta;
1903  phi_matchcode -= iphi;
1904 
1905  Amg::Vector2D lposnew(etapold->localParameters()[Trk::locX], (*lpos)[Trk::locY]);
1906 
1907  // calculate the corresponding global position
1908  const Amg::Vector3D gposnew = surf.localToGlobal(lposnew);
1909  const Amg::Vector3D& gdirnew = rsg.globalDirection();
1910 
1911  // create the new rots using the ROT creator
1912  const Trk::RIO_OnTrack* etaRot = m_rotCreator->createRIO_OnTrack(*etapold->prepRawData(), gposnew, gdirnew);
1913  rios.push_back(etaRot);
1914 
1915  const Trk::RIO_OnTrack* phiRot = m_rotCreator->createRIO_OnTrack(*phipold->prepRawData(), gposnew);
1916  rios.push_back(phiRot);
1917 
1918  // debug output, to be removed
1919  const Trk::Surface& surfPhi = phipold->associatedSurface();
1920  std::optional<Amg::Vector2D> lposEta = surf.globalToLocal(etaRot->globalPosition());
1921  std::optional<Amg::Vector2D> lposPhi = surfPhi.globalToLocal(phiRot->globalPosition());
1922  ATH_MSG_VERBOSE(" eta gp " << etapold->globalPosition() << " new " << etaRot->globalPosition() << " loc " << *lposEta);
1923  ATH_MSG_VERBOSE(" phi gp " << phipold->globalPosition() << " new " << phiRot->globalPosition() << " loc " << *lposPhi);
1924 
1925  } // end loop over phi
1926  } // end loop over eta
1927 
1928  // Handle unmatched hits here.
1929  int eta_single = maxeta - eta_match;
1930  int phi_single = maxphi - phi_match;
1931 
1932  ATH_MSG_DEBUG(" eta_single, phi_single: " << eta_single << " " << phi_single << " eta_matchcode, phi_matchcode: " << eta_matchcode
1933  << " " << phi_matchcode);
1934 
1935  // Add unmatched eta hit to the segment.
1936  if (maxeta > 2) { // only if there are at least 3 hits, if it's a 2-layer segment require both hits to be matched
1937  if (eta_single == 1) {
1938  if (eta_matchcode >= 0 && eta_matchcode < 4) {
1939  const Trk::RIO_OnTrack* pold = etarios[eta_matchcode];
1940  Trk::RIO_OnTrack* pnew = pold->clone();
1941  rios.push_back(pnew);
1942  }
1943  }
1944  // This should never happen
1945  else if (eta_single > 1) {
1946  ATH_MSG_WARNING("More than one unmatched eta hit: " << eta_single);
1947  }
1948  } else {
1949  if (eta_single != 0) {
1950  ATH_MSG_DEBUG("eta hit in a 2-layer segment not matched, bailing");
1951  delete psrf;
1952  return nullptr;
1953  }
1954  }
1955 
1956  // Add unmatched phi hit to the segment.
1957  if (maxphi > 2) {
1958  if (phi_single == 1) {
1959  if (phi_matchcode >= 0 && phi_matchcode < 4) {
1960  const Trk::RIO_OnTrack* pold = phirios[phi_matchcode];
1961  Trk::RIO_OnTrack* pnew = pold->clone();
1962  rios.push_back(pnew);
1963  }
1964  }
1965  // This should never happen
1966  else if (phi_single > 1) {
1967  ATH_MSG_WARNING("More than one unmatched phi hit: " << phi_single);
1968  }
1969  } else {
1970  if (phi_single != 0) {
1971  ATH_MSG_DEBUG("phi hit in a 2-layer segment not matched, bailing");
1972  delete psrf;
1973  return nullptr;
1974  }
1975  }
1976  } else { // We should never get here!
1977  ATH_MSG_WARNING("Unexpected input RIO counts: " << etarios.size() << " " << phirios.size());
1978  for (ICscSegmentFinder::RioList::const_iterator irio = etarios.begin(); irio != etarios.end(); ++irio) {
1979  const Trk::RIO_OnTrack* pold = *irio;
1980  Trk::RIO_OnTrack* pnew = pold->clone();
1981  rios.push_back(pnew);
1982  }
1983  for (ICscSegmentFinder::RioList::const_iterator irio = phirios.begin(); irio != phirios.end(); ++irio) {
1984  const Trk::RIO_OnTrack* pold = *irio;
1985  Trk::RIO_OnTrack* pnew = pold->clone();
1986  rios.push_back(pnew);
1987  }
1988  }
1989 
1990  unsigned int nMinRIOsTot = 5;
1991  if (use2LaySegsEta || use2LaySegsPhi) nMinRIOsTot = 4;
1992  if (rios.size() < nMinRIOsTot) {
1993  ATH_MSG_WARNING("too few CSC hits collected, not making segment: rios " << rios.size());
1994  delete psrf;
1995  return nullptr;
1996  }
1997  int ndof = rfq.numberDoF() + phifq.numberDoF();
1998  Trk::FitQuality* pfq = new Trk::FitQuality(chsq, ndof);
1999  // Build 4D segment.
2000 
2001  ATH_MSG_DEBUG("Segment " << rios.size() << " : ");
2002  MuonSegment* pseg = new MuonSegment(pos,
2003  pdir,
2004  std::move(cov),
2005  psrf,
2006  std::move(rios),
2007  pfq,
2008  Trk::Segment::Csc4dSegmentMaker);
2009  pseg->setT0Error(rtime, rerrorTime);
2010 
2011  return pseg;
2012 }
2013 
2015 void CscSegmentUtilTool::get2dSegments(Identifier eta_id, Identifier phi_id, ICscSegmentFinder::ChamberTrkClusters& eta_clus,
2016  ICscSegmentFinder::ChamberTrkClusters& phi_clus, ICscSegmentFinder::Segments& eta_segs,
2017  ICscSegmentFinder::Segments& phi_segs, const Amg::Vector3D& lpos000, const EventContext& ctx,
2018  int etaStat, int phiStat) const {
2019  if (!eta_id.is_valid() && !phi_id.is_valid()) {
2020  ATH_MSG_WARNING("in get2dSegments: got two invalid identifiers");
2021  return;
2022  }
2023  // check whether Id is valid
2024  Identifier chId = eta_id.is_valid() ? eta_id : phi_id;
2025  int col_station = m_idHelperSvc->cscIdHelper().stationName(chId) - 49;
2026  int col_eta = m_idHelperSvc->cscIdHelper().stationEta(chId);
2027  int col_phisec = m_idHelperSvc->cscIdHelper().stationPhi(chId);
2028 
2029  ATH_MSG_DEBUG("get2dSegments called " << eta_id << " " << phi_id << " " << col_station << " " << col_eta << " " << col_phisec
2030  << " " << etaStat << " " << phiStat << " "
2031  << " Counts: " << eta_clus[0].size() << " " << eta_clus[1].size() << " " << eta_clus[2].size()
2032  << " " << eta_clus[3].size());
2033 
2034  ICscSegmentFinder::Segments eta_segs3hit, phi_segs3hit;
2035 
2036  double pos_eta = -999;
2037  double slope_eta = -999;
2038  double pos_phi = -999;
2039  double slope_phi = -999;
2040 
2041  // Find 2D segments.
2042  find_2dsegments(false, col_station, col_eta, col_phisec, eta_clus, lpos000, eta_segs, pos_eta, slope_eta, ctx);
2043  find_2dsegments(true, col_station, col_eta, col_phisec, phi_clus, lpos000, phi_segs, pos_phi, slope_phi, ctx);
2044 
2045  // Find 3-hit 2D segments.
2046  find_2dseg3hit(false, col_station, col_eta, col_phisec, eta_clus, lpos000, eta_segs3hit, eta_segs, pos_eta, slope_eta, ctx);
2047  find_2dseg3hit(true, col_station, col_eta, col_phisec, phi_clus, lpos000, phi_segs3hit, phi_segs, pos_phi, slope_phi, ctx);
2048 
2049  // Add 3-hit segments to 4-hit segments.
2050  add_2dsegments(eta_segs, eta_segs3hit);
2051  add_2dsegments(phi_segs, phi_segs3hit);
2052 
2053  int nGoodEta = 0, nGoodPhi = 0;
2054  for (int i = 0; i < 4; ++i) {
2055  if ((etaStat % (int)std::pow(10, i + 1)) / (int)std::pow(10, i) == 0) nGoodEta++;
2056  if ((phiStat % (int)std::pow(10, i + 1)) / (int)std::pow(10, i) == 0) nGoodPhi++;
2057  }
2058  if (nGoodEta == 2) {
2059  // Find 2-hit 2D segments for eta.
2060  ICscSegmentFinder::Segments eta_segs2hit;
2061 
2062  ATH_MSG_VERBOSE(" start find_2dseg2hit eta ");
2063  find_2dseg2hit(false, col_station, col_eta, col_phisec, etaStat, eta_clus, lpos000, eta_segs2hit, pos_eta, slope_eta, ctx);
2064 
2065  // store 2-hit segments
2066  ATH_MSG_VERBOSE(" store 2hit eta segments");
2067  add_2dseg2hits(eta_segs, eta_segs2hit, etaStat);
2068  }
2069  if (nGoodPhi == 2) {
2070  // Find 2-hit 2D segments for eta.
2071  ICscSegmentFinder::Segments phi_segs2hit;
2072  ATH_MSG_VERBOSE(" start find_2dseg2hit phi ");
2073  find_2dseg2hit(true, col_station, col_eta, col_phisec, phiStat, phi_clus, lpos000, phi_segs2hit, pos_phi, slope_phi, ctx);
2074  ATH_MSG_VERBOSE(" store 2hit phi segments");
2075  add_2dseg2hits(phi_segs, phi_segs2hit, phiStat);
2076  }
2077 }
2078 
2079 //***********************************************************************
2080 // By Elliott
2082 bool CscSegmentUtilTool::unique_hits(ICscSegmentFinder::TrkClusters& fitclus,
2083  ICscSegmentFinder::Segments& segs) const { // segs => 4hitsegs...
2084 
2085  int nclus = fitclus.size();
2086  ATH_MSG_VERBOSE(" unique_hits nclus " << nclus << " segs size " << segs.size());
2087  // Loop over segments
2088  ICscSegmentFinder::Segments::const_iterator iseg;
2089  for (iseg = segs.begin(); iseg != segs.end(); ++iseg) {
2090  int nhits_common = 0;
2091  // Loop over hits in the 4-hit segment
2092  for (int iclus = 0; iclus < iseg->nclus; iclus++) {
2093  // Loop over hits in fitclus
2094  ICscSegmentFinder::TrkClusters::const_iterator ihit;
2095  for (ihit = fitclus.begin(); ihit != fitclus.end(); ++ihit) {
2096  // Increment number of common hits if the ids are the same.
2097  const Muon::CscClusterOnTrack* cot = iseg->clus[iclus].cl;
2098  const Muon::CscClusterOnTrack* hit = ihit->cl;
2099  if (cot->identify() == hit->identify()) nhits_common++;
2100  }
2101  // if (nhits_common == nclus) return false;
2102  if (nhits_common == nclus || (nclus >= 2 && nhits_common > m_max_3hitseg_sharehit)) return false; // WP 4/5/10
2103  }
2104  }
2105  // If we get here, then the hits from the 3-hit segment are unique.
2106  return true;
2107 }
2108 
2110 // if hit is in outlier, error is estimated in width/sqrt(12)
2111 void CscSegmentUtilTool::getRios(const ICscSegmentFinder::Segment& seg, ICscSegmentFinder::MbaseList* prios, bool measphi,
2112  const EventContext& ctx) const {
2113  // Getting CscPrepData
2114  for (int iclu = 0; iclu < seg.nclus; ++iclu) {
2115  const CscClusterOnTrack* pclu = seg.clus[iclu].cl;
2116 
2117  if (pclu != nullptr) {
2118  // if this cluster is picked up as outlier cluster
2119  if (seg.outlierid == iclu) {
2120  Identifier id = pclu->identify();
2121 
2122  // get prep raw data.
2123  const CscPrepData* prd = pclu->prepRawData();
2124  double err = getDefaultError(id, measphi, prd, ctx);
2125 
2126  Amg::MatrixX cov(1, 1);
2127  cov(0, 0) = err * err;
2128 
2129  // Create new calibrated hit with new error matrix.
2130  Amg::Vector2D lpos = prd->localPosition();
2131  Trk::DefinedParameter locPar(lpos.x(), Trk::locX);
2132  Trk::LocalParameters ppars(locPar);
2133 
2134  Trk::RIO_OnTrack* pclu2 = new CscClusterOnTrack(prd, std::move(ppars), std::move(cov), 0.0, prd->status(), prd->timeStatus(), prd->time());
2135  prios->push_back(pclu2);
2136  }
2137  // Just put into list of RIO_OnTrack
2138  else {
2139  Trk::RIO_OnTrack* pclu2 = pclu->clone();
2140  prios->push_back(pclu2);
2141  }
2142  }
2143  }
2144 }
2145 
2147 double CscSegmentUtilTool::getDefaultError(Identifier id, bool measphi, const CscPrepData* prd, const EventContext& ctx) const {
2148  const std::vector<Identifier>& strip_ids = prd->rdoList();
2149  unsigned int nstrip = strip_ids.size();
2150  SG::ReadCondHandle<MuonGM::MuonDetectorManager> DetectorManagerHandle{m_DetectorManagerKey, ctx};
2151  const MuonGM::MuonDetectorManager* MuonDetMgr{*DetectorManagerHandle};
2152  if (MuonDetMgr == nullptr) {
2153  ATH_MSG_ERROR("Null pointer to the read MuonDetectorManager conditions object");
2154  return 0.;
2155  }
2156  const CscReadoutElement* pro = MuonDetMgr->getCscReadoutElement(id);
2157  double pitch = pro->cathodeReadoutPitch(0, measphi);
2158  // Assign position error.
2159  double wmeas = pitch * nstrip;
2160  double weff = wmeas - 20;
2161  double weffmin = 0.5 * wmeas;
2162  if (weff < weffmin) weff = weffmin;
2163  double err = weff / std::sqrt(12.0);
2164  return err;
2165 }
2166 
2167 // Likelihood calculation for xy matching.
2168 double CscSegmentUtilTool::matchLikelihood(const MuonSegment& rsg, const MuonSegment& psg) const {
2169  // Loop over eta and phi segments.
2170  ICscSegmentFinder::RioList etarios;
2171  for (unsigned int irot = 0; irot < rsg.numberOfContainedROTs(); irot++) etarios.push_back(rsg.rioOnTrack(irot));
2172  ICscSegmentFinder::RioList phirios;
2173  for (unsigned int irot = 0; irot < psg.numberOfContainedROTs(); irot++) phirios.push_back(psg.rioOnTrack(irot)); // CHECK THIS
2174  int maxeta = etarios.size();
2175  int maxphi = phirios.size();
2176 
2177  double prob_real = 1, prob_ghost = 1;
2178  bool match = false;
2179  for (int ieta = 0; ieta < maxeta; ieta++) {
2180  for (int iphi = 0; iphi < maxphi; iphi++) {
2181  const Trk::RIO_OnTrack* etapold = etarios[ieta];
2182  const Trk::RIO_OnTrack* phipold = phirios[iphi];
2183  double qratio;
2184 
2185  const Muon::CscClusterOnTrack* eta_clu = dynamic_cast<const Muon::CscClusterOnTrack*>(etapold);
2186  if (!eta_clu) {
2187  ATH_MSG_DEBUG(" continuing because not a CSC eta cluster");
2188  continue;
2189  }
2190 
2191  const Muon::CscClusterOnTrack* phi_clu = dynamic_cast<const Muon::CscClusterOnTrack*>(phipold);
2192  if (!phi_clu) {
2193  ATH_MSG_DEBUG(" continuing because not a CSC phi cluster");
2194  continue;
2195  }
2196 
2197  // Get prep raw data
2198  const Muon::CscPrepData* eta_prd = eta_clu->prepRawData();
2199  const Muon::CscPrepData* phi_prd = phi_clu->prepRawData();
2200 
2201  // Figure out the layer for each view
2202  Identifier eta_sid = eta_prd->identify();
2203  Identifier phi_sid = phi_prd->identify();
2204  int eta_wlay = m_idHelperSvc->cscIdHelper().wireLayer(eta_sid);
2205  int phi_wlay = m_idHelperSvc->cscIdHelper().wireLayer(phi_sid);
2206 
2207  // Calculate the charge ratio for same eta & phi layers
2208  if (eta_wlay == phi_wlay) {
2209  // Get charge.
2210  double eta_chg = eta_prd->charge();
2211  double phi_chg = phi_prd->charge();
2212 
2213  if (eta_chg > 0 && phi_chg > 0) {
2214  qratio = eta_chg / phi_chg;
2215  prob_real *= pdf_sig(qratio);
2216  prob_ghost *= pdf_bkg(qratio);
2217  match = true;
2218  }
2219  } // if (eta_wlay)
2220  } // Loop over phi hits
2221  } // Loop over eta hits
2222 
2223  // Determine likelihood
2224  double like = qratio_like(prob_real, prob_ghost);
2225  if (m_allEtaPhiMatches) {
2226  if (match) ATH_MSG_VERBOSE(" MATCH allEtaPhiMatches " << like);
2227  if (!match) ATH_MSG_VERBOSE(" NO MATCH allEtaPhiMatches " << like);
2228  if (match) return 0.5;
2229  if (!match) return 0.;
2230  }
2231  return like;
2232 }
2233 
2234 // Function to return pdf value for signal distribution in xy matching.
2235 double CscSegmentUtilTool::pdf_sig(const double x) {
2236  double f1, f2;
2237  double par[6] = {1.25049, 1.02934, 0.0517436, 0.0229711, 0.900799, 0.374422};
2238 
2239  // Overflow bin.
2240  if (x > 10) return 0.015619;
2241 
2242  // Landau
2243  f1 = par[0] * TMath::Landau(x, par[1], par[2]);
2244 
2245  // gaussian
2246  f2 = par[3] * TMath::Gaus(x, par[4], par[5]);
2247 
2248  // function is the sum of a landau plus gaussian distribution.
2249  return f1 + f2;
2250 }
2251 
2252 // function to return pdf value for background distribution.
2253 double CscSegmentUtilTool::pdf_bkg(const double x) {
2254  double e1, e2;
2255  double par[8] = {0.0394188, 0.0486057, 0.0869231, 1.16153, -0.109998, 0.009729, 0.36183, 0.228344};
2256 
2257  // Overflow bin
2258  if (x > 10) return 0.0895146;
2259 
2260  // Treat values below 0.2 as flat bins.
2261  if (x < 0.1)
2262  e1 = par[0];
2263  else if (x < 0.2)
2264  e1 = par[1];
2265  else
2266  e1 = par[2] * exp(-1 * par[3] * (x - par[4]));
2267 
2268  // 2nd exponential
2269  e2 = par[5] * exp(-1 * par[6] * (x - par[7]));
2270 
2271  // value is the sum of two functions.
2272  return e1 + e2;
2273 }
2274 
2275 // likelihood function for charge ratio
2276 double CscSegmentUtilTool::qratio_like(const double pdf_sig, const double pdf_bkg) {
2277  double like = 0;
2278 
2279  // return zero if both probability distribution functions are zero.
2280  if (pdf_sig != 0 && pdf_bkg != 0) like = pdf_sig / (pdf_sig + pdf_bkg);
2281 
2282  return like;
2283 }
2284 
2285 bool CscSegmentUtilTool::isGood(const uint32_t stripHashId, const EventContext& ctx) const {
2286  unsigned int status = stripStatusBit(stripHashId, ctx);
2287  bool is_good = !((status & 0x1) || ((status >> 1) & 0x1)); // test for hot/dead channel
2288  return is_good;
2289 }
2290 
2291 int CscSegmentUtilTool::stripStatusBit(const uint32_t stripHashId, const EventContext& ctx) const {
2292  SG::ReadCondHandle<CscCondDbData> readHandle{m_readKey, ctx};
2293  const CscCondDbData* readCdo{*readHandle};
2294 
2295  int status = 0;
2296  if (!readCdo->readChannelStatus(stripHashId, status).isSuccess())
2297  ATH_MSG_WARNING(" failed to access CSC conditions database - status - "
2298  << "strip hash id = " << stripHashId);
2299  return status;
2300 }
Trk::SpaceTimePointBase::errorTime
float errorTime() const
access to the error on the measured time
Definition: SpaceTimePointBase.h:50
TrapezoidBounds.h
IDTPM::ndof
float ndof(const U &p)
Definition: TrackParametersHelper.h:142
Trk::LocalParameters
Definition: LocalParameters.h:98
CscSegmentUtilTool::stripStatusBit
int stripStatusBit(const uint32_t stripHashId, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:2291
MuonGM::MuonClusterReadoutElement::transform
virtual const Amg::Transform3D & transform() const override
Return local to global transform.
Definition: MuonClusterReadoutElement.h:124
ReadCellNoiseFromCoolCompare.s1
s1
Definition: ReadCellNoiseFromCoolCompare.py:378
ICscSegmentFinder::Segment::d1
double d1
Definition: ICscSegmentFinder.h:95
Muon::MuonSegmentCombination
Definition: MuonSegmentCombination.h:30
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
CscSegmentUtilTool::m_max_3hitseg_sharehit
Gaudi::Property< int > m_max_3hitseg_sharehit
Definition: CscSegmentUtilTool.h:93
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:29
yodamerge_tmp.dim
dim
Definition: yodamerge_tmp.py:239
CscSegmentUtilTool::m_DetectorManagerKey
SG::ReadCondHandleKey< MuonGM::MuonDetectorManager > m_DetectorManagerKey
Definition: CscSegmentUtilTool.h:112
Muon::MuonSegmentCombination::setNGoodCscLayers
void setNGoodCscLayers(int nEta, int nPhi)
Definition: MuonSegmentCombination.h:64
CscSegmentUtilTool::get2dMuonSegmentCombination
Muon::MuonSegmentCombination * get2dMuonSegmentCombination(Identifier eta_id, Identifier phi_id, ICscSegmentFinder::ChamberTrkClusters &eta_clus, ICscSegmentFinder::ChamberTrkClusters &phi_clus, const Amg::Vector3D &lpos000, const EventContext &ctx, int etaStat=0, int phiStat=0) const
Definition: CscSegmentUtilTool.cxx:148
CscSegmentUtilTool::get4dMuonSegmentCombination
Muon::MuonSegmentCombination * get4dMuonSegmentCombination(Identifier eta_id, Identifier phi_id, ICscSegmentFinder::ChamberTrkClusters &eta_clus, ICscSegmentFinder::ChamberTrkClusters &phi_clus, const Amg::Vector3D &lpos000, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:1738
python.rtime.rtime
def rtime(fn, *args, **kw)
Definition: rtime.py:15
ICscSegmentFinder::Segment::nunspoil
int nunspoil
Definition: ICscSegmentFinder.h:101
Muon::MuonClusterOnTrack::globalPosition
virtual const Amg::Vector3D & globalPosition() const override
Returns global position.
Definition: MuonClusterOnTrack.cxx:93
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
ICscSegmentFinder::Cluster
Definition: ICscSegmentFinder.h:59
Trk::locX
@ locX
Definition: ParamDefs.h:43
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:44
CscSegmentUtilTool::CscSegmentUtilTool
CscSegmentUtilTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: CscSegmentUtilTool.cxx:78
CscSegmentUtilTool::find_2dseg3hit
void find_2dseg3hit(bool measphi, int station, int eta, int phi, const ICscSegmentFinder::ChamberTrkClusters &clus, const Amg::Vector3D &lpos000, ICscSegmentFinder::Segments &segs, ICscSegmentFinder::Segments &segs4hit, double localPos, double localSlope, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:1403
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
ICscSegmentFinder::Segment::outlierid
int outlierid
Definition: ICscSegmentFinder.h:105
ICscSegmentFinder::Segment
Definition: ICscSegmentFinder.h:90
Muon::CscPrepData::charge
int charge() const
Returns the charge.
Definition: CscPrepData.h:153
Trk::RIO_OnTrack::clone
virtual RIO_OnTrack * clone() const override=0
Pseudo-constructor, needed to avoid excessive RTTI.
xAOD::uint32_t
setEventNumber uint32_t
Definition: EventInfo_v1.cxx:127
egammaEnergyPositionAllSamples::e1
double e1(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 1st sampling
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
Muon::MuonSegmentCombination::useStripsInSegment
bool useStripsInSegment(int isEta) const
Definition: MuonSegmentCombination.h:92
hist_file_dump.d
d
Definition: hist_file_dump.py:137
Trk::SurfaceBounds
Definition: SurfaceBounds.h:47
IsUnspoiled
bool IsUnspoiled(Muon::CscClusterStatus status)
Definition: CscDataPreparator.cxx:165
Trk::RotatedTrapezoidBounds
Definition: RotatedTrapezoidBounds.h:45
CscSegmentUtilTool::get2dSegments
void get2dSegments(Identifier eta_id, Identifier phi_id, ICscSegmentFinder::ChamberTrkClusters &eta_clus, ICscSegmentFinder::ChamberTrkClusters &phi_clus, ICscSegmentFinder::Segments &etasegs, ICscSegmentFinder::Segments &phisegs, const Amg::Vector3D &lpos000, const EventContext &ctx, int etaStat=0, int phiStat=0) const
Definition: CscSegmentUtilTool.cxx:2015
Muon::MuonSegment::rioOnTrack
const Trk::RIO_OnTrack * rioOnTrack(unsigned int) const
returns the RIO_OnTrack (also known as ROT) objects depending on the integer
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:187
Trk::ParamDefs
ParamDefs
Definition: ParamDefs.h:38
CscSegmentUtilTool::fit_rio_segment
void fit_rio_segment(const Trk::PlaneSurface &ssrf, bool dump, const ICscSegmentFinder::RioList &clus, double &s0, double &s1, double &d0, double &d1, double &d01, double &chsq, double &zshift) const
Definition: CscSegmentUtilTool.cxx:504
EventPrimitivesHelpers.h
ICscSegmentFinder::ChamberTrkClusters
std::vector< Cluster > ChamberTrkClusters[4]
Definition: ICscSegmentFinder.h:83
ICscSegmentFinder::Segment::d01
double d01
Definition: ICscSegmentFinder.h:96
Trk::RIO_OnTrack::globalPosition
virtual const Amg::Vector3D & globalPosition() const override=0
Interface method to get the global Position.
ICscSegmentFinder::Segment::s1
double s1
Definition: ICscSegmentFinder.h:93
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
MuonGM::MuonReadoutElement::GlobalToAmdbLRSTransform
virtual Amg::Transform3D GlobalToAmdbLRSTransform() const
Definition: MuonDetDescr/MuonReadoutGeometry/src/MuonReadoutElement.cxx:153
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
ICscStripFitter::Result::charge
double charge
Definition: ICscStripFitter.h:31
CscSegmentUtilTool::m_IPerror
Gaudi::Property< double > m_IPerror
Definition: CscSegmentUtilTool.h:88
Trk::PrepRawData::rdoList
const std::vector< Identifier > & rdoList() const
return the List of rdo identifiers (pointers)
CscStripPrepData.h
CscSegmentUtilTool::m_max_seg_per_chamber
Gaudi::Property< unsigned int > m_max_seg_per_chamber
Definition: CscSegmentUtilTool.h:95
Muon::MuonSegmentCombination::stationSegments
SegmentVec * stationSegments(unsigned int index) const
Access to segments in a given station.
Definition: MuonSegmentCombination.h:114
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
dq_defect_virtual_defect_validation.d1
d1
Definition: dq_defect_virtual_defect_validation.py:79
CscSegmentUtilTool::matchLikelihood
double matchLikelihood(const Muon::MuonSegment &rsg, const Muon::MuonSegment &psg) const
Definition: CscSegmentUtilTool.cxx:2168
CscSegmentUtilTool::initialize
virtual StatusCode initialize()
Definition: CscSegmentUtilTool.cxx:84
Trk::loc2
@ loc2
generic first and second local coordinate
Definition: ParamDefs.h:41
CscSegmentUtilTool::isGood
bool isGood(const uint32_t stripHashId, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:2285
CscSegmentUtilTool::fit_rio_residual
void fit_rio_residual(const Trk::PlaneSurface &ssrf, bool dump, const ICscSegmentFinder::RioList &clus, unsigned int irclu, double &res, double &dres, double &rs, double &drs) const
Definition: CscSegmentUtilTool.cxx:603
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
Trk::RIO_OnTrack
Definition: RIO_OnTrack.h:70
CscSegmentUtilTool::m_cluster_error_scaler
Gaudi::Property< float > m_cluster_error_scaler
Definition: CscSegmentUtilTool.h:90
CscSegmentUtilTool::m_TightenChi2
Gaudi::Property< bool > m_TightenChi2
Definition: CscSegmentUtilTool.h:101
CscSegmentUtilTool::m_max_slope_r
Gaudi::Property< double > m_max_slope_r
Definition: CscSegmentUtilTool.h:84
python.TurnDataReader.dr
dr
Definition: TurnDataReader.py:112
ICscSegmentFinder::MbaseList
DataVector< const Trk::MeasurementBase > MbaseList
Definition: ICscSegmentFinder.h:86
MuonGM::CscReadoutElement
Definition: CscReadoutElement.h:56
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
CscSegmentUtilTool.h
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
xAOD::MuonSegment
MuonSegment_v1 MuonSegment
Reference the current persistent version:
Definition: Event/xAOD/xAODMuon/xAODMuon/MuonSegment.h:13
empty
bool empty(TH1 *h)
Definition: computils.cxx:294
Muon::MuonSegment::numberOfContainedROTs
unsigned int numberOfContainedROTs() const
number of RIO_OnTracks
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:199
CscSegmentUtilTool::fit_detailCalcPart2
void fit_detailCalcPart2(double q0, double q1, double q2, double q01, double q11, double q02, double &s0, double &s1, double &d0, double &d1, double &d01, double &chsq) const
Definition: CscSegmentUtilTool.cxx:1597
ParamDefs.h
Identifier::is_valid
bool is_valid() const
Check if id is in a valid state.
ICscSegmentFinder::Segment::dtime
double dtime
Definition: ICscSegmentFinder.h:98
CscSegmentUtilTool::m_fitsegment_tantheta_tolerance
Gaudi::Property< double > m_fitsegment_tantheta_tolerance
Definition: CscSegmentUtilTool.h:87
ICscSegmentFinder::Segment::nclus
int nclus
Definition: ICscSegmentFinder.h:102
Trk::DefinedParameter
std::pair< double, ParamDefs > DefinedParameter
Definition: DefinedParameter.h:27
CscSegmentUtilTool::m_remove3Overlap
Gaudi::Property< bool > m_remove3Overlap
Definition: CscSegmentUtilTool.h:103
Trk::LocalDirection::angleYZ
double angleYZ() const
access method for angle of local YZ projection
Definition: LocalDirection.h:106
MuonSegmentCombination.h
read_hist_ntuple.f2
f2
Definition: read_hist_ntuple.py:20
CscSegmentUtilTool::build_segment
Muon::MuonSegment * build_segment(const ICscSegmentFinder::Segment &seg, bool measphi, Identifier chid, bool use2Lay, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:762
Muon::CscPrepData::timeStatus
CscTimeStatus timeStatus() const
Returns the Csc time status flag.
Definition: CscPrepData.h:168
CscSegmentUtilTool::find_2dsegments
void find_2dsegments(bool measphi, int station, int eta, int phi, const ICscSegmentFinder::ChamberTrkClusters &clus, const Amg::Vector3D &lpos000, ICscSegmentFinder::Segments &segs, double localPos, double localSlope, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:1046
PURW_create_actual_mu_profile.rdir
rdir
Definition: PURW_create_actual_mu_profile.py:62
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
CscSegmentUtilTool::getRios
void getRios(const ICscSegmentFinder::Segment &seg, ICscSegmentFinder::MbaseList *prios, bool measphi, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:2111
ICscSegmentFinder::sortByNunspoilAndChsq
Definition: ICscSegmentFinder.h:134
CscSegmentUtilTool::fit_residual
void fit_residual(const ICscSegmentFinder::TrkClusters &clus, const Amg::Vector3D &lpos000, unsigned int irclu, double &res, double &dres, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:450
Muon::CscTimeLate
@ CscTimeLate
Not successful time measurement but samples indicates late time above 200ns in RAW time.
Definition: CscTimeStatus.h:33
ICscClusterFitter.h
ICscSegmentFinder::Segment::d0
double d0
Definition: ICscSegmentFinder.h:94
CscSegmentUtilTool::find_2dseg2hit
void find_2dseg2hit(bool measphi, int station, int eta, int phi, int layStat, const ICscSegmentFinder::ChamberTrkClusters &clus, const Amg::Vector3D &lpos000, ICscSegmentFinder::Segments &segs, double localPos, double localSlope, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:1504
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CscPrepData.h
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:193
ICscSegmentFinder::Segment::clus
TrkClusters clus
Definition: ICscSegmentFinder.h:104
CscSegmentUtilTool::pdf_bkg
static double pdf_bkg(const double x)
Definition: CscSegmentUtilTool.cxx:2253
lumiFormat.i
int i
Definition: lumiFormat.py:92
CscSegmentUtilTool::pdf_sig
static double pdf_sig(const double x)
Definition: CscSegmentUtilTool.cxx:2235
ReadCalibFromCool.keep
keep
Definition: ReadCalibFromCool.py:85
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
Muon::CscPrepData
Class representing clusters from the CSC.
Definition: CscPrepData.h:39
Trk::theta
@ theta
Definition: ParamDefs.h:72
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
Muon::CscClusterOnTrack::clone
virtual CscClusterOnTrack * clone() const override final
Clone this ROT.
Definition: CscClusterOnTrack.h:149
MuonSegmentCombinationCollection.h
Muon::MuonSegment::localDirection
const Trk::LocalDirection & localDirection() const
local direction
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:169
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
createCablingJSON.station_name
int station_name
Simple script to generate a BIS78 cabling map as used for the Monte Carlo processing.
Definition: createCablingJSON.py:8
Muon::CscPrepData::status
CscClusterStatus status() const
Returns the Csc status (position measurement) flag.
Definition: CscPrepData.h:163
CscSegmentUtilTool::m_allEtaPhiMatches
Gaudi::Property< bool > m_allEtaPhiMatches
Definition: CscSegmentUtilTool.h:100
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:14
CscClusterOnTrack.h
CscSegmentUtilTool::fit_segment
void fit_segment(const ICscSegmentFinder::TrkClusters &clus, const Amg::Vector3D &lpos000, double &s0, double &s1, double &d0, double &d1, double &d01, double &chsq, double &time, double &dtime, double &zshift, const EventContext &ctx, int outlierLayer=-1) const
Definition: CscSegmentUtilTool.cxx:209
CscSegmentUtilTool::m_max_chisquare
Gaudi::Property< double > m_max_chisquare
Definition: CscSegmentUtilTool.h:83
CscReadoutElement.h
TRT::Track::d0
@ d0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:62
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Muon::MuonSegment::setT0Error
void setT0Error(float t0, float t0Error)
set the fitted time and error on the time
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:205
ICscSegmentFinder::TrkClusters
std::vector< Cluster > TrkClusters
Definition: ICscSegmentFinder.h:82
Trk::Surface::normal
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
ICscSegmentFinder::SegmentVec
Muon::MuonSegmentCombination::SegmentVec SegmentVec
Definition: ICscSegmentFinder.h:88
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Muon::CscPrepData::time
double time() const
Returns the time.
Definition: CscPrepData.h:158
Trk::Segment::Csc2dSegmentMaker
@ Csc2dSegmentMaker
Definition: TrkEvent/TrkSegment/TrkSegment/Segment.h:71
Trk::FitQuality
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition: FitQuality.h:97
ICscStripFitter::Result
Definition: ICscStripFitter.h:25
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
ICscSegmentFinder::Segments
std::vector< Segment > Segments
Definition: ICscSegmentFinder.h:132
python.ChapPy.dump
def dump(buf, stdout=sys.stdout)
Definition: ChapPy.py:25
DataVector< const Trk::MeasurementBase >
Trk::LocalDirection
represents the three-dimensional global direction with respect to a planar surface frame.
Definition: LocalDirection.h:81
CscSegmentUtilTool::m_max_slope_phi
Gaudi::Property< double > m_max_slope_phi
Definition: CscSegmentUtilTool.h:85
Trk::MeasurementBase::localCovariance
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
Definition: MeasurementBase.h:138
ICscSegmentFinder::Segment::s0
double s0
Definition: ICscSegmentFinder.h:92
CscSegmentUtilTool::m_max_chisquare_loose
Gaudi::Property< double > m_max_chisquare_loose
Definition: CscSegmentUtilTool.h:82
Muon::CscClusterOnTrack::time
float time() const
Return the time(ns)
Definition: CscClusterOnTrack.h:188
Muon::CscTimeSuccess
@ CscTimeSuccess
Time measurement is successful to use.
Definition: CscTimeStatus.h:27
CscSegmentUtilTool::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: CscSegmentUtilTool.h:107
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
Trk::Segment::Csc4dSegmentMaker
@ Csc4dSegmentMaker
Definition: TrkEvent/TrkSegment/TrkSegment/Segment.h:72
CscSegmentUtilTool::m_nunspoil
Gaudi::Property< int > m_nunspoil
Definition: CscSegmentUtilTool.h:92
CscPrepDataContainer.h
CscSegmentUtilTool::getDefaultError
double getDefaultError(Identifier id, bool measphi, const Muon::CscPrepData *prd, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:2147
ICscSegmentFinder::Segment::chsq
double chsq
Definition: ICscSegmentFinder.h:99
Trk::SurfaceBounds::clone
virtual SurfaceBounds * clone() const =0
clone() method to make deep copy in Surface copy constructor and for assigment operator of the Surfac...
CscSegmentUtilTool::qratio_like
static double qratio_like(const double pdf_sig, const double pdf_bkg)
Definition: CscSegmentUtilTool.cxx:2276
CscSegmentUtilTool::add_2dseg2hits
void add_2dseg2hits(ICscSegmentFinder::Segments &segs, ICscSegmentFinder::Segments &segs2, int layStat) const
Definition: CscSegmentUtilTool.cxx:1279
Muon::CscStatusUnspoiled
@ CscStatusUnspoiled
Clean cluster with precision fit.
Definition: CscClusterStatus.h:26
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
Amg::error
double error(const Amg::MatrixX &mat, int index)
return diagonal error of the matrix caller should ensure the matrix is symmetric and the index is in ...
Definition: EventPrimitivesHelpers.h:40
Trk::PrepRawData::localPosition
const Amg::Vector2D & localPosition() const
return the local position reference
threshold
Definition: chainparser.cxx:74
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
CscSegmentUtilTool::m_eventInfo
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo
Definition: CscSegmentUtilTool.h:110
Trk::Surface::insideBounds
virtual bool insideBounds(const Amg::Vector2D &locpos, double tol1=0., double tol2=0.) const =0
virtual methods to be overwritten by the inherited surfaces
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::RIO_OnTrack::prepRawData
virtual const Trk::PrepRawData * prepRawData() const =0
returns the PrepRawData (also known as RIO) object to which this RIO_OnTrack is associated.
Trk::MeasurementBaseType::RIO_OnTrack
@ RIO_OnTrack
Definition: MeasurementBase.h:49
Muon::MuonSegmentCombination::use2LayerSegments
bool use2LayerSegments(int isEta) const
Definition: MuonSegmentCombination.h:86
Muon::CscClusterOnTrack::prepRawData
virtual const CscPrepData * prepRawData() const override final
Returns the CscPrepData - is a CscPrepData in this scope.
Definition: CscClusterOnTrack.h:154
Trk::TrapezoidBounds
Definition: TrapezoidBounds.h:43
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
CscSegmentUtilTool::m_zshift
Gaudi::Property< bool > m_zshift
Definition: CscSegmentUtilTool.h:98
Trk::MeasurementBase::localParameters
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
Definition: MeasurementBase.h:132
Muon::MuonSegmentCombination::addSegments
bool addSegments(std::unique_ptr< SegmentVec >)
Add a set of Segments for a give station.
Definition: MuonSegmentCombination.h:102
CscSegmentUtilTool::m_readKey
SG::ReadCondHandleKey< CscCondDbData > m_readKey
Definition: CscSegmentUtilTool.h:109
Muon::MuonSegment::clone
virtual MuonSegment * clone() const override final
needed to avoid excessive RTTI
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:181
CscSegmentUtilTool::unique_hits
bool unique_hits(ICscSegmentFinder::TrkClusters &fitclus, ICscSegmentFinder::Segments &segs) const
Method for checking whether three hit segments are already part of 4 hit segments.
Definition: CscSegmentUtilTool.cxx:2082
CscSegmentUtilTool::m_min_xylike
Gaudi::Property< double > m_min_xylike
Definition: CscSegmentUtilTool.h:86
CscSegmentUtilTool::getMuonSegments
std::unique_ptr< std::vector< std::unique_ptr< Muon::MuonSegment > > > getMuonSegments(Identifier eta_id, Identifier phi_id, ICscSegmentFinder::ChamberTrkClusters &eta_clus, ICscSegmentFinder::ChamberTrkClusters &phi_clus, const Amg::Vector3D &lpos000, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:116
Muon::CscClusterOnTrack
Class to represent the calibrated clusters created from CSC strips.
Definition: CscClusterOnTrack.h:47
y
#define y
Muon::CscStatusSplitUnspoiled
@ CscStatusSplitUnspoiled
Clean cluster with precision fit after split cluster.
Definition: CscClusterStatus.h:59
CscSegmentUtilTool::fit_detailCalcPart1
void fit_detailCalcPart1(const ICscSegmentFinder::TrkClusters &clus, const Amg::Vector3D &lpos000, double &s0, double &s1, double &d0, double &d1, double &d01, double &chsq, bool &measphi, double &time, double &dtime, double &zshift, bool IsSlopeGive, int outlierHitLayer, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:226
Trk::Surface::globalToLocal
virtual bool globalToLocal(const Amg::Vector3D &glob, const Amg::Vector3D &mom, Amg::Vector2D &loc) const =0
Specified by each surface type: GlobalToLocal method without dynamic memory allocation - boolean chec...
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
egammaEnergyPositionAllSamples::e2
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
trigbs_pickEvents.cnt
cnt
Definition: trigbs_pickEvents.py:71
MuonGM::MuonDetectorManager
The MuonDetectorManager stores the transient representation of the Muon Spectrometer geometry and pro...
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonDetectorManager.h:49
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
PlaneSurface.h
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
Trk::PlaneSurface::bounds
virtual const SurfaceBounds & bounds() const override final
This method returns the bounds by reference, static NoBounds in case of no boundaries.
CscSegmentUtilTool::make_4dMuonSegment
Muon::MuonSegment * make_4dMuonSegment(const Muon::MuonSegment &rsg, const Muon::MuonSegment &psg, bool use2LaySegsEta, bool use2LaySegsPhi) const
Definition: CscSegmentUtilTool.cxx:1754
DEBUG
#define DEBUG
Definition: page_access.h:11
Trk::RIO_OnTrack::identify
virtual Identifier identify() const final
return the identifier -extends MeasurementBase
Definition: RIO_OnTrack.h:155
CscSegmentUtilTool::m_remove4Overlap
Gaudi::Property< bool > m_remove4Overlap
Definition: CscSegmentUtilTool.h:102
ICscStripFitter.h
Trk::FitQuality::chiSquared
double chiSquared() const
returns the of the overall track fit
Definition: FitQuality.h:56
MuonSegment.h
CscSegmentUtilTool::m_x5data
Gaudi::Property< bool > m_x5data
Definition: CscSegmentUtilTool.h:97
Trk::phi
@ phi
Definition: ParamDefs.h:81
Trk::FitQuality::numberDoF
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
Definition: FitQuality.h:60
merge.status
status
Definition: merge.py:17
Muon::MuonSegment::globalPosition
virtual const Amg::Vector3D & globalPosition() const override final
global position
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:157
Trk::RIO_OnTrack::associatedSurface
virtual const Surface & associatedSurface() const override=0
returns the surface for the local to global transformation
LocalDirection.h
CscSegmentUtilTool::add_2dsegments
void add_2dsegments(ICscSegmentFinder::Segments &segs4, ICscSegmentFinder::Segments &segs3) const
Adds 3-hit segments to 4-hit segments.
Definition: CscSegmentUtilTool.cxx:1146
Trk::SpaceTimePointBase::time
float time() const
access to the measured time
Definition: SpaceTimePointBase.h:47
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
ICscSegmentFinder::Segment::time
double time
Definition: ICscSegmentFinder.h:97
CscSegmentUtilTool::m_max_chisquare_tight
Gaudi::Property< double > m_max_chisquare_tight
Definition: CscSegmentUtilTool.h:81
Muon::MuonSegment
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:45
Muon::CscClusterStatus
CscClusterStatus
Enum to represent the cluster status - see the specific enum values for more details.
Definition: CscClusterStatus.h:23
AthAlgTool
Definition: AthAlgTool.h:26
Trk::LocalDirection::angleXZ
double angleXZ() const
access method for angle of local XZ projection
Definition: LocalDirection.h:103
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Trk::Segment::fitQuality
const FitQuality * fitQuality() const
return the FitQuality object, returns NULL if no FitQuality is defined
Definition: TrkEvent/TrkSegment/TrkSegment/Segment.h:160
FitQuality.h
RotatedTrapezoidBounds.h
Trk::loc1
@ loc1
Definition: ParamDefs.h:40
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
CscSegmentUtilTool::m_rotCreator
ToolHandle< Muon::ICscClusterOnTrackCreator > m_rotCreator
Definition: CscSegmentUtilTool.h:105
MuonGM::MuonClusterReadoutElement::bounds
virtual const Trk::SurfaceBounds & bounds() const override
Return the boundaries of the element.
Definition: MuonClusterReadoutElement.h:127
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
NSWL1::globalToLocal
Polygon globalToLocal(const Polygon &pol, float z, const Trk::PlaneSurface &surf)
Definition: GeoUtils.cxx:103
dq_make_web_display.cl
cl
print [x.__class__ for x in toList(dqregion.getSubRegions()) ]
Definition: dq_make_web_display.py:26
ICscSegmentFinder::RioList
std::vector< const Trk::RIO_OnTrack * > RioList
Definition: ICscSegmentFinder.h:85
CscSegmentUtilTool::m_IPconstraint
Gaudi::Property< bool > m_IPconstraint
Definition: CscSegmentUtilTool.h:99
Trk::PlaneSurface::clone
virtual PlaneSurface * clone() const override
Virtual constructor.
Muon::CscTimeEarly
@ CscTimeEarly
Not successful time measurement but samples indicates early time below -50ns in RAW time.
Definition: CscTimeStatus.h:30
readCCLHist.float
float
Definition: readCCLHist.py:83
Muon::MuonSegmentCombination::getNGoodCscLayers
int getNGoodCscLayers(int isEta) const
Definition: MuonSegmentCombination.h:97
ICscSegmentFinder::Segment::zshift
double zshift
Definition: ICscSegmentFinder.h:100
Trk::Surface::localToGlobal
virtual void localToGlobal(const Amg::Vector2D &locp, const Amg::Vector3D &mom, Amg::Vector3D &glob) const =0
Specified by each surface type: LocalToGlobal method without dynamic memory allocation.
ICscSegmentFinder.h
python.Dumpers.FitQuality
FitQuality
Definition: Dumpers.py:63
match
bool match(std::string s1, std::string s2)
match the individual directories of two strings
Definition: hcg.cxx:356
Muon::MuonSegment::associatedSurface
virtual const Trk::PlaneSurface & associatedSurface() const override final
returns the surface for the local to global transformation
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:175
CscSegmentUtilTool::find_outlier_cluster
int find_outlier_cluster(const ICscSegmentFinder::TrkClusters &clus, const Amg::Vector3D &lpos000, double &returned_chsq, const EventContext &ctx) const
Definition: CscSegmentUtilTool.cxx:390
Muon::MuonSegment::globalDirection
const Amg::Vector3D & globalDirection() const
global direction
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:163
CscCondDbData
Definition: CscCondDbData.h:24
CscSegmentUtilTool::spoiled_count
void spoiled_count(const ICscSegmentFinder::RioList &rios, double threshold, int &nspoil, int &nunspoil)
Definition: CscSegmentUtilTool.cxx:694
read_hist_ntuple.f1
f1
Definition: read_hist_ntuple.py:4
MuonGM::CscReadoutElement::cathodeReadoutPitch
double cathodeReadoutPitch(int chLayer, int measuresPhi) const
Definition: CscReadoutElement.cxx:147
Generated on Sun Jun 30 2024 21:13:27 for ATLAS Offline Software by   doxygen 1.8.18