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 |
|
8 |
#include "DataFormats/TrackReco/interface/TrackFwd.h"
|
9 |
#include "DataFormats/TrackReco/interface/Track.h"
|
10 |
|
11 |
|
12 |
#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
|
13 |
#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
|
14 |
|
15 |
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
|
16 |
#include "MagneticField/Engine/interface/MagneticField.h"
|
17 |
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
|
18 |
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
|
19 |
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
|
20 |
|
21 |
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
|
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), hPullX(0),hDistX(0)
|
33 |
{ }
|
34 |
|
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("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 |
// 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::checkUniqueRecHitMatching( const TrajectoryStateOnSurface & tsos, const RPCDetId & det, const edm::Event&ev, const edm::EventSetup& es)
|
80 |
{
|
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 |
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 |
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);
|
99 |
for (RPCRecHitCollection::const_iterator ih = recHitsInDet.first; ih != recHitsInDet.second; ++ih) {
|
100 |
LocalPoint hitPoint = ih->localPosition();
|
101 |
LocalError hitError = ih->localPositionError();
|
102 |
float distX = hitPoint.x()-trackAtRPCPoint.x();
|
103 |
float pullX = distX/ sqrt( trackAtRPCError.xx()+hitError.xx());
|
104 |
|
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;
|
110 |
}
|
111 |
if (hPullX) hPullX->Fill(pullX);
|
112 |
if (hDistX) hDistX->Fill(distX);
|
113 |
}
|
114 |
|
115 |
return (matched && unique);
|
116 |
}
|