CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Private Types | Private Member Functions | Private Attributes | Friends
RPCSeedPattern Class Reference

#include <RPCSeedPattern.h>

Public Types

typedef std::pair< ConstMuonRecHitPointer, ConstMuonRecHitPointerRPCSegment
 
typedef std::pair< TrajectorySeed, double > weightedTrajectorySeed
 

Public Member Functions

void add (const ConstMuonRecHitPointer &hit)
 
void clear ()
 
void configure (const edm::ParameterSet &iConfig)
 
unsigned int nrhit () const
 
 RPCSeedPattern ()
 
 ~RPCSeedPattern ()
 

Private Types

typedef MuonTransientTrackingRecHit::ConstMuonRecHitContainer ConstMuonRecHitContainer
 
typedef MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer
 
typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer
 
typedef MuonTransientTrackingRecHit::MuonRecHitPointer MuonRecHitPointer
 

Private Member Functions

ConstMuonRecHitPointer BestRefRecHit () const
 
bool checkSegment () const
 
void checkSegmentAlgorithmSpecial (const MagneticField &Field, const RPCGeometry &rpcGeometry)
 
void checkSimplePattern (const MagneticField &Field)
 
bool checkStraightwithSegment (const RPCSegment &Segment1, const RPCSegment &Segment2, double MinDeltaPhi) const
 
bool checkStraightwithThreerecHits (ConstMuonRecHitPointer(&precHit)[3], double MinDeltaPhi) const
 
bool checkStraightwithThreerecHits (double(&x)[3], double(&y)[3], double MinDeltaPhi) const
 
GlobalVector computePtwithSegment (const RPCSegment &Segment1, const RPCSegment &Segment2) const
 
GlobalVector computePtwithThreerecHits (double &pt, double &pt_err, ConstMuonRecHitPointer(&precHit)[3]) const
 
GlobalVector computePtWithThreerecHits (double &pt, double &pt_err, double(&x)[3], double(&y)[3]) const
 
weightedTrajectorySeed createFakeSeed (int &isGoodSeed, const MagneticField &)
 
weightedTrajectorySeed createSeed (int &isGoodSeed, const MagneticField &)
 
double extropolateStep (const GlobalPoint &startPosition, const GlobalVector &startMomentum, ConstMuonRecHitContainer::const_iterator iter, const int ClockwiseDirection, double &tracklength, const MagneticField &Field, const RPCGeometry &rpcGeometry)
 
ConstMuonRecHitPointer FirstRecHit () const
 
double getDistance (const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
 
LocalTrajectoryError getSpecialAlgorithmErrorMatrix (const ConstMuonRecHitPointer &first, const ConstMuonRecHitPointer &best)
 
void MiddlePointsAlgorithm ()
 
weightedTrajectorySeed seed (const MagneticField &Field, const RPCGeometry &rpcGeom, int &isGoodSeed)
 
void SegmentAlgorithm ()
 
void SegmentAlgorithmSpecial (const MagneticField &Field)
 
void ThreePointsAlgorithm ()
 

Private Attributes

unsigned int AlgorithmType
 
bool autoAlgorithmChoose
 
GlobalVector Center
 
GlobalVector Center2
 
int Charge
 
double deltaBz
 
double deltaRThreshold
 
GlobalPoint entryPosition
 
int isClockwise
 
bool isConfigured
 
int isGoodPattern
 
int isParralZ
 
bool isPatternChecked
 
bool isStraight
 
bool isStraight2
 
double lastPhi
 
GlobalPoint leavePosition
 
double MagnecticFieldThreshold
 
double MaxRSD
 
double meanBz
 
GlobalVector meanMagneticField2
 
double meanPt
 
double meanRadius
 
double meanRadius2
 
double meanSpt
 
double MinDeltaPhi
 
GlobalVector Momentum
 
double S
 
unsigned int sampleCount
 
RPCSegment SegmentRB [2]
 
double stepLength
 
ConstMuonRecHitContainer theRecHits
 
double ZError
 

Friends

class RPCSeedFinder
 

Detailed Description

Author
Haiyun.Teng - Peking University

Definition at line 30 of file RPCSeedPattern.h.

Member Typedef Documentation

◆ ConstMuonRecHitContainer

Definition at line 34 of file RPCSeedPattern.h.

◆ ConstMuonRecHitPointer

Definition at line 32 of file RPCSeedPattern.h.

◆ MuonRecHitContainer

Definition at line 33 of file RPCSeedPattern.h.

◆ MuonRecHitPointer

Definition at line 31 of file RPCSeedPattern.h.

◆ RPCSegment

Definition at line 37 of file RPCSeedPattern.h.

◆ weightedTrajectorySeed

Definition at line 38 of file RPCSeedPattern.h.

Constructor & Destructor Documentation

◆ RPCSeedPattern()

RPCSeedPattern::RPCSeedPattern ( )

Definition at line 30 of file RPCSeedPattern.cc.

30  {
31  isPatternChecked = false;
32  isConfigured = false;
34 }
double MagnecticFieldThreshold

◆ ~RPCSeedPattern()

RPCSeedPattern::~RPCSeedPattern ( )

Definition at line 36 of file RPCSeedPattern.cc.

36 {}

Member Function Documentation

◆ add()

void RPCSeedPattern::add ( const ConstMuonRecHitPointer hit)
inline

Definition at line 45 of file RPCSeedPattern.h.

References theRecHits.

Referenced by counter.Counter::register().

45 { theRecHits.push_back(hit); }
ConstMuonRecHitContainer theRecHits

◆ BestRefRecHit()

MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::BestRefRecHit ( ) const
private

Definition at line 467 of file RPCSeedPattern.cc.

467  {
469  int index = 0;
470  // Use the last one for recHit on last layer has minmum delta Z for barrel or delta R for endcap while calculating the momentum
471  // But for Algorithm 3 we use the 4th recHit on the 2nd segment for more accurate
472  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
473  if (AlgorithmType != 3)
474  best = (*iter);
475  else if (index < 4)
476  best = (*iter);
477  index++;
478  }
479  return best;
480 }
ConstMuonRecHitContainer theRecHits
unsigned int AlgorithmType
MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer

◆ checkSegment()

bool RPCSeedPattern::checkSegment ( ) const
private

Definition at line 437 of file RPCSeedPattern.cc.

References submitPVResolutionJobs::count, gather_cfg::cout, align::Detector, RPCChamber::id(), RPCDetId::region(), g4SimHits_cfi::Region, and RPCDetId::station().

437  {
438  bool isFit = true;
439  unsigned int count = 0;
440  // first 4 recHits should be located in RB1 and RB2
441  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
442  count++;
443  const GeomDet* Detector = (*iter)->det();
444  if (dynamic_cast<const RPCChamber*>(Detector) != nullptr) {
445  const RPCChamber* RPCCh = dynamic_cast<const RPCChamber*>(Detector);
446  RPCDetId RPCId = RPCCh->id();
447  int Region = RPCId.region();
448  unsigned int Station = RPCId.station();
449  //int Layer = RPCId.layer();
450  if (count <= 4) {
451  if (Region != 0)
452  isFit = false;
453  if (Station > 2)
454  isFit = false;
455  }
456  }
457  }
458  // more than 4 recHits for pattern building
459  if (count <= 4)
460  isFit = false;
461  cout << "Check for segment fit: " << isFit << endl;
462  return isFit;
463 }
RPCDetId id() const
Return the RPCChamberId of this chamber.
Definition: RPCChamber.cc:26
ConstMuonRecHitContainer theRecHits
int station() const
Definition: RPCDetId.h:78
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:53

◆ checkSegmentAlgorithmSpecial()

void RPCSeedPattern::checkSegmentAlgorithmSpecial ( const MagneticField Field,
const RPCGeometry rpcGeometry 
)
private

Definition at line 741 of file RPCSeedPattern.cc.

References gather_cfg::cout, PbPb_ZMuSkimMuonDPG_cff::deltaR, PV3DBase< T, PVType, FrameType >::perp(), mathSSE::sqrt(), upper_limit_pt, and relativeConstraints::value.

