CRackTrajectoryBuilder.h

Go to the documentation of this file.

#ifndef CRackTrajectoryBuilder_h
#define CRackTrajectoryBuilder_h

//
// Package:         RecoTracker/SingleTrackPattern
// Class:           CRackTrajectoryBuilder
// Original Author:  Michele Pioppi-INFN perugia

#include <string>

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/EventSetup.h"

#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
#include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
#include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
#include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"

#include "FWCore/Framework/interface/ESHandle.h"

//to sort hits by the det position
class CompareDetY_plus {
public:
  CompareDetY_plus(const TrackerGeometry& tracker) : _tracker(tracker) {}
  bool operator()(const TrackingRecHit* rh1, const TrackingRecHit* rh2) {
    const GeomDet* detPos1 = _tracker.idToDet(rh1->geographicalId());
    const GeomDet* detPos2 = _tracker.idToDet(rh2->geographicalId());

    const GlobalPoint& gp1 = detPos1->position();
    const GlobalPoint& gp2 = detPos2->position();

    if (gp1.y() > gp2.y())
      return true;
    if (gp1.y() < gp2.y())
      return false;
    //  if (gp1.y()== gp2.y())
    //
    return (rh1->geographicalId() < rh2->geographicalId());
  };

private:
  //   edm::ESHandle<TrackerGeometry> _tracker;
  const TrackerGeometry& _tracker;
};

class CompareDetY_minus {
public:
  CompareDetY_minus(const TrackerGeometry& tracker) : _tracker(tracker) {}
  bool operator()(const TrackingRecHit* rh1, const TrackingRecHit* rh2) {
    const GeomDet* detPos1 = _tracker.idToDet(rh1->geographicalId());
    const GeomDet* detPos2 = _tracker.idToDet(rh2->geographicalId());

    const GlobalPoint& gp1 = detPos1->position();
    const GlobalPoint& gp2 = detPos2->position();

    if (gp1.y() < gp2.y())
      return true;
    if (gp1.y() > gp2.y())
      return false;
    //  if (gp1.y()== gp2.y())
    //
    return (rh1->geographicalId() < rh2->geographicalId());
  };

private:
  //   edm::ESHandle<TrackerGeometry> _tracker;
  const TrackerGeometry& _tracker;
};

#ifndef TrajectoryBuilder_CompareHitY
#define TrajectoryBuilder_CompareHitY

class CompareHitY {
public:
  CompareHitY(const TrackerGeometry& tracker) : _tracker(tracker) {}
  bool operator()(const TrackingRecHit* rh1, const TrackingRecHit* rh2) {
    GlobalPoint gp1 = _tracker.idToDet(rh1->geographicalId())->surface().toGlobal(rh1->localPosition());
    GlobalPoint gp2 = _tracker.idToDet(rh2->geographicalId())->surface().toGlobal(rh2->localPosition());
    return gp1.y() < gp2.y();
  };

private:
  //   edm::ESHandle<TrackerGeometry> _tracker;
  const TrackerGeometry& _tracker;
};

class CompareHitY_plus {
public:
  CompareHitY_plus(const TrackerGeometry& tracker) : _tracker(tracker) {}
  bool operator()(const TrackingRecHit* rh1, const TrackingRecHit* rh2) {
    GlobalPoint gp1 = _tracker.idToDet(rh1->geographicalId())->surface().toGlobal(rh1->localPosition());
    GlobalPoint gp2 = _tracker.idToDet(rh2->geographicalId())->surface().toGlobal(rh2->localPosition());
    return gp1.y() > gp2.y();
  };

private:
  //   edm::ESHandle<TrackerGeometry> _tracker;
  const TrackerGeometry& _tracker;
};

#endif

class CRackTrajectoryBuilder {
  //  using namespace std;

  typedef TrajectoryStateOnSurface TSOS;
  typedef TrajectoryMeasurement TM;

  typedef std::vector<const TrackingRecHit*>::iterator TrackingRecHitIterator;

  typedef std::pair<TrackingRecHitIterator, TrackingRecHitIterator> TrackingRecHitRange;
  typedef std::vector<TrackingRecHitRange>::iterator TrackingRecHitRangeIterator;

  //  typedef std::pair<TrackingRecHitIterator, TSOS> PairTrackingRecHitTsos;
  typedef std::pair<TrackingRecHitRangeIterator, TSOS> PairTrackingRecHitTsos;

public:
  class CompareDetByTraj;
  friend class CompareDetByTraj;

