432 {
433
434 TrigL2MuonSA::MdtHits::iterator itMdtHit;
435
436 unsigned int i_layer_max = 0;
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
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;
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
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
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) {
676
677 int hit_index = std::distance(mdtSegment.begin(),itMdtHit);
678
679 if(mdtSegment.at(hit_index).isOutlier>1) continue;
680
682
684
685 for (unsigned int j=0; j<hitarray.size(); j++) {
686
687 if (hitarray[j]==hit_index) {
689 break;
690 }
691 }
692
693 if (fd==0) continue;
694
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
711 pbFitResult.
RILIN[
count] = itMdtHit->DriftSpace;
714
717
718 sum_Z_used = sum_Z_used + itMdtHit->Z;
719 sum_R_used = sum_R_used + itMdtHit->R;
720 } else {
722 }
723 }
724
726 for (
int i=0;
i<pbFitResult.
NPOI;
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) {
736
737
738 for (int cand=0; cand<6; cand++) {
739
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);
749 }
750 }
751 }
752 }
753 }
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){
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];
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) {
806
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];
819
821
822 for (int cand=0; cand<6; cand++) {
826 superPoint.
Chi2Cand[cand] = t_Chi2[cand];
827 }
828 }
829 }
830
831 Chbest=t_Chi2[0];
832
833 if (real_layer>3) {
835
840
841 for (int cand=0; cand<6; cand++) {
845 }
846 return;
847 }
848 }
849
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;
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];
875 superPoint.
Chi2Cand[cand] = Maxlayers_Chi2[cand];
876 }
877 }
878 }
879 }
880 }
881
882 return;
883}
constexpr int pow(int base, int exp) noexcept
Gaudi::Property< double > m_rwidth_Endcapmid_second
void findLayerCombination(std::vector< unsigned int > &a, int n, int r, std::vector< std::vector< unsigned int > > &c, int &nr) const
Gaudi::Property< double > m_rwidth_Endcapout_second
Gaudi::Property< double > m_endcapmid_mdt_chi2_limit
Gaudi::Property< double > m_rwidth_Endcapinn_second
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
int count(std::string s, const std::string ®x)
count how many occurances of a regx are in a string
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
1.14.0