File indexing completed on 2024-04-06 11:57:21
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
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
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
0096
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 }