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
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, bool hitsAreReverse, const MagneticField *magField, MaterialEffects materialEffects, PropagationDirection propDir, double mass, bool useBeamSpot, const reco::BeamSpot &beamSpot)
 
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< 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, double mass, MaterialEffects materialEffects, const PropagationDirection propDir, const MagneticField *magField, bool useBeamSpot, 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 PropagationDirection propDir, const MagneticField *magField) const
 
 ReferenceTrajectory (unsigned int nPar, unsigned int nHits, MaterialEffects materialEffects)
 
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)
 

Additional Inherited Members

- Protected Attributes inherited from ReferenceTrajectoryBase
AlgebraicMatrix theDerivatives
 
TMatrixD theGblExtDerivatives
 
TVectorD theGblExtMeasurements
 
TVectorD theGblExtPrecisions
 
std::vector< std::pair
< std::vector< 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,
bool  hitsAreReverse,
const MagneticField *  magField,
MaterialEffects  materialEffects,
PropagationDirection  propDir,
double  mass,
bool  useBeamSpot,
const reco::BeamSpot &  beamSpot 
)

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

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

Referenced by clone().

56  : ReferenceTrajectoryBase(
57  (materialEffects >= brokenLinesCoarse) ? 1 : refTsos.localParameters().mixedFormatVector().kSize,
58  (useBeamSpot == true) ? recHits.size()+1 : recHits.size(),
59  (materialEffects >= brokenLinesCoarse) ?
60  2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size()) :
61  ( (materialEffects == breakPoints) ? 2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())-2 : 0) ,
62  (materialEffects >= brokenLinesCoarse) ?
63  2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())-4 :
64  ( (materialEffects == breakPoints) ? 2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())-2 : 0) )
65 {
66  // no check against magField == 0
67  theParameters = asHepVector<5>( refTsos.localParameters().mixedFormatVector() );
68 
69  if (hitsAreReverse) {
70  TransientTrackingRecHit::ConstRecHitContainer fwdRecHits;
71  fwdRecHits.reserve(recHits.size());
72  for (TransientTrackingRecHit::ConstRecHitContainer::const_reverse_iterator it=recHits.rbegin();
73  it != recHits.rend(); ++it) {
74  fwdRecHits.push_back(*it);
75  }
76  theValidityFlag = this->construct(refTsos, fwdRecHits, mass, materialEffects,
77  propDir, magField,
78  useBeamSpot, beamSpot);
79  } else {
80  theValidityFlag = this->construct(refTsos, recHits, mass, materialEffects,
81  propDir, magField,
82  useBeamSpot, beamSpot);
83  }
84 }
ReferenceTrajectory::construct
virtual bool construct(const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, double mass, MaterialEffects materialEffects, const PropagationDirection propDir, const MagneticField *magField, bool useBeamSpot, const reco::BeamSpot &beamSpot)
Definition: ReferenceTrajectory.cc:101
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:187
TrackingRecHit::ConstRecHitContainer
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: TrackingRecHit.h:35
HLT_25ns14e33_v1_cff.useBeamSpot
tuple useBeamSpot
Definition: HLT_25ns14e33_v1_cff.py:7517
virtual ReferenceTrajectory::~ReferenceTrajectory ( )
inlinevirtual

Definition at line 79 of file ReferenceTrajectory.h.

79 {}
ReferenceTrajectory::ReferenceTrajectory ( unsigned int  nPar,
unsigned int  nHits,
MaterialEffects  materialEffects 
)
protected

Definition at line 89 of file ReferenceTrajectory.cc.

91  : ReferenceTrajectoryBase(
92  (materialEffects >= brokenLinesCoarse) ? 1 : nPar,
93  nHits,
94  (materialEffects >= brokenLinesCoarse) ? 2*nHits : ( (materialEffects == breakPoints) ? 2*nHits-2 : 0 ),
95  (materialEffects >= brokenLinesCoarse) ? 2*nHits-4 : ( (materialEffects == breakPoints) ? 2*nHits-2 : 0 ) )
96 {}
ReferenceTrajectoryBase::ReferenceTrajectoryBase
ReferenceTrajectoryBase(unsigned int nPar, unsigned int nHits, unsigned int nVirtualPar, unsigned int nVirtualMeas)
Definition: ReferenceTrajectoryBase.cc:3

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

