DualReferenceTrajectory.cc

Go to the documentation of this file.

#include "Alignment/ReferenceTrajectories/interface/DualReferenceTrajectory.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"

#include "DataFormats/CLHEP/interface/AlgebraicObjects.h" 
#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"

#include "Alignment/ReferenceTrajectories/interface/ReferenceTrajectory.h"

DualReferenceTrajectory::DualReferenceTrajectory(const TrajectoryStateOnSurface& tsos,
                                                 const ConstRecHitContainer& forwardRecHits,
                                                 const ConstRecHitContainer& backwardRecHits,
                                                 const MagneticField* magField,
                                                 const reco::BeamSpot& beamSpot,
                                                 const ReferenceTrajectoryBase::Config& config) :
  ReferenceTrajectoryBase(tsos.localParameters().mixedFormatVector().kSize,
                          numberOfUsedRecHits(forwardRecHits) + numberOfUsedRecHits(backwardRecHits) - 1,
                          0, 0),
  mass_(config.mass),
  materialEffects_(config.materialEffects),
  propDir_(config.propDir),
  useBeamSpot_(config.useBeamSpot)
{
  theValidityFlag = this->construct(tsos, forwardRecHits, backwardRecHits, magField, beamSpot);
}



DualReferenceTrajectory::DualReferenceTrajectory(unsigned int nPar, unsigned int nHits,
                         const ReferenceTrajectoryBase::Config& config)
  : ReferenceTrajectoryBase(nPar, nHits, 0, 0),
  mass_(config.mass),
  materialEffects_(config.materialEffects),
  propDir_(config.propDir),
  useBeamSpot_(config.useBeamSpot)
{}