741  {
742  if (isPatternChecked == true)
743  return;
744 
745  if (!checkSegment()) {
746  isPatternChecked = true;
747  isGoodPattern = -1;
748  return;
749  }
750 
751  // Set isGoodPattern to default true and check the failure
752  isGoodPattern = 1;
753 
754  // Print the recHit's position
755  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
756  cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
757 
758  // Check the Z direction
759  if (fabs((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError) {
760  if (((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
761  isParralZ = 1;
762  else
763  isParralZ = -1;
764  } else
765  isParralZ = 0;
766 
767  cout << " Check isParralZ is :" << isParralZ << endl;
768  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
769  if (isParralZ == 0) {
770  if (fabs((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) > ZError) {
771  cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
772  isGoodPattern = 0;
773  }
774  } else {
775  if ((int)(((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) / ZError) * isParralZ < 0) {
776  cout << "Pattern find error in Z direction: wrong Z direction" << endl;
777  isGoodPattern = 0;
778  }
779  }
780  }
781 
782  // Check the pattern
783  if (isStraight2 == true) {
784  // Set pattern to be straight
785  isClockwise = 0;
787  // Set the straight pattern with lowest priority among good pattern
789 
790  // Extrapolate to other recHits and check deltaR
791  GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
792  GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
793  GlobalVector startMomentum = startSegment * (meanPt / startSegment.perp());
794  unsigned int index = 0;
795 
796  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
797  if (index < 4) {
798  index++;
799  continue;
800  }
801  double tracklength = 0;
802  cout << "Now checking recHit " << index << endl;
803  double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, Field, rpcGeom);
804  cout << "Final distance is " << Distance << endl;
805  if (Distance > MaxRSD) {
806  cout << "Pattern find error in distance for other recHits: " << Distance << endl;
807  isGoodPattern = 0;
808  }
809  index++;
810  }
811  } else {
812  // Get clockwise direction
813  GlobalVector vec[2];
814  vec[0] = GlobalVector(
815  (entryPosition.x() - Center2.x()), (entryPosition.y() - Center2.y()), (entryPosition.z() - Center2.z()));
816  vec[1] = GlobalVector(
817  (leavePosition.x() - Center2.x()), (leavePosition.y() - Center2.y()), (leavePosition.z() - Center2.z()));
818  isClockwise = 0;
819  if ((vec[1].phi() - vec[0].phi()).value() > 0)
820  isClockwise = -1;
821  else
822  isClockwise = 1;
823 
824  cout << "Check isClockwise is : " << isClockwise << endl;
825 
826  // Get meanPt
827  meanPt = 0.01 * meanRadius2 * meanMagneticField2.z() * 0.3 * isClockwise;
828  //meanPt = 0.01 * meanRadius2[0] * (-3.8) * 0.3 * isClockwise;
829  cout << " meanRadius is " << meanRadius2 << ", with meanBz " << meanMagneticField2.z() << endl;
830  cout << " meanPt is " << meanPt << endl;
831  if (fabs(meanPt) > upper_limit_pt)
832  meanPt = upper_limit_pt * meanPt / fabs(meanPt);
833 
834  // Check the initial 3 segments
835  cout << "entryPosition: " << entryPosition << endl;
836  cout << "leavePosition: " << leavePosition << endl;
837  cout << "Center2 is : " << Center2 << endl;
838  double R1 = vec[0].perp();
839  double R2 = vec[1].perp();
840  double deltaR = sqrt(((R1 - meanRadius2) * (R1 - meanRadius2) + (R2 - meanRadius2) * (R2 - meanRadius2)) / 2);
841  meanSpt = deltaR;
842  cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
843  cout << "Delta R for the initial 3 segments is " << deltaR << endl;
844  if (deltaR > deltaRThreshold) {
845  cout << "Pattern find error in delta R for the initial 3 segments" << endl;
846  isGoodPattern = 0;
847  }
848 
849  // Extrapolate to other recHits and check deltaR
850  GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
851  GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
852  GlobalVector startMomentum = startSegment * (fabs(meanPt) / startSegment.perp());
853  unsigned int index = 0;
854  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
855  if (index < 4) {
856  index++;
857  continue;
858  }
859  double tracklength = 0;
860  cout << "Now checking recHit " << index << endl;
861  double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, Field, rpcGeom);
862  cout << "Final distance is " << Distance << endl;
863  if (Distance > MaxRSD) {
864  cout << "Pattern find error in distance for other recHits: " << Distance << endl;
865  isGoodPattern = 0;
866  }
867  index++;
868  }
869  }
870 
871  cout << "Checking finish, isGoodPattern now is " << isGoodPattern << endl;
872  // Set the pattern estimation signal
873  isPatternChecked = true;
874 }
double extropolateStep(const GlobalPoint &startPosition, const GlobalVector &startMomentum, ConstMuonRecHitContainer::const_iterator iter, const int ClockwiseDirection, double &tracklength, const MagneticField &Field, const RPCGeometry &rpcGeometry)
T perp() const
Definition: PV3DBase.h:69
GlobalPoint entryPosition
T z() const
Definition: PV3DBase.h:61
double deltaRThreshold
RPCSegment SegmentRB[2]
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
bool checkSegment() const
ConstMuonRecHitContainer theRecHits
GlobalVector Center2
GlobalVector meanMagneticField2
GlobalPoint leavePosition
#define upper_limit_pt
Global3DVector GlobalVector
Definition: GlobalVector.h:10

◆ checkSimplePattern()

void RPCSeedPattern::checkSimplePattern ( const MagneticField Field)
private

Definition at line 604 of file RPCSeedPattern.cc.

References funct::abs(), gather_cfg::cout, PbPb_ZMuSkimMuonDPG_cff::deltaR, PVValHelper::dx, PVValHelper::dy, PVValHelper::dz, MagneticField::inTesla(), mathSSE::sqrt(), upper_limit_pt, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

604  {
605  if (isPatternChecked == true)
606  return;
607 
608  unsigned int NumberofHitsinSeed = nrhit();
609 
610  // Print the recHit's position
611  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
612  cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
613 
614  // Get magnetice field information
615  std::vector<double> BzValue;
616  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
617  GlobalPoint gpFirst = (*iter)->globalPosition();
618  GlobalPoint gpLast = (*(iter + 1))->globalPosition();
620  double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
621  double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
622  double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
623  for (unsigned int index = 0; index < sampleCount; index++) {
624  gp[index] = GlobalPoint(
625  (gpFirst.x() + dx * (index + 1)), (gpFirst.y() + dy * (index + 1)), (gpFirst.z() + dz * (index + 1)));
626  GlobalVector MagneticVec_temp = Field.inTesla(gp[index]);
627  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
628  BzValue.push_back(MagneticVec_temp.z());
629  }
630  delete[] gp;
631  }
632  meanBz = 0.;
633  for (unsigned int index = 0; index < BzValue.size(); index++)
634  meanBz += BzValue[index];
635  meanBz /= BzValue.size();
636  cout << "Mean Bz is " << meanBz << endl;
637  deltaBz = 0.;
638  for (unsigned int index = 0; index < BzValue.size(); index++)
639  deltaBz += (BzValue[index] - meanBz) * (BzValue[index] - meanBz);
640  deltaBz /= BzValue.size();
641  deltaBz = sqrt(deltaBz);
642  cout << "delata Bz is " << deltaBz << endl;
643 
644  // Set isGoodPattern to default true and check the failure
645  isGoodPattern = 1;
646 
647  // Check the Z direction
648  if (fabs((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError) {
649  if (((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
650  isParralZ = 1;
651  else
652  isParralZ = -1;
653  } else
654  isParralZ = 0;
655 
656  cout << " Check isParralZ is :" << isParralZ << endl;
657  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
658  if (isParralZ == 0) {
659  if (fabs((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) > ZError) {
660  cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
661  isGoodPattern = 0;
662  }
663  } else {
664  if ((int)(((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) / ZError) * isParralZ < 0) {
665  cout << "Pattern find error in Z direction: wrong Z direction" << endl;
666  isGoodPattern = 0;
667  }
668  }
669  }
670 
671  // Check pattern
672  if (isStraight == false) {
673  // Check clockwise direction
674  GlobalVector* vec = new GlobalVector[NumberofHitsinSeed];
675  unsigned int index = 0;
676  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
677  GlobalVector vec_temp(((*iter)->globalPosition().x() - Center.x()),
678  ((*iter)->globalPosition().y() - Center.y()),
679  ((*iter)->globalPosition().z() - Center.z()));
680  vec[index] = vec_temp;
681  index++;
682  }
683  isClockwise = 0;
684  for (unsigned int index = 0; index < (NumberofHitsinSeed - 1); index++) {
685  // Check phi direction, all sub-dphi direction should be the same
686  if ((vec[index + 1].phi() - vec[index].phi()) > 0)
687  isClockwise--;
688  else
689  isClockwise++;
690  cout << "Current isClockwise is : " << isClockwise << endl;
691  }
692  cout << "Check isClockwise is : " << isClockwise << endl;
693  if ((unsigned int)abs(isClockwise) != (NumberofHitsinSeed - 1)) {
694  cout << "Pattern find error in Phi direction" << endl;
695  isGoodPattern = 0;
696  isClockwise = 0;
697  } else
699  delete[] vec;
700 
701  // Get meanPt and meanSpt
702  double deltaRwithBz = fabs(deltaBz * meanRadius / meanBz);
703  cout << "deltaR with Bz is " << deltaRwithBz << endl;
704 
705  if (isClockwise == 0)
707  else
708  meanPt = 0.01 * meanRadius * meanBz * 0.3 * isClockwise;
709  if (fabs(meanPt) > upper_limit_pt)
710  meanPt = upper_limit_pt * meanPt / fabs(meanPt);
711  cout << " meanRadius is " << meanRadius << endl;
712  cout << " meanPt is " << meanPt << endl;
713 
714  double deltaR = 0.;
715  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
716  deltaR += (getDistance(*iter, Center) - meanRadius) * (getDistance(*iter, Center) - meanRadius);
717  }
718  deltaR = deltaR / NumberofHitsinSeed;
719  deltaR = sqrt(deltaR);
720  //meanSpt = 0.01 * deltaR * meanBz * 0.3;
721  meanSpt = deltaR;
722  cout << "DeltaR is " << deltaR << endl;
723  if (deltaR > deltaRThreshold) {
724  cout << "Pattern find error: deltaR over threshold" << endl;
725  isGoodPattern = 0;
726  }
727  } else {
728  // Just set pattern to be straight
729  isClockwise = 0;
731  // Set the straight pattern with lowest priority among good pattern
733  }
734  cout << "III--> Seed Pt : " << meanPt << endl;
735  cout << "III--> Pattern is: " << isGoodPattern << endl;
736 
737  // Set the pattern estimation signal
738  isPatternChecked = true;
739 }
T z() const
Definition: PV3DBase.h:61
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
double deltaRThreshold
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
ConstMuonRecHitContainer theRecHits
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
GlobalVector Center
unsigned int nrhit() const
unsigned int sampleCount
#define upper_limit_pt

◆ checkStraightwithSegment()

bool RPCSeedPattern::checkStraightwithSegment ( const RPCSegment Segment1,
const RPCSegment Segment2,
double  MinDeltaPhi 
) const
private

Definition at line 525 of file RPCSeedPattern.cc.

References gather_cfg::cout, PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, PV3DBase< T, PVType, FrameType >::phi(), and relativeConstraints::value.

527  {
528  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
529  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
530  GlobalVector segvec3 = (Segment2.first)->globalPosition() - (Segment1.first)->globalPosition();
531  // compare segvec 1&2 for paralle, 1&3 for straight
532  double dPhi1 = (segvec2.phi() - segvec1.phi()).value();
533  double dPhi2 = (segvec3.phi() - segvec1.phi()).value();
534  cout << "Checking straight with 2 segments. dPhi1: " << dPhi1 << ", dPhi2: " << dPhi2 << endl;
535  cout << "Checking straight with 2 segments. dPhi1 in degree: " << dPhi1 * 180 / 3.1415926
536  << ", dPhi2 in degree: " << dPhi2 * 180 / 3.1415926 << endl;
537  if (fabs(dPhi1) > MinDeltaPhi || fabs(dPhi2) > MinDeltaPhi) {
538  cout << "Segment is estimate to be not straight" << endl;
539  return false;
540  } else {
541  cout << "Segment is estimate to be straight" << endl;
542  return true;
543  }
544 }
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66

◆ checkStraightwithThreerecHits() [1/2]

bool RPCSeedPattern::checkStraightwithThreerecHits ( ConstMuonRecHitPointer(&)  precHit[3],
double  MinDeltaPhi 
) const
private

Definition at line 487 of file RPCSeedPattern.cc.

References gather_cfg::cout, PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, PV3DBase< T, PVType, FrameType >::phi(), and relativeConstraints::value.

487  {
488  GlobalVector segvec1 = precHit[1]->globalPosition() - precHit[0]->globalPosition();
489  GlobalVector segvec2 = precHit[2]->globalPosition() - precHit[1]->globalPosition();
490  double dPhi = (segvec2.phi() - segvec1.phi()).value();
491  if (fabs(dPhi) > MinDeltaPhi) {
492  cout << "Part is estimate to be not straight" << endl;
493  return false;
494  } else {
495  cout << "Part is estimate to be straight" << endl;
496  return true;
497  }
498 }
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66

◆ checkStraightwithThreerecHits() [2/2]

bool RPCSeedPattern::checkStraightwithThreerecHits ( double(&)  x[3],
double(&)  y[3],
double  MinDeltaPhi 
) const
private

Definition at line 572 of file RPCSeedPattern.cc.

References gather_cfg::cout, PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, PV3DBase< T, PVType, FrameType >::phi(), relativeConstraints::value, and x.

572  {
573  GlobalVector segvec1((x[1] - x[0]), (y[1] - y[0]), 0);
574  GlobalVector segvec2((x[2] - x[1]), (y[2] - y[1]), 0);
575  double dPhi = (segvec2.phi() - segvec1.phi()).value();
576  if (fabs(dPhi) > MinDeltaPhi) {
577  cout << "Part is estimate to be not straight" << endl;
578  return false;
579  } else {
580  cout << "Part is estimate to be straight" << endl;
581  return true;
582  }
583 }

◆ clear()

void RPCSeedPattern::clear ( void  )
inline

Definition at line 44 of file RPCSeedPattern.h.

References theRecHits.

44 { theRecHits.clear(); }
ConstMuonRecHitContainer theRecHits

◆ computePtwithSegment()

GlobalVector RPCSeedPattern::computePtwithSegment ( const RPCSegment Segment1,
const RPCSegment Segment2 
) const
private

Definition at line 546 of file RPCSeedPattern.cc.

References Vector3DBase< T, FrameTag >::cross(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

546  {
547  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
548  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
549  GlobalPoint Point1(((Segment1.second)->globalPosition().x() + (Segment1.first)->globalPosition().x()) / 2,
550  ((Segment1.second)->globalPosition().y() + (Segment1.first)->globalPosition().y()) / 2,
551  ((Segment1.second)->globalPosition().z() + (Segment1.first)->globalPosition().z()) / 2);
552  GlobalPoint Point2(((Segment2.second)->globalPosition().x() + (Segment2.first)->globalPosition().x()) / 2,
553  ((Segment2.second)->globalPosition().y() + (Segment2.first)->globalPosition().y()) / 2,
554  ((Segment2.second)->globalPosition().z() + (Segment2.first)->globalPosition().z()) / 2);
555  GlobalVector vecZ(0, 0, 1);
556  GlobalVector gvec1 = segvec1.cross(vecZ);
557  GlobalVector gvec2 = segvec2.cross(vecZ);
558  double A1 = gvec1.x();
559  double B1 = gvec1.y();
560  double A2 = gvec2.x();
561  double B2 = gvec2.y();
562  double X1 = Point1.x();
563  double Y1 = Point1.y();
564  double X2 = Point2.x();
565  double Y2 = Point2.y();
566  double XO = (A1 * A2 * (Y2 - Y1) + A2 * B1 * X1 - A1 * B2 * X2) / (A2 * B1 - A1 * B2);
567  double YO = (B1 * B2 * (X2 - X1) + B2 * A1 * Y1 - B1 * A2 * Y2) / (B2 * A1 - B1 * A2);
568  GlobalVector Center(XO, YO, 0);
569  return Center;
570 }
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:110
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
GlobalVector Center

◆ computePtwithThreerecHits()

GlobalVector RPCSeedPattern::computePtwithThreerecHits ( double &  pt,
double &  pt_err,
ConstMuonRecHitPointer(&)  precHit[3] 
) const
private

Definition at line 500 of file RPCSeedPattern.cc.

References A, gather_cfg::cout, DiDispStaMuonMonitor_cfi::pt, mathSSE::sqrt(), and x.

502  {
503  double x[3], y[3];
504  x[0] = precHit[0]->globalPosition().x();
505  y[0] = precHit[0]->globalPosition().y();
506  x[1] = precHit[1]->globalPosition().x();
507  y[1] = precHit[1]->globalPosition().y();
508  x[2] = precHit[2]->globalPosition().x();
509  y[2] = precHit[2]->globalPosition().y();
510  double A = (y[2] - y[1]) / (x[2] - x[1]) - (y[1] - y[0]) / (x[1] - x[0]);
511  double TYO = (x[2] - x[0]) / A + (y[2] * y[2] - y[1] * y[1]) / ((x[2] - x[1]) * A) -
512  (y[1] * y[1] - y[0] * y[0]) / ((x[1] - x[0]) * A);
513  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]);
514  double XO = 0.5 * TXO;
515  double YO = 0.5 * TYO;
516  double R2 = (x[0] - XO) * (x[0] - XO) + (y[0] - YO) * (y[0] - YO);
517  cout << "R2 is " << R2 << endl;
518  // How this algorithm get the pt without magnetic field??
519  pt = 0.01 * sqrt(R2) * 2 * 0.3;
520  cout << "pt is " << pt << endl;
521  GlobalVector Center(XO, YO, 0);
522  return Center;
523 }
T sqrt(T t)
Definition: SSEVec.h:19
GlobalVector Center
Definition: APVGainStruct.h:7

◆ computePtWithThreerecHits()

GlobalVector RPCSeedPattern::computePtWithThreerecHits ( double &  pt,
double &  pt_err,
double(&)  x[3],
double(&)  y[3] 
) const
private

Definition at line 585 of file RPCSeedPattern.cc.

References A, gather_cfg::cout, DiDispStaMuonMonitor_cfi::pt, mathSSE::sqrt(), and x.

588  {
589  double A = (y[2] - y[1]) / (x[2] - x[1]) - (y[1] - y[0]) / (x[1] - x[0]);
590  double TYO = (x[2] - x[0]) / A + (y[2] * y[2] - y[1] * y[1]) / ((x[2] - x[1]) * A) -
591  (y[1] * y[1] - y[0] * y[0]) / ((x[1] - x[0]) * A);
592  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]);
593  double XO = 0.5 * TXO;
594  double YO = 0.5 * TYO;
595  double R2 = (x[0] - XO) * (x[0] - XO) + (y[0] - YO) * (y[0] - YO);
596  cout << "R2 is " << R2 << endl;
597  // How this algorithm get the pt without magnetic field??
598  pt = 0.01 * sqrt(R2) * 2 * 0.3;
599  cout << "pt is " << pt << endl;
600  GlobalVector Center(XO, YO, 0);
601  return Center;
602 }
T sqrt(T t)
Definition: SSEVec.h:19
GlobalVector Center
Definition: APVGainStruct.h:7

◆ configure()

void RPCSeedPattern::configure ( const edm::ParameterSet iConfig)

Definition at line 38 of file RPCSeedPattern.cc.

References edm::ParameterSet::getParameter(), and PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi.

38  {
39  MaxRSD = iConfig.getParameter<double>("MaxRSD");
40  deltaRThreshold = iConfig.getParameter<double>("deltaRThreshold");
41  AlgorithmType = iConfig.getParameter<unsigned int>("AlgorithmType");
42  autoAlgorithmChoose = iConfig.getParameter<bool>("autoAlgorithmChoose");
43  ZError = iConfig.getParameter<double>("ZError");
44  MinDeltaPhi = iConfig.getParameter<double>("MinDeltaPhi");
45  MagnecticFieldThreshold = iConfig.getParameter<double>("MagnecticFieldThreshold");
46  stepLength = iConfig.getParameter<double>("stepLength");
47  sampleCount = iConfig.getParameter<unsigned int>("sampleCount");
48  isConfigured = true;
49 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
double MagnecticFieldThreshold
bool autoAlgorithmChoose
double deltaRThreshold
unsigned int AlgorithmType
unsigned int sampleCount

◆ createFakeSeed()

RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createFakeSeed ( int &  isGoodSeed,
const MagneticField Field 
)
private

Definition at line 955 of file RPCSeedPattern.cc.

References alongMomentum, PixelTestBeamValidation_cfi::Charge, gather_cfg::cout, relativeConstraints::error, trajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), nano_mu_digi_cff::rawId, and upper_limit_pt.

955  {
956  // Create a fake seed and return
957  cout << "Now create a fake seed" << endl;
958  isPatternChecked = true;
959  isGoodPattern = -1;
960  isStraight = true;
962  meanSpt = 0;
963  Charge = 0;
964  isClockwise = 0;
965  isParralZ = 0;
966  meanRadius = -1;
967  //return createSeed(isGoodSeed, eSetup);
968 
969  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
971 
972  Momentum = GlobalVector(0, 0, 0);
973  LocalPoint segPos = best->localPosition();
974  LocalVector segDirFromPos = best->det()->toLocal(Momentum);
975  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
976 
977  //AlgebraicVector t(4);
978  AlgebraicSymMatrix mat(5, 0);
979  mat = best->parametersError().similarityT(best->projectionMatrix());
980  mat[0][0] = meanSpt;
981  LocalTrajectoryError error(asSMatrix<5>(mat));
982 
983  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &Field);
984 
985  DetId id = best->geographicalId();
986 
988 
990  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
991  container.push_back((*iter)->hit()->clone());
992 
993  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
994  weightedTrajectorySeed theweightedSeed;
995  theweightedSeed.first = theSeed;
996  theweightedSeed.second = meanSpt;
997  isGoodSeed = isGoodPattern;
998 
999  return theweightedSeed;
1000 }
ConstMuonRecHitPointer BestRefRecHit() const
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
GlobalVector Momentum
void push_back(D *&d)
Definition: OwnVector.h:326
ConstMuonRecHitContainer theRecHits
Definition: DetId.h:17
RPCSeedPattern::weightedTrajectorySeed weightedTrajectorySeed
CLHEP::HepSymMatrix AlgebraicSymMatrix
#define upper_limit_pt
Global3DVector GlobalVector
Definition: GlobalVector.h:10
MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer

◆ createSeed()

RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createSeed ( int &  isGoodSeed,
const MagneticField Field 
)
private

Definition at line 1002 of file RPCSeedPattern.cc.

References alongMomentum, PixelTestBeamValidation_cfi::Charge, gather_cfg::cout, cross(), debug, SiPixelRawToDigiRegional_cfi::deltaPhi, l1tTrackerHTMiss_cfi::deltaZ, relativeConstraints::error, dqmdumpme::first, createfilelist::int, trajectoryStateTransform::persistentState(), PV3DBase< T, PVType, FrameType >::phi(), edm::OwnVector< T, P >::push_back(), nano_mu_digi_cff::rawId, Vector3DBase< T, FrameTag >::unit(), and relativeConstraints::value.

1002  {
1003  if (isPatternChecked == false || isGoodPattern == -1) {
1004  cout << "Pattern is not yet checked! Create a fake seed instead!" << endl;
1005  return createFakeSeed(isGoodSeed, Field);
1006  }
1007 
1009 
1010  //double theMinMomentum = 3.0;
1011  //if(fabs(meanPt) < lower_limit_pt)
1012  //meanPt = lower_limit_pt * meanPt / fabs(meanPt);
1013 
1014  // For pattern we use is Clockwise other than isStraight to estimate charge
1015  if (isClockwise == 0)
1016  Charge = 0;
1017  else
1018  Charge = (int)(meanPt / fabs(meanPt));
1019 
1020  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
1023 
1024  if (isClockwise != 0) {
1025  if (AlgorithmType != 3) {
1026  // Get the momentum on reference recHit
1027  GlobalVector vecRef1((first->globalPosition().x() - Center.x()),
1028  (first->globalPosition().y() - Center.y()),
1029  (first->globalPosition().z() - Center.z()));
1030  GlobalVector vecRef2((best->globalPosition().x() - Center.x()),
1031  (best->globalPosition().y() - Center.y()),
1032  (best->globalPosition().z() - Center.z()));
1033 
1034  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1035  double deltaS = meanRadius * fabs(deltaPhi);
1036  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1037 
1038  GlobalVector vecZ(0, 0, 1);
1039  GlobalVector vecPt = (vecRef2.unit()).cross(vecZ);
1040  if (isClockwise == -1)
1041  vecPt *= -1;
1042  vecPt *= deltaS;
1043  Momentum = GlobalVector(0, 0, deltaZ);
1044  Momentum += vecPt;
1045  Momentum *= fabs(meanPt / deltaS);
1046  } else {
1047  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1049  Momentum *= fabs(meanPt / S);
1050  }
1051  } else {
1052  Momentum = best->globalPosition() - first->globalPosition();
1053  double deltaS = Momentum.perp();
1054  Momentum *= fabs(meanPt / deltaS);
1055  }
1056  LocalPoint segPos = best->localPosition();
1057  LocalVector segDirFromPos = best->det()->toLocal(Momentum);
1058  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
1059 
1061 
1062  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &Field);
1063  cout << "Trajectory State on Surface before the extrapolation" << endl;
1064  cout << debug.dumpTSOS(tsos);
1065  DetId id = best->geographicalId();
1066  cout << "The RecSegment relies on: " << endl;
1067  cout << debug.dumpMuonId(id);
1068  cout << debug.dumpTSOS(tsos);
1069 
1071 
1073  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
1074  // This casting withou clone will cause memory overflow when used in push_back
1075  // Since container's deconstructor functiion free the pointer menber!
1076  //TrackingRecHit* pt = dynamic_cast<TrackingRecHit*>(&*(*iter));
1077  //cout << "Push recHit type " << pt->getType() << endl;
1078  container.push_back((*iter)->hit()->clone());
1079  }
1080 
1081  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
1082  weightedTrajectorySeed theweightedSeed;
1083  theweightedSeed.first = theSeed;
1084  theweightedSeed.second = meanSpt;
1085  isGoodSeed = isGoodPattern;
1086 
1087  return theweightedSeed;
1088 }
Basic3DVector cross(const Basic3DVector &v) const
Vector product, or "cross" product, with a vector of same type.
T perp() const
Definition: PV3DBase.h:69
ConstMuonRecHitPointer BestRefRecHit() const
T z() const
Definition: PV3DBase.h:61
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
GlobalVector Momentum
LocalTrajectoryError getSpecialAlgorithmErrorMatrix(const ConstMuonRecHitPointer &first, const ConstMuonRecHitPointer &best)
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
void push_back(D *&d)
Definition: OwnVector.h:326
ConstMuonRecHitContainer theRecHits
GlobalVector Center
ConstMuonRecHitPointer FirstRecHit() const
Definition: DetId.h:17
#define debug
Definition: HDRShower.cc:19
RPCSeedPattern::weightedTrajectorySeed weightedTrajectorySeed
unsigned int AlgorithmType
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const MagneticField &)
Global3DVector GlobalVector
Definition: GlobalVector.h:10
MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer

◆ extropolateStep()

double RPCSeedPattern::extropolateStep ( const GlobalPoint startPosition,
const GlobalVector startMomentum,
ConstMuonRecHitContainer::const_iterator  iter,
const int  ClockwiseDirection,
double &  tracklength,
const MagneticField Field,
const RPCGeometry rpcGeometry 
)
private

Definition at line 876 of file RPCSeedPattern.cc.

References RPCGeometry::chamber(), gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), SiPixelRawToDigiRegional_cfi::deltaPhi, l1tTrackerHTMiss_cfi::deltaZ, MagneticField::inTesla(), Plane::localZ(), PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), VtxSmearedParameters_cfi::Phi, DetId::rawId(), GeomDet::surface(), Vector3DBase< T, FrameTag >::unit(), and PV3DBase< T, PVType, FrameType >::z().

882  {
883  cout << "Extrapolating the track to check the pattern" << endl;
884  tracklength = 0;
885  // Get the iter recHit's detector geometry
886  DetId hitDet = (*iter)->hit()->geographicalId();
887  RPCDetId RPCId = RPCDetId(hitDet.rawId());
888  //const RPCChamber* hitRPC = dynamic_cast<const RPCChamber*>(hitDet);
889 
890  const BoundPlane RPCSurface = rpcGeometry.chamber(RPCId)->surface();
891  double startSide = RPCSurface.localZ(startPosition);
892  cout << "Start side : " << startSide;
893 
894  GlobalPoint currentPosition = startPosition;
895  GlobalVector currentMomentum = startMomentum;
896  GlobalVector ZDirection(0, 0, 1);
897 
898  // Use the perp other than mag, since initial segment might have small value while final recHit have large difference value at Z direction
899  double currentDistance = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
900  cout << "Start current position is : " << currentPosition << endl;
901  cout << "Start current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
902  cout << "Start current distance is " << currentDistance << endl;
903  cout << "Start current radius is " << currentPosition.perp() << endl;
904  cout << "Destination radius is " << (*iter)->globalPosition().perp() << endl;
905 
906  // Judge roughly if the stepping cross the Det surface of the recHit
907  //while((currentPosition.perp() < ((*iter)->globalPosition().perp())))
908  double currentDistance_next = currentDistance;
909  do {
910  currentDistance = currentDistance_next;
911  if (ClockwiseDirection == 0) {
912  currentPosition += currentMomentum.unit() * stepLength;
913  } else {
914  double Bz = Field.inTesla(currentPosition).z();
915  double Radius = currentMomentum.perp() / fabs(Bz * 0.01 * 0.3);
916  double deltaPhi = (stepLength * currentMomentum.perp() / currentMomentum.mag()) / Radius;
917 
918  // Get the center for current step
919  GlobalVector currentPositiontoCenter = currentMomentum.unit().cross(ZDirection);
920  currentPositiontoCenter *= Radius;
921  // correction of ClockwiseDirection correction
922  currentPositiontoCenter *= ClockwiseDirection;
923  // continue to get the center for current step
924  GlobalPoint currentCenter = currentPosition;
925  currentCenter += currentPositiontoCenter;
926 
927  // Get the next step position
928  GlobalVector CentertocurrentPosition = (GlobalVector)(currentPosition - currentCenter);
929  double Phi = CentertocurrentPosition.phi().value();
930  Phi += deltaPhi * (-1) * ClockwiseDirection;
931  double deltaZ = stepLength * currentMomentum.z() / currentMomentum.mag();
932  GlobalVector CentertonewPosition(GlobalVector::Cylindrical(CentertocurrentPosition.perp(), Phi, deltaZ));
933  double PtPhi = currentMomentum.phi().value();
934  PtPhi += deltaPhi * (-1) * ClockwiseDirection;
935  currentMomentum = GlobalVector(GlobalVector::Cylindrical(currentMomentum.perp(), PtPhi, currentMomentum.z()));
936  currentPosition = currentCenter + CentertonewPosition;
937  }
938 
939  // count the total step length
940  tracklength += stepLength * currentMomentum.perp() / currentMomentum.mag();
941 
942  // Get the next step distance
943  double currentSide = RPCSurface.localZ(currentPosition);
944  cout << "Stepping current side : " << currentSide << endl;
945  cout << "Stepping current position is: " << currentPosition << endl;
946  cout << "Stepping current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
947  currentDistance_next = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
948  cout << "Stepping current distance is " << currentDistance << endl;
949  cout << "Stepping current radius is " << currentPosition.perp() << endl;
950  } while (currentDistance_next < currentDistance);
951 
952  return currentDistance;
953 }
T perp() const
Definition: PV3DBase.h:69
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:110
T z() const
Definition: PV3DBase.h:61
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
T mag() const
Definition: PV3DBase.h:64
float localZ(const GlobalPoint &gp) const
Definition: Plane.h:45
Definition: DetId.h:17
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
const RPCChamber * chamber(RPCDetId id) const
Definition: RPCGeometry.cc:46
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
Vector3DBase unit() const
Definition: Vector3DBase.h:54
Global3DVector GlobalVector
Definition: GlobalVector.h:10