References relval_steps::k, prof2calltree::l, ReferenceTrajectoryBase::nMeasPerHit, ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theNumberOfPars.

Referenced by construct().

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

References relval_steps::k, prof2calltree::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().

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

References delta, plotBeamSpotDB::first, i, relval_steps::k, prof2calltree::l, prof2calltree::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().

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

References relval_steps::k, prof2calltree::l, ReferenceTrajectoryBase::nMeasPerHit, and ReferenceTrajectoryBase::theMeasurementsCov.

Referenced by construct().

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

References gbl::GblPoint::addMeasurement(), gbl::GblPoint::addScatterer(), clhep2root(), relval_steps::k, ReferenceTrajectoryBase::theGblInput, ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

1017 {
1018 //CHK: add material effects using general broken lines
1019 // local track parameters are defined in the curvilinear system
1020 
1021  const double minPrec = 1.0; // minimum precision to use measurement (reject measurements in strip direction)
1022  int ierr = 0;
1023 
1024  AlgebraicMatrix OffsetToCurv(5,2); // dCurv/dU
1025  OffsetToCurv[3][0] = 1.; // dxt/du1
1026  OffsetToCurv[4][1] = 1.; // dyt/du2
1027 
1028  AlgebraicMatrix JacOffsetToMeas, tempMSCov;
1029 
1030  // GBL uses ROOT matrices as interface
1031  TMatrixDSym covariance(2), measPrecision(2);
1032  TMatrixD jacPointToPoint(5,5), firstLocalToCurv(5,5), proLocalToMeas(2,2);
1033  TVectorD measurement(2), scatterer(2), measPrecDiag(2), scatPrecDiag(2);
1034  scatterer.Zero();
1035 
1036 // measurements and scatterers from hits
1037  unsigned int numHits = allCurvlinJacobians.size();
1038  std::vector<GblPoint> GblPointList;
1039  GblPointList.reserve(numHits);
1040  for (unsigned int k = 0; k < numHits; ++k) {
1041 // (dMeas/dU) = (dMeas/dLoc) * (dLoc/dCurv) * (dCurv/dU)
1042  JacOffsetToMeas = (allProjections[k] * allLocalToCurv[k].inverse(ierr) ) * OffsetToCurv;
1043  if (ierr) {
1044  edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsGbl"
1045  << "Inversion 1 for general broken lines failed: " << ierr;
1046  return false;
1047  }
1048 
1049  // GBL point to point jacobian
1050  clhep2root(allCurvlinJacobians[k] * allCurvatureChanges[k], jacPointToPoint);
1051 
1052  // GBL point
1053  GblPoint aGblPoint( jacPointToPoint );
1054 
1055  // GBL projection from local to measurement system
1056  clhep2root(JacOffsetToMeas, proLocalToMeas);
1057 
1058  // GBL measurement (residuum to initial trajectory)
1059  clhep2root(theMeasurements.sub(2*k+1,2*k+2) - theTrajectoryPositions.sub(2*k+1,2*k+2), measurement);
1060 
1061  // GBL measurement covariance matrix
1062  clhep2root(theMeasurementsCov.sub(2*k+1,2*k+2), covariance);
1063 
1064  // GBL add measurement to point
1065  if (covariance(0,1) == 0.) {
1066  // covariance matrix is diagonal, independent measurements
1067  for (unsigned int row = 0; row < 2; ++row) {
1068  measPrecDiag(row) = ( 0. < covariance(row,row) ? 1.0/covariance(row,row) : 0. );
1069  }
1070  aGblPoint.addMeasurement(proLocalToMeas, measurement, measPrecDiag, minPrec);
1071  } else
1072  {
1073  // covariance matrix needs diagonalization
1074  measPrecision = covariance; measPrecision.InvertFast();
1075  aGblPoint.addMeasurement(proLocalToMeas, measurement, measPrecision, minPrec);
1076  }
1077 
1078  // GBL multiple scattering (diagonal matrix in curvilinear system)
1079  tempMSCov = allDeltaParameterCovs[k].similarity(allLocalToCurv[k]);
1080  for (unsigned int row = 0; row < 2; ++row) {
1081  scatPrecDiag(row) = 1.0/tempMSCov[row+1][row+1];
1082  }
1083 
1084  // GBL add scatterer to point
1085  aGblPoint.addScatterer(scatterer, scatPrecDiag);
1086 
1087  // add point to list
1088  GblPointList.push_back( aGblPoint );
1089  }
1090  // add list of points and transformation local to fit (=curvilinear) system at first hit
1091  clhep2root(allLocalToCurv[0], firstLocalToCurv);
1092  theGblInput.push_back(std::make_pair(GblPointList, firstLocalToCurv));
1093 
1094  return true;
1095 }
ReferenceTrajectoryBase::theGblInput
std::vector< std::pair< std::vector< GblPoint >, TMatrixD > > theGblInput
Definition: ReferenceTrajectoryBase.h:233
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
ReferenceTrajectoryBase::theMeasurementsCov
AlgebraicSymMatrix theMeasurementsCov
Definition: ReferenceTrajectoryBase.h:218
ReferenceTrajectory::clhep2root
void clhep2root(const AlgebraicVector &in, TVectorD &out)
Definition: ReferenceTrajectory.cc:1099
relval_steps.k
list k
Definition: relval_steps.py:1518
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 934 of file ReferenceTrajectory.cc.

