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

ReferenceTrajectory * clone () const override
 
 ReferenceTrajectory (const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, const MagneticField *magField, const reco::BeamSpot &beamSpot, const ReferenceTrajectoryBase::Config &config)
 
 ~ReferenceTrajectory () override
 
- 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
 
 ~ReferenceTrajectoryBase () override
 

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 55 of file ReferenceTrajectory.h.

Member Typedef Documentation

typedef SurfaceSideDefinition::SurfaceSide ReferenceTrajectory::SurfaceSide

Definition at line 57 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 52 of file ReferenceTrajectory.cc.

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

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

Definition at line 72 of file ReferenceTrajectory.h.

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

Definition at line 95 of file ReferenceTrajectory.cc.

98  : ReferenceTrajectoryBase((config.materialEffects >= brokenLinesCoarse) ? 1 : nPar,
99  nHits,
100  (config.materialEffects >= brokenLinesCoarse)
101  ? 2 * nHits
102  : ((config.materialEffects == breakPoints) ? 2 * nHits - 2 : 0),
103  (config.materialEffects >= brokenLinesCoarse)
104  ? 2 * nHits - 4
105  : ((config.materialEffects == breakPoints) ? 2 * nHits - 2 : 0)),
106  mass_(config.mass),
107  materialEffects_(config.materialEffects),
108  propDir_(config.propDir),
109  useBeamSpot_(config.useBeamSpot),
110  includeAPEs_(config.includeAPEs),
111  allowZeroMaterial_(config.allowZeroMaterial) {}
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:192
ReferenceTrajectory::materialEffects_
const MaterialEffects materialEffects_
Definition: ReferenceTrajectory.h:191

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 558 of file ReferenceTrajectory.cc.

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

Referenced by construct().

562  {
563  //CHK: add material effects using break points
564  int offsetPar = theNumberOfPars;
565  int offsetMeas = nMeasPerHit * allJacobians.size();
566  int ierr = 0;
567 
568  AlgebraicMatrix tempJacobian;
569  AlgebraicMatrix MSprojection(2, 5);
570  MSprojection[0][1] = 1;
571  MSprojection[1][2] = 1;
572  AlgebraicSymMatrix tempMSCov;
573  AlgebraicSymMatrix tempMSCovProj;
574  AlgebraicMatrix tempMSJacProj;
575 
576  for (unsigned int k = 1; k < allJacobians.size(); ++k) {
577  // CHK
578  int kbp = k - 1;
579  tempJacobian = allJacobians[k] * allCurvatureChanges[k];
580  tempMSCov = allDeltaParameterCovs[k - 1].similarity(allLocalToCurv[k - 1]);
581  tempMSCovProj = tempMSCov.similarity(MSprojection);
582  theMeasurementsCov[offsetMeas + nMeasPerHit * kbp][offsetMeas + nMeasPerHit * kbp] = tempMSCovProj[0][0];
583  theMeasurementsCov[offsetMeas + nMeasPerHit * kbp + 1][offsetMeas + nMeasPerHit * kbp + 1] = tempMSCovProj[1][1];
584  theDerivatives[offsetMeas + nMeasPerHit * kbp][offsetPar + 2 * kbp] = 1.0;
585  theDerivatives[offsetMeas + nMeasPerHit * kbp + 1][offsetPar + 2 * kbp + 1] = 1.0;
586 
587  tempMSJacProj = (allProjections[k] * (tempJacobian * allLocalToCurv[k - 1].inverse(ierr))) * MSprojection.T();
588  if (ierr) {
589  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBp"
590  << "Inversion 1 for break points failed: " << ierr;
591  return false;
592  }
593  theDerivatives[nMeasPerHit * k][offsetPar + 2 * kbp] = tempMSJacProj[0][0];
594  theDerivatives[nMeasPerHit * k][offsetPar + 2 * kbp + 1] = tempMSJacProj[0][1];
595  theDerivatives[nMeasPerHit * k + 1][offsetPar + 2 * kbp] = tempMSJacProj[1][0];
596  theDerivatives[nMeasPerHit * k + 1][offsetPar + 2 * kbp + 1] = tempMSJacProj[1][1];
597 
598  for (unsigned int l = k + 1; l < allJacobians.size(); ++l) {
599  // CHK
600  int kbp = k - 1;
601  tempJacobian = allJacobians[l] * allCurvatureChanges[l] * tempJacobian;
602  tempMSJacProj = (allProjections[l] * (tempJacobian * allLocalToCurv[k - 1].inverse(ierr))) * MSprojection.T();
603  if (ierr) {
604  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBp"
605  << "Inversion 2 for break points failed: " << ierr;
606  return false;
607  }
608  theDerivatives[nMeasPerHit * l][offsetPar + 2 * kbp] = tempMSJacProj[0][0];
609  theDerivatives[nMeasPerHit * l][offsetPar + 2 * kbp + 1] = tempMSJacProj[0][1];
610  theDerivatives[nMeasPerHit * l + 1][offsetPar + 2 * kbp] = tempMSJacProj[1][0];
611  theDerivatives[nMeasPerHit * l + 1][offsetPar + 2 * kbp + 1] = tempMSJacProj[1][1];
612  }
613  }
614 
615  return true;
616 }
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:247
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 620 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().

