CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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,
ConstMuonRecHitPointer
RPCSegment
 
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 35 of file RPCSeedPattern.cc.

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

Definition at line 42 of file RPCSeedPattern.cc.

42  {
43 
44 }

Member Function Documentation

void RPCSeedPattern::add ( const ConstMuonRecHitPointer hit)
inline

Definition at line 46 of file RPCSeedPattern.h.

References theRecHits.

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

Definition at line 526 of file RPCSeedPattern.cc.

References getHLTprescales::index.

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

Definition at line 489 of file RPCSeedPattern.cc.

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

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

Definition at line 825 of file RPCSeedPattern.cc.

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

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

Definition at line 667 of file RPCSeedPattern.cc.

References abs, gather_cfg::cout, deltaR(), edm::EventSetup::get(), getHLTprescales::index, phi, mathSSE::sqrt(), upper_limit_pt, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

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

Definition at line 589 of file RPCSeedPattern.cc.

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

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

Definition at line 549 of file RPCSeedPattern.cc.

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

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

Definition at line 634 of file RPCSeedPattern.cc.

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

635 {
636  GlobalVector segvec1((x[1]-x[0]), (y[1]-y[0]), 0);
637  GlobalVector segvec2((x[2]-x[1]), (y[2]-y[1]), 0);
638  double dPhi = (segvec2.phi() - segvec1.phi()).value();
639  if(fabs(dPhi) > MinDeltaPhi)
640  {
641  cout << "Part is estimate to be not straight" << endl;
642  return false;
643  }
644  else
645  {
646  cout << "Part is estimate to be straight" << endl;
647  return true;
648  }
649 }
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
tuple cout
Definition: gather_cfg.py:121
Definition: DDAxes.h:10
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 611 of file RPCSeedPattern.cc.

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

612 {
613  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
614  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
615  GlobalPoint Point1(((Segment1.second)->globalPosition().x() + (Segment1.first)->globalPosition().x()) / 2, ((Segment1.second)->globalPosition().y() + (Segment1.first)->globalPosition().y()) / 2, ((Segment1.second)->globalPosition().z() + (Segment1.first)->globalPosition().z()) / 2);
616  GlobalPoint Point2(((Segment2.second)->globalPosition().x() + (Segment2.first)->globalPosition().x()) / 2, ((Segment2.second)->globalPosition().y() + (Segment2.first)->globalPosition().y()) / 2, ((Segment2.second)->globalPosition().z() + (Segment2.first)->globalPosition().z()) / 2);
617  GlobalVector vecZ(0, 0, 1);
618  GlobalVector gvec1 = segvec1.cross(vecZ);
619  GlobalVector gvec2 = segvec2.cross(vecZ);
620  double A1 = gvec1.x();
621  double B1 = gvec1.y();
622  double A2 = gvec2.x();
623  double B2 = gvec2.y();
624  double X1 = Point1.x();
625  double Y1 = Point1.y();
626  double X2 = Point2.x();
627  double Y2 = Point2.y();
628  double XO = (A1*A2*(Y2-Y1)+A2*B1*X1-A1*B2*X2)/(A2*B1-A1*B2);
629  double YO = (B1*B2*(X2-X1)+B2*A1*Y1-B1*A2*Y2)/(B2*A1-B1*A2);
630  GlobalVector Center(XO, YO, 0);
631  return Center;
632 }
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 566 of file RPCSeedPattern.cc.

References funct::A, gather_cfg::cout, mathSSE::sqrt(), x, and detailsBasic3DVector::y.

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

Definition at line 651 of file RPCSeedPattern.cc.

References funct::A, gather_cfg::cout, mathSSE::sqrt(), x, and detailsBasic3DVector::y.

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

Definition at line 46 of file RPCSeedPattern.cc.

References edm::ParameterSet::getParameter().

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

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

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

Definition at line 1114 of file RPCSeedPattern.cc.

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

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

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

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

Definition at line 521 of file RPCSeedPattern.cc.

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

Definition at line 544 of file RPCSeedPattern.cc.

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

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

Definition at line 1210 of file RPCSeedPattern.cc.

References beta, gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, dPhi(), dttmaxenums::L, N, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), edm::second(), 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().

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

Definition at line 195 of file RPCSeedPattern.cc.

References gather_cfg::cout, i, n, x, X, and detailsBasic3DVector::y.

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

References gather_cfg::cout.

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

Definition at line 262 of file RPCSeedPattern.cc.

References gather_cfg::cout, i, and n.

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

Definition at line 323 of file RPCSeedPattern.cc.

References MagneticField_38T_polyFit2D_cff::BValue, gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), SiPixelRawToDigiRegional_cfi::deltaPhi, deltaR(), edm::EventSetup::get(), getHLTprescales::index, n, PV3DBase< T, PVType, FrameType >::perp(), phi, 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().

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

Definition at line 125 of file RPCSeedPattern.cc.

References gather_cfg::cout, i, j, gen::k, n, and upper_limit_pt.

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

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.