  class CompareDetByTraj {
  public:
    CompareDetByTraj(const TSOS& tSos) : _tSos(tSos){};
    bool operator()(const std::pair<TrackingRecHitRangeIterator, TSOS> rh1,
                    const std::pair<TrackingRecHitRangeIterator, TSOS> rh2) {
      GlobalPoint gp1 = rh1.second.globalPosition();
      GlobalPoint gp2 = rh2.second.globalPosition();

      GlobalPoint gpT = _tSos.globalPosition();
      GlobalVector gpDiff1 = gp1 - gpT;
      GlobalVector gpDiff2 = gp2 - gpT;

      //this might have a better performance ...
      //       float dist1 = ( gp1.x()-gpT.x() ) * ( gp1.x()-gpT.x() ) + ( gp1.y()-gpT.y() ) * ( gp1.y()-gpT.y() ) + ( gp1.z()-gpT.z() ) * ( gp1.z()-gpT.z() );
      //       float dist2 = ( gp2.x()-gpT.x() ) * ( gp2.x()-gpT.x() ) + ( gp2.y()-gpT.y() ) * ( gp2.y()-gpT.y() ) + ( gp2.z()-gpT.z() ) * ( gp2.z()-gpT.z() );
      //if ( dist1<dist2 )

      //     if ( gpDiff1.mag2() < gpDiff2.mag2() )

      float dist1 = gpDiff1 * _tSos.globalDirection();
      float dist2 = gpDiff2 * _tSos.globalDirection();

      if (dist1 < 0)
        return false;
      if (dist1 < dist2)
        return true;

      return false;
    };

  private:
    const TrajectoryStateOnSurface& _tSos;
  };

public:
  CRackTrajectoryBuilder(const edm::ParameterSet& conf);
  ~CRackTrajectoryBuilder();


  void run(const TrajectorySeedCollection& collseed,
           const SiStripRecHit2DCollection& collstereo,
           const SiStripRecHit2DCollection& collrphi,
           const SiStripMatchedRecHit2DCollection& collmatched,
           const SiPixelRecHitCollection& collpixel,
           const edm::EventSetup& es,
           edm::Event& e,
           std::vector<Trajectory>& trajoutput);

  void init(const edm::EventSetup& es, bool);
  Trajectory createStartingTrajectory(const TrajectorySeed& seed) const;

  const TransientTrackingRecHitBuilder* hitBuilder() const { return RHBuilder; }

private:
  std::vector<TrajectoryMeasurement> seedMeasurements(const TrajectorySeed& seed) const;

  std::vector<const TrackingRecHit*> SortHits(const SiStripRecHit2DCollection& collstereo,
                                              const SiStripRecHit2DCollection& collrphi,
                                              const SiStripMatchedRecHit2DCollection& collmatched,
                                              const SiPixelRecHitCollection& collpixel,
                                              const TrajectorySeed& seed,
                                              const bool bAddSeedHits);

  //    std::vector<TrackingRecHitRange> SortByTrajectory (const std::vector<TrackingRecHitRange>& inputHits);

  TSOS startingTSOS(const TrajectorySeed& seed) const;
  void updateTrajectory(Trajectory& traj, const TM& tm, const TransientTrackingRecHit& hit) const;

  void AddHit(Trajectory& traj, const std::vector<const TrackingRecHit*>& Hits, Propagator* currPropagator);
  //        edm::OwnVector<TransientTrackingRecHit> hits);
  bool qualityFilter(const Trajectory& traj);

  bool isDifferentStripReHit2D(const SiStripRecHit2D& hitA, const SiStripRecHit2D& hitB);

  std::pair<TrajectoryStateOnSurface, const GeomDet*> innerState(const Trajectory& traj) const;

private:
  edm::ESHandle<MagneticField> magfield;
  edm::ESHandle<TrackerGeometry> tracker;

  PropagatorWithMaterial* thePropagator;
  PropagatorWithMaterial* thePropagatorOp;

  //   AnalyticalPropagator *thePropagator;
  //   AnalyticalPropagator *thePropagatorOp;

  KFUpdator* theUpdator;
  Chi2MeasurementEstimator* theEstimator;
  const TransientTrackingRecHitBuilder* RHBuilder;
  const KFTrajectorySmoother* theSmoother;
  const KFTrajectoryFitter* theFitter;
  //   const KFTrajectoryFitter * theFitterOp;

  bool debug_info;
  bool fastPropagation;
  bool useMatchedHits;

  int theMinHits;
  double chi2cut;
  std::vector<Trajectory> trajFit;
  //RC edm::OwnVector<const TransientTrackingRecHit> hits;
  TransientTrackingRecHit::RecHitContainer hits;
  bool seed_plus;
  std::string geometry;
  std::string theBuilderName;
  //   TransientInitialStateEstimator*  theInitialState;
};

#endif