ME0SegmentMatcher.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h" 
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/MuonReco/interface/ME0Muon.h"
#include "DataFormats/MuonDetId/interface/ME0DetId.h"
#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/GEMRecHit/interface/ME0SegmentCollection.h"

#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Geometry/GEMGeometry/interface/ME0Geometry.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"

class ME0SegmentMatcher : public edm::stream::EDProducer<> {
public:
  explicit ME0SegmentMatcher(const edm::ParameterSet&);
  ~ME0SegmentMatcher() override;
  void produce(edm::Event&, const edm::EventSetup&) override;

  void beginRun(edm::Run const&, edm::EventSetup const&) override;

  FreeTrajectoryState getFTS(const GlobalVector& , const GlobalVector& , 
                 int , const AlgebraicSymMatrix66& ,
                 const MagneticField* );

  FreeTrajectoryState getFTS(const GlobalVector& , const GlobalVector& , 
                 int , const AlgebraicSymMatrix55& ,
                 const MagneticField* );

  void getFromFTS(const FreeTrajectoryState& ,
          GlobalVector& , GlobalVector& , 
          int& , AlgebraicSymMatrix66& );

private:

  double theX_RESIDUAL_CUT, theX_PULL_CUT, theY_RESIDUAL_CUT, theY_PULL_CUT, thePHIDIR_RESIDUAL_CUT;
  edm::InputTag OurSegmentsTag, generalTracksTag;
  edm::EDGetTokenT<ME0SegmentCollection> OurSegmentsToken_;
  edm::EDGetTokenT<reco::TrackCollection> generalTracksToken_;
  
};

ME0SegmentMatcher::ME0SegmentMatcher(const edm::ParameterSet& pas) {
  produces<std::vector<reco::ME0Muon> >();  
  theX_PULL_CUT   = pas.getParameter<double>("maxPullX");
  theX_RESIDUAL_CUT   = pas.getParameter<double>("maxDiffX");
  theY_PULL_CUT   = pas.getParameter<double>("maxPullY");
  theY_RESIDUAL_CUT   = pas.getParameter<double>("maxDiffY");
  thePHIDIR_RESIDUAL_CUT   = pas.getParameter<double>("maxDiffPhiDirection");
  //Might need to replace "OurSegments" with an edm::InputTag of "OurSegments"
  OurSegmentsTag = pas.getParameter<edm::InputTag>("me0SegmentTag");
  generalTracksTag = pas.getParameter<edm::InputTag>("tracksTag");
  OurSegmentsToken_ = consumes<ME0SegmentCollection>(OurSegmentsTag);
  generalTracksToken_ = consumes<reco::TrackCollection>(generalTracksTag);

}

ME0SegmentMatcher::~ME0SegmentMatcher() {}

