#include <ReferenceTrajectory.h>

Inheritance diagram for ReferenceTrajectory:

Public Types
typedef SurfaceSideDefinition::SurfaceSide	SurfaceSide

Public Types inherited from ReferenceTrajectoryBase
enum	MaterialEffects { none, multipleScattering, energyLoss, combined, breakPoints, brokenLinesCoarse, brokenLinesFine }

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

bool	isValid ()

const AlgebraicMatrix &	localToTrajectory () const

const AlgebraicSymMatrix &	measurementErrors () const

const AlgebraicVector &	measurements () 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	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

Additional Inherited Members
Protected Attributes inherited from ReferenceTrajectoryBase
AlgebraicMatrix	theDerivatives

AlgebraicMatrix	theInnerLocalToTrajectory

AlgebraicMatrix	theInnerTrajectoryToCurvilinear

AlgebraicVector	theMeasurements

AlgebraicSymMatrix	theMeasurementsCov

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

Member Typedef Documentation

typedef SurfaceSideDefinition::SurfaceSide ReferenceTrajectory::SurfaceSide

Definition at line 56 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().

  : ReferenceTrajectoryBase( 
    (materialEffects >= brokenLinesCoarse) ? 1 : refTsos.localParameters().mixedFormatVector().kSize, 
    (useBeamSpot == true) ? recHits.size()+1 : recHits.size(),
    (materialEffects >= brokenLinesCoarse) ? 
        2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())   :
    ( (materialEffects == breakPoints) ? 2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())-2 : 0) , 
    (materialEffects >= brokenLinesCoarse) ? 
        2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())-4 : 
    ( (materialEffects == breakPoints) ? 2*((useBeamSpot == true) ? recHits.size()+1 : recHits.size())-2 : 0) )
 {
   // no check against magField == 0  
   theParameters = asHepVector<5>( refTsos.localParameters().mixedFormatVector() );
   
   if (hitsAreReverse) {
     TransientTrackingRecHit::ConstRecHitContainer fwdRecHits;
     fwdRecHits.reserve(recHits.size());
     for (TransientTrackingRecHit::ConstRecHitContainer::const_reverse_iterator it=recHits.rbegin();
      it != recHits.rend(); ++it) {
       fwdRecHits.push_back(*it);
     }
     theValidityFlag = this->construct(refTsos, fwdRecHits, mass, materialEffects,
                       propDir, magField,
                       useBeamSpot, beamSpot);
   } else {
     theValidityFlag = this->construct(refTsos, recHits, mass, materialEffects,
                       propDir, magField,
                       useBeamSpot, beamSpot);
   }
 }

virtual ReferenceTrajectory::~ReferenceTrajectory ( )

inlinevirtual

Definition at line 74 of file ReferenceTrajectory.h.

74 {}

ReferenceTrajectory::ReferenceTrajectory	(	unsigned int	nPar,
		unsigned int	nHits,
		MaterialEffects	materialEffects
	)

protected

Definition at line 89 of file ReferenceTrajectory.cc.

  : ReferenceTrajectoryBase( 
    (materialEffects >= brokenLinesCoarse) ? 1 : nPar, 
    nHits, 
    (materialEffects >= brokenLinesCoarse) ? 2*nHits   : ( (materialEffects == breakPoints) ? 2*nHits-2 : 0 ), 
    (materialEffects >= brokenLinesCoarse) ? 2*nHits-4 : ( (materialEffects == breakPoints) ? 2*nHits-2 : 0 ) )
 {}

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

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

Referenced by construct().

 {
 //CHK: add material effects using break points
   int offsetPar = theNumberOfPars; 
   int offsetMeas = nMeasPerHit * allJacobians.size();
   int ierr = 0; 
 
   AlgebraicMatrix tempJacobian;
   AlgebraicMatrix MSprojection(2,5);
   MSprojection[0][1] = 1;
   MSprojection[1][2] = 1;
   AlgebraicSymMatrix tempMSCov; 
   AlgebraicSymMatrix tempMSCovProj;
   AlgebraicMatrix tempMSJacProj;   
 
   for (unsigned int k = 1; k < allJacobians.size(); ++k) {
 // CHK 
     int kbp = k-1;
     tempJacobian = allJacobians[k] * allCurvatureChanges[k];
     tempMSCov = allDeltaParameterCovs[k-1].similarity(allLocalToCurv[k-1]);
     tempMSCovProj = tempMSCov.similarity(MSprojection);
     theMeasurementsCov[offsetMeas+nMeasPerHit*kbp  ][offsetMeas+nMeasPerHit*kbp  ] = tempMSCovProj[0][0];
     theMeasurementsCov[offsetMeas+nMeasPerHit*kbp+1][offsetMeas+nMeasPerHit*kbp+1]=  tempMSCovProj[1][1];
     theDerivatives[offsetMeas+nMeasPerHit*kbp  ][offsetPar+2*kbp  ] = 1.0;
     theDerivatives[offsetMeas+nMeasPerHit*kbp+1][offsetPar+2*kbp+1] = 1.0 ; 
 
     tempMSJacProj = (allProjections[k] * ( tempJacobian * allLocalToCurv[k-1].inverse(ierr) )) * MSprojection.T();
     if (ierr) {
       edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBp"
                              << "Inversion 1 for break points failed: " << ierr;
       return false;
     }
     theDerivatives[nMeasPerHit*k  ][offsetPar+2*kbp  ] =  tempMSJacProj[0][0];  
     theDerivatives[nMeasPerHit*k  ][offsetPar+2*kbp+1] =  tempMSJacProj[0][1]; 
     theDerivatives[nMeasPerHit*k+1][offsetPar+2*kbp  ] =  tempMSJacProj[1][0];  
     theDerivatives[nMeasPerHit*k+1][offsetPar+2*kbp+1] =  tempMSJacProj[1][1];
               
     for (unsigned int l = k+1; l < allJacobians.size(); ++l) {
 // CHK    
       int kbp = k-1;
       tempJacobian = allJacobians[l] * allCurvatureChanges[l] * tempJacobian; 
       tempMSJacProj = (allProjections[l] * ( tempJacobian * allLocalToCurv[k-1].inverse(ierr) )) * MSprojection.T();
       if (ierr) {
     edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBp"
                    << "Inversion 2 for break points failed: " << ierr;
         return false;
       }
       theDerivatives[nMeasPerHit*l  ][offsetPar+2*kbp  ] =  tempMSJacProj[0][0];  
       theDerivatives[nMeasPerHit*l  ][offsetPar+2*kbp+1] =  tempMSJacProj[0][1]; 
       theDerivatives[nMeasPerHit*l+1][offsetPar+2*kbp  ] =  tempMSJacProj[1][0];  
       theDerivatives[nMeasPerHit*l+1][offsetPar+2*kbp+1] =  tempMSJacProj[1][1]; 
 
     }
 
   }
 
   return true;
 }

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

