ATLAS Offline Software
Loading...
Searching...
No Matches
MuFastStationFitter.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
7#include <math.h>
8
9#include "TMath.h"
10
12
14
15namespace{
16 constexpr double ZERO_LIMIT = 1.e-6;
17}
18
19
20// --------------------------------------------------------------------------------
21// --------------------------------------------------------------------------------
22
24 const std::string& name,
25 const IInterface* parent):
26 AthAlgTool(type,name,parent)
27{
28}
29
30// --------------------------------------------------------------------------------
31// --------------------------------------------------------------------------------
32
34{
35 // BackExtrapolator services
36 ATH_CHECK( m_backExtrapolator.retrieve() );
37
38 ATH_CHECK( m_alphaBetaEstimate.retrieve() );
39 ATH_MSG_DEBUG("Retrieved service " << m_alphaBetaEstimate);
40
41 ATH_CHECK( m_ptFromAlphaBeta.retrieve() );
42 ATH_MSG_DEBUG("Retrieved service " << m_ptFromAlphaBeta);
43
44 ATH_CHECK( m_nswStationFitter.retrieve(DisableTool{m_nswStationFitter.empty()}) );
45
46 return StatusCode::SUCCESS;
47}
48
49// --------------------------------------------------------------------------------
50// --------------------------------------------------------------------------------
52{
53 m_use_mcLUT = use_mcLUT;
54
55 if (m_use_mcLUT) {
56 const ServiceHandle<TrigL2MuonSA::PtEndcapLUTSvc> ptEndcapLUTSvc("PtEndcapLUTSvc_MC", name());
57 if ( ptEndcapLUTSvc.retrieve().isFailure() ) {
58 ATH_MSG_ERROR("Could not find PtEndcaplLUTSvc");
59 return StatusCode::FAILURE;
60 }
61 // Calculation of alpha and beta
62 m_alphaBetaEstimate->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
63 // conversion: alpha, beta -> pT
64 m_ptFromAlphaBeta->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
65 } else {
66 const ServiceHandle<TrigL2MuonSA::PtEndcapLUTSvc> ptEndcapLUTSvc("PtEndcapLUTSvc", name());
67 if ( ptEndcapLUTSvc.retrieve().isFailure() ) {
68 ATH_MSG_ERROR("Could not find PtEndcaplLUTSvc");
69 return StatusCode::FAILURE;
70 }
71 // Calculation of alpha and beta
72 m_alphaBetaEstimate->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
73 // conversion: alpha, beta -> pT
74 m_ptFromAlphaBeta->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
75 }
76
77 ATH_MSG_DEBUG( "Completed tp set " << (m_use_mcLUT?"MC":"not MC") << " flag" );
78
79 return StatusCode::SUCCESS;
80}
81
82//--------------------------------------------------------------------------------
83// --------------------------------------------------------------------------------
84
86 TrigL2MuonSA::RpcFitResult& rpcFitResult,
87 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns) const
88{
89
90 //
91 for (TrigL2MuonSA::TrackPattern& itTrack : v_trackPatterns) { // loop for track candidates
92
93 if (!rpcFitResult.isSuccess and std::abs(muonRoad.extFtfMiddlePhi) > ZERO_LIMIT) { //inside-out
94 itTrack.phiMSDir = (std::abs(std::cos(muonRoad.extFtfMiddlePhi)) > ZERO_LIMIT)? std::tan(muonRoad.extFtfMiddlePhi): 0;
95 }
96 else{
97 itTrack.phiMSDir = (std::abs(std::cos(rpcFitResult.phi)) > ZERO_LIMIT)? std::tan(rpcFitResult.phi): 0;
98 }
99 itTrack.isRpcFailure = !rpcFitResult.isSuccess;
100
101 ATH_CHECK( superPointFitter(itTrack) );
102 }
103 //
104
105 return StatusCode::SUCCESS;
106}
107
108// --------------------------------------------------------------------------------
109// --------------------------------------------------------------------------------
110
112 const TrigL2MuonSA::MuonRoad& muonRoad,
113 TrigL2MuonSA::TgcFitResult& tgcFitResult,
114 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns,
115 TrigL2MuonSA::StgcHits& stgcHits,
116 TrigL2MuonSA::MmHits& mmHits) const
117{
118
119 for (TrigL2MuonSA::TrackPattern& itTrack : v_trackPatterns) { // loop for track candidates
120
121 if (tgcFitResult.isSuccess) {
122 itTrack.phiMSDir = tgcFitResult.phiDir;
123 } else {
124 if ( std::abs(muonRoad.extFtfMiddlePhi) > ZERO_LIMIT ) { //insideout
125 itTrack.phiMSDir = (std::abs(std::cos(muonRoad.extFtfMiddlePhi)) > ZERO_LIMIT)? std::tan(muonRoad.extFtfMiddlePhi): 0;
126 } else {
127 itTrack.phiMSDir = (std::abs(std::cos(p_roids->phi())) > ZERO_LIMIT)? std::tan(p_roids->phi()): 0;
128 }
129 itTrack.isTgcFailure = true;
130 }
131
132 ATH_CHECK( superPointFitter(itTrack) );
133
134 if(!m_nswStationFitter.empty())
135 ATH_CHECK( m_nswStationFitter->superPointFitter(p_roids, itTrack, stgcHits, mmHits) );
136 }
137 //
138 return StatusCode::SUCCESS;
139}
140
141// --------------------------------------------------------------------------------
142// --------------------------------------------------------------------------------
143
145 const TrigL2MuonSA::MuonRoad& muonRoad,
146 TrigL2MuonSA::TgcFitResult& tgcFitResult,
147 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns,
148 TrigL2MuonSA::StgcHits& stgcHits,
149 TrigL2MuonSA::MmHits& mmHits) const
150{
151
152 for (TrigL2MuonSA::TrackPattern& itTrack : v_trackPatterns) { // loop for track candidates
153
154 if (tgcFitResult.isSuccess) {
155 itTrack.phiMSDir = tgcFitResult.phiDir;
156 } else {
157 itTrack.phiMSDir = (std::abs(std::cos(p_roids->phi())) > ZERO_LIMIT)? std::tan(p_roids->phi()): 0;
158 itTrack.isTgcFailure = true;
159 }
160
161 ATH_CHECK( superPointFitter(itTrack, muonRoad) );
162
163 makeReferenceLine(itTrack, muonRoad);
164 ATH_CHECK( m_alphaBetaEstimate->setAlphaBeta(p_roids, tgcFitResult, itTrack, muonRoad) );
165
166 if ( itTrack.etaBin < -1 ) {
167 itTrack.etaBin = (int)((std::abs(muonRoad.extFtfMiddleEta)-1.)/0.05); // eta binning is the same as AlphaBetaEstimate
168 if(itTrack.etaBin <= -1) itTrack.etaBin = 0;
169 }
170
171 ATH_CHECK( m_ptFromAlphaBeta->setPt(itTrack,tgcFitResult) );
172
173 double exInnerA = fromAlphaPtToInn(tgcFitResult,itTrack);
174 double bw = muonRoad.bw[3][0];
175 double aw = muonRoad.aw[3][0];
176 if(std::abs(exInnerA) > ZERO_LIMIT) updateInnSP(itTrack, exInnerA, aw,bw);
177
178 if(!m_nswStationFitter.empty())
179 ATH_CHECK( m_nswStationFitter->superPointFitter(p_roids, itTrack, stgcHits, mmHits) );
180
181 }
182 //
183 return StatusCode::SUCCESS;
184}
185
186// --------------------------------------------------------------------------------
187
188// --------------------------------------------------------------------------------
189
191{
192 const unsigned int MAX_STATION = 10; // no BMG(Backup=10)
193 const float SIGMA = 0.0080;
194 const float DRIFTSPACE_LIMIT = 16.;
195 const int MIN_MDT_FOR_FIT = 3;
196
197 for (unsigned int chamber=0; chamber<MAX_STATION; chamber++) { // loop for station
198 if (chamber==9) continue;//skip BME chamber
199
200 const TrigL2MuonSA::MdtHits& mdtSegment {trackPattern.mdtSegments[chamber]};
201 if (mdtSegment.size()==0) continue;
202
203 int count {0};
204 float sigma {0.}, phim {0.}, Xor {0.}, Yor {0.}, Ymid {0.};
205 TrigL2MuonSA::PBFitResult pbFitResult;
206
207 for (const TrigL2MuonSA::MdtHitData& itMdtHit : mdtSegment) { // loop for MDT hit
208
209 if (count >= NMEAMX) continue;
210 if (itMdtHit.isOutlier) continue;
211
212 if (!count) {
213 Ymid = itMdtHit.cYmid;
214 }
215 if (!Xor) {
216 Xor = itMdtHit.R;
217 Yor = itMdtHit.Z;
218 }
219
220 phim = itMdtHit.cPhip;
221 sigma = (std::abs(itMdtHit.DriftSigma) > ZERO_LIMIT)? itMdtHit.DriftSigma: SIGMA;
222
223 if ( std::abs(itMdtHit.DriftSpace) > ZERO_LIMIT &&
224 std::abs(itMdtHit.DriftSpace) < DRIFTSPACE_LIMIT &&
225 std::abs(itMdtHit.DriftTime) > ZERO_LIMIT ) {
226
227 pbFitResult.XILIN[count] = itMdtHit.R - Xor;
228 pbFitResult.YILIN[count] = itMdtHit.Z - Yor;
229 pbFitResult.RILIN[count] = itMdtHit.DriftSpace;
230 pbFitResult.WILIN[count] = 1/(sigma*sigma);
231 pbFitResult.RESI[count] = 0.;
232
233 count++;
234 }
235 } // end loop for MDT hits
236 pbFitResult.NPOI = count;
237
238 ATH_MSG_DEBUG("... MDT hit used in fit #=" << pbFitResult.NPOI);
239 if(msgLevel(MSG::DEBUG)){
240 for(int i=0;i<pbFitResult.NPOI;i++) {
241 ATH_MSG_DEBUG("i/XILIN[i]/YILIN[i]/RILIN[i]/WILIN[i] = "
242 << i << "/" << pbFitResult.XILIN[i] << "/" << pbFitResult.YILIN[i]
243 << "/" << pbFitResult.RILIN[i] << "/" << pbFitResult.WILIN[i]);
244 }
245 }
246
247 TrigL2MuonSA::SuperPoint& superPoint {trackPattern.superPoints[chamber]};
248 superPoint.Ndigi = pbFitResult.NPOI;
249 if (count >= MIN_MDT_FOR_FIT) {
250
251 Evlfit(pbFitResult);
252
253 superPoint.Npoint = pbFitResult.NPOI;
254
255 if ((trackPattern.s_address == -1 and chamber != 0) or // Endcap
256 (trackPattern.s_address != -1 and chamber == 3)){
257
258 if (std::abs(pbFitResult.ALIN) > ZERO_LIMIT) {
259 superPoint.Z = Ymid;
260 superPoint.R = (Ymid - Yor - pbFitResult.BLIN)/pbFitResult.ALIN + Xor;
261 superPoint.Alin = 1./pbFitResult.ALIN;
262 superPoint.Blin = -pbFitResult.BLIN/pbFitResult.ALIN;
263 }
264 }
265 else{ // Barrel
266 superPoint.R = Ymid;
267 superPoint.Z = pbFitResult.ALIN*(Ymid - Xor) + pbFitResult.BLIN + Yor;
268 superPoint.Alin = pbFitResult.ALIN;
269 superPoint.Blin = pbFitResult.BLIN;
270 }
271
272 superPoint.Phim = phim;
273 superPoint.Xor = Xor;
274 superPoint.Yor = Yor;
275 superPoint.Chi2 = pbFitResult.CHI2;
276 superPoint.PChi2 = pbFitResult.PCHI2;
277 for(int i=0;i<pbFitResult.NPOI;i++) superPoint.Residual[i] = pbFitResult.RESI[i];
278 }
279
280 ATH_MSG_DEBUG("... Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
281 << chamber << "/" << trackPattern.s_address << "/" << count << "/"
282 << superPoint.R << "/" << superPoint.Z << "/" << superPoint.Alin << "/"
283 << superPoint.Blin << "/" << superPoint.Phim << "/" << superPoint.Xor << "/"
284 << superPoint.Yor << "/" << superPoint.Chi2 << "/" << superPoint.PChi2);
285 } // end loop for stations
286
287 return StatusCode::SUCCESS;
288}
289
290// --------------------------------------------------------------------------------
291//---------------------------------------------------------------------------------
293 const TrigL2MuonSA::MuonRoad& muonRoad) const
294{
295 const unsigned int MAX_STATION = 10;
296 const float SIGMA = 0.0080;
297 const float DRIFTSPACE_LIMIT = 16.;
298 const int MIN_MDT_FOR_FIT = 3;
299
300 for (unsigned int chamber=0; chamber<MAX_STATION; chamber++) { // loop for station
301 ATH_MSG_DEBUG(" superpoint fit station "<<chamber);
302
303 if(chamber== 1 || chamber == 2 || chamber ==7 || chamber == 9) continue; // only loop for endcap Inner/Middle/Outer/EE/barrel inn
304
305 TrigL2MuonSA::MdtHits& mdtSegment {trackPattern.mdtSegments[chamber]};
306 if (mdtSegment.size()==0) continue;
307
308 TrigL2MuonSA::SuperPoint& superPoint {trackPattern.superPoints[chamber]};
309 TrigL2MuonSA::PBFitResult pbFitResult{};
310
311 if (chamber==0 || chamber == 6 || chamber==8){
312 int count=0;
313 float sigma {0.}, phim {0.}, Xor {0.}, Yor {0.}, Ymid {0.};
314
315 for (const TrigL2MuonSA::MdtHitData& itMdtHit : mdtSegment) { // loop for MDT hit
316
317 if (count >= NMEAMX) continue;
318 if (itMdtHit.isOutlier) continue;
319
320 superPoint.Ndigi++;
321 if (!count) {
322 Ymid = itMdtHit.cYmid;
323 }
324 if (!Xor) {
325 Xor = itMdtHit.R;
326 Yor = itMdtHit.Z;
327 }
328
329 phim = itMdtHit.cPhip;
330 sigma = (std::abs(itMdtHit.DriftSigma) > ZERO_LIMIT)? itMdtHit.DriftSigma: SIGMA;
331
332 if ( std::abs(itMdtHit.DriftSpace) > ZERO_LIMIT &&
333 std::abs(itMdtHit.DriftSpace) < DRIFTSPACE_LIMIT &&
334 std::abs(itMdtHit.DriftTime) > ZERO_LIMIT ) {
335
336 pbFitResult.XILIN[count] = itMdtHit.R - Xor;
337 pbFitResult.YILIN[count] = itMdtHit.Z - Yor;
338 pbFitResult.RILIN[count] = itMdtHit.DriftSpace;
339 pbFitResult.WILIN[count] = 1/(sigma*sigma);
340 pbFitResult.RESI[count] = 0.;
341 count++;
342 pbFitResult.NPOI = count;
343 } else {
344 superPoint.Ndigi--;
345 }
346 } // end loop for MDT hits
347
348 ATH_MSG_DEBUG("... MDT hit used in fit #=" << pbFitResult.NPOI);
349 for(int i=0;i<pbFitResult.NPOI;i++) {
350 ATH_MSG_DEBUG("i/XILIN[i]/YILIN[i]/RILIN[i]/WILIN[i] = "
351 << i << "/" << pbFitResult.XILIN[i] << "/" << pbFitResult.YILIN[i]
352 << "/" << pbFitResult.RILIN[i] << "/" << pbFitResult.WILIN[i]);
353 }
354 if (count >= MIN_MDT_FOR_FIT) {
355 Evlfit( pbFitResult);
356
357 float bc = (Ymid - Xor);
358 float X = (pbFitResult.ALIN*bc)+pbFitResult.BLIN ;
359
360 superPoint.Npoint = pbFitResult.NPOI;
361 if (trackPattern.s_address == -1) { // Endcap
362 if (std::abs(pbFitResult.ALIN) > ZERO_LIMIT) {
363 superPoint.Z = Ymid;
364 superPoint.R = (Ymid-Yor)/pbFitResult.ALIN - pbFitResult.BLIN/pbFitResult.ALIN + Xor;
365 superPoint.Alin = 1./pbFitResult.ALIN;
366 superPoint.Blin = -pbFitResult.BLIN/pbFitResult.ALIN;
367 if (chamber==0 || chamber==8){//endcap barrel inner or BEE
368 superPoint.R = bc + Xor;
369 superPoint.Z = X + Yor;
370 superPoint.Alin = pbFitResult.ALIN;
371 superPoint.Blin = pbFitResult.BLIN;
372 }
373 }
374 }
375 superPoint.Phim = phim;
376 superPoint.Xor = Xor;
377 superPoint.Yor = Yor;
378 superPoint.Chi2 = pbFitResult.CHI2;
379 superPoint.PChi2 = pbFitResult.PCHI2;
380 for(int i=0;i<pbFitResult.NPOI;i++) superPoint.Residual[i] = pbFitResult.RESI[i];
381
382 }
383 ATH_MSG_DEBUG("...Special Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
384 << chamber << "/" << trackPattern.s_address << "/" << count << "/"
385 << superPoint.R << "/" << superPoint.Z << "/" << superPoint.Alin << "/"
386 << superPoint.Blin << "/" << superPoint.Phim << "/" << superPoint.Xor << "/"
387 << superPoint.Yor << "/" << superPoint.Chi2 << "/" << superPoint.PChi2);
388 continue;
389 }
390
391 double aw = muonRoad.aw[chamber][0];
392 double bw = muonRoad.bw[chamber][0];
393 double nrWidth = 0.;
394 unsigned int sumN = 0;
395 //chamber=3/4/5 => Endcap Inner/Middle/Outer
396 if(chamber==3) {nrWidth = m_rwidth_Endcapinn_first;}
397 if(chamber==4) {nrWidth = m_rwidth_Endcapmid_first;}
398 if(chamber==5) {nrWidth = m_rwidth_Endcapout_first;}
399
400 for (TrigL2MuonSA::MdtHitData& itMdtHit : mdtSegment) { // loop for MDT hit
401 if (std::abs(itMdtHit.DriftSpace) < m_mdt_driftspace_downlimit ||
402 std::abs(itMdtHit.DriftSpace) > m_mdt_driftspace_uplimit){
403 itMdtHit.isOutlier = 2;
404 continue;
405 }
406
407 if(itMdtHit.isOutlier > 1)continue;
408 double Z = itMdtHit.Z;
409 double R = itMdtHit.R;
410 double nbw = aw*Z + bw;
411 if (R>(nbw-nrWidth) && R<(nbw+nrWidth)){
412 itMdtHit.isOutlier = 0;
413 sumN++;
414 }
415 else {
416 itMdtHit.isOutlier = 2;
417 continue;
418 }
419 }
420 if(sumN==0) continue;
421
422 stationSPFit(mdtSegment, superPoint,pbFitResult, trackPattern.s_address,chamber,aw);
423
424 } // end loop for stations
425
426 return StatusCode::SUCCESS;
427}
428// --------------------------------------------------------------------------------
429// --------------------------------------------------------------------------------
431 TrigL2MuonSA::SuperPoint& superPoint,
432 TrigL2MuonSA::PBFitResult& pbFitResult,int s_address, int i_station,double aw) const{
433
434 TrigL2MuonSA::MdtHits::iterator itMdtHit;
435
436 unsigned int i_layer_max = 0;
437 int count;
438 float Xor, Yor, sigma, phim=0;
439
440 const unsigned int MAX_STATION = 8;
441 const float SIGMA = 0.0080;
442 const float DRIFTSPACE_LIMIT = 16.;
443 const int MIN_MDT_FOR_FIT = 3;
444
445 const unsigned int MAX_LAYER = 12;
446
447 unsigned int MdtLayerHits[MAX_STATION][MAX_LAYER];
448 std:: vector<unsigned int> MdtLayerHits_index[MAX_STATION][MAX_LAYER];
449
450 for (unsigned int i_st=0; i_st<MAX_STATION; i_st++) {
451 for (unsigned int i_ly=0; i_ly<MAX_LAYER; i_ly++) {
452 MdtLayerHits[i_st][i_ly] = 0;
453 MdtLayerHits_index[i_st][i_ly].clear();
454 }
455 }
456
457 float Chbest = 1.e25;
458 double avZ[8];
459 double avR[8];
460 double sumZ[8];
461 double sumR[8];
462 unsigned int sumN[8];
463 double nsWidth=0.;
464
465 for (unsigned int i_st=0; i_st<8; i_st++) {
466 avZ[i_st] = 0.;
467 avR[i_st] = 0.;
468 sumZ[i_st] = 0.;
469 sumR[i_st] = 0.;
470 sumN[i_st] = 0.;
471 }
472
473 for (unsigned int i_hit=0; i_hit<mdtSegment.size(); i_hit++){
474
475 if (mdtSegment.at(i_hit).isOutlier>1) continue;
476
477 double Z = mdtSegment.at(i_hit).Z;
478 double R = mdtSegment.at(i_hit).R;
479
480 sumZ[i_station] = sumZ[i_station] + Z;
481 sumR[i_station] = sumR[i_station] + R;
482 sumN[i_station]++;
483
484 if (sumN[i_station]!=0) {
485 avZ[i_station] = sumZ[i_station]/sumN[i_station];
486 avR[i_station] = sumR[i_station]/sumN[i_station];
487 }
488 }
489
490 if (sumN[i_station]==0) return;
491
492 for (unsigned int i_hit=0; i_hit<mdtSegment.size(); i_hit++) {
493
494 if (mdtSegment.at(i_hit).isOutlier>1) continue;
495
496 double Z = mdtSegment.at(i_hit).Z;
497 double R = mdtSegment.at(i_hit).R;
498
499 if (i_station==3) nsWidth = m_rwidth_Endcapinn_second;
500 if (i_station==4) nsWidth = m_rwidth_Endcapmid_second;
501 if (i_station==5) nsWidth = m_rwidth_Endcapout_second;
502
503 double nbw = aw*Z+(avR[i_station]-aw*avZ[i_station]);
504
505 if ( R>(nbw-nsWidth) && R<(nbw+nsWidth) ) {
506 mdtSegment.at(i_hit).isOutlier = 0;
507 } else {
508 mdtSegment.at(i_hit).isOutlier = 2;
509 continue;
510 }
511 }
512
513 for (unsigned int i_hit=0; i_hit<mdtSegment.size(); i_hit++) {
514
515 unsigned int i_layer =mdtSegment.at(i_hit).Layer;
516 if (mdtSegment.at(i_hit).isOutlier>1) continue;
517 if (i_layer > i_layer_max) i_layer_max = i_layer;
518
519 MdtLayerHits[i_station][i_layer]++;
520 MdtLayerHits_index[i_station][i_layer].push_back(i_hit);
521 }
522
523 std::vector<unsigned int> Ly_1st{};
524 std::vector<float> Line_A{};
525 std::vector<float> Line_B{};
526 std::vector<float> Line_Chi2{};
527 std::vector<unsigned int> Line_count{};
528 std::vector<float> Line_Xor{};
529 std::vector<float>Line_Yor{};
530 std::vector<float> Line_sum_Z{};
531 std::vector<float> Line_phim{};
532 std::vector<float> Maxlayers_A{};
533 std::vector<float> Maxlayers_B{};
534 std::vector<float> Maxlayers_Chi2{};
535 std::vector<float> Maxlayers_RESI{};
536 float Maxlayers_Phim = 0;
537 float Maxlayers_R = 0;
538 float Maxlayers_Z = 0;
539 float Maxlayers_Xor = 0;
540 float Maxlayers_Yor = 0;
541 float Maxlayers_PChi2 = 0;
542 int Maxlayers_N = 0;
543
544 for (unsigned int i_layer=0; i_layer<=i_layer_max; i_layer++) {
545
546 if (MdtLayerHits[i_station][i_layer]==0) continue;
547
548 Ly_1st.push_back(i_layer);
549 }
550
551 const int real_layer= Ly_1st.size();
552 std::vector<std::vector<unsigned int> > Ly_flg{};
553
554 for (int pr=real_layer; pr>=3; pr--) {
555
556 Ly_flg.clear();
557 Line_A.clear();
558 Line_B.clear();
559 Line_Chi2.clear();
560 Line_count.clear();
561 Line_Xor.clear();
562 Line_Yor .clear();
563 Line_sum_Z.clear();
564 Line_phim.clear();
565
566 int total_cp = 0;
567 findLayerCombination(Ly_1st, real_layer, pr,Ly_flg, total_cp);
568
569 for (unsigned int i=0;i<Ly_flg.size(); i++) {
570
571 std::vector<std::vector<int> >tID{};
572 std::vector<std::vector<int> >tIndex{};
573
574 int tube_ID[NMEAMX][2];
575 int tubeindex[NMEAMX][2];
576
577 for (int ti=0; ti<8; ti++) {
578 for (int tj=0; tj<2; tj++) {
579 tube_ID[ti][tj] = 0;
580 tubeindex[ti][tj] = 0;
581 }
582 }
583
584 for (unsigned int j=0; j<Ly_flg[i].size(); j++) {
585
586 int i_layer = Ly_flg[i][j];
587 std::vector<int> tid{};
588 std::vector<int> tindex{};
589
590 if (MdtLayerHits[i_station][i_layer]==0) continue;
591
592 float tube_1st = 999999.;
593 float tube_2nd = 999999.;
594 int layer_1st= 9999;
595 int layer_2nd = 9999;
596
597 for (unsigned int i_hit=0; i_hit< MdtLayerHits[i_station][i_layer]; i_hit++) {
598
599 unsigned int i_index = MdtLayerHits_index[i_station][i_layer].at(i_hit);
600
601 if (mdtSegment.at(i_index).isOutlier>1) continue;
602
603 float nbw3 = (mdtSegment.at(i_index).Z)*(aw) + (avR[i_station]-(aw)*avZ[i_station]) ;
604 float dis_tube = std::abs(std::abs(nbw3-mdtSegment.at(i_index).R)- mdtSegment.at(i_index).DriftSpace);
605
606 if (dis_tube<tube_1st) {
607 tube_2nd = tube_1st;
608 layer_2nd = layer_1st;
609 tube_1st = dis_tube;
610 layer_1st = i_index;
611 } else if (dis_tube<tube_2nd) {
612 tube_2nd = dis_tube;
613 layer_2nd = i_index;
614 }
615 }
616
617 if ( layer_1st != 9999 ) {
618 mdtSegment.at(layer_1st).isOutlier = 0;
619 tid.push_back(1);
620 tindex.push_back(layer_1st);
621 }
622
623 if ( layer_2nd != 9999 ) {
624 mdtSegment.at(layer_2nd).isOutlier = 1;
625 tid.push_back(1);
626 tindex.push_back(layer_2nd);
627 }
628
629 tID.push_back(tid);
630 tIndex.push_back(tindex);
631 }
632
633 for (unsigned int ti=0; ti<tID.size();ti++) {
634 for (unsigned int tj=0; tj<tID[ti].size();tj++) {
635 tube_ID[ti][tj] = tID[ti][tj];
636 tubeindex[ti][tj] = tIndex[ti][tj];
637 }
638 }
639
640 std::vector<int> isg;
641 std::vector<int> hitarray;
642 int sumid;
643 int ntry = (int)floor(pow(2.,pr))-1;
644
645 for (int ntryi=0; ntryi<=ntry; ntryi++) {
646
647 isg.clear();
648 hitarray.clear();
649 sumid = 1;
650
651 for (int ntryj=1; ntryj<=pr; ntryj++) {
652 int yhit = 0;
653 int Isg = (ntryi&(int)pow(2.,ntryj-1))? 1 : 0;
654 isg.push_back(Isg);
655 if (tube_ID[ntryj-1][Isg] != 0) yhit = 1;
656 sumid = sumid * yhit;
657 }
658
659 if (sumid==1) {
660 for (unsigned int tt=0;tt<isg.size(); tt++) {
661 int tindex = tubeindex[tt][isg[tt]];
662 hitarray.push_back(tindex);
663 }
664 }
665
666 count = 0;
667 Xor = 0.;
668 Yor = 0.;
669
670 float sum_Z_used = 0.;
671 float sum_R_used = 0.;
672
673 if (hitarray.size()==0) continue;
674
675 for (itMdtHit=mdtSegment.begin(); itMdtHit!=mdtSegment.end(); ++itMdtHit) { // loop for MDT hit
676
677 int hit_index = std::distance(mdtSegment.begin(),itMdtHit);
678
679 if(mdtSegment.at(hit_index).isOutlier>1) continue;
680
681 if (count >= NMEAMX) continue;
682
683 int fd=0;
684
685 for (unsigned int j=0; j<hitarray.size(); j++) {
686
687 if (hitarray[j]==hit_index) {
688 fd=1;
689 break;
690 }
691 }
692
693 if (fd==0) continue;
694
695 superPoint.Ndigi++;
696
697 if (!Xor) {
698 Xor = itMdtHit->R;
699 Yor = itMdtHit->Z;
700 }
701
702 phim = itMdtHit->cPhip;
703 sigma = (std::abs(itMdtHit->DriftSigma) > ZERO_LIMIT)? itMdtHit->DriftSigma: SIGMA;
704
705 if ( std::abs(itMdtHit->DriftSpace) > ZERO_LIMIT &&
706 std::abs(itMdtHit->DriftSpace) < DRIFTSPACE_LIMIT &&
707 std::abs(itMdtHit->DriftTime) > ZERO_LIMIT ) {
708
709 pbFitResult.XILIN[count] = itMdtHit->R - Xor;
710 pbFitResult.YILIN[count] = itMdtHit->Z - Yor;
711 pbFitResult.RILIN[count] = itMdtHit->DriftSpace;
712 pbFitResult.WILIN[count] = 1/(sigma*sigma);
713 pbFitResult.RESI[count] = 0.;
714
715 count++;
716 pbFitResult.NPOI = count;
717
718 sum_Z_used = sum_Z_used + itMdtHit->Z;
719 sum_R_used = sum_R_used + itMdtHit->R;
720 } else {
721 superPoint.Ndigi--;
722 }
723 } // end loop for MDT hits
724
725 ATH_MSG_DEBUG("... MDT hit used in fit #=" << pbFitResult.NPOI);
726 for (int i=0;i<pbFitResult.NPOI;i++) {
727 ATH_MSG_DEBUG("i/XILIN[i]/YILIN[i]/RILIN[i]/WILIN[i] = "
728 << i << "/" << pbFitResult.XILIN[i] << "/" << pbFitResult.YILIN[i]
729 << "/" << pbFitResult.RILIN[i] << "/" << pbFitResult.WILIN[i]);
730 }
731
732 if (count >= MIN_MDT_FOR_FIT) {
733 Circles(pbFitResult.NPOI,pbFitResult.XILIN,pbFitResult.YILIN,pbFitResult.RILIN,pbFitResult.WILIN,
734 pbFitResult.ALIN,pbFitResult.BLIN,pbFitResult.CHI2,
735 pbFitResult.PCHI2, pbFitResult.SlopeCand, pbFitResult.InterceptCand, pbFitResult.Chi2Cand);
736
737
738 for (int cand=0; cand<6; cand++) {
739
740 if (std::abs(pbFitResult.SlopeCand[cand]) > ZERO_LIMIT) {
741 Line_A.push_back(1/pbFitResult.SlopeCand[cand]);
742 Line_B.push_back(-pbFitResult.InterceptCand[cand]/pbFitResult.SlopeCand[cand]-Yor/pbFitResult.SlopeCand[cand]+Xor);
743 Line_Chi2.push_back(pbFitResult.Chi2Cand[cand]);
744 Line_count.push_back(pbFitResult.NPOI);
745 Line_Xor.push_back(Xor);
746 Line_Yor .push_back(Yor);
747 Line_sum_Z.push_back(sum_Z_used/count);
748 Line_phim.push_back(phim);
749 }
750 }
751 }
752 }
753 }//end one of cp
754
755 if (Line_Chi2.size()==0) continue;
756
757 std::multimap<float, int>chi_map{};
758 std::vector<float> t_A;
759 std::vector<float> t_B;
760 std::vector<float> t_Chi2;
761 std::vector<float> t_count;
762 std::vector<float> t_Xor;
763 std::vector<float> t_Yor;
764 std::vector<float> t_sum_Z;
765 std::vector<float> t_phim;
766
767 t_A.clear();
768 t_B.clear();
769 t_Chi2.clear();
770 t_count.clear();
771 t_Xor.clear();
772 t_Yor.clear();
773 t_sum_Z.clear();
774 t_phim.clear();
775
776 for (unsigned int ir=0; ir<Line_Chi2.size(); ir++) chi_map.insert(std::make_pair(Line_Chi2.at(ir), ir));
777
778 for (std::multimap<float, int>::iterator jt = chi_map.begin(); jt != chi_map.end(); ++jt) {
779 t_A.push_back(Line_A.at(jt->second));
780 t_B.push_back(Line_B.at(jt->second));
781 t_Chi2.push_back(Line_Chi2.at(jt->second));
782 t_count.push_back(Line_count.at(jt->second));
783 t_Xor.push_back(Line_Xor.at(jt->second));
784 t_Yor.push_back(Line_Yor.at(jt->second));
785 t_sum_Z.push_back(Line_sum_Z.at(jt->second));
786 t_phim.push_back(Line_phim.at(jt->second));
787 if(pr==real_layer){//save max layers information
788 Maxlayers_A.push_back(Line_A.at(jt->second));
789 Maxlayers_B.push_back(Line_B.at(jt->second));
790 Maxlayers_Chi2.push_back(Line_Chi2.at(jt->second));
791 }
792 }
793
794 superPoint.Npoint = t_count[0];//pbFitResult.NPOI;
795 if(i_station==4 && pr==real_layer){
796 Maxlayers_Z = t_sum_Z[0];
797 Maxlayers_R = t_A[0]*t_sum_Z[0]+t_B[0];
798 Maxlayers_Phim = t_phim[0];
799 Maxlayers_Xor = t_Xor[0];
800 Maxlayers_Yor = t_Yor[0];
801 Maxlayers_PChi2 = pbFitResult.PCHI2;
802 Maxlayers_N = t_count[0];
803 }
804
805 if (s_address == -1) { // Endcap
806
807 if (std::abs(t_A[0]) > ZERO_LIMIT ) {
808 superPoint.Z = t_sum_Z[0];
809 superPoint.R = t_A[0]*t_sum_Z[0]+t_B[0];
810 superPoint.Alin =t_A[0];
811 superPoint.Blin =t_B[0];
812 }
813
814 superPoint.Phim = t_phim[0];
815 superPoint.Xor = t_Xor[0];
816 superPoint.Yor = t_Yor[0];
817 superPoint.Chi2 = t_Chi2[0];
818 superPoint.PChi2 = pbFitResult.PCHI2;
819
820 for (int i=0;i<pbFitResult.NPOI;i++) superPoint.Residual[i] = pbFitResult.RESI[i];
821
822 for (int cand=0; cand<6; cand++) {
823 if (std::abs(t_A[cand]) > ZERO_LIMIT ) {
824 superPoint.SlopeCand[cand] = t_A[cand];
825 superPoint.InterceptCand[cand] = t_B[cand];
826 superPoint.Chi2Cand[cand] = t_Chi2[cand];
827 }
828 }
829 }
830
831 Chbest=t_Chi2[0];
832
833 if (real_layer>3) {
834 if ((i_station == 3 || i_station == 5) && pr==4 && Chbest > m_endcapmid_mdt_chi2_limit) {
835
836 superPoint.Z =0.;
837 superPoint.R =0.;
838 superPoint.Alin=0.;
839 superPoint.Blin=0.;
840
841 for (int cand=0; cand<6; cand++) {
842 superPoint.SlopeCand[cand] = 0.;
843 superPoint.InterceptCand[cand] = 0.;
844 superPoint.Chi2Cand[cand] = 0.;
845 }
846 return;
847 }
848 }
849
850 if (Chbest<m_endcapmid_mdt_chi2_limit){
851
852 ATH_MSG_DEBUG("... Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
853 << i_station << "/" << s_address << "/" << count << "/"
854 << superPoint.R << "/" << superPoint.Z << "/" << superPoint.Alin << "/"
855 << superPoint.Blin << "/" << superPoint.Phim << "/" << superPoint.Xor << "/"
856 << superPoint.Yor << "/" << superPoint.Chi2 << "/" << superPoint.PChi2);
857
858 break;//jump out all cp
859 }else{
860 if(i_station==4 && Maxlayers_A.size()>0){
861 superPoint.Npoint = Maxlayers_N;
862 superPoint.Z = Maxlayers_Z;
863 superPoint.R = Maxlayers_R;
864 superPoint.Alin =Maxlayers_A[0];
865 superPoint.Blin =Maxlayers_B[0];
866 superPoint.Phim = Maxlayers_Phim;
867 superPoint.Xor = Maxlayers_Xor;
868 superPoint.Yor = Maxlayers_Yor;
869 superPoint.Chi2 = Maxlayers_Chi2[0];
870 superPoint.PChi2 = Maxlayers_PChi2;
871 for (int cand=0; cand<6; cand++) {
872 if (std::abs(Maxlayers_A[cand]) > ZERO_LIMIT ) {
873 superPoint.SlopeCand[cand] = Maxlayers_A[cand];
874 superPoint.InterceptCand[cand] = Maxlayers_B[cand];
875 superPoint.Chi2Cand[cand] = Maxlayers_Chi2[cand];
876 }
877 }
878 }
879 }
880 }//end all cp
881
882 return;
883}
884
885// --------------------------------------------------------------------------------
886//-----------------------------------------------------------------------------------
888 const TrigL2MuonSA::MuonRoad& muonRoad) const{
889
890 const unsigned int MAX_STATION = 8;
891 float aw[8];
892 float spZ[8];
893 float spR[8];
894 float A_cand[8][6];
895 float B_cand[8][6];
896 float Chi2_cand[8][6];
897 float spA[8];
898 float spB[8];
899 float spChi2[8];
900 const int middle = 4;
901 const int outer = 5;
902
903 for (unsigned int i_station=4; i_station<MAX_STATION; i_station++) {
904
905 aw[i_station]=0.;
906 spZ[i_station]=0.;
907 spR[i_station]=0.;
908 spA[i_station]=0.;
909 spB[i_station]=0.;
910 spChi2[i_station]=0.;
911
912 for (unsigned int ci=0; ci<NCAND; ci++) {
913 A_cand[i_station][ci] =0.;
914 B_cand[i_station][ci] =0.;
915 Chi2_cand[i_station][ci] =0.;
916 }
917
918 if (i_station<4 || i_station>5) continue; // only loop for endcap Inner/Middle/Outer
919
920 TrigL2MuonSA::SuperPoint& superPoint {trackPattern.superPoints[i_station]};
921 aw[i_station] = muonRoad.aw[i_station][0];
922 spZ[i_station] = superPoint.Z;
923 spR[i_station] = superPoint.R;
924 spA[i_station] = superPoint.Alin;
925 spB[i_station] = superPoint.Blin;
926
927 for (unsigned int cand=0; cand<NCAND; cand++) {
928 A_cand[i_station][cand] = superPoint.SlopeCand[cand];
929 B_cand[i_station][cand] = superPoint.InterceptCand[cand];
930 Chi2_cand[i_station][cand] = superPoint.Chi2Cand[cand];
931 }
932 }
933
934 float test_diff = 1.e25;
935 float best_diff = 1.e25;
936 float rmatched = 0.;
937 float match_midA = 0.;
938 float match_outA = 0.;
939
940 if (std::abs(A_cand[middle][0]) < ZERO_LIMIT && std::abs(A_cand[middle][1]) < ZERO_LIMIT) {
941 spZ[middle] = 0.;
942 spR[middle] = 0.;
943 spChi2[middle] = 0.;
944 }
945
946 if (std::abs(A_cand[outer][0]) < ZERO_LIMIT && std::abs(A_cand[outer][1]) < ZERO_LIMIT) {
947 spZ[outer] = 0.;
948 spR[outer] = 0.;
949 spChi2[outer] = 0.;
950 }
951
952 if (std::abs(A_cand[middle][0]) > ZERO_LIMIT && std::abs(A_cand[outer][0]) < ZERO_LIMIT) {
953
954 best_diff = 1.e25;
955 test_diff = 1.e25;
956
957 for (int m=0; m<6; m++) {
958
959 test_diff = std::abs(A_cand[middle][m] - aw[middle]);
960
961 if (test_diff<best_diff) {
962 best_diff = test_diff;
963 rmatched = A_cand[middle][m];
964 spB[middle] = B_cand[middle][m];
965 spChi2[middle] = Chi2_cand[middle][m];
966 spR[middle] = rmatched * spZ[middle] + spB[middle];
967
968 }
969 }
970 }
971
972 if(std::abs(A_cand[outer][1]) > ZERO_LIMIT && std::abs(A_cand[outer][0]) > ZERO_LIMIT && std::abs(spA[outer]) > ZERO_LIMIT && std::abs(spZ[outer]) > ZERO_LIMIT && std::abs(spR[outer]) > ZERO_LIMIT &&
973 std::abs(A_cand[middle][1]) > ZERO_LIMIT && std::abs(A_cand[middle][0]) > ZERO_LIMIT && std::abs(spA[middle]) > ZERO_LIMIT && std::abs(spZ[middle]) > ZERO_LIMIT && std::abs(spR[middle]) > ZERO_LIMIT){
974
975 float sp_line = 0.;
976 if(std::abs(spZ[outer]-spZ[middle]) > ZERO_LIMIT) sp_line = (spR[outer]-spR[middle])/(spZ[outer]-spZ[middle]);
977
978 for (int t=0; t<2; ++t) {
979 best_diff = 1.e25;
980 test_diff = 1.e25;
981 if (std::abs(sp_line)> ZERO_LIMIT) {
982 for (int i=0; i<6; ++i) {
983 if (t==0) test_diff = std::abs(A_cand[middle][i] - sp_line);
984 else if (t==1) test_diff = std::abs(A_cand[outer][i] - sp_line);
985 if (test_diff<best_diff) {
986 best_diff = test_diff;
987 if (t==0) {
988 match_midA = A_cand[middle][i];
989 spB[middle] = B_cand[middle][i];
990 spChi2[middle] = Chi2_cand[middle][i];
991 spR[middle] = match_midA * spZ[middle] + spB[middle];
992 } else if(t==1) {
993 match_outA = A_cand[outer][i];
994 spB[outer] = B_cand[outer][i];
995 spChi2[outer] = Chi2_cand[outer][i];
996 spR[outer] = match_outA * spZ[outer] + spB[outer];
997 }
998 }
999 }
1000 }
1001 }
1002 }
1003
1004 if (std::abs(spA[middle]) > ZERO_LIMIT) {
1005 if (std::abs(match_midA) > ZERO_LIMIT) {
1006 spA[middle] = match_midA;
1007 } else if (std::abs(rmatched) > ZERO_LIMIT) {
1008 spA[middle] = rmatched;
1009 }
1010
1011 if (std::abs(match_outA) > ZERO_LIMIT) spA[outer] = match_outA;
1012 }
1013
1014 for (unsigned int i_station=4; i_station<MAX_STATION; i_station++) {
1015 if (i_station<4 || i_station>5) continue; // only loop for endcap Inner/Middle/Outer
1016 TrigL2MuonSA::SuperPoint& superPoint {trackPattern.superPoints[i_station]};
1017 if(std::abs(spA[i_station]) > ZERO_LIMIT){
1018 superPoint.Alin =spA[i_station];
1019 superPoint.Blin =spB[i_station];
1020 superPoint.Chi2 =spChi2[i_station];
1021 }
1022 }
1023 return;
1024}
1025
1026//--------------------------------------------------------------------------------------
1027
1029 const TrigL2MuonSA::TrackPattern& trackPattern) const
1030{
1031 float MiddleSlope = 0;
1032 float OuterSlope = 0;
1033
1034 for (int i_station=4; i_station<6; i_station++) {
1035
1036 int chamberID = -1;
1037
1038 if ( i_station == 4 ) chamberID = xAOD::L2MuonParameters::Chamber::EndcapMiddle;
1039 if ( i_station == 5 ) chamberID = xAOD::L2MuonParameters::Chamber::EndcapOuter;
1040
1041 const TrigL2MuonSA::SuperPoint& superPoint {trackPattern.superPoints[chamberID]};
1042
1043 if ( superPoint.Npoint > 2 && superPoint.R > 0.) {
1044 if ( i_station==4 ) {
1045 MiddleSlope = superPoint.Alin;
1046 }
1047 if ( i_station==5 ) {
1048 OuterSlope = superPoint.Alin;
1049 }
1050 }
1051 }
1052
1053 double mdtpT = std::abs(tgcFitResult.tgcPT);
1054 double alpha_pt = std::abs(trackPattern.ptEndcapAlpha);
1055
1056 if (std::abs(MiddleSlope) > ZERO_LIMIT && std::abs(OuterSlope) > ZERO_LIMIT) {
1057 mdtpT = alpha_pt;
1058 } else if (std::abs(tgcFitResult.tgcPT)>=8.0 && std::abs(MiddleSlope) > ZERO_LIMIT) {
1059 mdtpT = alpha_pt;
1060 }
1061
1062 mdtpT = (std::abs(tgcFitResult.tgcPT)>1e-5)? mdtpT*(tgcFitResult.tgcPT/std::abs(tgcFitResult.tgcPT)) : 0;
1063 double etaMiddle = (tgcFitResult.tgcMid1[3])? tgcFitResult.tgcMid1[0] : tgcFitResult.tgcMid2[0];
1064 double phiMiddle = (tgcFitResult.tgcMid1[3])? tgcFitResult.tgcMid1[1] : tgcFitResult.tgcMid2[1];
1065 double eta;
1066 double sigma_eta;
1067 double extrInnerEta = 0;
1068 double naw = 0;
1070 muonSA->makePrivateStore();
1071 muonSA->setSAddress(-1);
1072 muonSA->setPt(mdtpT);
1073 muonSA->setEtaMS(etaMiddle);
1074 muonSA->setPhiMS(phiMiddle);
1075 muonSA->setRMS(0.);
1076 muonSA->setZMS(0.);
1077 double phi;
1078 double sigma_phi;
1079 double theta = 0.;
1080 StatusCode sc = m_backExtrapolator->give_eta_phi_at_vertex(muonSA, eta,sigma_eta,phi,sigma_phi,0.);
1081
1082 if (sc.isSuccess() ){
1083 extrInnerEta = eta;
1084 } else {
1085 extrInnerEta = etaMiddle;
1086 }
1087
1088 delete muonSA;
1089
1090 if (std::abs(extrInnerEta) > ZERO_LIMIT) {
1091 theta = std::atan(std::exp(-std::abs(extrInnerEta)))*2.;
1092 naw = std::tan(theta)*(std::abs(extrInnerEta)/extrInnerEta);
1093 }
1094
1095 return naw;
1096
1097}
1098// --------------------------------------------------------------------------------
1099//--------------------------------------------------------------------------------
1101 double &aw,
1102 double &tgc_aw,
1103 double &bw) const
1104{
1105 double nrWidth = m_rwidth_Endcapinn_first;
1106 unsigned int sumN[8];
1107
1108 for (unsigned int i_st=0; i_st<8;i_st++) {
1109 sumN[i_st] = 0;
1110 }
1111
1112 TrigL2MuonSA::PBFitResult pbFitResult;
1113 const int i_station = 3;//endcap inner
1114 int chamberID = i_station;
1115
1116 TrigL2MuonSA::MdtHits& mdtSegment {trackPattern.mdtSegments[chamberID]};
1117 TrigL2MuonSA::SuperPoint& superPoint {trackPattern.superPoints[chamberID]};
1118
1119 if (mdtSegment.size()==0) return;
1120
1121 for (TrigL2MuonSA::MdtHitData& itMdtHit : mdtSegment) { // loop for MDT hit
1122 if (std::abs(itMdtHit.DriftSpace) < m_mdt_driftspace_downlimit ||
1123 std::abs(itMdtHit.DriftSpace) > m_mdt_driftspace_uplimit){
1124 itMdtHit.isOutlier = 2;
1125 continue;
1126 }
1127
1128 if (itMdtHit.isOutlier > 1) continue;
1129
1130 double Z = itMdtHit.Z;
1131 double R = itMdtHit.R;
1132 double nbw = tgc_aw*Z + bw;
1133
1134 if (R>(nbw-nrWidth) && R<(nbw+nrWidth)) {
1135 itMdtHit.isOutlier = 0;
1136 sumN[i_station]++;
1137 } else {
1138 itMdtHit.isOutlier = 2;
1139 continue;
1140 }
1141 }
1142
1143 if (sumN[i_station]==0) return;
1144
1145 stationSPFit(mdtSegment, superPoint,pbFitResult, trackPattern.s_address,i_station, aw);
1146
1147 float df=1.e25;
1148
1149 for (int cand=0; cand<NCAND; cand++) {
1150 float ds=std::abs(superPoint.SlopeCand[cand]-aw);
1151 if (ds<df) {
1152 df=ds;
1153 superPoint.Alin = superPoint.SlopeCand[cand];
1154 superPoint.Blin = superPoint.InterceptCand[cand];
1155 superPoint.Chi2 = superPoint.Chi2Cand[cand];
1156 }
1157 }
1158}
1159
1160// --------------------------------------------------------------------------------
1161// --------------------------------------------------------------------------------
1162
1164 int n,
1165 int r,
1166 std::vector<std::vector<unsigned int> > &c,
1167 int &nr) const
1168{
1169 std::vector<unsigned int> b(r,0);
1170
1171 findSubLayerCombination(a,n,r,b,0,r,c,nr);
1172
1173 return;
1174}
1175
1176// --------------------------------------------------------------------------------
1177// --------------------------------------------------------------------------------
1178
1180 int n,
1181 int r,
1182 std::vector<unsigned int> &b,
1183 int index,
1184 int num,
1185 std::vector<std::vector<unsigned int> > &c,
1186 int &nr) const
1187{
1188 for (int i=index; i<n-num+1; i++) {
1189
1190 b[r-num] = a[i];
1191 std::vector<unsigned int> t;
1192 t.clear();
1193
1194 if (num==1) {
1195
1196 for (int j = 0;j<r; j++) t.push_back(b[j]);
1197
1198 c.push_back(t);
1199 nr++;
1200
1201 } else {
1202
1203 findSubLayerCombination(a,n,r,b,i+1,num-1,c,nr);
1204 }
1205 }
1206
1207 return;
1208}
1209
1210// --------------------------------------------------------------------------------
1211// --------------------------------------------------------------------------------
1212
1213/*==============================================================================
1214
1215 Call CIRCLES to fit a line, then rejects the worst point(s)
1216
1217 IFIT = 0 ok, = 1 failure
1218 IFLA = 0 normal, = 1 do not try to exclude points
1219
1220==============================================================================*/
1221
1223{
1224 Circles(pbFitResult.NPOI,pbFitResult.XILIN,pbFitResult.YILIN,pbFitResult.RILIN,pbFitResult.WILIN,
1225 pbFitResult.ALIN,pbFitResult.BLIN,pbFitResult.CHI2,pbFitResult.PCHI2);
1226
1227 if(pbFitResult.CHI2<=ZERO_LIMIT) return;
1228
1229 float Xnor = 1. / std::sqrt(1. + pbFitResult.ALIN * pbFitResult.ALIN);
1230
1231 for(int j=0;j<pbFitResult.NPOI;j++) {
1232
1233 float distj = (pbFitResult.ALIN * pbFitResult.XILIN[j] + pbFitResult.BLIN - pbFitResult.YILIN[j]) * Xnor;
1234 float rlin = (distj>=0.) ? pbFitResult.RILIN[j] : -pbFitResult.RILIN[j];
1235 pbFitResult.RESI[j] = distj - rlin;
1236
1237 }
1238}
1239
1240// --------------------------------------------------------------------------------
1241// --------------------------------------------------------------------------------
1242
1243/*==============================================================================
1244
1245 Computes the best line tangent to many circles or points
1246
1247 Version 1 - P.B. 960715 - First attempt
1248 Version 2 - P.B. 960910 - New algorithm for choice of the signs
1249
1250 y = A x + b
1251
1252 Input :
1253 =======
1254
1255 NMEAS : number of meas
1256 Xi, Yi : x, y of the centers (i=1...NMEAS) (i.e. wire of a mdt)
1257 RRi : distance line-center (i.e. measured drift distance)
1258 Wi : weight for the i-th meas (= 1/sigma**2)
1259
1260 Output :
1261 ========
1262
1263 A,B : line coefficients (y = Ax + b)
1264 CHI2 : Chi square = Sum Wi * (RRi - dist) ** 2 (error if .lt.0)
1265 PCHI2 : P(chisquare,d.o.f.)
1266==============================================================================*/
1267
1269 const std::array<float, NMEAMX>& XI,
1270 const std::array<float, NMEAMX>& YI,
1271 const std::array<float, NMEAMX>& RI,
1272 const std::array<float, NMEAMX>& WI,
1273 float& A,float& B, float& Chi2,float& Pchi2) const
1274{
1275 if (Nmeas<=2||Nmeas>=NMEAMX+1) return;
1276
1277 //------ First attempt, try with a line through the centers --------
1278 float SAA,SBB,SAB; //SAA,SBB,SAB, not used here
1279 Xline(XI,YI,WI,Nmeas,A,B,SAA,SBB,SAB,Chi2);
1280
1281 //------ Then choose 4 best hits and try all possible combinations (+Ri/-Ri) to identify the best segment candidate
1282 const float WIlim = 0.1 * (*std::max_element(WI.begin(), WI.begin() + Nmeas));
1283 const int Ngood = std::count_if(WI.begin(), WI.begin() + Nmeas, [WIlim](const float& w) {
1284 return w >= WIlim;
1285 }); //we count the good hits
1286
1287 std::vector<int> bestIdx{};
1288 for (int j=0;j<Nmeas;j++) { //we save the index of the 3/4 best hits
1289 if (WI[j]>=WIlim || Ngood<=3) { //if we have at least 4 good hits we use them, otherwise we use all
1290 if (bestIdx.size() < 4){
1291 bestIdx.push_back(j);
1292 }
1293 else{
1294 bestIdx.at(bestIdx.size()-2) = bestIdx.back();
1295 bestIdx.back() = j;
1296 }
1297 }
1298 }
1299 const int Nbest = bestIdx.size(); // is 3 or 4
1300
1301 const int Ntry = (1 << Nbest) - 1; // all possible combinations (+Ri/-Ri) of the best hits, i.e. Nbest=4 -> Ntry=15.
1302 std::array<float,NMEAMX> RRi{};
1303 float Abest{0.}, Bbest {0.}, CHbest{1.e25};
1304 for (int j=0;j<=Ntry;j++) { //loop over all combinations E.g. 0010 meas 2-th has -R and the other +R
1305 for (int k=0;k<Nbest;k++) { //loop over best hits
1306 const int Isig = (j & (1<<k)) ? 1 : 0; //is 1 if the k-th hit contributes as -R, otherwise +R
1307 RRi[bestIdx[k]] = (Isig==1) ? -RI[bestIdx[k]] : RI[bestIdx[k]];
1308 }
1309 float Aj = A;
1310 float Bj = B;
1311 float chi2j = -1;
1312 Circfit(Nmeas,XI,YI,RRi,WI,Aj,Bj,chi2j,&bestIdx);
1313
1314 if (chi2j>=0.0&&chi2j<=CHbest) {
1315 Abest = Aj;
1316 Bbest = Bj;
1317 CHbest = chi2j;
1318 }
1319 }
1320
1321 // ... and finally with all the points
1322 Chi2 = -1.;
1323 A=Abest;
1324 B=Bbest;
1325 Circfit(Nmeas,XI,YI,RI,WI,A,B,Chi2);
1326
1327 if (Chi2>=0.0) {
1328 Pchi2 = TMath::Prob(Chi2, Nmeas - 2);
1329 }
1330}
1331
1332// --------------------------------------------------------------------------------
1333// --------------------------------------------------------------------------------
1334
1335/*==============================================================================
1336
1337 Auxiliary to circles (see) - It makes the actual computations
1338
1339 NB - A and B MUST contain on input a first approximation of
1340 their values
1341
1342==============================================================================*/
1343
1345 const std::array<float, NMEAMX>& XI,
1346 const std::array<float, NMEAMX>& YI,
1347 const std::array<float, NMEAMX>& RI,
1348 const std::array<float, NMEAMX>& WI,
1349 float& A,float& B,float& Chi2, std::vector<int>* idx_vec) const
1350{
1351 const bool use_all {!idx_vec}; //if idx_vec=nullptr we use all hits
1352 std::vector<int> temp{};
1353 if(use_all){
1354 temp.resize(Nmeas);
1355 std::iota(temp.begin(), temp.end(), 0);
1356 idx_vec = &temp; //if we use all hits, idx_vec = 0,1,2,..Nmeas-1
1357 }
1358
1359 std::array<float,NMEAMX> XX{},YY{};
1360 int Niter{0};
1361 float Test, Toll {.1};
1362
1363 // Many iterations ...
1364 do {
1365 Niter++;
1366 float Xnor = 1. / std::sqrt(1. + A * A);
1367 float Aold = A;
1368 float Bold = B;
1369
1370 if(use_all){
1371 for(const int& j : *idx_vec) {
1372 const int Epsi {(A * XI[j] + B - YI[j]>=0.) ? 1 : -1};
1373 XX[j] = XI[j] - Epsi * Xnor * std::abs(RI[j]) * A;
1374 YY[j] = YI[j] + Epsi * Xnor * std::abs(RI[j]);
1375 }
1376 }
1377 else {
1378 for(const int& j : *idx_vec){
1379 XX[j] = XI[j] - Xnor * RI[j] * A;
1380 YY[j] = YI[j] + Xnor * RI[j];
1381 }
1382 }
1383 float SAA,SAB,SBB;
1384 Xline(XX,YY,WI,Nmeas,A,B,SAA,SBB,SAB,Chi2, idx_vec);
1385 if(Chi2<=0.) break;
1386 Test = ((Aold-A)*(Aold-A))/ SAA + ((Bold-B)*(Bold-B))/ SBB;
1387
1388 } while(Test>=Toll&&Niter<=20);
1389}
1390
1391// --------------------------------------------------------------------------------
1392// --------------------------------------------------------------------------------
1393
1394/*==============================================================================
1395 A simple linear fit : y = A x + B (see PDG 94, 17.20-25)
1396
1397 W = weights ( = 1./err**2)
1398==============================================================================*/
1399
1400void TrigL2MuonSA::MuFastStationFitter::Xline (const std::array<float, NMEAMX>& X,
1401 const std::array<float, NMEAMX>& Y,
1402 const std::array<float, NMEAMX>& W,
1403 const int NP,
1404 float& A,float& B,float& SAA,float& SBB,float& SAB,float& Square, std::vector<int>* idx_vec) const
1405{
1406 std::vector<int> temp{};
1407 if(!idx_vec){
1408 temp.resize(NP);
1409 std::iota(temp.begin(), temp.end(), 0);
1410 idx_vec = &temp; //if we use all hits, idx_vec = 0,1,2,..Nmeas-1
1411 }
1412
1413 float S1{0.},SX{0.},SY{0.},SXX{0.},SXY{0.},SYY{0.};
1414
1415 for(const int& j : *idx_vec) {
1416 S1 = S1 + W[j];
1417 SX = SX + W[j] * X[j];
1418 SY = SY + W[j] * Y[j];
1419 SXX = SXX + W[j] * X[j] * X[j];
1420 SXY = SXY + W[j] * X[j] * Y[j];
1421 SYY = SYY + W[j] * Y[j] * Y[j];
1422 }
1423
1424 float Deter = S1 * SXX - SX * SX;
1425
1426 if (std::abs(Deter) > ZERO_LIMIT) {
1427 A = (S1 * SXY - SX * SY) / Deter;
1428 B = (SY * SXX - SX * SXY) / Deter;
1429 SAA = S1 / Deter;
1430 SBB = SXX / Deter;
1431 SAB = - SX / Deter;
1432 }
1433 else {
1434 A= (S1 * SXY - SX * SY > 0.) ? 9.e+5 : -9.e+5;
1435 B = SY/S1 - SX/S1 * A;
1436 SAA = A;
1437 SBB = A;
1438 SAB = - A;
1439 }
1440 Square = 0.;
1441 for(const int& j : *idx_vec) {
1442 float DY =(Y[j] - A * X[j] - B)/std::sqrt(1 + A * A);
1443 Square = Square + W[j] * DY * DY;
1444 }
1445
1446}
1447
1448// --------------------------------------------------------------------------------
1449// --------------------------------------------------------------------------------
1451 const std::array<float, NMEAMX>& XI,
1452 const std::array<float, NMEAMX>& YI,
1453 const std::array<float, NMEAMX>& RI,
1454 const std::array<float, NMEAMX>& WI,
1455 float& A,float& B,float& Chi2,float& Pchi2,
1456 std::array<float, NCAND>& SlopeCand,
1457 std::array<float, NCAND>& InterceptCand,
1458 std::array<float, NCAND>& Chi2Cand) const
1459{
1460 if (Nmeas<=2||Nmeas>=NMEAMX+1) return;
1461
1462 //------ First attempt, try with a line through the centers --------
1463 float A0,B0, SAA,SBB,SAB,Square; //SAA,SBB,SAB,Square not used here
1464 Xline(XI,YI,WI,Nmeas,A0,B0,SAA,SBB,SAB,Square);
1465
1466 //------ Then choose 4 best hits and try all possible combinations (+Ri/-Ri) to identify the best segment candidate
1467 const float WIlim = 0.1 * (*std::max_element(WI.begin(), WI.begin() + Nmeas));
1468 const int Ngood = std::count_if(WI.begin(), WI.begin() + Nmeas, [WIlim](const float& w) {
1469 return w >= WIlim;
1470 }); //we count the good hits
1471
1472 std::vector<int> bestIdx{};
1473 for (int j=0;j<Nmeas;j++) { //we save the index of the 3/4 best hits
1474 if (WI[j]>=WIlim || Ngood<=3) { //if we have at least 4 good hits we use them, otherwise we use all
1475 if (bestIdx.size() < 4){
1476 bestIdx.push_back(j);
1477 }
1478 else{
1479 bestIdx.at(bestIdx.size()-2) = bestIdx.back();
1480 bestIdx.back() = j;
1481 }
1482 }
1483 }
1484 const int Nbest = bestIdx.size(); // is 3 or 4
1485
1486 std::vector<float> st_chi2{};
1487 std::vector<float> st_A{};
1488 std::vector<float> st_B{};
1489
1490 for (int i=0; i<NCAND; i++) {
1491 SlopeCand[i] = 0.;
1492 InterceptCand[i] = 0.;
1493 Chi2Cand[i] = 0.;
1494 }
1495
1496 const int Ntry = (1 << Nbest) - 1; // all possible combinations (+Ri/-Ri) of the best hits, i.e. Nbest=4 -> Ntry=15.
1497 std::array<float,NMEAMX> RRi{};
1498 float CHbest{1.e25}, Abest{0.}, Bbest{0.};
1499 Chi2 = -1.;
1500 for (int j=0;j<=Ntry;j++) { //loop over all combinations E.g. 0010 meas 2-th has -R and the other +R
1501 for (int k=0;k<Nbest;k++) { //loop over best hits
1502 int Isig = (j & (1<<k)) ? 1 : 0;
1503 RRi[bestIdx[k]] = (Isig==1) ? -RI[bestIdx[k]] : RI[bestIdx[k]];
1504 }
1505
1506 float Aj = A0;
1507 float Bj = B0;
1508 Circfit(Nmeas,XI,YI,RRi,WI,Aj,Bj,Chi2,&bestIdx);
1509 Circfit(Nmeas,XI,YI,RI,WI,Aj,Bj,Chi2);
1510 st_A.push_back(Aj); st_B.push_back(Bj); st_chi2.push_back(Chi2);
1511
1512 if (Chi2>=0.0&&Chi2<=CHbest) {
1513 Abest = Aj;
1514 Bbest = Bj;
1515 CHbest = Chi2;
1516 }
1517 }
1518
1519 std::multimap<float, int>chi_map {};
1520 std::vector<float> t_A {};
1521 std::vector<float> t_B {};
1522 std::vector<float> t_chi2 {};
1523
1524 for (size_t ir=0; ir<st_chi2.size(); ir++) chi_map.insert(std::make_pair(st_chi2.at(ir), ir));
1525
1526 for (std::multimap<float, int>::iterator jt = chi_map.begin(); jt != chi_map.end(); ++jt) {
1527 t_A.push_back(st_A.at(jt->second));
1528 t_B.push_back(st_B.at(jt->second));
1529 t_chi2.push_back(st_chi2.at(jt->second));
1530 }
1531
1532 for (int nv=0; nv<6; nv++) {
1533 SlopeCand[nv] = t_A[nv];
1534 InterceptCand[nv] = t_B[nv];
1535 Chi2Cand[nv] = t_chi2[nv];
1536 }
1537
1538 // ... and finally with all the points
1539 A = Abest;
1540 B = Bbest;
1541 Circfit(Nmeas,XI,YI,RI,WI,A,B,Chi2);
1542
1543 if (Chi2>=0.0) {
1544 Pchi2 = TMath::Prob(Chi2, Nmeas - 2);
1545 }
1546
1547 return;
1548}
1549
1550// --------------------------------------------------------------------------------
1551// --------------------------------------------------------------------------------
1552
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
Scalar theta() const
theta method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_DEBUG(x)
static const int B0
Definition AtlasPID.h:122
static Double_t a
static Double_t sc
#define NMEAMX
#define NCAND
struct TBPatternUnitContext S1
const float ZERO_LIMIT
constexpr int pow(int base, int exp) noexcept
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
virtual double phi() const override final
Methods to retrieve data members.
void makePrivateStore()
Create a new (empty) private store for this object.
Gaudi::Property< double > m_rwidth_Endcapmid_second
void Circfit(const int, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, float &, float &, float &, std::vector< int > *idx_vec=nullptr) const
virtual StatusCode initialize() override
ToolHandle< AlphaBetaEstimate > m_alphaBetaEstimate
void Evlfit(TrigL2MuonSA::PBFitResult &fitres) const
Gaudi::Property< double > m_rwidth_Endcapinn_first
StatusCode setMCFlag(bool use_mcLUT)
StatusCode superPointFitter(TrigL2MuonSA::TrackPattern &trackPattern) const
double fromAlphaPtToInn(const TrigL2MuonSA::TgcFitResult &tgcFitResult, const TrigL2MuonSA::TrackPattern &trackPattern) const
ToolHandle< PtFromAlphaBeta > m_ptFromAlphaBeta
StatusCode findSuperPointsSimple(const TrigRoiDescriptor *p_roids, const TrigL2MuonSA::MuonRoad &muonRoad, TrigL2MuonSA::TgcFitResult &tgcFitResult, std::vector< TrigL2MuonSA::TrackPattern > &v_trackPatterns, TrigL2MuonSA::StgcHits &stgcHits, TrigL2MuonSA::MmHits &mmHits) const
void updateInnSP(TrigL2MuonSA::TrackPattern &trackPattern, double &aw, double &tgc_aw, double &bw) const
Gaudi::Property< double > m_rwidth_Endcapmid_first
void findLayerCombination(std::vector< unsigned int > &a, int n, int r, std::vector< std::vector< unsigned int > > &c, int &nr) const
void Circles(const int, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, float &, float &, float &, float &) const
Gaudi::Property< double > m_rwidth_Endcapout_first
StatusCode findSuperPoints(const TrigL2MuonSA::MuonRoad &muonRoad, TrigL2MuonSA::RpcFitResult &rpcFitResult, std::vector< TrigL2MuonSA::TrackPattern > &v_trackPatterns) const
ToolHandle< NswStationFitter > m_nswStationFitter
MuFastStationFitter(const std::string &type, const std::string &name, const IInterface *parent)
Gaudi::Property< double > m_rwidth_Endcapout_second
Gaudi::Property< double > m_endcapmid_mdt_chi2_limit
void findSubLayerCombination(std::vector< unsigned int > &a, int n, int r, std::vector< unsigned int > &b, int index, int num, std::vector< std::vector< unsigned int > > &c, int &nr) const
Gaudi::Property< double > m_mdt_driftspace_downlimit
Gaudi::Property< double > m_mdt_driftspace_uplimit
void stationSPFit(TrigL2MuonSA::MdtHits &mdtSegment, TrigL2MuonSA::SuperPoint &superPoint, TrigL2MuonSA::PBFitResult &pbFitResult, int s_address, int i_station, double aw) const
Gaudi::Property< double > m_rwidth_Endcapinn_second
void Xline(const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, const std::array< float, NMEAMX > &, const int, float &, float &, float &, float &, float &, float &, std::vector< int > *idx_vec=nullptr) const
void makeReferenceLine(TrigL2MuonSA::TrackPattern &trackPattern, const TrigL2MuonSA::MuonRoad &muonRoad) const
ToolHandle< ITrigMuonBackExtrapolator > m_backExtrapolator
double aw[N_STATION][N_SECTOR]
Definition MuonRoad.h:83
double bw[N_STATION][N_SECTOR]
Definition MuonRoad.h:84
std::array< float, NMEAMX > YILIN
std::array< float, NMEAMX > RESI
std::array< float, NMEAMX > RILIN
std::array< float, NMEAMX > XILIN
std::array< float, NCAND > Chi2Cand
std::array< float, NCAND > InterceptCand
std::array< float, NMEAMX > WILIN
std::array< float, NCAND > SlopeCand
TrigL2MuonSA::MdtHits mdtSegments[s_NCHAMBER]
Definition TrackData.h:57
TrigL2MuonSA::SuperPoint superPoints[s_NCHAMBER]
Definition TrackData.h:60
nope - should be used for standalone also, perhaps need to protect the class def bits ifndef XAOD_ANA...
void setEtaMS(float value)
Set the eta at muon spectrometer.
void setRMS(float value)
Set the R at muon spectrometer.
void setSAddress(int value)
Set the station address of the muon.
void setPhiMS(float value)
Set the phi at muon spectrometer.
void setPt(float pt)
Set the transverse momentum ( ) of the muon.
void setZMS(float value)
Set the Z at muon spectrometer.
int ir
counter of the current depth
Definition fastadd.cxx:49
int r
Definition globals.cxx:22
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:146
std::vector< StgcHitData > StgcHits
Definition StgcData.h:49
std::vector< MdtHitData > MdtHits
Definition MdtData.h:56
std::vector< MmHitData > MmHits
Definition MmData.h:47
Definition index.py:1
@ EndcapOuter
Outer station in the endcap spectrometer.
@ EndcapMiddle
Middle station in the endcap spectrometer.
L2StandAloneMuon_v2 L2StandAloneMuon
Define the latest version of the muon SA class.
hold the test vectors and ease the comparison
Definition dumpNPs.cxx:34