MuonTrajectoryUpdator.cc

Go to the documentation of this file.

#include "RecoMuon/TrackingTools/interface/MuonTrajectoryUpdator.h"
#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"

#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"

#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include <algorithm>

using namespace edm;
using namespace std;

MuonTrajectoryUpdator::MuonTrajectoryUpdator(const edm::ParameterSet& par, NavigationDirection fitDirection)
    : theFitDirection(fitDirection) {
  // The max allowed chi2 to accept a rechit in the fit
  theMaxChi2 = par.getParameter<double>("MaxChi2");
  theEstimator = new Chi2MeasurementEstimator(theMaxChi2);

  // The KF updator
  theUpdator = new KFUpdator();

  // The granularity
  theGranularity = par.getParameter<int>("Granularity");

  // Rescale the error of the first state?
  theRescaleErrorFlag = par.getParameter<bool>("RescaleError");

  if (theRescaleErrorFlag)
    // The rescale factor
    theRescaleFactor = par.getParameter<double>("RescaleErrorFactor");

  // Use the invalid hits?
  useInvalidHits = par.getParameter<bool>("UseInvalidHits");

  // Flag needed for the rescaling
  theFirstTSOSFlag = true;

  // Exlude RPC from the fit?
  theRPCExFlag = par.getParameter<bool>("ExcludeRPCFromFit");
}

MuonTrajectoryUpdator::MuonTrajectoryUpdator(NavigationDirection fitDirection, double chi2, int granularity)
    : theMaxChi2(chi2), theGranularity(granularity), theFitDirection(fitDirection) {
  theEstimator = new Chi2MeasurementEstimator(theMaxChi2);

  // The KF updator
  theUpdator = new KFUpdator();
}

MuonTrajectoryUpdator::~MuonTrajectoryUpdator() {
  delete theEstimator;
  delete theUpdator;
}

void MuonTrajectoryUpdator::makeFirstTime() { theFirstTSOSFlag = true; }

