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