bool DualReferenceTrajectory::construct(const TrajectoryStateOnSurface &refTsos,
                                        const ConstRecHitContainer &forwardRecHits,
                                        const ConstRecHitContainer &backwardRecHits,
                                        const MagneticField *magField,
                                        const reco::BeamSpot &beamSpot)
{
  if (materialEffects_ >= breakPoints)  throw cms::Exception("BadConfig")
    << "[DualReferenceTrajectory::construct] Wrong MaterialEffects: " << materialEffects_;
    
  ReferenceTrajectoryBase* fwdTraj = construct(refTsos, forwardRecHits, magField, beamSpot);

  // set flag for opposite direction to true
  ReferenceTrajectoryBase* bwdTraj = construct(refTsos, backwardRecHits, magField, beamSpot, true);

  if ( !( fwdTraj->isValid() && bwdTraj->isValid() ) )
  {
    delete fwdTraj;
    delete bwdTraj;
    return false;
  }

  //
  // Combine both reference trajactories to a dual reference trajectory
  //

  const std::vector<TrajectoryStateOnSurface>& fwdTsosVec = fwdTraj->trajectoryStates();
  const std::vector<TrajectoryStateOnSurface>& bwdTsosVec = bwdTraj->trajectoryStates();
  theTsosVec.insert( theTsosVec.end(), fwdTsosVec.begin(), fwdTsosVec.end() );
  theTsosVec.insert( theTsosVec.end(), ++bwdTsosVec.begin(), bwdTsosVec.end() );

  const ConstRecHitContainer &fwdRecHits = fwdTraj->recHits();
  const ConstRecHitContainer &bwdRecHits = bwdTraj->recHits();
  theRecHits.insert( theRecHits.end(), fwdRecHits.begin(), fwdRecHits.end() );
  theRecHits.insert( theRecHits.end(), ++bwdRecHits.begin(), bwdRecHits.end() );

  theParameters = extractParameters( refTsos );
  
  unsigned int nParam   = theNumberOfPars;
  unsigned int nFwdMeas = fwdTraj->numberOfHitMeas();
  unsigned int nBwdMeas = bwdTraj->numberOfHitMeas();
  unsigned int nFwdBP   = fwdTraj->numberOfVirtualMeas();
  unsigned int nBwdBP   = bwdTraj->numberOfVirtualMeas();
  unsigned int nMeas    = nFwdMeas+nBwdMeas-nMeasPerHit; 
       
  theMeasurements.sub( 1, fwdTraj->measurements().sub( 1, nFwdMeas ) );
  theMeasurements.sub( nFwdMeas+1, bwdTraj->measurements().sub( nMeasPerHit+1, nBwdMeas ) );
    
  theMeasurementsCov.sub( 1, fwdTraj->measurementErrors().sub( 1, nFwdMeas ) );
  theMeasurementsCov.sub( nFwdMeas+1, bwdTraj->measurementErrors().sub( nMeasPerHit+1, nBwdMeas ) );
  
  theTrajectoryPositions.sub( 1, fwdTraj->trajectoryPositions() );
  theTrajectoryPositions.sub( nFwdMeas+1, bwdTraj->trajectoryPositions().sub( nMeasPerHit+1, nBwdMeas ) );

  theTrajectoryPositionCov.sub( 1, fwdTraj->trajectoryPositionErrors() );
  theTrajectoryPositionCov.sub( nFwdMeas+1, bwdTraj->trajectoryPositionErrors().sub( nMeasPerHit+1, nBwdMeas ) );

  theDerivatives.sub( 1, 1, fwdTraj->derivatives().sub( 1, nFwdMeas, 1, nParam ) );
  theDerivatives.sub( nFwdMeas+1, 1, bwdTraj->derivatives().sub( nMeasPerHit+1, nBwdMeas, 1, nParam ) );
  
// for the break points 
// DUAL with break points makes no sense: (MS) correlations between the two parts are lost ! 
  if (nFwdBP>0 )
  {
    theMeasurements.sub( nMeas+1, fwdTraj->measurements().sub( nFwdMeas+1, nFwdMeas+nFwdBP ) );
    theMeasurementsCov.sub( nMeas+1, fwdTraj->measurementErrors().sub( nFwdMeas+1, nFwdMeas+nFwdBP ) );  
    theDerivatives.sub( 1, nParam+1, fwdTraj->derivatives().sub( 1, nFwdMeas, nParam+1, nParam+nFwdBP ) );
    theDerivatives.sub( nMeas+1, nParam+1, fwdTraj->derivatives().sub( nFwdMeas+1, nFwdMeas+nFwdBP, nParam+1, nParam+nFwdBP ) );
  }  
  if (nBwdBP>0 )
  {
    theMeasurements.sub( nMeas+nFwdBP+1, bwdTraj->measurements().sub( nBwdMeas+1, nBwdMeas+nBwdBP ) );  
    theMeasurementsCov.sub( nMeas+nFwdBP+1, bwdTraj->measurementErrors().sub( nBwdMeas+1, nBwdMeas+nBwdBP ) );  
    theDerivatives.sub( nFwdMeas+1, nParam+nFwdBP+1, bwdTraj->derivatives().sub( nMeasPerHit+1, nBwdMeas, nParam+1, nParam+nBwdBP ) );
    theDerivatives.sub( nMeas+nFwdBP+1, nParam+nFwdBP+1, bwdTraj->derivatives().sub( nBwdMeas+1, nBwdMeas+nBwdBP, nParam+1, nParam+nBwdBP ) );
  }
    
  delete fwdTraj;
  delete bwdTraj;

  return true; 
}


ReferenceTrajectory*
DualReferenceTrajectory::construct(const TrajectoryStateOnSurface &referenceTsos, 
                   const ConstRecHitContainer &recHits,
                   const MagneticField *magField,
                   const reco::BeamSpot &beamSpot,
                   const bool revertDirection) const
{
  if (materialEffects_ >= breakPoints)  throw cms::Exception("BadConfig")
    << "[DualReferenceTrajectory::construct] Wrong MaterialEffects: " << materialEffects_;
  
  ReferenceTrajectoryBase::Config config(materialEffects_,
                                         (revertDirection ? oppositeDirection(propDir_): propDir_),
                                         mass_);
  config.useBeamSpot = useBeamSpot_;
  config.hitsAreReverse = false;
  return new ReferenceTrajectory(referenceTsos, recHits, magField, beamSpot, config);
}


AlgebraicVector
DualReferenceTrajectory::extractParameters(const TrajectoryStateOnSurface &referenceTsos) const
{
  return asHepVector<5>( referenceTsos.localParameters().mixedFormatVector() );
}