References gen::k, ConfigFiles::l, GlobalTrajectoryParameters::momentum(), 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().

 {
 //CHK: add material effects using broken lines
 //fine: use exact Jacobians, all detectors  
 //broken lines: pair of offsets (u1,u2) = (xt,yt) (in curvilinear frame (q/p,lambda,phi,xt,yt)) at each layer
 //              scattering angles (alpha1,alpha2) = (cosLambda*dPhi, dLambda) (cosLambda cancels in Chi2)
 //              DU' = (dU'/dU)*DU + (dU'/dAlpha)*DAlpha + (dU'/dQbyp)*DQbyp (propagation of U)
 //                  = J*DU + S*DAlpha + d*DQbyp
 //           => DAlpha = S^-1 (DU' - J*DU - d*DQbyp)
 
   int offsetPar = theNumberOfPars;
   int offsetMeas = nMeasPerHit*allCurvlinJacobians.size();
   int ierr = 0;
 
   GlobalVector p = gtp.momentum();
   double cosLambda = sqrt((p.x()*p.x()+p.y()*p.y())/(p.x()*p.x()+p.y()*p.y()+p.z()*p.z()));
 
 // transformations Curvilinear <-> BrokenLines
   AlgebraicMatrix QbypToCurv(5,1);   // dCurv/dQbyp
   QbypToCurv[0][0] = 1.;             // dQbyp/dQbyp
   AlgebraicMatrix AngleToCurv(5,2);  // dCurv/dAlpha
   AngleToCurv[1][1] = 1.;            // dlambda/dalpha2
   AngleToCurv[2][0] = 1./cosLambda;  // dphi/dalpha1
   AlgebraicMatrix CurvToAngle(2,5);  // dAlpha/dCurv
   CurvToAngle[1][1] = 1.;            // dalpha2/dlambda
   CurvToAngle[0][2] = cosLambda;     // dalpha1/dphi
   AlgebraicMatrix OffsetToCurv(5,2); // dCurv/dU
   OffsetToCurv[3][0] = 1.;           // dxt/du1
   OffsetToCurv[4][1] = 1.;           // dyt/du2
   AlgebraicMatrix CurvToOffset(2,5); // dU/dCurv
   CurvToOffset[0][3] = 1.;           // du1/dxt
   CurvToOffset[1][4] = 1.;           // du2/dyt
 
 // transformations  trajectory to components (Qbyp, U1, U2)
   AlgebraicMatrix TrajToQbyp(1,5);
   TrajToQbyp[0][0] = 1.;
   AlgebraicMatrix TrajToOff1(2,5);
   TrajToOff1[0][1] = 1.;
   TrajToOff1[1][2] = 1.;
   AlgebraicMatrix TrajToOff2(2,5);
   TrajToOff2[0][3] = 1.;
   TrajToOff2[1][4] = 1.;
 
   AlgebraicMatrix JacOffsetToAngleC, JacQbypToAngleC;
   AlgebraicMatrix JacCurvToOffsetL, JacOffsetToOffsetL, JacAngleToOffsetL, JacQbypToOffsetL, JacOffsetToAngleL;
   AlgebraicMatrix JacCurvToOffsetN, JacOffsetToOffsetN, JacAngleToOffsetN, JacQbypToOffsetN, JacOffsetToAngleN;
 
 // transformation from trajectory to curvilinear parameters
 
   JacCurvToOffsetN = CurvToOffset * allCurvlinJacobians[1]; // (dU'/dCurv') * (dCurv'/dCurv) @ 2nd point
   JacOffsetToOffsetN = JacCurvToOffsetN * OffsetToCurv; // J: (dU'/dU)     = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
   JacAngleToOffsetN  = JacCurvToOffsetN * AngleToCurv;  // S: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
   JacQbypToOffsetN   = JacCurvToOffsetN * QbypToCurv;   // d: (dU'/dQbyp)  = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
   JacOffsetToAngleN  = JacAngleToOffsetN.inverse(ierr); // W
   if (ierr) {
      edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                                 << "Inversion 1 for fine broken lines failed: " << ierr;
      return false;
   }
   JacOffsetToAngleC = -(JacOffsetToAngleN * JacOffsetToOffsetN); // (dAlpha/dU)
   JacQbypToAngleC   = -(JacOffsetToAngleN * JacQbypToOffsetN);   // (dAlpha/dQbyp)
   // (dAlpha/dTraj) = (dAlpha/dQbyp) * (dQbyp/dTraj) + (dAlpha/dU1) * (dU1/dTraj) + (dAlpha/dU2) * (dU2/dTraj)
   AlgebraicMatrix JacTrajToAngle = JacQbypToAngleC * TrajToQbyp + JacOffsetToAngleC * TrajToOff1 + JacOffsetToAngleN * TrajToOff2;
   // (dCurv/dTraj) = (dCurv/dQbyp) * (dQbyp/dTraj) + (dCurv/dAlpha) * (dAlpha/dTraj) + (dCurv/dU) * (dU/dTraj)
   theInnerTrajectoryToCurvilinear = QbypToCurv * TrajToQbyp + AngleToCurv * JacTrajToAngle + OffsetToCurv * TrajToOff1;
   theInnerLocalToTrajectory = theInnerTrajectoryToCurvilinear.inverse(ierr) * allLocalToCurv[0];
   if (ierr) {
     edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                                << "Inversion 2 for fine broken lines failed: " << ierr;
     return false;
   }
 
   AlgebraicMatrix tempJacobian(allCurvatureChanges[0]);
   AlgebraicSymMatrix tempMSCov;
   AlgebraicSymMatrix tempMSCovProj;
   AlgebraicMatrix tempJacL, tempJacN;
   AlgebraicMatrix JacOffsetToMeas;
 
 // measurements from hits  
   for (unsigned int k = 0; k < allCurvlinJacobians.size(); ++k) {
 //  (dMeas/dU) = (dMeas/dLoc) * (dLoc/dCurv) * (dCurv/dU)
     JacOffsetToMeas = (allProjections[k] * allLocalToCurv[k].inverse(ierr) ) * OffsetToCurv;
     if (ierr) {
        edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                    << "Inversion 3 for fine broken lines failed: " << ierr;
        return false;
     }
     theDerivatives[nMeasPerHit*k  ][offsetPar+2*k  ] =  JacOffsetToMeas[0][0];
     theDerivatives[nMeasPerHit*k  ][offsetPar+2*k+1] =  JacOffsetToMeas[0][1];
     theDerivatives[nMeasPerHit*k+1][offsetPar+2*k  ] =  JacOffsetToMeas[1][0];
     theDerivatives[nMeasPerHit*k+1][offsetPar+2*k+1] =  JacOffsetToMeas[1][1];
   }
 
 // measurement of MS kink  
   for (unsigned int k = 1; k < allCurvlinJacobians.size()-1; ++k) {
 // CHK 
     int iMsMeas = k-1; 
     int l    = k-1; // last hit
     int n    = k+1; // next hit
 
 // amount of multiple scattering in layer k  (angular error perp to direction)  
     tempMSCov = allDeltaParameterCovs[k].similarity(allLocalToCurv[k]);
     tempMSCovProj = tempMSCov.similarity(CurvToAngle);
     theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas  ][offsetMeas+nMeasPerHit*iMsMeas  ] = tempMSCovProj[1][1];
     theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas+1][offsetMeas+nMeasPerHit*iMsMeas+1] = tempMSCovProj[0][0];
 
 // transformation matices for offsets ( l <- k -> n )
     tempJacL = allCurvlinJacobians[k] * tempJacobian;
     JacCurvToOffsetL = CurvToOffset * tempJacL.inverse(ierr); // (dU'/dCurv') * (dCurv'/dCurv) @ last point
 
     if (ierr) {
        edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                    << "Inversion 4 for fine broken lines failed: " << ierr;
        return false;
     }
     JacOffsetToOffsetL = JacCurvToOffsetL * OffsetToCurv; // J-: (dU'/dU)     = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
     JacAngleToOffsetL  = JacCurvToOffsetL * AngleToCurv;  // S-: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
     JacQbypToOffsetL   = JacCurvToOffsetL * QbypToCurv;   // d-: (dU'/dQbyp)  = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
     JacOffsetToAngleL  =-JacAngleToOffsetL.inverse(ierr); // W-
     if (ierr) {
        edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                    << "Inversion 5 for fine broken lines failed: " << ierr;
        return false;
     }
     tempJacobian = tempJacobian * allCurvatureChanges[k];
     tempJacN = allCurvlinJacobians[n] * tempJacobian;
     JacCurvToOffsetN = CurvToOffset * tempJacN; // (dU'/dCurv') * (dCurv'/dCurv) @ next point
     JacOffsetToOffsetN = JacCurvToOffsetN * OffsetToCurv; // J+: (dU'/dU)     = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dU)
     JacAngleToOffsetN  = JacCurvToOffsetN * AngleToCurv;  // S+: (dU'/dAlpha) = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dAlpha)
     JacQbypToOffsetN   = JacCurvToOffsetN * QbypToCurv;   // d+: (dU'/dQbyp)  = (dU'/dCurv') * (dCurv'/dCurv) * (dCurv/dQbyp)
     JacOffsetToAngleN  = JacAngleToOffsetN.inverse(ierr); // W+
     if (ierr) {
        edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                    << "Inversion 6 for fine broken lines failed: " << ierr;
        return false;
     }
     JacOffsetToAngleC = -(JacOffsetToAngleL * JacOffsetToOffsetL + JacOffsetToAngleN * JacOffsetToOffsetN);
     JacQbypToAngleC   = -(JacOffsetToAngleL * JacQbypToOffsetL   + JacOffsetToAngleN * JacQbypToOffsetN);
 
     // bending    
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][              0] = JacQbypToAngleC[0][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][              0] = JacQbypToAngleC[1][0];
     // last layer
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*l  ] = JacOffsetToAngleL[0][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*l+1] = JacOffsetToAngleL[0][1];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*l  ] = JacOffsetToAngleL[1][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*l+1] = JacOffsetToAngleL[1][1];
     // current layer
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*k  ] = JacOffsetToAngleC[0][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*k+1] = JacOffsetToAngleC[0][1];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*k  ] = JacOffsetToAngleC[1][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*k+1] = JacOffsetToAngleC[1][1];
 
     // next layer
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*n  ] = JacOffsetToAngleN[0][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*n+1] = JacOffsetToAngleN[0][1];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*n  ] = JacOffsetToAngleN[1][0];
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*n+1] = JacOffsetToAngleN[1][1];
 
   }
 
   return true;
 }

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

 {
 //CHK: add material effects using broken lines
 //BrokenLinesCoarse: combine close by detectors, 
 //                   use approximate Jacobians from Steps (limit Qbyp -> 0),
 //                   bending only in RPhi (B=(0,0,Bz)), no energy loss correction
 
   int offsetPar = theNumberOfPars; 
   int offsetMeas = nMeasPerHit*allSteps.size();
   int ierr = 0;
 
   GlobalVector p = gtp.momentum();
   double cosLambda = sqrt((p.x()*p.x()+p.y()*p.y())/(p.x()*p.x()+p.y()*p.y()+p.z()*p.z()));
   double bFac = -gtp.magneticFieldInInverseGeV(gtp.position()).mag();
 
   // transformation from trajectory to curvilinear parameters at refTsos
   double delta (1.0/allSteps[1]);    
   theInnerTrajectoryToCurvilinear[0][0] = 1;
   theInnerTrajectoryToCurvilinear[1][2] = -delta;  
   theInnerTrajectoryToCurvilinear[1][4] =  delta;    
   theInnerTrajectoryToCurvilinear[2][0] = -0.5*bFac/delta;
   theInnerTrajectoryToCurvilinear[2][1] = -delta/cosLambda;
   theInnerTrajectoryToCurvilinear[2][3] =  delta/cosLambda;
   theInnerTrajectoryToCurvilinear[3][1] = 1;
   theInnerTrajectoryToCurvilinear[4][2] = 1;
   theInnerLocalToTrajectory = theInnerTrajectoryToCurvilinear.inverse(ierr) * allLocalToCurv[0];
   if (ierr) {
     edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                                << "Inversion 1 for coarse broken lines failed: " << ierr;
     return false;
   }
 
   AlgebraicMatrix CurvToAngle(2,5);  // dAlpha/dCurv
   CurvToAngle[1][1] = 1.;            // dalpha2/dlambda
   CurvToAngle[0][2] = cosLambda;     // dalpha1/dphi
   AlgebraicMatrix OffsetToCurv(5,2); // dCurv/dU
   OffsetToCurv[3][0] = 1.;           // dxt/du1
   OffsetToCurv[4][1] = 1.;           // dyt/du2
 
   AlgebraicSymMatrix tempMSCov;
   AlgebraicSymMatrix tempMSCovProj;
   AlgebraicMatrix JacOffsetToMeas;
 
   // combine closeby detectors into single plane
   std::vector<unsigned int> first(allSteps.size());
   std::vector<unsigned int> last (allSteps.size());
   std::vector<unsigned int> plane(allSteps.size());
   std::vector<double> sPlane(allSteps.size());
   unsigned int nPlane = 0;
   double sTot = 0;
 
   for (unsigned int k = 1; k < allSteps.size(); ++k) {
     sTot += allSteps[k];
     if (fabs(allSteps[k])>minStep) { nPlane += 1; first[nPlane] = k; }
     last[nPlane] = k;
     plane[k] = nPlane;
     sPlane[nPlane] += sTot;
   }
   if (nPlane < 2) return false; // pathological cases: need at least 2 planes
 
   theNumberOfVirtualPars = 2*(nPlane+1);
   theNumberOfVirtualMeas = 2*(nPlane-1);// unsigned underflow for nPlane == 0...
   for (unsigned int k = 0; k <= nPlane; ++k) { sPlane[k] /= (double) (last[k]-first[k]+1); }
 
   // measurements from hits
   sTot = 0; 
   for (unsigned int k = 0; k < allSteps.size(); ++k) {
     sTot += allSteps[k];
 //  (dMeas/dU) = (dMeas/dLoc) * (dLoc/dCurv) * (dCurv/dU)
     JacOffsetToMeas = (allProjections[k] * allLocalToCurv[k].inverse(ierr) ) * OffsetToCurv;
     if (ierr) {
       edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::addMaterialEffectsBrl"
                                  << "Inversion 2 for coarse broken lines failed: " << ierr;
       return false;
     }
 
     unsigned int iPlane = plane[k];
     if (last[iPlane] == first[iPlane])
     { // single plane
       theDerivatives[nMeasPerHit*k  ][offsetPar+2*iPlane  ] =  JacOffsetToMeas[0][0];
       theDerivatives[nMeasPerHit*k  ][offsetPar+2*iPlane+1] =  JacOffsetToMeas[0][1];
       theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane  ] =  JacOffsetToMeas[1][0];
       theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane+1] =  JacOffsetToMeas[1][1];
     } else
     { // combined plane: (linear) interpolation
       unsigned int jPlane; // neighbor plane for interpolation
       if (fabs(sTot) < fabs(sPlane[iPlane])) { jPlane = (iPlane>0) ? iPlane - 1 : 1; }
       else  { jPlane = (iPlane<nPlane) ? iPlane + 1 : nPlane -1 ;}
       // interpolation weights
       double sDiff = sPlane[iPlane] - sPlane[jPlane];
       double iFrac = (sTot - sPlane[jPlane]) / sDiff;
       double jFrac = 1.0 - iFrac;
       theDerivatives[nMeasPerHit*k  ][offsetPar+2*iPlane  ] =  JacOffsetToMeas[0][0]*iFrac;
       theDerivatives[nMeasPerHit*k  ][offsetPar+2*iPlane+1] =  JacOffsetToMeas[0][1]*iFrac;
       theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane  ] =  JacOffsetToMeas[1][0]*iFrac;
       theDerivatives[nMeasPerHit*k+1][offsetPar+2*iPlane+1] =  JacOffsetToMeas[1][1]*iFrac;
       theDerivatives[nMeasPerHit*k  ][offsetPar+2*jPlane  ] =  JacOffsetToMeas[0][0]*jFrac;
       theDerivatives[nMeasPerHit*k  ][offsetPar+2*jPlane+1] =  JacOffsetToMeas[0][1]*jFrac;
       theDerivatives[nMeasPerHit*k+1][offsetPar+2*jPlane  ] =  JacOffsetToMeas[1][0]*jFrac;
       theDerivatives[nMeasPerHit*k+1][offsetPar+2*jPlane+1] =  JacOffsetToMeas[1][1]*jFrac;
       // 2nd order neglected
       // theDerivatives[nMeasPerHit*k  ][                   0] = -0.5*bFac*sDiff*iFrac*sDiff*jFrac*cosLambda;
     }
   }
 
 // measurement of MS kink
   for (unsigned int i = 1; i < nPlane; ++i) {
 // CHK 
     int iMsMeas = i-1;
     int l    = i-1; // last hit
     int n    = i+1; // next hit
 
 // amount of multiple scattering in plane k
     for (unsigned int k = first[i]; k <= last[i]; ++k) {
       tempMSCov = allDeltaParameterCovs[k].similarity(allLocalToCurv[k]);
       tempMSCovProj = tempMSCov.similarity(CurvToAngle);
       theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas  ][offsetMeas+nMeasPerHit*iMsMeas  ] += tempMSCovProj[0][0];
       theMeasurementsCov[offsetMeas+nMeasPerHit*iMsMeas+1][offsetMeas+nMeasPerHit*iMsMeas+1] += tempMSCovProj[1][1];
     }
 // broken line measurements for layer k, correlations between both planes neglected
     double stepK = sPlane[i] - sPlane[l];
     double stepN = sPlane[n] - sPlane[i];
     double deltaK (1.0/stepK);
     double deltaN (1.0/stepN);
     // bending (only in RPhi)
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][              0] = -0.5*bFac*(stepK+stepN)*cosLambda;
     // last layer
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*l  ] = deltaK;
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*l+1] = deltaK;
     // current layer
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*i  ] = -(deltaK + deltaN);
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*i+1] = -(deltaK + deltaN);
     // next layer
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas  ][offsetPar+2*n  ] = deltaN;
     theDerivatives[offsetMeas+nMeasPerHit*iMsMeas+1][offsetPar+2*n+1] = deltaN;
   }
 
   return true;
 }

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

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

