26 const std::string&
name,
38 ATH_CHECK( m_backExtrapolator.retrieve() );
40 ATH_CHECK( m_alphaBetaEstimate.retrieve() );
43 ATH_CHECK( m_ptFromAlphaBeta.retrieve() );
46 ATH_CHECK( m_nswStationFitter.retrieve(DisableTool{m_nswStationFitter.empty()}) );
48 return StatusCode::SUCCESS;
55 m_use_mcLUT = use_mcLUT;
59 if ( ptEndcapLUTSvc.retrieve().isFailure() ) {
61 return StatusCode::FAILURE;
64 m_alphaBetaEstimate->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
66 m_ptFromAlphaBeta->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
69 if ( ptEndcapLUTSvc.retrieve().isFailure() ) {
71 return StatusCode::FAILURE;
74 m_alphaBetaEstimate->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
76 m_ptFromAlphaBeta->setMCFlag(m_use_mcLUT, &*ptEndcapLUTSvc);
79 ATH_MSG_DEBUG(
"Completed tp set " << (m_use_mcLUT?
"MC":
"not MC") <<
" flag" );
81 return StatusCode::SUCCESS;
90 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns)
const
105 itTrack.isRpcFailure =
true;
112 return StatusCode::SUCCESS;
121 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns,
129 itTrack.phiMSDir = tgcFitResult.
phiDir;
136 itTrack.isTgcFailure =
true;
141 if(!m_nswStationFitter.empty())
142 ATH_CHECK( m_nswStationFitter->superPointFitter(p_roids, itTrack, stgcHits, mmHits) );
145 return StatusCode::SUCCESS;
154 std::vector<TrigL2MuonSA::TrackPattern>& v_trackPatterns,
162 itTrack.phiMSDir = tgcFitResult.
phiDir;
165 itTrack.isTgcFailure =
true;
168 ATH_CHECK( superPointFitter(itTrack, muonRoad) );
170 makeReferenceLine(itTrack, muonRoad);
171 ATH_CHECK( m_alphaBetaEstimate->setAlphaBeta(p_roids, tgcFitResult, itTrack, muonRoad) );
173 if ( itTrack.etaBin < -1 ) {
175 if(itTrack.etaBin <= -1) itTrack.etaBin = 0;
178 ATH_CHECK( m_ptFromAlphaBeta->setPt(itTrack,tgcFitResult) );
180 double exInnerA = fromAlphaPtToInn(tgcFitResult,itTrack);
181 double bw = muonRoad.
bw[3][0];
182 double aw = muonRoad.
aw[3][0];
183 if(std::abs(exInnerA) >
ZERO_LIMIT) updateInnSP(itTrack, exInnerA, aw,bw);
185 if(!m_nswStationFitter.empty())
186 ATH_CHECK( m_nswStationFitter->superPointFitter(p_roids, itTrack, stgcHits, mmHits) );
190 return StatusCode::SUCCESS;
201 float Xor, Yor,
sigma, rm=0., phim=0;
202 float Ymid, Xmid, Amid;
204 const unsigned int MAX_STATION = 10;
205 const float SIGMA = 0.0080;
206 const float DRIFTSPACE_LIMIT = 16.;
207 const int MIN_MDT_FOR_FIT = 3;
226 if (mdtSegment->size()==0)
continue;
231 if (itMdtHit.isOutlier)
continue;
237 Amid = itMdtHit.cAmid;
238 Xmid = itMdtHit.cXmid;
239 Ymid = itMdtHit.cYmid;
246 phim = itMdtHit.cPhip;
247 sigma = (std::abs(itMdtHit.DriftSigma) >
ZERO_LIMIT)? itMdtHit.DriftSigma: SIGMA;
253 if (
chamber == 6 ) station = 3;
254 if ( std::abs(itMdtHit.DriftSpace) >
ZERO_LIMIT &&
255 std::abs(itMdtHit.DriftSpace) < DRIFTSPACE_LIMIT &&
262 itMdtHit.DriftSpace: SetDriftSpace(itMdtHit.DriftTime, itMdtHit.R, itMdtHit.Z, phim, trackPattern.
phiMSDir);
265 pbFitResult.
IDMEA[
count] = station*10 + itMdtHit.Layer;
280 for(
int i=0;
i<pbFitResult.
NPOI;
i++) {
282 <<
i <<
"/" << pbFitResult.
XILIN[
i] <<
"/" << pbFitResult.
YILIN[
i]
283 <<
"/" << pbFitResult.
RILIN[
i] <<
"/" << pbFitResult.
WILIN[
i]);
286 if (
count >= MIN_MDT_FOR_FIT) {
288 FitFlag = Evlfit(1, pbFitResult);
292 float bc = (Ymid - Xor) -ac*(Xmid - Yor);
293 float X = ( (pbFitResult.
ALIN*bc)+pbFitResult.
BLIN )/( 1-ac*pbFitResult.
ALIN );
301 superPoint->
R = (rm-Yor)/pbFitResult.
ALIN - pbFitResult.
BLIN/pbFitResult.
ALIN + Xor;
302 superPoint->
Alin = 1./pbFitResult.
ALIN;
306 superPoint->
R = ac*
X + bc + Xor;
307 superPoint->
Z =
X + Yor;
308 superPoint->
Alin = pbFitResult.
ALIN;
309 superPoint->
Blin = pbFitResult.
BLIN;
315 superPoint->
R = ac*
X + bc + Xor;
316 superPoint->
Z =
X + Yor;
317 superPoint->
Alin = pbFitResult.
ALIN;
318 superPoint->
Blin = pbFitResult.
BLIN;
321 superPoint->
R = (rm-Yor)/pbFitResult.
ALIN - pbFitResult.
BLIN/pbFitResult.
ALIN + Xor;
322 superPoint->
Alin = 1./pbFitResult.
ALIN;
327 superPoint->
Phim = phim;
328 superPoint->
Xor = Xor;
329 superPoint->
Yor = Yor;
330 superPoint->
Chi2 = pbFitResult.
CHI2;
336 ATH_MSG_DEBUG(
"... Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
338 << superPoint->
R <<
"/" << superPoint->
Z <<
"/" << superPoint->
Alin <<
"/"
339 << superPoint->
Blin <<
"/" << superPoint->
Phim <<
"/" << superPoint->
Xor <<
"/"
340 << superPoint->
Yor <<
"/" << superPoint->
Chi2 <<
"/" << superPoint->
PChi2);
345 return StatusCode::SUCCESS;
353 const unsigned int MAX_STATION = 10;
365 if (mdtSegment->size()==0)
continue;
378 const float SIGMA = 0.0080;
379 const float DRIFTSPACE_LIMIT = 16.;
380 const int MIN_MDT_FOR_FIT = 3;
385 if (itMdtHit.isOutlier)
continue;
390 Amid = itMdtHit.cAmid;
391 Xmid = itMdtHit.cXmid;
392 Ymid = itMdtHit.cYmid;
399 phim = itMdtHit.cPhip;
400 sigma = (std::abs(itMdtHit.DriftSigma) >
ZERO_LIMIT)? itMdtHit.DriftSigma: SIGMA;
404 if (
chamber == 0 ) station = 0;
405 if ( std::abs(itMdtHit.DriftSpace) >
ZERO_LIMIT &&
406 std::abs(itMdtHit.DriftSpace) < DRIFTSPACE_LIMIT &&
413 itMdtHit.DriftSpace: SetDriftSpace(itMdtHit.DriftTime, itMdtHit.R, itMdtHit.Z, phim, trackPattern.
phiMSDir);
416 pbFitResult.
IDMEA[
count] = station*10 + itMdtHit.Layer;
427 for(
int i=0;
i<pbFitResult.
NPOI;
i++) {
429 <<
i <<
"/" << pbFitResult.
XILIN[
i] <<
"/" << pbFitResult.
YILIN[
i]
430 <<
"/" << pbFitResult.
RILIN[
i] <<
"/" << pbFitResult.
WILIN[
i]);
432 if (
count >= MIN_MDT_FOR_FIT) {
433 FitFlag = Evlfit(1, pbFitResult);
437 float bc = (Ymid - Xor) -ac*(Xmid - Yor);
438 float X = ( (pbFitResult.
ALIN*bc)+pbFitResult.
BLIN )/( 1-ac*pbFitResult.
ALIN );
444 superPoint->
R = (rm-Yor)/pbFitResult.
ALIN - pbFitResult.
BLIN/pbFitResult.
ALIN + Xor;
445 superPoint->
Alin = 1./pbFitResult.
ALIN;
448 superPoint->
R = ac*
X + bc + Xor;
449 superPoint->
Z =
X + Yor;
450 superPoint->
Alin = pbFitResult.
ALIN;
451 superPoint->
Blin = pbFitResult.
BLIN;
455 superPoint->
Phim = phim;
456 superPoint->
Xor = Xor;
457 superPoint->
Yor = Yor;
458 superPoint->
Chi2 = pbFitResult.
CHI2;
463 ATH_MSG_DEBUG(
"...Special Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
465 << superPoint->
R <<
"/" << superPoint->
Z <<
"/" << superPoint->
Alin <<
"/"
466 << superPoint->
Blin <<
"/" << superPoint->
Phim <<
"/" << superPoint->
Xor <<
"/"
467 << superPoint->
Yor <<
"/" << superPoint->
Chi2 <<
"/" << superPoint->
PChi2);
474 unsigned int sumN = 0;
476 if(
chamber==3) {nrWidth = m_rwidth_Endcapinn_first;}
477 if(
chamber==4) {nrWidth = m_rwidth_Endcapmid_first;}
478 if(
chamber==5) {nrWidth = m_rwidth_Endcapout_first;}
481 if (std::abs(itMdtHit.DriftSpace) < m_mdt_driftspace_downlimit ||
482 std::abs(itMdtHit.DriftSpace) > m_mdt_driftspace_uplimit){
483 itMdtHit.isOutlier = 2;
487 if(itMdtHit.isOutlier > 1)
continue;
488 double Z = itMdtHit.Z;
489 double R = itMdtHit.R;
490 double nbw = aw*Z + bw;
491 if (R>(nbw-nrWidth) && R<(nbw+nrWidth)){
492 itMdtHit.isOutlier = 0;
495 itMdtHit.isOutlier = 2;
499 if(sumN==0)
continue;
505 return StatusCode::SUCCESS;
515 unsigned int i_layer_max = 0;
518 float Xor, Yor,
sigma, rm=0., phim=0;
519 float Ymid, Xmid, Amid;
521 const unsigned int MAX_STATION = 8;
522 const float SIGMA = 0.0080;
523 const float DRIFTSPACE_LIMIT = 16.;
524 const int MIN_MDT_FOR_FIT = 3;
526 const unsigned int MAX_LAYER = 12;
529 std:: vector<unsigned int> MdtLayerHits_index[MAX_STATION][MAX_LAYER];
531 for (
unsigned int i_st=0; i_st<MAX_STATION; i_st++) {
532 for (
unsigned int i_ly=0; i_ly<MAX_LAYER; i_ly++) {
534 MdtLayerHits_index[i_st][i_ly].clear();
538 float Chbest = 1.e25;
543 unsigned int sumN[8];
546 for (
unsigned int i_st=0; i_st<8; i_st++) {
554 for (
unsigned int i_hit=0; i_hit<mdtSegment->size(); i_hit++){
557 if (mdtSegment->at(i_hit).isOutlier>1)
continue;
559 double Z = mdtSegment->at(i_hit).Z;
560 double R = mdtSegment->at(i_hit).R;
562 sumZ[i_station] = sumZ[i_station] + Z;
563 sumR[i_station] = sumR[i_station] + R;
566 if (sumN[i_station]!=0) {
567 avZ[i_station] = sumZ[i_station]/sumN[i_station];
568 avR[i_station] = sumR[i_station]/sumN[i_station];
572 if (sumN[i_station]==0)
return;
574 for (
unsigned int i_hit=0; i_hit<mdtSegment->size(); i_hit++) {
576 if (mdtSegment->at(i_hit).isOutlier>1)
continue;
578 double Z = mdtSegment->at(i_hit).Z;
579 double R = mdtSegment->at(i_hit).R;
581 if (i_station==3) nsWidth = m_rwidth_Endcapinn_second;
582 if (i_station==4) nsWidth = m_rwidth_Endcapmid_second;
583 if (i_station==5) nsWidth = m_rwidth_Endcapout_second;
585 double nbw = aw*Z+(avR[i_station]-aw*avZ[i_station]);
587 if ( R>(nbw-nsWidth) && R<(nbw+nsWidth) ) {
588 mdtSegment->at(i_hit).isOutlier = 0;
590 mdtSegment->at(i_hit).isOutlier = 2;
595 for (
unsigned int i_hit=0; i_hit<mdtSegment->size(); i_hit++) {
597 unsigned int i_layer =mdtSegment->at(i_hit).Layer;
598 if (mdtSegment->at(i_hit).isOutlier>1)
continue;
599 if (i_layer > i_layer_max) i_layer_max = i_layer;
602 MdtLayerHits_index[i_station][i_layer].push_back(i_hit);
605 std::vector<unsigned int> Ly_1st;
607 std::vector<float> Line_A;
608 std::vector<float> Line_B;
609 std::vector<float> Line_Chi2;
610 std::vector<unsigned int> Line_count;
611 std::vector<float> Line_Xor;
612 std::vector<float>Line_Yor ;
613 std::vector<float>Line_Amid;
614 std::vector<float>Line_Xmid;
615 std::vector<float>Line_Ymid;
616 std::vector<float> Line_sum_Z;
617 std::vector<float> Line_sum_R;
618 std::vector<float> Line_rm;
619 std::vector<float> Line_phim;
620 std::vector<float> Maxlayers_A;
621 std::vector<float> Maxlayers_B;
622 std::vector<float> Maxlayers_Chi2;
623 std::vector<float> Maxlayers_RESI;
624 float Maxlayers_Phim = 0;
625 float Maxlayers_R = 0;
626 float Maxlayers_Z = 0;
627 float Maxlayers_Xor = 0;
628 float Maxlayers_Yor = 0;
629 float Maxlayers_PChi2 = 0;
633 Maxlayers_Chi2.clear();
634 Maxlayers_RESI.clear();
636 for (
unsigned int i_layer=0; i_layer<=i_layer_max; i_layer++) {
640 Ly_1st.push_back(i_layer);
643 const int real_layer= Ly_1st.size();
644 std::vector<std::vector<unsigned int> > Ly_flg;
646 for (
int pr=real_layer;
pr>=3;
pr--) {
664 findLayerCombination(Ly_1st, real_layer,
pr,Ly_flg, total_cp);
666 for (
unsigned int i=0;
i<Ly_flg.size();
i++) {
668 std::vector<std::vector<int> >tID;
670 std::vector<std::vector<int> >tIndex;
676 for (
int ti=0; ti<8; ti++) {
677 for (
int tj=0; tj<2; tj++) {
679 tubeindex[ti][tj] = 0;
683 for (
unsigned int j=0; j<Ly_flg[
i].size(); j++) {
685 int i_layer = Ly_flg[
i][j];
686 std::vector<int> tid;
688 std::vector<int> tindex;
692 float tube_1st = 999999.;
693 float tube_2nd = 999999.;
695 int layer_2nd = 9999;
697 for (
unsigned int i_hit=0; i_hit<
MdtLayerHits[i_station][i_layer]; i_hit++) {
699 unsigned int i_index = MdtLayerHits_index[i_station][i_layer].at(i_hit);
701 if (mdtSegment->at(i_index).isOutlier>1)
continue;
703 float nbw3 = (mdtSegment->at(i_index).Z)*(aw) + (avR[i_station]-(aw)*avZ[i_station]) ;
704 float dis_tube = std::abs(std::abs(nbw3-mdtSegment->at(i_index).R)- mdtSegment->at(i_index).DriftSpace);
706 if (dis_tube<tube_1st) {
708 layer_2nd = layer_1st;
711 }
else if (dis_tube<tube_2nd) {
717 if ( layer_1st != 9999 ) {
718 mdtSegment->at(layer_1st).isOutlier = 0;
720 tindex.push_back(layer_1st);
723 if ( layer_2nd != 9999 ) {
724 mdtSegment->at(layer_2nd).isOutlier = 1;
726 tindex.push_back(layer_2nd);
730 tIndex.push_back(tindex);
733 for (
unsigned int ti=0; ti<tID.size();ti++) {
734 for (
unsigned int tj=0; tj<tID[ti].size();tj++) {
735 tube_ID[ti][tj] = tID[ti][tj];
736 tubeindex[ti][tj] = tIndex[ti][tj];
740 std::vector<int> isg;
741 std::vector<int> hitarray;
743 int ntry = (
int)floor(
pow(2.,
pr))-1;
745 for (
int ntryi=0; ntryi<=ntry; ntryi++) {
751 for (
int ntryj=1; ntryj<=
pr; ntryj++) {
753 int Isg = (ntryi&(
int)
pow(2.,ntryj-1))? 1 : 0;
755 if (tube_ID[ntryj-1][Isg] != 0) yhit = 1;
756 sumid = sumid * yhit;
760 for (
unsigned int tt=0;
tt<isg.size();
tt++) {
761 int tindex = tubeindex[
tt][isg[
tt]];
762 hitarray.push_back(tindex);
773 float sum_Z_used = 0.;
774 float sum_R_used = 0.;
776 if (hitarray.size()==0)
continue;
778 for (itMdtHit=mdtSegment->begin(); itMdtHit!=mdtSegment->end(); ++itMdtHit) {
782 if(mdtSegment->at(hit_index).isOutlier>1)
continue;
788 for (
unsigned int j=0; j<hitarray.size(); j++) {
790 if (hitarray[j]==hit_index) {
801 rm = itMdtHit->cYmid;
802 Amid = itMdtHit->cAmid;
803 Xmid = itMdtHit->cXmid;
804 Ymid = itMdtHit->cYmid;
812 phim = itMdtHit->cPhip;
813 sigma = (std::abs(itMdtHit->DriftSigma) >
ZERO_LIMIT)? itMdtHit->DriftSigma: SIGMA;
816 std::abs(itMdtHit->DriftSpace) < DRIFTSPACE_LIMIT &&
817 std::abs(itMdtHit->DriftTime) >
ZERO_LIMIT ) {
823 itMdtHit->DriftSpace: SetDriftSpace(itMdtHit->DriftTime, itMdtHit->R, itMdtHit->Z, phim, phiDir);
828 if (i_station==3) i_st = 0;
829 if (i_station==4) i_st = 1;
830 if (i_station==5) i_st = 2;
831 pbFitResult.
IDMEA[
count] = i_st*10 + itMdtHit->Layer;
838 sum_Z_used = sum_Z_used + itMdtHit->Z;
839 sum_R_used = sum_R_used + itMdtHit->R;
846 for (
int i=0;
i<pbFitResult.
NPOI;
i++) {
848 <<
i <<
"/" << pbFitResult.
XILIN[
i] <<
"/" << pbFitResult.
YILIN[
i]
849 <<
"/" << pbFitResult.
RILIN[
i] <<
"/" << pbFitResult.
WILIN[
i]);
852 if (
count >= MIN_MDT_FOR_FIT) {
859 for (
int cand=0; cand<6; cand++) {
862 Line_A.push_back(1/pbFitResult.
SlopeCand[cand]);
864 Line_Chi2.push_back(pbFitResult.
Chi2Cand[cand]);
865 Line_count.push_back(pbFitResult.
NPOI);
866 Line_Xor.push_back(Xor);
867 Line_Yor .push_back(Yor);
868 Line_Amid.push_back(Amid);
869 Line_Xmid.push_back(Xmid);
870 Line_Ymid.push_back(Ymid);
871 Line_sum_Z.push_back(sum_Z_used/
count);
872 Line_sum_R.push_back(sum_R_used/
count);
873 Line_rm.push_back(rm);
874 Line_phim.push_back(phim);
881 if (Line_Chi2.size()==0)
continue;
883 std::multimap<float, int>chi_map;
885 std::vector<float> t_A;
886 std::vector<float> t_B;
887 std::vector<float> t_Chi2;
888 std::vector<float> t_count;
889 std::vector<float> t_Xor;
890 std::vector<float> t_Yor;
891 std::vector<float> t_Amid;
892 std::vector<float> t_Xmid;
893 std::vector<float> t_Ymid;
894 std::vector<float> t_sum_Z;
895 std::vector<float> t_sum_R;
896 std::vector<float> t_rm;
897 std::vector<float> t_phim;
913 for (
unsigned int ir=0;
ir<Line_Chi2.size();
ir++) chi_map.insert(std::make_pair(Line_Chi2.at(
ir),
ir));
916 t_A.push_back(Line_A.at(jt->second));
917 t_B.push_back(Line_B.at(jt->second));
918 t_Chi2.push_back(Line_Chi2.at(jt->second));
919 t_count.push_back(Line_count.at(jt->second));
920 t_Xor.push_back(Line_Xor.at(jt->second));
921 t_Yor.push_back(Line_Yor.at(jt->second));
922 t_Amid.push_back(Line_Amid.at(jt->second));
923 t_Xmid.push_back(Line_Xmid.at(jt->second));
924 t_Ymid.push_back(Line_Ymid.at(jt->second));
925 t_sum_Z.push_back(Line_sum_Z.at(jt->second));
926 t_sum_R.push_back(Line_sum_R.at(jt->second));
927 t_rm.push_back(Line_rm.at(jt->second));
928 t_phim.push_back(Line_phim.at(jt->second));
930 Maxlayers_A.push_back(Line_A.at(jt->second));
931 Maxlayers_B.push_back(Line_B.at(jt->second));
932 Maxlayers_Chi2.push_back(Line_Chi2.at(jt->second));
936 superPoint->
Npoint = t_count[0];
937 if(i_station==4 &&
pr==real_layer){
938 Maxlayers_Z = t_sum_Z[0];
939 Maxlayers_R = t_A[0]*t_sum_Z[0]+t_B[0];
940 Maxlayers_Phim = t_phim[0];
941 Maxlayers_Xor = t_Xor[0];
942 Maxlayers_Yor = t_Yor[0];
943 Maxlayers_PChi2 = pbFitResult.
PCHI2;
944 Maxlayers_N = t_count[0];
947 if (s_address == -1) {
950 superPoint->
Z = t_sum_Z[0];
951 superPoint->
R = t_A[0]*t_sum_Z[0]+t_B[0];
952 superPoint->
Alin =t_A[0];
953 superPoint->
Blin =t_B[0];
956 superPoint->
Phim = t_phim[0];
957 superPoint->
Xor = t_Xor[0];
958 superPoint->
Yor = t_Yor[0];
959 superPoint->
Chi2 = t_Chi2[0];
964 for (
int cand=0; cand<6; cand++) {
968 superPoint->
Chi2Cand[cand] = t_Chi2[cand];
976 if ((i_station == 3 || i_station == 5) &&
pr==4 && Chbest > m_endcapmid_mdt_chi2_limit) {
983 for (
int cand=0; cand<6; cand++) {
991 if (Chbest<m_endcapmid_mdt_chi2_limit){
993 ATH_MSG_DEBUG(
"... Superpoint chamber/s_address/count/R/Z/Alin/Blin/Phim/Xor/Yor/Chi2/PChi2="
994 << i_station <<
"/" << s_address <<
"/" <<
count <<
"/"
995 << superPoint->
R <<
"/" << superPoint->
Z <<
"/" << superPoint->
Alin <<
"/"
996 << superPoint->
Blin <<
"/" << superPoint->
Phim <<
"/" << superPoint->
Xor <<
"/"
997 << superPoint->
Yor <<
"/" << superPoint->
Chi2 <<
"/" << superPoint->
PChi2);
1001 if(i_station==4 && Maxlayers_A.size()>0){
1002 superPoint->
Npoint = Maxlayers_N;
1003 superPoint->
Z = Maxlayers_Z;
1004 superPoint->
R = Maxlayers_R;
1005 superPoint->
Alin =Maxlayers_A[0];
1006 superPoint->
Blin =Maxlayers_B[0];
1007 superPoint->
Phim = Maxlayers_Phim;
1008 superPoint->
Xor = Maxlayers_Xor;
1009 superPoint->
Yor = Maxlayers_Yor;
1010 superPoint->
Chi2 = Maxlayers_Chi2[0];
1011 superPoint->
PChi2 = Maxlayers_PChi2;
1012 for (
int cand=0; cand<6; cand++) {
1013 if (std::abs(Maxlayers_A[cand]) >
ZERO_LIMIT ) {
1014 superPoint->
SlopeCand[cand] = Maxlayers_A[cand];
1016 superPoint->
Chi2Cand[cand] = Maxlayers_Chi2[cand];
1032 const unsigned int MAX_STATION = 8;
1039 float Chi2_cand[8][6];
1043 const int middle = 4;
1044 const int outer = 5;
1046 for (
unsigned int i_station=4; i_station<MAX_STATION; i_station++) {
1054 spChi2[i_station]=0.;
1056 for (
unsigned int ci=0; ci<
NCAND; ci++) {
1057 A_cand[i_station][ci] =0.;
1058 B_cand[i_station][ci] =0.;
1059 Chi2_cand[i_station][ci] =0.;
1062 if (i_station<4 || i_station>5)
continue;
1064 superPoint = &(trackPattern.
superPoints[i_station]);
1065 aw[i_station] = muonRoad.
aw[i_station][0];
1067 spZ[i_station] = superPoint->
Z;
1068 spR[i_station] = superPoint->
R;
1069 spA[i_station] = superPoint->
Alin;
1070 spB[i_station] = superPoint->
Blin;
1072 for (
unsigned int cand=0; cand<
NCAND; cand++) {
1073 A_cand[i_station][cand] = superPoint->
SlopeCand[cand];
1075 Chi2_cand[i_station][cand] = superPoint->
Chi2Cand[cand];
1080 float test_diff = 1.e25;
1081 float best_diff = 1.e25;
1082 float rmatched = 0.;
1083 float match_midA = 0.;
1084 float match_outA = 0.;
1089 spChi2[middle] = 0.;
1103 for (
int m=0;
m<6;
m++) {
1105 test_diff = std::abs(A_cand[middle][
m] - aw[middle]);
1107 if (test_diff<best_diff) {
1108 best_diff = test_diff;
1109 rmatched = A_cand[middle][
m];
1110 spB[middle] = B_cand[middle][
m];
1111 spChi2[middle] = Chi2_cand[middle][
m];
1112 spR[middle] = rmatched * spZ[middle] + spB[middle];
1122 if(std::abs(spZ[outer]-spZ[middle]) >
ZERO_LIMIT) sp_line = (spR[outer]-spR[middle])/(spZ[outer]-spZ[middle]);
1124 for (
int t=0;
t<2; ++
t) {
1128 for (
int i=0;
i<6; ++
i) {
1129 if (
t==0) test_diff = std::abs(A_cand[middle][
i] - sp_line);
1130 else if (
t==1) test_diff = std::abs(A_cand[outer][
i] - sp_line);
1131 if (test_diff<best_diff) {
1132 best_diff = test_diff;
1134 match_midA = A_cand[middle][
i];
1135 spB[middle] = B_cand[middle][
i];
1136 spChi2[middle] = Chi2_cand[middle][
i];
1137 spR[middle] = match_midA * spZ[middle] + spB[middle];
1139 match_outA = A_cand[outer][
i];
1140 spB[outer] = B_cand[outer][
i];
1141 spChi2[outer] = Chi2_cand[outer][
i];
1142 spR[outer] = match_outA * spZ[outer] + spB[outer];
1152 spA[middle] = match_midA;
1153 }
else if (std::abs(rmatched) >
ZERO_LIMIT) {
1154 spA[middle] = rmatched;
1157 if (std::abs(match_outA) >
ZERO_LIMIT) spA[outer] = match_outA;
1161 for (
unsigned int i_station=4; i_station<MAX_STATION; i_station++) {
1162 if (i_station<4 || i_station>5)
continue;
1163 superPoint = &(trackPattern.
superPoints[i_station]);
1165 superPoint->
Alin =spA[i_station];
1166 superPoint->
Blin =spB[i_station];
1167 superPoint->
Chi2 =spChi2[i_station];
1180 float MiddleSlope = 0;
1181 float OuterSlope = 0;
1185 for (
int i_station=4; i_station<6; i_station++) {
1192 superPoint = &(trackPattern.
superPoints[chamberID]);
1194 if ( superPoint->
Npoint > 2 && superPoint->
R > 0.) {
1195 if ( i_station==4 ) {
1196 MiddleSlope = superPoint->
Alin;
1200 }
if ( i_station==5 ) {
1201 OuterSlope = superPoint->
Alin;
1209 double mdtpT = std::abs(tgcFitResult.
tgcPT);
1214 }
else if (std::abs(tgcFitResult.
tgcPT)>=8.0 && std::abs(MiddleSlope) >
ZERO_LIMIT) {
1218 mdtpT = (std::abs(tgcFitResult.
tgcPT)>1
e-5)? mdtpT*(tgcFitResult.
tgcPT/std::abs(tgcFitResult.
tgcPT)) : 0;
1223 double extrInnerEta = 0;
1228 muonSA->
setPt(mdtpT);
1236 StatusCode sc = m_backExtrapolator->give_eta_phi_at_vertex(muonSA, eta,sigma_eta,phi,sigma_phi,0.);
1238 if (
sc.isSuccess() ){
1241 extrInnerEta = etaMiddle;
1248 naw =
std::tan(theta)*(std::abs(extrInnerEta)/extrInnerEta);
1262 double nrWidth = m_rwidth_Endcapinn_first;
1265 unsigned int sumN[8];
1267 for (
unsigned int i_st=0; i_st<8;i_st++) {
1274 const int i_station = 3;
1275 int chamberID = i_station;
1277 mdtSegment = &(trackPattern.
mdtSegments[chamberID]);
1278 superPoint = &(trackPattern.
superPoints[chamberID]);
1280 if (mdtSegment->size()==0)
return;
1283 if (std::abs(itMdtHit.DriftSpace) < m_mdt_driftspace_downlimit ||
1284 std::abs(itMdtHit.DriftSpace) > m_mdt_driftspace_uplimit){
1285 itMdtHit.isOutlier = 2;
1289 if (itMdtHit.isOutlier > 1)
continue;
1291 double Z = itMdtHit.Z;
1292 double R = itMdtHit.R;
1293 double nbw = tgc_aw*Z + bw;
1295 if (R>(nbw-nrWidth) && R<(nbw+nrWidth)) {
1296 itMdtHit.isOutlier = 0;
1299 itMdtHit.isOutlier = 2;
1304 if (sumN[i_station]==0)
return;
1306 stationSPFit(mdtSegment, superPoint,pbFitResult, trackPattern.
s_address,i_station, aw, trackPattern.
phiMSDir);
1310 for (
int cand=0; cand<
NCAND; cand++) {
1311 float ds=std::abs(superPoint->
SlopeCand[cand]-aw);
1329 std::vector<std::vector<unsigned int> > &
c,
1332 std::vector<unsigned int>
b(
r,0);
1334 findSubLayerCombination(
a,
n,
r,
b,0,
r,
c,nr);
1345 std::vector<unsigned int> &
b,
1348 std::vector<std::vector<unsigned int> > &
c,
1354 std::vector<unsigned int>
t;
1359 for (
int j = 0;j<
r; j++)
t.push_back(
b[j]);
1366 findSubLayerCombination(
a,
n,
r,
b,
i+1,
num-1,
c,nr);
1379 const float CSPEED = 2.99979e+10;
1380 const float MDT_RED = 0.7;
1381 const float DRIFTVE = 0.00209;
1383 float flyspa,
sy,
x,
y,phip,phis,prt,tof;
1387 if(phim>=
M_PI) {phim = phim - 2*
M_PI;
sy = -1.;}
1389 y =
sy*std::abs(phiDir*
rad*std::sqrt(1./(1.+phiDir*phiDir)));
1391 phip = std::atan2(
y,
x);
1394 if(phim>=phip) prt = - phis/(CSPEED*MDT_RED);
1395 else prt = + phis/(CSPEED*MDT_RED);
1397 flyspa = std::sqrt(
rad*
rad+zeta*zeta+phis*phis);
1398 tof = flyspa/CSPEED;
1400 return (tdr - tof*1.
e+9 - prt*1.
e+9)*DRIFTVE;
1418 int i,j,
k,Ifit,Ntry,IGcur,Jbad;
1419 float Xnor,rlin,Xbad,
test;
1421 pbFitResult.
NDOF = -2;
1425 for(j=0;j<pbFitResult.
NPOI;j++)
if(pbFitResult.
IGLIN[j]>=1) pbFitResult.
NDOF++;
1427 while(pbFitResult.
NDOF>=1) {
1433 &pbFitResult.
PCHI2);
1437 Xnor = 1. / std::sqrt(1. + pbFitResult.
ALIN * pbFitResult.
ALIN);
1439 for(
i=0;
i<pbFitResult.
NPOI;
i++) pbFitResult.
RESI[
i] = 0.;
1441 for(j=0;j<pbFitResult.
NPOI;j++) {
1443 pbFitResult.
DISTJ[j] = (pbFitResult.
ALIN * pbFitResult.
XILIN[j] + pbFitResult.
BLIN - pbFitResult.
YILIN[j]) * Xnor;
1444 IGcur = std::abs(pbFitResult.
IGLIN[j])%100;
1447 pbFitResult.
RESI[j] = pbFitResult.
DISTJ[j];
1448 }
else if (IGcur==2) {
1449 rlin = (pbFitResult.
DISTJ[j]>=0.) ? pbFitResult.
RILIN[j] : -pbFitResult.
RILIN[j];
1450 pbFitResult.
RESI[j] = pbFitResult.
DISTJ[j] - rlin;
1451 }
else if (IGcur==3) {
1452 pbFitResult.
RESI[j] = pbFitResult.
DISTJ[j] - pbFitResult.
RILIN[j];
1456 if(pbFitResult.
PCHI2>=0.01||Ifla==1)
return Ifit;
1458 if (pbFitResult.
NDOF<=1||Ntry>=6) {
1461 pbFitResult.
NDOF = 0;
1462 pbFitResult.
ALIN = 0.;
1463 pbFitResult.
BLIN = 0.;
1464 pbFitResult.
CHI2 = -1.;
1465 pbFitResult.
PCHI2 = -.5;
1473 for(j=0;j<pbFitResult.
NPOI;j++) pbFitResult.
IGLIN[j] = - std::abs(pbFitResult.
IGLIN[j])%100;
1481 for (j=0;j<pbFitResult.
NPOI;j++) {
1482 if (pbFitResult.
IGLIN[j]>=1) {
1493 pbFitResult.
IGLIN[Jbad] = - pbFitResult.
IGLIN[Jbad] - 100;
1494 pbFitResult.
NDOF = pbFitResult.
NDOF - 1;
1547 float *
A,
float *B,
float DAB[2][2],
float *Chi2,
float *Pchi2)
const
1550 float RRi[
NMEAMX],WIlim,CHbest,Abest,Bbest;
1551 float A0,B0,SAA,SBB,SAB,Square,Aj,Bj;
1552 int i,j,jj,Ntrue,Ngood,Nbest,Ntry,NOgo,Isig,Iflg,Ksign[4];
1558 if (Nmeas<=2||Nmeas>=
NMEAMX+1)
return;
1567 for(
i=0;
i<4;
i++) Ksign[
i] = -1;
1569 for (j=0;j<Nmeas;j++) {
1572 WIlim = (WIlim>=WI[j])? WIlim : WI[j];
1576 if(Ntrue<=2)
return;
1578 WIlim = 0.1 * WIlim;
1580 for(j=0;j<Nmeas;j++)
if(IG[j]>=1&&WI[j]>=WIlim) Ngood++;
1582 for (j=0;j<Nmeas;j++) {
1583 if (IG[j]>=1&&(WI[j]>=WIlim||Ngood<=3)) {
1586 }
else if (Ksign[1]==-1) {
1588 }
else if(Ksign[2]==-1) {
1590 }
else if(Ksign[3]==-1) {
1594 Ksign[2] = Ksign[3];
1602 Xline(XI,YI,WI,IG,Nmeas,&A0,&B0,&SAA,&SBB,&SAB,&Square);
1609 Ntry = (
int)floor(
pow(2.,Nbest)) - 1;
1611 for (j=0;j<=Ntry;j++) {
1613 for (jj=1;jj<=Nbest;jj++) {
1614 Isig = (j&(
int)
pow(2.,jj-1))? 1 : 0;
1616 Iflg = IG[Ksign[jj-1]];
1617 Igg[Ksign[jj-1]] = Iflg;
1618 RRi[Ksign[jj-1]] = RI[Ksign[jj-1]];
1620 Igg[Ksign[jj-1]] = 3;
1621 if (Isig==1) RRi[Ksign[jj-1]] = - RI[Ksign[jj-1]];
1622 }
else if (Isig==1) {
1630 Circfit(Nmeas,XI,YI,RRi,WI,Igg,&Aj,&Bj,DAB,Chi2);
1632 if (*Chi2>=0.0&&*Chi2<=CHbest) {
1643 Circfit(Nmeas,XI,YI,RI,WI,IG,
A,B,DAB,Chi2);
1647 *Pchi2 = TMath::Prob(*Chi2, Ntrue2);
1666 float *
A,
float *B,
float DAB[2][2],
float *Chi2)
const
1670 float SAA,SAB,SBB,Square;
1686 Xnor = 1. / std::sqrt(1. + *
A * *
A);
1689 for(j=0;j<Nmeas;j++) {
1695 }
else if(IG[j]==2) {
1696 if(*
A * XI[j] + *B - YI[j]>=0.) Epsi = 1.0;
1698 XX[j] = XI[j] - Epsi * Xnor * std::abs(RI[j]) * *
A;
1699 YY[j] = YI[j] + Epsi * Xnor * std::abs(RI[j]);
1700 }
else if(IG[j]==3) {
1701 XX[j] = XI[j] - Xnor * RI[j] * *
A;
1702 YY[j] = YI[j] + Xnor * RI[j];
1706 Xline(XX,YY,WI,IG,Nmeas,
A,B,&SAA,&SBB,&SAB,&Square);
1707 if(Square<=0.)
break;
1708 Test = ((Aold-*
A)*(Aold-*
A))/ SAA + ((Bold-*B)*(Bold-*B))/ SBB;
1710 }
while(
Test>=Toll&&Niter<=20);
1731 float *
A,
float *B,
float *SAA,
float *SBB,
float *SAB,
float *Square)
const
1735 float S1,SX,SY,SXX,SXY,SYY,Deter,DY;
1748 SX = SX +
W[j] *
X[j];
1749 SY = SY +
W[j] *
Y[j];
1750 SXX = SXX +
W[j] *
X[j] *
X[j];
1751 SXY = SXY +
W[j] *
X[j] *
Y[j];
1752 SYY = SYY +
W[j] *
Y[j] *
Y[j];
1756 Deter =
S1 * SXX - SX * SX;
1759 *
A = (
S1 * SXY - SX * SY) / Deter;
1760 *B = (SY * SXX - SX * SXY) / Deter;
1763 *SAB = - SX / Deter;
1766 if(
S1 * SXY - SX * SY > 0.) {
1771 *B = SY/
S1 - SX/
S1 * *
A;
1779 DY =(
Y[j] - *
A *
X[j] - *B)/std::sqrt(1 + *
A * *
A);
1781 *Square = *Square +
W[j] * DY * DY;
1790 float *
A,
float *B,
float DAB[2][2],
float *Chi2,
float *Pchi2,
1791 float *SlopeCand,
float *InterceptCand,
float *Chi2Cand)
const
1793 float RRi[
NMEAMX],WIlim,CHbest,Abest,Bbest;
1794 float A0,B0,SAA,SBB,SAB,Square,Aj,Bj;
1795 int i,j,jj,Ntrue,Ngood,Nbest,Ntry,NOgo,Isig,Iflg,Ksign[4];
1798 std::vector<float> st_chi2;
1799 std::vector<float> st_A;
1800 std::vector<float> st_B;
1802 for (
int ii=0; ii<
NCAND; ii++) {
1804 InterceptCand[ii] = 0.;
1811 if (Nmeas<=2||Nmeas>=
NMEAMX+1)
return;
1820 for (
i=0;
i<4;
i++) Ksign[
i] = -1;
1822 for (j=0;j<Nmeas;j++) {
1825 WIlim = (WIlim>=WI[j])? WIlim : WI[j];
1829 if (Ntrue<=2)
return;
1831 WIlim = 0.1 * WIlim;
1833 for(j=0;j<Nmeas;j++)
if(IG[j]>=1&&WI[j]>=WIlim) Ngood++;
1835 for (j=0;j<Nmeas;j++) {
1836 if (IG[j]>=1&&(WI[j]>=WIlim||Ngood<=3)) {
1840 }
else if (Ksign[1]==-1) {
1842 }
else if (Ksign[2]==-1) {
1844 }
else if(Ksign[3]==-1) {
1848 Ksign[2] = Ksign[3];
1856 Xline(XI,YI,WI,IG,Nmeas,&A0,&B0,&SAA,&SBB,&SAB,&Square);
1859 st_A.clear(); st_B.clear(); st_chi2.clear();
1864 Ntry = (
int)floor(
pow(2.,Nbest)) - 1;
1866 for (j=0;j<=Ntry;j++) {
1870 for (jj=1;jj<=Nbest;jj++) {
1871 Isig = (j&(
int)
pow(2.,jj-1))? 1 : 0;
1873 Iflg = IG[Ksign[jj-1]];
1874 Igg[Ksign[jj-1]] = Iflg;
1875 RRi[Ksign[jj-1]] = RI[Ksign[jj-1]];
1879 Igg[Ksign[jj-1]] = 3;
1881 if (Isig==1) RRi[Ksign[jj-1]] = - RI[Ksign[jj-1]];
1883 }
else if (Isig==1) {
1894 Circfit(Nmeas,XI,YI,RRi,WI,Igg,&Aj,&Bj,DAB,Chi2);
1895 Circfit(Nmeas,XI,YI,RI,WI,IG,&Aj,&Bj,DAB,Chi2);
1896 st_A.push_back(Aj); st_B.push_back(Bj); st_chi2.push_back(*Chi2);
1898 if (*Chi2>=0.0&&*Chi2<=CHbest) {
1906 std::multimap<float, int>chi_map;
1908 std::vector<float> t_A;
1909 std::vector<float> t_B;
1910 std::vector<float> t_chi2;
1915 for (
unsigned int ir=0;
ir<st_chi2.size();
ir++) chi_map.insert(std::make_pair(st_chi2.at(
ir),
ir));
1918 t_A.push_back(st_A.at(jt->second));
1919 t_B.push_back(st_B.at(jt->second));
1920 t_chi2.push_back(st_chi2.at(jt->second));
1923 for (
int nv=0; nv<6; nv++) {
1924 SlopeCand[nv] = t_A[nv];
1925 InterceptCand[nv] = t_B[nv];
1926 Chi2Cand[nv] = t_chi2[nv];
1932 Circfit(Nmeas,XI,YI,RI,WI,IG,
A,B,DAB,Chi2);
1936 *Pchi2 = TMath::Prob(*Chi2, Ntrue2);