CosmicMuonSmoother.cc

Go to the documentation of this file.

#include "RecoMuon/CosmicMuonProducer/interface/CosmicMuonSmoother.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

using namespace edm;
using namespace std;

//
// constructor
//
CosmicMuonSmoother::CosmicMuonSmoother(const ParameterSet& par, const MuonServiceProxy *service) : theService(service) {

  theUpdator     = new KFUpdator;
  theUtilities   = new CosmicMuonUtilities; 
  theEstimator   = new Chi2MeasurementEstimator(200.0);
  thePropagatorAlongName = par.getParameter<string>("PropagatorAlong");
  thePropagatorOppositeName = par.getParameter<string>("PropagatorOpposite");
  theErrorRescaling = par.getParameter<double>("RescalingFactor");

  category_ = "Muon|RecoMuon|CosmicMuon|CosmicMuonSmoother";
}

//
// destructor
//
CosmicMuonSmoother::~CosmicMuonSmoother() {

  if ( theUpdator ) delete theUpdator;
  if ( theUtilities ) delete theUtilities;
  if ( theEstimator ) delete theEstimator;

}


//
// fit and smooth trajectory
//
Trajectory CosmicMuonSmoother::trajectory(const Trajectory& t) const {
   std::vector<Trajectory> && fitted = fit(t);
   if (fitted.empty()) return  Trajectory();
   std::vector<Trajectory> && smoothed = smooth(fitted);
   return  smoothed.empty() ? Trajectory() : smoothed.front(); 
}

//
// fit and smooth trajectory 
//
std::vector<Trajectory> CosmicMuonSmoother::trajectories(const TrajectorySeed& seed,
                                               const ConstRecHitContainer& hits, 
                                               const TrajectoryStateOnSurface& firstPredTsos) const {

  if ( hits.empty() ||!firstPredTsos.isValid() ) return vector<Trajectory>();

  LogTrace(category_)<< "trajectory begin (seed hits tsos)";

  TrajectoryStateOnSurface firstTsos = firstPredTsos;
  firstTsos.rescaleError(theErrorRescaling);

  LogTrace(category_)<< "first TSOS: "<<firstTsos;

  vector<Trajectory> fitted = fit(seed, hits, firstTsos);
  LogTrace(category_)<< "fitted: "<<fitted.size();
  vector<Trajectory> smoothed = smooth(fitted);
  LogTrace(category_)<< "smoothed: "<<smoothed.size();

  return smoothed;

}

//
// fit trajectory
//
vector<Trajectory> CosmicMuonSmoother::fit(const Trajectory& t) const {

  if ( t.empty() ) return vector<Trajectory>();

  LogTrace(category_)<< "fit begin (trajectory) ";

  TrajectoryStateOnSurface firstTsos = initialState(t); 
  if ( !firstTsos.isValid() ) {
    LogTrace(category_)<< "Error: firstTsos invalid. ";
    return vector<Trajectory>();
  }

  LogTrace(category_)<< "firstTsos: "<<firstTsos;

  ConstRecHitContainer hits = t.recHits();
  LogTrace(category_)<< "hits: "<<hits.size();
  LogTrace(category_)<<"hit front" <<hits.front()->globalPosition()<< " hit back" 
    <<hits.back()->globalPosition();

  sortHitsAlongMom(hits, firstTsos);

  LogTrace(category_)<<"after sorting hit front" <<hits.front()->globalPosition()<< " hit back" 
    <<hits.back()->globalPosition();

  return fit(t.seed(), hits, firstTsos);

}