◆ FirstRecHit()

MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::FirstRecHit ( ) const
private

Definition at line 465 of file RPCSeedPattern.cc.

465 { return theRecHits.front(); }
ConstMuonRecHitContainer theRecHits

◆ getDistance()

double RPCSeedPattern::getDistance ( const ConstMuonRecHitPointer precHit,
const GlobalVector Center 
) const
private

Definition at line 482 of file RPCSeedPattern.cc.

References mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

482  {
483  return sqrt((precHit->globalPosition().x() - Center.x()) * (precHit->globalPosition().x() - Center.x()) +
484  (precHit->globalPosition().y() - Center.y()) * (precHit->globalPosition().y() - Center.y()));
485 }
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
GlobalVector Center

◆ getSpecialAlgorithmErrorMatrix()

LocalTrajectoryError RPCSeedPattern::getSpecialAlgorithmErrorMatrix ( const ConstMuonRecHitPointer first,
const ConstMuonRecHitPointer best 
)
private

Definition at line 1090 of file RPCSeedPattern.cc.

References HLT_2023v12_cff::beta, gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, dqmdumpme::first, dttmaxenums::L, N, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), edm::second(), MuonME0Digis_cfi::sigma_x, mathSSE::sqrt(), GeomDet::toGlobal(), GeomDet::toLocal(), upper_limit_pt, relativeConstraints::value, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

