CMS 3D CMS Logo
List of all members | Public Types | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes
ReferenceTrajectory Class Reference

#include <ReferenceTrajectory.h>

Inheritance diagram for ReferenceTrajectory:
ReferenceTrajectoryBase ReferenceCounted BzeroReferenceTrajectory

Public Types

typedef SurfaceSideDefinition::SurfaceSide SurfaceSide
 
- Public Types inherited from ReferenceTrajectoryBase
enum  MaterialEffects {
  none, multipleScattering, energyLoss, combined,
  breakPoints, brokenLinesCoarse, brokenLinesFine, localGBL,
  curvlinGBL
}
 
typedef ReferenceCountingPointer< ReferenceTrajectoryBase > ReferenceTrajectoryPtr
 

Public Member Functions

virtual ReferenceTrajectory * clone () const
 
 ReferenceTrajectory (const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, const MagneticField *magField, const reco::BeamSpot &beamSpot, const ReferenceTrajectoryBase::Config &config)
 
virtual ~ReferenceTrajectory ()
 
- Public Member Functions inherited from ReferenceTrajectoryBase
const AlgebraicMatrix & derivatives () const
 
const Eigen::MatrixXd & gblExtDerivatives () const
 
const Eigen::VectorXd & gblExtMeasurements () const
 
const Eigen::VectorXd & gblExtPrecisions () const
 
std::vector< std::pair< std::vector< gbl::GblPoint >, Eigen::MatrixXd > > & gblInput ()
 
bool isValid ()
 
const AlgebraicMatrix & localToTrajectory () const
 
const AlgebraicSymMatrix & measurementErrors () const
 
const AlgebraicVector & measurements () const
 
int nominalField () const
 
unsigned int numberOfHitMeas () const
 
unsigned int numberOfHits () const
 
unsigned int numberOfPar () const
 
unsigned int numberOfVirtualMeas () const
 
unsigned int numberOfVirtualPar () const
 
const AlgebraicSymMatrix & parameterErrors () const
 
bool parameterErrorsAvailable () const
 
const AlgebraicVector & parameters () const
 
const TransientTrackingRecHit::ConstRecHitContainer & recHits () const
 
void setParameterErrors (const AlgebraicSymMatrix &error)
 
const AlgebraicSymMatrix & trajectoryPositionErrors () const
 
const AlgebraicVector & trajectoryPositions () const
 
const std::vector< TrajectoryStateOnSurface > & trajectoryStates () const
 
const AlgebraicMatrix & trajectoryToCurv () const
 
virtual ~ReferenceTrajectoryBase ()
 

Protected Member Functions

virtual bool addMaterialEffectsBp (const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv)
 
virtual bool addMaterialEffectsBrl (const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv, const GlobalTrajectoryParameters &gtp)
 
virtual bool addMaterialEffectsBrl (const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv, const std::vector< double > &allSteps, const GlobalTrajectoryParameters &gtp, const double minStep=1.0)
 
virtual bool addMaterialEffectsCov (const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs)
 
virtual bool addMaterialEffectsCurvlinGbl (const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv)
 
virtual bool addMaterialEffectsLocalGbl (const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvatureChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParameterCovs)
 
virtual bool construct (const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, const MagneticField *magField, const reco::BeamSpot &beamSpot)
 
MaterialEffectsUpdator * createUpdator (MaterialEffects materialEffects, double mass) const
 
virtual void fillDerivatives (const AlgebraicMatrix &projection, const AlgebraicMatrix &fullJacobian, unsigned int iRow)
 
virtual void fillMeasurementAndError (const TransientTrackingRecHit::ConstRecHitPointer &hitPtr, unsigned int iRow, const TrajectoryStateOnSurface &updatedTsos)
 
virtual void fillTrajectoryPositions (const AlgebraicMatrix &projection, const AlgebraicVector &mixedLocalParams, unsigned int iRow)
 
AlgebraicMatrix getHitProjectionMatrix (const TransientTrackingRecHit::ConstRecHitPointer &recHit) const
 
template<unsigned int N>
AlgebraicMatrix getHitProjectionMatrixT (const TransientTrackingRecHit::ConstRecHitPointer &recHit) const
 
virtual bool propagate (const Plane &previousSurface, const TrajectoryStateOnSurface &previousTsos, const Plane &newSurface, TrajectoryStateOnSurface &newTsos, AlgebraicMatrix &newJacobian, AlgebraicMatrix &newCurvlinJacobian, double &nextStep, const MagneticField *magField) const
 
 ReferenceTrajectory (unsigned int nPar, unsigned int nHits, const ReferenceTrajectoryBase::Config &config)
 
SurfaceSide surfaceSide (const PropagationDirection dir) const
 
- Protected Member Functions inherited from ReferenceTrajectoryBase
unsigned int numberOfUsedRecHits (const TransientTrackingRecHit::ConstRecHitContainer &recHits) const
 
 ReferenceTrajectoryBase (unsigned int nPar, unsigned int nHits, unsigned int nVirtualPar, unsigned int nVirtualMeas)
 
bool useRecHit (const TransientTrackingRecHit::ConstRecHitPointer &hitPtr) const
 

Private Member Functions

template<typename Derived >
void clhep2eigen (const AlgebraicVector &in, Eigen::MatrixBase< Derived > &out)
 
template<typename Derived >
void clhep2eigen (const AlgebraicMatrix &in, Eigen::MatrixBase< Derived > &out)
 
template<typename Derived >
void clhep2eigen (const AlgebraicSymMatrix &in, Eigen::MatrixBase< Derived > &out)
 

Private Attributes

const bool allowZeroMaterial_
 
const bool includeAPEs_
 
const double mass_
 
const MaterialEffects materialEffects_
 
const PropagationDirection propDir_
 
const bool useBeamSpot_
 

Additional Inherited Members

- Protected Attributes inherited from ReferenceTrajectoryBase
AlgebraicMatrix theDerivatives
 
Eigen::MatrixXd theGblExtDerivatives
 
Eigen::VectorXd theGblExtMeasurements
 
Eigen::VectorXd theGblExtPrecisions
 
std::vector< std::pair< std::vector< gbl::GblPoint >, Eigen::MatrixXd > > theGblInput
 
AlgebraicMatrix theInnerLocalToTrajectory
 
AlgebraicMatrix theInnerTrajectoryToCurvilinear
 
AlgebraicVector theMeasurements
 
AlgebraicSymMatrix theMeasurementsCov
 
int theNomField
 
unsigned int theNumberOfHits
 
unsigned int theNumberOfPars
 
unsigned int theNumberOfVirtualMeas
 
unsigned int theNumberOfVirtualPars
 
bool theParamCovFlag
 
AlgebraicSymMatrix theParameterCov
 
AlgebraicVector theParameters
 
TransientTrackingRecHit::ConstRecHitContainer theRecHits
 
AlgebraicSymMatrix theTrajectoryPositionCov
 
AlgebraicVector theTrajectoryPositions
 
std::vector< TrajectoryStateOnSurface > theTsosVec
 
bool theValidityFlag
 
- Static Protected Attributes inherited from ReferenceTrajectoryBase
static unsigned int nMeasPerHit {2}
 

Detailed Description

Definition at line 53 of file ReferenceTrajectory.h.

Member Typedef Documentation

typedef SurfaceSideDefinition::SurfaceSide ReferenceTrajectory::SurfaceSide

Definition at line 58 of file ReferenceTrajectory.h.

Constructor & Destructor Documentation

ReferenceTrajectory::ReferenceTrajectory ( const TrajectoryStateOnSurface &  referenceTsos,
const TransientTrackingRecHit::ConstRecHitContainer &  recHits,
const MagneticField *  magField,
const reco::BeamSpot &  beamSpot,
const ReferenceTrajectoryBase::Config &  config 
)

Constructor with Tsos at first hit (in physical order) and list of hits [if (hitsAreReverse) ==> order of hits is in opposite direction compared to the flight of particle, but note that ReferenceTrajectory::recHits() returns the hits always in order of flight], the material effects to be considered and a particle mass, the magnetic field and beamSpot of the event are needed for propagations etc.

Definition at line 53 of file ReferenceTrajectory.cc.

References construct(), ReferenceTrajectoryBase::Config::hitsAreReverse, TrajectoryStateOnSurface::localParameters(), LocalTrajectoryParameters::mixedFormatVector(), ReferenceTrajectoryBase::theParameters, and ReferenceTrajectoryBase::theValidityFlag.

57  :
58  ReferenceTrajectoryBase(
59  (config.materialEffects >= brokenLinesCoarse) ? 1 : refTsos.localParameters().mixedFormatVector().kSize,
60  (config.useBeamSpot) ? recHits.size()+1 : recHits.size(),
61  (config.materialEffects >= brokenLinesCoarse) ?
62  2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size()) :
63  ( (config.materialEffects == breakPoints) ? 2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size())-2 : 0) ,
64  (config.materialEffects >= brokenLinesCoarse) ?
65  2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size())-4 :
66  ( (config.materialEffects == breakPoints) ? 2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size())-2 : 0) ),
67  mass_(config.mass),
68  materialEffects_(config.materialEffects),
69  propDir_(config.propDir),
70  useBeamSpot_(config.useBeamSpot),
71  includeAPEs_(config.includeAPEs),
72  allowZeroMaterial_(config.allowZeroMaterial)
73 {
74  // no check against magField == 0
75  theParameters = asHepVector<5>( refTsos.localParameters().mixedFormatVector() );
76 
77  if (config.hitsAreReverse) {
78  TransientTrackingRecHit::ConstRecHitContainer fwdRecHits;
79  fwdRecHits.reserve(recHits.size());
80  for (TransientTrackingRecHit::ConstRecHitContainer::const_reverse_iterator it=recHits.rbegin();
81  it != recHits.rend(); ++it) {
82  fwdRecHits.push_back(*it);
83  }
84  theValidityFlag = this->construct(refTsos, fwdRecHits, magField, beamSpot);
85  } else {
86  theValidityFlag = this->construct(refTsos, recHits, magField, beamSpot);
87  }
88 }
ReferenceTrajectoryBase::ReferenceTrajectoryBase
ReferenceTrajectoryBase(unsigned int nPar, unsigned int nHits, unsigned int nVirtualPar, unsigned int nVirtualMeas)
Definition: ReferenceTrajectoryBase.cc:3
ReferenceTrajectoryBase::recHits
const TransientTrackingRecHit::ConstRecHitContainer & recHits() const
Definition: ReferenceTrajectoryBase.h:209
ReferenceTrajectory::construct
virtual bool construct(const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, const MagneticField *magField, const reco::BeamSpot &beamSpot)
Definition: ReferenceTrajectory.cc:111
ReferenceTrajectory::propDir_
const PropagationDirection propDir_
Definition: ReferenceTrajectory.h:194
ReferenceTrajectory::materialEffects_
const MaterialEffects materialEffects_
Definition: ReferenceTrajectory.h:193
TrackingRecHit::ConstRecHitContainer
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: TrackingRecHit.h:35
virtual ReferenceTrajectory::~ReferenceTrajectory ( )
inlinevirtual

Definition at line 73 of file ReferenceTrajectory.h.

73 {}
ReferenceTrajectory::ReferenceTrajectory ( unsigned int  nPar,
unsigned int  nHits,
const ReferenceTrajectoryBase::Config &  config 
)
protected

Definition at line 93 of file ReferenceTrajectory.cc.

