CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Static Public Member Functions | Static Private Member Functions
MuonIdTruthInfo Class Reference

#include <MuonIdTruthInfo.h>

Public Member Functions

void registerConsumes (edm::ConsumesCollector &iC)
 

Static Public Member Functions

static void truthMatchMuon (const edm::Event &iEvent, const edm::EventSetup &iSetup, reco::Muon &aMuon)
 

Static Private Member Functions

static void checkSimHitForBestMatch (reco::MuonSegmentMatch &segmentMatch, double &distance, const PSimHit &hit, const DetId &chamberId, const edm::ESHandle< GlobalTrackingGeometry > &geometry)
 
static double matchChi2 (const reco::Track &recoTrk, const SimTrack &simTrk)
 

Detailed Description

Definition at line 21 of file MuonIdTruthInfo.h.

Member Function Documentation

void MuonIdTruthInfo::checkSimHitForBestMatch ( reco::MuonSegmentMatch segmentMatch,
double &  distance,
const PSimHit hit,
const DetId chamberId,
const edm::ESHandle< GlobalTrackingGeometry > &  geometry 
)
staticprivate

Definition at line 104 of file MuonIdTruthInfo.cc.

References PSimHit::detUnitId(), reco::MuonSegmentMatch::dXdZ, reco::MuonSegmentMatch::dXdZErr, reco::MuonSegmentMatch::dYdZ, reco::MuonSegmentMatch::dYdZErr, PSimHit::entryPoint(), PSimHit::exitPoint(), LogTrace, DetId::rawId(), GeomDet::toGlobal(), GeomDet::toLocal(), reco::MuonSegmentMatch::x, PV3DBase< T, PVType, FrameType >::x(), reco::MuonSegmentMatch::xErr, reco::MuonSegmentMatch::y, PV3DBase< T, PVType, FrameType >::y(), reco::MuonSegmentMatch::yErr, and PV3DBase< T, PVType, FrameType >::z().

Referenced by truthMatchMuon().

109 {
110  printf("DONT FORGET TO CALL REGISTERCONSUMES()\n");
111 
112  // find the hit position projection at the reference surface of the chamber:
113  // first get entry and exit point of the hit in the global coordinates, then
114  // get local coordinates of these points wrt the chamber and then find the
115  // projected X-Y coordinates
116 
117  const GeomDet* chamberGeometry = geometry->idToDet( chamberId );
118  const GeomDet* simUnitGeometry = geometry->idToDet( DetId(hit.detUnitId()) );
119 
120  if (chamberGeometry && simUnitGeometry ) {
121  LocalPoint entryPoint = chamberGeometry->toLocal( simUnitGeometry->toGlobal( hit.entryPoint() ) );
122  LocalPoint exitPoint = chamberGeometry->toLocal( simUnitGeometry->toGlobal( hit.exitPoint() ) );
123  LocalVector direction = exitPoint - entryPoint;
124  if ( fabs(direction.z()) > 0.001) {
125  LocalPoint projection = entryPoint - direction*(entryPoint.z()/direction.z());
126  if ( fabs(projection.z()) > 0.001 )
127  edm::LogWarning("MuonIdentification") << "z coordinate of the hit projection must be zero and it's not!\n";
128 
129  double new_distance = 99999;
130  if( entryPoint.z()*exitPoint.z() < -1 ) // changed sign, so the reference point is inside
131  new_distance = 0;
132  else {
133  if ( fabs(entryPoint.z()) < fabs(exitPoint.z()) )
134  new_distance = fabs(entryPoint.z());
135  else
136  new_distance = fabs(exitPoint.z());
137  }
138 
139  if (new_distance < distance) {
140  // find a SimHit closer to the reference surface, update segmentMatch
141  segmentMatch.x = projection.x();
142  segmentMatch.y = projection.y();
143  segmentMatch.xErr = 0;
144  segmentMatch.yErr = 0;
145  segmentMatch.dXdZ = direction.x()/direction.z();
146  segmentMatch.dYdZ = direction.y()/direction.z();
147  segmentMatch.dXdZErr = 0;
148  segmentMatch.dYdZErr = 0;
149  distance = new_distance;
150  LogTrace("MuonIdentificationVerbose") << "Better truth matched segment found:\n" <<
151  "\tDetId: " << chamberId.rawId() << "\n" <<
152  "\tentry point: ( " << entryPoint.x() << ", " << entryPoint.y() << ", " << entryPoint.z() << " )\n" <<
153  "\texit point: ( " << exitPoint.x() << ", " << exitPoint.y() << ", " << exitPoint.z() << " )\n" <<
154  "\tprojection: ( " << projection.x() << ", " << projection.y() << ", " << projection.z() << " )\n";
155  }
156  }
157  }else{
158  if ( ! chamberGeometry ) edm::LogWarning("MuonIdentification") << "Cannot get chamber geomtry for DetId: " << chamberId.rawId();
159  if ( ! simUnitGeometry ) edm::LogWarning("MuonIdentification") << "Cannot get detector unit geomtry for DetId: " << hit.detUnitId();
160  }
161 }
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:47
T y() const
Definition: PV3DBase.h:63
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:62
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
Local3DPoint exitPoint() const
Exit point in the local Det frame.
Definition: PSimHit.h:38
T z() const
Definition: PV3DBase.h:64
#define LogTrace(id)
Definition: DetId.h:18
T x() const
Definition: PV3DBase.h:62
Local3DPoint entryPoint() const
Entry point in the local Det frame.
Definition: PSimHit.h:35
unsigned int detUnitId() const
Definition: PSimHit.h:93
double MuonIdTruthInfo::matchChi2 ( const reco::Track recoTrk,
const SimTrack simTrk 
)
staticprivate

