Me0 global reco selectors

DataFormats/MuonReco/interface/ME0Muon.h

@@ -8,6 +8,9 @@
  */
 #include "DataFormats/TrackReco/interface/TrackFwd.h"
 #include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "DataFormats/GeometryVector/interface/GlobalVector.h"
+#include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
 
 #include <DataFormats/GEMRecHit/interface/ME0SegmentCollection.h>
 
@@ -53,12 +56,42 @@ namespace reco {
     double phi() const { return innerTrack_.get()->phi(); }
     /// pseudorapidity of momentum vector
     double eta() const { return innerTrack_.get()->eta(); }
+
+    //functions for easy variation of me0muon criteria
+    /* double xpull() const { return xpull_; } */
+    /* double xdiff() const { return xdiff_; } */
+    /* double ypull() const { return ypull_; } */
+    /* double ydiff() const { return ydiff_; } */
+    /* double phidirdiff() const { return phidirdiff_; } */
+
+    /* void setXpull( const double xpull ) { xpull_ = xpull; } */
+    /* void setXdiff( const double xdiff ) { xdiff_ = xdiff; } */
+    /* void setYpull( const double ypull ) { ypull_ = ypull; } */
+    /* void setYdiff( const double ydiff ) { ydiff_ = ydiff; } */
+    /* void setPhidirdiff( const double phidirdiff ) { phidirdiff_ = phidirdiff; } */
+
+    GlobalPoint globalTrackPosAtSurface() const { return globalTrackPosAtSurface_; }
+    GlobalVector globalTrackMomAtSurface() const { return globalTrackMomAtSurface_; }
+    int trackCharge() const { return trackCharge_; }
+    AlgebraicSymMatrix66 trackCov() const { return trackCov_; }
+
+    void setGlobalTrackPosAtSurface(const GlobalPoint globalTrackPosAtSurface) { globalTrackPosAtSurface_ = globalTrackPosAtSurface; }
+    void setGlobalTrackMomAtSurface(const GlobalVector globalTrackMomAtSurface) { globalTrackMomAtSurface_ = globalTrackMomAtSurface; }
+    void setTrackCharge(const int trackCharge) { trackCharge_ = trackCharge; }
+    void setTrackCov(const AlgebraicSymMatrix66 trackCov) { trackCov_ = trackCov; }
 
   private:
     /// reference to Track reconstructed in the tracker only
     TrackRef innerTrack_;
     ME0Segment me0Segment_;
     int me0segid_;
+
+    GlobalPoint globalTrackPosAtSurface_;
+    GlobalVector globalTrackMomAtSurface_;
+    int trackCharge_;
+    AlgebraicSymMatrix66 trackCov_;
+
+    //double xpull_,ypull_,xdiff_,ydiff_,phidirdiff_;
   };
 
 }

DataFormats/MuonReco/src/classes_def.xml

@@ -184,8 +184,8 @@
     <class name="edm::Ref<std::vector<EmulatedME0Segment>,EmulatedME0Segment,edm::refhelper::FindUsingAdvance<std::vector<EmulatedME0Segment>,EmulatedME0Segment> >"/>
 
     <class name="std::vector<reco::ME0Muon>"/>
-    <class name="reco::ME0Muon" ClassVersion="15">
-      <version ClassVersion="15" checksum="3167678990"/>
+    <class name="reco::ME0Muon" ClassVersion="17">
+      <version ClassVersion="17" checksum="3425786499"/>
     </class>
     <class name="edm::Wrapper<std::vector<reco::ME0Muon> >"/>
 

RecoMuon/MuonIdentification/interface/ME0MuonSelector.h

@@ -0,0 +1,42 @@
+#ifndef RecoMuon_ME0MuonSelectors_h
+#define RecoMuon_ME0MuonSelectors_h
+//
+
+#include "DataFormats/MuonReco/interface/ME0Muon.h"
+
+#include "Geometry/GEMGeometry/interface/ME0Geometry.h"
+#include <Geometry/GEMGeometry/interface/ME0EtaPartition.h>
+#include <Geometry/Records/interface/MuonGeometryRecord.h>
+#include <DataFormats/MuonDetId/interface/ME0DetId.h>
+#include "FWCore/Framework/interface/ESHandle.h"
+
+#include "TMath.h"
+#include <string>
+
+
+
+
+namespace muon {
+   /// Selector type
+   enum SelectionType {
+      All = 0,                      // dummy options - always true
+      VeryLoose = 1,           //
+      Loose = 2,           //
+      Tight = 3,           //
+   };
+
+   /// a lightweight "map" for selection type string label and enum value
+   struct SelectionTypeStringToEnum { const char *label; SelectionType value; };
+   SelectionType selectionTypeFromString( const std::string &label );
+
+   /// main GoodMuon wrapper call
+   bool isGoodMuon( edm::ESHandle <ME0Geometry> me0Geom, const reco::ME0Muon& me0muon, SelectionType type );
+
+
+   /// Specialized isGoodMuon function called from main wrapper
+
+   bool isGoodMuon( edm::ESHandle <ME0Geometry> me0Geom, const reco::ME0Muon& me0muon, double MaxPullX, double MaxDiffX, double MaxPullY, double MaxDiffY, double MaxDiffPhiDir  );
+
+
+}
+#endif

RecoMuon/MuonIdentification/plugins/ME0MuonSelector.cc

@@ -0,0 +1,121 @@
+/** 
+ *Filter to select me0Muons based on pulls and differences w.r.t. me0Segments
+ *
+ *
+ */
+
+#include "RecoMuon/MuonIdentification/interface/ME0MuonSelector.h"
+
+#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
+#include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
+#include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCartesian.h"
+
+
+namespace muon {
+SelectionType selectionTypeFromString( const std::string &label )
+{
+   static SelectionTypeStringToEnum selectionTypeStringToEnumMap[] = {
+      { "All", All },
+      { "VeryLoose", VeryLoose },
+      { "Loose", Loose },
+      { "Tight", Tight },
+      { 0, (SelectionType)-1 }
+   };
+
+   SelectionType value = (SelectionType)-1;
+   bool found = false;
+   for(int i = 0; selectionTypeStringToEnumMap[i].label && (! found); ++i)
+      if (! strcmp(label.c_str(), selectionTypeStringToEnumMap[i].label)) {
+         found = true;
+         value = selectionTypeStringToEnumMap[i].value;
+      }
+
+   // in case of unrecognized selection type
+   if (! found) throw cms::Exception("MuonSelectorError") << label << " is not a recognized SelectionType";
+   return value;
+}
+}
+
+
+bool muon::isGoodMuon(edm::ESHandle<ME0Geometry> me0Geom, const reco::ME0Muon& me0muon, SelectionType type)
+{
+  switch (type)
+    {
+    case muon::All:
+      return true;
+      break;
+    case muon::VeryLoose:
+      return isGoodMuon(me0Geom, me0muon,3,4,20,20,3.14); 
+      break;
+    case muon::Loose:
+      return isGoodMuon(me0Geom, me0muon,3,2,3,2,0.5);
+      break;
+    case muon::Tight:
+      return isGoodMuon(me0Geom, me0muon,3,2,3,2,0.15);
+      break;
+    default:
+      return false;
+    }
+}
+
+
+
+
+bool muon::isGoodMuon(edm::ESHandle<ME0Geometry> me0Geom, const reco::ME0Muon& me0muon, double MaxPullX, double MaxDiffX, double MaxPullY, double MaxDiffY, double MaxDiffPhiDir )
+{
+  ME0Segment thisSegment = me0muon.me0segment();
+
+  ME0DetId id = thisSegment.me0DetId();
+
+  auto roll = me0Geom->etaPartition(id); 
+
+  float zSign  = me0muon.globalTrackMomAtSurface().z()/fabs(me0muon.globalTrackMomAtSurface().z());
+  if ( zSign * roll->toGlobal(thisSegment.localPosition()).z() < 0 ) return false;
+
+  //GlobalPoint r3FinalReco_glob(r3FinalReco_globv.x(),r3FinalReco_globv.y(),r3FinalReco_globv.z());
+
+  LocalPoint r3FinalReco = roll->toLocal(me0muon.globalTrackPosAtSurface());
+  LocalVector p3FinalReco=roll->toLocal(me0muon.globalTrackMomAtSurface());
+
+  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,me0muon.trackCharge());
+  JacobianCartesianToLocal jctl(roll->surface(),ltp);
+  AlgebraicMatrix56 jacobGlbToLoc = jctl.jacobian(); 
+
+  AlgebraicMatrix55 Ctmp =  (jacobGlbToLoc * me0muon.trackCov()) * 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 ( ( (fabs(thisPosition.x()-r3FinalReco.x())/sigmax ) < MaxPullX ) || (fabs(thisPosition.x()-r3FinalReco.x()) < MaxDiffX ) ) X_MatchFound = true;
+  if ( ( (fabs(thisPosition.y()-r3FinalReco.y())/sigmay ) < MaxPullY ) || (fabs(thisPosition.y()-r3FinalReco.y()) < MaxDiffY ) ) Y_MatchFound = true;
+
+  if ( fabs(me0muon.globalTrackMomAtSurface().phi()-roll->toGlobal(thisSegment.localDirection()).phi()) < MaxDiffPhiDir) Dir_MatchFound = true;
+
+  return (X_MatchFound && Y_MatchFound && Dir_MatchFound);
+
+}
+

RecoMuon/MuonIdentification/plugins/ME0SegmentMatcher.cc

@@ -39,6 +39,11 @@
 
 ME0SegmentMatcher::ME0SegmentMatcher(const edm::ParameterSet& pas) : iev(0){
   produces<std::vector<reco::ME0Muon> >();  
+  X_PULL_CUT   = pas.getParameter<double>("maxPullX");
+  X_RESIDUAL_CUT   = pas.getParameter<double>("maxDiffX");
+  Y_PULL_CUT   = pas.getParameter<double>("maxPullY");
+  Y_RESIDUAL_CUT   = pas.getParameter<double>("maxDiffY");
+  PHIDIR_RESIDUAL_CUT   = pas.getParameter<double>("maxDiffPhiDirection");
 }
 
 ME0SegmentMatcher::~ME0SegmentMatcher() {}
@@ -167,10 +172,16 @@ void ME0SegmentMatcher::produce(edm::Event& ev, const edm::EventSetup& setup) {
 	bool X_MatchFound = false, Y_MatchFound = false, Dir_MatchFound = false;
 
 
-	 if ( (fabs(thisPosition.x()-r3FinalReco.x()) < (3.0 * sigmax)) || (fabs(thisPosition.x()-r3FinalReco.x()) < 2.0 ) ) X_MatchFound = true;
-	 if ( (fabs(thisPosition.y()-r3FinalReco.y()) < (3.0 * sigmay)) || (fabs(thisPosition.y()-r3FinalReco.y()) < 2.0 ) ) Y_MatchFound = true;
+	 // if ( (fabs(thisPosition.x()-r3FinalReco.x()) < (3.0 * sigmax)) || (fabs(thisPosition.x()-r3FinalReco.x()) < 2.0 ) ) X_MatchFound = true;
+	 // if ( (fabs(thisPosition.y()-r3FinalReco.y()) < (3.0 * sigmay)) || (fabs(thisPosition.y()-r3FinalReco.y()) < 2.0 ) ) Y_MatchFound = true;
 
-	 if ( fabs(p3FinalReco_glob.phi()-roll->toGlobal(thisSegment->localDirection()).phi()) < 0.15) Dir_MatchFound = true;
+	 // if ( fabs(p3FinalReco_glob.phi()-roll->toGlobal(thisSegment->localDirection()).phi()) < 0.15) Dir_MatchFound = true;
+
+
+	 if ( (fabs(thisPosition.x()-r3FinalReco.x()) < (X_PULL_CUT * sigmax)) || (fabs(thisPosition.x()-r3FinalReco.x()) < X_RESIDUAL_CUT ) ) X_MatchFound = true;
+	 if ( (fabs(thisPosition.y()-r3FinalReco.y()) < (Y_PULL_CUT * sigmay)) || (fabs(thisPosition.y()-r3FinalReco.y()) < Y_RESIDUAL_CUT ) ) Y_MatchFound = true;
+
+	 if ( fabs(p3FinalReco_glob.phi()-roll->toGlobal(thisSegment->localDirection()).phi()) < PHIDIR_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) {
@@ -184,6 +195,20 @@ void ME0SegmentMatcher::produce(edm::Event& ev, const edm::EventSetup& setup) {
 	   if (thisDelR < ClosestDelR){
 	     ClosestDelR = thisDelR;
 	     MuonCandidate = reco::ME0Muon(thisTrackRef,(*thisSegment),SegmentNumber);
+
+	     //Setting the variables for easy me0muon selection later on
+	     // MuonCandidate.setXpull( ( fabs(thisPosition.x()-r3FinalReco.x()) / sigmax) );
+	     // MuonCandidate.setYpull( ( fabs(thisPosition.y()-r3FinalReco.y()) / sigmay) );
+
+	     // MuonCandidate.setXpull( fabs(thisPosition.x()-r3FinalReco.x()) );
+	     // MuonCandidate.setYpull( fabs(thisPosition.y()-r3FinalReco.y()) );
+
+	     // MuonCandidate.setPhidirdiff(fabs(p3FinalReco_glob.phi()-roll->toGlobal(thisSegment->localDirection()).phi()));
+
+	     MuonCandidate.setGlobalTrackPosAtSurface(r3FinalReco_glob);
+	     MuonCandidate.setGlobalTrackMomAtSurface(p3FinalReco_glob);
+	     MuonCandidate.setTrackCharge(chargeReco);
+	     MuonCandidate.setTrackCov(covFinalReco);
 	   }
 	 }
       }//End loop for (auto thisSegment = OurSegments->begin(); thisSegment != OurSegments->end(); ++thisSegment,++SegmentNumber)

RecoMuon/MuonIdentification/plugins/ME0SegmentMatcher.h

@@ -70,7 +70,7 @@ class ME0SegmentMatcher : public edm::EDProducer {
     int iev; // events through
 
     edm::ESHandle<ME0Geometry> me0Geom;
-
+    double X_PULL_CUT, Y_PULL_CUT,X_RESIDUAL_CUT,Y_RESIDUAL_CUT, PHIDIR_RESIDUAL_CUT;
 
 };
 

RecoMuon/MuonIdentification/python/me0SegmentMatcher_cfi.py

@@ -1,4 +1,10 @@
 import FWCore.ParameterSet.Config as cms
 
 
-me0SegmentMatching = cms.EDProducer("ME0SegmentMatcher")
+me0SegmentMatching = cms.EDProducer("ME0SegmentMatcher",
+                                    maxPullX = cms.double (3.0),
+                                    maxDiffX = cms.double (4.0),
+                                    maxPullY = cms.double (20.0),
+                                    maxDiffY = cms.double (20.0),
+                                    maxDiffPhiDirection = cms.double (3.14),
+                                    )