95  : ReferenceTrajectoryBase(
96  (config.materialEffects >= brokenLinesCoarse) ? 1 : nPar,
97  nHits,
98  (config.materialEffects >= brokenLinesCoarse) ? 2*nHits : ( (config.materialEffects == breakPoints) ? 2*nHits-2 : 0 ),
99  (config.materialEffects >= brokenLinesCoarse) ? 2*nHits-4 : ( (config.materialEffects == breakPoints) ? 2*nHits-2 : 0 ) ),
100  mass_(config.mass),
101  materialEffects_(config.materialEffects),
102  propDir_(config.propDir),
103  useBeamSpot_(config.useBeamSpot),
104  includeAPEs_(config.includeAPEs),
105  allowZeroMaterial_(config.allowZeroMaterial)
106 {}
ReferenceTrajectoryBase::ReferenceTrajectoryBase
ReferenceTrajectoryBase(unsigned int nPar, unsigned int nHits, unsigned int nVirtualPar, unsigned int nVirtualMeas)
Definition: ReferenceTrajectoryBase.cc:3
ReferenceTrajectory::propDir_
const PropagationDirection propDir_
Definition: ReferenceTrajectory.h:194
ReferenceTrajectory::materialEffects_
const MaterialEffects materialEffects_
Definition: ReferenceTrajectory.h:193

Member Function Documentation

bool ReferenceTrajectory::addMaterialEffectsBp ( const std::vector< AlgebraicMatrix > &  allJacobians,
const std::vector< AlgebraicMatrix > &  allProjections,
const std::vector< AlgebraicSymMatrix > &  allCurvChanges,
const std::vector< AlgebraicSymMatrix > &  allDeltaParaCovs,
const std::vector< AlgebraicMatrix > &  allLocalToCurv 
)
protectedvirtual

internal method to add material effects using break points

Definition at line 556 of file ReferenceTrajectory.cc.

References gen::k, checklumidiff::l, ReferenceTrajectoryBase::nMeasPerHit, ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theNumberOfPars.

Referenced by construct().

561 {
562 //CHK: add material effects using break points
563  int offsetPar = theNumberOfPars;
564  int offsetMeas = nMeasPerHit * allJacobians.size();
565  int ierr = 0;
566 
567  AlgebraicMatrix tempJacobian;
568  AlgebraicMatrix MSprojection(2,5);
569  MSprojection[0][1] = 1;
570  MSprojection[1][2] = 1;
571  AlgebraicSymMatrix tempMSCov;
572  AlgebraicSymMatrix tempMSCovProj;
573  AlgebraicMatrix tempMSJacProj;
574 
575  for (unsigned int k = 1; k < allJacobians.size(); ++k) {
576 // CHK
577  int kbp = k-1;
578  tempJacobian = allJacobians[k] * allCurvatureChanges[k];
579  tempMSCov = allDeltaParameterCovs[k-1].similarity(allLocalToCurv[k-1]);
580  tempMSCovProj = tempMSCov.similarity(MSprojection);
581  theMeasurementsCov[offsetMeas+nMeasPerHit*kbp ][offsetMeas+nMeasPerHit*kbp ] = tempMSCovProj[0][0];
582  theMeasurementsCov[offsetMeas+nMeasPerHit*kbp+1][offsetMeas+nMeasPerHit*kbp+1]= tempMSCovProj[1][1];
583  theDerivatives[offsetMeas+nMeasPerHit*kbp ][offsetPar+2*kbp ] = 1.0;
584  theDerivatives[offsetMeas+nMeasPerHit*kbp+1][offsetPar+2*kbp+1] = 1.0 ;
585 
586  tempMSJacProj = (allProjections[k] * ( tempJacobian * allLocalToCurv[k-1].inverse(ierr) )) * MSprojection.T();
587  if (ierr) {
588  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBp"
589  << "Inversion 1 for break points failed: " << ierr;
590  return false;
591  }
592  theDerivatives[nMeasPerHit*k ][offsetPar+2*kbp ] = tempMSJacProj[0][0];
593  theDerivatives[nMeasPerHit*k ][offsetPar+2*kbp+1] = tempMSJacProj[0][1];
594  theDerivatives[nMeasPerHit*k+1][offsetPar+2*kbp ] = tempMSJacProj[1][0];
595  theDerivatives[nMeasPerHit*k+1][offsetPar+2*kbp+1] = tempMSJacProj[1][1];
596 
597  for (unsigned int l = k+1; l < allJacobians.size(); ++l) {
598 // CHK
599  int kbp = k-1;
600  tempJacobian = allJacobians[l] * allCurvatureChanges[l] * tempJacobian;
601  tempMSJacProj = (allProjections[l] * ( tempJacobian * allLocalToCurv[k-1].inverse(ierr) )) * MSprojection.T();
602  if (ierr) {
603  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBp"
604  << "Inversion 2 for break points failed: " << ierr;
605  return false;
606  }
607  theDerivatives[nMeasPerHit*l ][offsetPar+2*kbp ] = tempMSJacProj[0][0];
608  theDerivatives[nMeasPerHit*l ][offsetPar+2*kbp+1] = tempMSJacProj[0][1];
609  theDerivatives[nMeasPerHit*l+1][offsetPar+2*kbp ] = tempMSJacProj[1][0];
610  theDerivatives[nMeasPerHit*l+1][offsetPar+2*kbp+1] = tempMSJacProj[1][1];
611 
612  }
613 
614  }
615 
616  return true;
617 }
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
bool ReferenceTrajectory::addMaterialEffectsBrl ( const std::vector< AlgebraicMatrix > &  allJacobians,
const std::vector< AlgebraicMatrix > &  allProjections,
const std::vector< AlgebraicSymMatrix > &  allCurvChanges,
const std::vector< AlgebraicSymMatrix > &  allDeltaParaCovs,
const std::vector< AlgebraicMatrix > &  allLocalToCurv,
const GlobalTrajectoryParameters &  gtp 
)
protectedvirtual

internal methods to add material effects using broken lines (fine version)

Definition at line 621 of file ReferenceTrajectory.cc.

References gen::k, checklumidiff::l, GlobalTrajectoryParameters::momentum(), gen::n, ReferenceTrajectoryBase::nMeasPerHit, AlCaHLTBitMon_ParallelJobs::p, mathSSE::sqrt(), ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theInnerLocalToTrajectory, ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear, ReferenceTrajectoryBase::theMeasurementsCov, ReferenceTrajectoryBase::theNumberOfPars, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by construct().

