SETFilter.h

Go to the documentation of this file.

#ifndef RecoMuon_MuonSeedGenerator_SETFilter_H
#define RecoMuon_MuonSeedGenerator_SETFilter_H

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/DetLayers/interface/NavigationDirection.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"

#include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"

//#include "CLHEP/Matrix/DiagMatrix.h"
#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Vector/ThreeVector.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "DataFormats/Common/interface/OwnVector.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"

#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"

#include "TrackingTools/PatternTools/interface/Trajectory.h"

#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
#include "FWCore/Framework/interface/ESHandle.h"


// used in the SET algorithm
struct SeedCandidate{
  MuonTransientTrackingRecHit::MuonRecHitContainer theSet;
  CLHEP::Hep3Vector momentum;
  int charge;
  double weight;
  Trajectory::DataContainer trajectoryMeasurementsInTheSet;
};

class DetLayer;
class Trajectory;
//class MuonServiceProxy;
class TrajectoryFitter;

namespace edm {class ParameterSet; class EventSetup; class Event;}

class SETFilter {

 public:
  SETFilter(const edm::ParameterSet& par, const MuonServiceProxy* service);

  virtual ~SETFilter();

  // Operations
  
  void refit(const TrajectoryStateOnSurface& initialState, const DetLayer*, Trajectory& trajectory);

  FreeTrajectoryState lastUpdatedFTS() const {return *theLastUpdatedTSOS.freeTrajectoryState();}

  TrajectoryStateOnSurface lastUpdatedTSOS() const {return theLastUpdatedTSOS;}

  bool fwfit_SET( std::vector < SeedCandidate> & validSegmentsSet_in,
                           std::vector < SeedCandidate> & validSegmentsSet_out);
  
  bool buildTrajectoryMeasurements( SeedCandidate * validSegmentsSet,
                    Trajectory::DataContainer & finalCandidate);


  bool transform(Trajectory::DataContainer &measurements_segments,
                 TransientTrackingRecHit::ConstRecHitContainer & hitContainer, 
         TrajectoryStateOnSurface & firstTSOS);

  bool transformLight(Trajectory::DataContainer &measurements_segments,
              TransientTrackingRecHit::ConstRecHitContainer & hitContainer, 
              TrajectoryStateOnSurface & firstTSOS);



  void reset();

  virtual void setEvent(const edm::Event& event);

  int getTotalChamberUsed() const {return totalChambers;}
  int getDTChamberUsed() const {return dtChambers;}
  int getCSCChamberUsed() const {return cscChambers;}
  int getRPCChamberUsed() const {return rpcChambers;}

  inline bool goodState() const {return totalChambers >= 2 && 
                   ((dtChambers + cscChambers) >0 );}
  
  std::vector<const DetLayer*> layers() const {return theDetLayers;}

  const DetLayer* lastDetLayer() const {return theDetLayers.back();}

  PropagationDirection propagationDirection() const;

protected:

private:

  void setLastUpdatedTSOS(TrajectoryStateOnSurface tsos) { theLastUpdatedTSOS = tsos;}
  
  void setLastButOneUpdatedTSOS(TrajectoryStateOnSurface tsos) { theLastButOneUpdatedTSOS = tsos;}
  
  void incrementChamberCounters(const DetLayer *layer);
 
  const Propagator *propagator() const;
  

  // chi2 functions (calculate chi2)
  double findChi2(double pX, double pY, double pZ,
                    const CLHEP::Hep3Vector& r3T,
                    SeedCandidate & muonCandidate,
                    TrajectoryStateOnSurface  &lastUpdatedTSOS,
                    Trajectory::DataContainer & trajectoryMeasurementsInTheSet,
                    bool detailedOutput);

  double findMinChi2(unsigned int iSet, const CLHEP::Hep3Vector& r3T,
                 SeedCandidate & muonCandidate,
                 std::vector < TrajectoryStateOnSurface > &lastUpdatedTSOS_Vect,
                 Trajectory::DataContainer & trajectoryMeasurementsInTheSet);

  double chi2AtSpecificStep(CLHEP::Hep3Vector &foot,
                            const CLHEP::Hep3Vector& r3T,
                            SeedCandidate & muonCandidate,
                            TrajectoryStateOnSurface  &lastUpdatedTSOS,
                            Trajectory::DataContainer & trajectoryMeasurementsInTheSet,
                            bool detailedOutput);

  // find initial points for the SIMPLEX minimization
  std::vector <CLHEP::Hep3Vector> find3MoreStartingPoints(CLHEP::Hep3Vector &key_foot,
                                                   const CLHEP::Hep3Vector& r3T,
                                                   SeedCandidate & muonCandidate);

  std::pair <double,double> findParabolaMinimum(std::vector <double> &quadratic_var,
                                                std::vector <double> &quadratic_chi2);

  // SIMPLEX minimization functions
  void pickElements(std::vector <double> &chi2Feet,
                   unsigned int & high, unsigned int & second_high, unsigned int & low);

  CLHEP::Hep3Vector reflectFoot(std::vector <CLHEP::Hep3Vector> & feet,
                         unsigned int key_foot, double scale );

  void nDimContract(std::vector <CLHEP::Hep3Vector> & feet, unsigned int low);
  //---- SET
  
  std::string thePropagatorName;

  NavigationDirection theFitDirection;

  TrajectoryStateOnSurface theLastUpdatedTSOS;
  TrajectoryStateOnSurface theLastButOneUpdatedTSOS;

  std::vector<const DetLayer*> theDetLayers;

  int totalChambers;
  int dtChambers;
  int cscChambers;
  int rpcChambers;

  bool useSegmentsInTrajectory;

  edm::ESHandle<TrajectoryFitter> theBWLightFitter;
  std::string theBWLightFitterName;

  const MuonServiceProxy *theService;
  //bool theOverlappingChambersFlag;
};
#endif