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 edm::EventSetup &eSetup)
 
void checkSimplePattern (const edm::EventSetup &eSetup)
 
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 edm::EventSetup &eSetup)
 
weightedTrajectorySeed createSeed (int &isGoodSeed, const edm::EventSetup &eSetup)
 
double extropolateStep (const GlobalPoint &startPosition, const GlobalVector &startMomentum, ConstMuonRecHitContainer::const_iterator iter, const int ClockwiseDirection, double &tracklength, const edm::EventSetup &eSetup)
 
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 edm::EventSetup &eSetup, int &isGoodSeed)
 
void SegmentAlgorithm ()
 
void SegmentAlgorithmSpecial (const edm::EventSetup &eSetup)
 
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 28 of file RPCSeedPattern.h.

Member Typedef Documentation

Definition at line 32 of file RPCSeedPattern.h.

Definition at line 30 of file RPCSeedPattern.h.

Definition at line 31 of file RPCSeedPattern.h.

Definition at line 29 of file RPCSeedPattern.h.

Definition at line 35 of file RPCSeedPattern.h.

Definition at line 36 of file RPCSeedPattern.h.

Constructor & Destructor Documentation

RPCSeedPattern::RPCSeedPattern ( )

Definition at line 32 of file RPCSeedPattern.cc.

32  {
33  isPatternChecked = false;
34  isConfigured = false;
36 }
double MagnecticFieldThreshold
RPCSeedPattern::~RPCSeedPattern ( )

Definition at line 38 of file RPCSeedPattern.cc.

38 {}

Member Function Documentation

void RPCSeedPattern::add ( const ConstMuonRecHitPointer hit)
inline

Definition at line 43 of file RPCSeedPattern.h.

References theRecHits.

Referenced by counter.Counter::register().

43 { theRecHits.push_back(hit); }
ConstMuonRecHitContainer theRecHits
MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::BestRefRecHit ( ) const
private

Definition at line 470 of file RPCSeedPattern.cc.

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

Definition at line 440 of file RPCSeedPattern.cc.

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

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

Definition at line 748 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.

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

Definition at line 607 of file RPCSeedPattern.cc.

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

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

Definition at line 528 of file RPCSeedPattern.cc.

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

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

Definition at line 490 of file RPCSeedPattern.cc.

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

490  {
491  GlobalVector segvec1 = precHit[1]->globalPosition() - precHit[0]->globalPosition();
492  GlobalVector segvec2 = precHit[2]->globalPosition() - precHit[1]->globalPosition();
493  double dPhi = (segvec2.phi() - segvec1.phi()).value();
494  if (fabs(dPhi) > MinDeltaPhi) {
495  cout << "Part is estimate to be not straight" << endl;
496  return false;
497  } else {
498  cout << "Part is estimate to be straight" << endl;
499  return true;
500  }
501 }
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
bool RPCSeedPattern::checkStraightwithThreerecHits ( double(&)  x[3],
double(&)  y[3],
double  MinDeltaPhi 
) const
private

Definition at line 575 of file RPCSeedPattern.cc.

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

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

Definition at line 42 of file RPCSeedPattern.h.

References theRecHits.

42 { theRecHits.clear(); }
ConstMuonRecHitContainer theRecHits
GlobalVector RPCSeedPattern::computePtwithSegment ( const RPCSegment Segment1,
const RPCSegment Segment2 
) const
private

Definition at line 549 of file RPCSeedPattern.cc.

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

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

Definition at line 503 of file RPCSeedPattern.cc.

References MaterialEffects_cfi::A, gather_cfg::cout, and mathSSE::sqrt().

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

Definition at line 588 of file RPCSeedPattern.cc.

References MaterialEffects_cfi::A, gather_cfg::cout, and mathSSE::sqrt().

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

Definition at line 40 of file RPCSeedPattern.cc.

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

40  {
41  MaxRSD = iConfig.getParameter<double>("MaxRSD");
42  deltaRThreshold = iConfig.getParameter<double>("deltaRThreshold");
43  AlgorithmType = iConfig.getParameter<unsigned int>("AlgorithmType");
44  autoAlgorithmChoose = iConfig.getParameter<bool>("autoAlgorithmChoose");
45  ZError = iConfig.getParameter<double>("ZError");
46  MinDeltaPhi = iConfig.getParameter<double>("MinDeltaPhi");
47  MagnecticFieldThreshold = iConfig.getParameter<double>("MagnecticFieldThreshold");
48  stepLength = iConfig.getParameter<double>("stepLength");
49  sampleCount = iConfig.getParameter<unsigned int>("sampleCount");
50  isConfigured = true;
51 }
T getParameter(std::string const &) const
double MagnecticFieldThreshold
bool autoAlgorithmChoose
double deltaRThreshold
unsigned int AlgorithmType
unsigned int sampleCount
RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createFakeSeed ( int &  isGoodSeed,
const edm::EventSetup eSetup 
)
private