625  {
626  //CHK: add material effects using broken lines
627  //fine: use exact Jacobians, all detectors
628  //broken lines: pair of offsets (u1,u2) = (xt,yt) (in curvilinear frame (q/p,lambda,phi,xt,yt)) at each layer
629  // scattering angles (alpha1,alpha2) = (cosLambda*dPhi, dLambda) (cosLambda cancels in Chi2)
630  // DU' = (dU'/dU)*DU + (dU'/dAlpha)*DAlpha + (dU'/dQbyp)*DQbyp (propagation of U)
631  // = J*DU + S*DAlpha + d*DQbyp
632  // => DAlpha = S^-1 (DU' - J*DU - d*DQbyp)
633 
634  int offsetPar = theNumberOfPars;
635  int offsetMeas = nMeasPerHit * allCurvlinJacobians.size();
636  int ierr = 0;
637 
638  GlobalVector p = gtp.momentum();
639  double cosLambda = sqrt((p.x() * p.x() + p.y() * p.y()) / (p.x() * p.x() + p.y() * p.y() + p.z() * p.z()));
640 
641  // transformations Curvilinear <-> BrokenLines
642  AlgebraicMatrix QbypToCurv(5, 1); // dCurv/dQbyp
643  QbypToCurv[0][0] = 1.; // dQbyp/dQbyp
644  AlgebraicMatrix AngleToCurv(5, 2); // dCurv/dAlpha
645  AngleToCurv[1][1] = 1.; // dlambda/dalpha2
646  AngleToCurv[2][0] = 1. / cosLambda; // dphi/dalpha1
647  AlgebraicMatrix CurvToAngle(2, 5); // dAlpha/dCurv
648  CurvToAngle[1][1] = 1.; // dalpha2/dlambda
649  CurvToAngle[0][2] = cosLambda; // dalpha1/dphi
650  AlgebraicMatrix OffsetToCurv(5, 2); // dCurv/dU
651  OffsetToCurv[3][0] = 1.; // dxt/du1
652  OffsetToCurv[4][1] = 1.; // dyt/du2
653  AlgebraicMatrix CurvToOffset(2, 5); // dU/dCurv
654  CurvToOffset[0][3] = 1.; // du1/dxt
655  CurvToOffset[1][4] = 1.; // du2/dyt
656 
657  // transformations trajectory to components (Qbyp, U1, U2)
658  AlgebraicMatrix TrajToQbyp(1, 5);
659  TrajToQbyp[0][0] = 1.;
660  AlgebraicMatrix TrajToOff1(2, 5);
661  TrajToOff1[0][1] = 1.;
662  TrajToOff1[1][2] = 1.;
663  AlgebraicMatrix TrajToOff2(2, 5);
664  TrajToOff2[0][3] = 1.;
665  TrajToOff2[1][4] = 1.;
666 
667  AlgebraicMatrix JacOffsetToAngleC, JacQbypToAngleC;
668  AlgebraicMatrix JacCurvToOffsetL, JacOffsetToOffsetL, JacAngleToOffsetL, JacQbypToOffsetL, JacOffsetToAngleL;
669  AlgebraicMatrix JacCurvToOffsetN, JacOffsetToOffsetN, JacAngleToOffsetN, JacQbypToOffsetN, JacOffsetToAngleN;
670 
671  // transformation from trajectory to curvilinear parameters
672 
673  JacCurvToOffsetN = CurvToOffset * allCurvlinJacobians[1]; // (dU'/dCurv') * (dCurv'/dCurv) @ 2nd point
674  JacOffsetToOffsetN = JacCurvToOffsetN * OffsetToCurv; // J: (dU'/dU) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
675  JacAngleToOffsetN =
676  JacCurvToOffsetN * AngleToCurv; // S: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
677  JacQbypToOffsetN = JacCurvToOffsetN * QbypToCurv; // d: (dU'/dQbyp) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
678  JacOffsetToAngleN = JacAngleToOffsetN.inverse(ierr); // W
679  if (ierr) {
680  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
681  << "Inversion 1 for fine broken lines failed: " << ierr;
682  return false;
683  }
684  JacOffsetToAngleC = -(JacOffsetToAngleN * JacOffsetToOffsetN); // (dAlpha/dU)
685  JacQbypToAngleC = -(JacOffsetToAngleN * JacQbypToOffsetN); // (dAlpha/dQbyp)
686  // (dAlpha/dTraj) = (dAlpha/dQbyp) * (dQbyp/dTraj) + (dAlpha/dU1) * (dU1/dTraj) + (dAlpha/dU2) * (dU2/dTraj)
687  AlgebraicMatrix JacTrajToAngle =
688  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] =
731  tempMSCovProj[0][0];
732 
733  // transformation matices for offsets ( l <- k -> n )
734  tempJacL = allCurvlinJacobians[k] * tempJacobian;
735  JacCurvToOffsetL = CurvToOffset * tempJacL.inverse(ierr); // (dU'/dCurv') * (dCurv'/dCurv) @ last point
736 
737  if (ierr) {
738  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
739  << "Inversion 4 for fine broken lines failed: " << ierr;
740  return false;
741  }
742  JacOffsetToOffsetL =
743  JacCurvToOffsetL * OffsetToCurv; // J-: (dU'/dU) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
744  JacAngleToOffsetL =
745  JacCurvToOffsetL * AngleToCurv; // S-: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
746  JacQbypToOffsetL =
747  JacCurvToOffsetL * QbypToCurv; // d-: (dU'/dQbyp) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
748  JacOffsetToAngleL = -JacAngleToOffsetL.inverse(ierr); // W-
749  if (ierr) {
750  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
751  << "Inversion 5 for fine broken lines failed: " << ierr;
752  return false;
753  }
754  tempJacobian = tempJacobian * allCurvatureChanges[k];
755  tempJacN = allCurvlinJacobians[n] * tempJacobian;
756  JacCurvToOffsetN = CurvToOffset * tempJacN; // (dU'/dCurv') * (dCurv'/dCurv) @ next point
757  JacOffsetToOffsetN =
758  JacCurvToOffsetN * OffsetToCurv; // J+: (dU'/dU) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
759  JacAngleToOffsetN =
760  JacCurvToOffsetN * AngleToCurv; // S+: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
761  JacQbypToOffsetN =
762  JacCurvToOffsetN * QbypToCurv; // d+: (dU'/dQbyp) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
763  JacOffsetToAngleN = JacAngleToOffsetN.inverse(ierr); // W+
764  if (ierr) {
765  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
766  << "Inversion 6 for fine broken lines failed: " << ierr;
767  return false;
768  }
769  JacOffsetToAngleC = -(JacOffsetToAngleL * JacOffsetToOffsetL + JacOffsetToAngleN * JacOffsetToOffsetN);
770  JacQbypToAngleC = -(JacOffsetToAngleL * JacQbypToOffsetL + JacOffsetToAngleN * JacQbypToOffsetN);
771 
772  // bending
773  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][0] = JacQbypToAngleC[0][0];
774  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][0] = JacQbypToAngleC[1][0];
775  // last layer
776  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][offsetPar + 2 * l] = JacOffsetToAngleL[0][0];
777  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][offsetPar + 2 * l + 1] = JacOffsetToAngleL[0][1];
778  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][offsetPar + 2 * l] = JacOffsetToAngleL[1][0];
779  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][offsetPar + 2 * l + 1] = JacOffsetToAngleL[1][1];
780  // current layer
781  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][offsetPar + 2 * k] = JacOffsetToAngleC[0][0];
782  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][offsetPar + 2 * k + 1] = JacOffsetToAngleC[0][1];
783  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][offsetPar + 2 * k] = JacOffsetToAngleC[1][0];
784  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][offsetPar + 2 * k + 1] = JacOffsetToAngleC[1][1];
785 
786  // next layer
787  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][offsetPar + 2 * n] = JacOffsetToAngleN[0][0];
788  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas][offsetPar + 2 * n + 1] = JacOffsetToAngleN[0][1];
789  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][offsetPar + 2 * n] = JacOffsetToAngleN[1][0];
790  theDerivatives[offsetMeas + nMeasPerHit * iMsMeas + 1][offsetPar + 2 * n + 1] = JacOffsetToAngleN[1][1];
791  }
792 
793  return true;
794 }
ReferenceTrajectoryBase::theInnerTrajectoryToCurvilinear
AlgebraicMatrix theInnerTrajectoryToCurvilinear
Definition: ReferenceTrajectoryBase.h:258
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
ReferenceTrajectoryBase::theInnerLocalToTrajectory
AlgebraicMatrix theInnerLocalToTrajectory
Definition: ReferenceTrajectoryBase.h:260
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:247
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 798 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().

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

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