//
// fit trajectory
//
vector<Trajectory> CosmicMuonSmoother::fit(const TrajectorySeed& seed,
                              const ConstRecHitContainer& hits, 
                          const TrajectoryStateOnSurface& firstPredTsos) const {

  LogTrace(category_)<< "fit begin (seed, hit, tsos).";

  if ( hits.empty() ) {
    LogTrace(category_)<< "Error: empty hits container.";
    return vector<Trajectory>();
  }

  Trajectory myTraj(seed, alongMomentum);

  TrajectoryStateOnSurface predTsos(firstPredTsos);
  LogTrace(category_)<< "first pred TSOS: "<<predTsos;

  if ( !predTsos.isValid() ) {
    LogTrace(category_)<< "Error: firstTsos invalid.";
    return vector<Trajectory>();
  }
  TrajectoryStateOnSurface currTsos;

  if ( hits.front()->isValid() ) {

    // FIXME  FIXME  CLONE !!!
   //  TrackingRecHit::RecHitPointer preciseHit = hits.front()->clone(predTsos);
   auto preciseHit = hits.front();
 
    LogTrace(category_)<<"first hit is at det "<< hits.front()->det()->surface().position();

    currTsos = theUpdator->update(predTsos, *preciseHit);
    if (!currTsos.isValid()){
      edm::LogWarning(category_)
    <<"an updated state is not valid. killing the trajectory.";
      return vector<Trajectory>();
    }
    myTraj.push(TrajectoryMeasurement(predTsos, currTsos, hits.front(),
                theEstimator->estimate(predTsos, *hits.front()).second));
    if ( currTsos.isValid() )  LogTrace(category_)<< "first curr TSOS: "<<currTsos;

  } else {

    currTsos = predTsos;
    myTraj.push(TrajectoryMeasurement(predTsos, hits.front()));
  }
  //const TransientTrackingRecHit& firsthit = *hits.front();

  for ( ConstRecHitContainer::const_iterator ihit = hits.begin() + 1; 
        ihit != hits.end(); ++ihit ) {

    if ( !(**ihit).isValid() ) {
      LogTrace(category_)<< "Error: invalid hit.";
      continue;
    }
    if (currTsos.isValid() && currTsos.globalMomentum().mag2() > 1e-18f)  {
     LogTrace(category_)<<"current pos "<<currTsos.globalPosition()
                       <<"mom "<<currTsos.globalMomentum();
    } else {
      LogTrace(category_)<<"current state invalid or momentum is too low";
      //logically, there's no way out: can't expect a valid result out of an invalid state
      LogTrace(category_)
    <<"Input state is not valid. This loop over hits is doomed: breaking out";
      break;
    }

    predTsos = propagatorAlong()->propagate(currTsos, (**ihit).det()->surface());
    LogTrace(category_)<<"predicted state propagate directly "<<predTsos.isValid();

    if ( !predTsos.isValid() ) {
      LogTrace(category_)<<"step-propagating from "<<currTsos.globalPosition() <<" to position: "<<(*ihit)->globalPosition();
      predTsos = theUtilities->stepPropagate(currTsos, (*ihit), *propagatorAlong());
    }
    if ( !predTsos.isValid() && (fabs(theService->magneticField()->inTesla(GlobalPoint(0,0,0)).z()) < 0.01) && (theService->propagator("StraightLinePropagator").isValid() ) ) {
      LogTrace(category_)<<"straight-line propagating from "<<currTsos.globalPosition() <<" to position: "<<(*ihit)->globalPosition();
      predTsos = theService->propagator("StraightLinePropagator")->propagate(currTsos, (**ihit).det()->surface());
    }
    if (predTsos.isValid())  {
      LogTrace(category_)<<"predicted pos "<<predTsos.globalPosition()
                         <<"mom "<<predTsos.globalMomentum();
    } else {
      LogTrace(category_)<<"predicted state invalid";
    }
    if ( !predTsos.isValid() ) {
      LogTrace(category_)<< "Error: predTsos is still invalid forward fit.";
//      return vector<Trajectory>();
      continue;
    } else if ( (**ihit).isValid() ) {
          // FIXME  FIXME  CLONE !!!
      // update  (FIXME!)
      // TransientTrackingRecHit::RecHitPointer preciseHit = (**ihit).clone(predTsos);
      auto preciseHit = *ihit;
 
      if ( !preciseHit->isValid() ) {
        currTsos = predTsos;
        myTraj.push(TrajectoryMeasurement(predTsos, *ihit));
      } else {
        currTsos = theUpdator->update(predTsos, *preciseHit);
    if (!currTsos.isValid()){
      edm::LogWarning(category_)
        <<"an updated state is not valid. killing the trajectory.";
      return vector<Trajectory>();
    }
        myTraj.push(TrajectoryMeasurement(predTsos, currTsos, preciseHit,
                       theEstimator->estimate(predTsos, *preciseHit).second));
      }
    } else {
      currTsos = predTsos;
      myTraj.push(TrajectoryMeasurement(predTsos, *ihit));
    }

  }

  std::vector<TrajectoryMeasurement> mytms = myTraj.measurements();
  LogTrace(category_)<<"fit result "<<mytms.size();
  for (std::vector<TrajectoryMeasurement>::const_iterator itm = mytms.begin();
       itm != mytms.end(); ++itm ) {
       LogTrace(category_)<<"updated pos "<<itm->updatedState().globalPosition()
                       <<"mom "<<itm->updatedState().globalMomentum();
       }


  return vector<Trajectory>(1, myTraj);

}


