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  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:41
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:65
int station() const
Definition: RPCDetId.h:98
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:66
GlobalPoint entryPosition
T y() const
Definition: PV3DBase.h:57
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:28
T z() const
Definition: PV3DBase.h:58
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:41
GlobalPoint leavePosition
#define upper_limit_pt
T x() const
Definition: PV3DBase.h:56
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:57
#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:28
T z() const
Definition: PV3DBase.h:58
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:41
unsigned int sampleCount
#define upper_limit_pt
T x() const
Definition: PV3DBase.h:56
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:63
tuple cout
Definition: gather_cfg.py:41
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:63
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
tuple cout
Definition: gather_cfg.py:41
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:41
Definition: DDAxes.h:10
void RPCSeedPattern::clear ( void  )
inline
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:57
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:56
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:28
GlobalVector Center
tuple cout
Definition: gather_cfg.py:41
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:28
GlobalVector Center
tuple cout
Definition: gather_cfg.py:41
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;
1091 
1093  eSetup.get<IdealMagneticFieldRecord>().get(Field);
1094 
1095  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &*Field);
1096 
1097  DetId id = best->geographicalId();
1098  TrajectoryStateTransform tsTransform;
1099  PTrajectoryStateOnDet *seedTSOS = tsTransform.persistentState(tsos, id.rawId());
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  delete seedTSOS;
1112  return theweightedSeed;
1113 }
GlobalVector Momentum
void push_back(D *&d)
Definition: OwnVector.h:290
ConstMuonRecHitContainer theRecHits
ConstMuonRecHitPointer BestRefRecHit() const
PTrajectoryStateOnDet * persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid) 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:41
#define upper_limit_pt
Global3DVector GlobalVector
Definition: GlobalVector.h:10
RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createSeed ( int &  isGoodSeed,
const edm::EventSetup eSetup 
)
private

Definition at line 1115 of file RPCSeedPattern.cc.

References alongMomentum, gather_cfg::cout, cross(), debug, Geom::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.