1091  {
1093  double dXdZ = 0;
1094  double dYdZ = 0;
1095  double dP = 0;
1096  AlgebraicSymMatrix mat(5, 0);
1097  mat = best->parametersError().similarityT(best->projectionMatrix());
1098  if (AlgorithmType != 3) {
1099  GlobalVector vecRef1((first->globalPosition().x() - Center.x()),
1100  (first->globalPosition().y() - Center.y()),
1101  (first->globalPosition().z() - Center.z()));
1102  GlobalVector vecRef2((best->globalPosition().x() - Center.x()),
1103  (best->globalPosition().y() - Center.y()),
1104  (best->globalPosition().z() - Center.z()));
1105  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1106  double L = meanRadius * fabs(deltaPhi);
1107  double N = nrhit();
1108  double A_N = 180 * N * N * N / ((N - 1) * (N + 1) * (N + 2) * (N + 3));
1109  double sigma_x = sqrt(mat[3][3]);
1110  double betaovergame = Momentum.mag() / 0.1066;
1111  double beta = sqrt((betaovergame * betaovergame) / (1 + betaovergame * betaovergame));
1112  double dPt = meanPt * (0.0136 * sqrt(1 / 100) * sqrt(4 * A_N / N) / (beta * 0.3 * meanBz * L) +
1113  sigma_x * meanPt * sqrt(4 * A_N) / (0.3 * meanBz * L * L));
1114  double dP = dPt * Momentum.mag() / meanPt;
1115  mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1116  mat[1][1] = dXdZ * dXdZ;
1117  mat[2][2] = dYdZ * dYdZ;
1118  Error = LocalTrajectoryError(asSMatrix<5>(mat));
1119  } else {
1120  AlgebraicSymMatrix mat0(5, 0);
1121  mat0 = (SegmentRB[0].first)->parametersError().similarityT((SegmentRB[0].first)->projectionMatrix());
1122  double dX0 = sqrt(mat0[3][3]);
1123  double dY0 = sqrt(mat0[4][4]);
1124  AlgebraicSymMatrix mat1(5, 0);
1125  mat1 = (SegmentRB[0].second)->parametersError().similarityT((SegmentRB[0].second)->projectionMatrix());
1126  double dX1 = sqrt(mat1[3][3]);
1127  double dY1 = sqrt(mat1[4][4]);
1128  AlgebraicSymMatrix mat2(5, 0);
1129  mat2 = (SegmentRB[1].first)->parametersError().similarityT((SegmentRB[1].first)->projectionMatrix());
1130  double dX2 = sqrt(mat2[3][3]);
1131  double dY2 = sqrt(mat2[4][4]);
1132  AlgebraicSymMatrix mat3(5, 0);
1133  mat3 = (SegmentRB[1].second)->parametersError().similarityT((SegmentRB[1].second)->projectionMatrix());
1134  double dX3 = sqrt(mat3[3][3]);
1135  double dY3 = sqrt(mat3[4][4]);
1136  cout << "Local error for 4 recHits are: " << dX0 << ", " << dY0 << ", " << dX1 << ", " << dY1 << ", " << dX2 << ", "
1137  << dY2 << ", " << dX3 << ", " << dY3 << endl;
1138  const GeomDetUnit* refRPC1 = (SegmentRB[0].second)->detUnit();
1139  LocalPoint recHit0 = refRPC1->toLocal((SegmentRB[0].first)->globalPosition());
1140  LocalPoint recHit1 = refRPC1->toLocal((SegmentRB[0].second)->globalPosition());
1141  LocalVector localSegment00 = (LocalVector)(recHit1 - recHit0);
1142  LocalVector localSegment01 = LocalVector(localSegment00.x() + dX0 + dX1, localSegment00.y(), localSegment00.z());
1143  LocalVector localSegment02 = LocalVector(localSegment00.x() - dX0 - dX1, localSegment00.y(), localSegment00.z());
1144  GlobalVector globalSegment00 = refRPC1->toGlobal(localSegment00);
1145  GlobalVector globalSegment01 = refRPC1->toGlobal(localSegment01);
1146  GlobalVector globalSegment02 = refRPC1->toGlobal(localSegment02);
1147 
1148  const GeomDetUnit* refRPC2 = (SegmentRB[1].first)->detUnit();
1149  LocalPoint recHit2 = refRPC2->toLocal((SegmentRB[1].first)->globalPosition());
1150  LocalPoint recHit3 = refRPC2->toLocal((SegmentRB[1].second)->globalPosition());
1151  LocalVector localSegment10 = (LocalVector)(recHit3 - recHit2);
1152  LocalVector localSegment11 = LocalVector(localSegment10.x() + dX2 + dX3, localSegment10.y(), localSegment10.z());
1153  LocalVector localSegment12 = LocalVector(localSegment10.x() - dX2 - dX3, localSegment10.y(), localSegment10.z());
1154  GlobalVector globalSegment10 = refRPC2->toGlobal(localSegment10);
1155  GlobalVector globalSegment11 = refRPC2->toGlobal(localSegment11);
1156  GlobalVector globalSegment12 = refRPC2->toGlobal(localSegment12);
1157 
1158  if (isClockwise != 0) {
1159  GlobalVector vec[2];
1160  vec[0] = GlobalVector(
1161  (entryPosition.x() - Center2.x()), (entryPosition.y() - Center2.y()), (entryPosition.z() - Center2.z()));
1162  vec[1] = GlobalVector(
1163  (leavePosition.x() - Center2.x()), (leavePosition.y() - Center2.y()), (leavePosition.z() - Center2.z()));
1164  double halfPhiCenter = fabs((vec[1].phi() - vec[0].phi()).value()) / 2;
1165  // dPhi0 shoule be the same clockwise direction, while dPhi1 should be opposite clockwise direction, w.r.t to the track clockwise
1166  double dPhi0 = (((globalSegment00.phi() - globalSegment01.phi()).value() * isClockwise) > 0)
1167  ? fabs((globalSegment00.phi() - globalSegment01.phi()).value())
1168  : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1169  double dPhi1 = (((globalSegment10.phi() - globalSegment11.phi()).value() * isClockwise) < 0)
1170  ? fabs((globalSegment10.phi() - globalSegment11.phi()).value())
1171  : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1172  // For the deltaR should be kept small, we assume the delta Phi0/Phi1 should be in a same limit value
1173  double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1174  cout << "DPhi for new Segment is about " << dPhi << endl;
1175  // Check the variance of halfPhiCenter
1176  double newhalfPhiCenter = ((halfPhiCenter - dPhi) > 0 ? (halfPhiCenter - dPhi) : 0);
1177  if (newhalfPhiCenter != 0) {
1178  double newmeanPt = meanPt * halfPhiCenter / newhalfPhiCenter;
1179  if (fabs(newmeanPt) > upper_limit_pt)
1180  newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1181  cout << "The error is inside range. Max meanPt could be " << newmeanPt << endl;
1182  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1183  } else {
1184  double newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1185  cout << "The error is outside range. Max meanPt could be " << newmeanPt << endl;
1186  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1187  }
1188  } else {
1189  double dPhi0 = (fabs((globalSegment00.phi() - globalSegment01.phi()).value()) <=
1190  fabs((globalSegment00.phi() - globalSegment02.phi()).value()))
1191  ? fabs((globalSegment00.phi() - globalSegment01.phi()).value())
1192  : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1193  double dPhi1 = (fabs((globalSegment10.phi() - globalSegment11.phi()).value()) <=
1194  fabs((globalSegment10.phi() - globalSegment12.phi()).value()))
1195  ? fabs((globalSegment10.phi() - globalSegment11.phi()).value())
1196  : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1197  double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1198  GlobalVector middleSegment = leavePosition - entryPosition;
1199  double halfDistance = middleSegment.perp() / 2;
1200  double newmeanPt = halfDistance / dPhi;
1201  cout << "The error is for straight. Max meanPt could be " << newmeanPt << endl;
1202  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1203  }
1204 
1205  double dXdZ1 = globalSegment11.x() / globalSegment11.z() - globalSegment10.x() / globalSegment10.z();
1206  double dXdZ2 = globalSegment12.x() / globalSegment12.z() - globalSegment10.x() / globalSegment10.z();
1207  dXdZ = (fabs(dXdZ1) >= fabs(dXdZ2)) ? dXdZ1 : dXdZ2;
1208 
1209  LocalVector localSegment13 = LocalVector(localSegment10.x(), localSegment10.y() + dY2 + dY3, localSegment10.z());
1210  LocalVector localSegment14 = LocalVector(localSegment10.x(), localSegment10.y() - dY2 - dY3, localSegment10.z());
1211  GlobalVector globalSegment13 = refRPC2->toGlobal(localSegment13);
1212  GlobalVector globalSegment14 = refRPC2->toGlobal(localSegment14);
1213  double dYdZ1 = globalSegment13.y() / globalSegment13.z() - globalSegment10.y() / globalSegment10.z();
1214  double dYdZ2 = globalSegment14.y() / globalSegment14.z() - globalSegment10.y() / globalSegment10.z();
1215  dYdZ = (fabs(dYdZ1) >= fabs(dYdZ2)) ? dYdZ1 : dYdZ2;
1216 
1217  mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1218  mat[1][1] = dXdZ * dXdZ;
1219  mat[2][2] = dYdZ * dYdZ;
1220  Error = LocalTrajectoryError(asSMatrix<5>(mat));
1221  }
1222  return Error;
1223 }
edm::ErrorSummaryEntry Error
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:69
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:58
GlobalPoint entryPosition
T z() const
Definition: PV3DBase.h:61
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
GlobalVector Momentum
RPCSegment SegmentRB[2]
U second(std::pair< T, U > const &p)
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
T mag() const
Definition: PV3DBase.h:64
GlobalVector Center
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
GlobalVector Center2
#define N
Definition: blowfish.cc:9
unsigned int AlgorithmType
unsigned int nrhit() const
CLHEP::HepSymMatrix AlgebraicSymMatrix
GlobalPoint leavePosition
#define upper_limit_pt
Global3DVector GlobalVector
Definition: GlobalVector.h:10