//
// smooth trajectory
//
vector<Trajectory> CosmicMuonSmoother::smooth(const vector<Trajectory>& tc) const {

  vector<Trajectory> result; 

  for ( vector<Trajectory>::const_iterator it = tc.begin(); it != tc.end(); ++it ) {
    vector<Trajectory> smoothed = smooth(*it);
    result.insert(result.end(), smoothed.begin(), smoothed.end());
  }

  return result;

}


//
// smooth trajectory
//
vector<Trajectory> CosmicMuonSmoother::smooth(const Trajectory& t) const {

  if ( t.empty() ) {
    LogTrace(category_)<< "Error: smooth: empty trajectory.";
    return vector<Trajectory>();
  }

  Trajectory myTraj(t.seed(), oppositeToMomentum);

  vector<TrajectoryMeasurement> avtm = t.measurements();

  if ( avtm.size() < 2 ) {
    LogTrace(category_)<< "Error: smooth: too little TM. ";
    return vector<Trajectory>();
  }

  TrajectoryStateOnSurface predTsos = avtm.back().forwardPredictedState();
  predTsos.rescaleError(theErrorRescaling);

  if ( !predTsos.isValid() ) {
    LogTrace(category_)<< "Error: smooth: first TSOS from back invalid. ";
    return vector<Trajectory>();
  }

  TrajectoryStateOnSurface currTsos;

  // first smoothed TrajectoryMeasurement is last fitted
  if ( avtm.back().recHit()->isValid() ) {
    currTsos = theUpdator->update(predTsos, (*avtm.back().recHit()));
    if (!currTsos.isValid()){
      edm::LogWarning(category_)
    <<"an updated state is not valid. killing the trajectory.";
      return vector<Trajectory>();
    }
    myTraj.push(TrajectoryMeasurement(avtm.back().forwardPredictedState(),
           predTsos,
           avtm.back().updatedState(),
           avtm.back().recHit(),
           avtm.back().estimate()//,
           /*avtm.back().layer()*/), 
            avtm.back().estimate());

  } else {
    currTsos = predTsos;
    myTraj.push(TrajectoryMeasurement(avtm.back().forwardPredictedState(),
           avtm.back().recHit()//,
           /*avtm.back().layer()*/));

  }

  TrajectoryStateCombiner combiner;


  for ( vector<TrajectoryMeasurement>::reverse_iterator itm = avtm.rbegin() + 1; 
        itm != avtm.rend() - 1; ++itm ) {

   if (currTsos.isValid())  {
     LogTrace(category_)<<"current pos "<<currTsos.globalPosition()
                       <<"mom "<<currTsos.globalMomentum();
    } else {
      LogTrace(category_)<<"current state invalid";
    }

    predTsos = propagatorOpposite()->propagate(currTsos,(*itm).recHit()->det()->surface());

    if ( !predTsos.isValid() ) {
      LogTrace(category_)<<"step-propagating from "<<currTsos.globalPosition() <<" to position: "<<(*itm).recHit()->globalPosition();
      predTsos = theUtilities->stepPropagate(currTsos, (*itm).recHit(), *propagatorOpposite());
    }
   if (predTsos.isValid())  {
      LogTrace(category_)<<"predicted pos "<<predTsos.globalPosition()
                       <<"mom "<<predTsos.globalMomentum();
    } else {
      LogTrace(category_)<<"predicted state invalid";
    }

    if ( !predTsos.isValid() ) {
      LogTrace(category_)<< "Error: predTsos is still invalid backward smooth.";
      return vector<Trajectory>();
    //  continue;
    } else if ( (*itm).recHit()->isValid() ) {
      //update
      currTsos = theUpdator->update(predTsos, (*(*itm).recHit()));
      if (!currTsos.isValid()){
    edm::LogWarning(category_)
      <<"an updated state is not valid. killing the trajectory.";
    return vector<Trajectory>();
      }
      TrajectoryStateOnSurface combTsos = combiner(predTsos, (*itm).forwardPredictedState());
      if ( !combTsos.isValid() ) {
         LogTrace(category_)<< "Error: smooth: combining pred TSOS failed. ";
         return vector<Trajectory>();
      }

      TrajectoryStateOnSurface smooTsos = combiner((*itm).updatedState(), predTsos);

      if ( !smooTsos.isValid() ) {
         LogTrace(category_)<< "Error: smooth: combining smooth TSOS failed. ";
         return vector<Trajectory>();
      }

      myTraj.push(TrajectoryMeasurement((*itm).forwardPredictedState(),
             predTsos,
             smooTsos,
             (*itm).recHit(),
             theEstimator->estimate(combTsos, (*(*itm).recHit())).second//,
             /*(*itm).layer()*/),
             (*itm).estimate());
    } else {
      currTsos = predTsos;
      TrajectoryStateOnSurface combTsos = combiner(predTsos, (*itm).forwardPredictedState());
      
      if ( !combTsos.isValid() ) {
         LogTrace(category_)<< "Error: smooth: combining TSOS failed. ";
         return vector<Trajectory>();
      }

      myTraj.push(TrajectoryMeasurement((*itm).forwardPredictedState(),
             predTsos,
             combTsos,
             (*itm).recHit()//,
             /*(*itm).layer()*/));
    }
  }

  // last smoothed TrajectoryMeasurement is last filtered
  predTsos = propagatorOpposite()->propagate(currTsos, avtm.front().recHit()->det()->surface());
  
  if ( !predTsos.isValid() ){
    LogTrace(category_)<< "Error: last predict TSOS failed, use original one. ";
 //    return vector<Trajectory>();
      myTraj.push(TrajectoryMeasurement(avtm.front().forwardPredictedState(),
                   avtm.front().recHit()));
  } else  {
    if ( avtm.front().recHit()->isValid() ) {
      //update
      currTsos = theUpdator->update(predTsos, (*avtm.front().recHit()));
      if (currTsos.isValid())
      myTraj.push(TrajectoryMeasurement(avtm.front().forwardPredictedState(),
           predTsos,
           currTsos,
           avtm.front().recHit(),
           theEstimator->estimate(predTsos, (*avtm.front().recHit())).second//,
           /*avtm.front().layer()*/),
            avtm.front().estimate());
    }
  }
  LogTrace(category_)<< "myTraj foundHits. "<<myTraj.foundHits();

  if (myTraj.foundHits() >= 3)
    return vector<Trajectory>(1, myTraj);
  else {
     LogTrace(category_)<< "Error: smooth: No enough hits in trajctory. ";
     return vector<Trajectory>();
  } 

}

