CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Types | Public Member Functions | 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 TMatrixD & gblExtDerivatives () const
 
const TVectorD & gblExtMeasurements () const
 
const TVectorD & gblExtPrecisions () const
 
std::vector< std::pair
< std::vector< gbl::GblPoint >
, TMatrixD > > & 		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

void clhep2root (const AlgebraicVector &in, TVectorD &out)
 
void clhep2root (const AlgebraicMatrix &in, TMatrixD &out)
 
void clhep2root (const AlgebraicSymMatrix &in, TMatrixDSym &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
 
TMatrixD theGblExtDerivatives
 
TVectorD theGblExtMeasurements
 
TVectorD theGblExtPrecisions
 
std::vector< std::pair
< std::vector< gbl::GblPoint >
, TMatrixD > > 		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 const unsigned int nMeasPerHit = 2
 

Detailed Description

Definition at line 56 of file ReferenceTrajectory.h.

Member Typedef Documentation

typedef SurfaceSideDefinition::SurfaceSide ReferenceTrajectory::SurfaceSide

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

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

Referenced by clone().

54  :
55  ReferenceTrajectoryBase(
56  (config.materialEffects >= brokenLinesCoarse) ? 1 : refTsos.localParameters().mixedFormatVector().kSize,
57  (config.useBeamSpot) ? recHits.size()+1 : recHits.size(),
58  (config.materialEffects >= brokenLinesCoarse) ?
59  2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size()) :
60  ( (config.materialEffects == breakPoints) ? 2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size())-2 : 0) ,
61  (config.materialEffects >= brokenLinesCoarse) ?
62  2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size())-4 :
63  ( (config.materialEffects == breakPoints) ? 2*((config.useBeamSpot) ? recHits.size()+1 : recHits.size())-2 : 0) ),
64  mass_(config.mass),
65  materialEffects_(config.materialEffects),
66  propDir_(config.propDir),
67  useBeamSpot_(config.useBeamSpot),
68  includeAPEs_(config.includeAPEs),
69  allowZeroMaterial_(config.allowZeroMaterial)
70 {
71  // no check against magField == 0
72  theParameters = asHepVector<5>( refTsos.localParameters().mixedFormatVector() );
73 
74  if (config.hitsAreReverse) {
75  TransientTrackingRecHit::ConstRecHitContainer fwdRecHits;
76  fwdRecHits.reserve(recHits.size());
77  for (TransientTrackingRecHit::ConstRecHitContainer::const_reverse_iterator it=recHits.rbegin();
78  it != recHits.rend(); ++it) {
79  fwdRecHits.push_back(*it);
80  }
81  theValidityFlag = this->construct(refTsos, fwdRecHits, magField, beamSpot);
82  } else {
83  theValidityFlag = this->construct(refTsos, recHits, magField, beamSpot);
84  }
85 }
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:108
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 76 of file ReferenceTrajectory.h.

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

Definition at line 90 of file ReferenceTrajectory.cc.

92  : ReferenceTrajectoryBase(
93  (config.materialEffects >= brokenLinesCoarse) ? 1 : nPar,
94  nHits,
95  (config.materialEffects >= brokenLinesCoarse) ? 2*nHits : ( (config.materialEffects == breakPoints) ? 2*nHits-2 : 0 ),
96  (config.materialEffects >= brokenLinesCoarse) ? 2*nHits-4 : ( (config.materialEffects == breakPoints) ? 2*nHits-2 : 0 ) ),
97  mass_(config.mass),
98  materialEffects_(config.materialEffects),
99  propDir_(config.propDir),
100  useBeamSpot_(config.useBeamSpot),
101  includeAPEs_(config.includeAPEs),
102  allowZeroMaterial_(config.allowZeroMaterial)
103 {}
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 553 of file ReferenceTrajectory.cc.

References relval_2017::k, cmsLHEtoEOSManager::l, ReferenceTrajectoryBase::nMeasPerHit, ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theNumberOfPars.

Referenced by construct().

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

