24 const std::string&
name,
36 ATH_CHECK( m_backExtrapolator.retrieve() );
38 ATH_CHECK( m_alphaBetaEstimate.retrieve() );
41 ATH_CHECK( m_ptFromAlphaBeta.retrieve() );
44 ATH_CHECK( m_nswStationFitter.retrieve(DisableTool{m_nswStationFitter.empty()}) );
46 return StatusCode::SUCCESS;
53 m_use_mcLUT = use_mcLUT;
57 if ( ptEndcapLUTSvc.retrieve().isFailure() ) {
59 return StatusCode::FAILURE;
62 m_alphaBetaEstimate->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
64 m_ptFromAlphaBeta->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
67 if ( ptEndcapLUTSvc.retrieve().isFailure() ) {
69 return StatusCode::FAILURE;
72 m_alphaBetaEstimate->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
74 m_ptFromAlphaBeta->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
77 ATH_MSG_DEBUG(
"Completed tp set " << (m_use_mcLUT?
"MC":
"not MC") <<
" flag" );
79 return StatusCode::SUCCESS;
87 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns)
const
99 itTrack.isRpcFailure = !rpcFitResult.
isSuccess;
105 return StatusCode::SUCCESS;
114 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns,
122 itTrack.phiMSDir = tgcFitResult.
phiDir;
129 itTrack.isTgcFailure =
true;
134 if(!m_nswStationFitter.empty())
135 ATH_CHECK( m_nswStationFitter->superPointFitter(p_roids, itTrack, stgcHits, mmHits) );
138 return StatusCode::SUCCESS;
147 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns,
155 itTrack.phiMSDir = tgcFitResult.
phiDir;
158 itTrack.isTgcFailure =
true;
161 ATH_CHECK( superPointFitter(itTrack, muonRoad) );
163 makeReferenceLine(itTrack, muonRoad);
164 ATH_CHECK( m_alphaBetaEstimate->setAlphaBeta(p_roids, tgcFitResult, itTrack, muonRoad) );
166 if ( itTrack.etaBin < -1 ) {
168 if(itTrack.etaBin <= -1) itTrack.etaBin = 0;
171 ATH_CHECK( m_ptFromAlphaBeta->setPt(itTrack,tgcFitResult) );
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);
178 if(!m_nswStationFitter.empty())
179 ATH_CHECK( m_nswStationFitter->superPointFitter(p_roids, itTrack, stgcHits, mmHits) );
183 return StatusCode::SUCCESS;
192 const unsigned int MAX_STATION = 10;
193 const float SIGMA = 0.0080;
194 const float DRIFTSPACE_LIMIT = 16.;
195 const int MIN_MDT_FOR_FIT = 3;
201 if (mdtSegment.size()==0)
continue;
204 float sigma {0.}, phim {0.}, Xor {0.}, Yor {0.}, Ymid {0.};
210 if (itMdtHit.isOutlier)
continue;
213 Ymid = itMdtHit.cYmid;
220 phim = itMdtHit.cPhip;
221 sigma = (std::abs(itMdtHit.DriftSigma) >
ZERO_LIMIT)? itMdtHit.DriftSigma: SIGMA;
224 std::abs(itMdtHit.DriftSpace) < DRIFTSPACE_LIMIT &&
229 pbFitResult.
RILIN[
count] = itMdtHit.DriftSpace;
240 for(
int i=0;
i<pbFitResult.
NPOI;
i++) {
242 <<
i <<
"/" << pbFitResult.
XILIN[
i] <<
"/" << pbFitResult.
YILIN[
i]
243 <<
"/" << pbFitResult.
RILIN[
i] <<
"/" << pbFitResult.
WILIN[
i]);
249 if (
count >= MIN_MDT_FOR_FIT) {
253 superPoint.Npoint = pbFitResult.
NPOI;
260 superPoint.R = (Ymid - Yor - pbFitResult.
BLIN)/pbFitResult.
ALIN + Xor;
261 superPoint.Alin = 1./pbFitResult.
ALIN;
262 superPoint.Blin = -pbFitResult.
BLIN/pbFitResult.
ALIN;
267 superPoint.Z = pbFitResult.
ALIN*(Ymid - Xor) + pbFitResult.
BLIN + Yor;
268 superPoint.Alin = pbFitResult.
ALIN;
269 superPoint.Blin = pbFitResult.
BLIN;
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];
280 ATH_MSG_DEBUG(
"... Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
282 << superPoint.R <<
"/" << superPoint.Z <<
"/" << superPoint.Alin <<
"/"
283 << superPoint.Blin <<
"/" << superPoint.Phim <<
"/" << superPoint.Xor <<
"/"
284 << superPoint.Yor <<
"/" << superPoint.Chi2 <<
"/" << superPoint.PChi2);
287 return StatusCode::SUCCESS;
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;
306 if (mdtSegment.size()==0)
continue;
313 float sigma {0.}, phim {0.}, Xor {0.}, Yor {0.}, Ymid {0.};
318 if (itMdtHit.isOutlier)
continue;
322 Ymid = itMdtHit.cYmid;
329 phim = itMdtHit.cPhip;
330 sigma = (std::abs(itMdtHit.DriftSigma) >
ZERO_LIMIT)? itMdtHit.DriftSigma: SIGMA;
333 std::abs(itMdtHit.DriftSpace) < DRIFTSPACE_LIMIT &&
336 pbFitResult.XILIN[
count] = itMdtHit.R - Xor;
337 pbFitResult.YILIN[
count] = itMdtHit.Z - Yor;
338 pbFitResult.RILIN[
count] = itMdtHit.DriftSpace;
340 pbFitResult.RESI[
count] = 0.;
342 pbFitResult.NPOI =
count;
348 ATH_MSG_DEBUG(
"... MDT hit used in fit #=" << pbFitResult.NPOI);
349 for(
int i=0;
i<pbFitResult.NPOI;
i++) {
351 <<
i <<
"/" << pbFitResult.XILIN[
i] <<
"/" << pbFitResult.YILIN[
i]
352 <<
"/" << pbFitResult.RILIN[
i] <<
"/" << pbFitResult.WILIN[
i]);
354 if (
count >= MIN_MDT_FOR_FIT) {
355 Evlfit( pbFitResult);
357 float bc = (Ymid - Xor);
358 float X = (pbFitResult.ALIN*bc)+pbFitResult.BLIN ;
360 superPoint.Npoint = pbFitResult.NPOI;
362 if (std::abs(pbFitResult.ALIN) >
ZERO_LIMIT) {
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;
368 superPoint.R = bc + Xor;
369 superPoint.Z =
X + Yor;
370 superPoint.Alin = pbFitResult.ALIN;
371 superPoint.Blin = pbFitResult.BLIN;
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];
383 ATH_MSG_DEBUG(
"...Special Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
385 << superPoint.R <<
"/" << superPoint.Z <<
"/" << superPoint.Alin <<
"/"
386 << superPoint.Blin <<
"/" << superPoint.Phim <<
"/" << superPoint.Xor <<
"/"
387 << superPoint.Yor <<
"/" << superPoint.Chi2 <<
"/" << superPoint.PChi2);
394 unsigned int sumN = 0;
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;}
401 if (std::abs(itMdtHit.DriftSpace) < m_mdt_driftspace_downlimit ||
402 std::abs(itMdtHit.DriftSpace) > m_mdt_driftspace_uplimit){
403 itMdtHit.isOutlier = 2;
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;
416 itMdtHit.isOutlier = 2;
420 if(sumN==0)
continue;
422 stationSPFit(mdtSegment, superPoint,pbFitResult, trackPattern.
s_address,
chamber,aw);
426 return StatusCode::SUCCESS;
436 unsigned int i_layer_max = 0;
438 float Xor, Yor,
sigma, phim=0;
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;
445 const unsigned int MAX_LAYER = 12;
448 std:: vector<unsigned int> MdtLayerHits_index[MAX_STATION][MAX_LAYER];
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++) {
453 MdtLayerHits_index[i_st][i_ly].clear();
457 float Chbest = 1.e25;
462 unsigned int sumN[8];
465 for (
unsigned int i_st=0; i_st<8; i_st++) {
473 for (
unsigned int i_hit=0; i_hit<mdtSegment.size(); i_hit++){
475 if (mdtSegment.at(i_hit).isOutlier>1)
continue;
477 double Z = mdtSegment.at(i_hit).Z;
478 double R = mdtSegment.at(i_hit).R;
480 sumZ[i_station] = sumZ[i_station] + Z;
481 sumR[i_station] = sumR[i_station] + R;
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];
490 if (sumN[i_station]==0)
return;
492 for (
unsigned int i_hit=0; i_hit<mdtSegment.size(); i_hit++) {
494 if (mdtSegment.at(i_hit).isOutlier>1)
continue;
496 double Z = mdtSegment.at(i_hit).Z;
497 double R = mdtSegment.at(i_hit).R;
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;
503 double nbw = aw*Z+(avR[i_station]-aw*avZ[i_station]);
505 if ( R>(nbw-nsWidth) && R<(nbw+nsWidth) ) {
506 mdtSegment.at(i_hit).isOutlier = 0;
508 mdtSegment.at(i_hit).isOutlier = 2;
513 for (
unsigned int i_hit=0; i_hit<mdtSegment.size(); i_hit++) {
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;
520 MdtLayerHits_index[i_station][i_layer].push_back(i_hit);
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;
544 for (
unsigned int i_layer=0; i_layer<=i_layer_max; i_layer++) {
548 Ly_1st.push_back(i_layer);
551 const int real_layer= Ly_1st.size();
552 std::vector<std::vector<unsigned int> > Ly_flg{};
554 for (
int pr=real_layer;
pr>=3;
pr--) {
567 findLayerCombination(Ly_1st, real_layer,
pr,Ly_flg, total_cp);
569 for (
unsigned int i=0;
i<Ly_flg.size();
i++) {
571 std::vector<std::vector<int> >tID{};
572 std::vector<std::vector<int> >tIndex{};
577 for (
int ti=0; ti<8; ti++) {
578 for (
int tj=0; tj<2; tj++) {
580 tubeindex[ti][tj] = 0;
584 for (
unsigned int j=0; j<Ly_flg[
i].size(); j++) {
586 int i_layer = Ly_flg[
i][j];
587 std::vector<int> tid{};
588 std::vector<int> tindex{};
592 float tube_1st = 999999.;
593 float tube_2nd = 999999.;
595 int layer_2nd = 9999;
597 for (
unsigned int i_hit=0; i_hit<
MdtLayerHits[i_station][i_layer]; i_hit++) {
599 unsigned int i_index = MdtLayerHits_index[i_station][i_layer].at(i_hit);
601 if (mdtSegment.at(i_index).isOutlier>1)
continue;
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);
606 if (dis_tube<tube_1st) {
608 layer_2nd = layer_1st;
611 }
else if (dis_tube<tube_2nd) {
617 if ( layer_1st != 9999 ) {
618 mdtSegment.at(layer_1st).isOutlier = 0;
620 tindex.push_back(layer_1st);
623 if ( layer_2nd != 9999 ) {
624 mdtSegment.at(layer_2nd).isOutlier = 1;
626 tindex.push_back(layer_2nd);
630 tIndex.push_back(tindex);
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];
640 std::vector<int> isg;
641 std::vector<int> hitarray;
643 int ntry = (
int)floor(
pow(2.,
pr))-1;
645 for (
int ntryi=0; ntryi<=ntry; ntryi++) {
651 for (
int ntryj=1; ntryj<=
pr; ntryj++) {
653 int Isg = (ntryi&(
int)
pow(2.,ntryj-1))? 1 : 0;
655 if (tube_ID[ntryj-1][Isg] != 0) yhit = 1;
656 sumid = sumid * yhit;
660 for (
unsigned int tt=0;
tt<isg.size();
tt++) {
661 int tindex = tubeindex[
tt][isg[
tt]];
662 hitarray.push_back(tindex);
670 float sum_Z_used = 0.;
671 float sum_R_used = 0.;
673 if (hitarray.size()==0)
continue;
675 for (itMdtHit=mdtSegment.begin(); itMdtHit!=mdtSegment.end(); ++itMdtHit) {
679 if(mdtSegment.at(hit_index).isOutlier>1)
continue;
685 for (
unsigned int j=0; j<hitarray.size(); j++) {
687 if (hitarray[j]==hit_index) {
702 phim = itMdtHit->cPhip;
703 sigma = (std::abs(itMdtHit->DriftSigma) >
ZERO_LIMIT)? itMdtHit->DriftSigma: SIGMA;
706 std::abs(itMdtHit->DriftSpace) < DRIFTSPACE_LIMIT &&
707 std::abs(itMdtHit->DriftTime) >
ZERO_LIMIT ) {
711 pbFitResult.
RILIN[
count] = itMdtHit->DriftSpace;
718 sum_Z_used = sum_Z_used + itMdtHit->Z;
719 sum_R_used = sum_R_used + itMdtHit->R;
726 for (
int i=0;
i<pbFitResult.
NPOI;
i++) {
728 <<
i <<
"/" << pbFitResult.
XILIN[
i] <<
"/" << pbFitResult.
YILIN[
i]
729 <<
"/" << pbFitResult.
RILIN[
i] <<
"/" << pbFitResult.
WILIN[
i]);
732 if (
count >= MIN_MDT_FOR_FIT) {
738 for (
int cand=0; cand<6; cand++) {
741 Line_A.push_back(1/pbFitResult.
SlopeCand[cand]);
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);
755 if (Line_Chi2.size()==0)
continue;
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;
776 for (
unsigned int ir=0;
ir<Line_Chi2.size();
ir++) chi_map.insert(std::make_pair(Line_Chi2.at(
ir),
ir));
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));
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));
794 superPoint.
Npoint = t_count[0];
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];
805 if (s_address == -1) {
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];
814 superPoint.
Phim = t_phim[0];
815 superPoint.
Xor = t_Xor[0];
816 superPoint.
Yor = t_Yor[0];
817 superPoint.
Chi2 = t_Chi2[0];
822 for (
int cand=0; cand<6; cand++) {
826 superPoint.
Chi2Cand[cand] = t_Chi2[cand];
834 if ((i_station == 3 || i_station == 5) &&
pr==4 && Chbest > m_endcapmid_mdt_chi2_limit) {
841 for (
int cand=0; cand<6; cand++) {
850 if (Chbest<m_endcapmid_mdt_chi2_limit){
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);
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];
875 superPoint.
Chi2Cand[cand] = Maxlayers_Chi2[cand];
890 const unsigned int MAX_STATION = 8;
896 float Chi2_cand[8][6];
900 const int middle = 4;
903 for (
unsigned int i_station=4; i_station<MAX_STATION; i_station++) {
910 spChi2[i_station]=0.;
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.;
918 if (i_station<4 || i_station>5)
continue;
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;
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];
934 float test_diff = 1.e25;
935 float best_diff = 1.e25;
937 float match_midA = 0.;
938 float match_outA = 0.;
957 for (
int m=0;
m<6;
m++) {
959 test_diff = std::abs(A_cand[middle][
m] - aw[middle]);
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];
976 if(std::abs(spZ[outer]-spZ[middle]) >
ZERO_LIMIT) sp_line = (spR[outer]-spR[middle])/(spZ[outer]-spZ[middle]);
978 for (
int t=0;
t<2; ++
t) {
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;
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];
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];
1006 spA[middle] = match_midA;
1007 }
else if (std::abs(rmatched) >
ZERO_LIMIT) {
1008 spA[middle] = rmatched;
1011 if (std::abs(match_outA) >
ZERO_LIMIT) spA[outer] = match_outA;
1014 for (
unsigned int i_station=4; i_station<MAX_STATION; i_station++) {
1015 if (i_station<4 || i_station>5)
continue;
1018 superPoint.
Alin =spA[i_station];
1019 superPoint.Blin =spB[i_station];
1020 superPoint.Chi2 =spChi2[i_station];
1031 float MiddleSlope = 0;
1032 float OuterSlope = 0;
1034 for (
int i_station=4; i_station<6; i_station++) {
1043 if ( superPoint.Npoint > 2 && superPoint.R > 0.) {
1044 if ( i_station==4 ) {
1045 MiddleSlope = superPoint.
Alin;
1047 if ( i_station==5 ) {
1048 OuterSlope = superPoint.Alin;
1053 double mdtpT = std::abs(tgcFitResult.
tgcPT);
1058 }
else if (std::abs(tgcFitResult.
tgcPT)>=8.0 && std::abs(MiddleSlope) >
ZERO_LIMIT) {
1062 mdtpT = (std::abs(tgcFitResult.
tgcPT)>1
e-5)? mdtpT*(tgcFitResult.
tgcPT/std::abs(tgcFitResult.
tgcPT)) : 0;
1067 double extrInnerEta = 0;
1072 muonSA->
setPt(mdtpT);
1080 StatusCode sc = m_backExtrapolator->give_eta_phi_at_vertex(muonSA, eta,sigma_eta,phi,sigma_phi,0.);
1082 if (
sc.isSuccess() ){
1085 extrInnerEta = etaMiddle;
1092 naw =
std::tan(theta)*(std::abs(extrInnerEta)/extrInnerEta);
1105 double nrWidth = m_rwidth_Endcapinn_first;
1106 unsigned int sumN[8];
1108 for (
unsigned int i_st=0; i_st<8;i_st++) {
1113 const int i_station = 3;
1114 int chamberID = i_station;
1119 if (mdtSegment.size()==0)
return;
1122 if (std::abs(itMdtHit.DriftSpace) < m_mdt_driftspace_downlimit ||
1123 std::abs(itMdtHit.DriftSpace) > m_mdt_driftspace_uplimit){
1124 itMdtHit.isOutlier = 2;
1128 if (itMdtHit.isOutlier > 1)
continue;
1130 double Z = itMdtHit.Z;
1131 double R = itMdtHit.R;
1132 double nbw = tgc_aw*Z + bw;
1134 if (R>(nbw-nrWidth) && R<(nbw+nrWidth)) {
1135 itMdtHit.isOutlier = 0;
1138 itMdtHit.isOutlier = 2;
1143 if (sumN[i_station]==0)
return;
1145 stationSPFit(mdtSegment, superPoint,pbFitResult, trackPattern.
s_address,i_station, aw);
1149 for (
int cand=0; cand<
NCAND; cand++) {
1150 float ds=std::abs(superPoint.SlopeCand[cand]-aw);
1153 superPoint.Alin = superPoint.SlopeCand[cand];
1154 superPoint.Blin = superPoint.InterceptCand[cand];
1155 superPoint.Chi2 = superPoint.Chi2Cand[cand];
1166 std::vector<std::vector<unsigned int> > &
c,
1169 std::vector<unsigned int>
b(
r,0);
1171 findSubLayerCombination(
a,
n,
r,
b,0,
r,
c,nr);
1182 std::vector<unsigned int> &
b,
1185 std::vector<std::vector<unsigned int> > &
c,
1191 std::vector<unsigned int>
t;
1196 for (
int j = 0;j<
r; j++)
t.push_back(
b[j]);
1203 findSubLayerCombination(
a,
n,
r,
b,
i+1,
num-1,
c,nr);
1229 float Xnor = 1. / std::sqrt(1. + pbFitResult.
ALIN * pbFitResult.
ALIN);
1231 for(
int j=0;j<pbFitResult.
NPOI;j++) {
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;
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
1275 if (Nmeas<=2||Nmeas>=
NMEAMX+1)
return;
1279 Xline(XI,YI,WI,Nmeas,
A,B,SAA,SBB,SAB,Chi2);
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) {
1287 std::vector<int> bestIdx{};
1288 for (
int j=0;j<Nmeas;j++) {
1289 if (WI[j]>=WIlim || Ngood<=3) {
1290 if (bestIdx.size() < 4){
1291 bestIdx.push_back(j);
1294 bestIdx.at(bestIdx.size()-2) = bestIdx.back();
1299 const int Nbest = bestIdx.size();
1301 const int Ntry = (1 << Nbest) - 1;
1302 std::array<float,NMEAMX> RRi{};
1303 float Abest{0.}, Bbest {0.}, CHbest{1.e25};
1304 for (
int j=0;j<=Ntry;j++) {
1305 for (
int k=0;
k<Nbest;
k++) {
1306 const int Isig = (j & (1<<
k)) ? 1 : 0;
1307 RRi[bestIdx[
k]] = (Isig==1) ? -RI[bestIdx[
k]] : RI[bestIdx[
k]];
1312 Circfit(Nmeas,XI,YI,RRi,WI,Aj,Bj,chi2j,&bestIdx);
1314 if (chi2j>=0.0&&chi2j<=CHbest) {
1325 Circfit(Nmeas,XI,YI,RI,WI,
A,B,Chi2);
1328 Pchi2 = TMath::Prob(Chi2, Nmeas - 2);
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
1351 const bool use_all {!idx_vec};
1352 std::vector<int> temp{};
1355 std::iota(temp.begin(), temp.end(), 0);
1359 std::array<float,NMEAMX> XX{},YY{};
1361 float Test, Toll {.1};
1366 float Xnor = 1. / std::sqrt(1. +
A *
A);
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]);
1378 for(
const int& j : *idx_vec){
1379 XX[j] = XI[j] - Xnor * RI[j] *
A;
1380 YY[j] = YI[j] + Xnor * RI[j];
1384 Xline(XX,YY,WI,Nmeas,
A,B,SAA,SBB,SAB,Chi2, idx_vec);
1386 Test = ((Aold-
A)*(Aold-
A))/ SAA + ((Bold-B)*(Bold-B))/ SBB;
1388 }
while(
Test>=Toll&&Niter<=20);
1401 const std::array<float, NMEAMX>&
Y,
1402 const std::array<float, NMEAMX>&
W,
1404 float&
A,
float& B,
float& SAA,
float& SBB,
float& SAB,
float& Square, std::vector<int>* idx_vec)
const
1406 std::vector<int> temp{};
1409 std::iota(temp.begin(), temp.end(), 0);
1413 float S1{0.},SX{0.},SY{0.},SXX{0.},SXY{0.},SYY{0.};
1415 for(
const int& j : *idx_vec) {
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];
1424 float Deter =
S1 * SXX - SX * SX;
1427 A = (
S1 * SXY - SX * SY) / Deter;
1428 B = (SY * SXX - SX * SXY) / Deter;
1434 A= (
S1 * SXY - SX * SY > 0.) ? 9.e+5 : -9.e+5;
1435 B = SY/
S1 - SX/
S1 *
A;
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;
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
1460 if (Nmeas<=2||Nmeas>=
NMEAMX+1)
return;
1463 float A0,B0, SAA,SBB,SAB,Square;
1464 Xline(XI,YI,WI,Nmeas,A0,B0,SAA,SBB,SAB,Square);
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) {
1472 std::vector<int> bestIdx{};
1473 for (
int j=0;j<Nmeas;j++) {
1474 if (WI[j]>=WIlim || Ngood<=3) {
1475 if (bestIdx.size() < 4){
1476 bestIdx.push_back(j);
1479 bestIdx.at(bestIdx.size()-2) = bestIdx.back();
1484 const int Nbest = bestIdx.size();
1486 std::vector<float> st_chi2{};
1487 std::vector<float> st_A{};
1488 std::vector<float> st_B{};
1492 InterceptCand[
i] = 0.;
1496 const int Ntry = (1 << Nbest) - 1;
1497 std::array<float,NMEAMX> RRi{};
1498 float CHbest{1.e25}, Abest{0.}, Bbest{0.};
1500 for (
int j=0;j<=Ntry;j++) {
1501 for (
int k=0;
k<Nbest;
k++) {
1502 int Isig = (j & (1<<
k)) ? 1 : 0;
1503 RRi[bestIdx[
k]] = (Isig==1) ? -RI[bestIdx[
k]] : RI[bestIdx[
k]];
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);
1512 if (Chi2>=0.0&&Chi2<=CHbest) {
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 {};
1524 for (
size_t ir=0;
ir<st_chi2.size();
ir++) chi_map.insert(std::make_pair(st_chi2.at(
ir),
ir));
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));
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];
1541 Circfit(Nmeas,XI,YI,RI,WI,
A,B,Chi2);
1544 Pchi2 = TMath::Prob(Chi2, Nmeas - 2);