Definition at line 967 of file RPCSeedPattern.cc.

References alongMomentum, gather_cfg::cout, relativeConstraints::error, edm::EventSetup::get(), trajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), and upper_limit_pt.

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

Definition at line 1017 of file RPCSeedPattern.cc.

References alongMomentum, gather_cfg::cout, cross(), debug, SiPixelRawToDigiRegional_cfi::deltaPhi, MuonPatternRecoDumper::dumpMuonId(), MuonPatternRecoDumper::dumpTSOS(), relativeConstraints::error, dqmdumpme::first, edm::EventSetup::get(), createfilelist::int, trajectoryStateTransform::persistentState(), PV3DBase< T, PVType, FrameType >::phi(), edm::OwnVector< T, P >::push_back(), Vector3DBase< T, FrameTag >::unit(), and relativeConstraints::value.

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

Definition at line 882 of file RPCSeedPattern.cc.

References RPCGeometry::chamber(), gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), SiPixelRawToDigiRegional_cfi::deltaPhi, edm::EventSetup::get(), MagneticField::inTesla(), PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), VtxSmearedParameters_cfi::Phi, JetPartonCorrections_cff::Radius, DetId::rawId(), GeomDet::surface(), Vector3DBase< T, FrameTag >::unit(), and PV3DBase< T, PVType, FrameType >::z().

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

Definition at line 468 of file RPCSeedPattern.cc.

468 { return theRecHits.front(); }
ConstMuonRecHitContainer theRecHits
double RPCSeedPattern::getDistance ( const ConstMuonRecHitPointer precHit,
const GlobalVector Center 
) const
private

Definition at line 485 of file RPCSeedPattern.cc.

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

485  {
486  return sqrt((precHit->globalPosition().x() - Center.x()) * (precHit->globalPosition().x() - Center.x()) +
487  (precHit->globalPosition().y() - Center.y()) * (precHit->globalPosition().y() - Center.y()));
488 }
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
T x() const
Definition: PV3DBase.h:59
LocalTrajectoryError RPCSeedPattern::getSpecialAlgorithmErrorMatrix ( const ConstMuonRecHitPointer first,
const ConstMuonRecHitPointer best 
)
private

Definition at line 1108 of file RPCSeedPattern.cc.

References zMuMuMuonUserData::beta, gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, HLT_2018_cff::dPhi, 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().

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

Definition at line 169 of file RPCSeedPattern.cc.

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

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

Definition at line 44 of file RPCSeedPattern.h.

References theRecHits.

44 { return theRecHits.size(); }
ConstMuonRecHitContainer theRecHits
RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::seed ( const edm::EventSetup eSetup,
int &  isGoodSeed 
)
private

Definition at line 53 of file RPCSeedPattern.cc.

References gather_cfg::cout.

53  {
54  if (isConfigured == false) {
55  cout << "Configuration not set yet" << endl;
56  return createFakeSeed(isGoodSeed, eSetup);
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, eSetup);
63  // If only three recHits, we don't have other choice
64  if (NumberofHitsinSeed == 3)
66 
67  if (NumberofHitsinSeed > 3) {
68  if (autoAlgorithmChoose == false) {
69  cout << "computePtWithmorerecHits" << endl;
70  if (AlgorithmType == 0)
72  if (AlgorithmType == 1)
74  if (AlgorithmType == 2)
76  if (AlgorithmType == 3) {
77  if (checkSegment())
79  else {
80  cout << "Not enough recHits for Special Segment Algorithm" << endl;
81  return createFakeSeed(isGoodSeed, eSetup);
82  }
83  }
84  } else {
85  if (checkSegment()) {
86  AlgorithmType = 3;
88  } else {
89  AlgorithmType = 2;
91  }
92  }
93  }
94 
95  // Check the pattern
96  if (isPatternChecked == false) {
97  if (AlgorithmType != 3) {
98  checkSimplePattern(eSetup);
99  } else {
101  }
102  }
103 
104  return createSeed(isGoodSeed, eSetup);
105 }
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
void checkSegmentAlgorithmSpecial(const edm::EventSetup &eSetup)
void MiddlePointsAlgorithm()
weightedTrajectorySeed createSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
bool autoAlgorithmChoose
unsigned int nrhit() const
void ThreePointsAlgorithm()
unsigned int AlgorithmType
void checkSimplePattern(const edm::EventSetup &eSetup)
void SegmentAlgorithmSpecial(const edm::EventSetup &eSetup)
bool checkSegment() const
void RPCSeedPattern::SegmentAlgorithm ( )
private