References relval_2017::k, cmsLHEtoEOSManager::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().

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

References delta, plotBeamSpotDB::first, i, relval_2017::k, cmsLHEtoEOSManager::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().

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

References relval_2017::k, cmsLHEtoEOSManager::l, ReferenceTrajectoryBase::nMeasPerHit, and ReferenceTrajectoryBase::theMeasurementsCov.

Referenced by construct().

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

References gbl::GblPoint::addMeasurement(), gbl::GblPoint::addScatterer(), allowZeroMaterial_, clhep2root(), Exception, infinity, relval_2017::k, ReferenceTrajectoryBase::theGblInput, ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

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

References gbl::GblPoint::addMeasurement(), gbl::GblPoint::addScatterer(), allowZeroMaterial_, edm::Exception::categoryCode(), clhep2root(), alignCSCRings::e, cms::Exception::explainSelf(), edm::errors::FatalRootError, relval_2017::k, ReferenceTrajectoryBase::theGblInput, ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

940 {
941 //CHK: add material effects using general broken lines, no initial kinks
942 // local track parameters are defined in the TSO system
943 
944  const double minPrec = 1.0; // minimum precision to use measurement (reject measurements in strip direction)
945 
946  AlgebraicMatrix OffsetToLocal(5,2); // dLocal/dU
947  OffsetToLocal[3][0] = 1.;
948  OffsetToLocal[4][1] = 1.;
949  AlgebraicMatrix SlopeToLocal(5,2); // dLocal/dU'
950  SlopeToLocal[1][0] = 1.;
951  SlopeToLocal[2][1] = 1.;
952 
953  // GBL uses ROOT matrices as interface
954  TMatrixDSym covariance(2), measPrecision(2), scatPrecision(2);
955  TMatrixD jacPointToPoint(5,5), identity(5,5), proLocalToMeas(2,2);
956  identity.UnitMatrix();
957  TVectorD measurement(2), scatterer(2), measPrecDiag(2);
958  scatterer.Zero();
959  //bool initialKinks = (allCurvlinKinks.size()>0);
960 
961 // measurements and scatterers from hits
962  unsigned int numHits = allJacobians.size();
963  std::vector<GblPoint> GblPointList;
964  GblPointList.reserve(numHits);
965  for (unsigned int k = 0; k < numHits; ++k) {
966 
967  // GBL point to point jacobian
968  clhep2root(allJacobians[k] * allCurvatureChanges[k], jacPointToPoint);
969 
970  // GBL point
971  GblPoint aGblPoint( jacPointToPoint );
972 
973  // GBL projection from local to measurement system
974  clhep2root(allProjections[k] * OffsetToLocal, proLocalToMeas);
975 
976  // GBL measurement (residuum to initial trajectory)
977  clhep2root(theMeasurements.sub(2*k+1,2*k+2) - theTrajectoryPositions.sub(2*k+1,2*k+2), measurement);
978 
979  // GBL measurement covariance matrix
980  clhep2root(theMeasurementsCov.sub(2*k+1,2*k+2), covariance);
981 
982  // GBL add measurement to point
983  if (covariance(0,1) == 0.) {
984  // covariance matrix is diagonal, independent measurements
985  for (unsigned int row = 0; row < 2; ++row) {
986  measPrecDiag(row) = ( 0. < covariance(row,row) ? 1.0/covariance(row,row) : 0. );
987  }
988  aGblPoint.addMeasurement(proLocalToMeas, measurement, measPrecDiag, minPrec);
989  } else
990  {
991  // covariance matrix needs diagonalization
992  measPrecision = covariance; measPrecision.InvertFast();
993  aGblPoint.addMeasurement(proLocalToMeas, measurement, measPrecision, minPrec);
994  }
995 
996  // GBL multiple scattering (full matrix in local system)
997  clhep2root(allDeltaParameterCovs[k].similarityT(SlopeToLocal), scatPrecision);
998  try {
999  scatPrecision.InvertFast();
1000 
1001  // GBL add scatterer to point
1002  aGblPoint.addScatterer(scatterer, scatPrecision);
1003  } catch (const edm::Exception& e) {
1004  if (e.categoryCode() == edm::errors::FatalRootError &&
1005  e.explainSelf().find("matrix is singular") != std::string::npos &&
1006  allowZeroMaterial_) {
1007  // allow for modules without material and do nothing
1008  } else {
1009  throw;
1010  }
1011  }
1012 
1013  // add point to list
1014  GblPointList.push_back( aGblPoint );
1015  }
1016  // add list of points and transformation local to fit (=local) system at first hit
1017  theGblInput.push_back(std::make_pair(GblPointList, identity));
1018 
1019  return true;
1020 }
edm::Exception::categoryCode
Code categoryCode() const
Definition: EDMException.h:93
cms::Exception::explainSelf
virtual std::string explainSelf() const
Definition: Exception.cc:146
relval_2017.k
list k
Definition: relval_2017.py:23
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
ReferenceTrajectoryBase::theGblInput
std::vector< std::pair< std::vector< gbl::GblPoint >, TMatrixD > > theGblInput
Definition: ReferenceTrajectoryBase.h:255
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:240
ReferenceTrajectory::clhep2root
void clhep2root(const AlgebraicVector &in, TVectorD &out)
Definition: ReferenceTrajectory.cc:1126
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:46
void ReferenceTrajectory::clhep2root ( const AlgebraicVector &  in,
TVectorD &  out 
)
private