References gbl::GblPoint::addMeasurement(), gbl::GblPoint::addScatterer(), clhep2root(), relval_steps::k, ReferenceTrajectoryBase::theGblInput, ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

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

Definition at line 1099 of file ReferenceTrajectory.cc.

Referenced by addMaterialEffectsCurvlinGbl(), and addMaterialEffectsLocalGbl().

1099  {
1100  // convert from CLHEP to ROOT matrix
1101  for (int row = 0; row < in.num_row(); ++row) {
1102  out[row] = in[row];
1103  }
1104 }
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
void ReferenceTrajectory::clhep2root ( const AlgebraicMatrix &  in,
TMatrixD &  out 
)
private

Definition at line 1106 of file ReferenceTrajectory.cc.

References cuy::col.

1106  {
1107  // convert from CLHEP to ROOT matrix
1108  for (int row = 0; row < in.num_row(); ++row) {
1109  for (int col = 0; col < in.num_col(); ++col) {
1110  out[row][col] = in[row][col];
1111  }
1112  }
1113 }
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
cuy.col
int col
Definition: cuy.py:1008
void ReferenceTrajectory::clhep2root ( const AlgebraicSymMatrix &  in,
TMatrixDSym &  out 
)
private

Definition at line 1115 of file ReferenceTrajectory.cc.

References cuy::col.

1115  {
1116  // convert from CLHEP to ROOT matrix
1117  for (int row = 0; row < in.num_row(); ++row) {
1118  for (int col = 0; col < in.num_col(); ++col) {
1119  out[row][col] = in[row][col];
1120  }
1121  }
1122 }
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
cuy.col
int col
Definition: cuy.py:1008
virtual ReferenceTrajectory* ReferenceTrajectory::clone ( void  ) const
inlinevirtual

Implements ReferenceTrajectoryBase.

Reimplemented in BzeroReferenceTrajectory.

Definition at line 81 of file ReferenceTrajectory.h.

References ReferenceTrajectory().

81 { return new ReferenceTrajectory(*this); }
ReferenceTrajectory::ReferenceTrajectory
ReferenceTrajectory(const TrajectoryStateOnSurface &referenceTsos, const TransientTrackingRecHit::ConstRecHitContainer &recHits, bool hitsAreReverse, const MagneticField *magField, MaterialEffects materialEffects, PropagationDirection propDir, double mass, bool useBeamSpot, const reco::BeamSpot &beamSpot)
Definition: ReferenceTrajectory.cc:48
bool ReferenceTrajectory::construct ( const TrajectoryStateOnSurface &  referenceTsos,
const TransientTrackingRecHit::ConstRecHitContainer &  recHits,
double  mass,
MaterialEffects  materialEffects,
const PropagationDirection  propDir,
const MagneticField *  magField,
bool  useBeamSpot,
const reco::BeamSpot &  beamSpot 
)
protectedvirtual

