CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

CosmicMuonSmoother Class Reference

#include <CosmicMuonSmoother.h>

Inheritance diagram for CosmicMuonSmoother:
TrajectorySmoother

List of all members.

Public Member Functions

virtual CosmicMuonSmootherclone () const
 CosmicMuonSmoother (const edm::ParameterSet &, const MuonServiceProxy *service)
Chi2MeasurementEstimatorestimator () const
std::vector< Trajectoryfit (const Trajectory &) const
std::vector< Trajectoryfit (const TrajectorySeed &seed, const ConstRecHitContainer &hits, const TrajectoryStateOnSurface &firstPredTsos) const
const PropagatorpropagatorAlong () const
const PropagatorpropagatorOpposite () const
virtual TrajectoryContainer trajectories (const Trajectory &traj) const
virtual TrajectoryContainer trajectories (const TrajectorySeed &seed, const ConstRecHitContainer &hits, const TrajectoryStateOnSurface &firstPredTsos) const
 refit trajectory
Trajectory trajectory (const Trajectory &) const
KFUpdatorupdator () const
CosmicMuonUtilitiesutilities () const
virtual ~CosmicMuonSmoother ()

Private Member Functions

TrajectoryStateOnSurface initialState (const Trajectory &) const
std::vector< Trajectorysmooth (const Trajectory &) const
std::vector< Trajectorysmooth (const std::vector< Trajectory > &) const
void sortHitsAlongMom (ConstRecHitContainer &hits, const TrajectoryStateOnSurface &) const

Private Attributes

std::string category_
double theErrorRescaling
Chi2MeasurementEstimatortheEstimator
std::string thePropagatorAlongName
std::string thePropagatorOppositeName
const MuonServiceProxytheService
KFUpdatortheUpdator
CosmicMuonUtilitiestheUtilities

Detailed Description

Definition at line 37 of file CosmicMuonSmoother.h.


Constructor & Destructor Documentation

CosmicMuonSmoother::CosmicMuonSmoother ( const edm::ParameterSet par,
const MuonServiceProxy *  service 
)
CosmicMuonSmoother::~CosmicMuonSmoother ( ) [virtual]

Definition at line 49 of file CosmicMuonSmoother.cc.

References theEstimator, theUpdator, and theUtilities.

                                        {

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

}

Member Function Documentation

virtual CosmicMuonSmoother* CosmicMuonSmoother::clone ( void  ) const [inline, virtual]

Implements TrajectorySmoother.

Definition at line 50 of file CosmicMuonSmoother.h.

References CosmicMuonSmoother().

                                            {
    return new CosmicMuonSmoother(*this);
  }
Chi2MeasurementEstimator* CosmicMuonSmoother::estimator ( void  ) const [inline]

Definition at line 68 of file CosmicMuonSmoother.h.

References theEstimator.

{return theEstimator;}
vector< Trajectory > CosmicMuonSmoother::fit ( const Trajectory t) const

Definition at line 96 of file CosmicMuonSmoother.cc.

References category_, Trajectory::empty(), initialState(), TrajectoryStateOnSurface::isValid(), LogTrace, Trajectory::recHits(), Trajectory::seed(), and sortHitsAlongMom().

Referenced by CosmicMuonTrajectoryBuilder::flipTrajectory(), GlobalCosmicMuonTrajectoryBuilder::trajectories(), trajectories(), and trajectory().

                                                                    {

  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);

}
vector< Trajectory > CosmicMuonSmoother::fit ( const TrajectorySeed seed,
const ConstRecHitContainer hits,
const TrajectoryStateOnSurface firstPredTsos 
) const

Definition at line 128 of file CosmicMuonSmoother.cc.

References alongMomentum, category_, Chi2MeasurementEstimator::estimate(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), LogTrace, Trajectory::measurements(), Propagator::propagate(), propagatorAlong(), Trajectory::push(), edm::second(), CosmicMuonUtilities::stepPropagate(), TrajectoryStateOnSurface::surface(), theEstimator, theService, theUpdator, theUtilities, and KFUpdator::update().

                                                                                               {

  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() ) {

    TransientTrackingRecHit::RecHitPointer preciseHit = hits.front()->clone(predTsos);
    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())  {
     LogTrace(category_)<<"current pos "<<currTsos.globalPosition()
                       <<"mom "<<currTsos.globalMomentum();
    } else {
      LogTrace(category_)<<"current state invalid";
    }

    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() ) {
      // update
      TransientTrackingRecHit::RecHitPointer preciseHit = (**ihit).clone(predTsos);

      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);

}
TrajectoryStateOnSurface CosmicMuonSmoother::initialState ( const Trajectory t) const [private]