Definition at line 1126 of file ReferenceTrajectory.cc.

Referenced by addMaterialEffectsCurvlinGbl(), and addMaterialEffectsLocalGbl().

1126  {
1127  // convert from CLHEP to ROOT matrix
1128  for (int row = 0; row < in.num_row(); ++row) {
1129  out[row] = in[row];
1130  }
1131 }
fileinputsource_cfi.out
tuple out
Definition: fileinputsource_cfi.py:62
void ReferenceTrajectory::clhep2root ( const AlgebraicMatrix &  in,
TMatrixD &  out 
)
private

Definition at line 1133 of file ReferenceTrajectory.cc.

References cuy::col.

1133  {
1134  // convert from CLHEP to ROOT matrix
1135  for (int row = 0; row < in.num_row(); ++row) {
1136  for (int col = 0; col < in.num_col(); ++col) {
1137  out[row][col] = in[row][col];
1138  }
1139  }
1140 }
fileinputsource_cfi.out
tuple out
Definition: fileinputsource_cfi.py:62
cuy.col
int col
Definition: cuy.py:1008
void ReferenceTrajectory::clhep2root ( const AlgebraicSymMatrix &  in,
TMatrixDSym &  out 
)
private

Definition at line 1142 of file ReferenceTrajectory.cc.

References cuy::col.

1142  {
1143  // convert from CLHEP to ROOT matrix
1144  for (int row = 0; row < in.num_row(); ++row) {
1145  for (int col = 0; col < in.num_col(); ++col) {
1146  out[row][col] = in[row][col];
1147  }
1148  }
1149 }
fileinputsource_cfi.out
tuple out
Definition: fileinputsource_cfi.py:62
cuy.col
int col
Definition: cuy.py:1008
virtual ReferenceTrajectory* ReferenceTrajectory::clone ( void  ) const
inlinevirtual

Implements ReferenceTrajectoryBase.

Reimplemented in BzeroReferenceTrajectory.

Definition at line 78 of file ReferenceTrajectory.h.

References ReferenceTrajectory().