627 {
628 //CHK: add material effects using broken lines
629 //fine: use exact Jacobians, all detectors
630 //broken lines: pair of offsets (u1,u2) = (xt,yt) (in curvilinear frame (q/p,lambda,phi,xt,yt)) at each layer
631 // scattering angles (alpha1,alpha2) = (cosLambda*dPhi, dLambda) (cosLambda cancels in Chi2)
632 // DU' = (dU'/dU)*DU + (dU'/dAlpha)*DAlpha + (dU'/dQbyp)*DQbyp (propagation of U)
633 // = J*DU + S*DAlpha + d*DQbyp
634 // => DAlpha = S^-1 (DU' - J*DU - d*DQbyp)
635 
636  int offsetPar = theNumberOfPars;
637  int offsetMeas = nMeasPerHit*allCurvlinJacobians.size();
638  int ierr = 0;
639 
640  GlobalVector p = gtp.momentum();
641  double cosLambda = sqrt((p.x()*p.x()+p.y()*p.y())/(p.x()*p.x()+p.y()*p.y()+p.z()*p.z()));
642 
643 // transformations Curvilinear <-> BrokenLines
644  AlgebraicMatrix QbypToCurv(5,1); // dCurv/dQbyp
645  QbypToCurv[0][0] = 1.; // dQbyp/dQbyp
646  AlgebraicMatrix AngleToCurv(5,2); // dCurv/dAlpha
647  AngleToCurv[1][1] = 1.; // dlambda/dalpha2
648  AngleToCurv[2][0] = 1./cosLambda; // dphi/dalpha1
649  AlgebraicMatrix CurvToAngle(2,5); // dAlpha/dCurv
650  CurvToAngle[1][1] = 1.; // dalpha2/dlambda
651  CurvToAngle[0][2] = cosLambda; // dalpha1/dphi
652  AlgebraicMatrix OffsetToCurv(5,2); // dCurv/dU
653  OffsetToCurv[3][0] = 1.; // dxt/du1
654  OffsetToCurv[4][1] = 1.; // dyt/du2
655  AlgebraicMatrix CurvToOffset(2,5); // dU/dCurv
656  CurvToOffset[0][3] = 1.; // du1/dxt
657  CurvToOffset[1][4] = 1.; // du2/dyt
658 
659 // transformations trajectory to components (Qbyp, U1, U2)
660  AlgebraicMatrix TrajToQbyp(1,5);
661  TrajToQbyp[0][0] = 1.;
662  AlgebraicMatrix TrajToOff1(2,5);
663  TrajToOff1[0][1] = 1.;
664  TrajToOff1[1][2] = 1.;
665  AlgebraicMatrix TrajToOff2(2,5);
666  TrajToOff2[0][3] = 1.;
667  TrajToOff2[1][4] = 1.;
668 
669  AlgebraicMatrix JacOffsetToAngleC, JacQbypToAngleC;
670  AlgebraicMatrix JacCurvToOffsetL, JacOffsetToOffsetL, JacAngleToOffsetL, JacQbypToOffsetL, JacOffsetToAngleL;
671  AlgebraicMatrix JacCurvToOffsetN, JacOffsetToOffsetN, JacAngleToOffsetN, JacQbypToOffsetN, JacOffsetToAngleN;
672 
673 // transformation from trajectory to curvilinear parameters
674 
675  JacCurvToOffsetN = CurvToOffset * allCurvlinJacobians[1]; // (dU'/dCurv') * (dCurv'/dCurv) @ 2nd point
676  JacOffsetToOffsetN = JacCurvToOffsetN * OffsetToCurv; // J: (dU'/dU) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
677  JacAngleToOffsetN = JacCurvToOffsetN * AngleToCurv; // S: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
678  JacQbypToOffsetN = JacCurvToOffsetN * QbypToCurv; // d: (dU'/dQbyp) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
679  JacOffsetToAngleN = JacAngleToOffsetN.inverse(ierr); // W
680  if (ierr) {
681  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
682  << "Inversion 1 for fine broken lines failed: " << ierr;
683  return false;
684  }
685  JacOffsetToAngleC = -(JacOffsetToAngleN * JacOffsetToOffsetN); // (dAlpha/dU)
686  JacQbypToAngleC = -(JacOffsetToAngleN * JacQbypToOffsetN); // (dAlpha/dQbyp)
687  // (dAlpha/dTraj) = (dAlpha/dQbyp) * (dQbyp/dTraj) + (dAlpha/dU1) * (dU1/dTraj) + (dAlpha/dU2) * (dU2/dTraj)
688  AlgebraicMatrix JacTrajToAngle = JacQbypToAngleC * TrajToQbyp + JacOffsetToAngleC * TrajToOff1 + JacOffsetToAngleN * TrajToOff2;
689  // (dCurv/dTraj) = (dCurv/dQbyp) * (dQbyp/dTraj) + (dCurv/dAlpha) * (dAlpha/dTraj) + (dCurv/dU) * (dU/dTraj)
690  theInnerTrajectoryToCurvilinear = QbypToCurv * TrajToQbyp + AngleToCurv * JacTrajToAngle + OffsetToCurv * TrajToOff1;
691  theInnerLocalToTrajectory = theInnerTrajectoryToCurvilinear.inverse(ierr) * allLocalToCurv[0];
692  if (ierr) {
693  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
694  << "Inversion 2 for fine broken lines failed: " << ierr;
695  return false;
696  }
697 
698  AlgebraicMatrix tempJacobian(allCurvatureChanges[0]);
699  AlgebraicSymMatrix tempMSCov;
700  AlgebraicSymMatrix tempMSCovProj;
701  AlgebraicMatrix tempJacL, tempJacN;
702  AlgebraicMatrix JacOffsetToMeas;
703 
704 // measurements from hits
705  for (unsigned int k = 0; k < allCurvlinJacobians.size(); ++k) {
706 // (dMeas/dU) = (dMeas/dLoc) * (dLoc/dCurv) * (dCurv/dU)
707  JacOffsetToMeas = (allProjections[k] * allLocalToCurv[k].inverse(ierr) ) * OffsetToCurv;
708  if (ierr) {
709  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
710  << "Inversion 3 for fine broken lines failed: " << ierr;
711  return false;
712  }
713  theDerivatives[nMeasPerHit*k ][offsetPar+2*k ] = JacOffsetToMeas[0][0];
714  theDerivatives[nMeasPerHit*k ][offsetPar+2*k+1] = JacOffsetToMeas[0][1];
715  theDerivatives[nMeasPerHit*k+1][offsetPar+2*k ] = JacOffsetToMeas[1][0];
716  theDerivatives[nMeasPerHit*k+1][offsetPar+2*k+1] = JacOffsetToMeas[1][1];
717  }
718 
719 // measurement of MS kink
720  for (unsigned int k = 1; k < allCurvlinJacobians.size()-1; ++k) {
721 // CHK
722  int iMsMeas = k-1;
723  int l = k-1; // last hit
724  int n = k+1; // next hit
725 
726 // amount of multiple scattering in layer k (angular error perp to direction)
727  tempMSCov = allDeltaParameterCovs[k].similarity(allLocalToCurv[k]);
728  tempMSCovProj = tempMSCov.similarity(CurvToAngle);
729  theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas ][offsetMeas+nMeasPerHit*iMsMeas ] = tempMSCovProj[1][1];
730  theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas+1][offsetMeas+nMeasPerHit*iMsMeas+1] = tempMSCovProj[0][0];
731 
732 // transformation matices for offsets ( l <- k -> n )
733  tempJacL = allCurvlinJacobians[k] * tempJacobian;
734  JacCurvToOffsetL = CurvToOffset * tempJacL.inverse(ierr); // (dU'/dCurv') * (dCurv'/dCurv) @ last point
735 
736  if (ierr) {
737  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
738  << "Inversion 4 for fine broken lines failed: " << ierr;
739  return false;
740  }
741  JacOffsetToOffsetL = JacCurvToOffsetL * OffsetToCurv; // J-: (dU'/dU) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
742  JacAngleToOffsetL = JacCurvToOffsetL * AngleToCurv; // S-: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
743  JacQbypToOffsetL = JacCurvToOffsetL * QbypToCurv; // d-: (dU'/dQbyp) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
744  JacOffsetToAngleL =-JacAngleToOffsetL.inverse(ierr); // W-
745  if (ierr) {
746  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
747  << "Inversion 5 for fine broken lines failed: " << ierr;
748  return false;
749  }
750  tempJacobian = tempJacobian * allCurvatureChanges[k];
751  tempJacN = allCurvlinJacobians[n] * tempJacobian;
752  JacCurvToOffsetN = CurvToOffset * tempJacN; // (dU'/dCurv') * (dCurv'/dCurv) @ next point
753  JacOffsetToOffsetN = JacCurvToOffsetN * OffsetToCurv; // J+: (dU'/dU) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
754  JacAngleToOffsetN = JacCurvToOffsetN * AngleToCurv; // S+: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
755  JacQbypToOffsetN = JacCurvToOffsetN * QbypToCurv; // d+: (dU'/dQbyp) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
756  JacOffsetToAngleN = JacAngleToOffsetN.inverse(ierr); // W+
757  if (ierr) {
758  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
759  << "Inversion 6 for fine broken lines failed: " << ierr;
760  return false;
761  }
762  JacOffsetToAngleC = -(JacOffsetToAngleL * JacOffsetToOffsetL + JacOffsetToAngleN * JacOffsetToOffsetN);
763  JacQbypToAngleC = -(JacOffsetToAngleL * JacQbypToOffsetL + JacOffsetToAngleN * JacQbypToOffsetN);
764 
765  // bending
766  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][ 0] = JacQbypToAngleC[0][0];
767  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][ 0] = JacQbypToAngleC[1][0];
768  // last layer
769  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*l ] = JacOffsetToAngleL[0][0];
770  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*l+1] = JacOffsetToAngleL[0][1];
771  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*l ] = JacOffsetToAngleL[1][0];
772  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*l+1] = JacOffsetToAngleL[1][1];
773  // current layer
774  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*k ] = JacOffsetToAngleC[0][0];
775  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*k+1] = JacOffsetToAngleC[0][1];
776  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*k ] = JacOffsetToAngleC[1][0];
777  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*k+1] = JacOffsetToAngleC[1][1];
778 
779  // next layer
780  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*n ] = JacOffsetToAngleN[0][0];
781  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*n+1] = JacOffsetToAngleN[0][1];
782  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*n ] = JacOffsetToAngleN[1][0];
783  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*n+1] = JacOffsetToAngleN[1][1];
784 
785  }
786 
787  return true;
788 }
ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear
AlgebraicMatrix theInnerTrajectoryToCurvilinear
Definition: ReferenceTrajectoryBase.h:251
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
ReferenceTrajectoryBase::theInnerLocalToTrajectory
AlgebraicMatrix theInnerLocalToTrajectory
Definition: ReferenceTrajectoryBase.h:253
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
bool ReferenceTrajectory::addMaterialEffectsBrl ( const std::vector< AlgebraicMatrix > &  allProjections,
const std::vector< AlgebraicSymMatrix > &  allDeltaParaCovs,
const std::vector< AlgebraicMatrix > &  allLocalToCurv,
const std::vector< double > &  allSteps,
const GlobalTrajectoryParameters &  gtp,
const double  minStep = 1.0 
)
protectedvirtual

internal methods to add material effects using broken lines (coarse version)

Definition at line 792 of file ReferenceTrajectory.cc.

References delta, plotBeamSpotDB::first, mps_fire::i, gen::k, checklumidiff::l, plotBeamSpotDB::last, mag(), GlobalTrajectoryParameters::magneticFieldInInverseGeV(), GlobalTrajectoryParameters::momentum(), gen::n, ReferenceTrajectoryBase::nMeasPerHit, AlCaHLTBitMon_ParallelJobs::p, GlobalTrajectoryParameters::position(), mathSSE::sqrt(), ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theInnerLocalToTrajectory, ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear, ReferenceTrajectoryBase::theMeasurementsCov, ReferenceTrajectoryBase::theNumberOfPars, ReferenceTrajectoryBase::theNumberOfVirtualMeas, ReferenceTrajectoryBase::theNumberOfVirtualPars, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

798 {
799 //CHK: add material effects using broken lines
800 //BrokenLinesCoarse: combine close by detectors,
801 // use approximate Jacobians from Steps (limit Qbyp -> 0),
802 // bending only in RPhi (B=(0,0,Bz)), no energy loss correction
803 
804  int offsetPar = theNumberOfPars;
805  int offsetMeas = nMeasPerHit*allSteps.size();
806  int ierr = 0;
807 
808  GlobalVector p = gtp.momentum();
809  double cosLambda = sqrt((p.x()*p.x()+p.y()*p.y())/(p.x()*p.x()+p.y()*p.y()+p.z()*p.z()));
810  double bFac = -gtp.magneticFieldInInverseGeV(gtp.position()).mag();
811 
812  // transformation from trajectory to curvilinear parameters at refTsos
813  double delta (1.0/allSteps[1]);
814  theInnerTrajectoryToCurvilinear[0][0] = 1;
815  theInnerTrajectoryToCurvilinear[1][2] = -delta;
816  theInnerTrajectoryToCurvilinear[1][4] = delta;
817  theInnerTrajectoryToCurvilinear[2][0] = -0.5*bFac/delta;
818  theInnerTrajectoryToCurvilinear[2][1] = -delta/cosLambda;
819  theInnerTrajectoryToCurvilinear[2][3] = delta/cosLambda;
820  theInnerTrajectoryToCurvilinear[3][1] = 1;
821  theInnerTrajectoryToCurvilinear[4][2] = 1;
822  theInnerLocalToTrajectory = theInnerTrajectoryToCurvilinear.inverse(ierr) * allLocalToCurv[0];
823  if (ierr) {
824  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
825  << "Inversion 1 for coarse broken lines failed: " << ierr;
826  return false;
827  }
828 
829  AlgebraicMatrix CurvToAngle(2,5); // dAlpha/dCurv
830  CurvToAngle[1][1] = 1.; // dalpha2/dlambda
831  CurvToAngle[0][2] = cosLambda; // dalpha1/dphi
832  AlgebraicMatrix OffsetToCurv(5,2); // dCurv/dU
833  OffsetToCurv[3][0] = 1.; // dxt/du1
834  OffsetToCurv[4][1] = 1.; // dyt/du2
835 
836  AlgebraicSymMatrix tempMSCov;
837  AlgebraicSymMatrix tempMSCovProj;
838  AlgebraicMatrix JacOffsetToMeas;
839 
840  // combine closeby detectors into single plane
841  std::vector<unsigned int> first(allSteps.size());
842  std::vector<unsigned int> last (allSteps.size());
843  std::vector<unsigned int> plane(allSteps.size());
844  std::vector<double> sPlane(allSteps.size());
845  unsigned int nPlane = 0;
846  double sTot = 0;
847 
848  for (unsigned int k = 1; k < allSteps.size(); ++k) {
849  sTot += allSteps[k];
850  if (fabs(allSteps[k])>minStep) { nPlane += 1; first[nPlane] = k; }
851  last[nPlane] = k;
852  plane[k] = nPlane;
853  sPlane[nPlane] += sTot;
854  }
855  if (nPlane < 2) return false; // pathological cases: need at least 2 planes
856 
857  theNumberOfVirtualPars = 2*(nPlane+1);
858  theNumberOfVirtualMeas = 2*(nPlane-1);// unsigned underflow for nPlane == 0...
859  for (unsigned int k = 0; k <= nPlane; ++k) { sPlane[k] /= (double) (last[k]-first[k]+1); }
860 
861  // measurements from hits
862  sTot = 0;
863  for (unsigned int k = 0; k < allSteps.size(); ++k) {
864  sTot += allSteps[k];
865 // (dMeas/dU) = (dMeas/dLoc) * (dLoc/dCurv) * (dCurv/dU)
866  JacOffsetToMeas = (allProjections[k] * allLocalToCurv[k].inverse(ierr) ) * OffsetToCurv;
867  if (ierr) {
868  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
869  << "Inversion 2 for coarse broken lines failed: " << ierr;
870  return false;
871  }
872 
873  unsigned int iPlane = plane[k];
874  if (last[iPlane] == first[iPlane])
875  { // single plane
876  theDerivatives[nMeasPerHit*k ][offsetPar+2*iPlane ] = JacOffsetToMeas[0][0];
877  theDerivatives[nMeasPerHit*k ][offsetPar+2*iPlane+1] = JacOffsetToMeas[0][1];
878  theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane ] = JacOffsetToMeas[1][0];
879  theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane+1] = JacOffsetToMeas[1][1];
880  } else
881  { // combined plane: (linear) interpolation
882  unsigned int jPlane; // neighbor plane for interpolation
883  if (fabs(sTot) < fabs(sPlane[iPlane])) { jPlane = (iPlane>0) ? iPlane - 1 : 1; }
884  else { jPlane = (iPlane<nPlane) ? iPlane + 1 : nPlane -1 ;}
885  // interpolation weights
886  double sDiff = sPlane[iPlane] - sPlane[jPlane];
887  double iFrac = (sTot - sPlane[jPlane]) / sDiff;
888  double jFrac = 1.0 - iFrac;
889  theDerivatives[nMeasPerHit*k ][offsetPar+2*iPlane ] = JacOffsetToMeas[0][0]*iFrac;
890  theDerivatives[nMeasPerHit*k ][offsetPar+2*iPlane+1] = JacOffsetToMeas[0][1]*iFrac;
891  theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane ] = JacOffsetToMeas[1][0]*iFrac;
892  theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane+1] = JacOffsetToMeas[1][1]*iFrac;
893  theDerivatives[nMeasPerHit*k ][offsetPar+2*jPlane ] = JacOffsetToMeas[0][0]*jFrac;
894  theDerivatives[nMeasPerHit*k ][offsetPar+2*jPlane+1] = JacOffsetToMeas[0][1]*jFrac;
895  theDerivatives[nMeasPerHit*k+1][offsetPar+2*jPlane ] = JacOffsetToMeas[1][0]*jFrac;
896  theDerivatives[nMeasPerHit*k+1][offsetPar+2*jPlane+1] = JacOffsetToMeas[1][1]*jFrac;
897  // 2nd order neglected
898  // theDerivatives[nMeasPerHit*k ][ 0] = -0.5*bFac*sDiff*iFrac*sDiff*jFrac*cosLambda;
899  }
900  }
901 
902 // measurement of MS kink
903  for (unsigned int i = 1; i < nPlane; ++i) {
904 // CHK
905  int iMsMeas = i-1;
906  int l = i-1; // last hit
907  int n = i+1; // next hit
908 
909 // amount of multiple scattering in plane k
910  for (unsigned int k = first[i]; k <= last[i]; ++k) {
911  tempMSCov = allDeltaParameterCovs[k].similarity(allLocalToCurv[k]);
912  tempMSCovProj = tempMSCov.similarity(CurvToAngle);
913  theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas ][offsetMeas+nMeasPerHit*iMsMeas ] += tempMSCovProj[0][0];
914  theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas+1][offsetMeas+nMeasPerHit*iMsMeas+1] += tempMSCovProj[1][1];
915  }
916 // broken line measurements for layer k, correlations between both planes neglected
917  double stepK = sPlane[i] - sPlane[l];
918  double stepN = sPlane[n] - sPlane[i];
919  double deltaK (1.0/stepK);
920  double deltaN (1.0/stepN);
921  // bending (only in RPhi)
922  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][ 0] = -0.5*bFac*(stepK+stepN)*cosLambda;
923  // last layer
924  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*l ] = deltaK;
925  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*l+1] = deltaK;
926  // current layer
927  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*i ] = -(deltaK + deltaN);
928  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*i+1] = -(deltaK + deltaN);
929  // next layer
930  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas ][offsetPar+2*n ] = deltaN;
931  theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*n+1] = deltaN;
932  }
933 
934  return true;
935 }
delta
dbl * delta
Definition: mlp_gen.cc:36
ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear
AlgebraicMatrix theInnerTrajectoryToCurvilinear
Definition: ReferenceTrajectoryBase.h:251
mps_fire.i
i
Definition: mps_fire.py:156
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
ReferenceTrajectoryBase::theInnerLocalToTrajectory
AlgebraicMatrix theInnerLocalToTrajectory
Definition: ReferenceTrajectoryBase.h:253
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:147
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
GlobalTrajectoryParameters::magneticFieldInInverseGeV
GlobalVector magneticFieldInInverseGeV(const GlobalPoint &x) const
Definition: GlobalTrajectoryParameters.cc:64
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
bool ReferenceTrajectory::addMaterialEffectsCov ( const std::vector< AlgebraicMatrix > &  allJacobians,
const std::vector< AlgebraicMatrix > &  allProjections,
const std::vector< AlgebraicSymMatrix > &  allCurvChanges,
const std::vector< AlgebraicSymMatrix > &  allDeltaParaCovs 
)
protectedvirtual