void ME0SegmentMatcher::produce(edm::Event& ev, const edm::EventSetup& setup) {

  //Getting the objects we'll need
  using namespace edm;

  ESHandle<ME0Geometry> me0Geom;
  setup.get<MuonGeometryRecord>().get(me0Geom);
  ESHandle<MagneticField> bField;
  setup.get<IdealMagneticFieldRecord>().get(bField);
  ESHandle<Propagator> ThisshProp;
  setup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAlong", ThisshProp);

  using namespace reco;

  Handle<ME0SegmentCollection> OurSegments;
  ev.getByToken(OurSegmentsToken_,OurSegments);

  auto oc = std::make_unique<std::vector<ME0Muon>>(); 
  std::vector<ME0Muon> TempStore; 

  Handle <TrackCollection > generalTracks;
  ev.getByToken(generalTracksToken_,generalTracks);

  int TrackNumber = 0;
    
  for (std::vector<Track>::const_iterator thisTrack = generalTracks->begin();
       thisTrack != generalTracks->end(); ++thisTrack,++TrackNumber){
    //Initializing our plane

    //Remove later
    if (std::abs(thisTrack->eta()) < 1.8) continue;

    //Getting the initial variables for propagation
    float zSign = thisTrack->pz() > 0 ? 1.0f : -1.0f;

    ReferenceCountingPointer<Plane> plane = Plane::build(Surface::PositionType(0,0,zSign*526.75),Surface::RotationType());

    int chargeReco = thisTrack->charge(); 
    GlobalVector p3reco, r3reco;

    p3reco = GlobalVector(thisTrack->outerPx(), thisTrack->outerPy(), thisTrack->outerPz());
    r3reco = GlobalVector(thisTrack->outerX(), thisTrack->outerY(), thisTrack->outerZ());

    AlgebraicSymMatrix66 covReco;
    //This is to fill the cov matrix correctly
    const AlgebraicSymMatrix55 & covReco_curv = thisTrack->outerStateCovariance();
    FreeTrajectoryState initrecostate = getFTS(p3reco, r3reco, chargeReco, covReco_curv, &*bField);
    getFromFTS(initrecostate, p3reco, r3reco, chargeReco, covReco);

    //Now we propagate and get the propagated variables from the propagated state
    TrajectoryStateOnSurface lastrecostate = ThisshProp->propagate(initrecostate,*plane);
    if (!lastrecostate.isValid()) continue;
    
    FreeTrajectoryState finalrecostate(*lastrecostate.freeTrajectoryState());

    AlgebraicSymMatrix66 covFinalReco;
    GlobalVector p3FinalReco_glob, r3FinalReco_globv;
    getFromFTS(finalrecostate, p3FinalReco_glob, r3FinalReco_globv, chargeReco, covFinalReco);

    //To transform the global propagated track to local coordinates
    int SegmentNumber = 0;

    reco::ME0Muon MuonCandidate;
    double ClosestDelR2 = 999.;

    for (auto thisSegment = OurSegments->begin(); thisSegment != OurSegments->end(); 
     ++thisSegment,++SegmentNumber){
      ME0DetId id = thisSegment->me0DetId();

      auto chamber = me0Geom->chamber(id);

      if ( zSign * chamber->toGlobal(thisSegment->localPosition()).z() < 0 ) continue;

      GlobalPoint r3FinalReco_glob(r3FinalReco_globv.x(),r3FinalReco_globv.y(),r3FinalReco_globv.z());

      LocalPoint r3FinalReco = chamber->toLocal(r3FinalReco_glob);
      LocalVector p3FinalReco=chamber->toLocal(p3FinalReco_glob);

      LocalPoint thisPosition(thisSegment->localPosition());
      LocalVector thisDirection(thisSegment->localDirection().x(),thisSegment->localDirection().y(),thisSegment->localDirection().z());  //FIXME

      //The same goes for the error
      AlgebraicMatrix thisCov(4,4,0);   
      for (int i = 1; i <=4; i++){
    for (int j = 1; j <=4; j++){
      thisCov(i,j) = thisSegment->parametersError()(i,j);
    }
      }

      LocalTrajectoryParameters ltp(r3FinalReco,p3FinalReco,chargeReco);
      JacobianCartesianToLocal jctl(chamber->surface(),ltp);
      const AlgebraicMatrix56& jacobGlbToLoc = jctl.jacobian(); 

      AlgebraicMatrix55 Ctmp =  (jacobGlbToLoc * covFinalReco) * ROOT::Math::Transpose(jacobGlbToLoc); 
      AlgebraicSymMatrix55 C;  // I couldn't find any other way, so I resort to the brute force
      for(int i=0; i<5; ++i) {
    for(int j=0; j<5; ++j) {
      C[i][j] = Ctmp[i][j]; 

    }
      }  

      Double_t sigmax = sqrt(C[3][3]+thisSegment->localPositionError().xx() );      
      Double_t sigmay = sqrt(C[4][4]+thisSegment->localPositionError().yy() );

      bool X_MatchFound = false, Y_MatchFound = false, Dir_MatchFound = false;
    
      if ( (std::abs(thisPosition.x()-r3FinalReco.x()) < (theX_PULL_CUT * sigmax)) || (std::abs(thisPosition.x()-r3FinalReco.x()) < theX_RESIDUAL_CUT ) ) X_MatchFound = true;
      if ( (std::abs(thisPosition.y()-r3FinalReco.y()) < (theY_PULL_CUT * sigmay)) || (std::abs(thisPosition.y()-r3FinalReco.y()) < theY_RESIDUAL_CUT ) ) Y_MatchFound = true;

      if ( std::abs(reco::deltaPhi(p3FinalReco_glob.barePhi(),chamber->toGlobal(thisSegment->localDirection()).barePhi())) < thePHIDIR_RESIDUAL_CUT) Dir_MatchFound = true;

      //Check for a Match, and if there is a match, check the delR from the segment, keeping only the closest in MuonCandidate
      if (X_MatchFound && Y_MatchFound && Dir_MatchFound) {
       
    TrackRef thisTrackRef(generalTracks,TrackNumber);
       
    GlobalPoint SegPos(chamber->toGlobal(thisSegment->localPosition()));
    GlobalPoint TkPos(r3FinalReco_globv.x(),r3FinalReco_globv.y(),r3FinalReco_globv.z());
       
    double thisDelR2 = reco::deltaR2(SegPos,TkPos);
    if (thisDelR2 < ClosestDelR2){
      ClosestDelR2 = thisDelR2;
      MuonCandidate = reco::ME0Muon(thisTrackRef,(*thisSegment),SegmentNumber,chargeReco);

      MuonCandidate.setGlobalTrackPosAtSurface(r3FinalReco_glob);
      MuonCandidate.setGlobalTrackMomAtSurface(p3FinalReco_glob);
      MuonCandidate.setLocalTrackPosAtSurface(r3FinalReco);
      MuonCandidate.setLocalTrackMomAtSurface(p3FinalReco);
      MuonCandidate.setGlobalTrackCov(covFinalReco);
      MuonCandidate.setLocalTrackCov(C);
    }
      }
    }//End loop for (auto thisSegment = OurSegments->begin(); thisSegment != OurSegments->end(); ++thisSegment,++SegmentNumber)

    //As long as the delR of the MuonCandidate is sensible, store the track-segment pair
    if (ClosestDelR2 < 500.) {
      oc->push_back(MuonCandidate);
    }
  }

  // put collection in event
  ev.put(std::move(oc));        
}

FreeTrajectoryState
ME0SegmentMatcher::getFTS(const GlobalVector& p3, const GlobalVector& r3, 
              int charge, const AlgebraicSymMatrix55& cov,
              const MagneticField* field){

  GlobalVector p3GV(p3.x(), p3.y(), p3.z());
  GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
  GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);

  CurvilinearTrajectoryError tCov(cov);
  
  return cov.kRows == 5 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars) ;
}

FreeTrajectoryState
ME0SegmentMatcher::getFTS(const GlobalVector& p3, const GlobalVector& r3, 
              int charge, const AlgebraicSymMatrix66& cov,
              const MagneticField* field){

  GlobalVector p3GV(p3.x(), p3.y(), p3.z());
  GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
  GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);

  CartesianTrajectoryError tCov(cov);
  
  return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars) ;
}

void ME0SegmentMatcher::getFromFTS(const FreeTrajectoryState& fts,
                   GlobalVector& p3, GlobalVector& r3, 
                   int& charge, AlgebraicSymMatrix66& cov){
  GlobalVector p3GV = fts.momentum();
  GlobalPoint r3GP = fts.position();

  GlobalVector p3T(p3GV.x(), p3GV.y(), p3GV.z());
  GlobalVector r3T(r3GP.x(), r3GP.y(), r3GP.z());
  p3 = p3T;
  r3 = r3T;  
  // p3.set(p3GV.x(), p3GV.y(), p3GV.z());
  // r3.set(r3GP.x(), r3GP.y(), r3GP.z());
  
  charge = fts.charge();
  cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
}

void ME0SegmentMatcher::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup){}

DEFINE_FWK_MODULE(ME0SegmentMatcher);