TrajectoryStateOnSurface CosmicMuonSmoother::initialState(const Trajectory& t) const {
  if ( t.empty() ) return TrajectoryStateOnSurface();

  if ( !t.firstMeasurement().updatedState().isValid() || !t.lastMeasurement().updatedState().isValid() )  return TrajectoryStateOnSurface();

  TrajectoryStateOnSurface result;

  bool beamhaloFlag = ( t.firstMeasurement().updatedState().globalMomentum().eta() > 4.0 || t.lastMeasurement().updatedState().globalMomentum().eta() > 4.0 ); 

  if ( !beamhaloFlag ) { //initialState is the top one
     if (t.firstMeasurement().updatedState().globalPosition().y() > t.lastMeasurement().updatedState().globalPosition().y()) {
     result = t.firstMeasurement().updatedState();
     } else {
       result = t.lastMeasurement().updatedState();
     } 
     if (result.globalMomentum().y()> 1.0 ) //top tsos should pointing down
       theUtilities->reverseDirection(result,&*theService->magneticField());
  } else {
     if ( t.firstMeasurement().updatedState().globalPosition().z() * t.lastMeasurement().updatedState().globalPosition().z() >0 ) { //same side
       if ( fabs(t.firstMeasurement().updatedState().globalPosition().z()) > fabs(t.lastMeasurement().updatedState().globalPosition().z()) ) {
         result = t.firstMeasurement().updatedState();
       } else {
         result = t.lastMeasurement().updatedState();
       }
     } else { //different sides

       if ( fabs(t.firstMeasurement().updatedState().globalPosition().eta()) > fabs(t.lastMeasurement().updatedState().globalPosition().eta()) ) {
         result = t.firstMeasurement().updatedState();
       } else {
         result = t.lastMeasurement().updatedState();
       }
     }

  }

  return result;

}

void CosmicMuonSmoother::sortHitsAlongMom(ConstRecHitContainer& hits, const TrajectoryStateOnSurface& tsos) const {

    if (hits.size() < 2) return;
    float dis1 = (hits.front()->globalPosition() - tsos.globalPosition()).mag();
    float dis2 = (hits.back()->globalPosition() - tsos.globalPosition()).mag();

    if ( dis1 > dis2 )
      std::reverse(hits.begin(),hits.end());

    return;
}