Referenced by construct().

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

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

Referenced by construct().

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

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

Referenced by construct().

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

Definition at line 1132 of file ReferenceTrajectory.cc.

References MillePedeFileConverter_cfg::out.

Referenced by addMaterialEffectsCurvlinGbl(), and addMaterialEffectsLocalGbl().

1132  {
1133  static_assert(Derived::ColsAtCompileTime == 1, "clhep2eigen: 'out' must be of vector type");
1134  for (int row = 0; row < in.num_row(); ++row) {
1135  out(row) = in[row];
1136  }
1137 }
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 1140 of file ReferenceTrajectory.cc.

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

1140  {
1141  for (int row = 0; row < in.num_row(); ++row) {
1142  for (int col = 0; col < in.num_col(); ++col) {
1143  out(row, col) = in[row][col];
1144  }
1145  }
1146 }
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
cuy.col
col
Definition: cuy.py:1010
template<typename Derived >
void ReferenceTrajectory::clhep2eigen ( const AlgebraicSymMatrix &  in,
Eigen::MatrixBase< Derived > &  out 
)
private

Definition at line 1149 of file ReferenceTrajectory.cc.

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

1149  {
1150  for (int row = 0; row < in.num_row(); ++row) {
1151  for (int col = 0; col < in.num_col(); ++col) {
1152  out(row, col) = in[row][col];
1153  }
1154  }
1155 }
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
cuy.col
col
Definition: cuy.py:1010
ReferenceTrajectory* ReferenceTrajectory::clone ( void  ) const
inlineoverridevirtual

