ATLAS Offline Software
MuFastPatternFinder.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "MuFastPatternFinder.h"
6 #include "xAODTrigMuon/TrigMuonDefs.h"
7 #include "MdtRegionDefiner.h"
8 #include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
9 
10 StatusCode TrigL2MuonSA::MuFastPatternFinder::initialize()
11 {
12  ATH_CHECK( m_nswPatternFinder.retrieve() );
13  ATH_CHECK( m_idHelperSvc.retrieve() );
14  return StatusCode::SUCCESS;
15 }
16 
17 // --------------------------------------------------------------------------------
18 // --------------------------------------------------------------------------------
19 
20 void TrigL2MuonSA::MuFastPatternFinder::doMdtCalibration(const EventContext& ctx, TrigL2MuonSA::MdtHitData& mdtHit, double track_phi, double phi0, bool isEndcap) const
21 {
22  int StationName = mdtHit.name;
23  int StationEta = mdtHit.StationEta;
24  int StationPhi = mdtHit.StationPhi;
25  int Multilayer = mdtHit.Multilayer;
26  int Layer = mdtHit.TubeLayer;
27  int Tube = mdtHit.Tube;
28 
29  ATH_MSG_DEBUG("...StationName/StationEta/StationPhi/Multilayer/Layer/Tube="
30  << StationName << "/" << StationEta << "/" << StationPhi << "/" << Multilayer << "/"
31  << Layer << "/" << Tube);
32 
33  Identifier id = ( mdtHit.Id.is_valid() ) ? mdtHit.Id : m_idHelperSvc->mdtIdHelper().channelID(StationName,StationEta,
34  StationPhi,Multilayer,Layer,Tube);
35  int tdcCounts = (int)mdtHit.DriftTime;
36  int adcCounts = mdtHit.Adc;
37 
38  double R = mdtHit.R;
39  const double cosDphi = std::cos(track_phi - phi0);
40  double InCo = cosDphi ? 1./ cosDphi: 0;
41  double X = (isEndcap)? R*std::cos(track_phi): R*InCo*std::cos(track_phi);
42  double Y = (isEndcap)? R*std::sin(track_phi): R*InCo*std::sin(track_phi);
43  double Z = mdtHit.Z;
44  const Amg::Vector3D point{X,Y,Z};
45 
46  MdtCalibInput calHit(id, adcCounts,tdcCounts, mdtHit.readEle);
47  calHit.setClosestApproach(point);
48  ATH_MSG_DEBUG("... MDT hit raw digit tdcCounts/adcCounts=" << tdcCounts << "/" << adcCounts);
49 
50  ATH_MSG_DEBUG("... MDT hit position X/Y/Z/track_phi/Multilayer/Layer/Tube="
51  << Amg::toString(point, 2) << "/" << track_phi << "/" << Multilayer << "/" << Layer << "/" << Tube);
52 
53  MdtCalibOutput calibOut = m_mdtCalibrationTool->calibrate(ctx, calHit, false);
54  double driftSpace = calibOut.driftRadius();
55  double driftSigma = calibOut.driftRadiusUncert();
56 
57  ATH_MSG_DEBUG("... MDT hit calibrated driftSpace/driftSigma=" << driftSpace << "/" << driftSigma);
58 
59  const double ZERO_LIMIT = 1e-4;
60 
61  if( std::abs(driftSpace) > ZERO_LIMIT ) {
62  mdtHit.DriftSpace = driftSpace;
63  mdtHit.DriftSigma = driftSigma;
64  }
65  else {
66  mdtHit.DriftSpace = 0;
67  mdtHit.DriftSigma = 0;
68  }
69 
70 }
71 
72 // --------------------------------------------------------------------------------
73 // --------------------------------------------------------------------------------
74 
75 
76 StatusCode TrigL2MuonSA::MuFastPatternFinder::findPatterns(const EventContext& ctx,
77  const TrigL2MuonSA::MuonRoad& muonRoad,
78  TrigL2MuonSA::MdtHits& mdtHits,
79  TrigL2MuonSA::StgcHits& stgcHits,
80  TrigL2MuonSA::MmHits& mmHits,
81  std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns) const
82 {
83  ATH_CHECK( findPatterns(ctx, muonRoad, mdtHits, v_trackPatterns) );
84  ATH_CHECK( m_nswPatternFinder->findPatterns(muonRoad, stgcHits, mmHits, v_trackPatterns.back()) );
85 
86  if (msgLvl(MSG::DEBUG)) {
87  for(unsigned int i_hit=0; i_hit<v_trackPatterns.back().stgcSegment.size(); i_hit++) {
88  ATH_MSG_DEBUG("PatternFinder: output sTGC hits global eta/phi/r/z/layer/channelType/isOutlier " << v_trackPatterns.back().stgcSegment[i_hit].eta << "/" << v_trackPatterns.back().stgcSegment[i_hit].phi << "/" << v_trackPatterns.back().stgcSegment[i_hit].r << "/" << v_trackPatterns.back().stgcSegment[i_hit].z << "/" << v_trackPatterns.back().stgcSegment[i_hit].layerNumber << "/" << v_trackPatterns.back().stgcSegment[i_hit].channelType << "/" << v_trackPatterns.back().stgcSegment[i_hit].isOutlier);
89  }
90 
91  for(unsigned int i_hit=0; i_hit<v_trackPatterns.back().mmSegment.size(); i_hit++) {
92  ATH_MSG_DEBUG("PatternFinder: output MM hits global eta/phi/r/z/layer/isOutlier " << v_trackPatterns.back().mmSegment[i_hit].eta << "/" << v_trackPatterns.back().mmSegment[i_hit].phi << "/" << v_trackPatterns.back().mmSegment[i_hit].r << "/" << v_trackPatterns.back().mmSegment[i_hit].z << "/" << v_trackPatterns.back().mmSegment[i_hit].layerNumber << "/" << v_trackPatterns.back().mmSegment[i_hit].isOutlier);
93  }
94  }
95  return StatusCode::SUCCESS;
96 }
97 
98 
99 StatusCode TrigL2MuonSA::MuFastPatternFinder::findPatterns(const EventContext& ctx,
100  const TrigL2MuonSA::MuonRoad& muonRoad,
101  TrigL2MuonSA::MdtHits& mdtHits,
102  std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns) const
103 {
104 
105  // find only 1 track pattern
106  v_trackPatterns.clear();
107  TrigL2MuonSA::TrackPattern trackPattern;
108  for (int i=0; i<xAOD::L2MuonParameters::Chamber::MaxChamber; i++) trackPattern.mdtSegments[i].clear();
109 
110  // Saddress = pos1/3 = 0:Large, 1:Large-SP, 2:Small, 3:Small-SP
111  trackPattern.s_address = (muonRoad.isEndcap)? -1: muonRoad.Special + 2*muonRoad.LargeSmall;
112 
113  const unsigned int MAX_STATION = 11;
114  const unsigned int MAX_LAYER = 12;
115 
116  TrigL2MuonSA::MdtLayerHits v_mdtLayerHits[MAX_STATION][MAX_LAYER];
117 
118  unsigned int i_layer_max = 0;
119  unsigned int chamber_max = 0;
120 
121  double trigger_phi = muonRoad.phiMiddle; // phi fit result to be used
122 
123  for(unsigned int i_hit=0; i_hit<mdtHits.size(); i_hit++) {
124 
125  unsigned int chamber = mdtHits[i_hit].Chamber;
126  if( chamber > chamber_max ) chamber_max = chamber;
127  if (chamber >= xAOD::L2MuonParameters::Chamber::MaxChamber) continue;
128 
129  int stationPhi = mdtHits[i_hit].StationPhi;
130  std::string name = m_idHelperSvc->mdtIdHelper().stationNameString(mdtHits[i_hit].name);
131  int chamber_this = 99;
132  int sector_this = 99;
133  bool isEndcap;
134  MdtRegionDefiner::find_station_sector(name, stationPhi, isEndcap, chamber_this, sector_this);
135  if( !isEndcap && sector_this != muonRoad.MDT_sector_trigger ) {
136  mdtHits[i_hit].isOutlier = 3;
137  continue;
138  }
139 
140  double aw = muonRoad.aw[chamber][0];
141  double bw = muonRoad.bw[chamber][0];
142  double Z = mdtHits[i_hit].Z;
143  double R = mdtHits[i_hit].R;
144  double residual = calc_residual(aw,bw,Z,R);
145 
146  mdtHits[i_hit].Residual = residual;
147  mdtHits[i_hit].isOutlier = 0;
148  unsigned int i_layer = mdtHits[i_hit].Layer;
149  double rWidth = muonRoad.rWidth[chamber][i_layer];
150 
151  ATH_MSG_DEBUG("... chamber/Z/R/aw/bw/residual/rWidth="
152  << chamber << "/" << Z << "/" << R << "/" << aw << "/" << bw << "/" << residual << "/" << rWidth);
153 
154  if( std::abs(residual) > rWidth ) {
155  mdtHits[i_hit].isOutlier = 2;
156  continue;
157  }
158 
159  ATH_MSG_DEBUG(" --> included at i_layer=" << i_layer);
160  if( mdtHits[i_hit].TubeLayer < 1 || i_layer > (MAX_LAYER-1) ) {
161  ATH_MSG_WARNING("strange i_layer=" << i_layer);
162  return StatusCode::FAILURE;
163  }
164  if( i_layer > i_layer_max ) i_layer_max = i_layer;
165 
166  v_mdtLayerHits[chamber][i_layer].indexes.push_back(i_hit);
167  v_mdtLayerHits[chamber][i_layer].ndigi++;
168  v_mdtLayerHits[chamber][i_layer].ndigi_all++;
169  v_mdtLayerHits[chamber][0].ntot++;
170  v_mdtLayerHits[chamber][0].ntot_all++;
171  v_mdtLayerHits[chamber][0].ResSum += residual;
172  }
173 
174  const double DeltaMin = 0.025;
175 
176  for(unsigned int chamber=0; chamber<=chamber_max; chamber++) {
177  ATH_MSG_DEBUG(" --- chamber=" << chamber);
178 
179  ATH_MSG_DEBUG("removing outliers...");
180  while(1) {
181  if (chamber==9) break;//BME skips this loop
182  if (chamber==10) break;//BMG skips this loop
183  unsigned int layer = 999999;
184  double DistMax = 0.;
185  double Residual = 0.;
186  unsigned int i_hit_max = 999999;
187  double ResMed = (v_mdtLayerHits[chamber][0].ntot!=0) ? v_mdtLayerHits[chamber][0].ResSum/v_mdtLayerHits[chamber][0].ntot : 0.;
188 
189  for(unsigned int i_layer=0; i_layer<=i_layer_max; i_layer++) {
190  for(unsigned int idigi=0; idigi<v_mdtLayerHits[chamber][i_layer].ndigi_all;idigi++) {
191 
192  unsigned int i_hit = v_mdtLayerHits[chamber][i_layer].indexes[idigi];
193  if(mdtHits[i_hit].isOutlier > 0) continue;
194  double DistMed = std::abs(mdtHits[i_hit].Residual - ResMed);
195  if(DistMed>=DistMax) {
196  DistMax = DistMed;
197  Residual = mdtHits[i_hit].Residual;
198  layer = i_layer;
199  i_hit_max = i_hit;
200  }
201  }
202  }
203  ATH_MSG_DEBUG("ResMed=" << ResMed << ": DistMax/layer/i_hit_max/ntot="
204  << DistMax << "/" << layer << "/" << i_hit_max << "/" << v_mdtLayerHits[chamber][0].ntot);
205 
206  // break conditions
207  if(layer == 999999) break;
208  if(v_mdtLayerHits[chamber][layer].ndigi==1) break;
209 
210  double Mednew = (v_mdtLayerHits[chamber][0].ResSum - Residual)/(v_mdtLayerHits[chamber][0].ntot - 1);
211  double Delta = 2.*std::abs((ResMed - Mednew)/(ResMed + Mednew));
212  ATH_MSG_DEBUG("Mednew/Delta/DeltaMin=" << Mednew << "/" << Delta << "/" << DeltaMin);
213  if(Delta<=DeltaMin) break;
214 
215  // if not, delete the maxRes and continue;
216  v_mdtLayerHits[chamber][0].ResSum -= Residual;
217  v_mdtLayerHits[chamber][0].ntot--;
218  v_mdtLayerHits[chamber][layer].ndigi--;
219  mdtHits[i_hit_max].isOutlier = 2;
220  }
221 
222  ATH_MSG_DEBUG("choosing one at each layer, and record it in segment...");
223  TrigL2MuonSA::MdtHits& mdtSegment {trackPattern.mdtSegments[chamber]};
224  mdtSegment.clear();
225 
226  double phi0 {0}; // phi at the first layer to be used
227  for(unsigned int i_layer=0; i_layer<=i_layer_max; i_layer++) {
228 
229  ATH_MSG_DEBUG("i_layer=" << i_layer << ": ndigi=" << v_mdtLayerHits[chamber][i_layer].ndigi);
230 
231  // choose one at each layer
232  while( v_mdtLayerHits[chamber][i_layer].ndigi>=2 ) {
233  double ResMax = 0.;
234  unsigned int i_hit_max = 999999;
235  for(unsigned int idigi=0;idigi<v_mdtLayerHits[chamber][i_layer].ndigi_all;idigi++) {
236  unsigned int i_hit = v_mdtLayerHits[chamber][i_layer].indexes[idigi];
237  if( mdtHits[i_hit].isOutlier > 0 ) continue;
238  if( std::abs(mdtHits[i_hit].Residual) >= ResMax ) {
239  ResMax = std::abs(mdtHits[i_hit].Residual);
240  i_hit_max = i_hit;
241  }
242  }
243  ATH_MSG_DEBUG("ResMax=" << ResMax << ": i_hit_max=" << i_hit_max);
244  if( i_hit_max == 999999 ) break;
245  v_mdtLayerHits[chamber][0].ResSum -= ResMax;
246  v_mdtLayerHits[chamber][0].ntot--;
247  v_mdtLayerHits[chamber][i_layer].ndigi--;
248  mdtHits[i_hit_max].isOutlier = 1;
249  }
250 
251  // record it into segement
252  for(unsigned int i_digi=0; i_digi< v_mdtLayerHits[chamber][i_layer].ndigi_all; i_digi++) {
253  unsigned int i_hit = v_mdtLayerHits[chamber][i_layer].indexes[i_digi];
254  if (i_layer==0) phi0 = mdtHits[i_hit].cPhi0;
255  doMdtCalibration(ctx, mdtHits[i_hit], trigger_phi, phi0, muonRoad.isEndcap);
256  if (mdtHits[i_hit].isOutlier > 1) continue;
257  mdtSegment.push_back(mdtHits[i_hit]);
258  }
259  }
260  ATH_MSG_DEBUG("nr of hits in segment=" << mdtSegment.size());
261  } // end loop on stations.
262 
263  v_trackPatterns.push_back(trackPattern);
264 
265  return StatusCode::SUCCESS;
266 }
267 
268 // --------------------------------------------------------------------------------
269 // --------------------------------------------------------------------------------
270 
271 double TrigL2MuonSA::MuFastPatternFinder::calc_residual(double aw,double bw,double x,double y) const
272 {
273  const double ZERO_LIMIT = 1e-4;
274  if( std::abs(aw) < ZERO_LIMIT ) return y-bw;
275  double ia = 1/aw;
276  double iaq = ia*ia;
277  double dz = x - (y-bw)*ia;
278  return dz/std::sqrt(1.+iaq);
279 }
MdtRegionDefiner.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
TrigL2MuonSA::MdtHitData::DriftSigma
double DriftSigma
Definition: MdtData.h:39
TrigL2MuonSA::MuFastPatternFinder::doMdtCalibration
void doMdtCalibration(const EventContext &ctx, TrigL2MuonSA::MdtHitData &mdtHit, double track_phi, double phi0, bool isEndcap) const
Definition: MuFastPatternFinder.cxx:20
TrigL2MuonSA::MdtHitData::readEle
const MuonGM::MdtReadoutElement * readEle
Definition: MdtData.h:50
MuFastPatternFinder.h
MdtCalibInput
Definition: MdtCalibInput.h:34
TrigL2MuonSA::MuonRoad::phiMiddle
double phiMiddle
Definition: MuonRoad.h:81
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
calibdata.chamber
chamber
Definition: calibdata.py:31
TrigL2MuonSA::MdtHitData::StationPhi
int StationPhi
Definition: MdtData.h:21
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
TrigL2MuonSA::MdtHits
std::vector< MdtHitData > MdtHits
Definition: MdtData.h:56
TrigL2MuonSA::MdtHitData
Definition: MdtData.h:18
TrigL2MuonSA::MuonRoad::bw
double bw[N_STATION][N_SECTOR]
Definition: MuonRoad.h:84
TrigL2MuonSA::TrackPattern::s_address
int s_address
Definition: TrackData.h:72
TrigL2MuonSA::MuonRoad::isEndcap
bool isEndcap
Definition: MuonRoad.h:74
TrigL2MuonSA::MdtHitData::R
double R
Definition: MdtData.h:36
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
TrigL2MuonSA::MdtHitData::name
unsigned int name
Definition: MdtData.h:19
x
#define x
TrigL2MuonSA::MdtHitData::StationEta
int StationEta
Definition: MdtData.h:20
Identifier::is_valid
bool is_valid() const
Check if id is in a valid state.
python.RingerConstants.Layer
Layer
Definition: RingerConstants.py:42
TrigL2MuonSA::MuonRoad::aw
double aw[N_STATION][N_SECTOR]
Definition: MuonRoad.h:83
TrigL2MuonSA::MdtHitData::Multilayer
int Multilayer
Definition: MdtData.h:22
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
TrigL2MuonSA::TrackPattern
Definition: TrackData.h:16
TrigL2MuonSA::MdtHitData::TubeLayer
int TubeLayer
Definition: MdtData.h:24
TrigL2MuonSA::MuFastPatternFinder::calc_residual
double calc_residual(double aw, double bw, double x, double y) const
Definition: MuFastPatternFinder.cxx:271
TrigL2MuonSA::MdtLayerHits
Definition: MuFastPatternFinder.h:26
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
TrigL2MuonSA::MdtHitData::Id
Identifier Id
Definition: MdtData.h:49
lumiFormat.i
int i
Definition: lumiFormat.py:85
TrigL2MuonSA::MuonRoad
Definition: MuonRoad.h:20
xAOD::L2MuonParameters::MaxChamber
@ MaxChamber
Number of measurement point definitions.
Definition: TrigMuonDefs.h:27
TrigL2MuonSA::MuFastPatternFinder::m_nswPatternFinder
ToolHandle< NswPatternFinder > m_nswPatternFinder
Definition: MuFastPatternFinder.h:65
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
Residual
Residual is a class that stores the residual, error, and type of residual.
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
TrigL2MuonSA::MuonRoad::MDT_sector_trigger
int MDT_sector_trigger
Definition: MuonRoad.h:88
TrigL2MuonSA::MuFastPatternFinder::initialize
virtual StatusCode initialize() override
Definition: MuFastPatternFinder.cxx:10
TrigL2MuonSA::MdtHitData::DriftTime
double DriftTime
Definition: MdtData.h:37
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TrigL2MuonSA::MdtHitData::Adc
uint16_t Adc
Definition: MdtData.h:41
MdtCalibInput::setClosestApproach
void setClosestApproach(const Amg::Vector3D &approach)
Sets the closest approach.
Definition: MdtCalibInput.cxx:106
TrigL2MuonSA::MuFastPatternFinder::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: MuFastPatternFinder.h:70
MdtCalibOutput
Definition: MdtCalibOutput.h:10
ZERO_LIMIT
const float ZERO_LIMIT
Definition: VP1TriggerHandleL2.cxx:37
TrigL2MuonSA::MdtRegionDefiner::find_station_sector
static void find_station_sector(const std::string &name, int phi, bool &endcap, int &chamber, int &sector)
Definition: MdtRegionDefiner.cxx:329
Monitored::Y
@ Y
Definition: HistogramFillerUtils.h:24
TrigL2MuonSA::MdtHitData::Tube
int Tube
Definition: MdtData.h:25
TrigL2MuonSA::MdtHitData::DriftSpace
double DriftSpace
Definition: MdtData.h:38
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
TrigMuonDefs.h
MdtCalibOutput::driftRadiusUncert
double driftRadiusUncert() const
Returns the uncertainty on the drift radius.
Definition: MdtCalibOutput.cxx:20
TrigL2MuonSA::MuonRoad::LargeSmall
int LargeSmall
Definition: MuonRoad.h:79
TrigL2MuonSA::MuFastPatternFinder::findPatterns
StatusCode findPatterns(const EventContext &ctx, const TrigL2MuonSA::MuonRoad &muonRoad, TrigL2MuonSA::MdtHits &mdtHits, std::vector< TrigL2MuonSA::TrackPattern > &v_trackPatterns) const
Definition: MuFastPatternFinder.cxx:99
TrigL2MuonSA::MmHits
std::vector< MmHitData > MmHits
Definition: MmData.h:47
y
#define y
TrigL2MuonSA::MuonRoad::rWidth
double rWidth[N_STATION][N_LAYER]
Definition: MuonRoad.h:86
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
TrigL2MuonSA::StgcHits
std::vector< StgcHitData > StgcHits
Definition: StgcData.h:49
DEBUG
#define DEBUG
Definition: page_access.h:11
TrigL2MuonSA::MdtHitData::Z
double Z
Definition: MdtData.h:35
GeoPrimitivesToStringConverter.h
TrigL2MuonSA::MuonRoad::Special
int Special
Definition: MuonRoad.h:80
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
MdtCalibOutput::driftRadius
double driftRadius() const
Returns the drift radius of the calibrated object.
Definition: MdtCalibOutput.cxx:19
TrigL2MuonSA::TrackPattern::mdtSegments
TrigL2MuonSA::MdtHits mdtSegments[s_NCHAMBER]
Definition: TrackData.h:57
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Sun Sep 28 2025 21:15:11 for ATLAS Offline Software by   doxygen 1.8.18