ReferenceTrajectories/src/TwoBodyDecayTrajectory.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0004 #include "DataFormats/GeometrySurface/interface/Surface.h"
0005 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectory.h"
0006 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0007 #include "DataFormats/Math/interface/Error.h"
0008 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayParameters.h"
0009
0010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0011
0012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0013
0014 TwoBodyDecayTrajectory::TwoBodyDecayTrajectory(const TwoBodyDecayTrajectoryState& tsos,
0015                                                const ConstRecHitCollection& recHits,
0016                                                const MagneticField* magField,
0017                                                const reco::BeamSpot& beamSpot,
0018                                                const ReferenceTrajectoryBase::Config& config)
0019     : ReferenceTrajectoryBase(
0020           TwoBodyDecayParameters::dimension,
0021           recHits.first.size() + recHits.second.size(),
0022           (config.materialEffects >= breakPoints) ? 2 * (recHits.first.size() + recHits.second.size()) - 4 : 0,
0023           (config.materialEffects >= breakPoints) ? 2 * (recHits.first.size() + recHits.second.size()) - 3 : 1),
0024       materialEffects_(config.materialEffects),
0025       propDir_(config.propDir),
0026       useRefittedState_(config.useRefittedState),
0027       constructTsosWithErrors_(config.constructTsosWithErrors)
0028
0029 {
0030   if (config.hitsAreReverse) {
0031     TransientTrackingRecHit::ConstRecHitContainer::const_reverse_iterator itRecHits;
0032     ConstRecHitCollection fwdRecHits;
0033
0034     fwdRecHits.first.reserve(recHits.first.size());
0035     for (itRecHits = recHits.first.rbegin(); itRecHits != recHits.first.rend(); ++itRecHits) {
0036       fwdRecHits.first.push_back(*itRecHits);
0037     }
0038
0039     fwdRecHits.second.reserve(recHits.second.size());
0040     for (itRecHits = recHits.second.rbegin(); itRecHits != recHits.second.rend(); ++itRecHits) {
0041       fwdRecHits.second.push_back(*itRecHits);
0042     }
0043
0044     theValidityFlag = this->construct(tsos, fwdRecHits, magField, beamSpot);
0045   } else {
0046     theValidityFlag = this->construct(tsos, recHits, magField, beamSpot);
0047   }
0048 }
0049
0050 TwoBodyDecayTrajectory::TwoBodyDecayTrajectory(void)
0051     : ReferenceTrajectoryBase(0, 0, 0, 0),
0052       materialEffects_(none),
0053       propDir_(anyDirection),
0054       useRefittedState_(false),
0055       constructTsosWithErrors_(false) {}
0056
0057 bool TwoBodyDecayTrajectory::construct(const TwoBodyDecayTrajectoryState& state,
0058                                        const ConstRecHitCollection& recHits,
0059                                        const MagneticField* field,
0060                                        const reco::BeamSpot& beamSpot) {
0061   const TwoBodyDecayTrajectoryState::TsosContainer& tsos = state.trajectoryStates(useRefittedState_);
0062   const TwoBodyDecayTrajectoryState::Derivatives& deriv = state.derivatives();
0063   double mass = state.particleMass();
0064
0065   //
0066   // first track
0067   //
0068
0069   // construct a trajectory (hits should be already in correct order)
0070   ReferenceTrajectoryBase::Config config(materialEffects_, propDir_, mass);
0071   config.useBeamSpot = false;
0072   config.hitsAreReverse = false;
0073
0074   ReferenceTrajectory trajectory1(tsos.first, recHits.first, field, beamSpot, config);
0075
0076   // check if construction of trajectory was successful
0077   if (!trajectory1.isValid())
0078     return false;
0079
0080   //
0081   // second track
0082   //
0083
0084   ReferenceTrajectory trajectory2(tsos.second, recHits.second, field, beamSpot, config);
0085
0086   if (!trajectory2.isValid())
0087     return false;
0088
0089   //
0090   // combine both tracks
0091   //
0092   unsigned int nLocal = deriv.first.num_row();
0093   unsigned int nTbd = deriv.first.num_col();
0094
0095   if (materialEffects_ >= localGBL) {
0096     // GBL trajectory inputs
0097     // convert to Eigen::MatrixXd
0098     Eigen::MatrixXd tbdToLocal1{nLocal, nTbd};
0099     for (unsigned int row = 0; row < nLocal; ++row) {
0100       for (unsigned int col = 0; col < nTbd; ++col) {
0101         tbdToLocal1(row, col) = deriv.first[row][col];
0102       }
0103     }
0104     // add first body
0105     theGblInput.push_back(
0106         std::make_pair(trajectory1.gblInput().front().first, trajectory1.gblInput().front().second * tbdToLocal1));
0107     // convert to Eigen::MatrixXd
0108     Eigen::MatrixXd tbdToLocal2{nLocal, nTbd};
0109     for (unsigned int row = 0; row < nLocal; ++row) {
0110       for (unsigned int col = 0; col < nTbd; ++col) {
0111         tbdToLocal2(row, col) = deriv.second[row][col];
0112       }
0113     }
0114     // add second body
0115     theGblInput.push_back(
0116         std::make_pair(trajectory2.gblInput().front().first, trajectory2.gblInput().front().second * tbdToLocal2));
0117     // add virtual mass measurement
0118     theGblExtDerivatives.resize(1, nTbd);
0119     theGblExtDerivatives.setZero();
0120     theGblExtDerivatives(0, TwoBodyDecayParameters::mass) = 1.0;
0121     theGblExtMeasurements.resize(1);
0122     theGblExtMeasurements(0) = state.primaryMass() - state.decayParameters()[TwoBodyDecayParameters::mass];
0123     theGblExtPrecisions.resize(1);
0124     theGblExtPrecisions(0) = 1.0 / (state.primaryWidth() * state.primaryWidth());
0125     // nominal field
0126     theNomField = trajectory1.nominalField();
0127   } else {
0128     unsigned int nHitMeas1 = trajectory1.numberOfHitMeas();
0129     unsigned int nVirtualMeas1 = trajectory1.numberOfVirtualMeas();
0130     unsigned int nPar1 = trajectory1.numberOfPar();
0131     unsigned int nVirtualPar1 = trajectory1.numberOfVirtualPar();
0132
0133     // derivatives of the trajectory w.r.t. to the decay parameters
0134     AlgebraicMatrix fullDeriv1 = trajectory1.derivatives().sub(1, nHitMeas1 + nVirtualMeas1, 1, nLocal) *
0135                                  trajectory1.localToTrajectory() * deriv.first;
0136
0137     unsigned int nHitMeas2 = trajectory2.numberOfHitMeas();
0138     unsigned int nVirtualMeas2 = trajectory2.numberOfVirtualMeas();
0139     unsigned int nPar2 = trajectory2.numberOfPar();
0140     unsigned int nVirtualPar2 = trajectory2.numberOfVirtualPar();
0141
0142     AlgebraicMatrix fullDeriv2 = trajectory2.derivatives().sub(1, nHitMeas2 + nVirtualMeas2, 1, nLocal) *
0143                                  trajectory2.localToTrajectory() * deriv.second;
0144
0145     theNumberOfRecHits.first = recHits.first.size();
0146     theNumberOfRecHits.second = recHits.second.size();
0147
0148     theNumberOfHits = trajectory1.numberOfHits() + trajectory2.numberOfHits();
0149     theNumberOfPars = nPar1 + nPar2;
0150     theNumberOfVirtualPars = nVirtualPar1 + nVirtualPar2;
0151     theNumberOfVirtualMeas = nVirtualMeas1 + nVirtualMeas2 + 1;  // add virtual mass measurement
0152
0153     // hit measurements from trajectory 1
0154     int rowOffset = 1;
0155     int colOffset = 1;
0156     theDerivatives.sub(rowOffset, colOffset, fullDeriv1.sub(1, nHitMeas1, 1, nTbd));
0157     colOffset += nTbd;
0158     theDerivatives.sub(
0159         rowOffset, colOffset, trajectory1.derivatives().sub(1, nHitMeas1, nLocal + 1, nPar1 + nVirtualPar1));
0160     // hit measurements from trajectory 2
0161     rowOffset += nHitMeas1;
0162     colOffset = 1;
0163     theDerivatives.sub(rowOffset, colOffset, fullDeriv2.sub(1, nHitMeas2, 1, nTbd));
0164     colOffset += (nPar1 + nVirtualPar1 + nTbd - nLocal);
0165     theDerivatives.sub(
0166         rowOffset, colOffset, trajectory2.derivatives().sub(1, nHitMeas2, nLocal + 1, nPar2 + nVirtualPar2));
0167     // MS measurements from trajectory 1
0168     rowOffset += nHitMeas2;
0169     colOffset = 1;
0170     theDerivatives.sub(rowOffset, colOffset, fullDeriv1.sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1, 1, nTbd));
0171     colOffset += nTbd;
0172     theDerivatives.sub(
0173         rowOffset,
0174         colOffset,
0175         trajectory1.derivatives().sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1, nLocal + 1, nPar1 + nVirtualPar1));
0176     // MS measurements from trajectory 2
0177     rowOffset += nVirtualMeas1;
0178     colOffset = 1;
0179     theDerivatives.sub(rowOffset, colOffset, fullDeriv2.sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2, 1, nTbd));
0180     colOffset += (nPar1 + nVirtualPar1 + nTbd - nLocal);
0181     theDerivatives.sub(
0182         rowOffset,
0183         colOffset,
0184         trajectory2.derivatives().sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2, nLocal + 1, nPar2 + nVirtualPar2));
0185
0186     theMeasurements.sub(1, trajectory1.measurements().sub(1, nHitMeas1));
0187     theMeasurements.sub(nHitMeas1 + 1, trajectory2.measurements().sub(1, nHitMeas2));
0188     theMeasurements.sub(nHitMeas1 + nHitMeas2 + 1,
0189                         trajectory1.measurements().sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1));
0190     theMeasurements.sub(nHitMeas1 + nHitMeas2 + nVirtualMeas1 + 1,
0191                         trajectory2.measurements().sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2));
0192
0193     theMeasurementsCov.sub(1, trajectory1.measurementErrors().sub(1, nHitMeas1));
0194     theMeasurementsCov.sub(nHitMeas1 + 1, trajectory2.measurementErrors().sub(1, nHitMeas2));
0195     theMeasurementsCov.sub(nHitMeas1 + nHitMeas2 + 1,
0196                            trajectory1.measurementErrors().sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1));
0197     theMeasurementsCov.sub(nHitMeas1 + nHitMeas2 + nVirtualMeas1 + 1,
0198                            trajectory2.measurementErrors().sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2));
0199
0200     theTrajectoryPositions.sub(1, trajectory1.trajectoryPositions());
0201     theTrajectoryPositions.sub(nHitMeas1 + 1, trajectory2.trajectoryPositions());
0202
0203     theTrajectoryPositionCov =
0204         state.decayParameters().covariance().similarity(theDerivatives.sub(1, nHitMeas1 + nHitMeas2, 1, 9));
0205
0206     theParameters = state.decayParameters().parameters();
0207
0208     // add virtual mass measurement
0209     rowOffset += nVirtualMeas2;
0210     int indMass = rowOffset - 1;
0211     theMeasurements[indMass] = state.primaryMass() - state.decayParameters()[TwoBodyDecayParameters::mass];
0212     theMeasurementsCov[indMass][indMass] = state.primaryWidth() * state.primaryWidth();
0213     theDerivatives[indMass][TwoBodyDecayParameters::mass] = 1.0;
0214   }
0215
0216   theRecHits.insert(theRecHits.end(), recHits.first.begin(), recHits.first.end());
0217   theRecHits.insert(theRecHits.end(), recHits.second.begin(), recHits.second.end());
0218
0219   if (constructTsosWithErrors_) {
0220     constructTsosVecWithErrors(trajectory1, trajectory2, field);
0221   } else {
0222     theTsosVec.insert(theTsosVec.end(), trajectory1.trajectoryStates().begin(), trajectory1.trajectoryStates().end());
0223
0224     theTsosVec.insert(theTsosVec.end(), trajectory2.trajectoryStates().begin(), trajectory2.trajectoryStates().end());
0225   }
0226
0227   return true;
0228 }
0229
0230 void TwoBodyDecayTrajectory::constructTsosVecWithErrors(const ReferenceTrajectory& traj1,
0231                                                         const ReferenceTrajectory& traj2,
0232                                                         const MagneticField* field) {
0233   int iTsos = 0;
0234
0235   std::vector<TrajectoryStateOnSurface>::const_iterator itTsos;
0236
0237   for (itTsos = traj1.trajectoryStates().begin(); itTsos != traj1.trajectoryStates().end(); itTsos++) {
0238     constructSingleTsosWithErrors(*itTsos, iTsos, field);
0239     iTsos++;
0240   }
0241
0242   for (itTsos = traj2.trajectoryStates().begin(); itTsos != traj2.trajectoryStates().end(); itTsos++) {
0243     constructSingleTsosWithErrors(*itTsos, iTsos, field);
0244     iTsos++;
0245   }
0246 }
0247
0248 void TwoBodyDecayTrajectory::constructSingleTsosWithErrors(const TrajectoryStateOnSurface& tsos,
0249                                                            int iTsos,
0250                                                            const MagneticField* field) {
0251   AlgebraicSymMatrix55 localErrors;
0252   const double coeff = 1e-2;
0253
0254   double invP = tsos.localParameters().signedInverseMomentum();
0255   LocalVector p = tsos.localParameters().momentum();
0256
0257   // rough estimate for the errors of q/p, dx/dz and dy/dz, assuming that
0258   // sigma(px) = sigma(py) = sigma(pz) = coeff*p.
0259   float dpinv = coeff * (fabs(p.x()) + fabs(p.y()) + fabs(p.z())) * invP * invP;
0260   float dxdir = coeff * (fabs(p.x()) + fabs(p.z())) / p.z() / p.z();
0261   float dydir = coeff * (fabs(p.y()) + fabs(p.z())) / p.z() / p.z();
0262   localErrors[0][0] = dpinv * dpinv;
0263   localErrors[1][1] = dxdir * dxdir;
0264   localErrors[2][2] = dydir * dydir;
0265
0266   // exact values for the errors on local x and y
0267   localErrors[3][3] = theTrajectoryPositionCov[nMeasPerHit * iTsos][nMeasPerHit * iTsos];
0268   localErrors[3][4] = theTrajectoryPositionCov[nMeasPerHit * iTsos][nMeasPerHit * iTsos + 1];
0269   localErrors[4][4] = theTrajectoryPositionCov[nMeasPerHit * iTsos + 1][nMeasPerHit * iTsos + 1];
0270
0271   // construct tsos with local errors
0272   theTsosVec[iTsos] = TrajectoryStateOnSurface(
0273       tsos.localParameters(), LocalTrajectoryError(localErrors), tsos.surface(), field, tsos.surfaceSide());
0274 }