Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-16 22:19:00

0001 #include "RecoMuon/TrackingTools/plugins/MuonErrorMatrixAdjuster.h"
0002 
0003 #include "TString.h"
0004 #include "TMath.h"
0005 #include <MagneticField/Records/interface/IdealMagneticFieldRecord.h>
0006 #include <MagneticField/Engine/interface/MagneticField.h>
0007 
0008 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
0009 #include "RecoMuon/TrackingTools/interface/MuonErrorMatrix.h"
0010 
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 
0013 #include <TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h>
0014 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0015 
0016 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0017 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0018 
0019 MuonErrorMatrixAdjuster::MuonErrorMatrixAdjuster(const edm::ParameterSet& iConfig)
0020     : theFieldToken{esConsumes()}, theHttopoToken{esConsumes()} {
0021   theCategory = "MuonErrorMatrixAdjuster";
0022   theInstanceName = iConfig.getParameter<std::string>("instanceName");
0023   //register your products
0024   produces<reco::TrackCollection>(theInstanceName);
0025   produces<TrackingRecHitCollection>();
0026   produces<reco::TrackExtraCollection>();
0027 
0028   theTrackLabel = iConfig.getParameter<edm::InputTag>("trackLabel");
0029   consumes<reco::TrackCollection>(theTrackLabel);
0030   theRescale = iConfig.getParameter<bool>("rescale");
0031 
0032   auto matrixProvider_pset = iConfig.getParameter<edm::ParameterSet>("errorMatrix_pset");
0033 
0034   theMatrixProvider = std::make_unique<MuonErrorMatrix>(matrixProvider_pset);
0035 }
0036 
0037 MuonErrorMatrixAdjuster::~MuonErrorMatrixAdjuster() {
0038   // do anything here that needs to be done at desctruction time
0039   // (e.g. close files, deallocate resources etc.)
0040 }
0041 
0042 //
0043 // member functions
0044 //
0045 
0046 //take the error matrix and rescale it or just replace it
0047 reco::TrackBase::CovarianceMatrix MuonErrorMatrixAdjuster::fix_cov_matrix(
0048     const reco::TrackBase::CovarianceMatrix& error_matrix, const GlobalVector& momentum) {
0049   //CovarianceMatrix is template for SMatrix
0050   reco::TrackBase::CovarianceMatrix revised_matrix(theMatrixProvider->get(momentum));
0051 
0052   if (theRescale) {
0053     //rescale old error matrix up by a factor
0054     multiply(revised_matrix, error_matrix);
0055   }
0056   return revised_matrix;
0057 }
0058 
0059 void MuonErrorMatrixAdjuster::multiply(reco::TrackBase::CovarianceMatrix& revised_matrix,
0060                                        const reco::TrackBase::CovarianceMatrix& scale_matrix) {
0061   //scale term by term the matrix
0062   // the true type of the matrix is such that [i][j] is the same memory object as [j][i]: looping i:0-5, j:0-5 double multiply the terms
0063   // need to loop only on i:0-5, j:i-5
0064   for (int i = 0; i != 5; i++) {
0065     for (int j = i; j != 5; j++) {
0066       revised_matrix(i, j) *= scale_matrix(i, j);
0067     }
0068   }
0069 }
0070 bool MuonErrorMatrixAdjuster::divide(reco::TrackBase::CovarianceMatrix& num_matrix,
0071                                      const reco::TrackBase::CovarianceMatrix& denom_matrix) {
0072   //divide term by term the matrix
0073   // the true type of the matrix is such that [i][j] is the same memory object as [j][i]: looping i:0-5, j:0-5 double multiply the terms
0074   // need to loop only on i:0-5, j:i-5
0075   for (int i = 0; i != 5; i++) {
0076     for (int j = i; j != 5; j++) {
0077       if (denom_matrix(i, j) == 0)
0078         return false;
0079       num_matrix(i, j) /= denom_matrix(i, j);
0080     }
0081   }
0082   return true;
0083 }
0084 
0085 reco::Track MuonErrorMatrixAdjuster::makeTrack(const reco::Track& recotrack_orig, const FreeTrajectoryState& PCAstate) {
0086   //get the parameters of the track so I can reconstruct it
0087   double chi2 = recotrack_orig.chi2();
0088   double ndof = recotrack_orig.ndof();
0089   const math::XYZPoint& refpos = recotrack_orig.referencePoint();
0090   const math::XYZVector& mom = recotrack_orig.momentum();
0091   int charge = recotrack_orig.charge();
0092 
0093   reco::TrackBase::CovarianceMatrix covariance_matrix =
0094       fix_cov_matrix(recotrack_orig.covariance(), PCAstate.momentum());
0095 
0096   LogDebug(theCategory) << "chi2: " << chi2 << "\n ndof: " << ndof << "\n refpos: " << refpos << "\n mom: " << mom
0097                         << "\n charge: " << charge << "\n covariance:\n"
0098                         << recotrack_orig.covariance() << "\n replaced by:\n"
0099                         << covariance_matrix;
0100 
0101   return reco::Track(chi2, ndof, refpos, mom, charge, covariance_matrix);
0102 }
0103 
0104 reco::TrackExtra* MuonErrorMatrixAdjuster::makeTrackExtra(const reco::Track& recotrack_orig,
0105                                                           reco::Track& recotrack,
0106                                                           reco::TrackExtraCollection& TEcol) {
0107   //get the 5x5 matrix of recotrack/recotrack_orig
0108   reco::TrackBase::CovarianceMatrix scale_matrix(recotrack.covariance());
0109   if (!divide(scale_matrix, recotrack_orig.covariance())) {
0110     edm::LogError(theCategory) << "original track error matrix has term ==0... skipping.";
0111     return nullptr;
0112   }
0113 
0114   const reco::TrackExtraRef& trackExtra_orig = recotrack_orig.extra();
0115   if (trackExtra_orig.isNull()) {
0116     edm::LogError(theCategory) << "original track has no track extra... skipping.";
0117     return nullptr;
0118   }
0119 
0120   //copy the outer state. rescaling the error matrix
0121   reco::TrackBase::CovarianceMatrix outerCov(trackExtra_orig->outerStateCovariance());
0122   multiply(outerCov, scale_matrix);
0123 
0124   //copy the inner state, rescaling the error matrix
0125   reco::TrackBase::CovarianceMatrix innerCov(trackExtra_orig->innerStateCovariance());
0126   multiply(innerCov, scale_matrix);
0127 
0128   //put the trackExtra
0129   TEcol.push_back(reco::TrackExtra(trackExtra_orig->outerPosition(),
0130                                    trackExtra_orig->outerMomentum(),
0131                                    true,
0132                                    trackExtra_orig->innerPosition(),
0133                                    trackExtra_orig->innerMomentum(),
0134                                    true,
0135                                    outerCov,
0136                                    trackExtra_orig->outerDetId(),
0137                                    innerCov,
0138                                    trackExtra_orig->innerDetId(),
0139                                    trackExtra_orig->seedDirection()));
0140 
0141   //add a reference to the trackextra on the track
0142   recotrack.setExtra(edm::Ref<reco::TrackExtraCollection>(theRefprodTE, theTEi++));
0143 
0144   //return the reference to the last inserted then
0145   return &(TEcol.back());
0146 }
0147 
0148 bool MuonErrorMatrixAdjuster::attachRecHits(const reco::Track& recotrack_orig,
0149                                             reco::Track& recotrack,
0150                                             reco::TrackExtra& trackextra,
0151                                             TrackingRecHitCollection& RHcol,
0152                                             const TrackerTopology& ttopo) {
0153   //loop over the hits of the original track
0154   trackingRecHit_iterator recHit = recotrack_orig.recHitsBegin();
0155   auto const firstHitIndex = theRHi;
0156   for (; recHit != recotrack_orig.recHitsEnd(); ++recHit) {
0157     //clone it. this is meandatory
0158     TrackingRecHit* hit = (*recHit)->clone();
0159 
0160     //put it on the new track
0161     recotrack.appendHitPattern(*hit, ttopo);
0162     //copy them in the new collection
0163     RHcol.push_back(hit);
0164     ++theRHi;
0165 
0166   }  //loop over original rechits
0167   //do something with the trackextra
0168   trackextra.setHits(theRefprodRH, firstHitIndex, theRHi - firstHitIndex);
0169 
0170   return true;  //if nothing fails
0171 }
0172 
0173 bool MuonErrorMatrixAdjuster::selectTrack(const reco::Track& recotrack_orig) { return true; }
0174 
0175 // ------------ method called to produce the data  ------------
0176 void MuonErrorMatrixAdjuster::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0177   using namespace edm;
0178 
0179   //open a collection of track
0180   edm::Handle<reco::TrackCollection> tracks;
0181   iEvent.getByLabel(theTrackLabel, tracks);
0182   LogDebug(theCategory) << "considering: " << tracks->size() << " uncorrected reco::Track from the event.("
0183                         << theTrackLabel << ")";
0184 
0185   //get the mag field
0186   theField = iSetup.getHandle(theFieldToken);
0187 
0188   const TrackerTopology& ttopo = iSetup.getData(theHttopoToken);
0189 
0190   //prepare the output collection
0191   auto Toutput = std::make_unique<reco::TrackCollection>();
0192   auto TRHoutput = std::make_unique<TrackingRecHitCollection>();
0193   auto TEoutput = std::make_unique<reco::TrackExtraCollection>();
0194   theRefprodTE = iEvent.getRefBeforePut<reco::TrackExtraCollection>();
0195   theTEi = 0;
0196   theRefprodRH = iEvent.getRefBeforePut<TrackingRecHitCollection>();
0197   theRHi = 0;
0198 
0199   for (unsigned int it = 0; it != tracks->size(); it++) {
0200     const reco::Track& recotrack_orig = (*tracks)[it];
0201     FreeTrajectoryState PCAstate = trajectoryStateTransform::initialFreeState(recotrack_orig, theField.product());
0202     if (PCAstate.position().mag() == 0) {
0203       edm::LogError(theCategory) << "invalid state from track initial state in " << theTrackLabel << ". skipping.";
0204       continue;
0205     }
0206 
0207     //create a reco::Track
0208     reco::Track recotrack = makeTrack(recotrack_orig, PCAstate);
0209 
0210     //make a selection on the create reco::Track
0211     if (!selectTrack(recotrack))
0212       continue;
0213 
0214     Toutput->push_back(recotrack);
0215     reco::Track& recotrackref = Toutput->back();
0216 
0217     //build the track extra
0218     reco::TrackExtra* extra = makeTrackExtra(recotrack_orig, recotrackref, *TEoutput);
0219     if (!extra) {
0220       edm::LogError(theCategory) << "cannot create the track extra for this track.";
0221       //pop the inserted track
0222       Toutput->pop_back();
0223       continue;
0224     }
0225 
0226     //attach the collection of rechits
0227     if (!attachRecHits(recotrack_orig, recotrackref, *extra, *TRHoutput, ttopo)) {
0228       edm::LogError(theCategory) << "cannot attach any rechits on this track";
0229       //pop the inserted track
0230       Toutput->pop_back();
0231       //pop the track extra
0232       TEoutput->pop_back();
0233       theTEi--;
0234       continue;
0235     }
0236 
0237   }  //loop over the original reco tracks
0238 
0239   LogDebug(theCategory) << "writing: " << Toutput->size() << " corrected reco::Track to the event.";
0240 
0241   iEvent.put(std::move(Toutput), theInstanceName);
0242   iEvent.put(std::move(TEoutput));
0243   iEvent.put(std::move(TRHoutput));
0244 }