27 #include "gsl/gsl_statistics.h"
37 isPatternChecked =
false;
39 MagnecticFieldThreshold = 0.5;
49 deltaRThreshold = iConfig.
getParameter<
double>(
"deltaRThreshold");
51 autoAlgorithmChoose = iConfig.
getParameter<
bool>(
"autoAlgorithmChoose");
53 MinDeltaPhi = iConfig.
getParameter<
double>(
"MinDeltaPhi");
54 MagnecticFieldThreshold = iConfig.
getParameter<
double>(
"MagnecticFieldThreshold");
56 sampleCount = iConfig.
getParameter<
unsigned int>(
"sampleCount");
62 if(isConfigured ==
false)
64 cout <<
"Configuration not set yet" << endl;
65 return createFakeSeed(isGoodSeed, eSetup);
69 unsigned int NumberofHitsinSeed = nrhit();
70 if(NumberofHitsinSeed < 3)
71 return createFakeSeed(isGoodSeed, eSetup);
73 if(NumberofHitsinSeed == 3)
74 ThreePointsAlgorithm();
76 if(NumberofHitsinSeed > 3)
78 if(autoAlgorithmChoose ==
false)
80 cout <<
"computePtWithmorerecHits" << endl;
82 ThreePointsAlgorithm();
84 MiddlePointsAlgorithm();
90 SegmentAlgorithmSpecial(eSetup);
93 cout <<
"Not enough recHits for Special Segment Algorithm" << endl;
94 return createFakeSeed(isGoodSeed, eSetup);
103 SegmentAlgorithmSpecial(eSetup);
114 if(isPatternChecked ==
false){
116 checkSimplePattern(eSetup);
118 checkSegmentAlgorithmSpecial(eSetup);
122 return createSeed(isGoodSeed, eSetup);
127 cout <<
"computePtWith3recHits" << endl;
128 unsigned int NumberofHitsinSeed = nrhit();
130 if(NumberofHitsinSeed < 3)
132 isPatternChecked =
true;
137 unsigned int NumberofPart = NumberofHitsinSeed * (NumberofHitsinSeed - 1) * (NumberofHitsinSeed - 2) / (3 * 2);;
138 double *pt =
new double[NumberofPart];
139 double *pt_err =
new double[NumberofPart];
143 unsigned int NumberofStraight = 0;
144 for(
unsigned int i = 0;
i < (NumberofHitsinSeed - 2);
i++)
145 for(
unsigned int j = (
i + 1);
j < (NumberofHitsinSeed - 1);
j++)
146 for(
unsigned int k = (
j + 1);
k < NumberofHitsinSeed;
k++)
148 precHit[0] = theRecHits[
i];
149 precHit[1] = theRecHits[
j];
150 precHit[2] = theRecHits[
k];
151 bool checkStraight = checkStraightwithThreerecHits(precHit, MinDeltaPhi);
154 GlobalVector Center_temp = computePtwithThreerecHits(pt[n], pt_err[n], precHit);
156 Center += Center_temp;
168 if(NumberofStraight == NumberofPart)
176 Center /= (NumberofPart - NumberofStraight);
178 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
179 meanR += getDistance(*iter, Center);
180 meanR /= NumberofHitsinSeed;
185 isPatternChecked =
false;
197 cout <<
"Using middle points algorithm" << endl;
198 unsigned int NumberofHitsinSeed = nrhit();
200 if(NumberofHitsinSeed < 4)
202 isPatternChecked =
true;
206 double *
X =
new double[NumberofHitsinSeed];
207 double *Y =
new double[NumberofHitsinSeed];
209 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
211 X[
n] = (*iter)->globalPosition().x();
212 Y[
n] = (*iter)->globalPosition().y();
213 cout <<
"X[" << n <<
"] = " << X[
n] <<
", Y[" << n <<
"]= " << Y[
n] << endl;
216 unsigned int NumberofPoints = NumberofHitsinSeed;
217 while(NumberofPoints > 3)
219 for(
unsigned int i = 0;
i <= (NumberofPoints - 2);
i++)
221 X[
i] = (X[
i] + X[
i+1]) / 2;
222 Y[
i] = (Y[
i] + Y[
i+1]) / 2;
227 for(
unsigned int i = 0;
i < 3;
i++)
234 bool checkStraight = checkStraightwithThreerecHits(x, y, MinDeltaPhi);
238 GlobalVector Center_temp = computePtWithThreerecHits(pt, pt_err, x, y);
240 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
241 meanR += getDistance(*iter, Center_temp);
242 meanR /= NumberofHitsinSeed;
245 Center = Center_temp;
256 isPatternChecked =
false;
264 cout <<
"Using segments algorithm" << endl;
265 unsigned int NumberofHitsinSeed = nrhit();
267 if(NumberofHitsinSeed < 4)
269 isPatternChecked =
true;
275 unsigned int NumberofSegment = NumberofHitsinSeed - 2;
278 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-2); iter++)
280 Segment[
n].first = (*iter);
281 Segment[
n].second = (*(iter + 2));
284 unsigned int NumberofStraight = 0;
285 for(
unsigned int i = 0;
i < NumberofSegment - 1;
i++)
287 bool checkStraight = checkStraightwithSegment(Segment[
i], Segment[i+1], MinDeltaPhi);
288 if(checkStraight ==
true)
295 GlobalVector Center_temp = computePtwithSegment(Segment[i], Segment[i+1]);
297 Center += Center_temp;
301 if((NumberofSegment-1-NumberofStraight) > 0)
304 Center /= (NumberofSegment - 1 - NumberofStraight);
306 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
307 meanR += getDistance(*iter, Center);
308 meanR /= NumberofHitsinSeed;
318 isPatternChecked =
false;
332 isPatternChecked =
true;
338 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
341 GlobalPoint gpLast = (*(iter+1))->globalPosition();
343 double dx = (gpLast.
x() - gpFirst.
x()) / (sampleCount + 1);
344 double dy = (gpLast.
y() - gpFirst.
y()) / (sampleCount + 1);
345 double dz = (gpLast.
z() - gpFirst.
z()) / (sampleCount + 1);
350 cout <<
"Sampling magnetic field : " << MagneticVec_temp << endl;
357 ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin();
358 for(
unsigned int n = 0;
n <= 1;
n++)
360 SegmentRB[
n].first = (*iter);
361 cout <<
"SegmentRB " <<
n <<
" recHit: " << (*iter)->globalPosition() << endl;
363 SegmentRB[
n].second = (*iter);
364 cout <<
"SegmentRB " <<
n <<
" recHit: " << (*iter)->globalPosition() << endl;
367 GlobalVector segvec1 = (SegmentRB[0].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
368 GlobalVector segvec2 = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
371 entryPosition = (SegmentRB[0].second)->globalPosition();
372 leavePosition = (SegmentRB[1].first)->globalPosition();
373 while(fabs(Field->inTesla(entryPosition).z()) < MagnecticFieldThreshold)
375 cout <<
"Entry position is : " << entryPosition <<
", and stepping into next point" << endl;
376 entryPosition += segvec1.
unit() * stepLength;
379 while(fabs(Field->inTesla(entryPosition).z()) >= MagnecticFieldThreshold)
381 cout <<
"Entry position is : " << entryPosition <<
", and stepping back into next point" << endl;
382 entryPosition -= segvec1.
unit() * stepLength / 10;
384 entryPosition += 0.5 * segvec1.
unit() * stepLength / 10;
385 cout <<
"Final entry position is : " << entryPosition << endl;
387 while(fabs(Field->inTesla(leavePosition).z()) < MagnecticFieldThreshold)
389 cout <<
"Leave position is : " << leavePosition <<
", and stepping into next point" << endl;
390 leavePosition -= segvec2.
unit() * stepLength;
393 while(fabs(Field->inTesla(leavePosition).z()) >= MagnecticFieldThreshold)
395 cout <<
"Leave position is : " << leavePosition <<
", and stepping back into next point" << endl;
396 leavePosition += segvec2.
unit() * stepLength / 10;
398 leavePosition -= 0.5 * segvec2.
unit() * stepLength / 10;
399 cout <<
"Final leave position is : " << leavePosition << endl;
403 double dx = (leavePosition.x() - entryPosition.x()) / (sampleCount + 1);
404 double dy = (leavePosition.y() - entryPosition.y()) / (sampleCount + 1);
405 double dz = (leavePosition.z() - entryPosition.z()) / (sampleCount + 1);
406 std::vector<GlobalVector>
BValue;
412 cout <<
"Sampling magnetic field : " << MagneticVec_temp << endl;
413 BValue.push_back(MagneticVec_temp);
417 for(std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
419 meanB2 /= BValue.size();
420 cout <<
"Mean B field is " << meanB2 << endl;
421 meanMagneticField2 = meanB2;
423 double meanBz2 = meanB2.
z();
424 double deltaBz2 = 0.;
425 for(std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
426 deltaBz2 += (BIter->z() - meanBz2) * (BIter->z() - meanBz2);;
427 deltaBz2 /= BValue.size();
428 deltaBz2 =
sqrt(deltaBz2);
429 cout<<
"delta Bz is " << deltaBz2 << endl;
433 bool checkStraight = checkStraightwithSegment(SegmentRB[0], SegmentRB[1], MinDeltaPhi);
434 if(checkStraight ==
true)
437 isStraight2 = checkStraight;
440 GlobalVector MomentumVec = (SegmentRB[1].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
441 S += MomentumVec.
perp();
442 lastPhi = MomentumVec.
phi().value();
446 GlobalVector seg1 = entryPosition - (SegmentRB[0].first)->globalPosition();
448 GlobalVector seg2 = (SegmentRB[1].second)->globalPosition() - leavePosition;
453 double A1 = gvec1.
x();
454 double B1 = gvec1.
y();
455 double A2 = gvec2.
x();
456 double B2 = gvec2.
y();
457 double X1 = entryPosition.x();
458 double Y1 = entryPosition.y();
459 double X2 = leavePosition.x();
460 double Y2 = leavePosition.y();
461 double XO = (A1*A2*(Y2-Y1)+A2*B1*X1-A1*B2*X2)/(A2*B1-A1*B2);
462 double YO = (B1*B2*(X2-X1)+B2*A1*Y1-B1*A2*Y2)/(B2*A1-B1*A2);
465 isStraight2 = checkStraight;
466 Center2 = Center_temp;
468 cout <<
"entryPosition: " << entryPosition << endl;
469 cout <<
"leavePosition: " << leavePosition << endl;
470 cout <<
"Center2 is : " << Center_temp << endl;
472 double R1 =
GlobalVector((entryPosition.x() - Center_temp.
x()), (entryPosition.y() - Center_temp.
y()), (entryPosition.z() - Center_temp.
z())).
perp();
473 double R2 =
GlobalVector((leavePosition.x() - Center_temp.
x()), (leavePosition.y() - Center_temp.
y()), (leavePosition.z() - Center_temp.
z())).
perp();
474 double meanR = (R1 + R2) / 2;
475 double deltaR =
sqrt(((R1-meanR)*(R1-meanR)+(R2-meanR)*(R2-meanR))/2);
477 cout <<
"R1 is " << R1 <<
", R2 is " << R2 << endl;
478 cout <<
"Mean radius is " << meanR << endl;
479 cout <<
"Delta R is " << deltaR << endl;
482 lastPhi = seg2.
phi().value();
486 isPatternChecked =
false;
492 unsigned int count = 0;
494 for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
498 if(dynamic_cast<const RPCChamber*>(Detector) != 0)
503 unsigned int Station = RPCId.
station();
517 cout <<
"Check for segment fit: " << isFit << endl;
523 return theRecHits.front();
532 for (ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
546 return sqrt((precHit->globalPosition().x()-Center.
x())*(precHit->globalPosition().x()-Center.
x())+(precHit->globalPosition().y()-Center.
y())*(precHit->globalPosition().y()-Center.
y()));
551 GlobalVector segvec1 = precHit[1]->globalPosition() - precHit[0]->globalPosition();
552 GlobalVector segvec2 = precHit[2]->globalPosition() - precHit[1]->globalPosition();
554 if(fabs(dPhi) > MinDeltaPhi)
556 cout <<
"Part is estimate to be not straight" << endl;
561 cout <<
"Part is estimate to be straight" << endl;
569 x[0] = precHit[0]->globalPosition().x();
570 y[0] = precHit[0]->globalPosition().y();
571 x[1] = precHit[1]->globalPosition().x();
572 y[1] = precHit[1]->globalPosition().y();
573 x[2] = precHit[2]->globalPosition().x();
574 y[2] = precHit[2]->globalPosition().y();
575 double A = (y[2]-y[1])/(x[2]-x[1]) - (y[1]-y[0])/(x[1]-x[0]);
576 double TYO = (x[2]-x[0])/A + (y[2]*y[2]-y[1]*y[1])/((x[2]-x[1])*A) - (y[1]*y[1]-y[0]*y[0])/((x[1]-x[0])*
A);
577 double TXO = (x[2]+x[1]) + (y[2]*y[2]-y[1]*y[1])/(x[2]-x[1]) - TYO*(y[2]-y[1])/(x[2]-x[1]);
578 double XO = 0.5 * TXO;
579 double YO = 0.5 * TYO;
580 double R2 = (x[0]-XO)*(x[0]-XO) + (y[0]-YO)*(y[0]-YO);
581 cout <<
"R2 is " << R2 << endl;
583 pt = 0.01 *
sqrt(R2) * 2 * 0.3;
584 cout <<
"pt is " << pt << endl;
591 GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
592 GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
593 GlobalVector segvec3 = (Segment2.first)->globalPosition() - (Segment1.first)->globalPosition();
595 double dPhi1 = (segvec2.
phi() - segvec1.
phi()).
value();
596 double dPhi2 = (segvec3.
phi() - segvec1.
phi()).
value();
597 cout <<
"Checking straight with 2 segments. dPhi1: " << dPhi1 <<
", dPhi2: " << dPhi2 << endl;
598 cout <<
"Checking straight with 2 segments. dPhi1 in degree: " << dPhi1*180/3.1415926 <<
", dPhi2 in degree: " << dPhi2*180/3.1415926 << endl;
599 if(fabs(dPhi1) > MinDeltaPhi || fabs(dPhi2) > MinDeltaPhi)
601 cout <<
"Segment is estimate to be not straight" << endl;
606 cout <<
"Segment is estimate to be straight" << endl;
613 GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
614 GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
615 GlobalPoint Point1(((Segment1.second)->globalPosition().x() + (Segment1.first)->globalPosition().x()) / 2, ((Segment1.second)->globalPosition().y() + (Segment1.first)->globalPosition().y()) / 2, ((Segment1.second)->globalPosition().z() + (Segment1.first)->globalPosition().z()) / 2);
616 GlobalPoint Point2(((Segment2.second)->globalPosition().x() + (Segment2.first)->globalPosition().x()) / 2, ((Segment2.second)->globalPosition().y() + (Segment2.first)->globalPosition().y()) / 2, ((Segment2.second)->globalPosition().z() + (Segment2.first)->globalPosition().z()) / 2);
620 double A1 = gvec1.
x();
621 double B1 = gvec1.
y();
622 double A2 = gvec2.
x();
623 double B2 = gvec2.
y();
624 double X1 = Point1.
x();
625 double Y1 = Point1.
y();
626 double X2 = Point2.
x();
627 double Y2 = Point2.
y();
628 double XO = (A1*A2*(Y2-Y1)+A2*B1*X1-A1*B2*X2)/(A2*B1-A1*B2);
629 double YO = (B1*B2*(X2-X1)+B2*A1*Y1-B1*A2*Y2)/(B2*A1-B1*A2);
639 if(fabs(dPhi) > MinDeltaPhi)
641 cout <<
"Part is estimate to be not straight" << endl;
646 cout <<
"Part is estimate to be straight" << endl;
653 double A = (
y[2]-
y[1])/(
x[2]-
x[1]) - (
y[1]-
y[0])/(
x[1]-
x[0]);
654 double TYO = (
x[2]-
x[0])/A + (
y[2]*
y[2]-
y[1]*
y[1])/((
x[2]-
x[1])*A) - (y[1]*y[1]-y[0]*y[0])/((
x[1]-
x[0])*
A);
655 double TXO = (
x[2]+
x[1]) + (y[2]*y[2]-y[1]*y[1])/(
x[2]-
x[1]) - TYO*(y[2]-y[1])/(
x[2]-
x[1]);
656 double XO = 0.5 * TXO;
657 double YO = 0.5 * TYO;
658 double R2 = (
x[0]-XO)*(
x[0]-XO) + (y[0]-YO)*(y[0]-YO);
659 cout <<
"R2 is " << R2 << endl;
661 pt = 0.01 *
sqrt(R2) * 2 * 0.3;
662 cout <<
"pt is " << pt << endl;
669 if(isPatternChecked ==
true)
676 unsigned int NumberofHitsinSeed = nrhit();
679 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
680 cout <<
"Position of recHit is: " << (*iter)->globalPosition() << endl;
683 std::vector<double> BzValue;
684 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
687 GlobalPoint gpLast = (*(iter+1))->globalPosition();
689 double dx = (gpLast.
x() - gpFirst.
x()) / (sampleCount + 1);
690 double dy = (gpLast.
y() - gpFirst.
y()) / (sampleCount + 1);
691 double dz = (gpLast.
z() - gpFirst.
z()) / (sampleCount + 1);
696 cout <<
"Sampling magnetic field : " << MagneticVec_temp << endl;
697 BzValue.push_back(MagneticVec_temp.z());
703 meanBz += BzValue[
index];
704 meanBz /= BzValue.size();
705 cout <<
"Mean Bz is " << meanBz << endl;
708 deltaBz += (BzValue[
index] - meanBz) * (BzValue[
index] - meanBz);
709 deltaBz /= BzValue.size();
710 deltaBz =
sqrt(deltaBz);
711 cout<<
"delata Bz is " << deltaBz << endl;
717 if(fabs((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
719 if(((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
727 cout <<
" Check isParralZ is :" << isParralZ << endl;
728 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
732 if(fabs((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z()) > ZError)
734 cout <<
"Pattern find error in Z direction: wrong perpendicular direction" << endl;
740 if((
int)(((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z())/ZError)*isParralZ < 0)
742 cout <<
"Pattern find error in Z direction: wrong Z direction" << endl;
749 if(isStraight ==
false)
753 unsigned int index = 0;
754 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
756 GlobalVector vec_temp(((*iter)->globalPosition().x()-Center.x()), ((*iter)->globalPosition().y()-Center.y()), ((*iter)->globalPosition().z()-Center.z()));
757 vec[
index] = vec_temp;
761 for(
unsigned int index = 0; index < (NumberofHitsinSeed-1); index++)
764 if((vec[index+1].
phi()-vec[index].
phi()) > 0)
768 cout <<
"Current isClockwise is : " << isClockwise << endl;
770 cout <<
"Check isClockwise is : " << isClockwise << endl;
771 if((
unsigned int)
abs(isClockwise) != (NumberofHitsinSeed-1))
773 cout <<
"Pattern find error in Phi direction" << endl;
778 isClockwise /=
abs(isClockwise);
782 double deltaRwithBz = fabs(deltaBz * meanRadius / meanBz);
783 cout <<
"deltaR with Bz is " << deltaRwithBz << endl;
788 meanPt = 0.01 * meanRadius * meanBz * 0.3 * isClockwise;
791 cout <<
" meanRadius is " << meanRadius << endl;
792 cout <<
" meanPt is " << meanPt << endl;
795 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
797 deltaR += (getDistance(*iter, Center) - meanRadius) * (getDistance(*iter, Center) - meanRadius);
799 deltaR = deltaR / NumberofHitsinSeed;
800 deltaR =
sqrt(deltaR);
803 cout <<
"DeltaR is " << deltaR << endl;
804 if(deltaR > deltaRThreshold)
806 cout <<
"Pattern find error: deltaR over threshold" << endl;
816 meanSpt = deltaRThreshold;
818 cout <<
"III--> Seed Pt : " << meanPt << endl;
819 cout <<
"III--> Pattern is: " << isGoodPattern << endl;
822 isPatternChecked =
true;
827 if(isPatternChecked ==
true)
832 isPatternChecked =
true;
841 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
842 cout <<
"Position of recHit is: " << (*iter)->globalPosition() << endl;
845 if(fabs((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
847 if(((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
855 cout <<
" Check isParralZ is :" << isParralZ << endl;
856 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
860 if(fabs((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z()) > ZError)
862 cout <<
"Pattern find error in Z direction: wrong perpendicular direction" << endl;
868 if((
int)(((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z())/ZError)*isParralZ < 0)
870 cout <<
"Pattern find error in Z direction: wrong Z direction" << endl;
877 if(isStraight2 ==
true)
883 meanSpt = deltaRThreshold;
886 GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
887 GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
889 unsigned int index = 0;
890 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
897 double tracklength = 0;
898 cout <<
"Now checking recHit " << index << endl;
899 double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, eSetup);
900 cout <<
"Final distance is " << Distance << endl;
901 if(Distance > MaxRSD)
903 cout <<
"Pattern find error in distance for other recHits: " << Distance << endl;
913 vec[0] =
GlobalVector((entryPosition.x()-Center2.x()), (entryPosition.y()-Center2.y()), (entryPosition.z()-Center2.z()));
914 vec[1] =
GlobalVector((leavePosition.x()-Center2.x()), (leavePosition.y()-Center2.y()), (leavePosition.z()-Center2.z()));
921 cout <<
"Check isClockwise is : " << isClockwise << endl;
924 meanPt = 0.01 * meanRadius2 * meanMagneticField2.z() * 0.3 * isClockwise;
926 cout <<
" meanRadius is " << meanRadius2 <<
", with meanBz " << meanMagneticField2.z() << endl;
927 cout <<
" meanPt is " << meanPt << endl;
932 cout <<
"entryPosition: " << entryPosition << endl;
933 cout <<
"leavePosition: " << leavePosition << endl;
934 cout <<
"Center2 is : " << Center2 << endl;
935 double R1 = vec[0].
perp();
936 double R2 = vec[1].
perp();
937 double deltaR =
sqrt(((R1-meanRadius2)*(R1-meanRadius2)+(R2-meanRadius2)*(R2-meanRadius2))/2);
939 cout <<
"R1 is " << R1 <<
", R2 is " << R2 << endl;
940 cout <<
"Delta R for the initial 3 segments is " << deltaR << endl;
941 if(deltaR > deltaRThreshold)
943 cout <<
"Pattern find error in delta R for the initial 3 segments" << endl;
948 GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
949 GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
950 GlobalVector startMomentum = startSegment*(fabs(meanPt)/startSegment.
perp());
951 unsigned int index = 0;
952 for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
959 double tracklength = 0;
960 cout <<
"Now checking recHit " << index << endl;
961 double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, eSetup);
962 cout <<
"Final distance is " << Distance << endl;
963 if(Distance > MaxRSD)
965 cout <<
"Pattern find error in distance for other recHits: " << Distance << endl;
972 cout <<
"Checking finish, isGoodPattern now is " << isGoodPattern << endl;
974 isPatternChecked =
true;
983 cout <<
"Extrapolating the track to check the pattern" << endl;
986 DetId hitDet = (*iter)->hit()->geographicalId();
994 double startSide = RPCSurface.localZ(startPosition);
995 cout <<
"Start side : " << startSide;
998 double currentSide = RPCSurface.localZ(currentPosition);
1003 double currentDistance = ((
GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
1004 cout <<
"Start current position is : " << currentPosition << endl;
1005 cout <<
"Start current Momentum is: " << currentMomentum.
mag() <<
", in vector: " << currentMomentum << endl;
1006 cout <<
"Start current distance is " << currentDistance << endl;
1007 cout <<
"Start current radius is " << currentPosition.
perp() << endl;
1008 cout <<
"Destination radius is " << (*iter)->globalPosition().perp() << endl;
1012 double currentDistance_next = currentDistance;
1015 currentDistance = currentDistance_next;
1016 if(ClockwiseDirection == 0)
1018 currentPosition += currentMomentum.
unit() * stepLength;
1022 double Bz = Field->inTesla(currentPosition).z();
1023 double Radius = currentMomentum.
perp()/fabs(Bz*0.01*0.3);
1024 double deltaPhi = (stepLength*currentMomentum.
perp()/currentMomentum.
mag())/Radius;
1028 currentPositiontoCenter *= Radius;
1030 currentPositiontoCenter *= ClockwiseDirection;
1033 currentCenter += currentPositiontoCenter;
1037 double Phi = CentertocurrentPosition.
phi().value();
1038 Phi += deltaPhi * (-1) * ClockwiseDirection;
1039 double deltaZ = stepLength*currentMomentum.
z()/currentMomentum.
mag();
1041 double PtPhi = currentMomentum.
phi().value();
1042 PtPhi += deltaPhi * (-1) * ClockwiseDirection;
1044 currentPosition = currentCenter + CentertonewPosition;
1048 tracklength += stepLength * currentMomentum.
perp() / currentMomentum.
mag();
1051 currentSide = RPCSurface.localZ(currentPosition);
1052 cout <<
"Stepping current side : " << currentSide << endl;
1053 cout <<
"Stepping current position is: " << currentPosition << endl;
1054 cout <<
"Stepping current Momentum is: " << currentMomentum.
mag() <<
", in vector: " << currentMomentum << endl;
1055 currentDistance_next = ((
GlobalVector)(currentPosition - (*iter)->globalPosition())).
perp();
1056 cout <<
"Stepping current distance is " << currentDistance << endl;
1057 cout <<
"Stepping current radius is " << currentPosition.
perp() << endl;
1058 }
while(currentDistance_next < currentDistance);
1060 return currentDistance;
1066 cout <<
"Now create a fake seed" << endl;
1067 isPatternChecked =
true;
1083 LocalVector segDirFromPos=best->det()->toLocal(Momentum);
1088 mat = best->parametersError().similarityT(best->projectionMatrix());
1089 mat[0][0] = meanSpt;
1097 DetId id = best->geographicalId();
1102 for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
1103 container.
push_back((*iter)->hit()->clone());
1107 theweightedSeed.first = theSeed;
1108 theweightedSeed.second = meanSpt;
1109 isGoodSeed = isGoodPattern;
1111 return theweightedSeed;
1116 if(isPatternChecked ==
false || isGoodPattern == -1)
1118 cout <<
"Pattern is not yet checked! Create a fake seed instead!" << endl;
1119 return createFakeSeed(isGoodSeed, eSetup);
1132 if(isClockwise == 0)
1135 Charge = (int)(meanPt / fabs(meanPt));
1141 if(isClockwise != 0)
1146 GlobalVector vecRef1((first->globalPosition().x()-Center.x()), (first->globalPosition().y()-Center.y()), (first->globalPosition().z()-Center.z()));
1147 GlobalVector vecRef2((best->globalPosition().x()-Center.x()), (best->globalPosition().y()-Center.y()), (best->globalPosition().z()-Center.z()));
1150 double deltaS = meanRadius * fabs(deltaPhi);
1151 double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1155 if(isClockwise == -1)
1160 Momentum *= fabs(meanPt / deltaS);
1164 double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1166 Momentum *= fabs(meanPt / S);
1171 Momentum = best->globalPosition() - first->globalPosition();
1172 double deltaS = Momentum.perp();
1173 Momentum *= fabs(meanPt / deltaS);
1176 LocalVector segDirFromPos = best->det()->toLocal(Momentum);
1182 cout <<
"Trajectory State on Surface before the extrapolation" << endl;
1184 DetId id = best->geographicalId();
1185 cout <<
"The RecSegment relies on: " << endl;
1192 for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
1198 container.
push_back((*iter)->hit()->clone());
1203 theweightedSeed.first = theSeed;
1204 theweightedSeed.second = meanSpt;
1205 isGoodSeed = isGoodPattern;
1207 return theweightedSeed;
1217 mat = best->parametersError().similarityT(best->projectionMatrix());
1219 GlobalVector vecRef1((first->globalPosition().x()-Center.x()), (first->globalPosition().y()-Center.y()), (first->globalPosition().z()-Center.z()));
1220 GlobalVector vecRef2((best->globalPosition().x()-Center.x()), (best->globalPosition().y()-Center.y()), (best->globalPosition().z()-Center.z()));
1222 double L = meanRadius * fabs(deltaPhi);
1224 double A_N = 180*N*N*N/((N-1)*(N+1)*(N+2)*(N+3));
1225 double sigma_x =
sqrt(mat[3][3]);
1226 double betaovergame = Momentum.mag()/0.1066;
1227 double beta =
sqrt((betaovergame*betaovergame)/(1+betaovergame*betaovergame));
1228 double dPt = meanPt*(0.0136*
sqrt(1/100)*
sqrt(4*A_N/N)/(beta*0.3*meanBz*
L) + sigma_x*meanPt*
sqrt(4*A_N)/(0.3*meanBz*L*
L));
1229 double dP = dPt * Momentum.mag() / meanPt;
1230 mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1231 mat[1][1] = dXdZ * dXdZ;
1232 mat[2][2] = dYdZ * dYdZ;
1237 mat0 = (SegmentRB[0].first)->parametersError().similarityT((SegmentRB[0].first)->projectionMatrix());
1238 double dX0 =
sqrt(mat0[3][3]);
1239 double dY0 =
sqrt(mat0[4][4]);
1241 mat1 = (SegmentRB[0].second)->parametersError().similarityT((SegmentRB[0].
second)->projectionMatrix());
1242 double dX1 =
sqrt(mat1[3][3]);
1243 double dY1 =
sqrt(mat1[4][4]);
1245 mat2 = (SegmentRB[1].first)->parametersError().similarityT((SegmentRB[1].first)->projectionMatrix());
1246 double dX2 =
sqrt(mat2[3][3]);
1247 double dY2 =
sqrt(mat2[4][4]);
1249 mat3 = (SegmentRB[1].second)->parametersError().similarityT((SegmentRB[1].second)->projectionMatrix());
1250 double dX3 =
sqrt(mat3[3][3]);
1251 double dY3 =
sqrt(mat3[4][4]);
1252 cout <<
"Local error for 4 recHits are: " << dX0 <<
", " << dY0 <<
", " << dX1 <<
", " << dY1 <<
", " << dX2 <<
", " << dY2 <<
", " << dX3 <<
", " << dY3 << endl;
1253 const GeomDetUnit* refRPC1 = (SegmentRB[0].second)->detUnit();
1263 const GeomDetUnit* refRPC2 = (SegmentRB[1].first)->detUnit();
1273 if(isClockwise != 0) {
1275 vec[0] =
GlobalVector((entryPosition.x()-Center2.x()), (entryPosition.y()-Center2.y()), (entryPosition.z()-Center2.z()));
1276 vec[1] =
GlobalVector((leavePosition.x()-Center2.x()), (leavePosition.y()-Center2.y()), (leavePosition.z()-Center2.z()));
1277 double halfPhiCenter = fabs((vec[1].
phi() - vec[0].
phi()).
value()) / 2;
1279 double dPhi0 = (((globalSegment00.phi() - globalSegment01.
phi()).
value()*isClockwise) > 0) ? fabs((globalSegment00.phi() - globalSegment01.
phi()).
value()) : fabs((globalSegment00.phi() - globalSegment02.
phi()).
value());
1280 double dPhi1 = (((globalSegment10.phi() - globalSegment11.
phi()).
value()*isClockwise) < 0) ? fabs((globalSegment10.phi() - globalSegment11.
phi()).
value()) : fabs((globalSegment10.phi() - globalSegment12.
phi()).
value());
1282 double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1283 cout <<
"DPhi for new Segment is about " << dPhi << endl;
1285 double newhalfPhiCenter = ((halfPhiCenter-
dPhi) > 0 ? (halfPhiCenter-dPhi) : 0);
1286 if(newhalfPhiCenter != 0) {
1287 double newmeanPt = meanPt * halfPhiCenter / newhalfPhiCenter;
1290 cout <<
"The error is inside range. Max meanPt could be " << newmeanPt << endl;
1291 dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1295 cout <<
"The error is outside range. Max meanPt could be " << newmeanPt << endl;
1296 dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1300 double dPhi0 = (fabs((globalSegment00.phi() - globalSegment01.
phi()).
value()) <= fabs((globalSegment00.phi() - globalSegment02.
phi()).
value())) ? fabs((globalSegment00.phi() - globalSegment01.
phi()).
value()) : fabs((globalSegment00.phi() - globalSegment02.
phi()).
value());
1301 double dPhi1 = (fabs((globalSegment10.phi() - globalSegment11.
phi()).
value()) <= fabs((globalSegment10.phi() - globalSegment12.
phi()).
value())) ? fabs((globalSegment10.phi() - globalSegment11.
phi()).
value()) : fabs((globalSegment10.phi() - globalSegment12.
phi()).
value());
1302 double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1303 GlobalVector middleSegment = leavePosition - entryPosition;
1304 double halfDistance = middleSegment.
perp() / 2;
1305 double newmeanPt = halfDistance /
dPhi;
1306 cout <<
"The error is for straight. Max meanPt could be " << newmeanPt << endl;
1307 dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1310 double dXdZ1 = globalSegment11.
x() / globalSegment11.
z() - globalSegment10.x() / globalSegment10.z();
1311 double dXdZ2 = globalSegment12.
x() / globalSegment12.
z() - globalSegment10.x() / globalSegment10.z();
1312 dXdZ = (fabs(dXdZ1) >= fabs(dXdZ2)) ? dXdZ1 : dXdZ2;
1318 double dYdZ1 = globalSegment13.
y() / globalSegment13.z() - globalSegment10.y() / globalSegment10.z();
1319 double dYdZ2 = globalSegment14.
y() / globalSegment14.
z() - globalSegment10.y() / globalSegment10.z();
1320 dYdZ = (fabs(dYdZ1) >= fabs(dYdZ2)) ? dYdZ1 : dYdZ2;
1322 mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1323 mat[1][1] = dXdZ * dXdZ;
1324 mat[2][2] = dYdZ * dYdZ;
T getParameter(std::string const &) const
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
void configure(const edm::ParameterSet &iConfig)
Local3DVector LocalVector
void checkSegmentAlgorithmSpecial(const edm::EventSetup &eSetup)
void MiddlePointsAlgorithm()
ConstMuonRecHitPointer FirstRecHit() const
bool checkStraightwithThreerecHits(ConstMuonRecHitPointer(&precHit)[3], double MinDeltaPhi) const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Geom::Phi< T > phi() const
weightedTrajectorySeed createSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
Global3DPoint GlobalPoint
double extropolateStep(const GlobalPoint &startPosition, const GlobalVector &startMomentum, ConstMuonRecHitContainer::const_iterator iter, const int ClockwiseDirection, double &tracklength, const edm::EventSetup &eSetup)
std::pair< ConstMuonRecHitPointer, ConstMuonRecHitPointer > RPCSegment
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
LocalTrajectoryError getSpecialAlgorithmErrorMatrix(const ConstMuonRecHitPointer &first, const ConstMuonRecHitPointer &best)
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
const Plane & surface() const
The nominal surface of the GeomDet.
std::string dumpMuonId(const DetId &id) const
uint32_t rawId() const
get the raw id
U second(std::pair< T, U > const &p)
const RPCChamber * chamber(RPCDetId id) const
std::string dumpTSOS(const TrajectoryStateOnSurface &tsos) const
double dPhi(double phi1, double phi2)
weightedTrajectorySeed seed(const edm::EventSetup &eSetup, int &isGoodSeed)
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
RPCDetId id() const
Return the RPCChamberId of this chamber.
void ThreePointsAlgorithm()
ConstMuonRecHitPointer BestRefRecHit() const
std::pair< TrajectorySeed, double > weightedTrajectorySeed
double deltaR(double eta1, double eta2, double phi1, double phi2)
Vector3DBase unit() const
void checkSimplePattern(const edm::EventSetup &eSetup)
void SegmentAlgorithmSpecial(const edm::EventSetup &eSetup)
GlobalVector computePtWithThreerecHits(double &pt, double &pt_err, double(&x)[3], double(&y)[3]) const
CLHEP::HepSymMatrix AlgebraicSymMatrix
bool checkSegment() const
bool checkStraightwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2, double MinDeltaPhi) const
GlobalVector computePtwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2) const
GlobalVector computePtwithThreerecHits(double &pt, double &pt_err, ConstMuonRecHitPointer(&precHit)[3]) const
Global3DVector GlobalVector
Basic3DVector cross(const Basic3DVector &v) const
Vector product, or "cross" product, with a vector of same type.
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.