ReferenceTrajectories/src/DualReferenceTrajectory.cc

0001
0002 #include "Alignment/ReferenceTrajectories/interface/DualReferenceTrajectory.h"
0003
0004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0005
0006 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0007 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0008
0009 #include "Alignment/ReferenceTrajectories/interface/ReferenceTrajectory.h"
0010
0011 DualReferenceTrajectory::DualReferenceTrajectory(const TrajectoryStateOnSurface& tsos,
0012                                                  const ConstRecHitContainer& forwardRecHits,
0013                                                  const ConstRecHitContainer& backwardRecHits,
0014                                                  const MagneticField* magField,
0015                                                  const reco::BeamSpot& beamSpot,
0016                                                  const ReferenceTrajectoryBase::Config& config)
0017     : ReferenceTrajectoryBase(tsos.localParameters().mixedFormatVector().kSize,
0018                               numberOfUsedRecHits(forwardRecHits) + numberOfUsedRecHits(backwardRecHits) - 1,
0019                               0,
0020                               0),
0021       mass_(config.mass),
0022       materialEffects_(config.materialEffects),
0023       propDir_(config.propDir),
0024       useBeamSpot_(config.useBeamSpot) {
0025   theValidityFlag = this->construct(tsos, forwardRecHits, backwardRecHits, magField, beamSpot);
0026 }
0027
0028 DualReferenceTrajectory::DualReferenceTrajectory(unsigned int nPar,
0029                                                  unsigned int nHits,
0030                                                  const ReferenceTrajectoryBase::Config& config)
0031     : ReferenceTrajectoryBase(nPar, nHits, 0, 0),
0032       mass_(config.mass),
0033       materialEffects_(config.materialEffects),
0034       propDir_(config.propDir),
0035       useBeamSpot_(config.useBeamSpot) {}
0036
0037 bool DualReferenceTrajectory::construct(const TrajectoryStateOnSurface& refTsos,
0038                                         const ConstRecHitContainer& forwardRecHits,
0039                                         const ConstRecHitContainer& backwardRecHits,
0040                                         const MagneticField* magField,
0041                                         const reco::BeamSpot& beamSpot) {
0042   if (materialEffects_ >= breakPoints)
0043     throw cms::Exception("BadConfig") << "[DualReferenceTrajectory::construct] Wrong MaterialEffects: "
0044                                       << materialEffects_;
0045
0046   ReferenceTrajectoryBase* fwdTraj = construct(refTsos, forwardRecHits, magField, beamSpot);
0047
0048   // set flag for opposite direction to true
0049   ReferenceTrajectoryBase* bwdTraj = construct(refTsos, backwardRecHits, magField, beamSpot, true);
0050
0051   if (!(fwdTraj->isValid() && bwdTraj->isValid())) {
0052     delete fwdTraj;
0053     delete bwdTraj;
0054     return false;
0055   }
0056
0057   //
0058   // Combine both reference trajactories to a dual reference trajectory
0059   //
0060
0061   const std::vector<TrajectoryStateOnSurface>& fwdTsosVec = fwdTraj->trajectoryStates();
0062   const std::vector<TrajectoryStateOnSurface>& bwdTsosVec = bwdTraj->trajectoryStates();
0063   theTsosVec.insert(theTsosVec.end(), fwdTsosVec.begin(), fwdTsosVec.end());
0064   theTsosVec.insert(theTsosVec.end(), ++bwdTsosVec.begin(), bwdTsosVec.end());
0065
0066   const ConstRecHitContainer& fwdRecHits = fwdTraj->recHits();
0067   const ConstRecHitContainer& bwdRecHits = bwdTraj->recHits();
0068   theRecHits.insert(theRecHits.end(), fwdRecHits.begin(), fwdRecHits.end());
0069   theRecHits.insert(theRecHits.end(), ++bwdRecHits.begin(), bwdRecHits.end());
0070
0071   theParameters = extractParameters(refTsos);
0072
0073   unsigned int nParam = theNumberOfPars;
0074   unsigned int nFwdMeas = fwdTraj->numberOfHitMeas();
0075   unsigned int nBwdMeas = bwdTraj->numberOfHitMeas();
0076   unsigned int nFwdBP = fwdTraj->numberOfVirtualMeas();
0077   unsigned int nBwdBP = bwdTraj->numberOfVirtualMeas();
0078   unsigned int nMeas = nFwdMeas + nBwdMeas - nMeasPerHit;
0079
0080   theMeasurements.sub(1, fwdTraj->measurements().sub(1, nFwdMeas));
0081   theMeasurements.sub(nFwdMeas + 1, bwdTraj->measurements().sub(nMeasPerHit + 1, nBwdMeas));
0082
0083   theMeasurementsCov.sub(1, fwdTraj->measurementErrors().sub(1, nFwdMeas));
0084   theMeasurementsCov.sub(nFwdMeas + 1, bwdTraj->measurementErrors().sub(nMeasPerHit + 1, nBwdMeas));
0085
0086   theTrajectoryPositions.sub(1, fwdTraj->trajectoryPositions());
0087   theTrajectoryPositions.sub(nFwdMeas + 1, bwdTraj->trajectoryPositions().sub(nMeasPerHit + 1, nBwdMeas));
0088
0089   theTrajectoryPositionCov.sub(1, fwdTraj->trajectoryPositionErrors());
0090   theTrajectoryPositionCov.sub(nFwdMeas + 1, bwdTraj->trajectoryPositionErrors().sub(nMeasPerHit + 1, nBwdMeas));
0091
0092   theDerivatives.sub(1, 1, fwdTraj->derivatives().sub(1, nFwdMeas, 1, nParam));
0093   theDerivatives.sub(nFwdMeas + 1, 1, bwdTraj->derivatives().sub(nMeasPerHit + 1, nBwdMeas, 1, nParam));
0094
0095   // for the break points
0096   // DUAL with break points makes no sense: (MS) correlations between the two parts are lost !
0097   if (nFwdBP > 0) {
0098     theMeasurements.sub(nMeas + 1, fwdTraj->measurements().sub(nFwdMeas + 1, nFwdMeas + nFwdBP));
0099     theMeasurementsCov.sub(nMeas + 1, fwdTraj->measurementErrors().sub(nFwdMeas + 1, nFwdMeas + nFwdBP));
0100     theDerivatives.sub(1, nParam + 1, fwdTraj->derivatives().sub(1, nFwdMeas, nParam + 1, nParam + nFwdBP));
0101     theDerivatives.sub(nMeas + 1,
0102                        nParam + 1,
0103                        fwdTraj->derivatives().sub(nFwdMeas + 1, nFwdMeas + nFwdBP, nParam + 1, nParam + nFwdBP));
0104   }
0105   if (nBwdBP > 0) {
0106     theMeasurements.sub(nMeas + nFwdBP + 1, bwdTraj->measurements().sub(nBwdMeas + 1, nBwdMeas + nBwdBP));
0107     theMeasurementsCov.sub(nMeas + nFwdBP + 1, bwdTraj->measurementErrors().sub(nBwdMeas + 1, nBwdMeas + nBwdBP));
0108     theDerivatives.sub(nFwdMeas + 1,
0109                        nParam + nFwdBP + 1,
0110                        bwdTraj->derivatives().sub(nMeasPerHit + 1, nBwdMeas, nParam + 1, nParam + nBwdBP));
0111     theDerivatives.sub(nMeas + nFwdBP + 1,
0112                        nParam + nFwdBP + 1,
0113                        bwdTraj->derivatives().sub(nBwdMeas + 1, nBwdMeas + nBwdBP, nParam + 1, nParam + nBwdBP));
0114   }
0115
0116   delete fwdTraj;
0117   delete bwdTraj;
0118
0119   return true;
0120 }
0121
0122 ReferenceTrajectory* DualReferenceTrajectory::construct(const TrajectoryStateOnSurface& referenceTsos,
0123                                                         const ConstRecHitContainer& recHits,
0124                                                         const MagneticField* magField,
0125                                                         const reco::BeamSpot& beamSpot,
0126                                                         const bool revertDirection) const {
0127   if (materialEffects_ >= breakPoints)
0128     throw cms::Exception("BadConfig") << "[DualReferenceTrajectory::construct] Wrong MaterialEffects: "
0129                                       << materialEffects_;
0130
0131   ReferenceTrajectoryBase::Config config(
0132       materialEffects_, (revertDirection ? oppositeDirection(propDir_) : propDir_), mass_);
0133   config.useBeamSpot = useBeamSpot_;
0134   config.hitsAreReverse = false;
0135   return new ReferenceTrajectory(referenceTsos, recHits, magField, beamSpot, config);
0136 }
0137
0138 AlgebraicVector DualReferenceTrajectory::extractParameters(const TrajectoryStateOnSurface& referenceTsos) const {
0139   return asHepVector<5>(referenceTsos.localParameters().mixedFormatVector());
0140 }