Definition at line 163 of file MuonIdTruthInfo.cc.

References SiPixelRawToDigiRegional_cfi::deltaPhi, CoreSimTrack::momentum(), reco::TrackBase::p(), reco::TrackBase::phi(), funct::pow(), and reco::TrackBase::theta().

Referenced by truthMatchMuon().

164 {
165  double deltaPhi = fabs(recoTrk.phi()-simTrk.momentum().phi());
166  if (deltaPhi>1.8*3.1416) deltaPhi = 2*3.1416-deltaPhi; // take care of phi discontinuity
167  return pow((recoTrk.p()-simTrk.momentum().rho())/simTrk.momentum().rho(),2)+
168  pow(deltaPhi/(2*3.1416),2)+pow((recoTrk.theta()-simTrk.momentum().theta())/3.1416,2);
169 }
double p() const
momentum vector magnitude
Definition: TrackBase.h:127
double theta() const
polar angle
Definition: TrackBase.h:115
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:137
const math::XYZTLorentzVectorD & momentum() const
Definition: CoreSimTrack.h:22
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void MuonIdTruthInfo::registerConsumes ( edm::ConsumesCollector iC)

Definition at line 11 of file MuonIdTruthInfo.cc.

References edm::ConsumesCollector::mayConsume().

11  {
12  iC.mayConsume<edm::SimTrackContainer>(edm::InputTag("g4SimHits",""));
13  iC.mayConsume<edm::PSimHitContainer>(edm::InputTag("g4SimHits","MuonDTHits"));
14  iC.mayConsume<edm::PSimHitContainer>(edm::InputTag("g4SimHits","MuonCSCHits"));
15 
16 }
EDGetTokenT< ProductType > mayConsume(edm::InputTag const &tag)
std::vector< PSimHit > PSimHitContainer
std::vector< SimTrack > SimTrackContainer
void MuonIdTruthInfo::truthMatchMuon ( const edm::Event iEvent,
const edm::EventSetup iSetup,
reco::Muon aMuon 
)
static

Definition at line 20 of file MuonIdTruthInfo.cc.

References checkSimHitForBestMatch(), MuonSubdetId::CSC, MuonSubdetId::DT, geometry, edm::EventSetup::get(), edm::Ref< C, T, F >::get(), edm::Event::getByLabel(), edm::ESHandleBase::isValid(), edm::HandleBase::isValid(), LogTrace, matchChi2(), reco::Muon::matches(), DetId::Muon, evf::evtn::offset(), trackerHits::simHits, reco::Muon::track(), reco::MuonSegmentMatch::x, and reco::MuonSegmentMatch::y.

Referenced by MuonIdProducer::produce().

