ATLAS Offline Software
NswStationFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "NswStationFitter.h"
6 #include "xAODTrigMuon/TrigMuonDefs.h"
7 
8 #include "RecMuonRoIUtils.h"
9 #include <cmath>
10 #include "AthenaBaseComps/AthMsgStreamMacros.h"
11 
12 #include <vector>
13 #include <array>
14 
15 namespace{
16  constexpr double ZERO_LIMIT = 1.e-9;
17 }
18 
19 static const double TanM1p5 = std::tan(-1.5/360.*2.*M_PI);
20 static const double TanP1p5 = std::tan(1.5/360.*2.*M_PI);
21 static const double CosM1p5 = std::cos(-1.5/360.*2.*M_PI);
22 static const double CosP1p5 = std::cos(1.5/360.*2.*M_PI);
23 static const double SinM1p5 = std::sin(-1.5/360.*2.*M_PI);
24 static const double SinP1p5 = std::sin(1.5/360.*2.*M_PI);
25 
26 TrigL2MuonSA::NswStationFitter::NswStationFitter(const std::string& type,
27  const std::string& name,
28  const IInterface* parent):
29  AthAlgTool(type,name,parent)
30 {
31 }
32 
33 StatusCode TrigL2MuonSA::NswStationFitter::superPointFitter(const TrigRoiDescriptor* p_roids,
34  TrigL2MuonSA::TrackPattern& trackPattern,
35  TrigL2MuonSA::StgcHits& stgcHits,
36  TrigL2MuonSA::MmHits& mmHits) const
37 {
38 
39  ATH_MSG_DEBUG("NswStationFitter::findSuperPoints() was called.");
40 
41  // selection for sTGC hits, RoI matching based
42  ATH_CHECK( selectStgcHits(p_roids,trackPattern.stgcSegment) );
43  ATH_CHECK( selectMmHits(p_roids,trackPattern.mmSegment) );
44 
45  ATH_CHECK( findStgcHitsInSegment(stgcHits) );
46  ATH_CHECK( findMmHitsInSegment(mmHits) );
47 
48  bool isLargeStgc = false;
49  bool isSmallStgc = false;
50  if(stgcHits.size() != 0){
51  for(unsigned int iHit = 0; iHit < stgcHits.size(); iHit++){
52  if(stgcHits.at(iHit).isOutlier == 0){
53  if(stgcHits.at(iHit).stationName == 58) isLargeStgc = true;
54  else if(stgcHits.at(iHit).stationName == 57) isSmallStgc = true;
55  if(isLargeStgc && isSmallStgc) continue;
56  }
57  }
58  }
59  if(mmHits.size() != 0){
60  for(unsigned int iHit = 0; iHit < mmHits.size(); iHit++){
61  if(mmHits.at(iHit).isOutlier == 0){
62  if(isLargeStgc){
63  if(mmHits.at(iHit).stationName == 55) mmHits.at(iHit).isOutlier = 1;
64  } else if(isSmallStgc) {
65  if(mmHits.at(iHit).stationName == 56) mmHits.at(iHit).isOutlier = 1;
66  }
67  }
68  }
69  }
70  ATH_CHECK( MakeSegment(trackPattern,stgcHits) );
71  ATH_CHECK( MakeSegment(trackPattern,mmHits) );
72 
73  ATH_MSG_DEBUG("Number of sTGC and MM hits for SPs " << trackPattern.stgcSegment.size() << " " << trackPattern.mmSegment.size());
74  if (trackPattern.stgcSegment.size() < 9 && trackPattern.mmSegment.size() < 6) {
75  ATH_MSG_DEBUG("Number of sTGC and MM hits for SPs is small " << trackPattern.stgcSegment.size() + trackPattern.mmSegment.size());
76  }
77  else {
78  ATH_CHECK( calcMergedHit(trackPattern) );
79  }
80 
81 
82  return StatusCode::SUCCESS;
83 
84 }
85 
86 StatusCode TrigL2MuonSA::NswStationFitter::selectStgcHits(const TrigRoiDescriptor* p_roids,
87  TrigL2MuonSA::StgcHits& stgcHits) const
88 {
89 
90  TrigL2MuonSA::StgcHits selectedStgcHits;
91  selectedStgcHits.clear();
92 
93  // define region where RoI is near
94  double etaMin = p_roids->eta() - 1.;
95  double etaMax = p_roids->eta() + 1.;
96  double phiMin = p_roids->phi() - 1.;
97  double phiMax = p_roids->phi() + 1.;
98  if( phiMin > M_PI ) phiMin -= 2*M_PI;
99  if( phiMax > M_PI ) phiMax -= 2*M_PI;
100  if( phiMin < -1.*M_PI ) phiMin += 2*M_PI;
101  if( phiMax < -1.*M_PI ) phiMax += 2*M_PI;
102  // boolian to check if sTGC hits are near RoI
103  bool inRoiEta, inRoiPhi, inRoi;
104 
105  // loop over sTGC digits.
106  unsigned int iHit;
107  if (stgcHits.size()>0) {
108  for (iHit = 0; iHit < stgcHits.size(); iHit++){
109  // Get the digit point.
110  TrigL2MuonSA::StgcHitData& hit = stgcHits[iHit];
111 
112  // check if sTGC hits are near RoI
113  inRoiEta = false; inRoiPhi = false; inRoi = false;
114  if ( etaMin<=hit.eta && hit.eta<=etaMax ) inRoiEta = true;
115  if ( phiMin<=phiMax && phiMin<=hit.phi && hit.phi<=phiMax ) inRoiPhi= true;
116  if ( phiMin>phiMax && (phiMin<=hit.phi || hit.phi<=phiMax)) inRoiPhi= true;
117  if ( inRoiEta && inRoiPhi ) inRoi = true;
118  ATH_MSG_DEBUG("sTGC hits eta = "<<hit.eta<<", phi = "<<hit.phi<<", r = "<<hit.r<<", z = "<<hit.z<<", stationEta = "<<hit.stationEta<<", stationPhi = "<<hit.stationPhi<<", channelType = "<<hit.channelType<<", stationName = "<<hit.stationName<<",layerNumber ="<<hit.layerNumber<<", matched RoI? = "<<inRoi);
119  // pushback if the hit is near RoI
120  if (!inRoi){
121  continue;
122  }
123  selectedStgcHits.push_back(hit);
124  }
125  }
126 
127  stgcHits.clear();
128  stgcHits = selectedStgcHits;
129 
130  return StatusCode::SUCCESS;
131 
132 }
133 
134 StatusCode TrigL2MuonSA::NswStationFitter::selectMmHits(const TrigRoiDescriptor* p_roids,
135  TrigL2MuonSA::MmHits& mmHits) const
136 {
137 
138  TrigL2MuonSA::MmHits selectedMmHits;
139  selectedMmHits.clear();
140 
141  // define region where RoI is near
142  double etaMin = p_roids->eta() - 1.;
143  double etaMax = p_roids->eta() + 1.;
144  double phiMin = p_roids->phi() - 1.;
145  double phiMax = p_roids->phi() + 1.;
146  if( phiMin > M_PI ) phiMin -= 2*M_PI;
147  if( phiMax > M_PI ) phiMax -= 2*M_PI;
148  if( phiMin < -1.*M_PI ) phiMin += 2*M_PI;
149  if( phiMax < -1.*M_PI ) phiMax += 2*M_PI;
150  // boolian to check if sTGC hits are near RoI
151  bool inRoiEta, inRoiPhi, inRoi;
152 
153  // loop over MM digits.
154  unsigned int iHit;
155  if (mmHits.size()>0) {
156  for (iHit = 0; iHit < mmHits.size(); iHit++){
157 
158  // Get the digit point.
159  TrigL2MuonSA::MmHitData& hit = mmHits[iHit];
160 
161  // check if MM hits are near RoI
162  inRoiEta = false; inRoiPhi = false; inRoi = false;
163  if ( etaMin<=hit.eta && hit.eta<=etaMax ) inRoiEta = true;
164  if ( phiMin<=phiMax && phiMin<=hit.phi && hit.phi<=phiMax ) inRoiPhi= true;
165  if ( phiMin>phiMax && (phiMin<=hit.phi || hit.phi<=phiMax)) inRoiPhi= true;
166  if ( inRoiEta && inRoiPhi ) inRoi = true;
167  ATH_MSG_DEBUG("MM hits eta = "<<hit.eta<<", phi = "<<hit.phi<<", r = "<<hit.r<<", z = "<<hit.z<<", stationEta = "<<hit.stationEta<<", stationPhi = "<<hit.stationPhi<<", stationName = "<<hit.stationName<<", matched RoI? = "<<inRoi);
168 
169  // pushback if the hit is near RoI
170  if (!inRoi){
171  continue;
172  }
173  selectedMmHits.push_back(hit);
174  }
175  }
176 
177  mmHits.clear();
178  mmHits = selectedMmHits;
179 
180  return StatusCode::SUCCESS;
181 
182 }
183 
184 StatusCode TrigL2MuonSA::NswStationFitter::calcWeightedSumHit(TrigL2MuonSA::TrackPattern& trackPattern) const
185 {
186 
187  TrigL2MuonSA::StgcHits stgcHits = trackPattern.stgcSegment;
188  TrigL2MuonSA::MmHits mmHits = trackPattern.mmSegment;
189 
190  if(stgcHits.size()==0 && mmHits.size()==0)
191  return StatusCode::SUCCESS;
192 
193  //superpoint information
194  float rWeightedCenter=0.,zWeightedCenter=0.,phiWeightedCenter=0.;
195  int totNumRWeightedCenter=0.,totNumZWeightedCenter=0.,totNumPhiWeightedCenter=0.;
196  float localPhiCenter=0., localPhi;
197 
198  // loop over sTGC digits.
199  unsigned int iHit;
200  if (stgcHits.size()>0) {
201  for (iHit = 0; iHit < stgcHits.size(); iHit++){
202 
203  // Get the digit point.
204  TrigL2MuonSA::StgcHitData& hit = stgcHits[iHit];
205 
206  if (iHit==0 && hit.stationPhi<=5) localPhiCenter = 0.25 * M_PI * ((float)hit.stationPhi-1.);
207  if (iHit==0 && hit.stationPhi> 5) localPhiCenter = 0.25 * M_PI * ((float)hit.stationPhi-9.);
208  if (hit.stationName == 57){
209  localPhiCenter += M_PI/8.; // small sTGC sectors
210  if (hit.stationPhi == 5) localPhiCenter -= 2 * M_PI;
211  }
212 
213  localPhi = hit.phi - localPhiCenter;
214  if (localPhi > M_PI) localPhi -= 2.0*M_PI;
215  if (localPhi < -1.*M_PI) localPhi += 2.0*M_PI;
216  // calculate Weighted Center
217  if (hit.channelType == 1) { // strip
218  rWeightedCenter += hit.r;
219  totNumRWeightedCenter += 1;
220  }
221  if (hit.channelType == 2) { // wire
222  phiWeightedCenter += localPhi;
223  totNumPhiWeightedCenter += 1;
224  }
225  zWeightedCenter += hit.z;
226  totNumZWeightedCenter += 1;
227  }
228  }
229 
230  // loop over MM digits.
231  if (mmHits.size()>0) {
232  for (iHit = 0; iHit < mmHits.size(); iHit++){
233 
234  // Get the digit point.
235  TrigL2MuonSA::MmHitData& hit = mmHits[iHit];
236 
237  if (hit.stationName == 55){
238  localPhiCenter += M_PI/8.; // small MM sectors
239  if (hit.stationPhi == 5) localPhiCenter -= 2 * M_PI;
240  }
241 
242  localPhi = hit.phi - localPhiCenter;
243  if (localPhi > M_PI) localPhi -= 2.0*M_PI;
244  if (localPhi < -1.*M_PI) localPhi += 2.0*M_PI;
245  }
246  }
247 
248  // calculate Weighted Center
249  if (totNumRWeightedCenter!=0) rWeightedCenter /= totNumRWeightedCenter;
250  if (totNumPhiWeightedCenter!=0) phiWeightedCenter /= totNumPhiWeightedCenter;
251  if (totNumZWeightedCenter!=0) zWeightedCenter /= totNumZWeightedCenter;
252 
253  // store superpoint info in TrackData
254  xAOD::L2MuonParameters::Chamber inner = xAOD::L2MuonParameters::Chamber::EndcapInner;
255  TrigL2MuonSA::SuperPoint* superPoint = &(trackPattern.superPoints[inner]);
256  superPoint->R = rWeightedCenter;
257  superPoint->Phim = phiWeightedCenter+localPhiCenter;
258  superPoint->Z = zWeightedCenter;
259  superPoint->Npoint = totNumZWeightedCenter;
260  if (zWeightedCenter!=0) superPoint->Alin = rWeightedCenter/zWeightedCenter;
261  superPoint->Blin = 0.;
262 
263  ATH_MSG_DEBUG("Nsw Super Point r/phi/z/slope = "<<superPoint->R<<"/"<<superPoint->Phim<<"/"<<superPoint->Z<<"/"<<superPoint->Alin);
264 
265  return StatusCode::SUCCESS;
266 
267 }
268 
269 StatusCode TrigL2MuonSA::NswStationFitter::findStgcHitsInSegment(TrigL2MuonSA::StgcHits& stgcHits) const
270 {
271  if(stgcHits.size() == 0) return StatusCode::SUCCESS;
272  int hitsInRoad = 0;
273  for(unsigned int iHit = 0; iHit < stgcHits.size(); iHit++){
274  if(stgcHits.at(iHit).isOutlier == 0){
275  hitsInRoad++;
276  stgcHits.at(iHit).isOutlier = 1;
277  }
278  }
279  if(hitsInRoad == 0) return StatusCode::SUCCESS;
280 
281  if(hitsInRoad < 9) {
282  ATH_MSG_DEBUG("Number of STGC hits is too small, at least 9 hits required : "<<hitsInRoad<<" hits");
283  return StatusCode::SUCCESS;
284  } else if(hitsInRoad > 100) {
285  ATH_MSG_WARNING("Number of STGC hits is too large, at most 100 hits allowed : "<<hitsInRoad<<" hits");
286  return StatusCode::SUCCESS;
287  }
288 
289  // cassifyDataEachLayer
290  std::array<std::vector<int>, 8> strHitIdByLayer;
291  std::array<std::vector<int>, 8> wireHitIdByLayer;
292  for(unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit){
293  if(stgcHits.at(iHit).isOutlier != 1) continue;
294  int layerNumber = stgcHits.at(iHit).layerNumber;
295  if (layerNumber > 7) {
296  ATH_MSG_WARNING("STGC hit layer number > 7");
297  continue;
298  }
299  if(stgcHits.at(iHit).channelType == 1){
300  strHitIdByLayer[layerNumber].push_back(iHit);
301  }else if(stgcHits.at(iHit).channelType == 2){
302  wireHitIdByLayer[layerNumber].push_back(iHit);
303  }
304  }
305  ATH_MSG_DEBUG("@@STGC@@ strip Nhits " << strHitIdByLayer[0].size()
306  << " " << strHitIdByLayer[1].size()
307  << " " << strHitIdByLayer[2].size()
308  << " " << strHitIdByLayer[3].size()
309  << " " << strHitIdByLayer[4].size()
310  << " " << strHitIdByLayer[5].size()
311  << " " << strHitIdByLayer[6].size()
312  << " " << strHitIdByLayer[7].size());
313  ATH_MSG_DEBUG("@@STGC@@ wire Nhits " << wireHitIdByLayer[0].size()
314  << " " << wireHitIdByLayer[1].size()
315  << " " << wireHitIdByLayer[2].size()
316  << " " << wireHitIdByLayer[3].size()
317  << " " << wireHitIdByLayer[4].size()
318  << " " << wireHitIdByLayer[5].size()
319  << " " << wireHitIdByLayer[6].size()
320  << " " << wireHitIdByLayer[7].size());
321 
322  std::vector< std::array<int, 8> > strHitIds; // Candidates' sets of strip hit ids in 8 layers
323  findSetOfStgcHitIds(stgcHits, strHitIdByLayer, strHitIds);
324  std::vector< std::array<int, 8> > wireHitIds; // Candidates' sets of wire hit ids in 8 layers
325  findSetOfStgcHitIds(stgcHits, wireHitIdByLayer, wireHitIds);
326  ATH_MSG_DEBUG("@@STGC@@ strip wire " << strHitIds.size() << " " << wireHitIds.size());
327 
328  bool isLargeStrip = false;
329  bool isSmallStrip = false;
330  for (unsigned int iHit = 0; iHit < strHitIds.size(); ++iHit) {
331  std::array<int, 8> hitIds = strHitIds.at(iHit);
332  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
333  if (hitIds[iLayer] != -1) {
334  if(stgcHits.at(hitIds[iLayer]).stationName == 58) isLargeStrip = true;
335  else if(stgcHits.at(hitIds[iLayer]).stationName == 57) isSmallStrip = true;
336  stgcHits.at(hitIds[iLayer]).isOutlier = 0;
337  }
338  }
339  }
340  for (unsigned int iHit = 0; iHit < wireHitIds.size(); ++iHit) {
341  std::array<int, 8> hitIds = wireHitIds.at(iHit);
342  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
343  if (hitIds[iLayer] != -1) {
344  if(isLargeStrip){
345  if(stgcHits.at(hitIds[iLayer]).stationName == 58){
346  stgcHits.at(hitIds[iLayer]).isOutlier = 0;
347  }
348  }
349  else if(isSmallStrip){
350  if(stgcHits.at(hitIds[iLayer]).stationName == 57){
351  stgcHits.at(hitIds[iLayer]).isOutlier = 0;
352  }
353  }
354  }
355  }
356  }
357  return StatusCode::SUCCESS;
358 
359 }
360 
361 void TrigL2MuonSA::NswStationFitter::findSetOfStgcHitIds(TrigL2MuonSA::StgcHits& stgcHits,
362  std::array<std::vector<int>,8> hitIdByLayer,
363  std::vector<std::array<int, 8>>& hitIdsCandidate) const
364 {
365  double NSWCenterZ = 7526.329;
366  int side = 0;
367 
368  bool isStrip = 0;
369  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
370  if ( hitIdByLayer[iLayer].size() > 0) {
371  side = std::abs(stgcHits.at(hitIdByLayer[iLayer].at(0)).z)/stgcHits.at(hitIdByLayer[iLayer].at(0)).z;
372  if ( stgcHits.at(hitIdByLayer[iLayer].at(0)).channelType == 1 ){
373  isStrip = 1;
374  break;
375  }
376  }
377  }
378  NSWCenterZ = NSWCenterZ * side;
379 
380  std::array<std::vector<unsigned long int>,4> hitIdsInTwo;
381  std::array<std::vector<double>,4> slopeInTwo;
382  std::array<std::vector<double>,4> interceptInTwo;
383 
384  // Loop over pairs of the i-th and the (i+4)-th layers
385  for(unsigned int iPair = 0; iPair < 4; ++iPair){
386  unsigned int nHitsInInner = hitIdByLayer[iPair].size();
387  unsigned int nHitsInOuter = hitIdByLayer[iPair+4].size();
388  if ( nHitsInInner > 0xffff-1 || nHitsInOuter > 0xffff-1) {
389  ATH_MSG_WARNING("Number of Stgc hits in layers exceeds the limit of (2^16 - 1) : Number of Stgc hits in "<<iPair<<"th layer = "<< nHitsInInner
390  <<", Number of Stgc hits in "<<iPair+4<<"th layer = "<<nHitsInOuter);
391  ATH_MSG_WARNING("Number of Stgc hits is limitted to (2^16 - 1) and hits with id more than (2^16 -1) will be trancated.");
392  if (nHitsInInner > 0xffff-1) {nHitsInInner = 0xffff-1;}
393  if (nHitsInOuter > 0xffff-1) {nHitsInOuter = 0xffff-1;}
394  }
395 
396  bool foundCounterparts[256] = {};
397  // Loop over hits in the i-th layer
398  for(unsigned int iHit = 0; iHit < nHitsInInner; ++iHit){
399  bool foundCounterpart = 0;
400 
401  double z[2] = {};
402  double r[2] = {};
403 
404  int iHitId = hitIdByLayer[iPair].at(iHit);
405  if(isStrip){
406  r[0] = stgcHits.at(iHitId).r;
407  z[0] = stgcHits.at(iHitId).z;
408  } else {
409  double localPhiCenter;
410  if (stgcHits.at(iHitId).stationPhi<=5) {
411  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHitId).stationPhi-1.);
412  } else {
413  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHitId).stationPhi-9.);
414  }
415 
416  if (stgcHits.at(iHitId).stationName == 57){
417  localPhiCenter += M_PI/8.; // small sTGC sectors
418  if (stgcHits.at(iHitId).stationPhi == 5) localPhiCenter -= 2. * M_PI;
419  }
420 
421  double phiProj = stgcHits.at(iHitId).phi - localPhiCenter;
422  if (phiProj > M_PI) phiProj -= 2.0*M_PI;
423  if (phiProj < -1.*M_PI) phiProj += 2.0*M_PI;
424  r[0] = stgcHits.at(iHitId).r;
425  z[0] = phiProj;
426  }
427 
428  // Loop over hits in the (i+4)-th layer
429  for(unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
430  int jHitId = hitIdByLayer[iPair+4].at(jHit);
431  double slope, intercept;
432  if(isStrip) {
433  r[1] = stgcHits.at(jHitId).r;
434  z[1] = stgcHits.at(jHitId).z;
435  slope = (r[1] - r[0]) / (z[1] - z[0]);
436  intercept = slope*(0. - z[0]) + r[0];
437  // select pairs whose slops in limited regions
438  if(std::abs(slope) < 0.14 || std::abs(slope) > 0.6 || std::abs(intercept) > 300.) continue;
439  } else {
440  double localPhiCenter;
441  if (stgcHits.at(jHitId).stationPhi<=5) {
442  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(jHitId).stationPhi-1.);
443  } else {
444  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(jHitId).stationPhi-9.);
445  }
446 
447  if (stgcHits.at(jHitId).stationName == 57){
448  localPhiCenter += M_PI/8.; // small sTGC sectors
449  if (stgcHits.at(jHitId).stationPhi == 5) localPhiCenter -= 2 * M_PI;
450  }
451 
452  double phiProj = stgcHits.at(jHitId).phi - localPhiCenter;
453  if (phiProj > M_PI) phiProj -= 2.0*M_PI;
454  if (phiProj < -1.*M_PI) phiProj += 2.0*M_PI;
455  r[1] = stgcHits.at(jHitId).r;
456  z[1] = phiProj;
457  slope = (z[0]+z[1])/2.;
458  intercept = (r[0]+r[1])/2.;
459  if(std::abs(r[0]*std::sin(z[0]) - r[1]*std::sin(z[1])) > 300.) continue;
460  }
461 
462  unsigned int encodedIds = (iHitId<<16) + jHitId;
463  hitIdsInTwo[iPair].push_back(encodedIds);
464  slopeInTwo[iPair].push_back(slope);
465  interceptInTwo[iPair].push_back(intercept);
466 
467  foundCounterpart = 1;
468  foundCounterparts[jHit] = 1;
469  }//end of jHit in the (i+4)-th layer
470  if(!foundCounterpart){ // in case of no counterpart in the (i+4)-th layer
471  unsigned int encodedIds = (iHitId<<16) + 0xffff; // fill all bits with 1 for hit id for the layer with no hit
472  hitIdsInTwo[iPair].push_back(encodedIds);
473  if(isStrip) {
474  slopeInTwo[iPair].push_back(r[0]/z[0]);
475  interceptInTwo[iPair].push_back(0.);
476  } else {
477  slopeInTwo[iPair].push_back(z[0]);
478  interceptInTwo[iPair].push_back(r[0]);
479  }
480  }
481  }//end of iHit in the i-th layer
482  // Loop over hits in the (i+4)-th layer
483  for(unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
484  if (!foundCounterparts[jHit]) {
485  int jHitId = hitIdByLayer[iPair+4].at(jHit);
486  unsigned int encodedIds = 0xffff0000 + jHitId; // fill all bits with 1 for hit id for the layer with no hit
487  hitIdsInTwo[iPair].push_back(encodedIds);
488  if(isStrip) {
489  slopeInTwo[iPair].push_back(stgcHits.at(jHitId).r/stgcHits.at(jHitId).z);
490  interceptInTwo[iPair].push_back(0.);
491  } else {
492  double localPhiCenter;
493  if (stgcHits.at(jHitId).stationPhi<=5) {
494  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(jHitId).stationPhi-1.);
495  } else {
496  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(jHitId).stationPhi-9.);
497  }
498  if (stgcHits.at(jHitId).stationName == 57){
499  localPhiCenter += M_PI/8.; // small sTGC sectors
500  if (stgcHits.at(jHitId).stationPhi == 5) localPhiCenter -= 2 * M_PI;
501  }
502 
503  double phiProj = stgcHits.at(jHitId).phi - localPhiCenter;
504  if (phiProj > M_PI) phiProj -= 2.0*M_PI;
505  if (phiProj < -1.*M_PI) phiProj += 2.0*M_PI;
506  slopeInTwo[iPair].push_back(phiProj);
507  interceptInTwo[iPair].push_back(stgcHits.at(jHitId).r);
508  }
509  }
510  }
511  }//end of pair loop
512 
513  ATH_MSG_DEBUG("@@STGC@@ isStrip= " << isStrip << " Npairs " << hitIdsInTwo[0].size() << " " << hitIdsInTwo[1].size() << " " << hitIdsInTwo[2].size() << " " << hitIdsInTwo[3].size());
514  for (unsigned int iLayer = 0; iLayer < 4; ++ iLayer) {
515  for (unsigned int iPair = 0; iPair < slopeInTwo[iLayer].size(); ++iPair) {
516  ATH_MSG_DEBUG("@@STGC@@ pair fit isStrip= " << isStrip << " slope= " << slopeInTwo[iLayer].at(iPair) << " intercept= " << interceptInTwo[iLayer].at(iPair));
517  }
518  }
519 
520  std::array<std::vector<unsigned long int>,2> hitIdsInFour;
521  std::array<std::vector<double>,2> slopeInFour;
522  std::array<std::vector<double>,2> interceptInFour;
523  for(unsigned int iQuad = 0; iQuad < 2; ++iQuad){
524  unsigned int nPairsInInner = hitIdsInTwo[iQuad].size();
525  unsigned int nPairsInOuter = hitIdsInTwo[iQuad+2].size();
526 
527  bool foundCounterparts[0xffff] = {};
528  for(unsigned int iPair = 0; iPair < nPairsInInner; ++iPair){
529  bool foundCounterpart = 0;
530  double slope[2];
531  double intercept[2];
532  slope[0] = slopeInTwo[iQuad].at(iPair);
533  intercept[0] = interceptInTwo[iQuad].at(iPair);
534 
535  for(unsigned int jPair = 0; jPair < nPairsInOuter; ++jPair){
536  unsigned int ihitIds = hitIdsInTwo[iQuad].at(iPair);
537  unsigned int jhitIds = hitIdsInTwo[iQuad+2].at(jPair);
538  if ( !(((ihitIds>>16 & 0xffff) != 0xffff || (ihitIds & 0xffff) != 0xffff) &&
539  ((jhitIds>>16 & 0xffff) != 0xffff || (jhitIds & 0xffff) != 0xffff )) ) continue; // require at least 2 hits in 4 layers
540 
541  slope[1] = slopeInTwo[iQuad+2].at(jPair);
542  intercept[1] = interceptInTwo[iQuad+2].at(jPair);
543 
544  if (isStrip) {
545  double spR0 = slope[0] * NSWCenterZ + intercept[0];
546  double spR1 = slope[1] * NSWCenterZ + intercept[1];
547  if(std::abs(spR1 - spR0) > 50.) continue; //OPTIMIZE ME!!!
548  } else {
549  if(std::abs(intercept[0]*std::sin(slope[0]) - intercept[1]*std::sin(slope[1])) > 100.) continue;
550  }
551 
552  foundCounterpart = 1;
553  foundCounterparts[jPair] = 1;
554 
555  unsigned long int encodedIds = (hitIdsInTwo[iQuad].at(iPair) << 32 ) + hitIdsInTwo[iQuad+2].at(jPair);
556  hitIdsInFour[iQuad].push_back(encodedIds);
557  slopeInFour[iQuad].push_back((slope[1] + slope[0])/2.);
558  interceptInFour[iQuad].push_back((intercept[1] + intercept[0])/2.);
559  }// end of iPair of the i-th and (i+4)-th layers, i=0,2
560 
561  if(foundCounterpart) continue; // in case of no counterpart of pairs in the inner layers
562  if((hitIdsInTwo[iQuad].at(iPair)>>16 & 0xffff) == 0xffff || (hitIdsInTwo[iQuad].at(iPair) & 0xffff) == 0xffff) continue;
563 
564  unsigned long int encodedIds = (hitIdsInTwo[iQuad].at(iPair) << 32 ) + 0xffffffff;
565  hitIdsInFour[iQuad].push_back(encodedIds);
566  slopeInFour[iQuad].push_back(slope[0]);
567  interceptInFour[iQuad].push_back(intercept[0]);
568  }// end of jPair of the j-th and (j+4)-th layers, j=1,3
569  for (unsigned int jPair = 0; jPair < nPairsInOuter; ++jPair) {
570  if(foundCounterparts[jPair]) continue; // in case of no counterpart of pairs in the outner layers
571  if((hitIdsInTwo[iQuad+2].at(jPair)>>16 & 0xffff) == 0xffff || (hitIdsInTwo[iQuad+2].at(jPair) & 0xffff) == 0xffff) continue;
572 // unsigned long int encodedIds = (0xffffffff << 32) + hitIdsInTwo[iQuad+2].at(jPair);
573  unsigned long int encodedIds = (0xffffffff00000000 ) + hitIdsInTwo[iQuad+2].at(jPair);
574  hitIdsInFour[iQuad].push_back(encodedIds);
575  slopeInFour[iQuad].push_back(slopeInTwo[iQuad+2].at(jPair));
576  interceptInFour[iQuad].push_back(interceptInTwo[iQuad+2].at(jPair));
577  }
578  }// end of quad loop
579 
580  ATH_MSG_DEBUG("@@STGC@@ isStrip= " << isStrip << " Nquads " << hitIdsInFour[0].size() << " " << hitIdsInFour[1].size());
581  for (unsigned int iLayer = 0; iLayer < 2; ++ iLayer) {
582  for (unsigned int iQuad = 0; iQuad < slopeInFour[iLayer].size(); ++iQuad) {
583  ATH_MSG_DEBUG("@@STGC@@ quad fit isStrip= " << isStrip << " slope= " << slopeInFour[iLayer].at(iQuad) << " intercept= " << interceptInFour[iLayer].at(iQuad));
584  }
585  }
586 
587  std::vector< std::array<int, 8> > hitIdsInEight;
588  std::vector<double> mseInEight;
589 
590  unsigned int nQuadInInner = hitIdsInFour[0].size();
591  unsigned int nQuadInOuter = hitIdsInFour[1].size();
592 
593  for(unsigned int iQuad = 0; iQuad < nQuadInInner; ++iQuad){
594  double slope[2];
595  double intercept[2];
596  slope[0] = slopeInFour[0].at(iQuad);
597  intercept[0] = interceptInFour[0].at(iQuad);
598 
599  for(unsigned int jQuad = 0; jQuad < nQuadInOuter; ++jQuad){
600  unsigned long int ihitIds = hitIdsInFour[0].at(iQuad);
601  unsigned long int jhitIds = hitIdsInFour[1].at(jQuad);
602  int nOfLayersWithNoHit = 0;
603  for (unsigned int iLayer = 0; iLayer < 4; ++iLayer) {
604  if ( (ihitIds>>(3-iLayer)*16 & 0xffff) == 0xffff ) {++nOfLayersWithNoHit;}
605  if ( (jhitIds>>(3-iLayer)*16 & 0xffff) == 0xffff ) {++nOfLayersWithNoHit;}
606  }
607  if (nOfLayersWithNoHit > 4) continue; // require at least 4 hits in 8 layers
608 
609  slope[1] = slopeInFour[1].at(jQuad);
610  intercept[1] = interceptInFour[1].at(jQuad);
611 
612  if(isStrip) {
613  double spR0 = slope[0] * NSWCenterZ + intercept[0];
614  double spR1 = slope[1] * NSWCenterZ + intercept[1];
615  if(std::abs(spR1 - spR0) > 10. ||
616  std::abs(intercept[1] + intercept[0]) / 2 > 100.) continue; // OPTIMIZE ME!!!
617  } else {
618  if(std::abs(intercept[0]*std::sin(slope[0]) - intercept[1]*std::sin(slope[1])) > 100. ) continue;
619  }
620 
621  std::array<int,8> setOfHitIds = {-1,-1,-1,-1,-1,-1,-1,-1};
622  std::vector<double> r, z;
623  for(unsigned int i = 0; i < 8; ++i) {
624  unsigned int iHitId, iLayer = 0;
625  if (i <= 3) {
626  iHitId = (unsigned int) ((ihitIds>>(3-i)*16) & 0xffff);
627  iLayer = i+3*(i%2);
628  } else {
629  iHitId = (unsigned int) ((jhitIds>>(3-i%4)*16) & 0xffff);
630  iLayer = (i-4)+3*(i%2)+1;
631  }
632  if ( iHitId != 0xffff ) {
633  if (isStrip) {
634  r.push_back(stgcHits.at(iHitId).r);
635  z.push_back(stgcHits.at(iHitId).z);
636  }
637  else {
638  double localPhiCenter;
639  if (stgcHits.at(iHitId).stationPhi<=5) {
640  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHitId).stationPhi-1.);
641  } else {
642  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHitId).stationPhi-9.);
643  }
644  if (stgcHits.at(iHitId).stationName == 57){
645  localPhiCenter += M_PI/8.; // small sTGC sectors
646  if (stgcHits.at(iHitId).stationPhi == 5) localPhiCenter -= 2 * M_PI;
647  }
648 
649  double phiProj = stgcHits.at(iHitId).phi - localPhiCenter;
650  if (phiProj > M_PI) phiProj -= 2.0*M_PI;
651  if (phiProj < -1.*M_PI) phiProj += 2.0*M_PI;
652  r.push_back(phiProj);
653  z.push_back(stgcHits.at(iHitId).r);
654  }
655  setOfHitIds[iLayer] = iHitId;
656  ATH_MSG_DEBUG("@@STGC@@ strip_pos iHitId " << iLayer << " " << iHitId);
657  }
658  }
659  double slopefit=0., interceptfit=99999., mse =-1.;
660  if (isStrip) {
661  LinearFit(z,r,&slopefit,&interceptfit,&mse);
662  } else {
663  double phiavg = 0.;
664  for (unsigned int iHit = 0; iHit < r.size(); ++iHit){
665  phiavg += r.at(iHit);
666  }
667  phiavg /= r.size();
668  mse = 0.;
669  for (unsigned int iHit = 0; iHit < r.size(); ++iHit){
670  mse += std::pow(r.at(iHit) - phiavg,2);
671  }
672  }
673  hitIdsInEight.push_back(setOfHitIds);
674  mseInEight.push_back(mse);
675  } // end of quad loop in the outer layers
676  } // end of quad loop in the inner layers
677  if(!hitIdsInEight.size()){
678  ATH_MSG_DEBUG("No candidate segment found in STGC");
679  return;
680  }
681 
682  ATH_MSG_DEBUG("@@STGC@@ isStrip= " << isStrip << " Noctets " << hitIdsInEight.size());
683  std::vector<int> nOctetSegments;
684  std::vector<int> patternStationName;
685 
686  for (unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet) {
687 
688  bool isFirstHit = true;
689  int hitStationName = 0;
690 
691  int nOctetSegment = 0;
692  ATH_MSG_DEBUG("@@STGC@@ octet fit isStrip= " << isStrip << " mse " << mseInEight.at(iOctet));
693  std::array<int, 8> tmpOctet = hitIdsInEight.at(iOctet);
694  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
695  if (tmpOctet[iLayer] != -1) {
696 
697  if(isFirstHit){
698  hitStationName = stgcHits.at(tmpOctet[iLayer]).stationName;
699  isFirstHit = false;
700 
701  ATH_MSG_DEBUG("@@STGC@@ octet pos isStrip= " << isStrip << " r= " << stgcHits.at(tmpOctet[iLayer]).r << " phi= " << stgcHits.at(tmpOctet[iLayer]).phi << " z= " << stgcHits.at(tmpOctet[iLayer]).z);
702  nOctetSegment++;
703  }
704  else if(stgcHits.at(tmpOctet[iLayer]).stationName == hitStationName){
705  ATH_MSG_DEBUG("@@STGC@@ octet pos isStrip= " << isStrip << " r= " << stgcHits.at(tmpOctet[iLayer]).r << " phi= " << stgcHits.at(tmpOctet[iLayer]).phi << " z= " << stgcHits.at(tmpOctet[iLayer]).z);
706  nOctetSegment++;
707  }
708  }
709  }
710  nOctetSegments.push_back(nOctetSegment);
711  patternStationName.push_back(hitStationName);
712  }
713  double nOcSegMax = 0;
714  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
715  if(nOctetSegments.at(iOctet) > nOcSegMax){
716  nOcSegMax = nOctetSegments.at(iOctet);
717  }
718  }
719 
720  double msemin = 1000000.;
721  double mseminWireL = 1000000.; // Large sector
722  double mseminWireS = 1000000.; // Small sector
723 
724  std::vector<int> octetIds(2,-1);
725 
726  if(isStrip){
727  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
728  if(nOctetSegments.at(iOctet) != nOcSegMax){
729  continue;
730  }
731  if( mseInEight.at(iOctet) < msemin) {
732  msemin = mseInEight.at(iOctet);
733  }
734  }// end of Octet loop
735  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
736  if(nOctetSegments.at(iOctet) != nOcSegMax) continue;
737  if(mseInEight.at(iOctet) != msemin){
738  continue;
739  }
740  octetIds.push_back(iOctet);
741  }
742  } else {
743  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
744  if(patternStationName.at(iOctet) == 58){
745  if( mseInEight.at(iOctet) < mseminWireL) {
746  mseminWireL = mseInEight.at(iOctet);
747  }
748  }
749  else if(patternStationName.at(iOctet) == 57){
750  if( mseInEight.at(iOctet) < mseminWireS) {
751  mseminWireS = mseInEight.at(iOctet);
752  }
753  }
754  }// end of Octet loop
755  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
756  if(patternStationName.at(iOctet) == 58){
757  if(mseInEight.at(iOctet) != mseminWireL){
758  continue;
759  }
760  }
761  else if(patternStationName.at(iOctet) == 57){
762  if(mseInEight.at(iOctet) != mseminWireS){
763  continue;
764  }
765  }
766  octetIds.push_back(iOctet);
767  }
768  }
769  for(unsigned int ids = 0; ids < octetIds.size(); ids++){
770  if (octetIds.at(ids) != -1) {
771  hitIdsCandidate.push_back(hitIdsInEight.at(octetIds.at(ids)));
772  }
773  }
774 }
775 
776 StatusCode TrigL2MuonSA::NswStationFitter::MakeSegment(TrigL2MuonSA::TrackPattern& trackPattern,
777  TrigL2MuonSA::StgcHits& stgcHits) const
778 {
779  TrigL2MuonSA::StgcHits selectedStgcHits;
780  selectedStgcHits.clear();
781  if(stgcHits.size() == 0) return StatusCode::SUCCESS;
782  for(unsigned int iHit = 0; iHit < stgcHits.size(); iHit++){
783  if(stgcHits.at(iHit).isOutlier != 0) continue;
784  selectedStgcHits.push_back(stgcHits.at(iHit));
785  }
786  trackPattern.stgcSegment.clear();
787  trackPattern.stgcSegment = selectedStgcHits;
788  return StatusCode::SUCCESS;
789 }
790 StatusCode TrigL2MuonSA::NswStationFitter::MakeSegment(TrigL2MuonSA::TrackPattern& trackPattern,
791  TrigL2MuonSA::MmHits& mmHits) const
792 {
793  TrigL2MuonSA::MmHits selectedMmHits;
794  selectedMmHits.clear();
795  if(mmHits.size() == 0) return StatusCode::SUCCESS;
796  for(unsigned int iHit = 0; iHit < mmHits.size(); iHit++){
797  if(mmHits.at(iHit).isOutlier != 0) continue;
798  selectedMmHits.push_back(mmHits.at(iHit));
799  }
800  trackPattern.mmSegment.clear();
801  trackPattern.mmSegment = selectedMmHits;
802  return StatusCode::SUCCESS;
803 }
804 
805 
806 void TrigL2MuonSA::NswStationFitter::LinearFit(std::vector<double>& x,std::vector<double>& y,
807  double* slope, double* intercept, double* mse) const
808 {
809  double sumX=0, sumY=0, sumXY=0, sumX2=0;
810  int nHits = x.size();
811  *mse = 0.;
812  for (unsigned int iHit = 0; iHit < x.size(); ++iHit){
813  sumX += x.at(iHit);
814  sumY += y.at(iHit);
815  sumXY += x.at(iHit)*y.at(iHit);
816  sumX2 += x.at(iHit)*x.at(iHit);
817  }
818 
819  if(nHits > 1) {
820  if(nHits * sumX2 - (sumX * sumX) > ZERO_LIMIT) {
821  *slope = (nHits * sumXY - sumX * sumY) / (nHits * sumX2 - (sumX * sumX));
822  *intercept = (sumX2 * sumY - sumXY * sumX) / (nHits * sumX2 - sumX * sumX);
823  } else {
824  *slope = 0.;
825  *intercept = sumY/nHits;
826  }
827  }
828  else if(nHits == 1) {
829  *slope = sumY/sumX;
830  *intercept = 0.;
831  }
832 
833  if(nHits > 2) {
834  for(unsigned int iHit = 0; iHit< x.size(); ++iHit){
835  *mse += std::pow(y.at(iHit) - (*slope * x.at(iHit) + *intercept), 2.0);
836  }
837  *mse = *mse / (nHits - 2);
838  }
839  else {
840  *mse = 1000.;
841  }
842 }
843 
844 void TrigL2MuonSA::NswStationFitter::LinearFitWeight(std::vector<double>& x,std::vector<double>& y,
845  std::vector<bool>& isStgc, double* slope, double* intercept, double* mse, double eta) const
846 {
847  double RmsDeltarEtaStgc[12] = {};
848  double RmsDeltarEtaMm[12]= {};
849  getNswResolution(RmsDeltarEtaStgc, RmsDeltarEtaMm, 12);
850 
851  double weightStgc[12] = {};
852  double weightMm[12] = {};
853  for(int i_weight=0; i_weight<12; i_weight++){
854  weightStgc[i_weight] = 1/std::pow(RmsDeltarEtaStgc[i_weight],2);
855  weightMm[i_weight] = 1/std::pow(RmsDeltarEtaMm[i_weight],2);
856  }
857  int weightBin = 0;
858  double minEta = 1.3;
859  double maxEta = 1.4;
860  for(int iBin=0; iBin<12; iBin++){
861  if(std::abs(eta) >= minEta && std::abs(eta) < maxEta){
862  weightBin = iBin;
863  break;
864  } else {
865  minEta += 0.1;
866  maxEta += 0.1;
867  }
868  }
869 
870  double sumX=0, sumY=0, sumXY=0, sumX2=0, sumW=0;
871  int nHits = x.size();
872  *mse = 0.;
873  for (unsigned int iHit = 0; iHit < x.size(); ++iHit){
874  if(isStgc.at(iHit)){
875  sumX += weightStgc[weightBin] * x.at(iHit);
876  sumY += weightStgc[weightBin] * y.at(iHit);
877  sumXY += weightStgc[weightBin] * x.at(iHit) * y.at(iHit);
878  sumX2 += weightStgc[weightBin] * x.at(iHit) * x.at(iHit);
879  sumW += weightStgc[weightBin];
880  } else {
881  sumX += weightMm[weightBin] * x.at(iHit);
882  sumY += weightMm[weightBin] * y.at(iHit);
883  sumXY += weightMm[weightBin] * x.at(iHit) * y.at(iHit);
884  sumX2 += weightMm[weightBin] * x.at(iHit) * x.at(iHit);
885  sumW += weightMm[weightBin];
886  }
887  }
888 
889  if(nHits > 1) {
890  if(nHits * sumX2 - (sumX * sumX) > ZERO_LIMIT) {
891  *slope = (sumW * sumXY - sumX * sumY) / (sumW * sumX2 - (sumX * sumX));
892  *intercept = (sumX2 * sumY - sumXY * sumX) / (sumW * sumX2 - sumX * sumX);
893  } else {
894  *slope = 0.;
895  *intercept = sumY/nHits;
896  }
897  }
898  else if(nHits == 1) {
899  *slope = sumY/sumX;
900  *intercept = 0.;
901  }
902 
903  if(nHits > 2) {
904  for(unsigned int iHit = 0; iHit< x.size(); ++iHit){
905  *mse += std::pow((y.at(iHit) - (*slope * x.at(iHit) + *intercept)), 2.0);
906  }
907  *mse = *mse / (nHits - 2);
908  }
909  else {
910  *mse = 1000.;
911  }
912 }
913 
914 void TrigL2MuonSA::NswStationFitter::getNswResolution(double *stgcDeltaR, double *mmDeltaR, unsigned int size) const
915 {
916  double RmsDeltarEtaStgc[12] = {1.43,1.53,1.53,1.56,1.59,1.54,1.70,1.69,1.76,1.81,1.83,1.84};
917  double RmsDeltarEtaMm[12] = {0.49,0.46,0.48,0.40,0.39,0.39,0.38,0.35,0.36,0.33,0.33,0.40};
918  for(unsigned int bin=0; bin < size; bin++){
919  stgcDeltaR[bin] = RmsDeltarEtaStgc[bin];
920  mmDeltaR[bin] = RmsDeltarEtaMm[bin];
921  }
922 }
923 
924 StatusCode TrigL2MuonSA::NswStationFitter::calcMergedHit(TrigL2MuonSA::TrackPattern& trackPattern) const
925 {
926  TrigL2MuonSA::StgcHits stgcHits = trackPattern.stgcSegment;
927  TrigL2MuonSA::MmHits mmHits = trackPattern.mmSegment;
928 
929  double side_mm = 0;
930  std::vector<double> r, z;
931  std::vector<bool> isStgc;
932  for(unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit) {
933  if (stgcHits.at(iHit).channelType == 1) {
934  r.push_back(stgcHits.at(iHit).r);
935  z.push_back(stgcHits.at(iHit).z);
936  isStgc.push_back(true);
937  }
938  }
939  for(unsigned int iHit = 0; iHit < mmHits.size(); ++iHit) {
940  if (mmHits.at(iHit).layerNumber < 2 || mmHits.at(iHit).layerNumber > 5) {
941  r.push_back(mmHits.at(iHit).r);
942  z.push_back(mmHits.at(iHit).z);
943  isStgc.push_back(false);
944  side_mm = std::abs(mmHits.at(iHit).z)/mmHits.at(iHit).z;
945  }
946  }
947  double slopefit=0., interceptfit=99999., mse=-1.;
948  LinearFit(z,r,&slopefit,&interceptfit,&mse);
949  ATH_MSG_DEBUG("@@Merge@@ stgc_mmX_fit slope= " << slopefit);
950  ATH_MSG_DEBUG("@@Merge@@ stgc_mmX_fit intercept= " << interceptfit);
951  ATH_MSG_DEBUG("@@Merge@@ stgc_mmX_fit mse= " << mse);
952 
953  std::vector<double> phiLocal;
954  double localPhiCenter = 3.*M_PI;
955  for (unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit) {
956  if (stgcHits.at(iHit).channelType != 2) {continue;}
957  if (stgcHits.at(iHit).stationPhi<=5) {
958  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHit).stationPhi-1.);
959  } else {
960  localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHit).stationPhi-9.);
961  }
962 
963  if (stgcHits.at(iHit).stationName == 57){
964  localPhiCenter += M_PI/8.; // small sTGC sectors
965  if (stgcHits.at(iHit).stationPhi == 5) localPhiCenter -= 2. * M_PI;
966  }
967 
968  double phiProj = stgcHits.at(iHit).phi - localPhiCenter;
969  if (phiProj > M_PI) phiProj -= 2.0*M_PI;
970  if (phiProj < -1.*M_PI) phiProj += 2.0*M_PI;
971  double rInterpolate = slopefit * stgcHits.at(iHit).z + interceptfit;
972  double rProj = stgcHits.at(iHit).r;
973  phiLocal.push_back( std::atan(rProj/rInterpolate*std::tan(phiProj)) );
974  ATH_MSG_DEBUG("@@Merge@@ philocalwire " << stgcHits.at(iHit).stationPhi << " " << stgcHits.at(iHit).stationName << " "
975  << localPhiCenter << " " << stgcHits.at(iHit).phi << " " << stgcHits.at(iHit).r << " " << stgcHits.at(iHit).z);
976  ATH_MSG_DEBUG("@@Merge@@ philocalwire " << rProj << " " << rInterpolate << " " << std::tan(phiProj) );
977  ATH_MSG_DEBUG("@@Merge@@ philocalwire " << std::atan(rProj/rInterpolate*std::tan(phiProj)) );
978  }
979  double tanTiltAngleU = 0,
980  tanTiltAngleV = 0;
981  double cosTiltAngleU = 0,
982  cosTiltAngleV = 0;
983  double sinTiltAngleU = 0,
984  sinTiltAngleV = 0;
985  if(side_mm > ZERO_LIMIT){
986  tanTiltAngleU = TanM1p5,
987  tanTiltAngleV = TanP1p5;
988  cosTiltAngleU = CosM1p5,
989  cosTiltAngleV = CosP1p5;
990  sinTiltAngleU = SinM1p5,
991  sinTiltAngleV = SinP1p5;
992  } else if(side_mm < -1.*ZERO_LIMIT){
993  tanTiltAngleU = TanP1p5,
994  tanTiltAngleV = TanM1p5;
995  cosTiltAngleU = CosP1p5,
996  cosTiltAngleV = CosM1p5;
997  sinTiltAngleU = SinP1p5,
998  sinTiltAngleV = SinM1p5;
999  } else {
1000  ATH_MSG_DEBUG("@@Merge@@ no U, V layer hits -> not consider tilt of U/V layers");
1001  }
1002  for (unsigned int iHit = 0; iHit < mmHits.size(); ++iHit) {
1003  if (localPhiCenter > 2.*M_PI) {
1004  if (mmHits.at(iHit).stationPhi<=5) {
1005  localPhiCenter = 0.25 * M_PI * ((double)mmHits.at(iHit).stationPhi-1.);
1006  } else {
1007  localPhiCenter = 0.25 * M_PI * ((double)mmHits.at(iHit).stationPhi-9.);
1008  }
1009  if (mmHits.at(iHit).stationName == 55){
1010  localPhiCenter += M_PI/8.; // small MM sectors
1011  if (mmHits.at(iHit).stationPhi == 5) localPhiCenter -= 2 * M_PI;
1012  }
1013  }
1014  if (mmHits.at(iHit).layerNumber >1 && mmHits.at(iHit).layerNumber < 6){
1015  double rInterpolate = slopefit * mmHits.at(iHit).z + interceptfit;
1016  double rProj = mmHits.at(iHit).r;
1017  if(std::abs(side_mm) < ZERO_LIMIT) {
1018  phiLocal.push_back(0);
1019  ATH_MSG_DEBUG("@@Merge@@ philocalmm 0");
1020  }
1021  else if ((mmHits.at(iHit).layerNumber)%2 == 0) { // layer U
1022  phiLocal.push_back( std::atan((rProj-rInterpolate)/tanTiltAngleU/rInterpolate));
1023  ATH_MSG_DEBUG("@@Merge@@ philocalmmU " << std::atan((rProj-rInterpolate)/tanTiltAngleU/rInterpolate));
1024  } else { // layer V
1025  phiLocal.push_back( std::atan((rProj-rInterpolate)/tanTiltAngleV/rInterpolate));
1026  ATH_MSG_DEBUG("@@Merge@@ philocalmmV " << std::atan((rProj-rInterpolate)/tanTiltAngleV/rInterpolate));
1027  }
1028  }
1029  }
1030  double sumPhiLocal = 0;
1031  for (unsigned int iHit = 0; iHit < phiLocal.size(); ++iHit ) {
1032  sumPhiLocal += phiLocal.at(iHit);
1033  }
1034  double phiLocalAvg = sumPhiLocal/phiLocal.size();
1035  ATH_MSG_DEBUG("@@Merge@@ philocalAvg " << phiLocalAvg);
1036 
1037  r.clear();
1038  z.clear();
1039  isStgc.clear();
1040  std::vector<double> r_stgc, z_stgc, r_mm, z_mm;
1041  std::vector<bool> isStgc_stgc, isStgc_mm;
1042  double side_stgc = 0;
1043  for(unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit) {
1044  if (stgcHits.at(iHit).stationPhi<=5) localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHit).stationPhi-1.);
1045  if (stgcHits.at(iHit).stationPhi> 5) localPhiCenter = 0.25 * M_PI * ((double)stgcHits.at(iHit).stationPhi-9.);
1046  if (stgcHits.at(iHit).channelType == 1) {
1047  r_stgc.push_back(stgcHits.at(iHit).r/std::cos(phiLocalAvg));
1048  z_stgc.push_back(stgcHits.at(iHit).z);
1049  isStgc_stgc.push_back(true);
1050  side_stgc = std::abs(stgcHits.at(iHit).z)/stgcHits.at(iHit).z;
1051 
1052  ATH_MSG_DEBUG("@@Merge@@ stgc strip_r " << phiLocalAvg << " " << stgcHits.at(iHit).z << " " << stgcHits.at(iHit).r/std::cos(phiLocalAvg));
1053  }
1054 
1055  }
1056  double slopefit_stgc=0., interceptfit_stgc=99999., mse_stgc=1.e20;
1057  if(r_stgc.size() == 0) {
1058  ATH_MSG_DEBUG("No STGC hit to calculate superoint");
1059  } else {
1060  LinearFit(z_stgc,r_stgc,&slopefit_stgc,&interceptfit_stgc,&mse_stgc);
1061  }
1062 
1063  for(unsigned int iHit = 0; iHit < mmHits.size(); ++iHit) {
1064  if (mmHits.at(iHit).layerNumber < 2 || mmHits.at(iHit).layerNumber > 5) {
1065  r_mm.push_back(mmHits.at(iHit).r/std::cos(phiLocalAvg));
1066  z_mm.push_back(mmHits.at(iHit).z);
1067  isStgc_mm.push_back(false);
1068  } else {
1069  z_mm.push_back(mmHits.at(iHit).z);
1070  isStgc_mm.push_back(false);
1071 
1072  double rProj = mmHits.at(iHit).r;
1073  if(std::abs(side_mm) < ZERO_LIMIT) { // no layer U/V
1074  r_mm.push_back(rProj);
1075  }
1076  else if ((mmHits.at(iHit).layerNumber)%2 == 0) { // layer U
1077  double rPrime = (rProj * cosTiltAngleU)/(cos(phiLocalAvg)*cosTiltAngleU + sin(phiLocalAvg)*sinTiltAngleU);
1078  r_mm.push_back(rPrime);
1079  } else { //layer V
1080  double rPrime = (rProj * cosTiltAngleV)/(cos(phiLocalAvg)*cosTiltAngleV + sin(phiLocalAvg)*sinTiltAngleV);
1081  r_mm.push_back(rPrime);
1082  }
1083  }
1084  }
1085  double slopefit_mm=0., interceptfit_mm=99999., mse_mm=1.e20;
1086  if(r_mm.size() == 0) {
1087  ATH_MSG_WARNING("No MM hit to calculate superoint");
1088  } else {
1089  LinearFit(z_mm,r_mm,&slopefit_mm,&interceptfit_mm,&mse_mm);
1090  }
1091 
1092  unsigned int fmerge = 0;
1093  slopefit=0., interceptfit=99999., mse=1.e20;
1094  double side = 0;
1095  double StgcSegZ = 7526.329;
1096  double StgcSegR = 0;
1097  double MmSegZ = 7526.329;
1098  double MmSegR = 0;
1099  if (mse_stgc < 1.e7 && mse_mm < 1.e7) {
1100  r = r_stgc;
1101  copy(r_mm.begin(), r_mm.end(), back_inserter(r));
1102  z = z_stgc;
1103  copy(z_mm.begin(), z_mm.end(), back_inserter(z));
1104  isStgc = isStgc_stgc;
1105  copy(isStgc_mm.begin(), isStgc_mm.end(), back_inserter(isStgc));
1106 
1107  if(side_stgc < -1.*ZERO_LIMIT){
1108  StgcSegZ = -7526.329;
1109  }
1110  StgcSegR = slopefit_stgc * StgcSegZ + interceptfit_stgc;
1111  double StgcSegOriginTheta = std::atan(StgcSegR / StgcSegZ);
1112  double StgcSegEta = side_stgc * (- std::log(std::abs(std::tan(StgcSegOriginTheta / 2))));
1113  if(side_mm < -1.*ZERO_LIMIT){
1114  MmSegZ = -7526.329;
1115  }
1116  MmSegR = slopefit_mm * MmSegZ + interceptfit_mm;
1117  double MmSegOriginTheta = std::atan(MmSegR / MmSegZ);
1118  double MmSegEta = side_stgc * (- std::log(std::abs(std::tan(MmSegOriginTheta / 2))));
1119 
1120  double SegEtaAve = 0;
1121  if(std::abs(side_stgc) > ZERO_LIMIT || std::abs(side_mm) > ZERO_LIMIT){
1122  if(side_stgc*side_mm > 0){
1123  side = side_stgc;
1124  SegEtaAve = (StgcSegEta + MmSegEta)/2;
1125  } else if(std::abs(side_stgc) < ZERO_LIMIT) {
1126  side = side_mm;
1127  SegEtaAve = MmSegEta;
1128  } else if(std::abs(side_mm) < ZERO_LIMIT) {
1129  side = side_stgc;
1130  SegEtaAve = StgcSegEta;
1131  }
1132  }
1133 
1134  LinearFitWeight(z,r,isStgc,&slopefit,&interceptfit,&mse,SegEtaAve);
1135  fmerge = 1;
1136  }
1137  if (mse > 1.e7) {
1138  if (mse_stgc < mse_mm) {
1139  slopefit = slopefit_stgc;
1140  interceptfit = interceptfit_stgc;
1141  mse = mse_stgc;
1142  z = z_stgc;
1143  side = side_stgc;
1144  fmerge = 2;
1145  } else {
1146  slopefit = slopefit_mm;
1147  interceptfit = interceptfit_mm;
1148  mse = mse_mm;
1149  z = z_mm;
1150  side = side_mm;
1151  fmerge = 3;
1152  }
1153  }
1154  if (mse > 1.e19) {
1155  ATH_MSG_WARNING("No sTGC and MM hit to calculate superoint");
1156  return StatusCode::SUCCESS;
1157  }
1158 
1159  // store superpoint info in TrackData
1160  xAOD::L2MuonParameters::Chamber inner = xAOD::L2MuonParameters::Chamber::EndcapInner;
1161  TrigL2MuonSA::SuperPoint* superPoint = &(trackPattern.superPoints[inner]);
1162  double NSWCenterZ = 7526.329;
1163  if(side < -1.*ZERO_LIMIT){
1164  NSWCenterZ = -7526.329;
1165  }
1166  superPoint->R = slopefit * NSWCenterZ + interceptfit;
1167  superPoint->Phim = phiLocalAvg+localPhiCenter;
1168  superPoint->Z = NSWCenterZ;
1169  superPoint->Npoint = z.size();
1170 
1171  if (NSWCenterZ != 0) superPoint->Alin = slopefit;
1172  superPoint->Blin = interceptfit;
1173 
1174  ATH_MSG_DEBUG("Nsw Super Point r/phi/z/slope = "<<superPoint->R<<"/"<<superPoint->Phim<<"/"<<superPoint->Z<<"/"<<superPoint->Alin);
1175 
1176  ATH_MSG_DEBUG("@@Merge@@ Nsw Super Point r/phi/z/slope = "<<superPoint->R<<"/"<<superPoint->Phim<<"/"<<superPoint->Z<<"/"<<superPoint->Alin);
1177  ATH_MSG_DEBUG("@@Merge@@ fit slope= " << slopefit << " " << slopefit_stgc << " " << slopefit_mm);
1178  ATH_MSG_DEBUG("@@Merge@@ fit intercept= " << interceptfit << " " << interceptfit_stgc << " " << interceptfit_mm);
1179  ATH_MSG_DEBUG("@@Merge@@ fit mse= " << mse << " " << mse_stgc << " " << mse_mm);
1180  ATH_MSG_DEBUG("@@Merge@@ fit tech= " << fmerge);
1181 
1182 
1183  return StatusCode::SUCCESS;
1184 
1185 }
1186 
1187 StatusCode TrigL2MuonSA::NswStationFitter::findMmHitsInSegment(TrigL2MuonSA::MmHits& mmHits) const
1188 {
1189  if(mmHits.size() == 0) return StatusCode::SUCCESS;
1190  int hitsInRoad = 0;
1191  for(unsigned int iHit = 0; iHit < mmHits.size(); iHit++){
1192  if(mmHits.at(iHit).isOutlier == 0){
1193  hitsInRoad++;
1194  mmHits.at(iHit).isOutlier = 1;
1195  }
1196  }
1197  if(hitsInRoad == 0) return StatusCode::SUCCESS;
1198 
1199  if(hitsInRoad < 6) {
1200  ATH_MSG_DEBUG("Number of MM hits is too small, at least 6 hits required : "<<hitsInRoad<<" hits");
1201  return StatusCode::SUCCESS;
1202  } else if(hitsInRoad > 100) {
1203  ATH_MSG_WARNING("Number of MM hits is too large, at most (2^16 - 1) hits allowed : "<<hitsInRoad<<" hits");
1204  return StatusCode::SUCCESS;
1205  }
1206 
1207  std::array< std::vector<int>, 8 > hitIdByLayer;
1208  for(unsigned int iHit = 0; iHit < mmHits.size(); ++iHit){
1209  if(mmHits.at(iHit).isOutlier != 1) continue;
1210  int layerNumber = mmHits.at(iHit).layerNumber;
1211  if (layerNumber > 7) {
1212  ATH_MSG_WARNING("MM hit layer number > 7");
1213  continue;
1214  }
1215  hitIdByLayer[layerNumber].push_back(iHit);
1216  }
1217  ATH_MSG_DEBUG("@@MM@@ Nhits " << hitIdByLayer[0].size()
1218  << " " << hitIdByLayer[1].size()
1219  << " " << hitIdByLayer[2].size()
1220  << " " << hitIdByLayer[3].size()
1221  << " " << hitIdByLayer[4].size()
1222  << " " << hitIdByLayer[5].size()
1223  << " " << hitIdByLayer[6].size()
1224  << " " << hitIdByLayer[7].size());
1225 
1226  std::vector< std::array<int, 8> > mmHitIds;
1227  findSetOfMmHitIds(mmHits, hitIdByLayer, mmHitIds);
1228  for (unsigned int iHit = 0; iHit < mmHitIds.size(); ++iHit) {
1229  std::array<int, 8> hitIds = mmHitIds.at(iHit);
1230  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
1231  if (hitIds[iLayer] != -1) {
1232  mmHits.at(hitIds[iLayer]).isOutlier = 0;
1233  }
1234  }
1235  }
1236  return StatusCode::SUCCESS;
1237 }
1238 
1239 void TrigL2MuonSA::NswStationFitter::findSetOfMmHitIds(TrigL2MuonSA::MmHits& mmHits,
1240  std::array<std::vector<int>,8> hitIdByLayer,
1241  std::vector<std::array<int, 8>>& hitIdsCandidate) const
1242 {
1243 
1244  double NSWCenterZ = 7526.329;
1245  int side = 0;
1246  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
1247  if ( hitIdByLayer[iLayer].size() > 0) {
1248  side = std::abs(mmHits.at(hitIdByLayer[iLayer].at(0)).z)/mmHits.at(hitIdByLayer[iLayer].at(0)).z;
1249  break;
1250  }
1251  }
1252  NSWCenterZ = NSWCenterZ * side;
1253 
1254  std::array<std::vector<unsigned long int>,4> hitIdsInTwo;
1255  std::array<std::vector<double>,4> slopeInTwo;
1256  std::array<std::vector<double>,4> interceptInTwo;
1257  // Loop over pairs of the i-th and the (i+6)-th layers
1258  for(unsigned int iPair = 0; iPair < 4; ++iPair){
1259 
1260  unsigned int nHitsInInner = hitIdByLayer[iPair].size();
1261  unsigned int nHitsInOuter;
1262  if (iPair < 2) { // paris of X layers
1263  nHitsInOuter = hitIdByLayer[iPair+6].size();
1264  } else { // pairs of U or V layers
1265  nHitsInOuter = hitIdByLayer[iPair+2].size();
1266  }
1267 
1268  if ( nHitsInInner > 0xffff-1 || nHitsInOuter > 0xffff-1) {
1269  ATH_MSG_WARNING("Number of Mm hits in layers exceeds the limit of (2^16 - 1) : Number of Mm hits in "<<iPair<<"th layer = "<< nHitsInInner
1270  <<", Number of Mm hits in "<<iPair+4<<"th layer = "<<nHitsInOuter);
1271  ATH_MSG_WARNING("Number of Mm hits is limitted to (2^16 - 1) and hits with id more than (2^16 -1) will be trancated.");
1272  if (nHitsInInner > 0xffff-1) {nHitsInInner = 0xffff-1;}
1273  if (nHitsInOuter > 0xffff-1) {nHitsInOuter = 0xffff-1;}
1274  }
1275 
1276  bool foundCounterparts[0xffff] = {};
1277  // Loop over hits in the i-th layer
1278  for(unsigned int iHit = 0; iHit < nHitsInInner; ++iHit){
1279 
1280  bool foundCounterpart = 0;
1281 
1282  double z[2] = {};
1283  double r[2] = {};
1284 
1285  int iHitId = hitIdByLayer[iPair].at(iHit);
1286  r[0] = mmHits.at(iHitId).r;
1287  z[0] = mmHits.at(iHitId).z;
1288 
1289  // Loop over hits in the (i+6)-th layer
1290  for(unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
1291 
1292  int jHitId;
1293  if (iPair < 2) {
1294  jHitId = hitIdByLayer[iPair+6].at(jHit);
1295  } else {
1296  jHitId = hitIdByLayer[iPair+2].at(jHit);
1297  }
1298  r[1] = mmHits.at(jHitId).r;
1299  z[1] = mmHits.at(jHitId).z;
1300 
1301  double slope = (r[1] - r[0]) / (z[1] - z[0]);
1302  double intercept = slope*(0. - z[0]) + r[0];
1303  // select pairs whose slops in limited regions
1304  if(std::abs(slope) < 0.1 || std::abs(slope) > 0.7 || std::abs(intercept) > 500.) continue;
1305 
1306  int encodedIds = (iHitId<<16) + jHitId;
1307  hitIdsInTwo[iPair].push_back(encodedIds);
1308  slopeInTwo[iPair].push_back(slope);
1309  interceptInTwo[iPair].push_back(intercept);
1310 
1311  foundCounterpart = 1;
1312  foundCounterparts[jHit] = 1;
1313  }//end of jHit in the (i+6)-th layer
1314  if(!foundCounterpart){ // in case of no counterpart in the (i+4)-th layer
1315  int encodedIds = (iHitId<<16) + 0xffff; // fill all bits with 1 for hit id for the layer with no hit
1316  hitIdsInTwo[iPair].push_back(encodedIds);
1317  slopeInTwo[iPair].push_back(r[0]/z[0]);
1318  interceptInTwo[iPair].push_back(0.);
1319  }
1320  }//end of iHit in the i-th layer
1321  // Loop over hits in the (i+4)-th layer
1322  for(unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
1323  if (!foundCounterparts[jHit]) {
1324  int jHitId;
1325  if (iPair < 2) {
1326  jHitId = hitIdByLayer[iPair+6].at(jHit);
1327  } else {
1328  jHitId = hitIdByLayer[iPair+2].at(jHit);
1329  }
1330  int encodedIds = (0xFFFFu<<16) + jHitId; // fill all bits with 1 for hit id for the layer with no hit
1331  hitIdsInTwo[iPair].push_back(encodedIds);
1332  slopeInTwo[iPair].push_back(mmHits.at(jHitId).r/mmHits.at(jHitId).z);
1333  interceptInTwo[iPair].push_back(0.);
1334  }
1335  }
1336  }//end of pair loop
1337  ATH_MSG_DEBUG("@@MM@@ Npairs " << hitIdsInTwo[0].size() << " " << hitIdsInTwo[1].size() << " " << hitIdsInTwo[2].size() << " " << hitIdsInTwo[3].size());
1338  for (unsigned int iLayer = 0; iLayer < 4; ++ iLayer) {
1339  for (unsigned int iPair = 0; iPair < slopeInTwo[iLayer].size(); ++iPair) {
1340  ATH_MSG_DEBUG("@@MM@@ pair fit slope= " << slopeInTwo[iLayer].at(iPair) << " intercept= " << interceptInTwo[iLayer].at(iPair));
1341  }
1342  }
1343 
1344  std::vector<std::array<int, 4>> hitIdsInFourX;
1345  std::vector<double> slopeInFourX;
1346  std::vector<double> interceptInFourX;
1347  std::vector<double> mseInFourX;
1348 
1349  unsigned int nPairsInInnerX = hitIdsInTwo[0].size();
1350  unsigned int nPairsInOuterX = hitIdsInTwo[1].size();
1351 
1352  for(unsigned int iPairX = 0; iPairX < nPairsInInnerX; ++iPairX){
1353 
1354  double slope[2];
1355  double intercept[2];
1356  double spR[2];
1357 
1358  slope[0] = slopeInTwo[0].at(iPairX);
1359  intercept[0] = interceptInTwo[0].at(iPairX);
1360  spR[0] = slope[0] * NSWCenterZ + intercept[0];
1361  for(unsigned int jPairX = 0; jPairX < nPairsInOuterX; ++jPairX){
1362  int ihitIds = hitIdsInTwo[0].at(iPairX);
1363  int jhitIds = hitIdsInTwo[1].at(jPairX);
1364  if ( ((ihitIds>>16 & 0xffff) == 0xffff || (ihitIds & 0xffff) == 0xffff) &&
1365  ((jhitIds>>16 & 0xffff) == 0xffff || (jhitIds & 0xffff) == 0xffff )) continue; // require at least 3 hits in 4 layers
1366 
1367  slope[1] = slopeInTwo[1].at(jPairX);
1368  intercept[1] = interceptInTwo[1].at(jPairX);
1369  spR[1] = slope[1] * NSWCenterZ + intercept[1];
1370 
1371  if(std::abs(spR[1] - spR[0]) > 50. ||
1372  std::abs((intercept[1] + intercept[0]) / 2) > 200.) continue;
1373 
1374  std::array<int, 4> setOfHitIds{};
1375  setOfHitIds[0] = (ihitIds>>16 & 0xffff);
1376  setOfHitIds[1] = (ihitIds & 0xffff);
1377  setOfHitIds[2] = (jhitIds>>16 & 0xffff);
1378  setOfHitIds[3] = (jhitIds & 0xffff);
1379  std::vector<double> r;
1380  std::vector<double> z;
1381  for(unsigned int iLayer = 0; iLayer < 4; ++iLayer){
1382  if(setOfHitIds[iLayer] == 0xffff) {
1383  continue;
1384  }
1385  double rhit = mmHits.at(setOfHitIds[iLayer]).r;
1386  double zhit = mmHits.at(setOfHitIds[iLayer]).z;
1387  r.push_back(rhit);
1388  z.push_back(zhit);
1389  }
1390  double slopefit=0., interceptfit=99999., mse=-1.;
1391  LinearFit(z,r,&slopefit, &interceptfit, &mse);
1392 
1393  hitIdsInFourX.push_back(setOfHitIds);
1394  slopeInFourX.push_back(slopefit);
1395  interceptInFourX.push_back(interceptfit);
1396  mseInFourX.push_back(mse);
1397 
1398  }// end of iPair of the i-th and (i+4)-th layers, i=0,2
1399  }// end of jPair of the j-th and (j+4)-th layers, j=1,3
1400  ATH_MSG_DEBUG("@@MM@@ X Nquads " << hitIdsInFourX.size());
1401  for (unsigned int iQuad = 0; iQuad < slopeInFourX.size(); ++iQuad) {
1402  ATH_MSG_DEBUG("@@MM@@ X quad fit slope= " << slopeInFourX.at(iQuad) << " intercept= " << interceptInFourX.at(iQuad) << " mse= " << mseInFourX.at(iQuad));
1403  }
1404 
1405  if(!hitIdsInFourX.size()){
1406  ATH_MSG_WARNING("No candidate segment found in MM X layers");
1407  return;
1408  }
1409 
1410  double tanTiltAngleU = 0,
1411  tanTiltAngleV = 0;
1412  if(side > ZERO_LIMIT){
1413  tanTiltAngleU = tan( 1.5/360.*2.*M_PI),
1414  tanTiltAngleV = tan(-1.5/360.*2.*M_PI);
1415  } else if(side < -1.*ZERO_LIMIT){
1416  tanTiltAngleU = tan(-1.5/360.*2.*M_PI),
1417  tanTiltAngleV = tan(1.5/360.*2.*M_PI);
1418  }
1419 
1420  std::vector< std::array<int, 8> > hitIdsInEight;
1421  std::vector<double> mseInEight;
1422 
1423  for(unsigned int iQuadX = 0; iQuadX < hitIdsInFourX.size(); ++iQuadX){
1424  if(mseInFourX.at(iQuadX) > 10) continue;
1425 
1426  double slopeX = slopeInFourX.at(iQuadX);
1427  double interceptX = interceptInFourX.at(iQuadX);
1428  std::array<int,4> hitIdsX{};
1429  hitIdsX = hitIdsInFourX.at(iQuadX);
1430 
1431  for (unsigned int iPairU = 0; iPairU < hitIdsInTwo[2].size(); ++iPairU) {
1432 
1433  int hitIdsU[2];
1434  hitIdsU[0] = hitIdsInTwo[2].at(iPairU)>>16 & 0xffff;
1435  hitIdsU[1] = hitIdsInTwo[2].at(iPairU) & 0xffff;
1436  double phiLocalU[2]={-99999,-99999};
1437  for(unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1438  if (hitIdsU[iLayer] == 0xffff) continue;
1439  if (hitIdsU[iLayer] < 0) {
1440  ATH_MSG_DEBUG("@@MM@@ hitIdsU[iLayer] iLayer= " << iLayer << " hitIdsU[iLayer]= " << hitIdsU[iLayer]);
1441  }
1442  double rInterpolate = slopeX * mmHits.at(hitIdsU[iLayer]).z + interceptX;
1443  double rProj = mmHits.at(hitIdsU[iLayer]).r;
1444  if(std::abs(tanTiltAngleU) < ZERO_LIMIT)
1445  phiLocalU[iLayer] = 0;
1446  else
1447  phiLocalU[iLayer] = std::atan((rProj-rInterpolate)/tanTiltAngleU/rInterpolate);
1448  }
1449 
1450  for(unsigned int iPairV = 0; iPairV < hitIdsInTwo[3].size(); ++iPairV) {
1451 
1452  int hitIdsV[2];
1453  hitIdsV[0] = hitIdsInTwo[3].at(iPairV)>>16 & 0xffff;
1454  hitIdsV[1] = hitIdsInTwo[3].at(iPairV) & 0xffff;
1455 
1456  if( (hitIdsU[0] == 0xffff || hitIdsU[1] == 0xffff) &&
1457  (hitIdsV[0] == 0xffff || hitIdsV[1] == 0xffff) ) continue; // require at least 3 UV layers having hits
1458 
1459  double phiLocalV[2]={-99999,-99999};
1460  for(unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1461  if (hitIdsV[iLayer] == 0xffff) continue;
1462  if (hitIdsV[iLayer] < 0) {
1463  ATH_MSG_DEBUG("@@MM@@ hitIdsV[iLayer] iLayer= " << iLayer << " hitIdsV[iLayer]= " << hitIdsV[iLayer]);
1464  }
1465  double rInterpolate = slopeX * mmHits.at(hitIdsV[iLayer]).z + interceptX;
1466  double rProj = mmHits.at(hitIdsV[iLayer]).r;
1467  if(std::abs(tanTiltAngleV) < ZERO_LIMIT)
1468  phiLocalV[iLayer] = 0;
1469  else
1470  phiLocalV[iLayer] = std::atan((rProj-rInterpolate)/tanTiltAngleV/rInterpolate);
1471  }
1472 
1473  if ( std::abs(phiLocalU[0]-phiLocalV[0]) > 0.05 &&
1474  std::abs(phiLocalU[1]-phiLocalV[1]) > 0.05) continue;
1475 
1476  // average of phis in 4 UV layers
1477  double phiLocalUV = 0;
1478  int nPhi = 0;
1479  for(unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1480  if(phiLocalU[iLayer] > -99999.) {
1481  phiLocalUV += phiLocalU[iLayer];
1482  ++nPhi;
1483  }
1484  if(phiLocalV[iLayer] > -99999.) {
1485  phiLocalUV += phiLocalV[iLayer];
1486  ++nPhi;
1487  }
1488  }
1489  phiLocalUV /= nPhi;
1490 
1491  std::array<int, 8> setOfHitIds = {-1,-1,-1,-1,-1,-1,-1,-1};
1492  std::vector<double> r, z;
1493  for (unsigned int iLayer = 0; iLayer < 4; ++iLayer) {
1494  if ( hitIdsX[iLayer] != 0xffff) {
1495  if (hitIdsX[iLayer] < 0) {
1496  ATH_MSG_DEBUG("@@MM@@ hitIdsX[iLayer] iLayer= " << iLayer << " hitIdsX[iLayer]= " << hitIdsX[iLayer]);
1497  }
1498  z.push_back(mmHits.at(hitIdsX[iLayer]).z);
1499  r.push_back(mmHits.at(hitIdsX[iLayer]).r / std::cos(phiLocalUV));
1500  setOfHitIds[iLayer] = hitIdsX[iLayer];
1501  }
1502  }
1503  for (unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1504  if ( hitIdsU[iLayer] != 0xffff) {
1505  if (hitIdsU[iLayer] < 0) {
1506  ATH_MSG_DEBUG("@@MM@@ 2 hitIdsU[iLayer] iLayer= " << iLayer << " hitIdsU[iLayer]= " << hitIdsU[iLayer]);
1507  }
1508  z.push_back(mmHits.at(hitIdsU[iLayer]).z);
1509  r.push_back(mmHits.at(hitIdsU[iLayer]).r*(std::cos(phiLocalUV) + 1/std::cos(phiLocalUV))/2.);
1510  setOfHitIds[iLayer+4] = hitIdsU[iLayer];
1511  }
1512  if ( hitIdsV[iLayer] != 0xffff) {
1513  if (hitIdsV[iLayer] < 0) {
1514  ATH_MSG_DEBUG("@@MM@@ 2 hitIdsV[iLayer] iLayer= " << iLayer << " hitIdsV[iLayer]= " << hitIdsV[iLayer]);
1515  }
1516  z.push_back(mmHits.at(hitIdsV[iLayer]).z);
1517  r.push_back(mmHits.at(hitIdsV[iLayer]).r*(std::cos(phiLocalUV) + 1/std::cos(phiLocalUV))/2.);
1518  setOfHitIds[iLayer+6] = hitIdsV[iLayer];
1519  }
1520  }
1521  double slopefit=0., interceptfit=99999., mse=-1.;
1522  LinearFit(z,r,&slopefit,&interceptfit,&mse);
1523 
1524  hitIdsInEight.push_back(setOfHitIds);
1525  mseInEight.push_back(mse);
1526  } // end of Pair loop of V layers
1527  } // end of Pair loop of U Layers
1528  }// end of Quad loop of X layers
1529  ATH_MSG_DEBUG("@@MM@@ Noctets " << hitIdsInEight.size());
1530 
1531  std::vector<int> nOctetSegments;
1532  std::vector<int> patternStationName;
1533  for (unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet) {
1534  bool isFirstHit = true;
1535  int hitStationName = 0;
1536  int nOctetSegment = 0;
1537  ATH_MSG_DEBUG("@@MM@@ octet fit mse " << mseInEight.at(iOctet));
1538  std::array<int, 8> tmpOctet = hitIdsInEight.at(iOctet);
1539  for (unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
1540  if (tmpOctet[iLayer] != -1) {
1541 
1542  if(isFirstHit){
1543  hitStationName = mmHits.at(tmpOctet[iLayer]).stationName;
1544  isFirstHit = false;
1545 
1546  ATH_MSG_DEBUG("@@MM@@ octet pos r= " << mmHits.at(tmpOctet[iLayer]).r << " phi= " << mmHits.at(tmpOctet[iLayer]).phi << " z= " << mmHits.at(tmpOctet[iLayer]).z);
1547  nOctetSegment++;
1548  }
1549  else if(mmHits.at(tmpOctet[iLayer]).stationName == hitStationName){
1550  ATH_MSG_DEBUG("@@MM@@ octet pos r= " << mmHits.at(tmpOctet[iLayer]).r << " phi= " << mmHits.at(tmpOctet[iLayer]).phi << " z= " << mmHits.at(tmpOctet[iLayer]).z);
1551  nOctetSegment++;
1552  }
1553 
1554  }
1555  }
1556  nOctetSegments.push_back(nOctetSegment);
1557  patternStationName.push_back(hitStationName);
1558  }
1559 
1560  double mseminL = 100000.;
1561  double mseminS = 100000.;
1562  std::vector<int> octetIds(2,-1);
1563  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
1564  if(patternStationName.at(iOctet) == 56){
1565  if( mseInEight.at(iOctet) < mseminL) {
1566  mseminL = mseInEight.at(iOctet);
1567  }
1568  }
1569  else if(patternStationName.at(iOctet) == 55){
1570  if( mseInEight.at(iOctet) < mseminS) {
1571  mseminS = mseInEight.at(iOctet);
1572  }
1573  }
1574  }// end of Octet loop
1575 
1576  for(unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
1577  if(patternStationName.at(iOctet) == 56){
1578  if( mseInEight.at(iOctet) != mseminL) {
1579  continue;
1580  }
1581  }
1582  else if(patternStationName.at(iOctet) == 55){
1583  if( mseInEight.at(iOctet) != mseminS) {
1584  continue;
1585  }
1586  }
1587  octetIds.push_back(iOctet);
1588  }
1589 
1590  for(unsigned int ids = 0; ids < octetIds.size(); ids++){
1591  if (octetIds.at(ids) != -1) {
1592  hitIdsCandidate.push_back(hitIdsInEight.at(octetIds.at(ids)));
1593  }
1594  }
1595 }
beamspotman.r
def r
Definition: beamspotman.py:676
TrigL2MuonSA::TrackPattern::stgcSegment
TrigL2MuonSA::StgcHits stgcSegment
Definition: TrackData.h:59
TrigL2MuonSA::TrackPattern::superPoints
TrigL2MuonSA::SuperPoint superPoints[s_NCHAMBER]
Definition: TrackData.h:60
TrigL2MuonSA::NswStationFitter::superPointFitter
StatusCode superPointFitter(const TrigRoiDescriptor *p_roids, TrigL2MuonSA::TrackPattern &trackPattern, TrigL2MuonSA::StgcHits &stgcHits, TrigL2MuonSA::MmHits &mmHits) const
Definition: NswStationFitter.cxx:33
xAOD::L2MuonParameters::Chamber
Chamber
Define chamber types and locations.
Definition: TrigMuonDefs.h:15
AthMsgStreamMacros.h
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TrigL2MuonSA::NswStationFitter::NswStationFitter
NswStationFitter(const std::string &type, const std::string &name, const IInterface *parent)
Definition: NswStationFitter.cxx:26
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
TrigL2MuonSA::StgcHitData::stationEta
int stationEta
Definition: StgcData.h:36
TrigL2MuonSA::TrackPattern::mmSegment
TrigL2MuonSA::MmHits mmSegment
Definition: TrackData.h:58
TrigL2MuonSA::MmHitData::eta
double eta
Definition: MmData.h:28
TrigL2MuonSA::NswStationFitter::MakeSegment
StatusCode MakeSegment(TrigL2MuonSA::TrackPattern &trackPattern, TrigL2MuonSA::StgcHits &stgcHits) const
Definition: NswStationFitter.cxx:776
TrigL2MuonSA::MmHitData::phi
double phi
Definition: MmData.h:32
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
TrigL2MuonSA::MmHitData::stationName
int stationName
Definition: MmData.h:37
M_PI
#define M_PI
Definition: ActiveFraction.h:11
bin
Definition: BinsDiffFromStripMedian.h:43
xAOD::etaMax
etaMax
Definition: HIEventShape_v2.cxx:46
TrigL2MuonSA::NswStationFitter::findStgcHitsInSegment
StatusCode findStgcHitsInSegment(TrigL2MuonSA::StgcHits &stgcHits) const
Definition: NswStationFitter.cxx:269
TrigRoiDescriptor
nope - should be used for standalone also, perhaps need to protect the class def bits #ifndef XAOD_AN...
Definition: TrigRoiDescriptor.h:56
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
TrigL2MuonSA::StgcHitData::phi
double phi
Definition: StgcData.h:33
x
#define x
TrigL2MuonSA::NswStationFitter::LinearFitWeight
void LinearFitWeight(std::vector< double > &x, std::vector< double > &y, std::vector< bool > &isStgc, double *slope, double *intercept, double *mse, double eta) const
Definition: NswStationFitter.cxx:844
TrigL2MuonSA::SuperPoint::Blin
float Blin
Definition: SuperPointData.h:106
TrigVSI::AlgConsts::nPhi
constexpr int nPhi
Default bin number of phi for vertex map.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:27
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:83
TrigL2MuonSA::SuperPoint::R
float R
Definition: SuperPointData.h:102
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
NswStationFitter.h
PUfitVar::maxEta
constexpr float maxEta
Definition: GepMETPufitAlg.cxx:13
TrigL2MuonSA::TrackPattern
Definition: TrackData.h:16
TRT::Hit::side
@ side
Definition: HitInfo.h:83
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
RecMuonRoIUtils.h
TrigL2MuonSA::SuperPoint::Phim
float Phim
Definition: SuperPointData.h:104
TrigL2MuonSA::MmHitData::r
double r
Definition: MmData.h:33
TrigL2MuonSA::MmHitData::stationPhi
int stationPhi
Definition: MmData.h:36
lumiFormat.i
int i
Definition: lumiFormat.py:92
TrigL2MuonSA::StgcHitData::r
double r
Definition: StgcData.h:34
z
#define z
TrigL2MuonSA::StgcHitData
Definition: StgcData.h:14
TrigL2MuonSA::NswStationFitter::findSetOfStgcHitIds
void findSetOfStgcHitIds(TrigL2MuonSA::StgcHits &stgcHits, std::array< std::vector< int >, 8 > hitIdByLayer, std::vector< std::array< int, 8 >> &hitIdsCandidate) const
Definition: NswStationFitter.cxx:361
TrigL2MuonSA::NswStationFitter::findSetOfMmHitIds
void findSetOfMmHitIds(TrigL2MuonSA::MmHits &mmHits, std::array< std::vector< int >, 8 > hitIdByLayer, std::vector< std::array< int, 8 >> &hitIdsCandidate) const
Definition: NswStationFitter.cxx:1239
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
test_pyathena.parent
parent
Definition: test_pyathena.py:15
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TrigL2MuonSA::NswStationFitter::selectMmHits
StatusCode selectMmHits(const TrigRoiDescriptor *p_roids, TrigL2MuonSA::MmHits &mmHits) const
Definition: NswStationFitter.cxx:134
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
TrigL2MuonSA::StgcHitData::channelType
int channelType
Definition: StgcData.h:43
ZERO_LIMIT
const float ZERO_LIMIT
Definition: VP1TriggerHandleL2.cxx:37
lumiFormat.array
array
Definition: lumiFormat.py:98
TrigL2MuonSA::NswStationFitter::calcWeightedSumHit
StatusCode calcWeightedSumHit(TrigL2MuonSA::TrackPattern &trackPattern) const
Definition: NswStationFitter.cxx:184
dumpTgcDigiThreshold.isStrip
list isStrip
Definition: dumpTgcDigiThreshold.py:33
TrigL2MuonSA::StgcHitData::stationPhi
int stationPhi
Definition: StgcData.h:37
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
python.subdetectors.mmg.ids
ids
Definition: mmg.py:8
TrigL2MuonSA::MmHitData
Definition: MmData.h:14
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
TrigL2MuonSA::StgcHitData::stationName
int stationName
Definition: StgcData.h:38
TrigL2MuonSA::NswStationFitter::getNswResolution
void getNswResolution(double *stgcDeltaR, double *mmDeltaR, unsigned int size) const
Definition: NswStationFitter.cxx:914
TrigMuonDefs.h
TrigL2MuonSA::StgcHitData::layerNumber
unsigned int layerNumber
Definition: StgcData.h:42
TrigL2MuonSA::NswStationFitter::calcMergedHit
StatusCode calcMergedHit(TrigL2MuonSA::TrackPattern &trackPattern) const
Definition: NswStationFitter.cxx:924
TrigL2MuonSA::MmHits
std::vector< MmHitData > MmHits
Definition: MmData.h:47
LArCellBinning.etaMin
etaMin
Definition: LArCellBinning.py:84
y
#define y
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
RoiDescriptor::phi
virtual double phi() const override final
Methods to retrieve data members.
Definition: RoiDescriptor.h:100
xAOD::L2MuonParameters::EndcapInner
@ EndcapInner
Inner station in the endcap spectrometer.
Definition: TrigMuonDefs.h:19
TrigL2MuonSA::SuperPoint::Npoint
int Npoint
Definition: SuperPointData.h:97
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
TrigL2MuonSA::StgcHits
std::vector< StgcHitData > StgcHits
Definition: StgcData.h:49
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
RoiDescriptor::eta
virtual double eta() const override final
Definition: RoiDescriptor.h:101
TrigL2MuonSA::NswStationFitter::findMmHitsInSegment
StatusCode findMmHitsInSegment(TrigL2MuonSA::MmHits &mmHits) const
Definition: NswStationFitter.cxx:1187
TrigL2MuonSA::SuperPoint::Z
float Z
Definition: SuperPointData.h:103
TrigL2MuonSA::StgcHitData::z
double z
Definition: StgcData.h:35
TrigL2MuonSA::SuperPoint
Definition: SuperPointData.h:74
calibdata.copy
bool copy
Definition: calibdata.py:27
TrigL2MuonSA::StgcHitData::eta
double eta
Definition: StgcData.h:29
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
TrigL2MuonSA::MmHitData::z
double z
Definition: MmData.h:34
TrigL2MuonSA::NswStationFitter::selectStgcHits
StatusCode selectStgcHits(const TrigRoiDescriptor *p_roids, TrigL2MuonSA::StgcHits &stgcHits) const
Definition: NswStationFitter.cxx:86
TrigL2MuonSA::SuperPoint::Alin
float Alin
Definition: SuperPointData.h:105
readCCLHist.float
float
Definition: readCCLHist.py:83
TrigL2MuonSA::NswStationFitter::LinearFit
void LinearFit(std::vector< double > &x, std::vector< double > &y, double *slope, double *intercept, double *mse) const
Definition: NswStationFitter.cxx:806
TrigL2MuonSA::MmHitData::stationEta
int stationEta
Definition: MmData.h:35
Generated on Sun Jun 30 2024 21:23:41 for ATLAS Offline Software by   doxygen 1.8.18