Definition at line 226 of file RPCSeedPattern.cc.

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

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

Definition at line 277 of file RPCSeedPattern.cc.

References volumeBasedMagneticField_1103l_cfi::BValue, gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), SiPixelRawToDigiRegional_cfi::deltaPhi, PbPb_ZMuSkimMuonDPG_cff::deltaR, PVValHelper::dx, PVValHelper::dy, PVValHelper::dz, edm::EventSetup::get(), runTauDisplay::gp, 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().

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

Definition at line 107 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.

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

Friends And Related Function Documentation

friend class RPCSeedFinder
friend

Definition at line 47 of file RPCSeedPattern.h.

Member Data Documentation

unsigned int RPCSeedPattern::AlgorithmType
private

Definition at line 84 of file RPCSeedPattern.h.

bool RPCSeedPattern::autoAlgorithmChoose
private

Definition at line 85 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::Center
private

Definition at line 107 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::Center2
private

Definition at line 98 of file RPCSeedPattern.h.

int RPCSeedPattern::Charge
private

Definition at line 116 of file RPCSeedPattern.h.

double RPCSeedPattern::deltaBz
private

Definition at line 110 of file RPCSeedPattern.h.

double RPCSeedPattern::deltaRThreshold
private

Definition at line 83 of file RPCSeedPattern.h.

GlobalPoint RPCSeedPattern::entryPosition
private

Definition at line 101 of file RPCSeedPattern.h.

int RPCSeedPattern::isClockwise
private

Definition at line 114 of file RPCSeedPattern.h.

bool RPCSeedPattern::isConfigured
private

Definition at line 91 of file RPCSeedPattern.h.

int RPCSeedPattern::isGoodPattern
private

Definition at line 113 of file RPCSeedPattern.h.

int RPCSeedPattern::isParralZ
private

Definition at line 115 of file RPCSeedPattern.h.

bool RPCSeedPattern::isPatternChecked
private

Definition at line 112 of file RPCSeedPattern.h.

bool RPCSeedPattern::isStraight
private

Definition at line 106 of file RPCSeedPattern.h.

bool RPCSeedPattern::isStraight2
private

Definition at line 97 of file RPCSeedPattern.h.

double RPCSeedPattern::lastPhi
private

Definition at line 103 of file RPCSeedPattern.h.

GlobalPoint RPCSeedPattern::leavePosition
private

Definition at line 102 of file RPCSeedPattern.h.

double RPCSeedPattern::MagnecticFieldThreshold
private

Definition at line 95 of file RPCSeedPattern.h.

double RPCSeedPattern::MaxRSD
private

Definition at line 82 of file RPCSeedPattern.h.

double RPCSeedPattern::meanBz
private

Definition at line 109 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::meanMagneticField2
private

Definition at line 96 of file RPCSeedPattern.h.

double RPCSeedPattern::meanPt
private

Definition at line 117 of file RPCSeedPattern.h.

double RPCSeedPattern::meanRadius
private

Definition at line 108 of file RPCSeedPattern.h.

double RPCSeedPattern::meanRadius2
private

Definition at line 99 of file RPCSeedPattern.h.

double RPCSeedPattern::meanSpt
private

Definition at line 118 of file RPCSeedPattern.h.

double RPCSeedPattern::MinDeltaPhi
private

Definition at line 87 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::Momentum
private

Definition at line 119 of file RPCSeedPattern.h.

double RPCSeedPattern::S
private

Definition at line 104 of file RPCSeedPattern.h.

unsigned int RPCSeedPattern::sampleCount
private

Definition at line 89 of file RPCSeedPattern.h.

RPCSegment RPCSeedPattern::SegmentRB[2]
private

Definition at line 100 of file RPCSeedPattern.h.

double RPCSeedPattern::stepLength
private

Definition at line 88 of file RPCSeedPattern.h.

ConstMuonRecHitContainer RPCSeedPattern::theRecHits
private

Definition at line 93 of file RPCSeedPattern.h.

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

double RPCSeedPattern::ZError
private

Definition at line 86 of file RPCSeedPattern.h.