◆ MiddlePointsAlgorithm()

void RPCSeedPattern::MiddlePointsAlgorithm ( )
private

Definition at line 170 of file RPCSeedPattern.cc.

References gather_cfg::cout, mps_fire::i, PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, dqmiodumpmetadata::n, DiDispStaMuonMonitor_cfi::pt, X, x, and beamSpotPI::Y.

170  {
171  cout << "Using middle points algorithm" << endl;
172  unsigned int NumberofHitsinSeed = nrhit();
173  // Check recHit number, if fail we set the pattern to "wrong"
174  if (NumberofHitsinSeed < 4) {
175  isPatternChecked = true;
176  isGoodPattern = -1;
177  return;
178  }
179  double* X = new double[NumberofHitsinSeed];
180  double* Y = new double[NumberofHitsinSeed];
181  unsigned int n = 0;
182  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
183  X[n] = (*iter)->globalPosition().x();
184  Y[n] = (*iter)->globalPosition().y();
185  cout << "X[" << n << "] = " << X[n] << ", Y[" << n << "]= " << Y[n] << endl;
186  n++;
187  }
188  unsigned int NumberofPoints = NumberofHitsinSeed;
189  while (NumberofPoints > 3) {
190  for (unsigned int i = 0; i <= (NumberofPoints - 2); i++) {
191  X[i] = (X[i] + X[i + 1]) / 2;
192  Y[i] = (Y[i] + Y[i + 1]) / 2;
193  }
194  NumberofPoints--;
195  }
196  double x[3], y[3];
197  for (unsigned int i = 0; i < 3; i++) {
198  x[i] = X[i];
199  y[i] = Y[i];
200  }
201  double pt = 0;
202  double pt_err = 0;
203  bool checkStraight = checkStraightwithThreerecHits(x, y, MinDeltaPhi);
204  if (!checkStraight) {
205  GlobalVector Center_temp = computePtWithThreerecHits(pt, pt_err, x, y);
206  double meanR = 0.;
207  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
208  meanR += getDistance(*iter, Center_temp);
209  meanR /= NumberofHitsinSeed;
210  // For simple pattern
211  isStraight = false;
212  Center = Center_temp;
213  meanRadius = meanR;
214  } else {
215  // For simple pattern
216  isStraight = true;
217  meanRadius = -1;
218  }
219 
220  // Unset the pattern estimation signa
221  isPatternChecked = false;
222 
223  delete[] X;
224  delete[] Y;
225 }
#define X(str)
Definition: MuonsGrabber.cc:38
GlobalVector computePtWithThreerecHits(double &pt, double &pt_err, double(&x)[3], double(&y)[3]) const
ConstMuonRecHitContainer theRecHits
bool checkStraightwithThreerecHits(ConstMuonRecHitPointer(&precHit)[3], double MinDeltaPhi) const
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
GlobalVector Center
unsigned int nrhit() const

◆ nrhit()

unsigned int RPCSeedPattern::nrhit ( ) const
inline

Definition at line 46 of file RPCSeedPattern.h.

References theRecHits.

46 { return theRecHits.size(); }
ConstMuonRecHitContainer theRecHits

◆ seed()

RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::seed ( const MagneticField Field,
const RPCGeometry rpcGeom,
int &  isGoodSeed 
)
private

Definition at line 51 of file RPCSeedPattern.cc.

References gather_cfg::cout.

53  {
54  if (isConfigured == false) {
55  cout << "Configuration not set yet" << endl;
56  return createFakeSeed(isGoodSeed, Field);
57  }
58 
59  // Check recHit number, if fail we return a fake seed and set pattern to "wrong"
60  unsigned int NumberofHitsinSeed = nrhit();
61  if (NumberofHitsinSeed < 3) {
62  return createFakeSeed(isGoodSeed, Field);
63  }
64  // If only three recHits, we don't have other choice
65  if (NumberofHitsinSeed == 3)
67 
68  if (NumberofHitsinSeed > 3) {
69  if (autoAlgorithmChoose == false) {
70  cout << "computePtWithmorerecHits" << endl;
71  if (AlgorithmType == 0)
73  if (AlgorithmType == 1)
75  if (AlgorithmType == 2)
77  if (AlgorithmType == 3) {
78  if (checkSegment()) {
80  } else {
81  cout << "Not enough recHits for Special Segment Algorithm" << endl;
82  return createFakeSeed(isGoodSeed, Field);
83  }
84  }
85  } else {
86  if (checkSegment()) {
87  AlgorithmType = 3;
89  } else {
90  AlgorithmType = 2;
92  }
93  }
94  }
95 
96  // Check the pattern
97  if (isPatternChecked == false) {
98  if (AlgorithmType != 3) {
99  checkSimplePattern(Field);
100  } else {
101  checkSegmentAlgorithmSpecial(Field, rpcGeom);
102  }
103  }
104 
105  return createSeed(isGoodSeed, Field);
106 }
weightedTrajectorySeed createSeed(int &isGoodSeed, const MagneticField &)
void MiddlePointsAlgorithm()
bool autoAlgorithmChoose
bool checkSegment() const
void ThreePointsAlgorithm()
void SegmentAlgorithmSpecial(const MagneticField &Field)
void checkSegmentAlgorithmSpecial(const MagneticField &Field, const RPCGeometry &rpcGeometry)
unsigned int AlgorithmType
unsigned int nrhit() const
void checkSimplePattern(const MagneticField &Field)
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const MagneticField &)