internal method to add material effects to measurments covariance matrix

Definition at line 488 of file ReferenceTrajectory.cc.

References gen::k, checklumidiff::l, ReferenceTrajectoryBase::nMeasPerHit, and ReferenceTrajectoryBase::theMeasurementsCov.

Referenced by construct().

492 {
493  // the uncertainty due to multiple scattering is 'transferred' to the error matrix of the hits
494 
495  // GF: Needs update once hit dimension is not hardcoded as nMeasPerHit!
496 
497  AlgebraicSymMatrix materialEffectsCov(nMeasPerHit * allJacobians.size(), 0);
498 
499  // additional covariance-matrix of the measurements due to material-effects (single measurement)
500  AlgebraicSymMatrix deltaMaterialEffectsCov;
501 
502  // additional covariance-matrix of the parameters due to material-effects
503  AlgebraicSymMatrix paramMaterialEffectsCov(allDeltaParameterCovs[0]); //initialization
504  // error-propagation to state after energy loss
505  //GFback paramMaterialEffectsCov = paramMaterialEffectsCov.similarity(allCurvatureChanges[0]);
506 
507  AlgebraicMatrix tempParameterCov;
508  AlgebraicMatrix tempMeasurementCov;
509 
510  for (unsigned int k = 1; k < allJacobians.size(); ++k) {
511  // error-propagation to next layer
512  paramMaterialEffectsCov = paramMaterialEffectsCov.similarity(allJacobians[k]);
513 
514  // get dependences for the measurements
515  deltaMaterialEffectsCov = paramMaterialEffectsCov.similarity(allProjections[k]);
516  materialEffectsCov[nMeasPerHit*k ][nMeasPerHit*k ] = deltaMaterialEffectsCov[0][0];
517  materialEffectsCov[nMeasPerHit*k ][nMeasPerHit*k+1] = deltaMaterialEffectsCov[0][1];
518  materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*k ] = deltaMaterialEffectsCov[1][0];
519  materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*k+1] = deltaMaterialEffectsCov[1][1];
520 
521  // GFback add uncertainties for the following layers due to scattering at this layer
522  paramMaterialEffectsCov += allDeltaParameterCovs[k];
523  // end GFback
524  tempParameterCov = paramMaterialEffectsCov;
525 
526  // compute "inter-layer-dependencies"
527  for (unsigned int l = k+1; l < allJacobians.size(); ++l) {
528  tempParameterCov = allJacobians[l] * allCurvatureChanges[l] * tempParameterCov;
529  tempMeasurementCov = allProjections[l] * tempParameterCov * allProjections[k].T();
530 
531  materialEffectsCov[nMeasPerHit*l][nMeasPerHit*k] = tempMeasurementCov[0][0];
532  materialEffectsCov[nMeasPerHit*k][nMeasPerHit*l] = tempMeasurementCov[0][0];
533 
534  materialEffectsCov[nMeasPerHit*l][nMeasPerHit*k+1] = tempMeasurementCov[0][1];
535  materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*l] = tempMeasurementCov[0][1];
536 
537  materialEffectsCov[nMeasPerHit*l+1][nMeasPerHit*k] = tempMeasurementCov[1][0];
538  materialEffectsCov[nMeasPerHit*k][nMeasPerHit*l+1] = tempMeasurementCov[1][0];
539 
540  materialEffectsCov[nMeasPerHit*l+1][nMeasPerHit*k+1] = tempMeasurementCov[1][1];
541  materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*l+1] = tempMeasurementCov[1][1];
542  }
543  // add uncertainties for the following layers due to scattering at this layer
544  // GFback paramMaterialEffectsCov += allDeltaParameterCovs[k];
545  // error-propagation to state after energy loss
546  paramMaterialEffectsCov = paramMaterialEffectsCov.similarity(allCurvatureChanges[k]);
547 
548  }
549  theMeasurementsCov += materialEffectsCov;
550 
551  return true; // cannot fail
552 }
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
bool ReferenceTrajectory::addMaterialEffectsCurvlinGbl ( const std::vector< AlgebraicMatrix > &  allJacobians,
const std::vector< AlgebraicMatrix > &  allProjections,
const std::vector< AlgebraicSymMatrix > &  allCurvChanges,
const std::vector< AlgebraicSymMatrix > &  allDeltaParaCovs,
const std::vector< AlgebraicMatrix > &  allLocalToCurv 
)
protectedvirtual

internal methods to add material effects using broken lines (fine version, curvilinear system)

Definition at line 1023 of file ReferenceTrajectory.cc.

References funct::abs(), gbl::GblPoint::addMeasurement(), gbl::GblPoint::addScatterer(), allowZeroMaterial_, clhep2eigen(), geometryDiff::epsilon, Exception, infinity, gen::k, ReferenceTrajectoryBase::theGblInput, ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

