434 TrigL2MuonSA::MdtHits::iterator itMdtHit;
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;
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--) {
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) {
677 int hit_index = std::distance(mdtSegment.begin(),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;
705 if ( std::abs(itMdtHit->DriftSpace) >
ZERO_LIMIT &&
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));
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));
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];
820 for (
int i=0;i<pbFitResult.
NPOI;i++) superPoint.
Residual[i] = pbFitResult.
RESI[i];
822 for (
int cand=0; cand<6; cand++) {
826 superPoint.
Chi2Cand[cand] = t_Chi2[cand];
841 for (
int cand=0; cand<6; cand++) {
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];
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];
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);
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{};
1490 for (
int i=0; i<
NCAND; i++) {
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));
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));
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);