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 30 of file RPCSeedPattern.h.

Member Typedef Documentation

Definition at line 35 of file RPCSeedPattern.h.

Definition at line 33 of file RPCSeedPattern.h.

Definition at line 34 of file RPCSeedPattern.h.

Definition at line 32 of file RPCSeedPattern.h.

Definition at line 38 of file RPCSeedPattern.h.

Definition at line 39 of file RPCSeedPattern.h.

Constructor & Destructor Documentation

RPCSeedPattern::RPCSeedPattern ( )

Definition at line 33 of file RPCSeedPattern.cc.

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

Definition at line 40 of file RPCSeedPattern.cc.

40  {
41 
42 }

Member Function Documentation

void RPCSeedPattern::add ( const ConstMuonRecHitPointer hit)
inline

Definition at line 46 of file RPCSeedPattern.h.

References theRecHits.

Referenced by counter.Counter::register().

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

Definition at line 524 of file RPCSeedPattern.cc.

525 {
527  int index = 0;
528  // Use the last one for recHit on last layer has minmum delta Z for barrel or delta R for endcap while calculating the momentum
529  // But for Algorithm 3 we use the 4th recHit on the 2nd segment for more accurate
530  for (ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
531  {
532  if(AlgorithmType != 3)
533  best = (*iter);
534  else
535  if(index < 4)
536  best = (*iter);
537  index++;
538  }
539  return best;
540 }
ConstMuonRecHitContainer theRecHits
unsigned int AlgorithmType
MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer
bool RPCSeedPattern::checkSegment ( ) const
private

Definition at line 487 of file RPCSeedPattern.cc.

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

488 {
489  bool isFit = true;
490  unsigned int count = 0;
491  // first 4 recHits should be located in RB1 and RB2
492  for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
493  {
494  count++;
495  const GeomDet* Detector = (*iter)->det();
496  if(dynamic_cast<const RPCChamber*>(Detector) != nullptr)
497  {
498  const RPCChamber* RPCCh = dynamic_cast<const RPCChamber*>(Detector);
499  RPCDetId RPCId = RPCCh->id();
500  int Region = RPCId.region();
501  unsigned int Station = RPCId.station();
502  //int Layer = RPCId.layer();
503  if(count <= 4)
504  {
505  if(Region != 0)
506  isFit = false;
507  if(Station > 2)
508  isFit = false;
509  }
510  }
511  }
512  // more than 4 recHits for pattern building
513  if(count <= 4)
514  isFit = false;
515  cout << "Check for segment fit: " << isFit << endl;
516  return isFit;
517 }
RPCDetId id() const
Return the RPCChamberId of this chamber.
Definition: RPCChamber.cc:33
ConstMuonRecHitContainer theRecHits
Region
Definition: Region.h:7
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:63
int station() const
Definition: RPCDetId.h:96
void RPCSeedPattern::checkSegmentAlgorithmSpecial ( const edm::EventSetup eSetup)
private

Definition at line 823 of file RPCSeedPattern.cc.

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

824 {
825  if(isPatternChecked == true)
826  return;
827 
828  if(!checkSegment())
829  {
830  isPatternChecked = true;
831  isGoodPattern = -1;
832  return;
833  }
834 
835  // Set isGoodPattern to default true and check the failure
836  isGoodPattern = 1;
837 
838  // Print the recHit's position
839  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
840  cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
841 
842  // Check the Z direction
843  if(fabs((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
844  {
845  if(((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
846  isParralZ = 1;
847  else
848  isParralZ = -1;
849  }
850  else
851  isParralZ = 0;
852 
853  cout << " Check isParralZ is :" << isParralZ << endl;
854  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
855  {
856  if(isParralZ == 0)
857  {
858  if(fabs((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z()) > ZError)
859  {
860  cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
861  isGoodPattern = 0;
862  }
863  }
864  else
865  {
866  if((int)(((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z())/ZError)*isParralZ < 0)
867  {
868  cout << "Pattern find error in Z direction: wrong Z direction" << endl;
869  isGoodPattern = 0;
870  }
871  }
872  }
873 
874  // Check the pattern
875  if(isStraight2 == true)
876  {
877  // Set pattern to be straight
878  isClockwise = 0;
880  // Set the straight pattern with lowest priority among good pattern
882 
883  // Extrapolate to other recHits and check deltaR
884  GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
885  GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
886  GlobalVector startMomentum = startSegment*(meanPt/startSegment.perp());
887  unsigned int index = 0;
888  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
889  {
890  if(index < 4)
891  {
892  index++;
893  continue;
894  }
895  double tracklength = 0;
896  cout << "Now checking recHit " << index << endl;
897  double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, eSetup);
898  cout << "Final distance is " << Distance << endl;
899  if(Distance > MaxRSD)
900  {
901  cout << "Pattern find error in distance for other recHits: " << Distance << endl;
902  isGoodPattern = 0;
903  }
904  index++;
905  }
906  }
907  else
908  {
909  // Get clockwise direction
910  GlobalVector vec[2];
913  isClockwise = 0;
914  if((vec[1].phi()-vec[0].phi()).value() > 0)
915  isClockwise = -1;
916  else
917  isClockwise = 1;
918 
919  cout << "Check isClockwise is : " << isClockwise << endl;
920 
921  // Get meanPt
922  meanPt = 0.01 * meanRadius2 * meanMagneticField2.z() * 0.3 * isClockwise;
923  //meanPt = 0.01 * meanRadius2[0] * (-3.8) * 0.3 * isClockwise;
924  cout << " meanRadius is " << meanRadius2 << ", with meanBz " << meanMagneticField2.z() << endl;
925  cout << " meanPt is " << meanPt << endl;
926  if(fabs(meanPt) > upper_limit_pt)
927  meanPt = upper_limit_pt * meanPt / fabs(meanPt);
928 
929  // Check the initial 3 segments
930  cout << "entryPosition: " << entryPosition << endl;
931  cout << "leavePosition: " << leavePosition << endl;
932  cout << "Center2 is : " << Center2 << endl;
933  double R1 = vec[0].perp();
934  double R2 = vec[1].perp();
935  double deltaR = sqrt(((R1-meanRadius2)*(R1-meanRadius2)+(R2-meanRadius2)*(R2-meanRadius2))/2);
936  meanSpt = deltaR;
937  cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
938  cout << "Delta R for the initial 3 segments is " << deltaR << endl;
939  if(deltaR > deltaRThreshold)
940  {
941  cout << "Pattern find error in delta R for the initial 3 segments" << endl;
942  isGoodPattern = 0;
943  }
944 
945  // Extrapolate to other recHits and check deltaR
946  GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
947  GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
948  GlobalVector startMomentum = startSegment*(fabs(meanPt)/startSegment.perp());
949  unsigned int index = 0;
950  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
951  {
952  if(index < 4)
953  {
954  index++;
955  continue;
956  }
957  double tracklength = 0;
958  cout << "Now checking recHit " << index << endl;
959  double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, eSetup);
960  cout << "Final distance is " << Distance << endl;
961  if(Distance > MaxRSD)
962  {
963  cout << "Pattern find error in distance for other recHits: " << Distance << endl;
964  isGoodPattern = 0;
965  }
966  index++;
967  }
968  }
969 
970  cout << "Checking finish, isGoodPattern now is " << isGoodPattern << endl;
971  // Set the pattern estimation signal
972  isPatternChecked = true;
973 }
T perp() const
Definition: PV3DBase.h:72
GlobalPoint entryPosition
T y() const
Definition: PV3DBase.h:63
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:18
T z() const
Definition: PV3DBase.h:64
ConstMuonRecHitContainer theRecHits
GlobalVector Center2
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
GlobalVector meanMagneticField2
bool checkSegment() const
GlobalPoint leavePosition
#define upper_limit_pt
T x() const
Definition: PV3DBase.h:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void RPCSeedPattern::checkSimplePattern ( const edm::EventSetup eSetup)
private

Definition at line 665 of file RPCSeedPattern.cc.

References funct::abs(), gather_cfg::cout, 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().

666 {
667  if(isPatternChecked == true)
668  return;
669 
670  // Get magnetic field
672  eSetup.get<IdealMagneticFieldRecord>().get(Field);
673 
674  unsigned int NumberofHitsinSeed = nrhit();
675 
676  // Print the recHit's position
677  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
678  cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
679 
680  // Get magnetice field information
681  std::vector<double> BzValue;
682  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
683  {
684  GlobalPoint gpFirst = (*iter)->globalPosition();
685  GlobalPoint gpLast = (*(iter+1))->globalPosition();
687  double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
688  double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
689  double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
690  for(unsigned int index = 0; index < sampleCount; index++)
691  {
692  gp[index] = GlobalPoint((gpFirst.x()+dx*(index+1)), (gpFirst.y()+dy*(index+1)), (gpFirst.z()+dz*(index+1)));
693  GlobalVector MagneticVec_temp = Field->inTesla(gp[index]);
694  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
695  BzValue.push_back(MagneticVec_temp.z());
696  }
697  delete [] gp;
698  }
699  meanBz = 0.;
700  for(unsigned int index = 0; index < BzValue.size(); index++)
701  meanBz += BzValue[index];
702  meanBz /= BzValue.size();
703  cout << "Mean Bz is " << meanBz << endl;
704  deltaBz = 0.;
705  for(unsigned int index = 0; index < BzValue.size(); index++)
706  deltaBz += (BzValue[index] - meanBz) * (BzValue[index] - meanBz);
707  deltaBz /= BzValue.size();
708  deltaBz = sqrt(deltaBz);
709  cout<< "delata Bz is " << deltaBz << endl;
710 
711  // Set isGoodPattern to default true and check the failure
712  isGoodPattern = 1;
713 
714  // Check the Z direction
715  if(fabs((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
716  {
717  if(((*(theRecHits.end()-1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
718  isParralZ = 1;
719  else
720  isParralZ = -1;
721  }
722  else
723  isParralZ = 0;
724 
725  cout << " Check isParralZ is :" << isParralZ << endl;
726  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
727  {
728  if(isParralZ == 0)
729  {
730  if(fabs((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z()) > ZError)
731  {
732  cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
733  isGoodPattern = 0;
734  }
735  }
736  else
737  {
738  if((int)(((*(iter+1))->globalPosition().z()-(*iter)->globalPosition().z())/ZError)*isParralZ < 0)
739  {
740  cout << "Pattern find error in Z direction: wrong Z direction" << endl;
741  isGoodPattern = 0;
742  }
743  }
744  }
745 
746  // Check pattern
747  if(isStraight == false)
748  {
749  // Check clockwise direction
750  GlobalVector *vec = new GlobalVector[NumberofHitsinSeed];
751  unsigned int index = 0;
752  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
753  {
754  GlobalVector vec_temp(((*iter)->globalPosition().x()-Center.x()), ((*iter)->globalPosition().y()-Center.y()), ((*iter)->globalPosition().z()-Center.z()));
755  vec[index] = vec_temp;
756  index++;
757  }
758  isClockwise = 0;
759  for(unsigned int index = 0; index < (NumberofHitsinSeed-1); index++)
760  {
761  // Check phi direction, all sub-dphi direction should be the same
762  if((vec[index+1].phi()-vec[index].phi()) > 0)
763  isClockwise--;
764  else
765  isClockwise++;
766  cout << "Current isClockwise is : " << isClockwise << endl;
767  }
768  cout << "Check isClockwise is : " << isClockwise << endl;
769  if((unsigned int)abs(isClockwise) != (NumberofHitsinSeed-1))
770  {
771  cout << "Pattern find error in Phi direction" << endl;
772  isGoodPattern = 0;
773  isClockwise = 0;
774  }
775  else
777  delete [] vec;
778 
779  // Get meanPt and meanSpt
780  double deltaRwithBz = fabs(deltaBz * meanRadius / meanBz);
781  cout << "deltaR with Bz is " << deltaRwithBz << endl;
782 
783  if(isClockwise == 0)
785  else
786  meanPt = 0.01 * meanRadius * meanBz * 0.3 * isClockwise;
787  if(fabs(meanPt) > upper_limit_pt)
788  meanPt = upper_limit_pt * meanPt / fabs(meanPt);
789  cout << " meanRadius is " << meanRadius << endl;
790  cout << " meanPt is " << meanPt << endl;
791 
792  double deltaR = 0.;
793  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
794  {
795  deltaR += (getDistance(*iter, Center) - meanRadius) * (getDistance(*iter, Center) - meanRadius);
796  }
797  deltaR = deltaR / NumberofHitsinSeed;
798  deltaR = sqrt(deltaR);
799  //meanSpt = 0.01 * deltaR * meanBz * 0.3;
800  meanSpt = deltaR;
801  cout << "DeltaR is " << deltaR << endl;
802  if(deltaR > deltaRThreshold)
803  {
804  cout << "Pattern find error: deltaR over threshold" << endl;
805  isGoodPattern = 0;
806  }
807  }
808  else
809  {
810  // Just set pattern to be straight
811  isClockwise =0;
813  // Set the straight pattern with lowest priority among good pattern
815  }
816  cout << "III--> Seed Pt : " << meanPt << endl;
817  cout << "III--> Pattern is: " << isGoodPattern << endl;
818 
819  // Set the pattern estimation signal
820  isPatternChecked = true;
821 }
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
unsigned int nrhit() const
double deltaRThreshold
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
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.
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
T get() const
Definition: EventSetup.h:62
unsigned int sampleCount
#define upper_limit_pt
T x() const
Definition: PV3DBase.h:62
bool RPCSeedPattern::checkStraightwithSegment ( const RPCSegment Segment1,
const RPCSegment Segment2,
double  MinDeltaPhi 
) const
private

Definition at line 587 of file RPCSeedPattern.cc.

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

588 {
589  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
590  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
591  GlobalVector segvec3 = (Segment2.first)->globalPosition() - (Segment1.first)->globalPosition();
592  // compare segvec 1&2 for paralle, 1&3 for straight
593  double dPhi1 = (segvec2.phi() - segvec1.phi()).value();
594  double dPhi2 = (segvec3.phi() - segvec1.phi()).value();
595  cout << "Checking straight with 2 segments. dPhi1: " << dPhi1 << ", dPhi2: " << dPhi2 << endl;
596  cout << "Checking straight with 2 segments. dPhi1 in degree: " << dPhi1*180/3.1415926 << ", dPhi2 in degree: " << dPhi2*180/3.1415926 << endl;
597  if(fabs(dPhi1) > MinDeltaPhi || fabs(dPhi2) > MinDeltaPhi)
598  {
599  cout << "Segment is estimate to be not straight" << endl;
600  return false;
601  }
602  else
603  {
604  cout << "Segment is estimate to be straight" << endl;
605  return true;
606  }
607 }
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
bool RPCSeedPattern::checkStraightwithThreerecHits ( ConstMuonRecHitPointer(&)  precHit[3],
double  MinDeltaPhi 
) const
private

Definition at line 547 of file RPCSeedPattern.cc.

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

548 {
549  GlobalVector segvec1 = precHit[1]->globalPosition() - precHit[0]->globalPosition();
550  GlobalVector segvec2 = precHit[2]->globalPosition() - precHit[1]->globalPosition();
551  double dPhi = (segvec2.phi() - segvec1.phi()).value();
552  if(fabs(dPhi) > MinDeltaPhi)
553  {
554  cout << "Part is estimate to be not straight" << endl;
555  return false;
556  }
557  else
558  {
559  cout << "Part is estimate to be straight" << endl;
560  return true;
561  }
562 }
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
bool RPCSeedPattern::checkStraightwithThreerecHits ( double(&)  x[3],
double(&)  y[3],
double  MinDeltaPhi 
) const
private

Definition at line 632 of file RPCSeedPattern.cc.

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

633 {
634  GlobalVector segvec1((x[1]-x[0]), (y[1]-y[0]), 0);
635  GlobalVector segvec2((x[2]-x[1]), (y[2]-y[1]), 0);
636  double dPhi = (segvec2.phi() - segvec1.phi()).value();
637  if(fabs(dPhi) > MinDeltaPhi)
638  {
639  cout << "Part is estimate to be not straight" << endl;
640  return false;
641  }
642  else
643  {
644  cout << "Part is estimate to be straight" << endl;
645  return true;
646  }
647 }
void RPCSeedPattern::clear ( void  )
inline

Definition at line 45 of file RPCSeedPattern.h.

References theRecHits.

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

Definition at line 609 of file RPCSeedPattern.cc.

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

610 {
611  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
612  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
613  GlobalPoint Point1(((Segment1.second)->globalPosition().x() + (Segment1.first)->globalPosition().x()) / 2, ((Segment1.second)->globalPosition().y() + (Segment1.first)->globalPosition().y()) / 2, ((Segment1.second)->globalPosition().z() + (Segment1.first)->globalPosition().z()) / 2);
614  GlobalPoint Point2(((Segment2.second)->globalPosition().x() + (Segment2.first)->globalPosition().x()) / 2, ((Segment2.second)->globalPosition().y() + (Segment2.first)->globalPosition().y()) / 2, ((Segment2.second)->globalPosition().z() + (Segment2.first)->globalPosition().z()) / 2);
615  GlobalVector vecZ(0, 0, 1);
616  GlobalVector gvec1 = segvec1.cross(vecZ);
617  GlobalVector gvec2 = segvec2.cross(vecZ);
618  double A1 = gvec1.x();
619  double B1 = gvec1.y();
620  double A2 = gvec2.x();
621  double B2 = gvec2.y();
622  double X1 = Point1.x();
623  double Y1 = Point1.y();
624  double X2 = Point2.x();
625  double Y2 = Point2.y();
626  double XO = (A1*A2*(Y2-Y1)+A2*B1*X1-A1*B2*X2)/(A2*B1-A1*B2);
627  double YO = (B1*B2*(X2-X1)+B2*A1*Y1-B1*A2*Y2)/(B2*A1-B1*A2);
628  GlobalVector Center(XO, YO, 0);
629  return Center;
630 }
T y() const
Definition: PV3DBase.h:63
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
GlobalVector Center
T x() const
Definition: PV3DBase.h:62
GlobalVector RPCSeedPattern::computePtwithThreerecHits ( double &  pt,
double &  pt_err,
ConstMuonRecHitPointer(&)  precHit[3] 
) const
private

Definition at line 564 of file RPCSeedPattern.cc.

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

565 {
566  double x[3], y[3];
567  x[0] = precHit[0]->globalPosition().x();
568  y[0] = precHit[0]->globalPosition().y();
569  x[1] = precHit[1]->globalPosition().x();
570  y[1] = precHit[1]->globalPosition().y();
571  x[2] = precHit[2]->globalPosition().x();
572  y[2] = precHit[2]->globalPosition().y();
573  double A = (y[2]-y[1])/(x[2]-x[1]) - (y[1]-y[0])/(x[1]-x[0]);
574  double TYO = (x[2]-x[0])/A + (y[2]*y[2]-y[1]*y[1])/((x[2]-x[1])*A) - (y[1]*y[1]-y[0]*y[0])/((x[1]-x[0])*A);
575  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]);
576  double XO = 0.5 * TXO;
577  double YO = 0.5 * TYO;
578  double R2 = (x[0]-XO)*(x[0]-XO) + (y[0]-YO)*(y[0]-YO);
579  cout << "R2 is " << R2 << endl;
580  // How this algorithm get the pt without magnetic field??
581  pt = 0.01 * sqrt(R2) * 2 * 0.3;
582  cout << "pt is " << pt << endl;
583  GlobalVector Center(XO, YO, 0);
584  return Center;
585 }
T sqrt(T t)
Definition: SSEVec.h:18
GlobalVector Center
GlobalVector RPCSeedPattern::computePtWithThreerecHits ( double &  pt,
double &  pt_err,
double(&)  x[3],
double(&)  y[3] 
) const
private

Definition at line 649 of file RPCSeedPattern.cc.

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

650 {
651  double A = (y[2]-y[1])/(x[2]-x[1]) - (y[1]-y[0])/(x[1]-x[0]);
652  double TYO = (x[2]-x[0])/A + (y[2]*y[2]-y[1]*y[1])/((x[2]-x[1])*A) - (y[1]*y[1]-y[0]*y[0])/((x[1]-x[0])*A);
653  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]);
654  double XO = 0.5 * TXO;
655  double YO = 0.5 * TYO;
656  double R2 = (x[0]-XO)*(x[0]-XO) + (y[0]-YO)*(y[0]-YO);
657  cout << "R2 is " << R2 << endl;
658  // How this algorithm get the pt without magnetic field??
659  pt = 0.01 * sqrt(R2) * 2 * 0.3;
660  cout << "pt is " << pt << endl;
661  GlobalVector Center(XO, YO, 0);
662  return Center;
663 }
T sqrt(T t)
Definition: SSEVec.h:18
GlobalVector Center
void RPCSeedPattern::configure ( const edm::ParameterSet iConfig)

Definition at line 44 of file RPCSeedPattern.cc.

References edm::ParameterSet::getParameter().

44  {
45 
46  MaxRSD = iConfig.getParameter<double>("MaxRSD");
47  deltaRThreshold = iConfig.getParameter<double>("deltaRThreshold");
48  AlgorithmType = iConfig.getParameter<unsigned int>("AlgorithmType");
49  autoAlgorithmChoose = iConfig.getParameter<bool>("autoAlgorithmChoose");
50  ZError = iConfig.getParameter<double>("ZError");
51  MinDeltaPhi = iConfig.getParameter<double>("MinDeltaPhi");
52  MagnecticFieldThreshold = iConfig.getParameter<double>("MagnecticFieldThreshold");
53  stepLength = iConfig.getParameter<double>("stepLength");
54  sampleCount = iConfig.getParameter<unsigned int>("sampleCount");
55  isConfigured = true;
56 }
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 1060 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.

1061 {
1062  // Create a fake seed and return
1063  cout << "Now create a fake seed" << endl;
1064  isPatternChecked = true;
1065  isGoodPattern = -1;
1066  isStraight = true;
1068  meanSpt = 0;
1069  Charge = 0;
1070  isClockwise = 0;
1071  isParralZ = 0;
1072  meanRadius = -1;
1073  //return createSeed(isGoodSeed, eSetup);
1074 
1075  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
1076  const ConstMuonRecHitPointer best = BestRefRecHit();
1077 
1078  Momentum = GlobalVector(0, 0, 0);
1079  LocalPoint segPos=best->localPosition();
1080  LocalVector segDirFromPos=best->det()->toLocal(Momentum);
1081  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
1082 
1083  //AlgebraicVector t(4);
1084  AlgebraicSymMatrix mat(5,0);
1085  mat = best->parametersError().similarityT(best->projectionMatrix());
1086  mat[0][0] = meanSpt;
1087  LocalTrajectoryError error(asSMatrix<5>(mat));
1088 
1090  eSetup.get<IdealMagneticFieldRecord>().get(Field);
1091 
1092  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &*Field);
1093 
1094  DetId id = best->geographicalId();
1095 
1097 
1099  for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
1100  container.push_back((*iter)->hit()->clone());
1101 
1102  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
1103  weightedTrajectorySeed theweightedSeed;
1104  theweightedSeed.first = theSeed;
1105  theweightedSeed.second = meanSpt;
1106  isGoodSeed = isGoodPattern;
1107 
1108  return theweightedSeed;
1109 }
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
GlobalVector Momentum
void push_back(D *&d)
Definition: OwnVector.h:290
ConstMuonRecHitContainer theRecHits
ConstMuonRecHitPointer BestRefRecHit() const
Definition: DetId.h:18
RPCSeedPattern::weightedTrajectorySeed weightedTrajectorySeed
T get() const
Definition: EventSetup.h:62
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 1111 of file RPCSeedPattern.cc.

References alongMomentum, gather_cfg::cout, cross(), debug, hiPixelPairStep_cff::deltaPhi, MuonPatternRecoDumper::dumpMuonId(), MuonPatternRecoDumper::dumpTSOS(), relativeConstraints::error, plotBeamSpotDB::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.

1112 {
1113  if(isPatternChecked == false || isGoodPattern == -1)
1114  {
1115  cout <<"Pattern is not yet checked! Create a fake seed instead!" << endl;
1116  return createFakeSeed(isGoodSeed, eSetup);
1117  }
1118 
1120  eSetup.get<IdealMagneticFieldRecord>().get(Field);
1121 
1123 
1124  //double theMinMomentum = 3.0;
1125  //if(fabs(meanPt) < lower_limit_pt)
1126  //meanPt = lower_limit_pt * meanPt / fabs(meanPt);
1127 
1128  // For pattern we use is Clockwise other than isStraight to estimate charge
1129  if(isClockwise == 0)
1130  Charge = 0;
1131  else
1132  Charge = (int)(meanPt / fabs(meanPt));
1133 
1134  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
1135  const ConstMuonRecHitPointer best = BestRefRecHit();
1137 
1138  if(isClockwise != 0)
1139  {
1140  if(AlgorithmType != 3)
1141  {
1142  // Get the momentum on reference recHit
1143  GlobalVector vecRef1((first->globalPosition().x()-Center.x()), (first->globalPosition().y()-Center.y()), (first->globalPosition().z()-Center.z()));
1144  GlobalVector vecRef2((best->globalPosition().x()-Center.x()), (best->globalPosition().y()-Center.y()), (best->globalPosition().z()-Center.z()));
1145 
1146  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1147  double deltaS = meanRadius * fabs(deltaPhi);
1148  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1149 
1150  GlobalVector vecZ(0, 0, 1);
1151  GlobalVector vecPt = (vecRef2.unit()).cross(vecZ);
1152  if(isClockwise == -1)
1153  vecPt *= -1;
1154  vecPt *= deltaS;
1155  Momentum = GlobalVector(0, 0, deltaZ);
1156  Momentum += vecPt;
1157  Momentum *= fabs(meanPt / deltaS);
1158  }
1159  else
1160  {
1161  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1163  Momentum *= fabs(meanPt / S);
1164  }
1165  }
1166  else
1167  {
1168  Momentum = best->globalPosition() - first->globalPosition();
1169  double deltaS = Momentum.perp();
1170  Momentum *= fabs(meanPt / deltaS);
1171  }
1172  LocalPoint segPos = best->localPosition();
1173  LocalVector segDirFromPos = best->det()->toLocal(Momentum);
1174  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
1175 
1177 
1178  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &*Field);
1179  cout << "Trajectory State on Surface before the extrapolation" << endl;
1180  cout << debug.dumpTSOS(tsos);
1181  DetId id = best->geographicalId();
1182  cout << "The RecSegment relies on: " << endl;
1183  cout << debug.dumpMuonId(id);
1184  cout << debug.dumpTSOS(tsos);
1185 
1186  PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
1187 
1189  for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
1190  {
1191  // This casting withou clone will cause memory overflow when used in push_back
1192  // Since container's deconstructor functiion free the pointer menber!
1193  //TrackingRecHit* pt = dynamic_cast<TrackingRecHit*>(&*(*iter));
1194  //cout << "Push recHit type " << pt->getType() << endl;
1195  container.push_back((*iter)->hit()->clone());
1196  }
1197 
1198  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
1199  weightedTrajectorySeed theweightedSeed;
1200  theweightedSeed.first = theSeed;
1201  theweightedSeed.second = meanSpt;
1202  isGoodSeed = isGoodPattern;
1203 
1204  return theweightedSeed;
1205 }
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
T perp() const
Definition: PV3DBase.h:72
ConstMuonRecHitPointer FirstRecHit() const
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
T y() const
Definition: PV3DBase.h:63
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:290
T z() const
Definition: PV3DBase.h:64
ConstMuonRecHitContainer theRecHits
ConstMuonRecHitPointer BestRefRecHit() const
GlobalVector Center
Definition: DetId.h:18
#define debug
Definition: HDRShower.cc:19
RPCSeedPattern::weightedTrajectorySeed weightedTrajectorySeed
unsigned int AlgorithmType
T get() const
Definition: EventSetup.h:62
T x() const
Definition: PV3DBase.h:62
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 975 of file RPCSeedPattern.cc.

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

976 {
977  // Get magnetic field
979  eSetup.get<IdealMagneticFieldRecord>().get(Field);
980 
981  cout << "Extrapolating the track to check the pattern" << endl;
982  tracklength = 0;
983  // Get the iter recHit's detector geometry
984  DetId hitDet = (*iter)->hit()->geographicalId();
985  RPCDetId RPCId = RPCDetId(hitDet.rawId());
986  //const RPCChamber* hitRPC = dynamic_cast<const RPCChamber*>(hitDet);
988  eSetup.get<MuonGeometryRecord>().get(pRPCGeom);
989  const RPCGeometry* rpcGeometry = (const RPCGeometry*)&*pRPCGeom;
990 
991  const BoundPlane RPCSurface = rpcGeometry->chamber(RPCId)->surface();
992  double startSide = RPCSurface.localZ(startPosition);
993  cout << "Start side : " << startSide;
994 
995  GlobalPoint currentPosition = startPosition;
996  GlobalVector currentMomentum = startMomentum;
997  GlobalVector ZDirection(0, 0, 1);
998 
999  // Use the perp other than mag, since initial segment might have small value while final recHit have large difference value at Z direction
1000  double currentDistance = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
1001  cout << "Start current position is : " << currentPosition << endl;
1002  cout << "Start current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
1003  cout << "Start current distance is " << currentDistance << endl;
1004  cout << "Start current radius is " << currentPosition.perp() << endl;
1005  cout << "Destination radius is " << (*iter)->globalPosition().perp() << endl;
1006 
1007  // Judge roughly if the stepping cross the Det surface of the recHit
1008  //while((currentPosition.perp() < ((*iter)->globalPosition().perp())))
1009  double currentDistance_next = currentDistance;
1010  do
1011  {
1012  currentDistance = currentDistance_next;
1013  if(ClockwiseDirection == 0)
1014  {
1015  currentPosition += currentMomentum.unit() * stepLength;
1016  }
1017  else
1018  {
1019  double Bz = Field->inTesla(currentPosition).z();
1020  double Radius = currentMomentum.perp()/fabs(Bz*0.01*0.3);
1021  double deltaPhi = (stepLength*currentMomentum.perp()/currentMomentum.mag())/Radius;
1022 
1023  // Get the center for current step
1024  GlobalVector currentPositiontoCenter = currentMomentum.unit().cross(ZDirection);
1025  currentPositiontoCenter *= Radius;
1026  // correction of ClockwiseDirection correction
1027  currentPositiontoCenter *= ClockwiseDirection;
1028  // continue to get the center for current step
1029  GlobalPoint currentCenter = currentPosition;
1030  currentCenter += currentPositiontoCenter;
1031 
1032  // Get the next step position
1033  GlobalVector CentertocurrentPosition = (GlobalVector)(currentPosition - currentCenter);
1034  double Phi = CentertocurrentPosition.phi().value();
1035  Phi += deltaPhi * (-1) * ClockwiseDirection;
1036  double deltaZ = stepLength*currentMomentum.z()/currentMomentum.mag();
1037  GlobalVector CentertonewPosition(GlobalVector::Cylindrical(CentertocurrentPosition.perp(), Phi, deltaZ));
1038  double PtPhi = currentMomentum.phi().value();
1039  PtPhi += deltaPhi * (-1) * ClockwiseDirection;
1040  currentMomentum = GlobalVector(GlobalVector::Cylindrical(currentMomentum.perp(), PtPhi, currentMomentum.z()));
1041  currentPosition = currentCenter + CentertonewPosition;
1042  }
1043 
1044  // count the total step length
1045  tracklength += stepLength * currentMomentum.perp() / currentMomentum.mag();
1046 
1047  // Get the next step distance
1048  double currentSide = RPCSurface.localZ(currentPosition);
1049  cout << "Stepping current side : " << currentSide << endl;
1050  cout << "Stepping current position is: " << currentPosition << endl;
1051  cout << "Stepping current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
1052  currentDistance_next = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
1053  cout << "Stepping current distance is " << currentDistance << endl;
1054  cout << "Stepping current radius is " << currentPosition.perp() << endl;
1055  }while(currentDistance_next < currentDistance);
1056 
1057  return currentDistance;
1058 }
T perp() const
Definition: PV3DBase.h:72
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:47
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
const RPCChamber * chamber(RPCDetId id) const
Definition: RPCGeometry.cc:71
T mag() const
Definition: PV3DBase.h:67
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:64
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
T value() const
Explicit access to value in case implicit conversion not OK.
Definition: Phi.h:38
Vector3DBase unit() const
Definition: Vector3DBase.h:57
Definition: DetId.h:18
T get() const
Definition: EventSetup.h:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::FirstRecHit ( ) const
private

Definition at line 519 of file RPCSeedPattern.cc.

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

Definition at line 542 of file RPCSeedPattern.cc.

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

543 {
544  return sqrt((precHit->globalPosition().x()-Center.x())*(precHit->globalPosition().x()-Center.x())+(precHit->globalPosition().y()-Center.y())*(precHit->globalPosition().y()-Center.y()));
545 }
T y() const
Definition: PV3DBase.h:63
T sqrt(T t)
Definition: SSEVec.h:18
T x() const
Definition: PV3DBase.h:62
LocalTrajectoryError RPCSeedPattern::getSpecialAlgorithmErrorMatrix ( const ConstMuonRecHitPointer first,
const ConstMuonRecHitPointer best 
)
private

Definition at line 1207 of file RPCSeedPattern.cc.

References pfBoostedDoubleSVAK8TagInfos_cfi::beta, gather_cfg::cout, hiPixelPairStep_cff::deltaPhi, particleFlow_cfi::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().

1207  {
1208 
1210  double dXdZ = 0;
1211  double dYdZ = 0;
1212  double dP = 0;
1213  AlgebraicSymMatrix mat(5, 0);
1214  mat = best->parametersError().similarityT(best->projectionMatrix());
1215  if(AlgorithmType != 3) {
1216  GlobalVector vecRef1((first->globalPosition().x()-Center.x()), (first->globalPosition().y()-Center.y()), (first->globalPosition().z()-Center.z()));
1217  GlobalVector vecRef2((best->globalPosition().x()-Center.x()), (best->globalPosition().y()-Center.y()), (best->globalPosition().z()-Center.z()));
1218  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1219  double L = meanRadius * fabs(deltaPhi);
1220  double N = nrhit();
1221  double A_N = 180*N*N*N/((N-1)*(N+1)*(N+2)*(N+3));
1222  double sigma_x = sqrt(mat[3][3]);
1223  double betaovergame = Momentum.mag()/0.1066;
1224  double beta = sqrt((betaovergame*betaovergame)/(1+betaovergame*betaovergame));
1225  double dPt = meanPt*(0.0136*sqrt(1/100)*sqrt(4*A_N/N)/(beta*0.3*meanBz*L) + sigma_x*meanPt*sqrt(4*A_N)/(0.3*meanBz*L*L));
1226  double dP = dPt * Momentum.mag() / meanPt;
1227  mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1228  mat[1][1] = dXdZ * dXdZ;
1229  mat[2][2] = dYdZ * dYdZ;
1230  Error = LocalTrajectoryError(asSMatrix<5>(mat));
1231  }
1232  else {
1233  AlgebraicSymMatrix mat0(5, 0);
1234  mat0 = (SegmentRB[0].first)->parametersError().similarityT((SegmentRB[0].first)->projectionMatrix());
1235  double dX0 = sqrt(mat0[3][3]);
1236  double dY0 = sqrt(mat0[4][4]);
1237  AlgebraicSymMatrix mat1(5, 0);
1238  mat1 = (SegmentRB[0].second)->parametersError().similarityT((SegmentRB[0].second)->projectionMatrix());
1239  double dX1 = sqrt(mat1[3][3]);
1240  double dY1 = sqrt(mat1[4][4]);
1241  AlgebraicSymMatrix mat2(5, 0);
1242  mat2 = (SegmentRB[1].first)->parametersError().similarityT((SegmentRB[1].first)->projectionMatrix());
1243  double dX2 = sqrt(mat2[3][3]);
1244  double dY2 = sqrt(mat2[4][4]);
1245  AlgebraicSymMatrix mat3(5, 0);
1246  mat3 = (SegmentRB[1].second)->parametersError().similarityT((SegmentRB[1].second)->projectionMatrix());
1247  double dX3 = sqrt(mat3[3][3]);
1248  double dY3 = sqrt(mat3[4][4]);
1249  cout << "Local error for 4 recHits are: " << dX0 << ", " << dY0 << ", " << dX1 << ", " << dY1 << ", " << dX2 << ", " << dY2 << ", " << dX3 << ", " << dY3 << endl;
1250  const GeomDetUnit* refRPC1 = (SegmentRB[0].second)->detUnit();
1251  LocalPoint recHit0 = refRPC1->toLocal((SegmentRB[0].first)->globalPosition());
1252  LocalPoint recHit1 = refRPC1->toLocal((SegmentRB[0].second)->globalPosition());
1253  LocalVector localSegment00 = (LocalVector)(recHit1 - recHit0);
1254  LocalVector localSegment01 = LocalVector(localSegment00.x()+dX0+dX1, localSegment00.y(), localSegment00.z());
1255  LocalVector localSegment02 = LocalVector(localSegment00.x()-dX0-dX1, localSegment00.y(), localSegment00.z());
1256  GlobalVector globalSegment00 = refRPC1->toGlobal(localSegment00);
1257  GlobalVector globalSegment01 = refRPC1->toGlobal(localSegment01);
1258  GlobalVector globalSegment02 = refRPC1->toGlobal(localSegment02);
1259 
1260  const GeomDetUnit* refRPC2 = (SegmentRB[1].first)->detUnit();
1261  LocalPoint recHit2 = refRPC2->toLocal((SegmentRB[1].first)->globalPosition());
1262  LocalPoint recHit3 = refRPC2->toLocal((SegmentRB[1].second)->globalPosition());
1263  LocalVector localSegment10 = (LocalVector)(recHit3 - recHit2);
1264  LocalVector localSegment11 = LocalVector(localSegment10.x()+dX2+dX3, localSegment10.y(), localSegment10.z());
1265  LocalVector localSegment12 = LocalVector(localSegment10.x()-dX2-dX3, localSegment10.y(), localSegment10.z());
1266  GlobalVector globalSegment10 = refRPC2->toGlobal(localSegment10);
1267  GlobalVector globalSegment11 = refRPC2->toGlobal(localSegment11);
1268  GlobalVector globalSegment12 = refRPC2->toGlobal(localSegment12);
1269 
1270  if(isClockwise != 0) {
1271  GlobalVector vec[2];
1274  double halfPhiCenter = fabs((vec[1].phi() - vec[0].phi()).value()) / 2;
1275  // dPhi0 shoule be the same clockwise direction, while dPhi1 should be opposite clockwise direction, w.r.t to the track clockwise
1276  double dPhi0 = (((globalSegment00.phi() - globalSegment01.phi()).value()*isClockwise) > 0) ? fabs((globalSegment00.phi() - globalSegment01.phi()).value()) : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1277  double dPhi1 = (((globalSegment10.phi() - globalSegment11.phi()).value()*isClockwise) < 0) ? fabs((globalSegment10.phi() - globalSegment11.phi()).value()) : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1278  // For the deltaR should be kept small, we assume the delta Phi0/Phi1 should be in a same limit value
1279  double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1280  cout << "DPhi for new Segment is about " << dPhi << endl;
1281  // Check the variance of halfPhiCenter
1282  double newhalfPhiCenter = ((halfPhiCenter-dPhi) > 0 ? (halfPhiCenter-dPhi) : 0);
1283  if(newhalfPhiCenter != 0) {
1284  double newmeanPt = meanPt * halfPhiCenter / newhalfPhiCenter;
1285  if(fabs(newmeanPt) > upper_limit_pt)
1286  newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1287  cout << "The error is inside range. Max meanPt could be " << newmeanPt << endl;
1288  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1289  }
1290  else {
1291  double newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1292  cout << "The error is outside range. Max meanPt could be " << newmeanPt << endl;
1293  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1294  }
1295  }
1296  else {
1297  double dPhi0 = (fabs((globalSegment00.phi() - globalSegment01.phi()).value()) <= fabs((globalSegment00.phi() - globalSegment02.phi()).value())) ? fabs((globalSegment00.phi() - globalSegment01.phi()).value()) : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1298  double dPhi1 = (fabs((globalSegment10.phi() - globalSegment11.phi()).value()) <= fabs((globalSegment10.phi() - globalSegment12.phi()).value())) ? fabs((globalSegment10.phi() - globalSegment11.phi()).value()) : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1299  double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1300  GlobalVector middleSegment = leavePosition - entryPosition;
1301  double halfDistance = middleSegment.perp() / 2;
1302  double newmeanPt = halfDistance / dPhi;
1303  cout << "The error is for straight. Max meanPt could be " << newmeanPt << endl;
1304  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1305  }
1306 
1307  double dXdZ1 = globalSegment11.x() / globalSegment11.z() - globalSegment10.x() / globalSegment10.z();
1308  double dXdZ2 = globalSegment12.x() / globalSegment12.z() - globalSegment10.x() / globalSegment10.z();
1309  dXdZ = (fabs(dXdZ1) >= fabs(dXdZ2)) ? dXdZ1 : dXdZ2;
1310 
1311  LocalVector localSegment13 = LocalVector(localSegment10.x(), localSegment10.y()+dY2+dY3, localSegment10.z());
1312  LocalVector localSegment14 = LocalVector(localSegment10.x(), localSegment10.y()-dY2-dY3, localSegment10.z());
1313  GlobalVector globalSegment13 = refRPC2->toGlobal(localSegment13);
1314  GlobalVector globalSegment14 = refRPC2->toGlobal(localSegment14);
1315  double dYdZ1 = globalSegment13.y() / globalSegment13.z() - globalSegment10.y() / globalSegment10.z();
1316  double dYdZ2 = globalSegment14.y() / globalSegment14.z() - globalSegment10.y() / globalSegment10.z();
1317  dYdZ = (fabs(dYdZ1) >= fabs(dYdZ2)) ? dYdZ1 : dYdZ2;
1318 
1319  mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1320  mat[1][1] = dXdZ * dXdZ;
1321  mat[2][2] = dYdZ * dYdZ;
1322  Error = LocalTrajectoryError(asSMatrix<5>(mat));
1323  }
1324  return Error;
1325 }
edm::ErrorSummaryEntry Error
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:72
GlobalPoint entryPosition
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
GlobalVector Momentum
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:69
unsigned int nrhit() const
RPCSegment SegmentRB[2]
U second(std::pair< T, U > const &p)
T mag() const
Definition: PV3DBase.h:67
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
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:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void RPCSeedPattern::MiddlePointsAlgorithm ( )
private

Definition at line 193 of file RPCSeedPattern.cc.

References gather_cfg::cout, mps_fire::i, gen::n, EnergyCorrector::pt, X, and DOFs::Y.

194 {
195  cout << "Using middle points algorithm" << endl;
196  unsigned int NumberofHitsinSeed = nrhit();
197  // Check recHit number, if fail we set the pattern to "wrong"
198  if(NumberofHitsinSeed < 4)
199  {
200  isPatternChecked = true;
201  isGoodPattern = -1;
202  return;
203  }
204  double *X = new double[NumberofHitsinSeed];
205  double *Y = new double[NumberofHitsinSeed];
206  unsigned int n = 0;
207  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
208  {
209  X[n] = (*iter)->globalPosition().x();
210  Y[n] = (*iter)->globalPosition().y();
211  cout << "X[" << n <<"] = " << X[n] << ", Y[" << n <<"]= " << Y[n] << endl;
212  n++;
213  }
214  unsigned int NumberofPoints = NumberofHitsinSeed;
215  while(NumberofPoints > 3)
216  {
217  for(unsigned int i = 0; i <= (NumberofPoints - 2); i++)
218  {
219  X[i] = (X[i] + X[i+1]) / 2;
220  Y[i] = (Y[i] + Y[i+1]) / 2;
221  }
222  NumberofPoints--;
223  }
224  double x[3], y[3];
225  for(unsigned int i = 0; i < 3; i++)
226  {
227  x[i] = X[i];
228  y[i] = Y[i];
229  }
230  double pt = 0;
231  double pt_err = 0;
232  bool checkStraight = checkStraightwithThreerecHits(x, y, MinDeltaPhi);
233  if(!checkStraight)
234  {
235 
236  GlobalVector Center_temp = computePtWithThreerecHits(pt, pt_err, x, y);
237  double meanR = 0.;
238  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
239  meanR += getDistance(*iter, Center_temp);
240  meanR /= NumberofHitsinSeed;
241  // For simple pattern
242  isStraight = false;
243  Center = Center_temp;
244  meanRadius = meanR;
245  }
246  else
247  {
248  // For simple pattern
249  isStraight = true;
250  meanRadius = -1;
251  }
252 
253  // Unset the pattern estimation signa
254  isPatternChecked = false;
255 
256  delete [] X;
257  delete [] Y;
258 }
bool checkStraightwithThreerecHits(ConstMuonRecHitPointer(&precHit)[3], double MinDeltaPhi) const
#define X(str)
Definition: MuonsGrabber.cc:48
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 47 of file RPCSeedPattern.h.

References theRecHits.

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

Definition at line 58 of file RPCSeedPattern.cc.

References gather_cfg::cout.

58  {
59 
60  if(isConfigured == false)
61  {
62  cout << "Configuration not set yet" << endl;
63  return createFakeSeed(isGoodSeed, eSetup);
64  }
65 
66  // Check recHit number, if fail we return a fake seed and set pattern to "wrong"
67  unsigned int NumberofHitsinSeed = nrhit();
68  if(NumberofHitsinSeed < 3)
69  return createFakeSeed(isGoodSeed, eSetup);
70  // If only three recHits, we don't have other choice
71  if(NumberofHitsinSeed == 3)
73 
74  if(NumberofHitsinSeed > 3)
75  {
76  if(autoAlgorithmChoose == false)
77  {
78  cout << "computePtWithmorerecHits" << endl;
79  if(AlgorithmType == 0)
81  if(AlgorithmType == 1)
83  if(AlgorithmType == 2)
85  if(AlgorithmType == 3)
86  {
87  if(checkSegment())
89  else
90  {
91  cout << "Not enough recHits for Special Segment Algorithm" << endl;
92  return createFakeSeed(isGoodSeed, eSetup);
93  }
94  }
95  }
96  else
97  {
98  if(checkSegment())
99  {
100  AlgorithmType = 3;
101  SegmentAlgorithmSpecial(eSetup);
102  }
103  else
104  {
105  AlgorithmType = 2;
107  }
108  }
109  }
110 
111  // Check the pattern
112  if(isPatternChecked == false){
113  if(AlgorithmType != 3){
114  checkSimplePattern(eSetup);
115  } else {
117  }
118  }
119 
120  return createSeed(isGoodSeed, eSetup);
121 }
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 260 of file RPCSeedPattern.cc.

References gather_cfg::cout, mps_fire::i, and gen::n.

261 {
262  cout << "Using segments algorithm" << endl;
263  unsigned int NumberofHitsinSeed = nrhit();
264  // Check recHit number, if fail we set the pattern to "wrong"
265  if(NumberofHitsinSeed < 4)
266  {
267  isPatternChecked = true;
268  isGoodPattern = -1;
269  return;
270  }
271 
272  RPCSegment* Segment;
273  unsigned int NumberofSegment = NumberofHitsinSeed - 2;
274  Segment = new RPCSegment[NumberofSegment];
275  unsigned int n = 0;
276  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-2); iter++)
277  {
278  Segment[n].first = (*iter);
279  Segment[n].second = (*(iter + 2));
280  n++;
281  }
282  unsigned int NumberofStraight = 0;
283  for(unsigned int i = 0; i < NumberofSegment - 1; i++)
284  {
285  bool checkStraight = checkStraightwithSegment(Segment[i], Segment[i+1], MinDeltaPhi);
286  if(checkStraight == true)
287  {
288  // For simple patterm
289  NumberofStraight++;
290  }
291  else
292  {
293  GlobalVector Center_temp = computePtwithSegment(Segment[i], Segment[i+1]);
294  // For simple patterm
295  Center += Center_temp;
296  }
297  }
298  // For simple pattern, only one general parameter for pattern
299  if((NumberofSegment-1-NumberofStraight) > 0)
300  {
301  isStraight = false;
302  Center /= (NumberofSegment - 1 - NumberofStraight);
303  double meanR = 0.;
304  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
305  meanR += getDistance(*iter, Center);
306  meanR /= NumberofHitsinSeed;
307  meanRadius = meanR;
308  }
309  else
310  {
311  isStraight = true;
312  meanRadius = -1;
313  }
314 
315  // Unset the pattern estimation signal
316  isPatternChecked = false;
317 
318  delete [] Segment;
319 }
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 321 of file RPCSeedPattern.cc.

References gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), hiPixelPairStep_cff::deltaPhi, deltaR(), PVValHelper::dx, PVValHelper::dy, PVValHelper::dz, edm::EventSetup::get(), runTauDisplay::gp, MagneticField::inTesla(), gen::n, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), mathSSE::sqrt(), Vector3DBase< T, FrameTag >::unit(), Geom::Phi< T >::value(), relativeConstraints::value, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

322 {
323  // Get magnetic field
325  eSetup.get<IdealMagneticFieldRecord>().get(Field);
326 
327  //unsigned int NumberofHitsinSeed = nrhit();
328  if(!checkSegment())
329  {
330  isPatternChecked = true;
331  isGoodPattern = -1;
332  return;
333  }
334 
335  // Get magnetice field sampling information, recHit's position is not the border of Chamber and Iron
336  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end()-1); iter++)
337  {
338  GlobalPoint gpFirst = (*iter)->globalPosition();
339  GlobalPoint gpLast = (*(iter+1))->globalPosition();
341  double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
342  double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
343  double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
344  for(unsigned int index = 0; index < sampleCount; index++)
345  {
346  gp[index] = GlobalPoint((gpFirst.x()+dx*(index+1)), (gpFirst.y()+dy*(index+1)), (gpFirst.z()+dz*(index+1)));
347  GlobalVector MagneticVec_temp = Field->inTesla(gp[index]);
348  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
349  //BValue.push_back(MagneticVec_temp);
350  }
351  delete [] gp;
352  }
353 
354  // form two segments
355  ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin();
356  for(unsigned int n = 0; n <= 1; n++)
357  {
358  SegmentRB[n].first = (*iter);
359  cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
360  iter++;
361  SegmentRB[n].second = (*iter);
362  cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
363  iter++;
364  }
365  GlobalVector segvec1 = (SegmentRB[0].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
366  GlobalVector segvec2 = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
367 
368  // extrapolate the segment to find the Iron border which magnetic field is at large value
369  entryPosition = (SegmentRB[0].second)->globalPosition();
370  leavePosition = (SegmentRB[1].first)->globalPosition();
371  while(fabs(Field->inTesla(entryPosition).z()) < MagnecticFieldThreshold)
372  {
373  cout << "Entry position is : " << entryPosition << ", and stepping into next point" << endl;
374  entryPosition += segvec1.unit() * stepLength;
375  }
376  // Loop back for more accurate by stepLength/10
377  while(fabs(Field->inTesla(entryPosition).z()) >= MagnecticFieldThreshold)
378  {
379  cout << "Entry position is : " << entryPosition << ", and stepping back into next point" << endl;
380  entryPosition -= segvec1.unit() * stepLength / 10;
381  }
382  entryPosition += 0.5 * segvec1.unit() * stepLength / 10;
383  cout << "Final entry position is : " << entryPosition << endl;
384 
385  while(fabs(Field->inTesla(leavePosition).z()) < MagnecticFieldThreshold)
386  {
387  cout << "Leave position is : " << leavePosition << ", and stepping into next point" << endl;
388  leavePosition -= segvec2.unit() * stepLength;
389  }
390  // Loop back for more accurate by stepLength/10
391  while(fabs(Field->inTesla(leavePosition).z()) >= MagnecticFieldThreshold)
392  {
393  cout << "Leave position is : " << leavePosition << ", and stepping back into next point" << endl;
394  leavePosition += segvec2.unit() * stepLength / 10;
395  }
396  leavePosition -= 0.5 * segvec2.unit() * stepLength / 10;
397  cout << "Final leave position is : " << leavePosition << endl;
398 
399  // Sampling magnetic field in Iron region
401  double dx = (leavePosition.x() - entryPosition.x()) / (sampleCount + 1);
402  double dy = (leavePosition.y() - entryPosition.y()) / (sampleCount + 1);
403  double dz = (leavePosition.z() - entryPosition.z()) / (sampleCount + 1);
404  std::vector<GlobalVector> BValue;
405  BValue.clear();
406  for(unsigned int index = 0; index < sampleCount; index++)
407  {
408  gp[index] = GlobalPoint((entryPosition.x()+dx*(index+1)), (entryPosition.y()+dy*(index+1)), (entryPosition.z()+dz*(index+1)));
409  GlobalVector MagneticVec_temp = Field->inTesla(gp[index]);
410  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
411  BValue.push_back(MagneticVec_temp);
412  }
413  delete [] gp;
414  GlobalVector meanB2(0, 0, 0);
415  for(std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
416  meanB2 += (*BIter);
417  meanB2 /= BValue.size();
418  cout << "Mean B field is " << meanB2 << endl;
419  meanMagneticField2 = meanB2;
420 
421  double meanBz2 = meanB2.z();
422  double deltaBz2 = 0.;
423  for(std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
424  deltaBz2 += (BIter->z() - meanBz2) * (BIter->z() - meanBz2);;
425  deltaBz2 /= BValue.size();
426  deltaBz2 = sqrt(deltaBz2);
427  cout<< "delta Bz is " << deltaBz2 << endl;
428 
429  // Distance of the initial 3 segment
430  S = 0;
431  bool checkStraight = checkStraightwithSegment(SegmentRB[0], SegmentRB[1], MinDeltaPhi);
432  if(checkStraight == true)
433  {
434  // Just for complex pattern
435  isStraight2 = checkStraight;
436  Center2 = GlobalVector(0, 0, 0);
437  meanRadius2 = -1;
438  GlobalVector MomentumVec = (SegmentRB[1].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
439  S += MomentumVec.perp();
440  lastPhi = MomentumVec.phi().value();
441  }
442  else
443  {
444  GlobalVector seg1 = entryPosition - (SegmentRB[0].first)->globalPosition();
445  S += seg1.perp();
446  GlobalVector seg2 = (SegmentRB[1].second)->globalPosition() - leavePosition;
447  S += seg2.perp();
448  GlobalVector vecZ(0, 0, 1);
449  GlobalVector gvec1 = seg1.cross(vecZ);
450  GlobalVector gvec2 = seg2.cross(vecZ);
451  double A1 = gvec1.x();
452  double B1 = gvec1.y();
453  double A2 = gvec2.x();
454  double B2 = gvec2.y();
455  double X1 = entryPosition.x();
456  double Y1 = entryPosition.y();
457  double X2 = leavePosition.x();
458  double Y2 = leavePosition.y();
459  double XO = (A1*A2*(Y2-Y1)+A2*B1*X1-A1*B2*X2)/(A2*B1-A1*B2);
460  double YO = (B1*B2*(X2-X1)+B2*A1*Y1-B1*A2*Y2)/(B2*A1-B1*A2);
461  GlobalVector Center_temp(XO, YO, 0);
462  // Just for complex pattern
463  isStraight2 = checkStraight;
464  Center2 = Center_temp;
465 
466  cout << "entryPosition: " << entryPosition << endl;
467  cout << "leavePosition: " << leavePosition << endl;
468  cout << "Center2 is : " << Center_temp << endl;
469 
470  double R1 = GlobalVector((entryPosition.x() - Center_temp.x()), (entryPosition.y() - Center_temp.y()), (entryPosition.z() - Center_temp.z())).perp();
471  double R2 = GlobalVector((leavePosition.x() - Center_temp.x()), (leavePosition.y() - Center_temp.y()), (leavePosition.z() - Center_temp.z())).perp();
472  double meanR = (R1 + R2) / 2;
473  double deltaR = sqrt(((R1-meanR)*(R1-meanR)+(R2-meanR)*(R2-meanR))/2);
474  meanRadius2 = meanR;
475  cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
476  cout << "Mean radius is " << meanR << endl;
477  cout << "Delta R is " << deltaR << endl;
478  double deltaPhi = fabs(((leavePosition-GlobalPoint(XO, YO, 0)).phi()-(entryPosition-GlobalPoint(XO, YO, 0)).phi()).value());
479  S += meanR * deltaPhi;
480  lastPhi = seg2.phi().value();
481  }
482 
483  // Unset the pattern estimation signa
484  isPatternChecked = false;
485 }
T perp() const
Definition: PV3DBase.h:72
double MagnecticFieldThreshold
GlobalPoint entryPosition
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
RPCSegment SegmentRB[2]
T sqrt(T t)
Definition: SSEVec.h:18
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:64
ConstMuonRecHitContainer theRecHits
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
GlobalVector Center2
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
T value() const
Explicit access to value in case implicit conversion not OK.
Definition: Phi.h:38
Vector3DBase unit() const
Definition: Vector3DBase.h:57
GlobalVector meanMagneticField2
T get() const
Definition: EventSetup.h:62
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:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void RPCSeedPattern::ThreePointsAlgorithm ( )
private

Definition at line 123 of file RPCSeedPattern.cc.

References gather_cfg::cout, mps_fire::i, gen::k, gen::n, EnergyCorrector::pt, and upper_limit_pt.

124 {
125  cout << "computePtWith3recHits" << endl;
126  unsigned int NumberofHitsinSeed = nrhit();
127  // Check recHit number, if fail we set the pattern to "wrong"
128  if(NumberofHitsinSeed < 3)
129  {
130  isPatternChecked = true;
131  isGoodPattern = -1;
132  return;
133  }
134  // Choose every 3 recHits to form a part
135  unsigned int NumberofPart = NumberofHitsinSeed * (NumberofHitsinSeed - 1) * (NumberofHitsinSeed - 2) / (3 * 2);;
136  double *pt = new double[NumberofPart];
137  double *pt_err = new double[NumberofPart];
138  // Loop for each three-recHits part
139  ConstMuonRecHitPointer precHit[3];
140  unsigned int n = 0;
141  unsigned int NumberofStraight = 0;
142  for(unsigned int i = 0; i < (NumberofHitsinSeed - 2); i++)
143  for(unsigned int j = (i + 1); j < (NumberofHitsinSeed - 1); j++)
144  for(unsigned int k = (j + 1); k < NumberofHitsinSeed; k++)
145  {
146  precHit[0] = theRecHits[i];
147  precHit[1] = theRecHits[j];
148  precHit[2] = theRecHits[k];
149  bool checkStraight = checkStraightwithThreerecHits(precHit, MinDeltaPhi);
150  if(!checkStraight)
151  {
152  GlobalVector Center_temp = computePtwithThreerecHits(pt[n], pt_err[n], precHit);
153  // For simple pattern
154  Center += Center_temp;
155  }
156  else
157  {
158  // For simple pattern
159  NumberofStraight++;
160  pt[n] = upper_limit_pt;
161  pt_err[n] = 0;
162  }
163  n++;
164  }
165  // For simple pattern, only one general parameter for pattern
166  if(NumberofStraight == NumberofPart)
167  {
168  isStraight = true;
169  meanRadius = -1;
170  }
171  else
172  {
173  isStraight = false;
174  Center /= (NumberofPart - NumberofStraight);
175  double meanR = 0.;
176  for(ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
177  meanR += getDistance(*iter, Center);
178  meanR /= NumberofHitsinSeed;
179  meanRadius = meanR;
180  }
181 
182  // Unset the pattern estimation signa
183  isPatternChecked = false;
184 
185  //double ptmean0 = 0;
186  //double sptmean0 = 0;
187  //computeBestPt(pt, pt_err, ptmean0, sptmean0, (NumberofPart - NumberofStraight));
188 
189  delete [] pt;
190  delete [] pt_err;
191 }
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
int k[5][pyjets_maxn]
#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 50 of file RPCSeedPattern.h.

Member Data Documentation

unsigned int RPCSeedPattern::AlgorithmType
private

Definition at line 81 of file RPCSeedPattern.h.

bool RPCSeedPattern::autoAlgorithmChoose
private

Definition at line 82 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::Center
private

Definition at line 104 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::Center2
private

Definition at line 95 of file RPCSeedPattern.h.

int RPCSeedPattern::Charge
private

Definition at line 113 of file RPCSeedPattern.h.

double RPCSeedPattern::deltaBz
private

Definition at line 107 of file RPCSeedPattern.h.

double RPCSeedPattern::deltaRThreshold
private

Definition at line 80 of file RPCSeedPattern.h.

GlobalPoint RPCSeedPattern::entryPosition
private

Definition at line 98 of file RPCSeedPattern.h.

int RPCSeedPattern::isClockwise
private

Definition at line 111 of file RPCSeedPattern.h.

bool RPCSeedPattern::isConfigured
private

Definition at line 88 of file RPCSeedPattern.h.

int RPCSeedPattern::isGoodPattern
private

Definition at line 110 of file RPCSeedPattern.h.

int RPCSeedPattern::isParralZ
private

Definition at line 112 of file RPCSeedPattern.h.

bool RPCSeedPattern::isPatternChecked
private

Definition at line 109 of file RPCSeedPattern.h.

bool RPCSeedPattern::isStraight
private

Definition at line 103 of file RPCSeedPattern.h.

bool RPCSeedPattern::isStraight2
private

Definition at line 94 of file RPCSeedPattern.h.

double RPCSeedPattern::lastPhi
private

Definition at line 100 of file RPCSeedPattern.h.

GlobalPoint RPCSeedPattern::leavePosition
private

Definition at line 99 of file RPCSeedPattern.h.

double RPCSeedPattern::MagnecticFieldThreshold
private

Definition at line 92 of file RPCSeedPattern.h.

double RPCSeedPattern::MaxRSD
private

Definition at line 79 of file RPCSeedPattern.h.

double RPCSeedPattern::meanBz
private

Definition at line 106 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::meanMagneticField2
private

Definition at line 93 of file RPCSeedPattern.h.

double RPCSeedPattern::meanPt
private

Definition at line 114 of file RPCSeedPattern.h.

double RPCSeedPattern::meanRadius
private

Definition at line 105 of file RPCSeedPattern.h.

double RPCSeedPattern::meanRadius2
private

Definition at line 96 of file RPCSeedPattern.h.

double RPCSeedPattern::meanSpt
private

Definition at line 115 of file RPCSeedPattern.h.

double RPCSeedPattern::MinDeltaPhi
private

Definition at line 84 of file RPCSeedPattern.h.

GlobalVector RPCSeedPattern::Momentum
private

Definition at line 116 of file RPCSeedPattern.h.

double RPCSeedPattern::S
private

Definition at line 101 of file RPCSeedPattern.h.

unsigned int RPCSeedPattern::sampleCount
private

Definition at line 86 of file RPCSeedPattern.h.

RPCSegment RPCSeedPattern::SegmentRB[2]
private

Definition at line 97 of file RPCSeedPattern.h.

double RPCSeedPattern::stepLength
private

Definition at line 85 of file RPCSeedPattern.h.

ConstMuonRecHitContainer RPCSeedPattern::theRecHits
private

Definition at line 90 of file RPCSeedPattern.h.

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

double RPCSeedPattern::ZError
private

Definition at line 83 of file RPCSeedPattern.h.