1028 {
1029 //CHK: add material effects using general broken lines
1030 // local track parameters are defined in the curvilinear system
1031 
1032  const double minPrec = 1.0; // minimum precision to use measurement (reject measurements in strip direction)
1033  int ierr = 0;
1034 
1035  AlgebraicMatrix OffsetToCurv(5,2); // dCurv/dU
1036  OffsetToCurv[3][0] = 1.; // dxt/du1
1037  OffsetToCurv[4][1] = 1.; // dyt/du2
1038 
1039  AlgebraicMatrix JacOffsetToMeas, tempMSCov;
1040 
1041  // GBL uses Eigen matrices as interface
1042  Eigen::Matrix2d covariance, proLocalToMeas;
1043  Matrix5d jacPointToPoint, firstLocalToCurv;
1044  Eigen::Vector2d measurement, measPrecDiag, scatPrecDiag;
1045  auto scatterer = Eigen::Vector2d::Zero();
1046 
1047 // measurements and scatterers from hits
1048  unsigned int numHits = allCurvlinJacobians.size();
1049  std::vector<GblPoint> GblPointList;
1050  GblPointList.reserve(numHits);
1051  for (unsigned int k = 0; k < numHits; ++k) {
1052 // (dMeas/dU) = (dMeas/dLoc) * (dLoc/dCurv) * (dCurv/dU)
1053  JacOffsetToMeas = (allProjections[k] * allLocalToCurv[k].inverse(ierr) ) * OffsetToCurv;
1054  if (ierr) {
1055  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsGbl"
1056  << "Inversion 1 for general broken lines failed: " << ierr;
1057  return false;
1058  }
1059 
1060  // GBL point to point jacobian
1061  clhep2eigen(allCurvlinJacobians[k] * allCurvatureChanges[k], jacPointToPoint);
1062 
1063  // GBL point
1064  GblPoint aGblPoint( jacPointToPoint );
1065 
1066  // GBL projection from local to measurement system
1067  clhep2eigen(JacOffsetToMeas, proLocalToMeas);
1068 
1069  // GBL measurement (residuum to initial trajectory)
1070  clhep2eigen(theMeasurements.sub(2*k+1,2*k+2) - theTrajectoryPositions.sub(2*k+1,2*k+2), measurement);
1071 
1072  // GBL measurement covariance matrix
1073  clhep2eigen(theMeasurementsCov.sub(2*k+1,2*k+2), covariance);
1074 
1075  // GBL add measurement to point
1076  if (std::abs(covariance(0,1)) < std::numeric_limits<double>::epsilon()) {
1077  // covariance matrix is diagonal, independent measurements
1078  for (unsigned int row = 0; row < 2; ++row) {
1079  measPrecDiag(row) = ( 0. < covariance(row,row) ? 1.0/covariance(row,row) : 0. );
1080  }
1081  aGblPoint.addMeasurement(proLocalToMeas, measurement, measPrecDiag, minPrec);
1082  } else
1083  {
1084  // covariance matrix needs diagonalization
1085  aGblPoint.addMeasurement(proLocalToMeas, measurement, covariance.inverse(), minPrec);
1086  }
1087 
1088  // GBL multiple scattering (diagonal matrix in curvilinear system)
1089  tempMSCov = allDeltaParameterCovs[k].similarity(allLocalToCurv[k]);
1090  for (unsigned int row = 0; row < 2; ++row) {
1091  scatPrecDiag(row) = 1.0/tempMSCov[row+1][row+1];
1092  }
1093 
1094  // check for singularity
1095  bool singularCovariance{false};
1096  for (int row = 0; row < scatPrecDiag.rows(); ++row) {
1097  if (!(scatPrecDiag[row] < std::numeric_limits<double>::infinity())) {
1098  singularCovariance = true;
1099  break;
1100  }
1101  }
1102  if (singularCovariance && !allowZeroMaterial_) {
1103  throw cms::Exception("Alignment")
1104  << "@SUB=ReferenceTrajectory::addMaterialEffectsCurvlinGbl"
1105  << "\nEncountered singular scatter covariance-matrix without allowing "
1106  << "for zero material.";
1107  }
1108 
1109  // GBL add scatterer to point
1110  aGblPoint.addScatterer(scatterer, scatPrecDiag);
1111 
1112  // add point to list
1113  GblPointList.push_back( aGblPoint );
1114  }
1115  // add list of points and transformation local to fit (=curvilinear) system at first hit
1116  clhep2eigen(allLocalToCurv[0], firstLocalToCurv);
1117  theGblInput.push_back(std::make_pair(GblPointList, firstLocalToCurv));
1118 
1119  return true;
1120 }
ReferenceTrajectory::clhep2eigen
void clhep2eigen(const AlgebraicVector &in, Eigen::MatrixBase< Derived > &out)
Definition: ReferenceTrajectory.cc:1124
ReferenceTrajectoryBase::theGblInput
std::vector< std::pair< std::vector< gbl::GblPoint >, Eigen::MatrixXd > > theGblInput
Definition: ReferenceTrajectoryBase.h:255
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
infinity
const double infinity
Definition: CSCChamberFitter.cc:10
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:66
gbl::Matrix5d
Eigen::Matrix< double, 5, 5 > Matrix5d
Definition: GblData.h:43
bool ReferenceTrajectory::addMaterialEffectsLocalGbl ( const std::vector< AlgebraicMatrix > &  allJacobians,
const std::vector< AlgebraicMatrix > &  allProjections,
const std::vector< AlgebraicSymMatrix > &  allCurvatureChanges,
const std::vector< AlgebraicSymMatrix > &  allDeltaParameterCovs 
)
protectedvirtual

internal methods to add material effects using broken lines (fine version, local system)

Definition at line 939 of file ReferenceTrajectory.cc.

References funct::abs(), gbl::GblPoint::addMeasurement(), gbl::GblPoint::addScatterer(), allowZeroMaterial_, clhep2eigen(), geometryDiff::epsilon, Exception, gen::k, ReferenceTrajectoryBase::theGblInput, ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

943 {
944 //CHK: add material effects using general broken lines, no initial kinks
945 // local track parameters are defined in the TSO system
946 
947  const double minPrec = 1.0; // minimum precision to use measurement (reject measurements in strip direction)
948 
949  AlgebraicMatrix OffsetToLocal(5,2); // dLocal/dU
950  OffsetToLocal[3][0] = 1.;
951  OffsetToLocal[4][1] = 1.;
952  AlgebraicMatrix SlopeToLocal(5,2); // dLocal/dU'
953  SlopeToLocal[1][0] = 1.;
954  SlopeToLocal[2][1] = 1.;
955 
956  // GBL uses Eigen matrices as interface
957  Eigen::Matrix2d covariance, scatPrecision, proLocalToMeas;
958  Matrix5d jacPointToPoint;
959  auto identity = Matrix5d::Identity();
960  Eigen::Vector2d measurement, measPrecDiag;
961  auto scatterer = Eigen::Vector2d::Zero();
962 
963  //bool initialKinks = (allCurvlinKinks.size()>0);
964 
965 // measurements and scatterers from hits
966  unsigned int numHits = allJacobians.size();
967  std::vector<GblPoint> GblPointList;
968  GblPointList.reserve(numHits);
969  for (unsigned int k = 0; k < numHits; ++k) {
970 
971  // GBL point to point jacobian
972  clhep2eigen(allJacobians[k] * allCurvatureChanges[k], jacPointToPoint);
973 
974  // GBL point
975  GblPoint aGblPoint( jacPointToPoint );
976 
977  // GBL projection from local to measurement system
978  clhep2eigen(allProjections[k] * OffsetToLocal, proLocalToMeas);
979 
980  // GBL measurement (residuum to initial trajectory)
981  clhep2eigen(theMeasurements.sub(2*k+1,2*k+2) - theTrajectoryPositions.sub(2*k+1,2*k+2), measurement);
982 
983  // GBL measurement covariance matrix
984  clhep2eigen(theMeasurementsCov.sub(2*k+1,2*k+2), covariance);
985 
986  // GBL add measurement to point
987  if (std::abs(covariance(0,1)) < std::numeric_limits<double>::epsilon()) {
988  // covariance matrix is diagonal, independent measurements
989  for (unsigned int row = 0; row < 2; ++row) {
990  measPrecDiag(row) = ( 0. < covariance(row,row) ? 1.0/covariance(row,row) : 0. );
991  }
992  aGblPoint.addMeasurement(proLocalToMeas, measurement, measPrecDiag, minPrec);
993  } else
994  {
995  // covariance matrix needs diagonalization
996  aGblPoint.addMeasurement(proLocalToMeas, measurement, covariance.inverse(), minPrec);
997  }
998 
999  // GBL multiple scattering (full matrix in local system)
1000  clhep2eigen(allDeltaParameterCovs[k].similarityT(SlopeToLocal), scatPrecision);
1001  if (!(scatPrecision.colPivHouseholderQr().isInvertible())) {
1002  if (!allowZeroMaterial_) {
1003  throw cms::Exception("Alignment")
1004  << "@SUB=ReferenceTrajectory::addMaterialEffectsLocalGbl"
1005  << "\nEncountered singular scatter covariance-matrix without allowing "
1006  << "for zero material.";
1007  }
1008  } else {
1009  // GBL add scatterer to point
1010  aGblPoint.addScatterer(scatterer, scatPrecision.inverse());
1011  }
1012  // add point to list
1013  GblPointList.push_back( aGblPoint );
1014  }
1015  // add list of points and transformation local to fit (=local) system at first hit
1016  theGblInput.push_back(std::make_pair(GblPointList, identity));
1017 
1018  return true;
1019 }
ReferenceTrajectory::clhep2eigen
void clhep2eigen(const AlgebraicVector &in, Eigen::MatrixBase< Derived > &out)
Definition: ReferenceTrajectory.cc:1124
ReferenceTrajectoryBase::theGblInput
std::vector< std::pair< std::vector< gbl::GblPoint >, Eigen::MatrixXd > > theGblInput
Definition: ReferenceTrajectoryBase.h:255
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:66
gbl::Matrix5d
Eigen::Matrix< double, 5, 5 > Matrix5d
Definition: GblData.h:43
template<typename Derived >
void ReferenceTrajectory::clhep2eigen ( const AlgebraicVector &  in,
Eigen::MatrixBase< Derived > &  out 
)
private

Definition at line 1124 of file ReferenceTrajectory.cc.

References MillePedeFileConverter_cfg::out.

Referenced by addMaterialEffectsCurvlinGbl(), and addMaterialEffectsLocalGbl().

1125  {
1126  static_assert(Derived::ColsAtCompileTime == 1,
1127  "clhep2eigen: 'out' must be of vector type");
1128  for (int row = 0; row < in.num_row(); ++row) {
1129  out(row) = in[row];
1130  }
1131 }
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
template<typename Derived >
void ReferenceTrajectory::clhep2eigen ( const AlgebraicMatrix &  in,
Eigen::MatrixBase< Derived > &  out 
)
private

Definition at line 1134 of file ReferenceTrajectory.cc.

References cuy::col, and MillePedeFileConverter_cfg::out.

1135  {
1136  for (int row = 0; row < in.num_row(); ++row) {
1137  for (int col = 0; col < in.num_col(); ++col) {
1138  out(row, col) = in[row][col];
1139  }
1140  }
1141 }
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
cuy.col
col
Definition: cuy.py:1008
template<typename Derived >
void ReferenceTrajectory::clhep2eigen ( const AlgebraicSymMatrix &  in,
Eigen::MatrixBase< Derived > &  out 
)
private

Definition at line 1144 of file ReferenceTrajectory.cc.

References cuy::col, getHitProjectionMatrix(), and MillePedeFileConverter_cfg::out.

1145  {
1146  for (int row = 0; row < in.num_row(); ++row) {
1147  for (int col = 0; col < in.num_col(); ++col) {
1148  out(row, col) = in[row][col];
1149  }
1150  }
1151 }
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
cuy.col
col
Definition: cuy.py:1008
virtual ReferenceTrajectory* ReferenceTrajectory::clone ( void  ) const
inlinevirtual

Implements ReferenceTrajectoryBase.

Reimplemented in BzeroReferenceTrajectory.

Definition at line 75 of file ReferenceTrajectory.h.

References ecalDrivenElectronSeedsParameters_cff::beamSpot, and ResonanceBuilder::mass.

75 { return new ReferenceTrajectory(*this); }
ReferenceTrajectory::ReferenceTrajectory
ReferenceTrajectory(const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, const MagneticField *magField, const reco::BeamSpot &beamSpot, const ReferenceTrajectoryBase::Config &config)
Definition: ReferenceTrajectory.cc:53
bool ReferenceTrajectory::construct ( const TrajectoryStateOnSurface &  referenceTsos,
const TransientTrackingRecHit::ConstRecHitContainer &  recHits,
const MagneticField *  magField,
const reco::BeamSpot &  beamSpot 
)
protectedvirtual

internal method to calculate members

Definition at line 111 of file ReferenceTrajectory.cc.