1116 {
1117  if(isPatternChecked == false || isGoodPattern == -1)
1118  {
1119  cout <<"Pattern is not yet checked! Create a fake seed instead!" << endl;
1120  return createFakeSeed(isGoodSeed, eSetup);
1121  }
1122 
1124  eSetup.get<IdealMagneticFieldRecord>().get(Field);
1125 
1127 
1128  //double theMinMomentum = 3.0;
1129  //if(fabs(meanPt) < lower_limit_pt)
1130  //meanPt = lower_limit_pt * meanPt / fabs(meanPt);
1131 
1132  // For pattern we use is Clockwise other than isStraight to estimate charge
1133  if(isClockwise == 0)
1134  Charge = 0;
1135  else
1136  Charge = (int)(meanPt / fabs(meanPt));
1137 
1138  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
1139  const ConstMuonRecHitPointer best = BestRefRecHit();
1141 
1142  if(isClockwise != 0)
1143  {
1144  if(AlgorithmType != 3)
1145  {
1146  // Get the momentum on reference recHit
1147  GlobalVector vecRef1((first->globalPosition().x()-Center.x()), (first->globalPosition().y()-Center.y()), (first->globalPosition().z()-Center.z()));
1148  GlobalVector vecRef2((best->globalPosition().x()-Center.x()), (best->globalPosition().y()-Center.y()), (best->globalPosition().z()-Center.z()));
1149 
1150  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1151  double deltaS = meanRadius * fabs(deltaPhi);
1152  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1153 
1154  GlobalVector vecZ(0, 0, 1);
1155  GlobalVector vecPt = (vecRef2.unit()).cross(vecZ);
1156  if(isClockwise == -1)
1157  vecPt *= -1;
1158  vecPt *= deltaS;
1159  Momentum = GlobalVector(0, 0, deltaZ);
1160  Momentum += vecPt;
1161  Momentum *= fabs(meanPt / deltaS);
1162  }
1163  else
1164  {
1165  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1167  Momentum *= fabs(meanPt / S);
1168  }
1169  }
1170  else
1171  {
1172  Momentum = best->globalPosition() - first->globalPosition();
1173  double deltaS = Momentum.perp();
1174  Momentum *= fabs(meanPt / deltaS);
1175  }
1176  LocalPoint segPos = best->localPosition();
1177  LocalVector segDirFromPos = best->det()->toLocal(Momentum);
1178  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
1179 
1181 
1182  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &*Field);
1183  cout << "Trajectory State on Surface before the extrapolation" << endl;
1184  cout << debug.dumpTSOS(tsos);
1185  DetId id = best->geographicalId();
1186  cout << "The RecSegment relies on: " << endl;
1187  cout << debug.dumpMuonId(id);
1188  cout << debug.dumpTSOS(tsos);
1189  TrajectoryStateTransform tsTransform;
1190  PTrajectoryStateOnDet *seedTSOS = tsTransform.persistentState(tsos, id.rawId());
1191 
1193  for(ConstMuonRecHitContainer::const_iterator iter=theRecHits.begin(); iter!=theRecHits.end(); iter++)
1194  {
1195  // This casting withou clone will cause memory overflow when used in push_back
1196  // Since container's deconstructor functiion free the pointer menber!
1197  //TrackingRecHit* pt = dynamic_cast<TrackingRecHit*>(&*(*iter));
1198  //cout << "Push recHit type " << pt->getType() << endl;
1199  container.push_back((*iter)->hit()->clone());
1200  }
1201 
1202  TrajectorySeed theSeed(*seedTSOS, container, alongMomentum);
1203  weightedTrajectorySeed theweightedSeed;
1204  theweightedSeed.first = theSeed;
1205  theweightedSeed.second = meanSpt;
1206  isGoodSeed = isGoodPattern;
1207 
1208  delete seedTSOS;
1209  return theweightedSeed;
1210 }
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
T perp() const
Definition: PV3DBase.h:66
ConstMuonRecHitPointer FirstRecHit() const
double deltaPhi(float phi1, float phi2)
Definition: VectorUtil.h:30
T y() const
Definition: PV3DBase.h:57
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:58
ConstMuonRecHitContainer theRecHits
ConstMuonRecHitPointer BestRefRecHit() const
PTrajectoryStateOnDet * persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid) const
GlobalVector Center
bool first
Definition: L1TdeRCT.cc:79
Definition: DetId.h:20
RPCSeedPattern::weightedTrajectorySeed weightedTrajectorySeed
const T & get() const
Definition: EventSetup.h:55
unsigned int AlgorithmType
tuple cout
Definition: gather_cfg.py:41
#define debug
Definition: MEtoEDMFormat.h:34
T x() const
Definition: PV3DBase.h:56
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(), Geom::deltaPhi(), edm::EventSetup::get(), Plane::localZ(), 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:66
double deltaPhi(float phi1, float phi2)
Definition: VectorUtil.h:30
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
float localZ(const GlobalPoint &gp) const
Fast access to distance from plane for a point.
Definition: Plane.h:52
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
const RPCChamber * chamber(RPCDetId id) const
Definition: RPCGeometry.cc:66
T mag() const
Definition: PV3DBase.h:61
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:58
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:41
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
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().

Referenced by python.Vispa.Gui.ConnectableWidget.ConnectableWidget::arrangePorts(), python.Vispa.Gui.WidgetContainer.WidgetContainer::autolayoutChildren(), python.Vispa.Views.LineDecayView.LineDecayContainer::autolayoutThreadFinished(), python.Vispa.Gui.ConnectableWidget.ConnectableWidget::centerSinglePortVertically(), python.Vispa.Views.LineDecayView.LineDecayContainer::childrenRect(), python.Vispa.Gui.WidgetContainer.WidgetContainer::contentStartX(), python.Vispa.Gui.WidgetContainer.WidgetContainer::contentStartY(), python.Vispa.Gui.MenuWidget.MenuWidget::drawMenuEntries(), python.Vispa.Gui.ConnectableWidget.ConnectableWidget::drawPortNames(), python.Vispa.Gui.ConnectableWidget.ConnectableWidget::dropArea(), python.Vispa.Gui.WidgetContainer.WidgetContainer::mouseDoubleClickEvent(), python.Vispa.Gui.WidgetContainer.WidgetContainer::mouseMoveEvent(), python.Vispa.Gui.WidgetContainer.WidgetContainer::sizeHint(), and python.Vispa.Gui.ConnectableWidget.ConnectableWidget::sizeHint().

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:57
T sqrt(T t)
Definition: SSEVec.h:28
T x() const
Definition: PV3DBase.h:56
LocalTrajectoryError RPCSeedPattern::getSpecialAlgorithmErrorMatrix ( const ConstMuonRecHitPointer first,
const ConstMuonRecHitPointer best 
)
private

Definition at line 1212 of file RPCSeedPattern.cc.

References beta, gather_cfg::cout, Geom::deltaPhi(), dPhi(), dttmaxenums::L, MultiGaussianStateTransform::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().

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