19 static const double TanM1p5 =
std::tan(-1.5/360.*2.*
M_PI);
20 static const double TanP1p5 =
std::tan(1.5/360.*2.*
M_PI);
21 static const double CosM1p5 =
std::cos(-1.5/360.*2.*
M_PI);
22 static const double CosP1p5 =
std::cos(1.5/360.*2.*
M_PI);
23 static const double SinM1p5 =
std::sin(-1.5/360.*2.*
M_PI);
24 static const double SinP1p5 =
std::sin(1.5/360.*2.*
M_PI);
27 const std::string&
name,
39 ATH_MSG_DEBUG(
"NswStationFitter::findSuperPoints() was called.");
45 ATH_CHECK( findStgcHitsInSegment(stgcHits) );
48 bool isLargeStgc =
false;
49 bool isSmallStgc =
false;
50 if(stgcHits.size() != 0){
51 for(
unsigned int iHit = 0; iHit < stgcHits.size(); iHit++){
52 if(stgcHits.at(iHit).isOutlier == 0){
53 if(stgcHits.at(iHit).stationName == 58) isLargeStgc =
true;
54 else if(stgcHits.at(iHit).stationName == 57) isSmallStgc =
true;
55 if(isLargeStgc && isSmallStgc)
continue;
59 if(mmHits.size() != 0){
60 for(
unsigned int iHit = 0; iHit < mmHits.size(); iHit++){
61 if(mmHits.at(iHit).isOutlier == 0){
63 if(mmHits.at(iHit).stationName == 55) mmHits.at(iHit).isOutlier = 1;
64 }
else if(isSmallStgc) {
65 if(mmHits.at(iHit).stationName == 56) mmHits.at(iHit).isOutlier = 1;
70 ATH_CHECK( MakeSegment(trackPattern,stgcHits) );
71 ATH_CHECK( MakeSegment(trackPattern,mmHits) );
82 return StatusCode::SUCCESS;
91 selectedStgcHits.clear();
96 double phiMin = p_roids->
phi() - 1.;
97 double phiMax = p_roids->
phi() + 1.;
98 if( phiMin >
M_PI ) phiMin -= 2*
M_PI;
99 if( phiMax >
M_PI ) phiMax -= 2*
M_PI;
100 if( phiMin < -1.*
M_PI ) phiMin += 2*
M_PI;
101 if( phiMax < -1.*
M_PI ) phiMax += 2*
M_PI;
103 bool inRoiEta, inRoiPhi, inRoi;
107 if (stgcHits.size()>0) {
108 for (iHit = 0; iHit < stgcHits.size(); iHit++){
113 inRoiEta =
false; inRoiPhi =
false; inRoi =
false;
115 if ( phiMin<=phiMax && phiMin<=hit.
phi && hit.
phi<=phiMax ) inRoiPhi=
true;
116 if ( phiMin>phiMax && (phiMin<=hit.
phi || hit.
phi<=phiMax)) inRoiPhi=
true;
117 if ( inRoiEta && inRoiPhi ) inRoi =
true;
123 selectedStgcHits.push_back(hit);
128 stgcHits = selectedStgcHits;
130 return StatusCode::SUCCESS;
139 selectedMmHits.clear();
144 double phiMin = p_roids->
phi() - 1.;
145 double phiMax = p_roids->
phi() + 1.;
146 if( phiMin >
M_PI ) phiMin -= 2*
M_PI;
147 if( phiMax >
M_PI ) phiMax -= 2*
M_PI;
148 if( phiMin < -1.*
M_PI ) phiMin += 2*
M_PI;
149 if( phiMax < -1.*
M_PI ) phiMax += 2*
M_PI;
151 bool inRoiEta, inRoiPhi, inRoi;
155 if (mmHits.size()>0) {
156 for (iHit = 0; iHit < mmHits.size(); iHit++){
162 inRoiEta =
false; inRoiPhi =
false; inRoi =
false;
164 if ( phiMin<=phiMax && phiMin<=hit.
phi && hit.
phi<=phiMax ) inRoiPhi=
true;
165 if ( phiMin>phiMax && (phiMin<=hit.
phi || hit.
phi<=phiMax)) inRoiPhi=
true;
166 if ( inRoiEta && inRoiPhi ) inRoi =
true;
173 selectedMmHits.push_back(hit);
178 mmHits = selectedMmHits;
180 return StatusCode::SUCCESS;
190 if(stgcHits.size()==0 && mmHits.size()==0)
191 return StatusCode::SUCCESS;
194 float rWeightedCenter=0.,zWeightedCenter=0.,phiWeightedCenter=0.;
195 int totNumRWeightedCenter=0.,totNumZWeightedCenter=0.,totNumPhiWeightedCenter=0.;
196 float localPhiCenter=0., localPhi;
200 if (stgcHits.size()>0) {
201 for (iHit = 0; iHit < stgcHits.size(); iHit++){
209 localPhiCenter +=
M_PI/8.;
213 localPhi = hit.
phi - localPhiCenter;
214 if (localPhi >
M_PI) localPhi -= 2.0*
M_PI;
215 if (localPhi < -1.*
M_PI) localPhi += 2.0*
M_PI;
218 rWeightedCenter += hit.
r;
219 totNumRWeightedCenter += 1;
222 phiWeightedCenter += localPhi;
223 totNumPhiWeightedCenter += 1;
225 zWeightedCenter += hit.
z;
226 totNumZWeightedCenter += 1;
231 if (mmHits.size()>0) {
232 for (iHit = 0; iHit < mmHits.size(); iHit++){
238 localPhiCenter +=
M_PI/8.;
242 localPhi = hit.
phi - localPhiCenter;
243 if (localPhi >
M_PI) localPhi -= 2.0*
M_PI;
244 if (localPhi < -1.*
M_PI) localPhi += 2.0*
M_PI;
249 if (totNumRWeightedCenter!=0) rWeightedCenter /= totNumRWeightedCenter;
250 if (totNumPhiWeightedCenter!=0) phiWeightedCenter /= totNumPhiWeightedCenter;
251 if (totNumZWeightedCenter!=0) zWeightedCenter /= totNumZWeightedCenter;
256 superPoint->
R = rWeightedCenter;
257 superPoint->
Phim = phiWeightedCenter+localPhiCenter;
258 superPoint->
Z = zWeightedCenter;
259 superPoint->
Npoint = totNumZWeightedCenter;
260 if (zWeightedCenter!=0) superPoint->
Alin = rWeightedCenter/zWeightedCenter;
261 superPoint->
Blin = 0.;
263 ATH_MSG_DEBUG(
"Nsw Super Point r/phi/z/slope = "<<superPoint->
R<<
"/"<<superPoint->
Phim<<
"/"<<superPoint->
Z<<
"/"<<superPoint->
Alin);
265 return StatusCode::SUCCESS;
271 if(stgcHits.size() == 0)
return StatusCode::SUCCESS;
273 for(
unsigned int iHit = 0; iHit < stgcHits.size(); iHit++){
274 if(stgcHits.at(iHit).isOutlier == 0){
276 stgcHits.at(iHit).isOutlier = 1;
279 if(hitsInRoad == 0)
return StatusCode::SUCCESS;
282 ATH_MSG_DEBUG(
"Number of STGC hits is too small, at least 9 hits required : "<<hitsInRoad<<
" hits");
283 return StatusCode::SUCCESS;
284 }
else if(hitsInRoad > 100) {
285 ATH_MSG_WARNING(
"Number of STGC hits is too large, at most 100 hits allowed : "<<hitsInRoad<<
" hits");
286 return StatusCode::SUCCESS;
290 std::array<std::vector<int>, 8> strHitIdByLayer;
291 std::array<std::vector<int>, 8> wireHitIdByLayer;
292 for(
unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit){
293 if(stgcHits.at(iHit).isOutlier != 1)
continue;
294 int layerNumber = stgcHits.at(iHit).layerNumber;
295 if (layerNumber > 7) {
299 if(stgcHits.at(iHit).channelType == 1){
300 strHitIdByLayer[layerNumber].push_back(iHit);
301 }
else if(stgcHits.at(iHit).channelType == 2){
302 wireHitIdByLayer[layerNumber].push_back(iHit);
306 <<
" " << strHitIdByLayer[1].
size()
307 <<
" " << strHitIdByLayer[2].
size()
308 <<
" " << strHitIdByLayer[3].
size()
309 <<
" " << strHitIdByLayer[4].
size()
310 <<
" " << strHitIdByLayer[5].
size()
311 <<
" " << strHitIdByLayer[6].
size()
312 <<
" " << strHitIdByLayer[7].
size());
314 <<
" " << wireHitIdByLayer[1].
size()
315 <<
" " << wireHitIdByLayer[2].
size()
316 <<
" " << wireHitIdByLayer[3].
size()
317 <<
" " << wireHitIdByLayer[4].
size()
318 <<
" " << wireHitIdByLayer[5].
size()
319 <<
" " << wireHitIdByLayer[6].
size()
320 <<
" " << wireHitIdByLayer[7].
size());
322 std::vector< std::array<int, 8> > strHitIds;
323 findSetOfStgcHitIds(stgcHits, strHitIdByLayer, strHitIds);
324 std::vector< std::array<int, 8> > wireHitIds;
325 findSetOfStgcHitIds(stgcHits, wireHitIdByLayer, wireHitIds);
326 ATH_MSG_DEBUG(
"@@STGC@@ strip wire " << strHitIds.size() <<
" " << wireHitIds.size());
328 bool isLargeStrip =
false;
329 bool isSmallStrip =
false;
330 for (
unsigned int iHit = 0; iHit < strHitIds.size(); ++iHit) {
331 std::array<int, 8> hitIds = strHitIds.at(iHit);
332 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
333 if (hitIds[iLayer] != -1) {
334 if(stgcHits.at(hitIds[iLayer]).stationName == 58) isLargeStrip =
true;
335 else if(stgcHits.at(hitIds[iLayer]).stationName == 57) isSmallStrip =
true;
336 stgcHits.at(hitIds[iLayer]).isOutlier = 0;
340 for (
unsigned int iHit = 0; iHit < wireHitIds.size(); ++iHit) {
341 std::array<int, 8> hitIds = wireHitIds.at(iHit);
342 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
343 if (hitIds[iLayer] != -1) {
345 if(stgcHits.at(hitIds[iLayer]).stationName == 58){
346 stgcHits.at(hitIds[iLayer]).isOutlier = 0;
349 else if(isSmallStrip){
350 if(stgcHits.at(hitIds[iLayer]).stationName == 57){
351 stgcHits.at(hitIds[iLayer]).isOutlier = 0;
357 return StatusCode::SUCCESS;
363 std::vector<std::array<int, 8>>& hitIdsCandidate)
const
365 double NSWCenterZ = 7526.329;
369 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
370 if ( hitIdByLayer[iLayer].
size() > 0) {
371 side = std::abs(stgcHits.at(hitIdByLayer[iLayer].at(0)).z)/stgcHits.at(hitIdByLayer[iLayer].at(0)).z;
372 if ( stgcHits.at(hitIdByLayer[iLayer].at(0)).channelType == 1 ){
378 NSWCenterZ = NSWCenterZ *
side;
380 std::array<std::vector<unsigned long int>,4> hitIdsInTwo;
381 std::array<std::vector<double>,4> slopeInTwo;
382 std::array<std::vector<double>,4> interceptInTwo;
385 for(
unsigned int iPair = 0; iPair < 4; ++iPair){
386 unsigned int nHitsInInner = hitIdByLayer[iPair].size();
387 unsigned int nHitsInOuter = hitIdByLayer[iPair+4].size();
388 if ( nHitsInInner > 0xffff-1 || nHitsInOuter > 0xffff-1) {
389 ATH_MSG_WARNING(
"Number of Stgc hits in layers exceeds the limit of (2^16 - 1) : Number of Stgc hits in "<<iPair<<
"th layer = "<< nHitsInInner
390 <<
", Number of Stgc hits in "<<iPair+4<<
"th layer = "<<nHitsInOuter);
391 ATH_MSG_WARNING(
"Number of Stgc hits is limitted to (2^16 - 1) and hits with id more than (2^16 -1) will be trancated.");
392 if (nHitsInInner > 0xffff-1) {nHitsInInner = 0xffff-1;}
393 if (nHitsInOuter > 0xffff-1) {nHitsInOuter = 0xffff-1;}
396 bool foundCounterparts[256] = {};
398 for(
unsigned int iHit = 0; iHit < nHitsInInner; ++iHit){
399 bool foundCounterpart = 0;
404 int iHitId = hitIdByLayer[iPair].at(iHit);
406 r[0] = stgcHits.at(iHitId).r;
407 z[0] = stgcHits.at(iHitId).z;
409 double localPhiCenter;
410 if (stgcHits.at(iHitId).stationPhi<=5) {
411 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHitId).stationPhi-1.);
413 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHitId).stationPhi-9.);
416 if (stgcHits.at(iHitId).stationName == 57){
417 localPhiCenter +=
M_PI/8.;
418 if (stgcHits.at(iHitId).stationPhi == 5) localPhiCenter -= 2. *
M_PI;
421 double phiProj = stgcHits.at(iHitId).phi - localPhiCenter;
422 if (phiProj >
M_PI) phiProj -= 2.0*
M_PI;
423 if (phiProj < -1.*
M_PI) phiProj += 2.0*
M_PI;
424 r[0] = stgcHits.at(iHitId).r;
429 for(
unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
430 int jHitId = hitIdByLayer[iPair+4].at(jHit);
431 double slope, intercept;
433 r[1] = stgcHits.at(jHitId).r;
434 z[1] = stgcHits.at(jHitId).z;
435 slope = (
r[1] -
r[0]) / (
z[1] -
z[0]);
436 intercept = slope*(0. -
z[0]) +
r[0];
438 if(std::abs(slope) < 0.14 || std::abs(slope) > 0.6 || std::abs(intercept) > 300.)
continue;
440 double localPhiCenter;
441 if (stgcHits.at(jHitId).stationPhi<=5) {
442 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(jHitId).stationPhi-1.);
444 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(jHitId).stationPhi-9.);
447 if (stgcHits.at(jHitId).stationName == 57){
448 localPhiCenter +=
M_PI/8.;
449 if (stgcHits.at(jHitId).stationPhi == 5) localPhiCenter -= 2 *
M_PI;
452 double phiProj = stgcHits.at(jHitId).phi - localPhiCenter;
453 if (phiProj >
M_PI) phiProj -= 2.0*
M_PI;
454 if (phiProj < -1.*
M_PI) phiProj += 2.0*
M_PI;
455 r[1] = stgcHits.at(jHitId).r;
457 slope = (
z[0]+
z[1])/2.;
458 intercept = (
r[0]+
r[1])/2.;
462 unsigned int encodedIds = (iHitId<<16) + jHitId;
463 hitIdsInTwo[iPair].push_back(encodedIds);
464 slopeInTwo[iPair].push_back(slope);
465 interceptInTwo[iPair].push_back(intercept);
467 foundCounterpart = 1;
468 foundCounterparts[jHit] = 1;
470 if(!foundCounterpart){
471 unsigned int encodedIds = (iHitId<<16) + 0xffff;
472 hitIdsInTwo[iPair].push_back(encodedIds);
474 slopeInTwo[iPair].push_back(
r[0]/
z[0]);
475 interceptInTwo[iPair].push_back(0.);
477 slopeInTwo[iPair].push_back(
z[0]);
478 interceptInTwo[iPair].push_back(
r[0]);
483 for(
unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
484 if (!foundCounterparts[jHit]) {
485 int jHitId = hitIdByLayer[iPair+4].at(jHit);
486 unsigned int encodedIds = 0xffff0000 + jHitId;
487 hitIdsInTwo[iPair].push_back(encodedIds);
489 slopeInTwo[iPair].push_back(stgcHits.at(jHitId).r/stgcHits.at(jHitId).z);
490 interceptInTwo[iPair].push_back(0.);
492 double localPhiCenter;
493 if (stgcHits.at(jHitId).stationPhi<=5) {
494 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(jHitId).stationPhi-1.);
496 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(jHitId).stationPhi-9.);
498 if (stgcHits.at(jHitId).stationName == 57){
499 localPhiCenter +=
M_PI/8.;
500 if (stgcHits.at(jHitId).stationPhi == 5) localPhiCenter -= 2 *
M_PI;
503 double phiProj = stgcHits.at(jHitId).phi - localPhiCenter;
504 if (phiProj >
M_PI) phiProj -= 2.0*
M_PI;
505 if (phiProj < -1.*
M_PI) phiProj += 2.0*
M_PI;
506 slopeInTwo[iPair].push_back(phiProj);
507 interceptInTwo[iPair].push_back(stgcHits.at(jHitId).r);
513 ATH_MSG_DEBUG(
"@@STGC@@ isStrip= " <<
isStrip <<
" Npairs " << hitIdsInTwo[0].
size() <<
" " << hitIdsInTwo[1].
size() <<
" " << hitIdsInTwo[2].
size() <<
" " << hitIdsInTwo[3].
size());
514 for (
unsigned int iLayer = 0; iLayer < 4; ++ iLayer) {
515 for (
unsigned int iPair = 0; iPair < slopeInTwo[iLayer].size(); ++iPair) {
516 ATH_MSG_DEBUG(
"@@STGC@@ pair fit isStrip= " <<
isStrip <<
" slope= " << slopeInTwo[iLayer].at(iPair) <<
" intercept= " << interceptInTwo[iLayer].at(iPair));
520 std::array<std::vector<unsigned long int>,2> hitIdsInFour;
521 std::array<std::vector<double>,2> slopeInFour;
522 std::array<std::vector<double>,2> interceptInFour;
523 for(
unsigned int iQuad = 0; iQuad < 2; ++iQuad){
524 unsigned int nPairsInInner = hitIdsInTwo[iQuad].size();
525 unsigned int nPairsInOuter = hitIdsInTwo[iQuad+2].size();
527 bool foundCounterparts[0xffff] = {};
528 for(
unsigned int iPair = 0; iPair < nPairsInInner; ++iPair){
529 bool foundCounterpart = 0;
532 slope[0] = slopeInTwo[iQuad].at(iPair);
533 intercept[0] = interceptInTwo[iQuad].at(iPair);
535 for(
unsigned int jPair = 0; jPair < nPairsInOuter; ++jPair){
536 unsigned int ihitIds = hitIdsInTwo[iQuad].at(iPair);
537 unsigned int jhitIds = hitIdsInTwo[iQuad+2].at(jPair);
538 if ( !(((ihitIds>>16 & 0xffff) != 0xffff || (ihitIds & 0xffff) != 0xffff) &&
539 ((jhitIds>>16 & 0xffff) != 0xffff || (jhitIds & 0xffff) != 0xffff )) )
continue;
541 slope[1] = slopeInTwo[iQuad+2].at(jPair);
542 intercept[1] = interceptInTwo[iQuad+2].at(jPair);
545 double spR0 = slope[0] * NSWCenterZ + intercept[0];
546 double spR1 = slope[1] * NSWCenterZ + intercept[1];
547 if(std::abs(spR1 - spR0) > 50.)
continue;
549 if(std::abs(intercept[0]*
std::sin(slope[0]) - intercept[1]*
std::sin(slope[1])) > 100.)
continue;
552 foundCounterpart = 1;
553 foundCounterparts[jPair] = 1;
555 unsigned long int encodedIds = (hitIdsInTwo[iQuad].at(iPair) << 32 ) + hitIdsInTwo[iQuad+2].at(jPair);
556 hitIdsInFour[iQuad].push_back(encodedIds);
557 slopeInFour[iQuad].push_back((slope[1] + slope[0])/2.);
558 interceptInFour[iQuad].push_back((intercept[1] + intercept[0])/2.);
561 if(foundCounterpart)
continue;
562 if((hitIdsInTwo[iQuad].at(iPair)>>16 & 0xffff) == 0xffff || (hitIdsInTwo[iQuad].at(iPair) & 0xffff) == 0xffff)
continue;
564 unsigned long int encodedIds = (hitIdsInTwo[iQuad].at(iPair) << 32 ) + 0xffffffff;
565 hitIdsInFour[iQuad].push_back(encodedIds);
566 slopeInFour[iQuad].push_back(slope[0]);
567 interceptInFour[iQuad].push_back(intercept[0]);
569 for (
unsigned int jPair = 0; jPair < nPairsInOuter; ++jPair) {
570 if(foundCounterparts[jPair])
continue;
571 if((hitIdsInTwo[iQuad+2].at(jPair)>>16 & 0xffff) == 0xffff || (hitIdsInTwo[iQuad+2].at(jPair) & 0xffff) == 0xffff)
continue;
573 unsigned long int encodedIds = (0xffffffff00000000 ) + hitIdsInTwo[iQuad+2].at(jPair);
574 hitIdsInFour[iQuad].push_back(encodedIds);
575 slopeInFour[iQuad].push_back(slopeInTwo[iQuad+2].at(jPair));
576 interceptInFour[iQuad].push_back(interceptInTwo[iQuad+2].at(jPair));
581 for (
unsigned int iLayer = 0; iLayer < 2; ++ iLayer) {
582 for (
unsigned int iQuad = 0; iQuad < slopeInFour[iLayer].size(); ++iQuad) {
583 ATH_MSG_DEBUG(
"@@STGC@@ quad fit isStrip= " <<
isStrip <<
" slope= " << slopeInFour[iLayer].at(iQuad) <<
" intercept= " << interceptInFour[iLayer].at(iQuad));
587 std::vector< std::array<int, 8> > hitIdsInEight;
588 std::vector<double> mseInEight;
590 unsigned int nQuadInInner = hitIdsInFour[0].size();
591 unsigned int nQuadInOuter = hitIdsInFour[1].size();
593 for(
unsigned int iQuad = 0; iQuad < nQuadInInner; ++iQuad){
596 slope[0] = slopeInFour[0].at(iQuad);
597 intercept[0] = interceptInFour[0].at(iQuad);
599 for(
unsigned int jQuad = 0; jQuad < nQuadInOuter; ++jQuad){
600 unsigned long int ihitIds = hitIdsInFour[0].at(iQuad);
601 unsigned long int jhitIds = hitIdsInFour[1].at(jQuad);
602 int nOfLayersWithNoHit = 0;
603 for (
unsigned int iLayer = 0; iLayer < 4; ++iLayer) {
604 if ( (ihitIds>>(3-iLayer)*16 & 0xffff) == 0xffff ) {++nOfLayersWithNoHit;}
605 if ( (jhitIds>>(3-iLayer)*16 & 0xffff) == 0xffff ) {++nOfLayersWithNoHit;}
607 if (nOfLayersWithNoHit > 4)
continue;
609 slope[1] = slopeInFour[1].at(jQuad);
610 intercept[1] = interceptInFour[1].at(jQuad);
613 double spR0 = slope[0] * NSWCenterZ + intercept[0];
614 double spR1 = slope[1] * NSWCenterZ + intercept[1];
615 if(std::abs(spR1 - spR0) > 10. ||
616 std::abs(intercept[1] + intercept[0]) / 2 > 100.)
continue;
618 if(std::abs(intercept[0]*
std::sin(slope[0]) - intercept[1]*
std::sin(slope[1])) > 100. )
continue;
621 std::array<int,8> setOfHitIds = {-1,-1,-1,-1,-1,-1,-1,-1};
622 std::vector<double>
r,
z;
623 for(
unsigned int i = 0;
i < 8; ++
i) {
624 unsigned int iHitId, iLayer = 0;
626 iHitId = (
unsigned int) ((ihitIds>>(3-
i)*16) & 0xffff);
629 iHitId = (
unsigned int) ((jhitIds>>(3-
i%4)*16) & 0xffff);
630 iLayer = (
i-4)+3*(
i%2)+1;
632 if ( iHitId != 0xffff ) {
634 r.push_back(stgcHits.at(iHitId).r);
635 z.push_back(stgcHits.at(iHitId).z);
638 double localPhiCenter;
639 if (stgcHits.at(iHitId).stationPhi<=5) {
640 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHitId).stationPhi-1.);
642 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHitId).stationPhi-9.);
644 if (stgcHits.at(iHitId).stationName == 57){
645 localPhiCenter +=
M_PI/8.;
646 if (stgcHits.at(iHitId).stationPhi == 5) localPhiCenter -= 2 *
M_PI;
649 double phiProj = stgcHits.at(iHitId).phi - localPhiCenter;
650 if (phiProj >
M_PI) phiProj -= 2.0*
M_PI;
651 if (phiProj < -1.*
M_PI) phiProj += 2.0*
M_PI;
652 r.push_back(phiProj);
653 z.push_back(stgcHits.at(iHitId).r);
655 setOfHitIds[iLayer] = iHitId;
656 ATH_MSG_DEBUG(
"@@STGC@@ strip_pos iHitId " << iLayer <<
" " << iHitId);
659 double slopefit=0., interceptfit=99999., mse =-1.;
661 LinearFit(
z,
r,&slopefit,&interceptfit,&mse);
664 for (
unsigned int iHit = 0; iHit <
r.size(); ++iHit){
665 phiavg +=
r.at(iHit);
669 for (
unsigned int iHit = 0; iHit <
r.size(); ++iHit){
673 hitIdsInEight.push_back(setOfHitIds);
674 mseInEight.push_back(mse);
677 if(!hitIdsInEight.size()){
683 std::vector<int> nOctetSegments;
684 std::vector<int> patternStationName;
686 for (
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet) {
688 bool isFirstHit =
true;
689 int hitStationName = 0;
691 int nOctetSegment = 0;
693 std::array<int, 8> tmpOctet = hitIdsInEight.at(iOctet);
694 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
695 if (tmpOctet[iLayer] != -1) {
698 hitStationName = stgcHits.at(tmpOctet[iLayer]).stationName;
701 ATH_MSG_DEBUG(
"@@STGC@@ octet pos isStrip= " <<
isStrip <<
" r= " << stgcHits.at(tmpOctet[iLayer]).r <<
" phi= " << stgcHits.at(tmpOctet[iLayer]).phi <<
" z= " << stgcHits.at(tmpOctet[iLayer]).z);
704 else if(stgcHits.at(tmpOctet[iLayer]).stationName == hitStationName){
705 ATH_MSG_DEBUG(
"@@STGC@@ octet pos isStrip= " <<
isStrip <<
" r= " << stgcHits.at(tmpOctet[iLayer]).r <<
" phi= " << stgcHits.at(tmpOctet[iLayer]).phi <<
" z= " << stgcHits.at(tmpOctet[iLayer]).z);
710 nOctetSegments.push_back(nOctetSegment);
711 patternStationName.push_back(hitStationName);
713 double nOcSegMax = 0;
714 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
715 if(nOctetSegments.at(iOctet) > nOcSegMax){
716 nOcSegMax = nOctetSegments.at(iOctet);
720 double msemin = 1000000.;
721 double mseminWireL = 1000000.;
722 double mseminWireS = 1000000.;
724 std::vector<int> octetIds(2,-1);
727 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
728 if(nOctetSegments.at(iOctet) != nOcSegMax){
731 if( mseInEight.at(iOctet) < msemin) {
732 msemin = mseInEight.at(iOctet);
735 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
736 if(nOctetSegments.at(iOctet) != nOcSegMax)
continue;
737 if(mseInEight.at(iOctet) != msemin){
740 octetIds.push_back(iOctet);
743 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
744 if(patternStationName.at(iOctet) == 58){
745 if( mseInEight.at(iOctet) < mseminWireL) {
746 mseminWireL = mseInEight.at(iOctet);
749 else if(patternStationName.at(iOctet) == 57){
750 if( mseInEight.at(iOctet) < mseminWireS) {
751 mseminWireS = mseInEight.at(iOctet);
755 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
756 if(patternStationName.at(iOctet) == 58){
757 if(mseInEight.at(iOctet) != mseminWireL){
761 else if(patternStationName.at(iOctet) == 57){
762 if(mseInEight.at(iOctet) != mseminWireS){
766 octetIds.push_back(iOctet);
769 for(
unsigned int ids = 0;
ids < octetIds.size();
ids++){
770 if (octetIds.at(
ids) != -1) {
771 hitIdsCandidate.push_back(hitIdsInEight.at(octetIds.at(
ids)));
780 selectedStgcHits.clear();
781 if(stgcHits.size() == 0)
return StatusCode::SUCCESS;
782 for(
unsigned int iHit = 0; iHit < stgcHits.size(); iHit++){
783 if(stgcHits.at(iHit).isOutlier != 0)
continue;
784 selectedStgcHits.push_back(stgcHits.at(iHit));
788 return StatusCode::SUCCESS;
794 selectedMmHits.clear();
795 if(mmHits.size() == 0)
return StatusCode::SUCCESS;
796 for(
unsigned int iHit = 0; iHit < mmHits.size(); iHit++){
797 if(mmHits.at(iHit).isOutlier != 0)
continue;
798 selectedMmHits.push_back(mmHits.at(iHit));
802 return StatusCode::SUCCESS;
807 double* slope,
double* intercept,
double* mse)
const
809 double sumX=0, sumY=0, sumXY=0, sumX2=0;
810 int nHits =
x.size();
812 for (
unsigned int iHit = 0; iHit <
x.size(); ++iHit){
815 sumXY +=
x.at(iHit)*
y.at(iHit);
816 sumX2 +=
x.at(iHit)*
x.at(iHit);
820 if(nHits * sumX2 - (sumX * sumX) >
ZERO_LIMIT) {
821 *slope = (nHits * sumXY - sumX * sumY) / (nHits * sumX2 - (sumX * sumX));
822 *intercept = (sumX2 * sumY - sumXY * sumX) / (nHits * sumX2 - sumX * sumX);
825 *intercept = sumY/nHits;
828 else if(nHits == 1) {
834 for(
unsigned int iHit = 0; iHit<
x.size(); ++iHit){
835 *mse +=
std::pow(
y.at(iHit) - (*slope *
x.at(iHit) + *intercept), 2.0);
837 *mse = *mse / (nHits - 2);
845 std::vector<bool>& isStgc,
double* slope,
double* intercept,
double* mse,
double eta)
const
847 double RmsDeltarEtaStgc[12] = {};
848 double RmsDeltarEtaMm[12]= {};
849 getNswResolution(RmsDeltarEtaStgc, RmsDeltarEtaMm, 12);
851 double weightStgc[12] = {};
852 double weightMm[12] = {};
853 for(
int i_weight=0; i_weight<12; i_weight++){
854 weightStgc[i_weight] = 1/
std::pow(RmsDeltarEtaStgc[i_weight],2);
855 weightMm[i_weight] = 1/
std::pow(RmsDeltarEtaMm[i_weight],2);
860 for(
int iBin=0; iBin<12; iBin++){
861 if(std::abs(eta) >= minEta && std::abs(eta) <
maxEta){
870 double sumX=0, sumY=0, sumXY=0, sumX2=0, sumW=0;
871 int nHits =
x.size();
873 for (
unsigned int iHit = 0; iHit <
x.size(); ++iHit){
875 sumX += weightStgc[weightBin] *
x.at(iHit);
876 sumY += weightStgc[weightBin] *
y.at(iHit);
877 sumXY += weightStgc[weightBin] *
x.at(iHit) *
y.at(iHit);
878 sumX2 += weightStgc[weightBin] *
x.at(iHit) *
x.at(iHit);
879 sumW += weightStgc[weightBin];
881 sumX += weightMm[weightBin] *
x.at(iHit);
882 sumY += weightMm[weightBin] *
y.at(iHit);
883 sumXY += weightMm[weightBin] *
x.at(iHit) *
y.at(iHit);
884 sumX2 += weightMm[weightBin] *
x.at(iHit) *
x.at(iHit);
885 sumW += weightMm[weightBin];
890 if(nHits * sumX2 - (sumX * sumX) >
ZERO_LIMIT) {
891 *slope = (sumW * sumXY - sumX * sumY) / (sumW * sumX2 - (sumX * sumX));
892 *intercept = (sumX2 * sumY - sumXY * sumX) / (sumW * sumX2 - sumX * sumX);
895 *intercept = sumY/nHits;
898 else if(nHits == 1) {
904 for(
unsigned int iHit = 0; iHit<
x.size(); ++iHit){
905 *mse +=
std::pow((
y.at(iHit) - (*slope *
x.at(iHit) + *intercept)), 2.0);
907 *mse = *mse / (nHits - 2);
916 double RmsDeltarEtaStgc[12] = {1.43,1.53,1.53,1.56,1.59,1.54,1.70,1.69,1.76,1.81,1.83,1.84};
917 double RmsDeltarEtaMm[12] = {0.49,0.46,0.48,0.40,0.39,0.39,0.38,0.35,0.36,0.33,0.33,0.40};
919 stgcDeltaR[
bin] = RmsDeltarEtaStgc[
bin];
920 mmDeltaR[
bin] = RmsDeltarEtaMm[
bin];
930 std::vector<double>
r,
z;
931 std::vector<bool> isStgc;
932 for(
unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit) {
933 if (stgcHits.at(iHit).channelType == 1) {
934 r.push_back(stgcHits.at(iHit).r);
935 z.push_back(stgcHits.at(iHit).z);
936 isStgc.push_back(
true);
939 for(
unsigned int iHit = 0; iHit < mmHits.size(); ++iHit) {
940 if (mmHits.at(iHit).layerNumber < 2 || mmHits.at(iHit).layerNumber > 5) {
941 r.push_back(mmHits.at(iHit).r);
942 z.push_back(mmHits.at(iHit).z);
943 isStgc.push_back(
false);
944 side_mm = std::abs(mmHits.at(iHit).z)/mmHits.at(iHit).z;
947 double slopefit=0., interceptfit=99999., mse=-1.;
948 LinearFit(
z,
r,&slopefit,&interceptfit,&mse);
950 ATH_MSG_DEBUG(
"@@Merge@@ stgc_mmX_fit intercept= " << interceptfit);
953 std::vector<double> phiLocal;
954 double localPhiCenter = 3.*
M_PI;
955 for (
unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit) {
956 if (stgcHits.at(iHit).channelType != 2) {
continue;}
957 if (stgcHits.at(iHit).stationPhi<=5) {
958 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHit).stationPhi-1.);
960 localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHit).stationPhi-9.);
963 if (stgcHits.at(iHit).stationName == 57){
964 localPhiCenter +=
M_PI/8.;
965 if (stgcHits.at(iHit).stationPhi == 5) localPhiCenter -= 2. *
M_PI;
968 double phiProj = stgcHits.at(iHit).phi - localPhiCenter;
969 if (phiProj >
M_PI) phiProj -= 2.0*
M_PI;
970 if (phiProj < -1.*
M_PI) phiProj += 2.0*
M_PI;
971 double rInterpolate = slopefit * stgcHits.at(iHit).z + interceptfit;
972 double rProj = stgcHits.at(iHit).r;
974 ATH_MSG_DEBUG(
"@@Merge@@ philocalwire " << stgcHits.at(iHit).stationPhi <<
" " << stgcHits.at(iHit).stationName <<
" "
975 << localPhiCenter <<
" " << stgcHits.at(iHit).phi <<
" " << stgcHits.at(iHit).r <<
" " << stgcHits.at(iHit).z);
976 ATH_MSG_DEBUG(
"@@Merge@@ philocalwire " << rProj <<
" " << rInterpolate <<
" " <<
std::tan(phiProj) );
979 double tanTiltAngleU = 0,
981 double cosTiltAngleU = 0,
983 double sinTiltAngleU = 0,
986 tanTiltAngleU = TanM1p5,
987 tanTiltAngleV = TanP1p5;
988 cosTiltAngleU = CosM1p5,
989 cosTiltAngleV = CosP1p5;
990 sinTiltAngleU = SinM1p5,
991 sinTiltAngleV = SinP1p5;
993 tanTiltAngleU = TanP1p5,
994 tanTiltAngleV = TanM1p5;
995 cosTiltAngleU = CosP1p5,
996 cosTiltAngleV = CosM1p5;
997 sinTiltAngleU = SinP1p5,
998 sinTiltAngleV = SinM1p5;
1000 ATH_MSG_DEBUG(
"@@Merge@@ no U, V layer hits -> not consider tilt of U/V layers");
1002 for (
unsigned int iHit = 0; iHit < mmHits.size(); ++iHit) {
1003 if (localPhiCenter > 2.*
M_PI) {
1004 if (mmHits.at(iHit).stationPhi<=5) {
1005 localPhiCenter = 0.25 *
M_PI * ((
double)mmHits.at(iHit).stationPhi-1.);
1007 localPhiCenter = 0.25 *
M_PI * ((
double)mmHits.at(iHit).stationPhi-9.);
1009 if (mmHits.at(iHit).stationName == 55){
1010 localPhiCenter +=
M_PI/8.;
1011 if (mmHits.at(iHit).stationPhi == 5) localPhiCenter -= 2 *
M_PI;
1014 if (mmHits.at(iHit).layerNumber >1 && mmHits.at(iHit).layerNumber < 6){
1015 double rInterpolate = slopefit * mmHits.at(iHit).z + interceptfit;
1016 double rProj = mmHits.at(iHit).r;
1018 phiLocal.push_back(0);
1021 else if ((mmHits.at(iHit).layerNumber)%2 == 0) {
1022 phiLocal.push_back(
std::atan((rProj-rInterpolate)/tanTiltAngleU/rInterpolate));
1025 phiLocal.push_back(
std::atan((rProj-rInterpolate)/tanTiltAngleV/rInterpolate));
1030 double sumPhiLocal = 0;
1031 for (
unsigned int iHit = 0; iHit < phiLocal.size(); ++iHit ) {
1032 sumPhiLocal += phiLocal.at(iHit);
1034 double phiLocalAvg = sumPhiLocal/phiLocal.size();
1040 std::vector<double> r_stgc, z_stgc, r_mm, z_mm;
1041 std::vector<bool> isStgc_stgc, isStgc_mm;
1042 double side_stgc = 0;
1043 for(
unsigned int iHit = 0; iHit < stgcHits.size(); ++iHit) {
1044 if (stgcHits.at(iHit).stationPhi<=5) localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHit).stationPhi-1.);
1045 if (stgcHits.at(iHit).stationPhi> 5) localPhiCenter = 0.25 *
M_PI * ((
double)stgcHits.at(iHit).stationPhi-9.);
1046 if (stgcHits.at(iHit).channelType == 1) {
1047 r_stgc.push_back(stgcHits.at(iHit).r/
std::cos(phiLocalAvg));
1048 z_stgc.push_back(stgcHits.at(iHit).z);
1049 isStgc_stgc.push_back(
true);
1050 side_stgc = std::abs(stgcHits.at(iHit).z)/stgcHits.at(iHit).z;
1052 ATH_MSG_DEBUG(
"@@Merge@@ stgc strip_r " << phiLocalAvg <<
" " << stgcHits.at(iHit).z <<
" " << stgcHits.at(iHit).r/
std::cos(phiLocalAvg));
1056 double slopefit_stgc=0., interceptfit_stgc=99999., mse_stgc=1.e20;
1057 if(r_stgc.size() == 0) {
1060 LinearFit(z_stgc,r_stgc,&slopefit_stgc,&interceptfit_stgc,&mse_stgc);
1063 for(
unsigned int iHit = 0; iHit < mmHits.size(); ++iHit) {
1064 if (mmHits.at(iHit).layerNumber < 2 || mmHits.at(iHit).layerNumber > 5) {
1065 r_mm.push_back(mmHits.at(iHit).r/
std::cos(phiLocalAvg));
1066 z_mm.push_back(mmHits.at(iHit).z);
1067 isStgc_mm.push_back(
false);
1069 z_mm.push_back(mmHits.at(iHit).z);
1070 isStgc_mm.push_back(
false);
1072 double rProj = mmHits.at(iHit).r;
1074 r_mm.push_back(rProj);
1076 else if ((mmHits.at(iHit).layerNumber)%2 == 0) {
1077 double rPrime = (rProj * cosTiltAngleU)/(
cos(phiLocalAvg)*cosTiltAngleU +
sin(phiLocalAvg)*sinTiltAngleU);
1078 r_mm.push_back(rPrime);
1080 double rPrime = (rProj * cosTiltAngleV)/(
cos(phiLocalAvg)*cosTiltAngleV +
sin(phiLocalAvg)*sinTiltAngleV);
1081 r_mm.push_back(rPrime);
1085 double slopefit_mm=0., interceptfit_mm=99999., mse_mm=1.e20;
1086 if(r_mm.size() == 0) {
1089 LinearFit(z_mm,r_mm,&slopefit_mm,&interceptfit_mm,&mse_mm);
1092 unsigned int fmerge = 0;
1093 slopefit=0., interceptfit=99999., mse=1.e20;
1095 double StgcSegZ = 7526.329;
1096 double StgcSegR = 0;
1097 double MmSegZ = 7526.329;
1099 if (mse_stgc < 1.e7 && mse_mm < 1.e7) {
1101 copy(r_mm.begin(), r_mm.end(), back_inserter(
r));
1103 copy(z_mm.begin(), z_mm.end(), back_inserter(
z));
1104 isStgc = isStgc_stgc;
1105 copy(isStgc_mm.begin(), isStgc_mm.end(), back_inserter(isStgc));
1108 StgcSegZ = -7526.329;
1110 StgcSegR = slopefit_stgc * StgcSegZ + interceptfit_stgc;
1111 double StgcSegOriginTheta =
std::atan(StgcSegR / StgcSegZ);
1112 double StgcSegEta = side_stgc * (-
std::log(std::abs(
std::tan(StgcSegOriginTheta / 2))));
1116 MmSegR = slopefit_mm * MmSegZ + interceptfit_mm;
1117 double MmSegOriginTheta =
std::atan(MmSegR / MmSegZ);
1118 double MmSegEta = side_stgc * (-
std::log(std::abs(
std::tan(MmSegOriginTheta / 2))));
1120 double SegEtaAve = 0;
1122 if(side_stgc*side_mm > 0){
1124 SegEtaAve = (StgcSegEta + MmSegEta)/2;
1125 }
else if(std::abs(side_stgc) <
ZERO_LIMIT) {
1127 SegEtaAve = MmSegEta;
1130 SegEtaAve = StgcSegEta;
1134 LinearFitWeight(
z,
r,isStgc,&slopefit,&interceptfit,&mse,SegEtaAve);
1138 if (mse_stgc < mse_mm) {
1139 slopefit = slopefit_stgc;
1140 interceptfit = interceptfit_stgc;
1146 slopefit = slopefit_mm;
1147 interceptfit = interceptfit_mm;
1156 return StatusCode::SUCCESS;
1162 double NSWCenterZ = 7526.329;
1164 NSWCenterZ = -7526.329;
1166 superPoint->
R = slopefit * NSWCenterZ + interceptfit;
1167 superPoint->
Phim = phiLocalAvg+localPhiCenter;
1168 superPoint->
Z = NSWCenterZ;
1169 superPoint->
Npoint =
z.size();
1171 if (NSWCenterZ != 0) superPoint->
Alin = slopefit;
1172 superPoint->
Blin = interceptfit;
1174 ATH_MSG_DEBUG(
"Nsw Super Point r/phi/z/slope = "<<superPoint->
R<<
"/"<<superPoint->
Phim<<
"/"<<superPoint->
Z<<
"/"<<superPoint->
Alin);
1176 ATH_MSG_DEBUG(
"@@Merge@@ Nsw Super Point r/phi/z/slope = "<<superPoint->
R<<
"/"<<superPoint->
Phim<<
"/"<<superPoint->
Z<<
"/"<<superPoint->
Alin);
1177 ATH_MSG_DEBUG(
"@@Merge@@ fit slope= " << slopefit <<
" " << slopefit_stgc <<
" " << slopefit_mm);
1178 ATH_MSG_DEBUG(
"@@Merge@@ fit intercept= " << interceptfit <<
" " << interceptfit_stgc <<
" " << interceptfit_mm);
1179 ATH_MSG_DEBUG(
"@@Merge@@ fit mse= " << mse <<
" " << mse_stgc <<
" " << mse_mm);
1183 return StatusCode::SUCCESS;
1189 if(mmHits.size() == 0)
return StatusCode::SUCCESS;
1191 for(
unsigned int iHit = 0; iHit < mmHits.size(); iHit++){
1192 if(mmHits.at(iHit).isOutlier == 0){
1194 mmHits.at(iHit).isOutlier = 1;
1197 if(hitsInRoad == 0)
return StatusCode::SUCCESS;
1199 if(hitsInRoad < 6) {
1200 ATH_MSG_DEBUG(
"Number of MM hits is too small, at least 6 hits required : "<<hitsInRoad<<
" hits");
1201 return StatusCode::SUCCESS;
1202 }
else if(hitsInRoad > 100) {
1203 ATH_MSG_WARNING(
"Number of MM hits is too large, at most (2^16 - 1) hits allowed : "<<hitsInRoad<<
" hits");
1204 return StatusCode::SUCCESS;
1207 std::array< std::vector<int>, 8 > hitIdByLayer;
1208 for(
unsigned int iHit = 0; iHit < mmHits.size(); ++iHit){
1209 if(mmHits.at(iHit).isOutlier != 1)
continue;
1210 int layerNumber = mmHits.at(iHit).layerNumber;
1211 if (layerNumber > 7) {
1215 hitIdByLayer[layerNumber].push_back(iHit);
1218 <<
" " << hitIdByLayer[1].
size()
1219 <<
" " << hitIdByLayer[2].
size()
1220 <<
" " << hitIdByLayer[3].
size()
1221 <<
" " << hitIdByLayer[4].
size()
1222 <<
" " << hitIdByLayer[5].
size()
1223 <<
" " << hitIdByLayer[6].
size()
1224 <<
" " << hitIdByLayer[7].
size());
1226 std::vector< std::array<int, 8> > mmHitIds;
1227 findSetOfMmHitIds(mmHits, hitIdByLayer, mmHitIds);
1228 for (
unsigned int iHit = 0; iHit < mmHitIds.size(); ++iHit) {
1229 std::array<int, 8> hitIds = mmHitIds.at(iHit);
1230 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
1231 if (hitIds[iLayer] != -1) {
1232 mmHits.at(hitIds[iLayer]).isOutlier = 0;
1236 return StatusCode::SUCCESS;
1241 std::vector<std::array<int, 8>>& hitIdsCandidate)
const
1244 double NSWCenterZ = 7526.329;
1246 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
1247 if ( hitIdByLayer[iLayer].
size() > 0) {
1248 side = std::abs(mmHits.at(hitIdByLayer[iLayer].at(0)).z)/mmHits.at(hitIdByLayer[iLayer].at(0)).z;
1252 NSWCenterZ = NSWCenterZ *
side;
1254 std::array<std::vector<unsigned long int>,4> hitIdsInTwo;
1255 std::array<std::vector<double>,4> slopeInTwo;
1256 std::array<std::vector<double>,4> interceptInTwo;
1258 for(
unsigned int iPair = 0; iPair < 4; ++iPair){
1260 unsigned int nHitsInInner = hitIdByLayer[iPair].size();
1261 unsigned int nHitsInOuter;
1263 nHitsInOuter = hitIdByLayer[iPair+6].size();
1265 nHitsInOuter = hitIdByLayer[iPair+2].size();
1268 if ( nHitsInInner > 0xffff-1 || nHitsInOuter > 0xffff-1) {
1269 ATH_MSG_WARNING(
"Number of Mm hits in layers exceeds the limit of (2^16 - 1) : Number of Mm hits in "<<iPair<<
"th layer = "<< nHitsInInner
1270 <<
", Number of Mm hits in "<<iPair+4<<
"th layer = "<<nHitsInOuter);
1271 ATH_MSG_WARNING(
"Number of Mm hits is limitted to (2^16 - 1) and hits with id more than (2^16 -1) will be trancated.");
1272 if (nHitsInInner > 0xffff-1) {nHitsInInner = 0xffff-1;}
1273 if (nHitsInOuter > 0xffff-1) {nHitsInOuter = 0xffff-1;}
1276 bool foundCounterparts[0xffff] = {};
1278 for(
unsigned int iHit = 0; iHit < nHitsInInner; ++iHit){
1280 bool foundCounterpart = 0;
1285 int iHitId = hitIdByLayer[iPair].at(iHit);
1286 r[0] = mmHits.at(iHitId).r;
1287 z[0] = mmHits.at(iHitId).z;
1290 for(
unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
1294 jHitId = hitIdByLayer[iPair+6].at(jHit);
1296 jHitId = hitIdByLayer[iPair+2].at(jHit);
1298 r[1] = mmHits.at(jHitId).r;
1299 z[1] = mmHits.at(jHitId).z;
1301 double slope = (
r[1] -
r[0]) / (
z[1] -
z[0]);
1302 double intercept = slope*(0. -
z[0]) +
r[0];
1304 if(std::abs(slope) < 0.1 || std::abs(slope) > 0.7 || std::abs(intercept) > 500.)
continue;
1306 int encodedIds = (iHitId<<16) + jHitId;
1307 hitIdsInTwo[iPair].push_back(encodedIds);
1308 slopeInTwo[iPair].push_back(slope);
1309 interceptInTwo[iPair].push_back(intercept);
1311 foundCounterpart = 1;
1312 foundCounterparts[jHit] = 1;
1314 if(!foundCounterpart){
1315 int encodedIds = (iHitId<<16) + 0xffff;
1316 hitIdsInTwo[iPair].push_back(encodedIds);
1317 slopeInTwo[iPair].push_back(
r[0]/
z[0]);
1318 interceptInTwo[iPair].push_back(0.);
1322 for(
unsigned int jHit = 0; jHit < nHitsInOuter; ++jHit){
1323 if (!foundCounterparts[jHit]) {
1326 jHitId = hitIdByLayer[iPair+6].at(jHit);
1328 jHitId = hitIdByLayer[iPair+2].at(jHit);
1330 int encodedIds = (0xFFFF
u<<16) + jHitId;
1331 hitIdsInTwo[iPair].push_back(encodedIds);
1332 slopeInTwo[iPair].push_back(mmHits.at(jHitId).r/mmHits.at(jHitId).z);
1333 interceptInTwo[iPair].push_back(0.);
1337 ATH_MSG_DEBUG(
"@@MM@@ Npairs " << hitIdsInTwo[0].
size() <<
" " << hitIdsInTwo[1].
size() <<
" " << hitIdsInTwo[2].
size() <<
" " << hitIdsInTwo[3].
size());
1338 for (
unsigned int iLayer = 0; iLayer < 4; ++ iLayer) {
1339 for (
unsigned int iPair = 0; iPair < slopeInTwo[iLayer].size(); ++iPair) {
1340 ATH_MSG_DEBUG(
"@@MM@@ pair fit slope= " << slopeInTwo[iLayer].at(iPair) <<
" intercept= " << interceptInTwo[iLayer].at(iPair));
1344 std::vector<std::array<int, 4>> hitIdsInFourX;
1345 std::vector<double> slopeInFourX;
1346 std::vector<double> interceptInFourX;
1347 std::vector<double> mseInFourX;
1349 unsigned int nPairsInInnerX = hitIdsInTwo[0].size();
1350 unsigned int nPairsInOuterX = hitIdsInTwo[1].size();
1352 for(
unsigned int iPairX = 0; iPairX < nPairsInInnerX; ++iPairX){
1355 double intercept[2];
1358 slope[0] = slopeInTwo[0].at(iPairX);
1359 intercept[0] = interceptInTwo[0].at(iPairX);
1360 spR[0] = slope[0] * NSWCenterZ + intercept[0];
1361 for(
unsigned int jPairX = 0; jPairX < nPairsInOuterX; ++jPairX){
1362 int ihitIds = hitIdsInTwo[0].at(iPairX);
1363 int jhitIds = hitIdsInTwo[1].at(jPairX);
1364 if ( ((ihitIds>>16 & 0xffff) == 0xffff || (ihitIds & 0xffff) == 0xffff) &&
1365 ((jhitIds>>16 & 0xffff) == 0xffff || (jhitIds & 0xffff) == 0xffff ))
continue;
1367 slope[1] = slopeInTwo[1].at(jPairX);
1368 intercept[1] = interceptInTwo[1].at(jPairX);
1369 spR[1] = slope[1] * NSWCenterZ + intercept[1];
1371 if(std::abs(spR[1] - spR[0]) > 50. ||
1372 std::abs((intercept[1] + intercept[0]) / 2) > 200.)
continue;
1374 std::array<int, 4> setOfHitIds{};
1375 setOfHitIds[0] = (ihitIds>>16 & 0xffff);
1376 setOfHitIds[1] = (ihitIds & 0xffff);
1377 setOfHitIds[2] = (jhitIds>>16 & 0xffff);
1378 setOfHitIds[3] = (jhitIds & 0xffff);
1379 std::vector<double>
r;
1380 std::vector<double>
z;
1381 for(
unsigned int iLayer = 0; iLayer < 4; ++iLayer){
1382 if(setOfHitIds[iLayer] == 0xffff) {
1385 double rhit = mmHits.at(setOfHitIds[iLayer]).r;
1386 double zhit = mmHits.at(setOfHitIds[iLayer]).z;
1390 double slopefit=0., interceptfit=99999., mse=-1.;
1391 LinearFit(
z,
r,&slopefit, &interceptfit, &mse);
1393 hitIdsInFourX.push_back(setOfHitIds);
1394 slopeInFourX.push_back(slopefit);
1395 interceptInFourX.push_back(interceptfit);
1396 mseInFourX.push_back(mse);
1401 for (
unsigned int iQuad = 0; iQuad < slopeInFourX.size(); ++iQuad) {
1402 ATH_MSG_DEBUG(
"@@MM@@ X quad fit slope= " << slopeInFourX.at(iQuad) <<
" intercept= " << interceptInFourX.at(iQuad) <<
" mse= " << mseInFourX.at(iQuad));
1405 if(!hitIdsInFourX.size()){
1410 double tanTiltAngleU = 0,
1413 tanTiltAngleU =
tan( 1.5/360.*2.*
M_PI),
1414 tanTiltAngleV =
tan(-1.5/360.*2.*
M_PI);
1416 tanTiltAngleU =
tan(-1.5/360.*2.*
M_PI),
1417 tanTiltAngleV =
tan(1.5/360.*2.*
M_PI);
1420 std::vector< std::array<int, 8> > hitIdsInEight;
1421 std::vector<double> mseInEight;
1423 for(
unsigned int iQuadX = 0; iQuadX < hitIdsInFourX.size(); ++iQuadX){
1424 if(mseInFourX.at(iQuadX) > 10)
continue;
1426 double slopeX = slopeInFourX.at(iQuadX);
1427 double interceptX = interceptInFourX.at(iQuadX);
1428 std::array<int,4> hitIdsX{};
1429 hitIdsX = hitIdsInFourX.at(iQuadX);
1431 for (
unsigned int iPairU = 0; iPairU < hitIdsInTwo[2].size(); ++iPairU) {
1434 hitIdsU[0] = hitIdsInTwo[2].at(iPairU)>>16 & 0xffff;
1435 hitIdsU[1] = hitIdsInTwo[2].at(iPairU) & 0xffff;
1436 double phiLocalU[2]={-99999,-99999};
1437 for(
unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1438 if (hitIdsU[iLayer] == 0xffff)
continue;
1439 if (hitIdsU[iLayer] < 0) {
1440 ATH_MSG_DEBUG(
"@@MM@@ hitIdsU[iLayer] iLayer= " << iLayer <<
" hitIdsU[iLayer]= " << hitIdsU[iLayer]);
1442 double rInterpolate = slopeX * mmHits.at(hitIdsU[iLayer]).z + interceptX;
1443 double rProj = mmHits.at(hitIdsU[iLayer]).r;
1445 phiLocalU[iLayer] = 0;
1447 phiLocalU[iLayer] =
std::atan((rProj-rInterpolate)/tanTiltAngleU/rInterpolate);
1450 for(
unsigned int iPairV = 0; iPairV < hitIdsInTwo[3].size(); ++iPairV) {
1453 hitIdsV[0] = hitIdsInTwo[3].at(iPairV)>>16 & 0xffff;
1454 hitIdsV[1] = hitIdsInTwo[3].at(iPairV) & 0xffff;
1456 if( (hitIdsU[0] == 0xffff || hitIdsU[1] == 0xffff) &&
1457 (hitIdsV[0] == 0xffff || hitIdsV[1] == 0xffff) )
continue;
1459 double phiLocalV[2]={-99999,-99999};
1460 for(
unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1461 if (hitIdsV[iLayer] == 0xffff)
continue;
1462 if (hitIdsV[iLayer] < 0) {
1463 ATH_MSG_DEBUG(
"@@MM@@ hitIdsV[iLayer] iLayer= " << iLayer <<
" hitIdsV[iLayer]= " << hitIdsV[iLayer]);
1465 double rInterpolate = slopeX * mmHits.at(hitIdsV[iLayer]).z + interceptX;
1466 double rProj = mmHits.at(hitIdsV[iLayer]).r;
1468 phiLocalV[iLayer] = 0;
1470 phiLocalV[iLayer] =
std::atan((rProj-rInterpolate)/tanTiltAngleV/rInterpolate);
1473 if ( std::abs(phiLocalU[0]-phiLocalV[0]) > 0.05 &&
1474 std::abs(phiLocalU[1]-phiLocalV[1]) > 0.05)
continue;
1477 double phiLocalUV = 0;
1479 for(
unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1480 if(phiLocalU[iLayer] > -99999.) {
1481 phiLocalUV += phiLocalU[iLayer];
1484 if(phiLocalV[iLayer] > -99999.) {
1485 phiLocalUV += phiLocalV[iLayer];
1491 std::array<int, 8> setOfHitIds = {-1,-1,-1,-1,-1,-1,-1,-1};
1492 std::vector<double>
r,
z;
1493 for (
unsigned int iLayer = 0; iLayer < 4; ++iLayer) {
1494 if ( hitIdsX[iLayer] != 0xffff) {
1495 if (hitIdsX[iLayer] < 0) {
1496 ATH_MSG_DEBUG(
"@@MM@@ hitIdsX[iLayer] iLayer= " << iLayer <<
" hitIdsX[iLayer]= " << hitIdsX[iLayer]);
1498 z.push_back(mmHits.at(hitIdsX[iLayer]).z);
1499 r.push_back(mmHits.at(hitIdsX[iLayer]).r /
std::cos(phiLocalUV));
1500 setOfHitIds[iLayer] = hitIdsX[iLayer];
1503 for (
unsigned int iLayer = 0; iLayer < 2; ++iLayer) {
1504 if ( hitIdsU[iLayer] != 0xffff) {
1505 if (hitIdsU[iLayer] < 0) {
1506 ATH_MSG_DEBUG(
"@@MM@@ 2 hitIdsU[iLayer] iLayer= " << iLayer <<
" hitIdsU[iLayer]= " << hitIdsU[iLayer]);
1508 z.push_back(mmHits.at(hitIdsU[iLayer]).z);
1509 r.push_back(mmHits.at(hitIdsU[iLayer]).r*(
std::cos(phiLocalUV) + 1/
std::cos(phiLocalUV))/2.);
1510 setOfHitIds[iLayer+4] = hitIdsU[iLayer];
1512 if ( hitIdsV[iLayer] != 0xffff) {
1513 if (hitIdsV[iLayer] < 0) {
1514 ATH_MSG_DEBUG(
"@@MM@@ 2 hitIdsV[iLayer] iLayer= " << iLayer <<
" hitIdsV[iLayer]= " << hitIdsV[iLayer]);
1516 z.push_back(mmHits.at(hitIdsV[iLayer]).z);
1517 r.push_back(mmHits.at(hitIdsV[iLayer]).r*(
std::cos(phiLocalUV) + 1/
std::cos(phiLocalUV))/2.);
1518 setOfHitIds[iLayer+6] = hitIdsV[iLayer];
1521 double slopefit=0., interceptfit=99999., mse=-1.;
1522 LinearFit(
z,
r,&slopefit,&interceptfit,&mse);
1524 hitIdsInEight.push_back(setOfHitIds);
1525 mseInEight.push_back(mse);
1531 std::vector<int> nOctetSegments;
1532 std::vector<int> patternStationName;
1533 for (
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet) {
1534 bool isFirstHit =
true;
1535 int hitStationName = 0;
1536 int nOctetSegment = 0;
1537 ATH_MSG_DEBUG(
"@@MM@@ octet fit mse " << mseInEight.at(iOctet));
1538 std::array<int, 8> tmpOctet = hitIdsInEight.at(iOctet);
1539 for (
unsigned int iLayer = 0; iLayer < 8; ++iLayer) {
1540 if (tmpOctet[iLayer] != -1) {
1543 hitStationName = mmHits.at(tmpOctet[iLayer]).stationName;
1546 ATH_MSG_DEBUG(
"@@MM@@ octet pos r= " << mmHits.at(tmpOctet[iLayer]).r <<
" phi= " << mmHits.at(tmpOctet[iLayer]).phi <<
" z= " << mmHits.at(tmpOctet[iLayer]).z);
1549 else if(mmHits.at(tmpOctet[iLayer]).stationName == hitStationName){
1550 ATH_MSG_DEBUG(
"@@MM@@ octet pos r= " << mmHits.at(tmpOctet[iLayer]).r <<
" phi= " << mmHits.at(tmpOctet[iLayer]).phi <<
" z= " << mmHits.at(tmpOctet[iLayer]).z);
1556 nOctetSegments.push_back(nOctetSegment);
1557 patternStationName.push_back(hitStationName);
1560 double mseminL = 100000.;
1561 double mseminS = 100000.;
1562 std::vector<int> octetIds(2,-1);
1563 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
1564 if(patternStationName.at(iOctet) == 56){
1565 if( mseInEight.at(iOctet) < mseminL) {
1566 mseminL = mseInEight.at(iOctet);
1569 else if(patternStationName.at(iOctet) == 55){
1570 if( mseInEight.at(iOctet) < mseminS) {
1571 mseminS = mseInEight.at(iOctet);
1576 for(
unsigned int iOctet = 0; iOctet < hitIdsInEight.size(); ++iOctet){
1577 if(patternStationName.at(iOctet) == 56){
1578 if( mseInEight.at(iOctet) != mseminL) {
1582 else if(patternStationName.at(iOctet) == 55){
1583 if( mseInEight.at(iOctet) != mseminS) {
1587 octetIds.push_back(iOctet);
1590 for(
unsigned int ids = 0;
ids < octetIds.size();
ids++){
1591 if (octetIds.at(
ids) != -1) {
1592 hitIdsCandidate.push_back(hitIdsInEight.at(octetIds.at(
ids)));