pair<bool, TrajectoryStateOnSurface> MuonTrajectoryUpdator::update(const TrajectoryMeasurement* measurement,
                                                                   Trajectory& trajectory,
                                                                   const Propagator* propagator) {
  const std::string metname = "Muon|RecoMuon|MuonTrajectoryUpdator";

  MuonPatternRecoDumper muonDumper;

  // Status of the updating
  bool updated = false;

  if (!measurement)
    return pair<bool, TrajectoryStateOnSurface>(updated, TrajectoryStateOnSurface());

  // measurement layer
  const DetLayer* detLayer = measurement->layer();

  // these are the 4D segment for the CSC/DT and a point for the RPC
  TransientTrackingRecHit::ConstRecHitPointer muonRecHit = measurement->recHit();

  // The KFUpdator takes TransientTrackingRecHits as arg.
  TransientTrackingRecHit::ConstRecHitContainer recHitsForFit =
      MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, theGranularity);

  // sort the container in agreement with the porpagation direction
  sort(recHitsForFit, detLayer);

  TrajectoryStateOnSurface lastUpdatedTSOS = measurement->predictedState();

  LogTrace(metname) << "Number of rechits for the fit: " << recHitsForFit.size() << endl;

  TransientTrackingRecHit::ConstRecHitContainer::iterator recHit;
  for (recHit = recHitsForFit.begin(); recHit != recHitsForFit.end(); ++recHit) {
    if ((*recHit)->isValid()) {
      // propagate the TSOS onto the rechit plane
      TrajectoryStateOnSurface propagatedTSOS = propagateState(lastUpdatedTSOS, measurement, *recHit, propagator);

      if (propagatedTSOS.isValid()) {
        pair<bool, double> thisChi2 = estimator()->estimate(propagatedTSOS, *((*recHit).get()));

        LogTrace(metname) << "Estimation for Kalman Fit. Chi2: " << thisChi2.second;

        // Is an RPC hit? Prepare the variable to possibly exluding it from the fit
        bool wantIncludeThisHit = true;
        if (theRPCExFlag && (*recHit)->geographicalId().det() == DetId::Muon &&
            (*recHit)->geographicalId().subdetId() == MuonSubdetId::RPC) {
          wantIncludeThisHit = false;
          LogTrace(metname)
              << "This is an RPC hit and the present configuration is such that it will be excluded from the fit";
        }

        // The Chi2 cut was already applied in the estimator, which
        // returns 0 if the chi2 is bigger than the cut defined in its
        // constructor
        if (thisChi2.first) {
          updated = true;
          if (wantIncludeThisHit) {  // This split is a trick to have the RPC hits counted as updatable (in used chamber counting), while are not actually included in the fit when the proper obtion is activated.

            LogTrace(metname) << endl
                              << "     Kalman Start"
                              << "\n"
                              << "\n";
            LogTrace(metname) << "  Meas. Position : " << (**recHit).globalPosition() << "\n"
                              << "  Pred. Position : " << propagatedTSOS.globalPosition()
                              << "  Pred Direction : " << propagatedTSOS.globalDirection() << endl;

            if (theRescaleErrorFlag && theFirstTSOSFlag) {
              propagatedTSOS.rescaleError(theRescaleFactor);
              theFirstTSOSFlag = false;
            }

            lastUpdatedTSOS = measurementUpdator()->update(propagatedTSOS, *((*recHit).get()));

            if (!lastUpdatedTSOS.isValid()) {
              edm::LogInfo(metname) << "Invalid last TSOS, will skip RecHit ";
              lastUpdatedTSOS = propagatedTSOS;  // Revert update
              continue;
            }

            LogTrace(metname) << "  Fit   Position : " << lastUpdatedTSOS.globalPosition()
                              << "  Fit  Direction : " << lastUpdatedTSOS.globalDirection() << "\n"
                              << "  Fit position radius : " << lastUpdatedTSOS.globalPosition().perp()
                              << "filter updated" << endl;

            LogTrace(metname) << muonDumper.dumpTSOS(lastUpdatedTSOS);

            LogTrace(metname) << "\n\n     Kalman End"
                              << "\n"
                              << "\n";

            TrajectoryMeasurement&& updatedMeasurement =
                updateMeasurement(propagatedTSOS, lastUpdatedTSOS, *recHit, thisChi2.second, detLayer, measurement);
            // FIXME: check!
            trajectory.push(std::move(updatedMeasurement), thisChi2.second);
          } else {
            LogTrace(metname) << "  Compatible RecHit with good chi2 but made with RPC when it was decided to not "
                                 "include it in the fit"
                              << "  --> trajectory NOT updated, invalid RecHit added." << endl;

            MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
                MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
            invalidRhPtr->invalidateHit();
            TrajectoryMeasurement invalidRhMeasurement(
                propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
            trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
          }
        } else {
          if (useInvalidHits) {
            LogTrace(metname) << "  Compatible RecHit with too large chi2"
                              << "  --> trajectory NOT updated, invalid RecHit added." << endl;

            MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
                MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
            invalidRhPtr->invalidateHit();
            TrajectoryMeasurement invalidRhMeasurement(
                propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
            trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
          }
        }
      }
    }
  }
  recHitsForFit.clear();
  return pair<bool, TrajectoryStateOnSurface>(updated, lastUpdatedTSOS);
}

TrajectoryStateOnSurface MuonTrajectoryUpdator::propagateState(
    const TrajectoryStateOnSurface& state,
    const TrajectoryMeasurement* measurement,
    const TransientTrackingRecHit::ConstRecHitPointer& current,
    const Propagator* propagator) const {
  const TransientTrackingRecHit::ConstRecHitPointer& mother = measurement->recHit();

  if (current->geographicalId() == mother->geographicalId())
    return measurement->predictedState();

  const TrajectoryStateOnSurface tsos = propagator->propagate(state, current->det()->surface());
  return tsos;
}

// FIXME: would I a different threatment for the two prop dirrections??
TrajectoryMeasurement MuonTrajectoryUpdator::updateMeasurement(const TrajectoryStateOnSurface& propagatedTSOS,
                                                               const TrajectoryStateOnSurface& lastUpdatedTSOS,
                                                               const TransientTrackingRecHit::ConstRecHitPointer& recHit,
                                                               const double& chi2,
                                                               const DetLayer* detLayer,
                                                               const TrajectoryMeasurement* initialMeasurement) {
  return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, recHit, chi2, detLayer);

  //   // FIXME: put a better check! One could fit in first out-in and then in - out
  //   if(propagator()->propagationDirection() == alongMomentum)
  //     return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS,
  //                 recHit,thisChi2.second,detLayer);

  //   // FIXME: Check this carefully!!
  //   else if(propagator()->propagationDirection() == oppositeToMomentum)
  //     return TrajectoryMeasurement(initialMeasurement->forwardPredictedState(),
  //                 propagatedTSOS, lastUpdatedTSOS,
  //                 recHit,thisChi2.second,detLayer);
  //   else{
  //     LogError("MuonTrajectoryUpdator::updateMeasurement") <<"Wrong propagation direction!!";
  //   }
}

void MuonTrajectoryUpdator::sort(TransientTrackingRecHit::ConstRecHitContainer& recHitsForFit,
                                 const DetLayer* detLayer) {
  if (detLayer->subDetector() == GeomDetEnumerators::DT) {
    if (fitDirection() == insideOut)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorInOut());
    else if (fitDirection() == outsideIn)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorOutIn());
    else
      LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
  }

  else if (detLayer->subDetector() == GeomDetEnumerators::CSC) {
    if (fitDirection() == insideOut)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
    else if (fitDirection() == outsideIn)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
    else
      LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
  }

  else if (detLayer->subDetector() == GeomDetEnumerators::GEM) {
    if (fitDirection() == insideOut)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
    else if (fitDirection() == outsideIn)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
    else
      LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
  } else if (detLayer->subDetector() == GeomDetEnumerators::ME0) {
    if (fitDirection() == insideOut)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
    else if (fitDirection() == outsideIn)
      stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
    else
      LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
  }
}