23 {
24  // get a list of simulated track and find a track with the best match to
25  // the muon.track(). Use its id and chamber id to localize hits
26  // If a hit has non-zero local z coordinate, it's position wrt
27  // to the center of a chamber is extrapolated by a straight line
28 
30  iEvent.getByLabel<edm::SimTrackContainer>("g4SimHits", "", simTracks);
31  if (! simTracks.isValid() ) {
32  LogTrace("MuonIdentification") <<"No tracks found";
33  return;
34  }
35 
36  // get the tracking Geometry
38  iSetup.get<GlobalTrackingGeometryRecord>().get(geometry);
39  if ( ! geometry.isValid() )
40  throw cms::Exception("FatalError") << "Unable to find GlobalTrackingGeometryRecord in event!\n";
41 
42  float bestMatchChi2 = 9999; //minimization creteria
43  unsigned int bestMatch = 0;
44  const unsigned int offset = 0; // kludge to fix a problem in trackId matching between tracks and hits.
45 
46  for( edm::SimTrackContainer::const_iterator simTrk = simTracks->begin();
47  simTrk != simTracks->end(); simTrk++ )
48  {
49  float chi2 = matchChi2(*aMuon.track().get(),*simTrk);
50  if (chi2>bestMatchChi2) continue;
51  bestMatchChi2 = chi2;
52  bestMatch = simTrk->trackId();
53  }
54 
55  bestMatch -= offset;
56 
57  std::vector<reco::MuonChamberMatch>& matches = aMuon.matches();
58  int numberOfTruthMatchedChambers = 0;
59 
60  // loop over chambers
61  for(std::vector<reco::MuonChamberMatch>::iterator chamberMatch = matches.begin();
62  chamberMatch != matches.end(); chamberMatch ++)
63  {
64  if (chamberMatch->id.det() != DetId::Muon ) {
65  edm::LogWarning("MuonIdentification") << "Detector id of a muon chamber corresponds to not a muon detector";
66  continue;
67  }
68  reco::MuonSegmentMatch bestSegmentMatch;
69  double distance = 99999;
70 
71  if ( chamberMatch->id.subdetId() == MuonSubdetId::DT) {
72  DTChamberId detId(chamberMatch->id.rawId());
73 
75  iEvent.getByLabel("g4SimHits", "MuonDTHits", simHits);
76  if ( simHits.isValid() ) {
77  for( edm::PSimHitContainer::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++)
78  if ( hit->trackId() == bestMatch ) checkSimHitForBestMatch(bestSegmentMatch, distance, *hit, detId, geometry );
79  }else LogTrace("MuonIdentification") <<"No DT simulated hits are found";
80  }
81 
82  if ( chamberMatch->id.subdetId() == MuonSubdetId::CSC) {
83  CSCDetId detId(chamberMatch->id.rawId());
84 
86  iEvent.getByLabel("g4SimHits", "MuonCSCHits", simHits);
87  if ( simHits.isValid() ) {
88  for( edm::PSimHitContainer::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++)
89  if ( hit->trackId() == bestMatch ) checkSimHitForBestMatch(bestSegmentMatch, distance, *hit, detId, geometry );
90  }else LogTrace("MuonIdentification") <<"No CSC simulated hits are found";
91  }
92  if (distance < 9999) {
93  chamberMatch->truthMatches.push_back( bestSegmentMatch );
94  numberOfTruthMatchedChambers++;
95  LogTrace("MuonIdentification") << "Best truth matched hit:" <<
96  "\tDetId: " << chamberMatch->id.rawId() << "\n" <<
97  "\tprojection: ( " << bestSegmentMatch.x << ", " << bestSegmentMatch.y << " )\n";
98  }
99  }
100  LogTrace("MuonIdentification") << "Truth matching summary:\n\tnumber of chambers: " << matches.size() <<
101  "\n\tnumber of truth matched chambers: " << numberOfTruthMatchedChambers << "\n";
102 }
virtual TrackRef track() const
reference to a Track
Definition: Muon.h:49
static void checkSimHitForBestMatch(reco::MuonSegmentMatch &segmentMatch, double &distance, const PSimHit &hit, const DetId &chamberId, const edm::ESHandle< GlobalTrackingGeometry > &geometry)
static const int CSC
Definition: MuonSubdetId.h:13
unsigned int offset(bool)
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
#define LogTrace(id)
const T & get() const
Definition: EventSetup.h:55
tuple simHits
Definition: trackerHits.py:16
std::vector< MuonChamberMatch > & matches()
get muon matching information
Definition: Muon.h:140
ESHandle< TrackerGeometry > geometry
static const int DT
Definition: MuonSubdetId.h:12
bool isValid() const
Definition: ESHandle.h:37
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
std::vector< SimTrack > SimTrackContainer
static double matchChi2(const reco::Track &recoTrk, const SimTrack &simTrk)