Implements ReferenceTrajectoryBase.

Definition at line 74 of file ReferenceTrajectory.h.

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

74 { 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:52
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 115 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().

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

350  {
351  switch (materialEffects) {
352  // MultipleScatteringUpdator doesn't change the trajectory-state
353  // during update and can therefore be used if material effects should be ignored:
354  case none:
355  case multipleScattering:
356  return new MultipleScatteringUpdator(mass);
357  case energyLoss:
358  return new EnergyLossUpdator(mass);
359  case combined:
360  return new CombinedMaterialEffectsUpdator(mass);
361  case breakPoints:
362  return new CombinedMaterialEffectsUpdator(mass);
363  case brokenLinesCoarse:
364  case brokenLinesFine:
365  case localGBL:
366  case curvlinGBL:
367  return new CombinedMaterialEffectsUpdator(mass);
368  }
369 
370  return nullptr;
371 }
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 466 of file ReferenceTrajectory.cc.

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

Referenced by construct().

468  {
469  // derivatives of the local coordinates of the reference track w.r.t. to the inital track-parameters
470  const AlgebraicMatrix projectedJacobian(projection * fullJacobian);
471  for (int i = 0; i < parameters().num_row(); ++i) {
472  for (int j = 0; j < projectedJacobian.num_row(); ++j) {
473  theDerivatives[iRow + j][i] = projectedJacobian[j][i];
474  }
475  }
476 }
mps_fire.i
i
Definition: mps_fire.py:338
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
ReferenceTrajectoryBase::parameters
const AlgebraicVector & parameters() const
Definition: ReferenceTrajectoryBase.h:191
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 428 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().

430  {
431  // get the measurements and their errors, use information updated with tsos if improving
432  // (GF: Also for measurements or only for errors or do the former not change?)
433  // GF 10/2008: I doubt that it makes sense to update the hit with the tsos here:
434  // That is an analytical extrapolation and not the best guess of the real
435  // track state on the module, but the latter should be better to get the best
436  // hit uncertainty estimate!
437 
438  // FIXME FIXME CLONE
439  const auto &newHitPtr = hitPtr;
440  // TransientTrackingRecHit::ConstRecHitPointer newHitPtr(hitPtr->canImproveWithTrack() ?
441  // hitPtr->clone(updatedTsos) : hitPtr);
442 
443  const LocalPoint localMeasurement = newHitPtr->localPosition();
444  const LocalError localMeasurementCov = newHitPtr->localPositionError(); // CPE+APE
445 
446  theMeasurements[iRow] = localMeasurement.x();
447  theMeasurements[iRow + 1] = localMeasurement.y();
448  theMeasurementsCov[iRow][iRow] = localMeasurementCov.xx();
449  theMeasurementsCov[iRow][iRow + 1] = localMeasurementCov.xy();
450  theMeasurementsCov[iRow + 1][iRow + 1] = localMeasurementCov.yy();
451 
452  if (!includeAPEs_) {
453  // subtract APEs (if existing) from covariance matrix
454  auto det = static_cast<const TrackerGeomDet *>(newHitPtr->det());
455  const auto localAPE = det->localAlignmentError();
456  if (localAPE.valid()) {
457  theMeasurementsCov[iRow][iRow] -= localAPE.xx();
458  theMeasurementsCov[iRow][iRow + 1] -= localAPE.xy();
459  theMeasurementsCov[iRow + 1][iRow + 1] -= localAPE.yy();
460  }
461  }
462 }
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:247
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 480 of file ReferenceTrajectory.cc.

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

Referenced by construct().

482  {
483  // get the local coordinates of the reference trajectory
484  const AlgebraicVector localPosition(projection * mixedLocalParams);
485  for (int i = 0; i < localPosition.num_row(); ++i) {
486  theTrajectoryPositions[iRow + i] = localPosition[i];
487  }
488 }
mps_fire.i
i
Definition: mps_fire.py:338
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 1159 of file ReferenceTrajectory.cc.

References Exception, and ReferenceTrajectoryBase::useRecHit().

Referenced by construct().

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

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 375 of file ReferenceTrajectory.cc.

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

Referenced by construct().

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

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

Member Data Documentation

const bool ReferenceTrajectory::allowZeroMaterial_
private

Definition at line 195 of file ReferenceTrajectory.h.

Referenced by addMaterialEffectsCurvlinGbl(), and addMaterialEffectsLocalGbl().

const bool ReferenceTrajectory::includeAPEs_
private

Definition at line 194 of file ReferenceTrajectory.h.

Referenced by fillMeasurementAndError().

const double ReferenceTrajectory::mass_
private

Definition at line 190 of file ReferenceTrajectory.h.

Referenced by construct().

const MaterialEffects ReferenceTrajectory::materialEffects_
private

Definition at line 191 of file ReferenceTrajectory.h.

Referenced by construct().

const PropagationDirection ReferenceTrajectory::propDir_
private

Definition at line 192 of file ReferenceTrajectory.h.

Referenced by construct(), and propagate().

const bool ReferenceTrajectory::useBeamSpot_
private

Definition at line 193 of file ReferenceTrajectory.h.

Referenced by construct().

Generated for CMSSW Reference Manual by doxygen 1.8.11