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.3 by konec, Sun Jun 13 22:16:47 2010 UTC

#	Line 4 | Line 4
4		#include "FWCore/Framework/interface/ESHandle.h"
5		#include "FWCore/Utilities/interface/InputTag.h"
6
7	–
7		#include "DataFormats/TrackReco/interface/TrackFwd.h"
8		#include "DataFormats/TrackReco/interface/Track.h"
9		#include "DataFormats/TrackReco/interface/print.h"
#	Line 33 | Line 32
32
33
34		#include "DataFormats/MuonDetId/interface/RPCDetId.h"
35	+	#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
36	+
37		#include "TH1.h"
38		#include <cmath>
39		using namespace edm;
40		using namespace std;
41
42		SynchroSelector::SynchroSelector(const edm::ParameterSet & cfg)
43	<	: theConfig(cfg), hDxy(0),hNum(0),hDeltaEta(0),hDeltaPhi(0)
43	>	: theConfig(cfg), hDxy(0),hNum(0),hDeltaEta(0),hDeltaPhi(0),
44	>	hPullX(0),hDistX(0)
45		{ }
46
47	<	bool SynchroSelector::checkMatching( const TrajectoryStateOnSurface & tsos,
46	<	float rpcEta, float rpcPhi, const edm::Event&ev, const edm::EventSetup& es)
47	<	{
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. ) ) );
56	<	edm::ESHandle<Propagator> propagator;
57	<	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);
68	<
69	<
70	<	double maxDeltaEta = theConfig.getParameter<double>("maxDeltaEta");
71	<	double maxDeltaPhi = theConfig.getParameter<double>("maxDeltaPhi");
72	<	return ( fabs(dphi) < maxDeltaPhi && fabs(deta) < maxDeltaEta) ? true : false;
73	<	}
74	<
75	<	bool SynchroSelector::checkTriggerMatching( const TrajectoryStateOnSurface & tsos,  const edm::Event&ev, const edm::EventSetup& es)
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");
#	Line 80 | Line 52 | bool SynchroSelector::checkTriggerMatchi
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;
#	Line 88 | Line 66 | bool SynchroSelector::checkTriggerMatchi
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	<	if (checkMatching(tsos, it->etaValue(), it->phiValue(),ev,es)) return true;
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("SteppingHelixPropagatorAlong", 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		}
#	Line 97 | Line 97 | bool SynchroSelector::checkTriggerMatchi
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");
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(theOption),trackCollection);
105	>	ev.getByLabel(InputTag( theConfig.getParameter<std::string>("collection")),trackCollection);
106
107		edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
108		es.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
#	Line 119 | Line 118 | bool SynchroSelector::takeIt(const RPCDe
118		if (bsHandle.isValid()) reference = math::XYZPoint(bsHandle->x0(), bsHandle->y0(), bsHandle->z0());
119		}
120
122	–
121		reco::TrackCollection tracks = *(trackCollection.product());
122		//cout <<"#RECONSTRUCTED tracks: " << tracks.size() << endl;
123		if (hNum) hNum->Fill(tracks.size());
124
127	–	bool inside = false;
125		mindeta = 100.;
126		mindphi = 100.;
127
128		typedef reco::TrackCollection::const_iterator IT;
132	–	GlobalPoint detPos, trackPos;
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	<
133	>	if (dxy > theMaxTIP) continue;
134	>
135		TrajectoryStateOnSurface aTSOS = TrajectoryStateTransform().outerStateOnSurface(track, *globalGeometry, magField.product());
136	<	if (theCheckL1Matching && !checkTriggerMatching(aTSOS, ev,es) ) continue;
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	<	edm::ESHandle<Propagator> propagator;
141	<	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	<	}
140	>	if (hDxy) hDxy->Fill(dxy);
141	>	result = true;
142		}
171	–	if (mindphi < 99.) hDeltaPhi->Fill(mindphi);
172	–	if (mindeta < 99.) hDeltaEta->Fill(mindeta);
143
144	<	if (inside) {
145	<	std::cout <<" detector: " << detPos <<" track: "<< trackPos << std::endl;
176	<	thePos.push_back(trackPos);
177	<	}
144	>	if (result && mindphi < 99.) hDeltaPhi->Fill(mindphi);
145	>	if (result && mindeta < 99.) hDeltaEta->Fill(mindeta);
146
147	<	return inside;
147	>	return result;
148		}
149	+
150	+	bool SynchroSelector::checkRpcDetMatching( const TrajectoryStateOnSurface & tsos,  const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es)
151	+	{
152	+	edm::ESHandle<Propagator> propagator;
153	+	es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAlong", propagator);
154	+
155	+	edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
156	+	es.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
157	+
158	+	GlobalPoint detPos, trackPos;
159	+	TrajectoryStateOnSurface trackAtRPC =  propagator->propagate(tsos, globalGeometry->idToDet(det)->surface());
160	+	if (!trackAtRPC.isValid()) return false;
161	+	//globalGeometry->idToDet(det)->specificSurface();
162	+	detPos = globalGeometry->idToDet(det)->position();
163	+	trackPos = trackAtRPC.globalPosition();
164	+
165	+	/*
166	+	std::cout <<" **** "<<theConfig.getParameter<std::string>("collection")<<" position at RPC det:"<<det.rawId()
167	+	//<< is :"<globalGeometry->idToDet(det)->position()
168	+	<<", r= "<<trackAtRPC.globalPosition().perp()
169	+	<<", z= "<<trackAtRPC.globalPosition().z()
170	+	<<", phi= "<<trackAtRPC.globalPosition().phi()
171	+	<<", eta= "<<trackAtRPC.globalPosition().eta()
172	+	//<<" localPosition: "<<trackAtRPC.localPosition()
173	+	//<<" localError: (xx:"<<trackAtRPC.localError().positionError()
174	+	<<" isInside: "<< globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition())
175	+	<< std::endl;
176	+	*/
177	+	float scale = 0.;
178	+	if (     trackAtRPC.localError().positionError().xx() > 2500 || trackAtRPC.localError().positionError().yy() > 2500) return false;
179	+	else if (trackAtRPC.localError().positionError().xx() > 1000 || trackAtRPC.localError().positionError().yy() > 1000)  scale = 0;
180	+	else if (trackAtRPC.localError().positionError().xx() > 500  || trackAtRPC.localError().positionError().yy() > 500)   scale = 1;
181	+	else if (trackAtRPC.localError().positionError().xx() > 100  || trackAtRPC.localError().positionError().yy() > 100)   scale = 2;
182	+	else scale = 3.;
183	+
184	+
185	+	bool inside = globalGeometry->idToDet(det)->surface().bounds().inside( trackAtRPC.localPosition(), trackAtRPC.localError().positionError(), scale);
186	+	//if (inside) { std::cout <<" detector: " << detPos <<" track: "<< trackPos << std::endl; thePos.push_back(trackPos);
187	+	return inside;
188	+	}
189	+
190	+	bool SynchroSelector::checkUniqueRecHitMatching( const TrajectoryStateOnSurface & tsos,  const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es)
191	+	{
192	+	//propagate to det
193	+	edm::ESHandle<Propagator> propagator;
194	+	es.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAlong", propagator);
195	+	edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
196	+	es.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
197	+	TrajectoryStateOnSurface trackAtRPC =  propagator->propagate(tsos, globalGeometry->idToDet(det)->surface());
198	+	if (!trackAtRPC.isValid()) return false;
199	+
200	+	LocalPoint trackAtRPCPoint = trackAtRPC.localPosition() ;
201	+	LocalError trackAtRPCError = trackAtRPC.localError().positionError();
202	+
203	+	edm::Handle<RPCRecHitCollection> rpcHits;
204	+	ev.getByType(rpcHits);
205	+	bool matched = false;
206	+	bool unique = true;
207	+	pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> recHitsInDet =  rpcHits->get(det);
208	+	for (RPCRecHitCollection::const_iterator ih = recHitsInDet.first; ih != recHitsInDet.second; ++ih) {
209	+	LocalPoint hitPoint = ih->localPosition();
210	+	LocalError hitError = ih->localPositionError();
211	+	float distX = hitPoint.x()-trackAtRPCPoint.x();
212	+	float pullX = distX/ sqrt( trackAtRPCError.xx()+hitError.xx());
213	+
214	+	if (fabs(pullX) < 2. && fabs(distX) < 20.) {
215	+	matched = true;
216	+	} else {
217	+	unique = false;
218	+	}
219	+	if (hPullX) hPullX->Fill(pullX);
220	+	if (hDistX) hDistX->Fill(distX);
221	+	}
222	+
223	+	return (matched && unique);
224	+	}

Diff Legend

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