ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/L1RpcTriggerAnalysis/src/SynchroSelector.cc

Comparing UserCode/L1RpcTriggerAnalysis/src/SynchroSelector.cc (file contents):
Revision 1.1 by konec, Wed May 26 06:26:45 2010 UTC vs.
Revision 1.5 by konec, Tue Oct 30 11:13:16 2012 UTC

#	Line 1 | Line 1
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"
6		#include "FWCore/Utilities/interface/InputTag.h"
7
7	–
8		#include "DataFormats/TrackReco/interface/TrackFwd.h"
9		#include "DataFormats/TrackReco/interface/Track.h"
10	–	#include "DataFormats/TrackReco/interface/print.h"
10
12	–	#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h"
13	–	#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTCand.h"
14	–	#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTExtendedCand.h"
15	–	#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
11
12		#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
13		#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
#	Line 24 | Line 19
19		#include "TrackingTools/GeomPropagators/interface/Propagator.h"
20
21		#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
27	–	#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
28	–	#include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
29	–	#include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
30	–	#include "DataFormats/GeometrySurface/interface/BoundDisk.h"
31	–	#include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
32	–	#include "DataFormats/BeamSpot/interface/BeamSpot.h"
33	–
22
23		#include "DataFormats/MuonDetId/interface/RPCDetId.h"
24	+	#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
25	+
26		#include "TH1.h"
27		#include <cmath>
28		using namespace edm;
29		using namespace std;
30
31		SynchroSelector::SynchroSelector(const edm::ParameterSet & cfg)
32	<	: theConfig(cfg), hDxy(0),hNum(0),hDeltaEta(0),hDeltaPhi(0)
32	>	: theConfig(cfg), hPullX(0),hDistX(0)
33		{ }
34
35	<	bool SynchroSelector::checkMatching( const TrajectoryStateOnSurface & tsos,
46	<	float rpcEta, float rpcPhi, const edm::Event&ev, const edm::EventSetup& es)
35	>	bool SynchroSelector::checkRpcDetMatching( const TrajectoryStateOnSurface & tsos,  const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es)
36		{
48	–	edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
49	–	es.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
50	–	edm::ESHandle<MagneticField> magField;
51	–	es.get<IdealMagneticFieldRecord>().get(magField);
52	–	ReferenceCountingPointer<Surface> rpc;
53	–	if (rpcEta < -1.04)      rpc = ReferenceCountingPointer<Surface>(new  BoundDisk( GlobalPoint(0.,0.,-800.), TkRotation<float>(), SimpleDiskBounds( 300., 710., -10., 10. ) ) );
54	–	else if (rpcEta < 1.04)  rpc = ReferenceCountingPointer<Surface>(new  BoundCylinder( GlobalPoint(0.,0.,0.), TkRotation<float>(), SimpleCylinderBounds( 500, 520, -700, 700 ) ) );
55	–	else                     rpc = ReferenceCountingPointer<Surface>(new  BoundDisk( GlobalPoint(0.,0.,800.), TkRotation<float>(), SimpleDiskBounds( 300., 710., -10., 10. ) ) );
37		edm::ESHandle<Propagator> propagator;
38	<	es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
58	<	TrajectoryStateOnSurface trackAtRPC =  propagator->propagate(tsos, *rpc);
59	<	if (!trackAtRPC.isValid()) return false;
60	<	float phi = trackAtRPC.globalPosition().phi();
61	<	float dphi = phi-rpcPhi;
62	<	if (dphi < -M_PI) dphi+=2*M_PI;
63	<	if (dphi > M_PI) dphi-=2*M_PI;
64	<	float eta = trackAtRPC.globalPosition().eta();
65	<	float deta = eta-rpcEta;
66	<	if (fabs(dphi) < mindphi ) mindphi = fabs(dphi);
67	<	if (fabs(deta) < mindeta ) mindeta = fabs(deta);
38	>	es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAlong", propagator);
39
40	+	edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
41	+	es.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
42
43	<	double maxDeltaEta = theConfig.getParameter<double>("maxDeltaEta");
44	<	double maxDeltaPhi = theConfig.getParameter<double>("maxDeltaPhi");
45	<	return ( fabs(dphi) < maxDeltaPhi && fabs(deta) < maxDeltaEta) ? true : false;
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	>	//  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	>
51	>	/*
52	>	std::cout <<" **** "<<theConfig.getParameter<std::string>("collection")<<" position at RPC det:"<<det.rawId()
53	>	//<< is :"<globalGeometry->idToDet(det)->position()
54	>	<<", r= "<<trackAtRPC.globalPosition().perp()
55	>	<<", z= "<<trackAtRPC.globalPosition().z()
56	>	<<", phi= "<<trackAtRPC.globalPosition().phi()
57	>	<<", eta= "<<trackAtRPC.globalPosition().eta()
58	>	//<<" localPosition: "<<trackAtRPC.localPosition()
59	>	//<<" localError: (xx:"<<trackAtRPC.localError().positionError()
60	>	<<" isInside: "<< globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition())
61	>	<< std::endl;
62	>	*/
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)  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	>	//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
79	<	bool SynchroSelector::checkTriggerMatching( const TrajectoryStateOnSurface & tsos,  const edm::Event&ev, const edm::EventSetup& es)
79	>	bool SynchroSelector::checkUniqueRecHitMatching( const TrajectoryStateOnSurface & tsos,  const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es)
80		{
81	<	edm::Handle<L1MuGMTReadoutCollection> pCollection;
82	<	InputTag l1MuReadout = theConfig.getParameter<edm::InputTag>("l1MuReadout");
83	<	ev.getByLabel(l1MuReadout,pCollection);
80	<	L1MuGMTReadoutCollection const* gmtrc = pCollection.product();
81	<	if (!gmtrc) return false;
82	<
83	<	vector<L1MuGMTReadoutRecord> gmt_records = gmtrc->getRecords();
84	<	vector<L1MuGMTReadoutRecord>::const_iterator RRItr;
85	<	typedef vector<L1MuRegionalCand>::const_iterator ITC;
86	<	for( RRItr = gmt_records.begin() ; RRItr != gmt_records.end() ; RRItr++ ) {
87	<	for (int i=1; i<=2; ++i) {
88	<	vector<L1MuRegionalCand> Cands = (i==1) ? RRItr->getBrlRPCCands() : RRItr->getFwdRPCCands();
89	<	for(ITC it = Cands.begin() ; it != Cands.end() ; ++it ) {
90	<	if (it->empty()) continue;
91	<	if (checkMatching(tsos, it->etaValue(), it->phiValue(),ev,es)) 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	<	std::string theOption = theConfig.getParameter<std::string>("collection");
101	<	bool theCheckL1Matching = theConfig.getParameter<bool>("matchL1RPC");
102	<	float theMaxTIP = theConfig.getParameter<double>("maxTIP");
103	<	float theMinPt = theConfig.getParameter<double>("minPt");
104	<
105	<	edm::Handle<reco::TrackCollection> trackCollection;
106	<	ev.getByLabel(InputTag(theOption),trackCollection);
107	<
81	>	//propagate to det
82	>	edm::ESHandle<Propagator> propagator;
83	>	es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
84		edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
85		es.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
86	<
87	<	edm::ESHandle<MagneticField> magField;
88	<	es.get<IdealMagneticFieldRecord>().get(magField);
89	<
90	<	math::XYZPoint reference(0.,0.,0.);
91	<	if (theConfig.exists("beamSpot")) {
92	<	edm::InputTag beamSpot =  theConfig.getParameter<edm::InputTag>("beamSpot");
93	<	edm::Handle<reco::BeamSpot> bsHandle;
94	<	ev.getByLabel( beamSpot, bsHandle);
95	<	if (bsHandle.isValid()) reference = math::XYZPoint(bsHandle->x0(), bsHandle->y0(), bsHandle->z0());
96	<	}
97	<
98	<
99	<	reco::TrackCollection tracks = *(trackCollection.product());
100	<	//cout <<"#RECONSTRUCTED tracks: " << tracks.size() << endl;
101	<	if (hNum) hNum->Fill(tracks.size());
102	<
103	<	bool inside = false;
104	<	mindeta = 100.;
105	<	mindphi = 100.;
106	<
107	<	typedef reco::TrackCollection::const_iterator IT;
108	<	GlobalPoint detPos, trackPos;
109	<	for (IT it = tracks.begin(); it !=tracks.end(); ++it) {
110	<	const reco::Track & track = *it;
135	<	if (track.pt() < theMinPt ) continue;
136	<	double dxy = track.dxy(reference);
137	<
138	<	TrajectoryStateOnSurface aTSOS = TrajectoryStateTransform().outerStateOnSurface(track, *globalGeometry, magField.product());
139	<	if (theCheckL1Matching && !checkTriggerMatching(aTSOS, ev,es) ) continue;
140	<
141	<	edm::ESHandle<Propagator> propagator;
142	<	es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
143	<
144	<	TrajectoryStateOnSurface trackAtRPC =  propagator->propagate(aTSOS, globalGeometry->idToDet(det)->surface());
145	<	if (!trackAtRPC.isValid()) continue;
146	<	globalGeometry->idToDet(det)->specificSurface();
147	<	detPos = globalGeometry->idToDet(det)->position();
148	<	trackPos = trackAtRPC.globalPosition();
149	<	std::cout <<" **** "<<theOption<<" position at RPC det:"<<det.rawId()
150	<	//<< is :"<globalGeometry->idToDet(det)->position()
151	<	<<", r= "<<trackAtRPC.globalPosition().perp()
152	<	<<", z= "<<trackAtRPC.globalPosition().z()
153	<	<<", phi= "<<trackAtRPC.globalPosition().phi()
154	<	<<", eta= "<<trackAtRPC.globalPosition().eta()
155	<	//<<" localPosition: "<<trackAtRPC.localPosition()
156	<	//<<" localError: (xx:"<<trackAtRPC.localError().positionError()
157	<	<<" isInside: "<< globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition())
158	<	<< std::endl;
159	<	float scale = 0.;
160	<	if (     trackAtRPC.localError().positionError().xx() > 2500 || trackAtRPC.localError().positionError().yy() > 2500) continue;
161	<	else if (trackAtRPC.localError().positionError().xx() > 1000 || trackAtRPC.localError().positionError().yy() > 1000)  scale = 0;
162	<	else if (trackAtRPC.localError().positionError().xx() > 500  || trackAtRPC.localError().positionError().yy() > 500)   scale = 1;
163	<	else if (trackAtRPC.localError().positionError().xx() > 100  || trackAtRPC.localError().positionError().yy() > 100)   scale = 2;
164	<	else scale = 3.;
165	<
166	<	if (globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition(), trackAtRPC.localError().positionError(), scale)) {
167	<	if (hDxy) hDxy->Fill(dxy);
168	<	if (dxy < theMaxTIP) inside = true;
169	<	}
170	<	}
171	<	if (mindphi < 99.) hDeltaPhi->Fill(mindphi);
172	<	if (mindeta < 99.) hDeltaEta->Fill(mindeta);
173	<
174	<	if (inside) {
175	<	std::cout <<" detector: " << detPos <<" track: "<< trackPos << std::endl;
176	<	thePos.push_back(trackPos);
86	>	TrajectoryStateOnSurface trackAtRPC =  propagator->propagate(tsos, globalGeometry->idToDet(det)->surface());
87	>	if (!trackAtRPC.isValid()) return false;
88	>
89	>	LocalPoint trackAtRPCPoint = trackAtRPC.localPosition() ;
90	>	LocalError trackAtRPCError = trackAtRPC.localError().positionError();
91	>
92	>	edm::Handle<RPCRecHitCollection> rpcHits;
93	>	ev.getByType(rpcHits);
94	>	bool matched = false;
95	>	bool unique = true;
96	>	pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> recHitsInDet =  rpcHits->get(det);
97	>	for (RPCRecHitCollection::const_iterator ih = recHitsInDet.first; ih != recHitsInDet.second; ++ih) {
98	>	LocalPoint hitPoint = ih->localPosition();
99	>	LocalError hitError = ih->localPositionError();
100	>	float distX = hitPoint.x()-trackAtRPCPoint.x();
101	>	float pullX = distX/ sqrt( trackAtRPCError.xx()+hitError.xx());
102	>
103	>	if (    fabs(pullX) < theConfig.getParameter<double>("checkUniqueRecHitMatching_maxPull")
104	>	&& fabs(distX) < theConfig.getParameter<double>("checkUniqueRecHitMatching_maxDist") ) {
105	>	matched = true;
106	>	} else {
107	>	unique = false;
108	>	}
109	>	if (hPullX) hPullX->Fill(pullX);
110	>	if (hDistX) hDistX->Fill(distX);
111		}
112	<
113	<	return inside;
114	<	}
112	>
113	>	return (matched && unique);
114	>	}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines