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;
447 unsigned int MdtLayerHits[MAX_STATION][MAX_LAYER];
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++) {
452 MdtLayerHits[i_st][i_ly] = 0;
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;
519 MdtLayerHits[i_station][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++) {
546 if (MdtLayerHits[i_station][i_layer]==0)
continue;
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{};
590 if (MdtLayerHits[i_station][i_layer]==0)
continue;
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];
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];