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