References addMaterialEffectsBp(), addMaterialEffectsBrl(), addMaterialEffectsCov(), addMaterialEffectsCurvlinGbl(), addMaterialEffectsLocalGbl(), alongMomentum, AnalyticalPropagator_cfi::AnalyticalPropagator, ReferenceTrajectoryBase::breakPoints, ReferenceTrajectoryBase::brokenLinesCoarse, ReferenceTrajectoryBase::brokenLinesFine, Plane::build(), ReferenceTrajectoryBase::combined, createUpdator(), Vector3DBase< T, FrameTag >::cross(), ReferenceTrajectoryBase::curvlinGBL, reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), ReferenceTrajectoryBase::energyLoss, fillDerivatives(), fillMeasurementAndError(), fillTrajectoryPositions(), TrajectoryStateOnSurface::freeState(), getHitProjectionMatrix(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::hasError(), TrajectoryStateOnSurface::localError(), ReferenceTrajectoryBase::localGBL, TrajectoryStateOnSurface::localParameters(), mass_, materialEffects_, LocalTrajectoryError::matrix(), FreeTrajectoryState::momentum(), ReferenceTrajectoryBase::multipleScattering, ReferenceTrajectoryBase::nMeasPerHit, MagneticField::nominalValue(), ReferenceTrajectoryBase::none, PV3DBase< T, PVType, FrameType >::phi(), propagate(), PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi::propagator, propDir_, idealTransformation::rotation, LocalTrajectoryParameters::signedInverseMomentum(), GeomDet::surface(), TrajectoryStateOnSurface::surface(), surfaceSide(), ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theInnerLocalToTrajectory, ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear, ReferenceTrajectoryBase::theNomField, ReferenceTrajectoryBase::theNumberOfHits, ReferenceTrajectoryBase::theNumberOfVirtualPars, ReferenceTrajectoryBase::theParameters, ReferenceTrajectoryBase::theRecHits, ReferenceTrajectoryBase::theTrajectoryPositionCov, ReferenceTrajectoryBase::theTsosVec, TSCBLBuilderNoMaterial_cfi::TSCBLBuilderNoMaterial, useBeamSpot_, ReferenceTrajectoryBase::useRecHit(), reco::BeamSpot::x0(), reco::BeamSpot::y0(), and reco::BeamSpot::z0().

Referenced by BzeroReferenceTrajectory::BzeroReferenceTrajectory(), and ReferenceTrajectory().

115 {
116  TrajectoryStateOnSurface theRefTsos = refTsos;
117 
118  const SurfaceSide surfaceSide = this->surfaceSide(propDir_);
119  // auto_ptr to avoid memory leaks in case of not reaching delete at end of method:
120  std::auto_ptr<MaterialEffectsUpdator> aMaterialEffectsUpdator
121  (this->createUpdator(materialEffects_, mass_));
122  if (!aMaterialEffectsUpdator.get()) return false; // empty auto_ptr
123 
124  AlgebraicMatrix fullJacobian(theParameters.num_row(), theParameters.num_row());
125  std::vector<AlgebraicMatrix> allJacobians;
126  allJacobians.reserve(theNumberOfHits);
127 
128  TransientTrackingRecHit::ConstRecHitPointer previousHitPtr;
129  TrajectoryStateOnSurface previousTsos;
130  AlgebraicSymMatrix previousChangeInCurvature(theParameters.num_row(), 1);
131  std::vector<AlgebraicSymMatrix> allCurvatureChanges;
132  allCurvatureChanges.reserve(theNumberOfHits);
133 
134  const LocalTrajectoryError zeroErrors(0., 0., 0., 0., 0.);
135 
136  std::vector<AlgebraicMatrix> allProjections;
137  allProjections.reserve(theNumberOfHits);
138  std::vector<AlgebraicSymMatrix> allDeltaParameterCovs;
139  allDeltaParameterCovs.reserve(theNumberOfHits);
140 
141  // CHK
142  std::vector<AlgebraicMatrix> allLocalToCurv;
143  allLocalToCurv.reserve(theNumberOfHits);
144  std::vector<double> allSteps;
145  allSteps.reserve(theNumberOfHits);
146  std::vector<AlgebraicMatrix> allCurvlinJacobians;
147  allCurvlinJacobians.reserve(theNumberOfHits);
148 
149  AlgebraicMatrix firstCurvlinJacobian(5, 5, 1);
150 
151  unsigned int iRow = 0;
152 
153  theNomField = magField->nominalValue(); // nominal magnetic field in kGauss
154  // local storage vector of all rechits (including rechit for beam spot in case it is used)
155  TransientTrackingRecHit::ConstRecHitContainer allRecHits;
156 
157  if (useBeamSpot_ && propDir_==alongMomentum) {
158 
159  GlobalPoint bs(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
160 
161  TrajectoryStateClosestToBeamLine tsctbl(TSCBLBuilderNoMaterial()(*(refTsos.freeState()), beamSpot));
162  if (!tsctbl.isValid()) {
163  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::construct"
164  << "TrajectoryStateClostestToBeamLine invalid. Skip track.";
165  return false;
166  }
167 
168  FreeTrajectoryState pcaFts = tsctbl.trackStateAtPCA();
169  GlobalVector bd(beamSpot.dxdz(), beamSpot.dydz(), 1.0);
170 
171  //propagation FIXME: Should use same propagator everywhere...
172  AnalyticalPropagator propagator(magField);
173  std::pair< TrajectoryStateOnSurface, double > tsosWithPath =
174  propagator.propagateWithPath(pcaFts, refTsos.surface());
175 
176  if (!tsosWithPath.first.isValid()) return false;
177 
178  GlobalVector momDir(pcaFts.momentum());
179  GlobalVector perpDir(bd.cross(momDir));
180  Plane::RotationType rotation(perpDir, bd);
181 
182  BeamSpotGeomDet * bsGeom = new BeamSpotGeomDet(Plane::build(bs, rotation));
183 
184  // There is also a constructor taking the magentic field. Use this one instead?
185  theRefTsos = TrajectoryStateOnSurface(pcaFts, bsGeom->surface());
186 
187  TransientTrackingRecHit::ConstRecHitPointer bsRecHit(new BeamSpotTransientTrackingRecHit(beamSpot,
188  bsGeom,
189  theRefTsos.freeState()->momentum().phi()));
190  allRecHits.push_back(bsRecHit);
191 
192  }
193 
194  // copy all rechits to the local storage vector
195  TransientTrackingRecHit::ConstRecHitContainer::const_iterator itRecHit;
196  for ( itRecHit = recHits.begin(); itRecHit != recHits.end(); ++itRecHit ) {
197  const TransientTrackingRecHit::ConstRecHitPointer &hitPtr = *itRecHit;
198  allRecHits.push_back(hitPtr);
199  }
200 
201  for ( itRecHit = allRecHits.begin(); itRecHit != allRecHits.end(); ++itRecHit ) {
202 
203  const TransientTrackingRecHit::ConstRecHitPointer &hitPtr = *itRecHit;
204  theRecHits.push_back(hitPtr);
205 
206  if (0 == iRow) {
207 
208  // compute the derivatives of the reference-track's parameters w.r.t. the initial ones
209  // derivative of the initial reference-track parameters w.r.t. themselves is of course the identity
210  fullJacobian = AlgebraicMatrix(theParameters.num_row(), theParameters.num_row(), 1);
211  allJacobians.push_back(fullJacobian);
212  theTsosVec.push_back(theRefTsos);
213  const JacobianLocalToCurvilinear startTrafo(hitPtr->det()->surface(), theRefTsos.localParameters(), *magField);
214  const AlgebraicMatrix localToCurvilinear = asHepMatrix<5>(startTrafo.jacobian());
215  if (materialEffects_ <= breakPoints) {
216  theInnerTrajectoryToCurvilinear = asHepMatrix<5>(startTrafo.jacobian());
217  theInnerLocalToTrajectory = AlgebraicMatrix(5, 5, 1);
218  }
219  allLocalToCurv.push_back(localToCurvilinear);
220  allSteps.push_back(0.);
221  allCurvlinJacobians.push_back(firstCurvlinJacobian);
222 
223  } else {
224 
225  AlgebraicMatrix nextJacobian;
226  AlgebraicMatrix nextCurvlinJacobian;
227  double nextStep = 0.;
228  TrajectoryStateOnSurface nextTsos;
229 
230  if (!this->propagate(previousHitPtr->det()->surface(), previousTsos,
231  hitPtr->det()->surface(), nextTsos,
232  nextJacobian, nextCurvlinJacobian, nextStep, magField)) {
233  return false; // stop if problem...// no delete aMaterialEffectsUpdator needed
234  }
235 
236  allJacobians.push_back(nextJacobian);
237  fullJacobian = nextJacobian * previousChangeInCurvature * fullJacobian;
238  theTsosVec.push_back(nextTsos);
239 
240  const JacobianLocalToCurvilinear startTrafo(hitPtr->det()->surface(), nextTsos.localParameters(), *magField);
241  const AlgebraicMatrix localToCurvilinear = asHepMatrix<5>(startTrafo.jacobian());
242  allLocalToCurv.push_back(localToCurvilinear);
243  if (nextStep == 0.) {
244  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::construct"
245  << "step 0. from id " << previousHitPtr->geographicalId()
246  << " to " << hitPtr->det()->geographicalId() << ".";
247  // brokenLinesFine will not work, brokenLinesCoarse combines close by layers
248  if (materialEffects_ == brokenLinesFine) {
249  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::construct" << "Skip track.";
250  return false;
251  }
252  }
253  allSteps.push_back(nextStep);
254  allCurvlinJacobians.push_back(nextCurvlinJacobian);
255 
256  }
257 
258  // take material effects into account. since trajectory-state is constructed with errors equal zero,
259  // the updated state contains only the uncertainties due to interactions in the current layer.
260  const TrajectoryStateOnSurface tmpTsos(theTsosVec.back().localParameters(), zeroErrors,
261  theTsosVec.back().surface(), magField, surfaceSide);
262  const TrajectoryStateOnSurface updatedTsos = aMaterialEffectsUpdator->updateState(tmpTsos, propDir_);
263 
264  if ( !updatedTsos.isValid() ) return false;// no delete aMaterialEffectsUpdator needed
265 
266  if ( theTsosVec.back().localParameters().charge() )
267  {
268  previousChangeInCurvature[0][0] = updatedTsos.localParameters().signedInverseMomentum()
269  / theTsosVec.back().localParameters().signedInverseMomentum();
270  }
271 
272  // get multiple-scattering covariance-matrix
273  allDeltaParameterCovs.push_back( asHepMatrix<5>(updatedTsos.localError().matrix()) );
274  allCurvatureChanges.push_back(previousChangeInCurvature);
275 
276  // projection-matrix tsos-parameters -> measurement-coordinates
277  allProjections.push_back(this->getHitProjectionMatrix(hitPtr));
278  // set start-parameters for next propagation. trajectory-state without error
279  // - no error propagation needed here.
280  previousHitPtr = hitPtr;
281  previousTsos = TrajectoryStateOnSurface(updatedTsos.globalParameters(),
282  updatedTsos.surface(), surfaceSide);
283 
284  if (materialEffects_ < brokenLinesCoarse) {
285  this->fillDerivatives(allProjections.back(), fullJacobian, iRow);
286  }
287 
288  AlgebraicVector mixedLocalParams = asHepVector<5>(theTsosVec.back().localParameters().mixedFormatVector());
289  this->fillTrajectoryPositions(allProjections.back(), mixedLocalParams, iRow);
290  if ( useRecHit( hitPtr ) ) this->fillMeasurementAndError(hitPtr, iRow, updatedTsos);
291 
292  iRow += nMeasPerHit;
293  } // end of loop on hits
294 
295  bool msOK = true;
296  switch (materialEffects_) {
297  case none:
298  break;
299  case multipleScattering:
300  case energyLoss:
301  case combined:
302  msOK = this->addMaterialEffectsCov(allJacobians, allProjections, allCurvatureChanges,
303  allDeltaParameterCovs);
304  break;
305  case breakPoints:
306  msOK = this->addMaterialEffectsBp(allJacobians, allProjections, allCurvatureChanges,
307  allDeltaParameterCovs, allLocalToCurv);
308  break;
309  case brokenLinesCoarse:
310  msOK = this->addMaterialEffectsBrl(allProjections, allDeltaParameterCovs, allLocalToCurv,
311  allSteps, refTsos.globalParameters());
312  break;
313  case brokenLinesFine:
314  msOK = this->addMaterialEffectsBrl(allCurvlinJacobians, allProjections, allCurvatureChanges,
315  allDeltaParameterCovs, allLocalToCurv, refTsos.globalParameters());
316  break;
317  case localGBL:
318  msOK = this->addMaterialEffectsLocalGbl(allJacobians, allProjections, allCurvatureChanges,
319  allDeltaParameterCovs);
320  break;
321  case curvlinGBL:
322  msOK = this->addMaterialEffectsCurvlinGbl(allCurvlinJacobians, allProjections, allCurvatureChanges,
323  allDeltaParameterCovs, allLocalToCurv);
324  }
325  if (!msOK) return false;
326 
327  if (refTsos.hasError()) {
328  AlgebraicSymMatrix parameterCov = asHepMatrix<5>(refTsos.localError().matrix());
329  AlgebraicMatrix parDeriv;
330  if (theNumberOfVirtualPars>0) {
331  parDeriv = theDerivatives.sub( 1, nMeasPerHit*allJacobians.size(), 1, theParameters.num_row() );
332  } else {
333  parDeriv = theDerivatives;
334  }
335  theTrajectoryPositionCov = parameterCov.similarity(parDeriv);
336  } else {
337  theTrajectoryPositionCov = AlgebraicSymMatrix(theDerivatives.num_row(), 1);
338  }
339 
340  return true;
341 }
ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear
AlgebraicMatrix theInnerTrajectoryToCurvilinear
Definition: ReferenceTrajectoryBase.h:251
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:68
ReferenceTrajectory::getHitProjectionMatrix
AlgebraicMatrix getHitProjectionMatrix(const TransientTrackingRecHit::ConstRecHitPointer &recHit) const
Definition: ReferenceTrajectory.cc:1157
ReferenceTrajectory::addMaterialEffectsBrl
virtual bool addMaterialEffectsBrl(const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv, const GlobalTrajectoryParameters &gtp)
Definition: ReferenceTrajectory.cc:621
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:107
ReferenceTrajectoryBase::theInnerLocalToTrajectory
AlgebraicMatrix theInnerLocalToTrajectory
Definition: ReferenceTrajectoryBase.h:253
MagneticField::nominalValue
int nominalValue() const
The nominal field value for this map in kGauss.
Definition: MagneticField.h:56
ReferenceTrajectory::fillMeasurementAndError
virtual void fillMeasurementAndError(const TransientTrackingRecHit::ConstRecHitPointer &hitPtr, unsigned int iRow, const TrajectoryStateOnSurface &updatedTsos)
Definition: ReferenceTrajectory.cc:421
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
ReferenceTrajectoryBase::recHits
const TransientTrackingRecHit::ConstRecHitContainer & recHits() const
Definition: ReferenceTrajectoryBase.h:209
ReferenceTrajectory::addMaterialEffectsLocalGbl
virtual bool addMaterialEffectsLocalGbl(const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvatureChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParameterCovs)
Definition: ReferenceTrajectory.cc:939
ReferenceTrajectory::surfaceSide
SurfaceSide surfaceSide(const PropagationDirection dir) const
Definition: ReferenceTrajectory.h:167
ReferenceTrajectory::propagate
virtual bool propagate(const Plane &previousSurface, const TrajectoryStateOnSurface &previousTsos, const Plane &newSurface, TrajectoryStateOnSurface &newTsos, AlgebraicMatrix &newJacobian, AlgebraicMatrix &newCurvlinJacobian, double &nextStep, const MagneticField *magField) const
Definition: ReferenceTrajectory.cc:372
ReferenceTrajectory::addMaterialEffectsBp
virtual bool addMaterialEffectsBp(const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv)
Definition: ReferenceTrajectory.cc:556
GeomDet::surface
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
LocalTrajectoryParameters::signedInverseMomentum
float signedInverseMomentum() const
Signed inverse momentum q/p (zero for neutrals).
Definition: LocalTrajectoryParameters.h:116
ReferenceTrajectory::propDir_
const PropagationDirection propDir_
Definition: ReferenceTrajectory.h:194
ReferenceTrajectoryBase::theTrajectoryPositionCov
AlgebraicSymMatrix theTrajectoryPositionCov
Definition: ReferenceTrajectoryBase.h:243
reco::BeamSpot::dydz
double dydz() const
dydz slope
Definition: BeamSpot.h:84
ReferenceTrajectory::addMaterialEffectsCurvlinGbl
virtual bool addMaterialEffectsCurvlinGbl(const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs, const std::vector< AlgebraicMatrix > &allLocalToCurv)
Definition: ReferenceTrajectory.cc:1023
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:122
ReferenceTrajectory::SurfaceSide
SurfaceSideDefinition::SurfaceSide SurfaceSide
Definition: ReferenceTrajectory.h:58
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
TrackingRecHit::ConstRecHitPointer
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
Definition: TrackingRecHit.h:28
Plane::build
static PlanePointer build(Args &&...args)
Definition: Plane.h:33
TrajectoryStateOnSurface::freeState
FreeTrajectoryState const * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:70
ReferenceTrajectory::createUpdator
MaterialEffectsUpdator * createUpdator(MaterialEffects materialEffects, double mass) const
Definition: ReferenceTrajectory.cc:346
ReferenceTrajectory::fillDerivatives
virtual void fillDerivatives(const AlgebraicMatrix &projection, const AlgebraicMatrix &fullJacobian, unsigned int iRow)
Definition: ReferenceTrajectory.cc:460
AnalyticalPropagator_cfi.AnalyticalPropagator
AnalyticalPropagator
Definition: AnalyticalPropagator_cfi.py:3
ReferenceTrajectory::materialEffects_
const MaterialEffects materialEffects_
Definition: ReferenceTrajectory.h:193
FreeTrajectoryState::momentum
GlobalVector momentum() const
Definition: FreeTrajectoryState.h:87
reco::BeamSpot::dxdz
double dxdz() const
dxdz slope
Definition: BeamSpot.h:82
ReferenceTrajectoryBase::theRecHits
TransientTrackingRecHit::ConstRecHitContainer theRecHits
Definition: ReferenceTrajectoryBase.h:237
TrackingRecHit::ConstRecHitContainer
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: TrackingRecHit.h:35
AlgebraicVector
CLHEP::HepVector AlgebraicVector
Definition: AlgebraicObjects.h:12
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
ReferenceTrajectoryBase::theTsosVec
std::vector< TrajectoryStateOnSurface > theTsosVec
Definition: ReferenceTrajectoryBase.h:236
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:66
ReferenceTrajectory::addMaterialEffectsCov
virtual bool addMaterialEffectsCov(const std::vector< AlgebraicMatrix > &allJacobians, const std::vector< AlgebraicMatrix > &allProjections, const std::vector< AlgebraicSymMatrix > &allCurvChanges, const std::vector< AlgebraicSymMatrix > &allDeltaParaCovs)
Definition: ReferenceTrajectory.cc:488
ReferenceTrajectory::fillTrajectoryPositions
virtual void fillTrajectoryPositions(const AlgebraicMatrix &projection, const AlgebraicVector &mixedLocalParams, unsigned int iRow)
Definition: ReferenceTrajectory.cc:475
PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi.propagator
propagator
Definition: PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi.py:41
ReferenceTrajectoryBase::useRecHit
bool useRecHit(const TransientTrackingRecHit::ConstRecHitPointer &hitPtr) const
Definition: ReferenceTrajectoryBase.cc:35
TSCBLBuilderNoMaterial_cfi.TSCBLBuilderNoMaterial
TSCBLBuilderNoMaterial
Definition: TSCBLBuilderNoMaterial_cfi.py:3
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:64
MaterialEffectsUpdator * ReferenceTrajectory::createUpdator ( MaterialEffects  materialEffects,
double  mass 
) const
protected

internal method to get apropriate updator

Definition at line 346 of file ReferenceTrajectory.cc.

References ReferenceTrajectoryBase::breakPoints, ReferenceTrajectoryBase::brokenLinesCoarse, ReferenceTrajectoryBase::brokenLinesFine, ReferenceTrajectoryBase::combined, ReferenceTrajectoryBase::curvlinGBL, ReferenceTrajectoryBase::energyLoss, GsfElectronMaterialEffects_cfi::EnergyLossUpdator, ReferenceTrajectoryBase::localGBL, ReferenceTrajectoryBase::multipleScattering, GsfElectronMaterialEffects_cfi::MultipleScatteringUpdator, and ReferenceTrajectoryBase::none.

Referenced by construct().

347 {
348  switch (materialEffects) {
349  // MultipleScatteringUpdator doesn't change the trajectory-state
350  // during update and can therefore be used if material effects should be ignored:
351  case none:
352  case multipleScattering:
353  return new MultipleScatteringUpdator(mass);
354  case energyLoss:
355  return new EnergyLossUpdator(mass);
356  case combined:
357  return new CombinedMaterialEffectsUpdator(mass);
358  case breakPoints:
359  return new CombinedMaterialEffectsUpdator(mass);
360  case brokenLinesCoarse:
361  case brokenLinesFine:
362  case localGBL:
363  case curvlinGBL:
364  return new CombinedMaterialEffectsUpdator(mass);
365 }
366 
367  return 0;
368 }
GsfElectronMaterialEffects_cfi.MultipleScatteringUpdator
MultipleScatteringUpdator
Definition: GsfElectronMaterialEffects_cfi.py:7
GsfElectronMaterialEffects_cfi.EnergyLossUpdator
EnergyLossUpdator
Definition: GsfElectronMaterialEffects_cfi.py:5
void ReferenceTrajectory::fillDerivatives ( const AlgebraicMatrix &  projection,
const AlgebraicMatrix &  fullJacobian,
unsigned int  iRow 
)
protectedvirtual

internal method to fill derivatives for hit iRow/2

Definition at line 460 of file ReferenceTrajectory.cc.

References mps_fire::i, ReferenceTrajectoryBase::parameters(), and ReferenceTrajectoryBase::theDerivatives.

Referenced by construct().

463 {
464  // derivatives of the local coordinates of the reference track w.r.t. to the inital track-parameters
465  const AlgebraicMatrix projectedJacobian(projection * fullJacobian);
466  for (int i = 0; i < parameters().num_row(); ++i) {
467  for (int j = 0; j < projectedJacobian.num_row(); ++j) {
468  theDerivatives[iRow+j][i] = projectedJacobian[j][i];
469  }
470  }
471 }
mps_fire.i
i
Definition: mps_fire.py:156
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
ReferenceTrajectoryBase::parameters
const AlgebraicVector & parameters() const
Definition: ReferenceTrajectoryBase.h:188
void ReferenceTrajectory::fillMeasurementAndError ( const TransientTrackingRecHit::ConstRecHitPointer &  hitPtr,
unsigned int  iRow,
const TrajectoryStateOnSurface &  updatedTsos 
)
protectedvirtual

internal method to fill measurement and error matrix for hit iRow/2

Definition at line 421 of file ReferenceTrajectory.cc.

References includeAPEs_, TrackerGeomDet::localAlignmentError(), ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), PV3DBase< T, PVType, FrameType >::y(), and LocalError::yy().

Referenced by construct().

424 {
425  // get the measurements and their errors, use information updated with tsos if improving
426  // (GF: Also for measurements or only for errors or do the former not change?)
427  // GF 10/2008: I doubt that it makes sense to update the hit with the tsos here:
428  // That is an analytical extrapolation and not the best guess of the real
429  // track state on the module, but the latter should be better to get the best
430  // hit uncertainty estimate!
431 
432  // FIXME FIXME CLONE
433  auto newHitPtr = hitPtr;
434 // TransientTrackingRecHit::ConstRecHitPointer newHitPtr(hitPtr->canImproveWithTrack() ?
435 // hitPtr->clone(updatedTsos) : hitPtr);
436 
437  const LocalPoint localMeasurement = newHitPtr->localPosition();
438  const LocalError localMeasurementCov = newHitPtr->localPositionError(); // CPE+APE
439 
440  theMeasurements[iRow] = localMeasurement.x();
441  theMeasurements[iRow+1] = localMeasurement.y();
442  theMeasurementsCov[iRow][iRow] = localMeasurementCov.xx();
443  theMeasurementsCov[iRow][iRow+1] = localMeasurementCov.xy();
444  theMeasurementsCov[iRow+1][iRow+1] = localMeasurementCov.yy();
445 
446  if (!includeAPEs_) {
447  // subtract APEs (if existing) from covariance matrix
448  auto det = static_cast<const TrackerGeomDet*>(newHitPtr->det());
449  const auto localAPE = det->localAlignmentError();
450  if (localAPE.valid()) {
451  theMeasurementsCov[iRow][iRow] -= localAPE.xx();
452  theMeasurementsCov[iRow][iRow+1] -= localAPE.xy();
453  theMeasurementsCov[iRow+1][iRow+1] -= localAPE.yy();
454  }
455  }
456 }
LocalError::xx
float xx() const
Definition: LocalError.h:24
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
LocalError::xy
float xy() const
Definition: LocalError.h:25
LocalError::yy
float yy() const
Definition: LocalError.h:26
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
TrackerGeomDet::localAlignmentError
LocalError const & localAlignmentError() const
Return local alligment error.
Definition: TrackerGeomDet.h:14
void ReferenceTrajectory::fillTrajectoryPositions ( const AlgebraicMatrix &  projection,
const AlgebraicVector &  mixedLocalParams,
unsigned int  iRow 
)
protectedvirtual

internal method to fill the trajectory positions for hit iRow/2

Definition at line 475 of file ReferenceTrajectory.cc.

References mps_fire::i, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

478 {
479  // get the local coordinates of the reference trajectory
480  const AlgebraicVector localPosition(projection * mixedLocalParams);
481  for (int i = 0; i < localPosition.num_row(); ++i) {
482  theTrajectoryPositions[iRow+i] = localPosition[i];
483  }
484 }
mps_fire.i
i
Definition: mps_fire.py:156
AlgebraicVector
CLHEP::HepVector AlgebraicVector
Definition: AlgebraicObjects.h:12
AlgebraicMatrix ReferenceTrajectory::getHitProjectionMatrix ( const TransientTrackingRecHit::ConstRecHitPointer &  recHit) const
protected

first (generic) helper to get the projection matrix

Definition at line 1157 of file ReferenceTrajectory.cc.

References Exception, getHitProjectionMatrixT(), and ReferenceTrajectoryBase::useRecHit().

Referenced by clhep2eigen(), and construct().

1158 {
1159  if (this->useRecHit(hitPtr)) {
1160  // check which templated non-member function to call:
1161  switch (hitPtr->dimension()) {
1162  case 1:
1163  return getHitProjectionMatrixT<1>(hitPtr);
1164  case 2:
1165  return getHitProjectionMatrixT<2>(hitPtr);
1166  case 3:
1167  return getHitProjectionMatrixT<3>(hitPtr);
1168  case 4:
1169  return getHitProjectionMatrixT<4>(hitPtr);
1170  case 5:
1171  return getHitProjectionMatrixT<5>(hitPtr);
1172  default:
1173  throw cms::Exception("ReferenceTrajectory::getHitProjectionMatrix")
1174  << "Unexpected hit dimension: " << hitPtr->dimension() << "\n";
1175  }
1176  }
1177  // invalid or (to please compiler) unknown dimension
1178  return AlgebraicMatrix(2, 5, 0); // get size from ???
1179 }
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
ReferenceTrajectoryBase::useRecHit
bool useRecHit(const TransientTrackingRecHit::ConstRecHitPointer &hitPtr) const
Definition: ReferenceTrajectoryBase.cc:35
template<unsigned int N>
AlgebraicMatrix ReferenceTrajectory::getHitProjectionMatrixT ( const TransientTrackingRecHit::ConstRecHitPointer &  recHit) const
protected

second helper (templated on the dimension) to get the projection matrix

Definition at line 1186 of file ReferenceTrajectory.cc.

References N, packedPFCandidateRefMixer_cfi::pf, KfComponentsHolder::projection(), alignCSCRings::r, and KfComponentsHolder::setup().

Referenced by getHitProjectionMatrix().

1187 {
1188  // define variables that will be used to setup the KfComponentsHolder
1189  // (their allocated memory is needed by 'hitPtr->getKfComponents(..)'
1190 
1191  ProjectMatrix<double,5,N> pf;
1192  typename AlgebraicROOTObject<N>::Vector r, rMeas;
1193  typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas;
1194  // input for the holder - but dummy is OK here to just get the projection matrix:
1195  const AlgebraicVector5 dummyPars;
1196  const AlgebraicSymMatrix55 dummyErr;
1197 
1198  // setup the holder with the correct dimensions and get the values
1199  KfComponentsHolder holder;
1200  holder.setup<N>(&r, &V, &pf, &rMeas, &VMeas, dummyPars, dummyErr);
1201  hitPtr->getKfComponents(holder);
1202 
1203  return asHepMatrix<N,5>(holder.projection<N>());
1204 }
KfComponentsHolder::setup
void setup(typename AlgebraicROOTObject< D >::Vector *params, typename AlgebraicROOTObject< D, D >::SymMatrix *errors, ProjectMatrix< double, 5, D > *projFunc, typename AlgebraicROOTObject< D >::Vector *measuredParams, typename AlgebraicROOTObject< D, D >::SymMatrix *measuredErrors, const AlgebraicVector5 &tsosLocalParameters, const AlgebraicSymMatrix55 &tsosLocalErrors)
Definition: KfComponentsHolder.h:108
AlgebraicROOTObject::SymMatrix
ROOT::Math::SMatrix< double, D1, D1, ROOT::Math::MatRepSym< double, D1 > > SymMatrix
Definition: AlgebraicROOTObjects.h:68
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
KfComponentsHolder::projection
AlgebraicROOTObject< D, 5 >::Matrix projection()
Definition: KfComponentsHolder.h:51
AlgebraicROOTObject::Vector
ROOT::Math::SVector< double, D1 > Vector
Definition: AlgebraicROOTObjects.h:67
N
#define N
Definition: blowfish.cc:9
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:14
bool ReferenceTrajectory::propagate ( const Plane &  previousSurface,
const TrajectoryStateOnSurface &  previousTsos,
const Plane &  newSurface,
TrajectoryStateOnSurface &  newTsos,
AlgebraicMatrix &  newJacobian,
AlgebraicMatrix &  newCurvlinJacobian,
double &  nextStep,
const MagneticField *  magField 
) const
protectedvirtual

internal method to calculate jacobian

From TrackingTools/ GeomPropagators/ interface/ AnalyticalPropagator.h NB: this propagator assumes constant, non-zero magnetic field parallel to the z-axis!

Definition at line 372 of file ReferenceTrajectory.cc.

References TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::localParameters(), Propagator::propagateWithPath(), defaultRKPropagator::Product::propagator, and propDir_.

Referenced by construct().

376 {
377  // propagate to next layer
381  //AnalyticalPropagator aPropagator(magField, propDir_);
382  // Hard coded RungeKutta instead Analytical (avoid bias in TEC), but
383  // work around TrackPropagation/RungeKutta/interface/RKTestPropagator.h and
384  // http://www.parashift.com/c++-faq-lite/strange-inheritance.html#faq-23.9
385  defaultRKPropagator::Product rkprod(magField, propDir_); //double tolerance = 5.e-5)
386  Propagator &aPropagator = rkprod.propagator;
387  const std::pair<TrajectoryStateOnSurface, double> tsosWithPath =
388  aPropagator.propagateWithPath(previousTsos, newSurface);
389 
390  // stop if propagation wasn't successful
391  if (!tsosWithPath.first.isValid()) return false;
392 
393  nextStep = tsosWithPath.second;
394  // calculate derivative of reference-track parameters on the actual layer w.r.t. the ones
395  // on the previous layer (both in global coordinates)
396  const AnalyticalCurvilinearJacobian aJacobian(previousTsos.globalParameters(),
397  tsosWithPath.first.globalPosition(),
398  tsosWithPath.first.globalMomentum(),
399  tsosWithPath.second);
400  const AlgebraicMatrix curvilinearJacobian = asHepMatrix<5,5>(aJacobian.jacobian());
401 
402  // jacobian of the track parameters on the previous layer for local->global transformation
403  const JacobianLocalToCurvilinear startTrafo(previousSurface, previousTsos.localParameters(), *magField);
404  const AlgebraicMatrix localToCurvilinear = asHepMatrix<5>(startTrafo.jacobian());
405 
406  // jacobian of the track parameters on the actual layer for global->local transformation
407  const JacobianCurvilinearToLocal endTrafo(newSurface, tsosWithPath.first.localParameters(), *magField);
408  const AlgebraicMatrix curvilinearToLocal = asHepMatrix<5>(endTrafo.jacobian());
409 
410  // compute derivative of reference-track parameters on the actual layer w.r.t. the ones on
411  // the previous layer (both in their local representation)
412  newCurvlinJacobian = curvilinearJacobian;
413  newJacobian = curvilinearToLocal * curvilinearJacobian * localToCurvilinear;
414  newTsos = tsosWithPath.first;
415 
416  return true;
417 }
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:107
ReferenceTrajectory::propDir_
const PropagationDirection propDir_
Definition: ReferenceTrajectory.h:194
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
Propagator::propagateWithPath
virtual std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &, const Surface &) const final
Definition: Propagator.cc:15
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:80
SurfaceSide ReferenceTrajectory::surfaceSide ( const PropagationDirection  dir) const
inlineprotected