Referenced by construct().

 {
   // the uncertainty due to multiple scattering is 'transferred' to the error matrix of the hits
 
   // GF: Needs update once hit dimension is not hardcoded as nMeasPerHit!
 
   AlgebraicSymMatrix materialEffectsCov(nMeasPerHit * allJacobians.size(), 0);
 
   // additional covariance-matrix of the measurements due to material-effects (single measurement)
   AlgebraicSymMatrix deltaMaterialEffectsCov;
 
   // additional covariance-matrix of the parameters due to material-effects
   AlgebraicSymMatrix paramMaterialEffectsCov(allDeltaParameterCovs[0]); //initialization
   // error-propagation to state after energy loss
   //GFback  paramMaterialEffectsCov = paramMaterialEffectsCov.similarity(allCurvatureChanges[0]);
 
   AlgebraicMatrix tempParameterCov;
   AlgebraicMatrix tempMeasurementCov;
 
   for (unsigned int k = 1; k < allJacobians.size(); ++k) {
     // error-propagation to next layer
     paramMaterialEffectsCov = paramMaterialEffectsCov.similarity(allJacobians[k]);
 
     // get dependences for the measurements
     deltaMaterialEffectsCov = paramMaterialEffectsCov.similarity(allProjections[k]);
     materialEffectsCov[nMeasPerHit*k  ][nMeasPerHit*k  ] = deltaMaterialEffectsCov[0][0];
     materialEffectsCov[nMeasPerHit*k  ][nMeasPerHit*k+1] = deltaMaterialEffectsCov[0][1];
     materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*k  ] = deltaMaterialEffectsCov[1][0];
     materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*k+1] = deltaMaterialEffectsCov[1][1];
 
     // GFback add uncertainties for the following layers due to scattering at this layer
     paramMaterialEffectsCov += allDeltaParameterCovs[k];
     // end GFback
     tempParameterCov = paramMaterialEffectsCov;
 
     // compute "inter-layer-dependencies"
     for (unsigned int l = k+1; l < allJacobians.size(); ++l) {
       tempParameterCov   = allJacobians[l]   * allCurvatureChanges[l] * tempParameterCov;
       tempMeasurementCov = allProjections[l] * tempParameterCov       * allProjections[k].T();
 
       materialEffectsCov[nMeasPerHit*l][nMeasPerHit*k] = tempMeasurementCov[0][0];
       materialEffectsCov[nMeasPerHit*k][nMeasPerHit*l] = tempMeasurementCov[0][0];
 
       materialEffectsCov[nMeasPerHit*l][nMeasPerHit*k+1] = tempMeasurementCov[0][1];
       materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*l] = tempMeasurementCov[0][1];
 
       materialEffectsCov[nMeasPerHit*l+1][nMeasPerHit*k] = tempMeasurementCov[1][0];
       materialEffectsCov[nMeasPerHit*k][nMeasPerHit*l+1] = tempMeasurementCov[1][0];
 
       materialEffectsCov[nMeasPerHit*l+1][nMeasPerHit*k+1] = tempMeasurementCov[1][1];
       materialEffectsCov[nMeasPerHit*k+1][nMeasPerHit*l+1] = tempMeasurementCov[1][1];
     } 
     // add uncertainties for the following layers due to scattering at this layer
     // GFback paramMaterialEffectsCov += allDeltaParameterCovs[k];    
     // error-propagation to state after energy loss
     paramMaterialEffectsCov = paramMaterialEffectsCov.similarity(allCurvatureChanges[k]);
      
   }
   theMeasurementsCov += materialEffectsCov;
 
   return true; // cannot fail
 }