◆ SegmentAlgorithm()

void RPCSeedPattern::SegmentAlgorithm ( )
private

Definition at line 227 of file RPCSeedPattern.cc.

References gather_cfg::cout, mps_fire::i, PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, and dqmiodumpmetadata::n.

227  {
228  cout << "Using segments algorithm" << endl;
229  unsigned int NumberofHitsinSeed = nrhit();
230  // Check recHit number, if fail we set the pattern to "wrong"
231  if (NumberofHitsinSeed < 4) {
232  isPatternChecked = true;
233  isGoodPattern = -1;
234  return;
235  }
236 
237  RPCSegment* Segment;
238  unsigned int NumberofSegment = NumberofHitsinSeed - 2;
239  Segment = new RPCSegment[NumberofSegment];
240  unsigned int n = 0;
241  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 2); iter++) {
242  Segment[n].first = (*iter);
243  Segment[n].second = (*(iter + 2));
244  n++;
245  }
246  unsigned int NumberofStraight = 0;
247  for (unsigned int i = 0; i < NumberofSegment - 1; i++) {
248  bool checkStraight = checkStraightwithSegment(Segment[i], Segment[i + 1], MinDeltaPhi);
249  if (checkStraight == true) {
250  // For simple patterm
251  NumberofStraight++;
252  } else {
253  GlobalVector Center_temp = computePtwithSegment(Segment[i], Segment[i + 1]);
254  // For simple patterm
255  Center += Center_temp;
256  }
257  }
258  // For simple pattern, only one general parameter for pattern
259  if ((NumberofSegment - 1 - NumberofStraight) > 0) {
260  isStraight = false;
261  Center /= (NumberofSegment - 1 - NumberofStraight);
262  double meanR = 0.;
263  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
264  meanR += getDistance(*iter, Center);
265  meanR /= NumberofHitsinSeed;
266  meanRadius = meanR;
267  } else {
268  isStraight = true;
269  meanRadius = -1;
270  }
271 
272  // Unset the pattern estimation signal
273  isPatternChecked = false;
274 
275  delete[] Segment;
276 }
std::pair< ConstMuonRecHitPointer, ConstMuonRecHitPointer > RPCSegment
GlobalVector computePtwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2) const
ConstMuonRecHitContainer theRecHits
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
GlobalVector Center
unsigned int nrhit() const
bool checkStraightwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2, double MinDeltaPhi) const

◆ SegmentAlgorithmSpecial()

void RPCSeedPattern::SegmentAlgorithmSpecial ( const MagneticField Field)
private

Definition at line 278 of file RPCSeedPattern.cc.

References volumeBasedMagneticField_160812_cfi::BValue, gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), SiPixelRawToDigiRegional_cfi::deltaPhi, PbPb_ZMuSkimMuonDPG_cff::deltaR, PVValHelper::dx, PVValHelper::dy, PVValHelper::dz, MagneticField::inTesla(), PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, dqmiodumpmetadata::n, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), mathSSE::sqrt(), Vector3DBase< T, FrameTag >::unit(), relativeConstraints::value, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