78 { 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:50
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 108 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(), HLT_FULL_cff::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().

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

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

Referenced by construct().

344 {
345  switch (materialEffects) {
346  // MultipleScatteringUpdator doesn't change the trajectory-state
347  // during update and can therefore be used if material effects should be ignored:
348  case none:
349  case multipleScattering:
350  return new MultipleScatteringUpdator(mass);
351  case energyLoss:
352  return new EnergyLossUpdator(mass);
353  case combined:
354  return new CombinedMaterialEffectsUpdator(mass);
355  case breakPoints:
356  return new CombinedMaterialEffectsUpdator(mass);
357  case brokenLinesCoarse:
358  case brokenLinesFine:
359  case localGBL:
360  case curvlinGBL:
361  return new CombinedMaterialEffectsUpdator(mass);
362 }
363 
364  return 0;
365 }
HLT_FULL_cff.MultipleScatteringUpdator
tuple MultipleScatteringUpdator
Definition: HLT_FULL_cff.py:4803
HLT_FULL_cff.EnergyLossUpdator
tuple EnergyLossUpdator
Definition: HLT_FULL_cff.py:4801
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 457 of file ReferenceTrajectory.cc.

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

Referenced by construct().

460 {
461  // derivatives of the local coordinates of the reference track w.r.t. to the inital track-parameters
462  const AlgebraicMatrix projectedJacobian(projection * fullJacobian);
463  for (int i = 0; i < parameters().num_row(); ++i) {
464  for (int j = 0; j < projectedJacobian.num_row(); ++j) {
465  theDerivatives[iRow+j][i] = projectedJacobian[j][i];
466  }
467  }
468 }
i
int i
Definition: DBlmapReader.cc:9
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
j
int j
Definition: DBlmapReader.cc:9
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 418 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().

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

References i, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

475 {
476  // get the local coordinates of the reference trajectory
477  const AlgebraicVector localPosition(projection * mixedLocalParams);
478  for (int i = 0; i < localPosition.num_row(); ++i) {
479  theTrajectoryPositions[iRow+i] = localPosition[i];
480  }
481 }
i
int i
Definition: DBlmapReader.cc:9
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 1155 of file ReferenceTrajectory.cc.

References Exception.

Referenced by construct().

1156 {
1157  if (this->useRecHit(hitPtr)) {
1158  // check which templated non-member function to call:
1159  switch (hitPtr->dimension()) {
1160  case 1:
1161  return getHitProjectionMatrixT<1>(hitPtr);
1162  case 2:
1163  return getHitProjectionMatrixT<2>(hitPtr);
1164  case 3:
1165  return getHitProjectionMatrixT<3>(hitPtr);
1166  case 4:
1167  return getHitProjectionMatrixT<4>(hitPtr);
1168  case 5:
1169  return getHitProjectionMatrixT<5>(hitPtr);
1170  default:
1171  throw cms::Exception("ReferenceTrajectory::getHitProjectionMatrix")
1172  << "Unexpected hit dimension: " << hitPtr->dimension() << "\n";
1173  }
1174  }
1175  // invalid or (to please compiler) unknown dimension
1176  return AlgebraicMatrix(2, 5, 0); // get size from ???
1177 }
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 1184 of file ReferenceTrajectory.cc.

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

1185 {
1186  // define variables that will be used to setup the KfComponentsHolder
1187  // (their allocated memory is needed by 'hitPtr->getKfComponents(..)'
1188 
1189  ProjectMatrix<double,5,N> pf;
1190  typename AlgebraicROOTObject<N>::Vector r, rMeas;
1191  typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas;
1192  // input for the holder - but dummy is OK here to just get the projection matrix:
1193  const AlgebraicVector5 dummyPars;
1194  const AlgebraicSymMatrix55 dummyErr;
1195 
1196  // setup the holder with the correct dimensions and get the values
1197  KfComponentsHolder holder;
1198  holder.setup<N>(&r, &V, &pf, &rMeas, &VMeas, dummyPars, dummyErr);
1199  hitPtr->getKfComponents(holder);
1200 
1201  return asHepMatrix<N,5>(holder.projection<N>());
1202 }
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
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 369 of file ReferenceTrajectory.cc.

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

Referenced by construct().

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

References SurfaceSideDefinition::afterSurface, alongMomentum, and SurfaceSideDefinition::beforeSurface.

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

171  {
172  return ( dir == alongMomentum ) ?
173  SurfaceSideDefinition::beforeSurface :
174  SurfaceSideDefinition::afterSurface;
175  }
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.5