internal method to calculate members

Definition at line 101 of file ReferenceTrajectory.cc.

References addMaterialEffectsBp(), addMaterialEffectsBrl(), addMaterialEffectsCov(), addMaterialEffectsCurvlinGbl(), addMaterialEffectsLocalGbl(), alongMomentum, AnalyticalPropagator_cfi::AnalyticalPropagator, ReferenceTrajectoryBase::breakPoints, ReferenceTrajectoryBase::brokenLinesCoarse, ReferenceTrajectoryBase::brokenLinesFine, newFWLiteAna::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(), LocalTrajectoryError::matrix(), FreeTrajectoryState::momentum(), ReferenceTrajectoryBase::multipleScattering, ReferenceTrajectoryBase::nMeasPerHit, MagneticField::nominalValue(), ReferenceTrajectoryBase::none, PV3DBase< T, PVType, FrameType >::phi(), propagate(), HLT_25ns14e33_v1_cff::propagator, 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, ReferenceTrajectoryBase::useRecHit(), reco::BeamSpot::x0(), reco::BeamSpot::y0(), and reco::BeamSpot::z0().

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

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

internal method to get apropriate updator

Definition at line 339 of file ReferenceTrajectory.cc.

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

Referenced by construct().

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

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

Referenced by construct().

452 {
453  // derivatives of the local coordinates of the reference track w.r.t. to the inital track-parameters
454  const AlgebraicMatrix projectedJacobian(projection * fullJacobian);
455  for (int i = 0; i < parameters().num_row(); ++i) {
456  theDerivatives[iRow ][i] = projectedJacobian[0][i];
457  theDerivatives[iRow+1][i] = projectedJacobian[1][i];
458  // GF: Should be a loop once the hit dimension is not hardcoded as nMeasPerHit (to be checked):
459  // for (int j = 0; j < projection.num_col(); ++j) {
460  // theDerivatives[iRow+j][i] = projectedJacobian[j][i];
461  // }
462  }
463 }
i
int i
Definition: DBlmapReader.cc:9
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:13
ReferenceTrajectoryBase::parameters
const AlgebraicVector & parameters() const
Definition: ReferenceTrajectoryBase.h:166
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 414 of file ReferenceTrajectory.cc.

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

Referenced by construct().

417 {
418  // get the measurements and their errors, use information updated with tsos if improving
419  // (GF: Also for measurements or only for errors or do the former not change?)
420  // GF 10/2008: I doubt that it makes sense to update the hit with the tsos here:
421  // That is an analytical extrapolation and not the best guess of the real
422  // track state on the module, but the latter should be better to get the best
423  // hit uncertainty estimate!
424 
425  // FIXME FIXME CLONE
426  auto newHitPtr = hitPtr;
427 // TransientTrackingRecHit::ConstRecHitPointer newHitPtr(hitPtr->canImproveWithTrack() ?
428 // hitPtr->clone(updatedTsos) : hitPtr);
429 
430  const LocalPoint localMeasurement = newHitPtr->localPosition();
431  const LocalError localMeasurementCov = newHitPtr->localPositionError();
432 
433  theMeasurements[iRow] = localMeasurement.x();
434  theMeasurements[iRow+1] = localMeasurement.y();
435  theMeasurementsCov[iRow][iRow] = localMeasurementCov.xx();
436  theMeasurementsCov[iRow][iRow+1] = localMeasurementCov.xy();
437  theMeasurementsCov[iRow+1][iRow+1] = localMeasurementCov.yy();
438  // GF: Should be a loop once the hit dimension is not hardcoded as nMeasPerHit (to be checked):
439  // for (int i = 0; i < hitPtr->dimension(); ++i) {
440  // theMeasurements[iRow+i] = hitPtr->parameters()[i]; // fixme: parameters() is by value!
441  // for (int j = i; j < hitPtr->dimension(); ++j) {
442  // theMeasurementsCov[iRow+i][iRow+j] = hitPtr->parametersError()[i][j];
443  // }
444  // }
445 }
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:218
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
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 467 of file ReferenceTrajectory.cc.

References ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

470 {
471  // get the local coordinates of the reference trajectory
472  const AlgebraicVector localPosition(projection * mixedLocalParams);
473  theTrajectoryPositions[iRow] = localPosition[0];
474  theTrajectoryPositions[iRow+1] = localPosition[1];
475  // GF: Should be a loop once the hit dimension is not hardcoded as nMeasPerHit (to be checked):
476  // for (int j = 0; j < projection.num_col(); ++j) {
477  // theTrajectoryPositions[iRow+j] = localPosition[j];
478  // }
479 }
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 1128 of file ReferenceTrajectory.cc.

References edm::hlt::Exception.

Referenced by construct().

1129 {
1130  if (this->useRecHit(hitPtr)) {
1131  // check which templated non-member function to call:
1132  switch (hitPtr->dimension()) {
1133  case 1:
1134  return getHitProjectionMatrixT<1>(hitPtr);
1135  case 2:
1136  return getHitProjectionMatrixT<2>(hitPtr);
1137  case 3:
1138  return getHitProjectionMatrixT<3>(hitPtr);
1139  case 4:
1140  return getHitProjectionMatrixT<4>(hitPtr);
1141  case 5:
1142  return getHitProjectionMatrixT<5>(hitPtr);
1143  default:
1144  throw cms::Exception("ReferenceTrajectory::getHitProjectionMatrix")
1145  << "Unexpected hit dimension: " << hitPtr->dimension() << "\n";
1146  }
1147  }
1148  // invalid or (to please compiler) unknown dimension
1149  return AlgebraicMatrix(2, 5, 0); // get size from ???
1150 }
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 1157 of file ReferenceTrajectory.cc.

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

1158 {
1159  // define variables that will be used to setup the KfComponentsHolder
1160  // (their allocated memory is needed by 'hitPtr->getKfComponents(..)'
1161 
1162  ProjectMatrix<double,5,N> pf;
1163  typename AlgebraicROOTObject<N>::Vector r, rMeas;
1164  typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas;
1165  // input for the holder - but dummy is OK here to just get the projection matrix:
1166  const AlgebraicVector5 dummyPars;
1167  const AlgebraicSymMatrix55 dummyErr;
1168 
1169  // setup the holder with the correct dimensions and get the values
1170  KfComponentsHolder holder;
1171  holder.setup<N>(&r, &V, &pf, &rMeas, &VMeas, dummyPars, dummyErr);
1172  hitPtr->getKfComponents(holder);
1173 
1174  return asHepMatrix<N,5>(holder.projection<N>());
1175 }
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:76
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:31
KfComponentsHolder::projection
AlgebraicROOTObject< D, 5 >::Matrix projection()
Definition: KfComponentsHolder.h:51
AlgebraicROOTObject::Vector
ROOT::Math::SVector< double, D1 > Vector
Definition: AlgebraicROOTObjects.h:75
N
#define N
Definition: blowfish.cc:9
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:22
bool ReferenceTrajectory::propagate ( const Plane &  previousSurface,
const TrajectoryStateOnSurface &  previousTsos,
const Plane &  newSurface,
TrajectoryStateOnSurface &  newTsos,
AlgebraicMatrix &  newJacobian,
AlgebraicMatrix &  newCurvlinJacobian,
double &  nextStep,
const PropagationDirection  propDir,
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 365 of file ReferenceTrajectory.cc.

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

Referenced by construct().

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

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

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

174  {
175  return ( dir == alongMomentum ) ?
176  SurfaceSideDefinition::beforeSurface :
177  SurfaceSideDefinition::afterSurface;
178  }
dir
dbl *** dir
Definition: mlp_gen.cc:35
Generated for CMSSW Reference Manual by doxygen 1.8.5