278  {
279  //unsigned int NumberofHitsinSeed = nrhit();
280  if (!checkSegment()) {
281  isPatternChecked = true;
282  isGoodPattern = -1;
283  return;
284  }
285 
286  // Get magnetice field sampling information, recHit's position is not the border of Chamber and Iron
287  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
288  GlobalPoint gpFirst = (*iter)->globalPosition();
289  GlobalPoint gpLast = (*(iter + 1))->globalPosition();
291  double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
292  double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
293  double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
294  for (unsigned int index = 0; index < sampleCount; index++) {
295  gp[index] = GlobalPoint(
296  (gpFirst.x() + dx * (index + 1)), (gpFirst.y() + dy * (index + 1)), (gpFirst.z() + dz * (index + 1)));
297  GlobalVector MagneticVec_temp = Field.inTesla(gp[index]);
298  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
299  //BValue.push_back(MagneticVec_temp);
300  }
301  delete[] gp;
302  }
303 
304  // form two segments
305  ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin();
306  for (unsigned int n = 0; n <= 1; n++) {
307  SegmentRB[n].first = (*iter);
308  cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
309  iter++;
310  SegmentRB[n].second = (*iter);
311  cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
312  iter++;
313  }
314  GlobalVector segvec1 = (SegmentRB[0].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
315  GlobalVector segvec2 = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
316 
317  // extrapolate the segment to find the Iron border which magnetic field is at large value
318  entryPosition = (SegmentRB[0].second)->globalPosition();
319  leavePosition = (SegmentRB[1].first)->globalPosition();
320  while (fabs(Field.inTesla(entryPosition).z()) < MagnecticFieldThreshold) {
321  cout << "Entry position is : " << entryPosition << ", and stepping into next point" << endl;
322  entryPosition += segvec1.unit() * stepLength;
323  }
324  // Loop back for more accurate by stepLength/10
325  while (fabs(Field.inTesla(entryPosition).z()) >= MagnecticFieldThreshold) {
326  cout << "Entry position is : " << entryPosition << ", and stepping back into next point" << endl;
327  entryPosition -= segvec1.unit() * stepLength / 10;
328  }
329  entryPosition += 0.5 * segvec1.unit() * stepLength / 10;
330  cout << "Final entry position is : " << entryPosition << endl;
331 
332  while (fabs(Field.inTesla(leavePosition).z()) < MagnecticFieldThreshold) {
333  cout << "Leave position is : " << leavePosition << ", and stepping into next point" << endl;
334  leavePosition -= segvec2.unit() * stepLength;
335  }
336  // Loop back for more accurate by stepLength/10
337  while (fabs(Field.inTesla(leavePosition).z()) >= MagnecticFieldThreshold) {
338  cout << "Leave position is : " << leavePosition << ", and stepping back into next point" << endl;
339  leavePosition += segvec2.unit() * stepLength / 10;
340  }
341  leavePosition -= 0.5 * segvec2.unit() * stepLength / 10;
342  cout << "Final leave position is : " << leavePosition << endl;
343 
344  // Sampling magnetic field in Iron region
346  double dx = (leavePosition.x() - entryPosition.x()) / (sampleCount + 1);
347  double dy = (leavePosition.y() - entryPosition.y()) / (sampleCount + 1);
348  double dz = (leavePosition.z() - entryPosition.z()) / (sampleCount + 1);
349  std::vector<GlobalVector> BValue;
350  BValue.clear();
351  for (unsigned int index = 0; index < sampleCount; index++) {
352  gp[index] = GlobalPoint((entryPosition.x() + dx * (index + 1)),
353  (entryPosition.y() + dy * (index + 1)),
354  (entryPosition.z() + dz * (index + 1)));
355  GlobalVector MagneticVec_temp = Field.inTesla(gp[index]);
356  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
357  BValue.push_back(MagneticVec_temp);
358  }
359  delete[] gp;
360  GlobalVector meanB2(0, 0, 0);
361  for (std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
362  meanB2 += (*BIter);
363  meanB2 /= BValue.size();
364  cout << "Mean B field is " << meanB2 << endl;
365  meanMagneticField2 = meanB2;
366 
367  double meanBz2 = meanB2.z();
368  double deltaBz2 = 0.;
369  for (std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
370  deltaBz2 += (BIter->z() - meanBz2) * (BIter->z() - meanBz2);
371  deltaBz2 /= BValue.size();
372  deltaBz2 = sqrt(deltaBz2);
373  cout << "delta Bz is " << deltaBz2 << endl;
374 
375  // Distance of the initial 3 segment
376  S = 0;
377  bool checkStraight = checkStraightwithSegment(SegmentRB[0], SegmentRB[1], MinDeltaPhi);
378  if (checkStraight == true) {
379  // Just for complex pattern
380  isStraight2 = checkStraight;
381  Center2 = GlobalVector(0, 0, 0);
382  meanRadius2 = -1;
383  GlobalVector MomentumVec = (SegmentRB[1].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
384  S += MomentumVec.perp();
385  lastPhi = MomentumVec.phi().value();
386  } else {
387  GlobalVector seg1 = entryPosition - (SegmentRB[0].first)->globalPosition();
388  S += seg1.perp();
389  GlobalVector seg2 = (SegmentRB[1].second)->globalPosition() - leavePosition;
390  S += seg2.perp();
391  GlobalVector vecZ(0, 0, 1);
392  GlobalVector gvec1 = seg1.cross(vecZ);
393  GlobalVector gvec2 = seg2.cross(vecZ);
394  double A1 = gvec1.x();
395  double B1 = gvec1.y();
396  double A2 = gvec2.x();
397  double B2 = gvec2.y();
398  double X1 = entryPosition.x();
399  double Y1 = entryPosition.y();
400  double X2 = leavePosition.x();
401  double Y2 = leavePosition.y();
402  double XO = (A1 * A2 * (Y2 - Y1) + A2 * B1 * X1 - A1 * B2 * X2) / (A2 * B1 - A1 * B2);
403  double YO = (B1 * B2 * (X2 - X1) + B2 * A1 * Y1 - B1 * A2 * Y2) / (B2 * A1 - B1 * A2);
404  GlobalVector Center_temp(XO, YO, 0);
405  // Just for complex pattern
406  isStraight2 = checkStraight;
407  Center2 = Center_temp;
408 
409  cout << "entryPosition: " << entryPosition << endl;
410  cout << "leavePosition: " << leavePosition << endl;
411  cout << "Center2 is : " << Center_temp << endl;
412 
413  double R1 = GlobalVector((entryPosition.x() - Center_temp.x()),
414  (entryPosition.y() - Center_temp.y()),
415  (entryPosition.z() - Center_temp.z()))
416  .perp();
417  double R2 = GlobalVector((leavePosition.x() - Center_temp.x()),
418  (leavePosition.y() - Center_temp.y()),
419  (leavePosition.z() - Center_temp.z()))
420  .perp();
421  double meanR = (R1 + R2) / 2;
422  double deltaR = sqrt(((R1 - meanR) * (R1 - meanR) + (R2 - meanR) * (R2 - meanR)) / 2);
423  meanRadius2 = meanR;
424  cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
425  cout << "Mean radius is " << meanR << endl;
426  cout << "Delta R is " << deltaR << endl;
427  double deltaPhi =
428  fabs(((leavePosition - GlobalPoint(XO, YO, 0)).phi() - (entryPosition - GlobalPoint(XO, YO, 0)).phi()).value());
429  S += meanR * deltaPhi;
430  lastPhi = seg2.phi().value();
431  }
432 
433  // Unset the pattern estimation signa
434  isPatternChecked = false;
435 }
T perp() const
Definition: PV3DBase.h:69
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:110
double MagnecticFieldThreshold
GlobalPoint entryPosition
T z() const
Definition: PV3DBase.h:61
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
RPCSegment SegmentRB[2]
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
bool checkSegment() const
ConstMuonRecHitContainer theRecHits
GlobalVector Center2
GlobalVector meanMagneticField2
bool checkStraightwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2, double MinDeltaPhi) const
GlobalPoint leavePosition
Vector3DBase unit() const
Definition: Vector3DBase.h:54
unsigned int sampleCount
Global3DVector GlobalVector
Definition: GlobalVector.h:10

◆ ThreePointsAlgorithm()

void RPCSeedPattern::ThreePointsAlgorithm ( )
private

Definition at line 108 of file RPCSeedPattern.cc.

References gather_cfg::cout, mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, PythiaFilterGammaJetWithBg_cfi::MinDeltaPhi, dqmiodumpmetadata::n, DiDispStaMuonMonitor_cfi::pt, and upper_limit_pt.

108  {
109  cout << "computePtWith3recHits" << endl;
110  unsigned int NumberofHitsinSeed = nrhit();
111  // Check recHit number, if fail we set the pattern to "wrong"
112  if (NumberofHitsinSeed < 3) {
113  isPatternChecked = true;
114  isGoodPattern = -1;
115  return;
116  }
117  // Choose every 3 recHits to form a part
118  unsigned int NumberofPart = NumberofHitsinSeed * (NumberofHitsinSeed - 1) * (NumberofHitsinSeed - 2) / (3 * 2);
119  ;
120  double* pt = new double[NumberofPart];
121  double* pt_err = new double[NumberofPart];
122  // Loop for each three-recHits part
123  ConstMuonRecHitPointer precHit[3];
124  unsigned int n = 0;
125  unsigned int NumberofStraight = 0;
126  for (unsigned int i = 0; i < (NumberofHitsinSeed - 2); i++)
127  for (unsigned int j = (i + 1); j < (NumberofHitsinSeed - 1); j++)
128  for (unsigned int k = (j + 1); k < NumberofHitsinSeed; k++) {
129  precHit[0] = theRecHits[i];
130  precHit[1] = theRecHits[j];
131  precHit[2] = theRecHits[k];
132  bool checkStraight = checkStraightwithThreerecHits(precHit, MinDeltaPhi);
133  if (!checkStraight) {
134  GlobalVector Center_temp = computePtwithThreerecHits(pt[n], pt_err[n], precHit);
135  // For simple pattern
136  Center += Center_temp;
137  } else {
138  // For simple pattern
139  NumberofStraight++;
140  pt[n] = upper_limit_pt;
141  pt_err[n] = 0;
142  }
143  n++;
144  }
145  // For simple pattern, only one general parameter for pattern
146  if (NumberofStraight == NumberofPart) {
147  isStraight = true;
148  meanRadius = -1;
149  } else {
150  isStraight = false;
151  Center /= (NumberofPart - NumberofStraight);
152  double meanR = 0.;
153  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
154  meanR += getDistance(*iter, Center);
155  meanR /= NumberofHitsinSeed;
156  meanRadius = meanR;
157  }
158 
159  // Unset the pattern estimation signa
160  isPatternChecked = false;
161 
162  //double ptmean0 = 0;
163  //double sptmean0 = 0;
164  //computeBestPt(pt, pt_err, ptmean0, sptmean0, (NumberofPart - NumberofStraight));
165 
166  delete[] pt;
167  delete[] pt_err;
168 }
ConstMuonRecHitContainer theRecHits
bool checkStraightwithThreerecHits(ConstMuonRecHitPointer(&precHit)[3], double MinDeltaPhi) const
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
GlobalVector Center
unsigned int nrhit() const
GlobalVector computePtwithThreerecHits(double &pt, double &pt_err, ConstMuonRecHitPointer(&precHit)[3]) const
#define upper_limit_pt
MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer

Friends And Related Function Documentation

◆ RPCSeedFinder

friend class RPCSeedFinder
friend

Definition at line 49 of file RPCSeedPattern.h.

Member Data Documentation

◆ AlgorithmType

unsigned int RPCSeedPattern::AlgorithmType
private

Definition at line 87 of file RPCSeedPattern.h.

◆ autoAlgorithmChoose

bool RPCSeedPattern::autoAlgorithmChoose
private

Definition at line 88 of file RPCSeedPattern.h.

◆ Center

GlobalVector RPCSeedPattern::Center
private

Definition at line 110 of file RPCSeedPattern.h.

◆ Center2

GlobalVector RPCSeedPattern::Center2
private

Definition at line 101 of file RPCSeedPattern.h.

◆ Charge

int RPCSeedPattern::Charge
private

Definition at line 119 of file RPCSeedPattern.h.

◆ deltaBz

double RPCSeedPattern::deltaBz
private

Definition at line 113 of file RPCSeedPattern.h.

◆ deltaRThreshold

double RPCSeedPattern::deltaRThreshold
private

Definition at line 86 of file RPCSeedPattern.h.

◆ entryPosition

GlobalPoint RPCSeedPattern::entryPosition
private

Definition at line 104 of file RPCSeedPattern.h.

◆ isClockwise

int RPCSeedPattern::isClockwise
private

Definition at line 117 of file RPCSeedPattern.h.

◆ isConfigured

bool RPCSeedPattern::isConfigured
private

Definition at line 94 of file RPCSeedPattern.h.

◆ isGoodPattern

int RPCSeedPattern::isGoodPattern
private

Definition at line 116 of file RPCSeedPattern.h.

◆ isParralZ

int RPCSeedPattern::isParralZ
private

Definition at line 118 of file RPCSeedPattern.h.

◆ isPatternChecked

bool RPCSeedPattern::isPatternChecked
private

Definition at line 115 of file RPCSeedPattern.h.

◆ isStraight

bool RPCSeedPattern::isStraight
private

Definition at line 109 of file RPCSeedPattern.h.

◆ isStraight2

bool RPCSeedPattern::isStraight2
private

Definition at line 100 of file RPCSeedPattern.h.

◆ lastPhi

double RPCSeedPattern::lastPhi
private

Definition at line 106 of file RPCSeedPattern.h.

◆ leavePosition

GlobalPoint RPCSeedPattern::leavePosition
private

Definition at line 105 of file RPCSeedPattern.h.

◆ MagnecticFieldThreshold

double RPCSeedPattern::MagnecticFieldThreshold
private

Definition at line 98 of file RPCSeedPattern.h.

◆ MaxRSD

double RPCSeedPattern::MaxRSD
private

Definition at line 85 of file RPCSeedPattern.h.

◆ meanBz

double RPCSeedPattern::meanBz
private

Definition at line 112 of file RPCSeedPattern.h.

◆ meanMagneticField2

GlobalVector RPCSeedPattern::meanMagneticField2
private

Definition at line 99 of file RPCSeedPattern.h.

◆ meanPt

double RPCSeedPattern::meanPt
private

Definition at line 120 of file RPCSeedPattern.h.

◆ meanRadius

double RPCSeedPattern::meanRadius
private

Definition at line 111 of file RPCSeedPattern.h.

◆ meanRadius2

double RPCSeedPattern::meanRadius2
private

Definition at line 102 of file RPCSeedPattern.h.

◆ meanSpt

double RPCSeedPattern::meanSpt
private

Definition at line 121 of file RPCSeedPattern.h.

◆ MinDeltaPhi

double RPCSeedPattern::MinDeltaPhi
private

Definition at line 90 of file RPCSeedPattern.h.

◆ Momentum

GlobalVector RPCSeedPattern::Momentum
private

Definition at line 122 of file RPCSeedPattern.h.

◆ S

double RPCSeedPattern::S
private

Definition at line 107 of file RPCSeedPattern.h.

◆ sampleCount

unsigned int RPCSeedPattern::sampleCount
private

Definition at line 92 of file RPCSeedPattern.h.

◆ SegmentRB

RPCSegment RPCSeedPattern::SegmentRB[2]
private

Definition at line 103 of file RPCSeedPattern.h.

◆ stepLength

double RPCSeedPattern::stepLength
private

Definition at line 91 of file RPCSeedPattern.h.

◆ theRecHits

ConstMuonRecHitContainer RPCSeedPattern::theRecHits
private

Definition at line 96 of file RPCSeedPattern.h.

Referenced by add(), clear(), and nrhit().

◆ ZError

double RPCSeedPattern::ZError
private

Definition at line 89 of file RPCSeedPattern.h.