virtual ReferenceTrajectory* ReferenceTrajectory::clone ( void ) const

inlinevirtual

Implements ReferenceTrajectoryBase.

Reimplemented in BzeroReferenceTrajectory.

Definition at line 76 of file ReferenceTrajectory.h.

References ReferenceTrajectory().

76 { 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.

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

 {   
   TrajectoryStateOnSurface theRefTsos = refTsos;
 
   const SurfaceSide surfaceSide = this->surfaceSide(propDir);
   // auto_ptr to avoid memory leaks in case of not reaching delete at end of method:
   std::auto_ptr<MaterialEffectsUpdator> aMaterialEffectsUpdator
     (this->createUpdator(materialEffects, mass));
   if (!aMaterialEffectsUpdator.get()) return false; // empty auto_ptr
 
   AlgebraicMatrix                 fullJacobian(theParameters.num_row(), theParameters.num_row());
   std::vector<AlgebraicMatrix>    allJacobians; 
   allJacobians.reserve(theNumberOfHits);
 
   TransientTrackingRecHit::ConstRecHitPointer  previousHitPtr;
   TrajectoryStateOnSurface                     previousTsos;
   AlgebraicSymMatrix              previousChangeInCurvature(theParameters.num_row(), 1);
   std::vector<AlgebraicSymMatrix> allCurvatureChanges; 
   allCurvatureChanges.reserve(theNumberOfHits);
   
   const LocalTrajectoryError zeroErrors(0., 0., 0., 0., 0.);
 
   std::vector<AlgebraicMatrix> allProjections;
   allProjections.reserve(theNumberOfHits);
   std::vector<AlgebraicSymMatrix> allDeltaParameterCovs;
   allDeltaParameterCovs.reserve(theNumberOfHits);
 
   // CHK
   std::vector<AlgebraicMatrix> allLocalToCurv;
   allLocalToCurv.reserve(theNumberOfHits); 
   std::vector<double> allSteps;
   allSteps.reserve(theNumberOfHits); 
   std::vector<AlgebraicMatrix>    allCurvlinJacobians; 
   allCurvlinJacobians.reserve(theNumberOfHits);
   
   AlgebraicMatrix firstCurvlinJacobian(5, 5, 1);
   
   unsigned int iRow = 0;
 
   // local storage vector of all rechits (including rechit for beam spot in case it is used)
   TransientTrackingRecHit::ConstRecHitContainer allRecHits;
 
   if (useBeamSpot && propDir==alongMomentum) {
     
     GlobalPoint bs(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
     
     TrajectoryStateClosestToBeamLine tsctbl(TSCBLBuilderNoMaterial()(*(refTsos.freeState()), beamSpot));
     if (!tsctbl.isValid()) {
       edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::construct" 
                  << "TrajectoryStateClostestToBeamLine invalid. Skip track.";
       return false;
     }
 
     FreeTrajectoryState pcaFts = tsctbl.trackStateAtPCA();
     GlobalVector bd(beamSpot.dxdz(), beamSpot.dydz(), 1.0);
     
     //propagation FIXME: Should use same propagator everywhere...
     AnalyticalPropagator propagator(magField);
     std::pair< TrajectoryStateOnSurface, double > tsosWithPath =
       propagator.propagateWithPath(pcaFts, refTsos.surface());
     
     if (!tsosWithPath.first.isValid()) return false;
     
     GlobalVector momDir(pcaFts.momentum());
     GlobalVector perpDir(bd.cross(momDir));
     Plane::RotationType rotation(perpDir, bd);
     
     BeamSpotGeomDet * bsGeom = new BeamSpotGeomDet(Plane::build(bs, rotation));
 
     // There is also a constructor taking the magentic field. Use this one instead?
     theRefTsos = TrajectoryStateOnSurface(pcaFts, bsGeom->surface());
     
     TransientTrackingRecHit::ConstRecHitPointer bsRecHit(new BeamSpotTransientTrackingRecHit(beamSpot,
                       bsGeom,
                       theRefTsos.freeState()->momentum().phi()));
     allRecHits.push_back(bsRecHit);
 
   }
   
   // copy all rechits to the local storage vector
   TransientTrackingRecHit::ConstRecHitContainer::const_iterator itRecHit;
   for ( itRecHit = recHits.begin(); itRecHit != recHits.end(); ++itRecHit ) { 
     const TransientTrackingRecHit::ConstRecHitPointer &hitPtr = *itRecHit;
     allRecHits.push_back(hitPtr);
   }
 
   for ( itRecHit = allRecHits.begin(); itRecHit != allRecHits.end(); ++itRecHit ) { 
     
     const TransientTrackingRecHit::ConstRecHitPointer &hitPtr = *itRecHit;
     theRecHits.push_back(hitPtr);
 
     if (0 == iRow) {
 
       // compute the derivatives of the reference-track's parameters w.r.t. the initial ones
       // derivative of the initial reference-track parameters w.r.t. themselves is of course the identity 
       fullJacobian = AlgebraicMatrix(theParameters.num_row(), theParameters.num_row(), 1);
       allJacobians.push_back(fullJacobian);
       theTsosVec.push_back(theRefTsos);
       const JacobianLocalToCurvilinear startTrafo(hitPtr->det()->surface(), theRefTsos.localParameters(), *magField);
       const AlgebraicMatrix localToCurvilinear =  asHepMatrix<5>(startTrafo.jacobian());
       if (materialEffects <= breakPoints) {
          theInnerTrajectoryToCurvilinear = asHepMatrix<5>(startTrafo.jacobian());
      theInnerLocalToTrajectory = AlgebraicMatrix(5, 5, 1);
       }  
       allLocalToCurv.push_back(localToCurvilinear);
       allSteps.push_back(0.);
       allCurvlinJacobians.push_back(firstCurvlinJacobian);
 
     } else {
 
       AlgebraicMatrix nextJacobian;
       AlgebraicMatrix nextCurvlinJacobian;
       double nextStep = 0.;
       TrajectoryStateOnSurface nextTsos;
 
       if (!this->propagate(previousHitPtr->det()->surface(), previousTsos,
                hitPtr->det()->surface(), nextTsos,
                nextJacobian, nextCurvlinJacobian, nextStep, propDir, magField)) {
     return false; // stop if problem...// no delete aMaterialEffectsUpdator needed
       }
       
       allJacobians.push_back(nextJacobian);
       fullJacobian = nextJacobian * previousChangeInCurvature * fullJacobian;
       theTsosVec.push_back(nextTsos);
       
       const JacobianLocalToCurvilinear startTrafo(hitPtr->det()->surface(), nextTsos.localParameters(), *magField);
       const AlgebraicMatrix localToCurvilinear =  asHepMatrix<5>(startTrafo.jacobian());
       allLocalToCurv.push_back(localToCurvilinear);
       if (nextStep == 0.) { 
     edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::construct"
                    << "step 0. from id " << previousHitPtr->geographicalId()
                    << " to " << hitPtr->det()->geographicalId() << ".";
     // brokenLinesFine will not work, brokenLinesCoarse combines close by layers
     if (materialEffects == brokenLinesFine) {
       edm::LogError("Alignment") << "@SUB=ReferenceTrajectory::construct" << "Skip track.";
       return false;
     }
       }
       allSteps.push_back(nextStep);
       allCurvlinJacobians.push_back(nextCurvlinJacobian);
 
     }
 
     // take material effects into account. since trajectory-state is constructed with errors equal zero,
     // the updated state contains only the uncertainties due to interactions in the current layer.
     const TrajectoryStateOnSurface tmpTsos(theTsosVec.back().localParameters(), zeroErrors,
                        theTsosVec.back().surface(), magField, surfaceSide);
     const TrajectoryStateOnSurface updatedTsos = aMaterialEffectsUpdator->updateState(tmpTsos, propDir);
 
     if ( !updatedTsos.isValid() ) return false;// no delete aMaterialEffectsUpdator needed
     
     if ( theTsosVec.back().localParameters().charge() )
     {
       previousChangeInCurvature[0][0] = updatedTsos.localParameters().signedInverseMomentum() 
     / theTsosVec.back().localParameters().signedInverseMomentum();
     }
     
     // get multiple-scattering covariance-matrix
     allDeltaParameterCovs.push_back( asHepMatrix<5>(updatedTsos.localError().matrix()) );
     allCurvatureChanges.push_back(previousChangeInCurvature);
     
     // projection-matrix tsos-parameters -> measurement-coordinates
     allProjections.push_back(this->getHitProjectionMatrix(hitPtr));
     // set start-parameters for next propagation. trajectory-state without error
     //  - no error propagation needed here.
     previousHitPtr = hitPtr;
     previousTsos   = TrajectoryStateOnSurface(updatedTsos.globalParameters(),
                                               updatedTsos.surface(), surfaceSide);
     
     if (materialEffects < brokenLinesCoarse) {
       this->fillDerivatives(allProjections.back(), fullJacobian, iRow);
     }
 
     AlgebraicVector mixedLocalParams = asHepVector<5>(theTsosVec.back().localParameters().mixedFormatVector());
     this->fillTrajectoryPositions(allProjections.back(), mixedLocalParams, iRow);
     if ( useRecHit( hitPtr ) ) this->fillMeasurementAndError(hitPtr, iRow, updatedTsos);
 
     iRow += nMeasPerHit;
   } // end of loop on hits
 
   bool msOK = true;
   switch (materialEffects) {
   case none:
     break;
   case multipleScattering:
   case energyLoss:
   case combined:
     msOK = this->addMaterialEffectsCov(allJacobians, allProjections, allCurvatureChanges,
                                        allDeltaParameterCovs);
     break;
   case breakPoints:
     msOK = this->addMaterialEffectsBp(allJacobians, allProjections, allCurvatureChanges,
                                       allDeltaParameterCovs, allLocalToCurv);
     break;
   case brokenLinesCoarse:
     msOK = this->addMaterialEffectsBrl(allProjections, allDeltaParameterCovs, allLocalToCurv,
                                        allSteps, refTsos.globalParameters());
     break;
   case brokenLinesFine:
     msOK = this->addMaterialEffectsBrl(allCurvlinJacobians, allProjections, allCurvatureChanges,
                                        allDeltaParameterCovs, allLocalToCurv, refTsos.globalParameters());
   }
   if (!msOK) return false;
  
   if (refTsos.hasError()) {
     AlgebraicSymMatrix parameterCov = asHepMatrix<5>(refTsos.localError().matrix());
     AlgebraicMatrix  parDeriv;
     if (theNumberOfVirtualPars>0) { 
       parDeriv = theDerivatives.sub( 1, nMeasPerHit*allJacobians.size(), 1, theParameters.num_row() ); 
     } else {
       parDeriv = theDerivatives; 
     }
     theTrajectoryPositionCov = parameterCov.similarity(parDeriv);
   } else {
     theTrajectoryPositionCov = AlgebraicSymMatrix(theDerivatives.num_row(), 1);
   }
 
   return true;
 }

MaterialEffectsUpdator * ReferenceTrajectory::createUpdator	(	MaterialEffects	materialEffects,
		double	mass
	)		const

protected

internal method to get apropriate updator

Definition at line 330 of file ReferenceTrajectory.cc.

References ReferenceTrajectoryBase::breakPoints, ReferenceTrajectoryBase::brokenLinesCoarse, ReferenceTrajectoryBase::brokenLinesFine, ReferenceTrajectoryBase::combined, ReferenceTrajectoryBase::energyLoss, ReferenceTrajectoryBase::multipleScattering, and ReferenceTrajectoryBase::none.

Referenced by construct().

 {
   switch (materialEffects) {
     // MultipleScatteringUpdator doesn't change the trajectory-state
     // during update and can therefore be used if material effects should be ignored:
   case none:
   case multipleScattering: 
     return new MultipleScatteringUpdator(mass);
   case energyLoss:
     return new EnergyLossUpdator(mass);
   case combined:
     return new CombinedMaterialEffectsUpdator(mass);
   case breakPoints:
     return new CombinedMaterialEffectsUpdator(mass);
   case brokenLinesCoarse:
   case brokenLinesFine:
     return new CombinedMaterialEffectsUpdator(mass);
 }
 
   return 0;
 }

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

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

Referenced by construct().

 {
   // derivatives of the local coordinates of the reference track w.r.t. to the inital track-parameters
   const AlgebraicMatrix projectedJacobian(projection * fullJacobian);
   for (int i = 0; i < parameters().num_row(); ++i) {
     theDerivatives[iRow  ][i] = projectedJacobian[0][i];
     theDerivatives[iRow+1][i] = projectedJacobian[1][i];
     // GF: Should be a loop once the hit dimension is not hardcoded as nMeasPerHit (to be checked):
     // for (int j = 0; j < projection.num_col(); ++j) {
     //   theDerivatives[iRow+j][i] = projectedJacobian[j][i];
     // }
   }
 }

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 403 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().

 {
   // get the measurements and their errors, use information updated with tsos if improving
   // (GF: Also for measurements or only for errors or do the former not change?)
   // GF 10/2008: I doubt that it makes sense to update the hit with the tsos here:
   //             That is an analytical extrapolation and not the best guess of the real 
   //             track state on the module, but the latter should be better to get the best
   //             hit uncertainty estimate!
 
    // FIXME FIXME  CLONE
   auto newHitPtr =  hitPtr;
 //  TransientTrackingRecHit::ConstRecHitPointer newHitPtr(hitPtr->canImproveWithTrack() ?
 //                          hitPtr->clone(updatedTsos) : hitPtr);
 
   const LocalPoint localMeasurement    = newHitPtr->localPosition();
   const LocalError localMeasurementCov = newHitPtr->localPositionError();
   
   theMeasurements[iRow]   = localMeasurement.x();
   theMeasurements[iRow+1] = localMeasurement.y();
   theMeasurementsCov[iRow][iRow]     = localMeasurementCov.xx();
   theMeasurementsCov[iRow][iRow+1]   = localMeasurementCov.xy();
   theMeasurementsCov[iRow+1][iRow+1] = localMeasurementCov.yy();
   // GF: Should be a loop once the hit dimension is not hardcoded as nMeasPerHit (to be checked):
   // for (int i = 0; i < hitPtr->dimension(); ++i) {
   //   theMeasurements[iRow+i]   = hitPtr->parameters()[i]; // fixme: parameters() is by value!
   //   for (int j = i; j < hitPtr->dimension(); ++j) {
   //     theMeasurementsCov[iRow+i][iRow+j] = hitPtr->parametersError()[i][j];
   //   }
   // }
 }

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

References ReferenceTrajectoryBase::theTrajectoryPositions.

Referenced by construct().

 {
   // get the local coordinates of the reference trajectory
   const AlgebraicVector localPosition(projection * mixedLocalParams);
   theTrajectoryPositions[iRow] = localPosition[0];
   theTrajectoryPositions[iRow+1] = localPosition[1];
   // GF: Should be a loop once the hit dimension is not hardcoded as nMeasPerHit (to be checked):
   // for (int j = 0; j < projection.num_col(); ++j) {
   //   theTrajectoryPositions[iRow+j] = localPosition[j];
   // }
 }

AlgebraicMatrix ReferenceTrajectory::getHitProjectionMatrix ( const TransientTrackingRecHit::ConstRecHitPointer & recHit ) const

protected

first (generic) helper to get the projection matrix

Definition at line 925 of file ReferenceTrajectory.cc.

References edm::hlt::Exception.

Referenced by construct().

 {
   if (this->useRecHit(hitPtr)) {
     // check which templated non-member function to call:
     switch (hitPtr->dimension()) {
     case 1:
       return getHitProjectionMatrixT<1>(hitPtr);
     case 2:
       return getHitProjectionMatrixT<2>(hitPtr);
     case 3:
       return getHitProjectionMatrixT<3>(hitPtr);
     case 4:
       return getHitProjectionMatrixT<4>(hitPtr);
     case 5:
       return getHitProjectionMatrixT<5>(hitPtr);
     default:
       throw cms::Exception("ReferenceTrajectory::getHitProjectionMatrix")
         << "Unexpected hit dimension: " << hitPtr->dimension() << "\n";
     }
   }
   // invalid or (to please compiler) unknown dimension
   return AlgebraicMatrix(2, 5, 0); // get size from ???
 }

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

References data-class-funcs::H, N, KfComponentsHolder::projection(), alignCSCRings::r, and KfComponentsHolder::setup().

 {
   // define variables that will be used to setup the KfComponentsHolder
   // (their allocated memory is needed by 'hitPtr->getKfComponents(..)'
   // ProjectMatrix<double,5,N>  pf; // not needed
   typename AlgebraicROOTObject<N,5>::Matrix H; 
   typename AlgebraicROOTObject<N>::Vector r, rMeas; 
   typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas;
   // input for the holder - but dummy is OK here to just get the projection matrix:
   const AlgebraicVector5 dummyPars;
   const AlgebraicSymMatrix55 dummyErr;
   ProjectMatrix<double,5,N> dummyProjFunc;
 
   // setup the holder with the correct dimensions and get the values
   KfComponentsHolder holder;
   holder.setup<N>(&r, &V, &H, &dummyProjFunc, &rMeas, &VMeas, dummyPars, dummyErr);
   hitPtr->getKfComponents(holder);
 
   return asHepMatrix<N,5>(holder.projection<N>());
 }

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

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

Referenced by construct().

 {
   // propagate to next layer
   //AnalyticalPropagator aPropagator(magField, propDir);
   // Hard coded RungeKutta instead Analytical (avoid bias in TEC), but
   // work around TrackPropagation/RungeKutta/interface/RKTestPropagator.h and
   // http://www.parashift.com/c++-faq-lite/strange-inheritance.html#faq-23.9
   defaultRKPropagator::Product  rkprod(magField, propDir); //double tolerance = 5.e-5)
   Propagator &aPropagator = rkprod.propagator;
   const std::pair<TrajectoryStateOnSurface, double> tsosWithPath =
     aPropagator.propagateWithPath(previousTsos, newSurface);
 
   // stop if propagation wasn't successful
   if (!tsosWithPath.first.isValid()) return false;
   
   nextStep = tsosWithPath.second;
   // calculate derivative of reference-track parameters on the actual layer w.r.t. the ones
   // on the previous layer (both in global coordinates)
   const AnalyticalCurvilinearJacobian aJacobian(previousTsos.globalParameters(), 
                         tsosWithPath.first.globalPosition(),
                         tsosWithPath.first.globalMomentum(),
                         tsosWithPath.second);
   const AlgebraicMatrix curvilinearJacobian = asHepMatrix<5,5>(aJacobian.jacobian());
 
   // jacobian of the track parameters on the previous layer for local->global transformation
   const JacobianLocalToCurvilinear startTrafo(previousSurface, previousTsos.localParameters(), *magField);
   const AlgebraicMatrix localToCurvilinear = asHepMatrix<5>(startTrafo.jacobian());
     
   // jacobian of the track parameters on the actual layer for global->local transformation
   const JacobianCurvilinearToLocal endTrafo(newSurface, tsosWithPath.first.localParameters(), *magField);
   const AlgebraicMatrix curvilinearToLocal = asHepMatrix<5>(endTrafo.jacobian());
   
   // compute derivative of reference-track parameters on the actual layer w.r.t. the ones on
   // the previous layer (both in their local representation)
   newCurvlinJacobian = curvilinearJacobian;
   newJacobian = curvilinearToLocal * curvilinearJacobian * localToCurvilinear;
   newTsos     = tsosWithPath.first;
 
   return true;
 }

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

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

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

   {
     return ( dir == alongMomentum ) ?
       SurfaceSideDefinition::beforeSurface :
       SurfaceSideDefinition::afterSurface;
   }

Public Types

Public Member Functions

Protected Member Functions

Additional Inherited Members

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation