| 1 |
|
#include "UserCode/L1RpcTriggerAnalysis/interface/SynchroSelector.h" |
| 2 |
+ |
|
| 3 |
|
#include "FWCore/Framework/interface/Event.h" |
| 4 |
|
#include "FWCore/Framework/interface/EventSetup.h" |
| 5 |
|
#include "FWCore/Framework/interface/ESHandle.h" |
| 7 |
|
|
| 8 |
|
#include "DataFormats/TrackReco/interface/TrackFwd.h" |
| 9 |
|
#include "DataFormats/TrackReco/interface/Track.h" |
| 9 |
– |
#include "DataFormats/TrackReco/interface/print.h" |
| 10 |
|
|
| 11 |
– |
#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h" |
| 12 |
– |
#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTCand.h" |
| 13 |
– |
#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTExtendedCand.h" |
| 14 |
– |
#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h" |
| 11 |
|
|
| 12 |
|
#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" |
| 13 |
|
#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h" |
| 19 |
|
#include "TrackingTools/GeomPropagators/interface/Propagator.h" |
| 20 |
|
|
| 21 |
|
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" |
| 26 |
– |
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" |
| 27 |
– |
#include "DataFormats/GeometrySurface/interface/BoundCylinder.h" |
| 28 |
– |
#include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h" |
| 29 |
– |
#include "DataFormats/GeometrySurface/interface/BoundDisk.h" |
| 30 |
– |
#include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h" |
| 31 |
– |
#include "DataFormats/BeamSpot/interface/BeamSpot.h" |
| 32 |
– |
|
| 22 |
|
|
| 23 |
|
#include "DataFormats/MuonDetId/interface/RPCDetId.h" |
| 24 |
|
#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h" |
| 29 |
|
using namespace std; |
| 30 |
|
|
| 31 |
|
SynchroSelector::SynchroSelector(const edm::ParameterSet & cfg) |
| 32 |
< |
: theConfig(cfg), hDxy(0),hNum(0),hDeltaEta(0),hDeltaPhi(0), |
| 44 |
< |
hPullX(0),hDistX(0) |
| 32 |
> |
: theConfig(cfg), hPullX(0),hDistX(0) |
| 33 |
|
{ } |
| 34 |
|
|
| 47 |
– |
bool SynchroSelector::checkL1RpcMatching( const TrajectoryStateOnSurface & tsos, const edm::Event&ev, const edm::EventSetup& es) |
| 48 |
– |
{ |
| 49 |
– |
edm::Handle<L1MuGMTReadoutCollection> pCollection; |
| 50 |
– |
InputTag l1MuReadout = theConfig.getParameter<edm::InputTag>("l1MuReadout"); |
| 51 |
– |
ev.getByLabel(l1MuReadout,pCollection); |
| 52 |
– |
L1MuGMTReadoutCollection const* gmtrc = pCollection.product(); |
| 53 |
– |
if (!gmtrc) return false; |
| 54 |
– |
|
| 55 |
– |
edm::ESHandle<GlobalTrackingGeometry> globalGeometry; |
| 56 |
– |
es.get<GlobalTrackingGeometryRecord>().get(globalGeometry); |
| 57 |
– |
edm::ESHandle<MagneticField> magField; |
| 58 |
– |
es.get<IdealMagneticFieldRecord>().get(magField); |
| 59 |
– |
ReferenceCountingPointer<Surface> rpc; |
| 60 |
– |
|
| 61 |
– |
vector<L1MuGMTReadoutRecord> gmt_records = gmtrc->getRecords(); |
| 62 |
– |
vector<L1MuGMTReadoutRecord>::const_iterator RRItr; |
| 63 |
– |
typedef vector<L1MuRegionalCand>::const_iterator ITC; |
| 64 |
– |
for( RRItr = gmt_records.begin() ; RRItr != gmt_records.end() ; RRItr++ ) { |
| 65 |
– |
for (int i=1; i<=2; ++i) { |
| 66 |
– |
vector<L1MuRegionalCand> Cands = (i==1) ? RRItr->getBrlRPCCands() : RRItr->getFwdRPCCands(); |
| 67 |
– |
for(ITC it = Cands.begin() ; it != Cands.end() ; ++it ) { |
| 68 |
– |
if (it->empty()) continue; |
| 69 |
– |
//propagate and check matching for candidate |
| 70 |
– |
float rpcEta = it->etaValue(); |
| 71 |
– |
float rpcPhi = it->phiValue(); |
| 72 |
– |
if (rpcEta < -1.04) rpc = ReferenceCountingPointer<Surface>(new BoundDisk( GlobalPoint(0.,0.,-800.), TkRotation<float>(), SimpleDiskBounds( 300., 710., -10., 10. ) ) ); |
| 73 |
– |
else if (rpcEta < 1.04) rpc = ReferenceCountingPointer<Surface>(new BoundCylinder( GlobalPoint(0.,0.,0.), TkRotation<float>(), SimpleCylinderBounds( 500, 520, -700, 700 ) ) ); |
| 74 |
– |
else rpc = ReferenceCountingPointer<Surface>(new BoundDisk( GlobalPoint(0.,0.,800.), TkRotation<float>(), SimpleDiskBounds( 300., 710., -10., 10. ) ) ); |
| 75 |
– |
edm::ESHandle<Propagator> propagator; |
| 76 |
– |
es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator); |
| 77 |
– |
TrajectoryStateOnSurface trackAtRPC = propagator->propagate(tsos, *rpc); |
| 78 |
– |
if (!trackAtRPC.isValid()) return false; |
| 79 |
– |
float phi = trackAtRPC.globalPosition().phi(); |
| 80 |
– |
float dphi = phi-rpcPhi; |
| 81 |
– |
if (dphi < -M_PI) dphi+=2*M_PI; |
| 82 |
– |
if (dphi > M_PI) dphi-=2*M_PI; |
| 83 |
– |
float eta = trackAtRPC.globalPosition().eta(); |
| 84 |
– |
float deta = eta-rpcEta; |
| 85 |
– |
if (fabs(dphi) < mindphi ) mindphi = fabs(dphi); |
| 86 |
– |
if (fabs(deta) < mindeta ) mindeta = fabs(deta); |
| 87 |
– |
|
| 88 |
– |
double maxDeltaEta = theConfig.getParameter<double>("maxDeltaEta"); |
| 89 |
– |
double maxDeltaPhi = theConfig.getParameter<double>("maxDeltaPhi"); |
| 90 |
– |
bool matching = ( fabs(dphi) < maxDeltaPhi && fabs(deta) < maxDeltaEta) ? true : false; |
| 91 |
– |
if (matching) return true; |
| 92 |
– |
} |
| 93 |
– |
} |
| 94 |
– |
} |
| 95 |
– |
return false; |
| 96 |
– |
} |
| 97 |
– |
|
| 98 |
– |
bool SynchroSelector::takeIt(const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es) |
| 99 |
– |
{ |
| 100 |
– |
bool result = false; |
| 101 |
– |
float theMaxTIP = theConfig.getParameter<double>("maxTIP"); |
| 102 |
– |
float theMinPt = theConfig.getParameter<double>("minPt"); |
| 103 |
– |
|
| 104 |
– |
edm::Handle<reco::TrackCollection> trackCollection; |
| 105 |
– |
ev.getByLabel(InputTag( theConfig.getParameter<std::string>("collection")),trackCollection); |
| 106 |
– |
|
| 107 |
– |
edm::ESHandle<GlobalTrackingGeometry> globalGeometry; |
| 108 |
– |
es.get<GlobalTrackingGeometryRecord>().get(globalGeometry); |
| 109 |
– |
|
| 110 |
– |
edm::ESHandle<MagneticField> magField; |
| 111 |
– |
es.get<IdealMagneticFieldRecord>().get(magField); |
| 112 |
– |
|
| 113 |
– |
math::XYZPoint reference(0.,0.,0.); |
| 114 |
– |
if (theConfig.exists("beamSpot")) { |
| 115 |
– |
edm::InputTag beamSpot = theConfig.getParameter<edm::InputTag>("beamSpot"); |
| 116 |
– |
edm::Handle<reco::BeamSpot> bsHandle; |
| 117 |
– |
ev.getByLabel( beamSpot, bsHandle); |
| 118 |
– |
if (bsHandle.isValid()) reference = math::XYZPoint(bsHandle->x0(), bsHandle->y0(), bsHandle->z0()); |
| 119 |
– |
} |
| 120 |
– |
|
| 121 |
– |
reco::TrackCollection tracks = *(trackCollection.product()); |
| 122 |
– |
//cout <<"#RECONSTRUCTED tracks: " << tracks.size() << endl; |
| 123 |
– |
if (hNum) hNum->Fill(tracks.size()); |
| 124 |
– |
|
| 125 |
– |
mindeta = 100.; |
| 126 |
– |
mindphi = 100.; |
| 127 |
– |
|
| 128 |
– |
typedef reco::TrackCollection::const_iterator IT; |
| 129 |
– |
for (IT it = tracks.begin(); it !=tracks.end(); ++it) { |
| 130 |
– |
const reco::Track & track = *it; |
| 131 |
– |
if (track.pt() < theMinPt ) continue; |
| 132 |
– |
double dxy = track.dxy(reference); |
| 133 |
– |
if (dxy > theMaxTIP) continue; |
| 134 |
– |
|
| 135 |
– |
TrajectoryStateOnSurface aTSOS = TrajectoryStateTransform().outerStateOnSurface(track, *globalGeometry, magField.product()); |
| 136 |
– |
if (theConfig.getParameter<bool>("matchL1RPC") && !checkL1RpcMatching(aTSOS, ev,es) ) continue; |
| 137 |
– |
if (!checkRpcDetMatching(aTSOS,det,ev,es)) continue; |
| 138 |
– |
if (theConfig.getParameter<bool>("checkUniqueRecHitMatching") && !checkUniqueRecHitMatching(aTSOS,det,ev,es)) continue; |
| 139 |
– |
|
| 140 |
– |
if (hDxy) hDxy->Fill(dxy); |
| 141 |
– |
result = true; |
| 142 |
– |
} |
| 143 |
– |
|
| 144 |
– |
if (result && mindphi < 99.) hDeltaPhi->Fill(mindphi); |
| 145 |
– |
if (result && mindeta < 99.) hDeltaEta->Fill(mindeta); |
| 146 |
– |
|
| 147 |
– |
return result; |
| 148 |
– |
} |
| 149 |
– |
|
| 35 |
|
bool SynchroSelector::checkRpcDetMatching( const TrajectoryStateOnSurface & tsos, const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es) |
| 36 |
|
{ |
| 37 |
|
edm::ESHandle<Propagator> propagator; |
| 38 |
< |
es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator); |
| 38 |
> |
es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAlong", propagator); |
| 39 |
|
|
| 40 |
|
edm::ESHandle<GlobalTrackingGeometry> globalGeometry; |
| 41 |
|
es.get<GlobalTrackingGeometryRecord>().get(globalGeometry); |
| 42 |
|
|
| 43 |
|
GlobalPoint detPos, trackPos; |
| 44 |
+ |
// std::cout <<"BEFORE PROPAGATION POSITION:"<< tsos.globalPosition().perp()<<" GLOBAL MOMENTUM:"<<tsos.globalMomentum().mag()<< std::endl; |
| 45 |
|
TrajectoryStateOnSurface trackAtRPC = propagator->propagate(tsos, globalGeometry->idToDet(det)->surface()); |
| 46 |
|
if (!trackAtRPC.isValid()) return false; |
| 47 |
< |
//globalGeometry->idToDet(det)->specificSurface(); |
| 47 |
> |
// std::cout <<"AFTER PROPAGATION POSITION:"<< trackAtRPC.globalPosition().perp()<<" GLOBAL MOMENTUM:"<<trackAtRPC.globalMomentum().mag()<< std::endl; |
| 48 |
|
detPos = globalGeometry->idToDet(det)->position(); |
| 49 |
|
trackPos = trackAtRPC.globalPosition(); |
| 50 |
|
|
| 60 |
|
<<" isInside: "<< globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition()) |
| 61 |
|
<< std::endl; |
| 62 |
|
*/ |
| 63 |
< |
float scale = 0.; |
| 63 |
> |
float propQuality = 0.; |
| 64 |
|
if ( trackAtRPC.localError().positionError().xx() > 2500 || trackAtRPC.localError().positionError().yy() > 2500) return false; |
| 65 |
< |
else if (trackAtRPC.localError().positionError().xx() > 1000 || trackAtRPC.localError().positionError().yy() > 1000) scale = 0; |
| 66 |
< |
else if (trackAtRPC.localError().positionError().xx() > 500 || trackAtRPC.localError().positionError().yy() > 500) scale = 1; |
| 67 |
< |
else if (trackAtRPC.localError().positionError().xx() > 100 || trackAtRPC.localError().positionError().yy() > 100) scale = 2; |
| 68 |
< |
else scale = 3.; |
| 65 |
> |
else if (trackAtRPC.localError().positionError().xx() > 1000 || trackAtRPC.localError().positionError().yy() > 1000) propQuality = 0; |
| 66 |
> |
else if (trackAtRPC.localError().positionError().xx() > 500 || trackAtRPC.localError().positionError().yy() > 500) propQuality = 1; |
| 67 |
> |
else if (trackAtRPC.localError().positionError().xx() > 100 || trackAtRPC.localError().positionError().yy() > 100) propQuality = 2; |
| 68 |
> |
else propQuality = 3.; |
| 69 |
|
|
| 70 |
+ |
if (propQuality < theConfig.getParameter<int>("checkRpcDetMatching_minPropagationQuality") ) return false; |
| 71 |
+ |
double scale = theConfig.getParameter<bool>("checkRpcDetMatching_matchingScaleAuto") ? propQuality : theConfig.getParameter<double>("checkRpcDetMatching_matchingScaleValue") ; |
| 72 |
|
|
| 73 |
|
bool inside = globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition(), trackAtRPC.localError().positionError(), scale); |
| 74 |
< |
//if (inside) { std::cout <<" detector: " << detPos <<" track: "<< trackPos << std::endl; thePos.push_back(trackPos); |
| 74 |
> |
//std::cout<<"In:"<<inside <<" detector r:" << detPos.perp()<<" phi:"<<detPos.phi()<<" z:"<<detPos.z() |
| 75 |
> |
// <<" track r:"<< trackPos.perp()<<" phi:"<<trackPos.phi()<<" z:"<<trackPos.z() <<"Error: "<<trackAtRPC.localError().positionError()<< std::endl; |
| 76 |
|
return inside; |
| 77 |
|
} |
| 78 |
|
|
| 90 |
|
LocalError trackAtRPCError = trackAtRPC.localError().positionError(); |
| 91 |
|
|
| 92 |
|
edm::Handle<RPCRecHitCollection> rpcHits; |
| 93 |
< |
ev.getByType(rpcHits); |
| 93 |
> |
std::vector<edm::Handle<RPCRecHitCollection> > handle_vec; |
| 94 |
> |
ev.getManyByType(handle_vec); |
| 95 |
> |
if(handle_vec.size()) rpcHits = handle_vec[0]; |
| 96 |
|
bool matched = false; |
| 97 |
|
bool unique = true; |
| 98 |
|
pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> recHitsInDet = rpcHits->get(det); |
| 102 |
|
float distX = hitPoint.x()-trackAtRPCPoint.x(); |
| 103 |
|
float pullX = distX/ sqrt( trackAtRPCError.xx()+hitError.xx()); |
| 104 |
|
|
| 105 |
< |
if (fabs(pullX) < 2. && fabs(distX) < 20.) { |
| 105 |
> |
if ( fabs(pullX) < theConfig.getParameter<double>("checkUniqueRecHitMatching_maxPull") |
| 106 |
> |
&& fabs(distX) < theConfig.getParameter<double>("checkUniqueRecHitMatching_maxDist") ) { |
| 107 |
|
matched = true; |
| 108 |
|
} else { |
| 109 |
|
unique = false; |