Definition at line 423 of file CosmicMuonSmoother.cc.

References Trajectory::empty(), PV3DBase< T, PVType, FrameType >::eta(), Trajectory::firstMeasurement(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), Trajectory::lastMeasurement(), query::result, CosmicMuonUtilities::reverseDirection(), theService, theUtilities, TrajectoryMeasurement::updatedState(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by fit().

                                                                                   {
  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;

}
const Propagator* CosmicMuonSmoother::propagatorAlong ( ) const [inline]
const Propagator* CosmicMuonSmoother::propagatorOpposite ( ) const [inline]
vector< Trajectory > CosmicMuonSmoother::smooth ( const Trajectory t) const [private]

Definition at line 265 of file CosmicMuonSmoother.cc.

References category_, heavyFlavorValidationHarvestingSequence_cff::combiner, Trajectory::empty(), Chi2MeasurementEstimator::estimate(), TrajectoryMeasurement::estimate(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), LogTrace, Trajectory::measurements(), oppositeToMomentum, Propagator::propagate(), propagatorOpposite(), TrajectoryStateOnSurface::rescaleError(), Trajectory::seed(), CosmicMuonUtilities::stepPropagate(), TrajectoryStateOnSurface::surface(), theErrorRescaling, theEstimator, theUpdator, theUtilities, and KFUpdator::update().

                                                                       {

  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>();
  } 

}
std::vector<Trajectory> CosmicMuonSmoother::smooth ( const std::vector< Trajectory > &  ) const [private]

Referenced by trajectories(), and trajectory().

void CosmicMuonSmoother::sortHitsAlongMom ( ConstRecHitContainer hits,
const TrajectoryStateOnSurface tsos 
) const [private]

Definition at line 462 of file CosmicMuonSmoother.cc.

References TrajectoryStateOnSurface::globalPosition(), and mag().

Referenced by fit().

                                                                                                                {

    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;
}
virtual TrajectoryContainer CosmicMuonSmoother::trajectories ( const Trajectory traj) const [inline, virtual]
std::vector< Trajectory > CosmicMuonSmoother::trajectories ( const TrajectorySeed seed,
const ConstRecHitContainer hits,
const TrajectoryStateOnSurface firstPredTsos 
) const [virtual]

refit trajectory

Definition at line 71 of file CosmicMuonSmoother.cc.

References category_, fit(), TrajectoryStateOnSurface::isValid(), LogTrace, TrajectoryStateOnSurface::rescaleError(), smooth(), and theErrorRescaling.

                                                                                                        {

  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;

}
Trajectory CosmicMuonSmoother::trajectory ( const Trajectory t) const [virtual]

Implements TrajectorySmoother.

Definition at line 61 of file CosmicMuonSmoother.cc.

References fit(), and smooth().

                                                                   {
   std::vector<Trajectory> && fitted = fit(t);
   if (fitted.empty()) return  Trajectory();
   std::vector<Trajectory> && smoothed = smooth(fitted);
   return  smoothed.empty() ? Trajectory() : smoothed.front(); 
}
KFUpdator* CosmicMuonSmoother::updator ( ) const [inline]

Definition at line 64 of file CosmicMuonSmoother.h.

References theUpdator.

{return theUpdator;}
CosmicMuonUtilities* CosmicMuonSmoother::utilities ( ) const [inline]

Member Data Documentation

std::string CosmicMuonSmoother::category_ [private]

Definition at line 92 of file CosmicMuonSmoother.h.

Referenced by CosmicMuonSmoother(), fit(), smooth(), and trajectories().

Definition at line 91 of file CosmicMuonSmoother.h.

Referenced by CosmicMuonSmoother(), smooth(), and trajectories().

Definition at line 84 of file CosmicMuonSmoother.h.

Referenced by CosmicMuonSmoother(), estimator(), fit(), smooth(), and ~CosmicMuonSmoother().

Definition at line 89 of file CosmicMuonSmoother.h.

Referenced by CosmicMuonSmoother(), and propagatorAlong().

Definition at line 90 of file CosmicMuonSmoother.h.

Referenced by CosmicMuonSmoother(), and propagatorOpposite().

Definition at line 87 of file CosmicMuonSmoother.h.

Referenced by fit(), initialState(), propagatorAlong(), and propagatorOpposite().

Definition at line 83 of file CosmicMuonSmoother.h.

Referenced by CosmicMuonSmoother(), fit(), smooth(), updator(), and ~CosmicMuonSmoother().