Don't care for propagation direction 'anyDirection' - in that case the material effects are anyway not updated ...

Definition at line 167 of file ReferenceTrajectory.h.

References SurfaceSideDefinition::afterSurface, alongMomentum, SurfaceSideDefinition::beforeSurface, recoMuon::in, MillePedeFileConverter_cfg::out, and rpcPointValidation_cfi::recHit.

Referenced by BzeroReferenceTrajectory::BzeroReferenceTrajectory(), and construct().

168  {
169  return ( dir == alongMomentum ) ?
170  SurfaceSideDefinition::beforeSurface :
171  SurfaceSideDefinition::afterSurface;
172  }
dir
dbl *** dir
Definition: mlp_gen.cc:35

Member Data Documentation

const bool ReferenceTrajectory::allowZeroMaterial_
private

Definition at line 197 of file ReferenceTrajectory.h.

Referenced by addMaterialEffectsCurvlinGbl(), and addMaterialEffectsLocalGbl().

const bool ReferenceTrajectory::includeAPEs_
private

Definition at line 196 of file ReferenceTrajectory.h.

Referenced by fillMeasurementAndError().

const double ReferenceTrajectory::mass_
private

Definition at line 192 of file ReferenceTrajectory.h.

Referenced by construct().

const MaterialEffects ReferenceTrajectory::materialEffects_
private

Definition at line 193 of file ReferenceTrajectory.h.

Referenced by construct().

const PropagationDirection ReferenceTrajectory::propDir_
private

Definition at line 194 of file ReferenceTrajectory.h.

Referenced by construct(), and propagate().

const bool ReferenceTrajectory::useBeamSpot_
private

Definition at line 195 of file ReferenceTrajectory.h.

Referenced by construct().

Generated for CMSSW Reference Manual by doxygen 1.8.11