File indexing completed on 2022-10-10 23:56:24
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include "RecoMuon/MuonIdentification/interface/MuonShowerInformationFiller.h"
0013
0014
0015 #include <memory>
0016 #include <algorithm>
0017 #include <iostream>
0018
0019
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/PluginManager/interface/PluginManager.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0025 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0026 #include "DataFormats/DetId/interface/DetId.h"
0027 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0028 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0029 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0030 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0031 #include "DataFormats/GeometrySurface/interface/Bounds.h"
0032 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0033 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0034 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0035 #include "MagneticField/Engine/interface/MagneticField.h"
0036 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0037 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0038 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0039 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0040 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0041 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0042 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0043 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0044 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0045 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0046 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0047 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
0048 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryBuilder.h"
0049 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0050 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0051 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0052
0053 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0054 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0055
0056 #include "DataFormats/MuonReco/interface/MuonShower.h"
0057 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0058 #include "DataFormats/MuonReco/interface/Muon.h"
0059 #include "DataFormats/TrackReco/interface/Track.h"
0060 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0061
0062 using namespace std;
0063 using namespace edm;
0064
0065
0066
0067
0068 MuonShowerInformationFiller::MuonShowerInformationFiller(const edm::ParameterSet& par, edm::ConsumesCollector& iC)
0069 : theService(nullptr),
0070 theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
0071 theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
0072 theCSCSegmentsLabel(par.getParameter<InputTag>("CSCSegmentLabel")),
0073 theDT4DRecSegmentLabel(par.getParameter<InputTag>("DT4DRecSegmentLabel")) {
0074 theDTRecHitToken = iC.consumes<DTRecHitCollection>(theDTRecHitLabel);
0075 theCSCRecHitToken = iC.consumes<CSCRecHit2DCollection>(theCSCRecHitLabel);
0076 theCSCSegmentsToken = iC.consumes<CSCSegmentCollection>(theCSCSegmentsLabel);
0077 theDT4DRecSegmentToken = iC.consumes<DTRecSegment4DCollection>(theDT4DRecSegmentLabel);
0078
0079 theTrackerToken = iC.esConsumes();
0080 theTrackingGeometryToken = iC.esConsumes();
0081 theFieldToken = iC.esConsumes();
0082 theCSCGeometryToken = iC.esConsumes();
0083 theDTGeometryToken = iC.esConsumes();
0084
0085 edm::ParameterSet serviceParameters = par.getParameter<edm::ParameterSet>("ServiceParameters");
0086 theService = new MuonServiceProxy(serviceParameters, edm::ConsumesCollector(iC));
0087
0088 auto trackerRecHitBuilderName = par.getParameter<string>("TrackerRecHitBuilder");
0089 theTrackerRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", trackerRecHitBuilderName));
0090 auto muonRecHitBuilderName = par.getParameter<string>("MuonRecHitBuilder");
0091 theMuonRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", muonRecHitBuilderName));
0092
0093 theCacheId_TRH = 0;
0094 theCacheId_MT = 0;
0095
0096 category_ = "MuonShowerInformationFiller";
0097
0098 for (int istat = 0; istat < 4; istat++) {
0099 theStationShowerDeltaR.push_back(0.);
0100 theStationShowerTSize.push_back(0.);
0101 theAllStationHits.push_back(0);
0102 theCorrelatedStationHits.push_back(0);
0103 }
0104 }
0105
0106
0107
0108
0109 MuonShowerInformationFiller::~MuonShowerInformationFiller() {
0110 if (theService)
0111 delete theService;
0112 }
0113
0114
0115
0116
0117 reco::MuonShower MuonShowerInformationFiller::fillShowerInformation(const reco::Muon& muon,
0118 const edm::Event& iEvent,
0119 const edm::EventSetup& iSetup) {
0120 reco::MuonShower returnShower;
0121
0122
0123 theService->update(iSetup);
0124 setEvent(iEvent);
0125 setServices(theService->eventSetup());
0126
0127 fillHitsByStation(muon);
0128 std::vector<int> nStationHits = theAllStationHits;
0129 std::vector<int> nStationCorrelatedHits = theCorrelatedStationHits;
0130 std::vector<float> stationShowerSizeT = theStationShowerTSize;
0131 std::vector<float> stationShowerDeltaR = theStationShowerDeltaR;
0132
0133 returnShower.nStationHits = nStationHits;
0134 returnShower.nStationCorrelatedHits = nStationCorrelatedHits;
0135 returnShower.stationShowerSizeT = stationShowerSizeT;
0136 returnShower.stationShowerDeltaR = stationShowerDeltaR;
0137
0138 return returnShower;
0139 }
0140
0141
0142
0143
0144 void MuonShowerInformationFiller::setEvent(const edm::Event& event) {
0145
0146 event.getByToken(theDTRecHitToken, theDTRecHits);
0147 event.getByToken(theCSCRecHitToken, theCSCRecHits);
0148 event.getByToken(theCSCSegmentsToken, theCSCSegments);
0149 event.getByToken(theDT4DRecSegmentToken, theDT4DRecSegments);
0150
0151 for (int istat = 0; istat < 4; istat++) {
0152 theStationShowerDeltaR.at(istat) = 0.;
0153 theStationShowerTSize.at(istat) = 0.;
0154 theAllStationHits.at(istat) = 0;
0155 theCorrelatedStationHits.at(istat) = 0;
0156 }
0157 }
0158
0159
0160
0161
0162 void MuonShowerInformationFiller::setServices(const EventSetup& setup) {
0163
0164 theTrackingGeometry = setup.getHandle(theTrackingGeometryToken);
0165 theField = setup.getHandle(theFieldToken);
0166 theTracker = setup.getHandle(theTrackerToken);
0167 theCSCGeometry = setup.getHandle(theCSCGeometryToken);
0168 theDTGeometry = setup.getHandle(theDTGeometryToken);
0169
0170
0171 unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
0172 if (newCacheId_TRH != theCacheId_TRH) {
0173 theTrackerRecHitBuilder = setup.getHandle(theTrackerRecHitBuilderToken);
0174 theMuonRecHitBuilder = setup.getHandle(theMuonRecHitBuilderToken);
0175 }
0176 }
0177
0178
0179
0180
0181 TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::hitsFromSegments(
0182 const GeomDet* geomDet,
0183 edm::Handle<DTRecSegment4DCollection> dtSegments,
0184 edm::Handle<CSCSegmentCollection> cscSegments) const {
0185 MuonTransientTrackingRecHit::MuonRecHitContainer segments;
0186
0187 DetId geoId = geomDet->geographicalId();
0188
0189 if (geoId.subdetId() == MuonSubdetId::DT) {
0190 DTChamberId chamberId(geoId.rawId());
0191
0192
0193 DTRecSegment4DCollection::id_iterator chamberIdIt;
0194 for (chamberIdIt = dtSegments->id_begin(); chamberIdIt != dtSegments->id_end(); ++chamberIdIt) {
0195 if (*chamberIdIt != chamberId)
0196 continue;
0197
0198
0199 DTRecSegment4DCollection::range range = dtSegments->get((*chamberIdIt));
0200
0201 for (DTRecSegment4DCollection::const_iterator iseg = range.first; iseg != range.second; ++iseg) {
0202 if (iseg->dimension() != 4)
0203 continue;
0204 segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*iseg));
0205 }
0206 }
0207 } else if (geoId.subdetId() == MuonSubdetId::CSC) {
0208 CSCDetId did(geoId.rawId());
0209
0210 for (CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin(); chamberId != cscSegments->id_end();
0211 ++chamberId) {
0212 if ((*chamberId).chamber() != did.chamber())
0213 continue;
0214
0215
0216 CSCSegmentCollection::range range = cscSegments->get((*chamberId));
0217
0218 for (CSCSegmentCollection::const_iterator iseg = range.first; iseg != range.second; ++iseg) {
0219 if (iseg->dimension() != 3)
0220 continue;
0221 segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*iseg));
0222 }
0223 }
0224 } else {
0225 LogTrace(category_) << "Segments are not built in RPCs" << endl;
0226 }
0227
0228 TransientTrackingRecHit::ConstRecHitContainer allhitscorrelated;
0229
0230 if (segments.empty())
0231 return allhitscorrelated;
0232
0233 TransientTrackingRecHit::ConstRecHitPointer muonRecHit(segments.front());
0234 allhitscorrelated = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, 2);
0235
0236 if (segments.size() == 1)
0237 return allhitscorrelated;
0238
0239 for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
0240 iseg != segments.end();
0241 ++iseg) {
0242 TransientTrackingRecHit::ConstRecHitPointer muonRecHit((*iseg));
0243 TransientTrackingRecHit::ConstRecHitContainer hits1 =
0244 MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, 2);
0245
0246 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin(); ihit1 != hits1.end();
0247 ++ihit1) {
0248 bool usedbefore = false;
0249
0250
0251 GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();
0252
0253 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
0254 ihit2 != allhitscorrelated.end();
0255 ++ihit2) {
0256
0257
0258 GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();
0259
0260 if ((gp1dinsegHit2 - gp1dinsegHit).mag() < 1.0)
0261 usedbefore = true;
0262 }
0263 if (!usedbefore)
0264 allhitscorrelated.push_back(*ihit1);
0265 }
0266 }
0267
0268 return allhitscorrelated;
0269 }
0270
0271
0272
0273
0274
0275 TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::findThetaCluster(
0276 TransientTrackingRecHit::ConstRecHitContainer& muonRecHits, const GlobalPoint& refpoint) const {
0277 if (muonRecHits.empty())
0278 return muonRecHits;
0279
0280
0281 float step = 0.05;
0282 TransientTrackingRecHit::ConstRecHitContainer result;
0283
0284 stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDTheta(refpoint));
0285
0286 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHits.begin();
0287 ihit != muonRecHits.end() - 1;
0288 ++ihit) {
0289 if (fabs((*(ihit + 1))->globalPosition().theta() - (*ihit)->globalPosition().theta()) < step) {
0290 result.push_back(*ihit);
0291 } else {
0292 break;
0293 }
0294 }
0295
0296 return result;
0297 }
0298
0299
0300
0301
0302 MuonTransientTrackingRecHit::MuonRecHitContainer MuonShowerInformationFiller::findPerpCluster(
0303 MuonTransientTrackingRecHit::MuonRecHitContainer& muonRecHits) const {
0304 if (muonRecHits.empty())
0305 return muonRecHits;
0306
0307 stable_sort(muonRecHits.begin(), muonRecHits.end(), LessPerp());
0308
0309 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator seedhit =
0310 min_element(muonRecHits.begin(), muonRecHits.end(), LessPerp());
0311
0312 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
0313 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;
0314
0315 float step = 0.1;
0316 while (ihigh != muonRecHits.end() - 1 &&
0317 (fabs((*(ihigh + 1))->globalPosition().perp() - (*ihigh)->globalPosition().perp()) < step)) {
0318 ihigh++;
0319 }
0320 while (ilow != muonRecHits.begin() &&
0321 (fabs((*ilow)->globalPosition().perp() - (*(ilow - 1))->globalPosition().perp()) < step)) {
0322 ilow--;
0323 }
0324
0325 MuonTransientTrackingRecHit::MuonRecHitContainer result(ilow, ihigh);
0326
0327 return result;
0328 }
0329
0330
0331
0332
0333 vector<const GeomDet*> MuonShowerInformationFiller::getCompatibleDets(const reco::Track& track) const {
0334 vector<const GeomDet*> total;
0335 total.reserve(1000);
0336
0337 LogTrace(category_) << "Consider a track " << track.p() << " eta: " << track.eta() << " phi " << track.phi() << endl;
0338
0339 TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(
0340 track, *theService->trackingGeometry(), &*theService->magneticField());
0341 TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(
0342 track, *theService->trackingGeometry(), &*theService->magneticField());
0343
0344 GlobalPoint innerPos = innerTsos.globalPosition();
0345 GlobalPoint outerPos = outerTsos.globalPosition();
0346
0347 vector<GlobalPoint> allCrossingPoints;
0348
0349 const vector<const DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();
0350
0351 for (auto iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {
0352
0353 GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const BarrelDetLayer*>(*iLayer));
0354
0355
0356 if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500) &&
0357 (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0)))
0358 allCrossingPoints.push_back(xPoint);
0359 }
0360
0361 stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos));
0362
0363 vector<const GeomDet*> tempDT;
0364
0365 for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
0366 tempDT = dtPositionToDets(*ipos);
0367 vector<const GeomDet*>::const_iterator begin = tempDT.begin();
0368 vector<const GeomDet*>::const_iterator end = tempDT.end();
0369
0370 for (; begin != end; ++begin) {
0371 total.push_back(*begin);
0372 }
0373 }
0374 allCrossingPoints.clear();
0375
0376 const vector<const DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
0377 for (auto iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {
0378 GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const ForwardDetLayer*>(*iLayer));
0379
0380
0381 if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500.0) &&
0382 (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0)))
0383 allCrossingPoints.push_back(xPoint);
0384 }
0385 stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos));
0386
0387 vector<const GeomDet*> tempCSC;
0388 for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
0389 tempCSC = cscPositionToDets(*ipos);
0390 vector<const GeomDet*>::const_iterator begin = tempCSC.begin();
0391 vector<const GeomDet*>::const_iterator end = tempCSC.end();
0392
0393 for (; begin != end; ++begin) {
0394 total.push_back(*begin);
0395 }
0396 }
0397
0398 return total;
0399 }
0400
0401
0402
0403
0404 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0405 const GlobalPoint& p2,
0406 const BarrelDetLayer* dl) const {
0407 const BoundCylinder& bc = dl->specificSurface();
0408 return crossingPoint(p1, p2, bc);
0409 }
0410
0411 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0412 const GlobalPoint& p2,
0413 const Cylinder& cyl) const {
0414 float radius = cyl.radius();
0415
0416 GlobalVector dp = p1 - p2;
0417 float slope = dp.x() / dp.y();
0418 float a = p1.x() - slope * p1.y();
0419
0420 float n2 = (1 + slope * slope);
0421 float n1 = 2 * a * slope;
0422 float n0 = a * a - radius * radius;
0423
0424 float y1 = 9999;
0425 float y2 = 9999;
0426 if (n1 * n1 - 4 * n2 * n0 > 0) {
0427 y1 = (-n1 + sqrt(n1 * n1 - 4 * n2 * n0)) / (2 * n2);
0428 y2 = (-n1 - sqrt(n1 * n1 - 4 * n2 * n0)) / (2 * n2);
0429 }
0430
0431 float x1 = p1.x() + slope * (y1 - p1.y());
0432 float x2 = p1.x() + slope * (y2 - p1.y());
0433
0434 float slopeZ = dp.z() / dp.y();
0435
0436 float z1 = p1.z() + slopeZ * (y1 - p1.y());
0437 float z2 = p1.z() + slopeZ * (y2 - p1.y());
0438
0439
0440 if ((p2.x() * x1 > 0) && (y1 * p2.y() > 0) && (z1 * p2.z() > 0)) {
0441 return GlobalPoint(x1, y1, z1);
0442 } else {
0443 return GlobalPoint(x2, y2, z2);
0444 }
0445 }
0446
0447
0448
0449
0450 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0451 const GlobalPoint& p2,
0452 const ForwardDetLayer* dl) const {
0453 const BoundDisk& bc = dl->specificSurface();
0454 return crossingPoint(p1, p2, bc);
0455 }
0456
0457 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0458 const GlobalPoint& p2,
0459 const BoundDisk& disk) const {
0460 float diskZ = disk.position().z();
0461 int endcap = diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
0462 diskZ = diskZ + endcap * dynamic_cast<const SimpleDiskBounds&>(disk.bounds()).thickness() / 2.;
0463
0464 GlobalVector dp = p1 - p2;
0465
0466 float slopeZ = dp.z() / dp.y();
0467 float y1 = diskZ / slopeZ;
0468
0469 float slopeX = dp.z() / dp.x();
0470 float x1 = diskZ / slopeX;
0471
0472 float z1 = diskZ;
0473
0474 if (p2.z() * z1 > 0) {
0475 return GlobalPoint(x1, y1, z1);
0476 } else {
0477 return GlobalPoint(0, 0, 0);
0478 }
0479 }
0480
0481
0482
0483
0484 vector<const GeomDet*> MuonShowerInformationFiller::dtPositionToDets(const GlobalPoint& gp) const {
0485 int minwheel = -3;
0486 int maxwheel = -3;
0487 if (gp.z() < -680.0) {
0488 minwheel = -3;
0489 maxwheel = -3;
0490 } else if (gp.z() < -396.0) {
0491 minwheel = -2;
0492 maxwheel = -1;
0493 } else if (gp.z() < -126.8) {
0494 minwheel = -2;
0495 maxwheel = 0;
0496 } else if (gp.z() < 126.8) {
0497 minwheel = -1;
0498 maxwheel = 1;
0499 } else if (gp.z() < 396.0) {
0500 minwheel = 0;
0501 maxwheel = 2;
0502 } else if (gp.z() < 680.0) {
0503 minwheel = 1;
0504 maxwheel = 2;
0505 } else {
0506 minwheel = 3;
0507 maxwheel = 3;
0508 }
0509
0510 int station = 5;
0511 if (gp.perp() > 680.0 && gp.perp() < 755.0)
0512 station = 4;
0513 else if (gp.perp() > 580.0)
0514 station = 3;
0515 else if (gp.perp() > 480.0)
0516 station = 2;
0517 else if (gp.perp() > 380.0)
0518 station = 1;
0519 else
0520 station = 0;
0521
0522 vector<int> sectors;
0523
0524 float phistep = M_PI / 6;
0525
0526 float phigp = (float)gp.barePhi();
0527
0528 if (fabs(deltaPhi(phigp, 0 * phistep)) < phistep)
0529 sectors.push_back(1);
0530 if (fabs(deltaPhi(phigp, phistep)) < phistep)
0531 sectors.push_back(2);
0532 if (fabs(deltaPhi(phigp, 2 * phistep)) < phistep)
0533 sectors.push_back(3);
0534 if (fabs(deltaPhi(phigp, 3 * phistep)) < phistep) {
0535 sectors.push_back(4);
0536 if (station == 4)
0537 sectors.push_back(13);
0538 }
0539 if (fabs(deltaPhi(phigp, 4 * phistep)) < phistep)
0540 sectors.push_back(5);
0541 if (fabs(deltaPhi(phigp, 5 * phistep)) < phistep)
0542 sectors.push_back(6);
0543 if (fabs(deltaPhi(phigp, 6 * phistep)) < phistep)
0544 sectors.push_back(7);
0545 if (fabs(deltaPhi(phigp, 7 * phistep)) < phistep)
0546 sectors.push_back(8);
0547 if (fabs(deltaPhi(phigp, 8 * phistep)) < phistep)
0548 sectors.push_back(9);
0549 if (fabs(deltaPhi(phigp, 9 * phistep)) < phistep) {
0550 sectors.push_back(10);
0551 if (station == 4)
0552 sectors.push_back(14);
0553 }
0554 if (fabs(deltaPhi(phigp, 10 * phistep)) < phistep)
0555 sectors.push_back(11);
0556 if (fabs(deltaPhi(phigp, 11 * phistep)) < phistep)
0557 sectors.push_back(12);
0558
0559 LogTrace(category_) << "DT position to dets" << endl;
0560 LogTrace(category_) << "number of sectors to consider: " << sectors.size() << endl;
0561 LogTrace(category_) << "station: " << station << " wheels: " << minwheel << " " << maxwheel << endl;
0562
0563 vector<const GeomDet*> result;
0564 if (station > 4 || station < 1)
0565 return result;
0566 if (minwheel > 2 || maxwheel < -2)
0567 return result;
0568
0569 for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
0570 for (int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
0571 DTChamberId chamberid(iwheel, station, (*isector));
0572 result.push_back(theService->trackingGeometry()->idToDet(chamberid));
0573 }
0574 }
0575
0576 LogTrace(category_) << "number of GeomDets for this track: " << result.size() << endl;
0577
0578 return result;
0579 }
0580
0581
0582
0583
0584 vector<const GeomDet*> MuonShowerInformationFiller::cscPositionToDets(const GlobalPoint& gp) const {
0585
0586 int endcap = 0;
0587 if (gp.z() > 0) {
0588 endcap = 1;
0589 } else {
0590 endcap = 2;
0591 }
0592
0593
0594 int station = 5;
0595
0596
0597 if (fabs(gp.z()) > 1000. && fabs(gp.z()) < 1055.0) {
0598 station = 4;
0599 } else if (fabs(gp.z()) > 910.0 && fabs(gp.z()) < 965.0) {
0600 station = 3;
0601 } else if (fabs(gp.z()) > 800.0 && fabs(gp.z()) < 860.0) {
0602 station = 2;
0603 } else if (fabs(gp.z()) > 570.0 && fabs(gp.z()) < 730.0) {
0604 station = 1;
0605 }
0606
0607 vector<int> sectors;
0608
0609 float phistep1 = M_PI / 18.;
0610 float phistep2 = M_PI / 9.;
0611 float phigp = (float)gp.barePhi();
0612
0613 int ring = -1;
0614
0615
0616 if (station == 1) {
0617
0618 if (gp.perp() > 100 && gp.perp() < 270)
0619 ring = 1;
0620 else if (gp.perp() > 270 && gp.perp() < 450)
0621 ring = 2;
0622 else if (gp.perp() > 450 && gp.perp() < 695)
0623 ring = 3;
0624 else if (gp.perp() > 100 && gp.perp() < 270)
0625 ring = 4;
0626
0627 } else if (station == 2) {
0628 if (gp.perp() > 140 && gp.perp() < 350)
0629 ring = 1;
0630 else if (gp.perp() > 350 && gp.perp() < 700)
0631 ring = 2;
0632
0633 } else if (station == 3) {
0634 if (gp.perp() > 160 && gp.perp() < 350)
0635 ring = 1;
0636 else if (gp.perp() > 350 && gp.perp() < 700)
0637 ring = 2;
0638
0639 } else if (station == 4) {
0640 if (gp.perp() > 175 && gp.perp() < 350)
0641 ring = 1;
0642 else if (gp.perp() > 350 && gp.perp() < 700)
0643 ring = 2;
0644 }
0645
0646 if (station > 1 && ring == 1) {
0647
0648 for (int i = 0; i < 18; i++) {
0649 if (fabs(deltaPhi(phigp, i * phistep2)) < phistep2)
0650 sectors.push_back(i + 1);
0651 }
0652
0653 } else {
0654
0655 for (int i = 0; i < 36; i++) {
0656 if (fabs(deltaPhi(phigp, i * phistep1)) < phistep1)
0657 sectors.push_back(i + 1);
0658 }
0659 }
0660
0661 LogTrace(category_) << "CSC position to dets" << endl;
0662 LogTrace(category_) << "ring: " << ring << endl;
0663 LogTrace(category_) << "endcap: " << endcap << endl;
0664 LogTrace(category_) << "station: " << station << endl;
0665 LogTrace(category_) << "CSC number of sectors to consider: " << sectors.size() << endl;
0666
0667
0668 vector<const GeomDet*> result;
0669 if (station > 4 || station < 1)
0670 return result;
0671 if (endcap == 0)
0672 return result;
0673 if (ring == -1)
0674 return result;
0675
0676 int minlayer = 1;
0677 int maxlayer = 6;
0678
0679 for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
0680 for (int ilayer = minlayer; ilayer != maxlayer + 1; ++ilayer) {
0681 CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
0682 result.push_back(theService->trackingGeometry()->idToDet(cscid));
0683 }
0684 }
0685
0686 return result;
0687 }
0688
0689
0690
0691
0692 void MuonShowerInformationFiller::fillHitsByStation(const reco::Muon& muon) {
0693 reco::TrackRef track;
0694 if (muon.isGlobalMuon())
0695 track = muon.globalTrack();
0696 else if (muon.isStandAloneMuon())
0697 track = muon.outerTrack();
0698 else
0699 return;
0700
0701
0702 vector<MuonRecHitContainer> muonRecHits(4);
0703
0704
0705 vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);
0706
0707
0708 vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);
0709
0710
0711 MuonRecHitContainer tmpCSC1;
0712 bool dtOverlapToCheck = false;
0713 bool cscOverlapToCheck = false;
0714
0715 for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++) {
0716
0717 DetId geoId = (*igd)->geographicalId();
0718
0719
0720 if (geoId.det() != DetId::Muon)
0721 continue;
0722
0723
0724 if (geoId.subdetId() == MuonSubdetId::DT) {
0725
0726 DTChamberId detid(geoId.rawId());
0727 int station = detid.station();
0728 int wheel = detid.wheel();
0729
0730
0731 TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp =
0732 hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
0733 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
0734 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
0735
0736 for (; hits_begin != hits_end; ++hits_begin) {
0737 muonCorrelatedHits.at(station - 1).push_back(*hits_begin);
0738 }
0739
0740
0741 if (abs(wheel) == 2 && station != 4 && station != 1)
0742 dtOverlapToCheck = true;
0743
0744
0745 for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1;
0746 ++isuperlayer) {
0747
0748 for (int ilayer = DTChamberId::minLayerId; ilayer != DTChamberId::maxLayerId + 1; ++ilayer) {
0749 DTLayerId lid(detid, isuperlayer, ilayer);
0750 DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
0751 for (DTRecHitCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
0752 vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
0753 for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end();
0754 ++irechit) {
0755 muonRecHits.at(station - 1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &**irechit));
0756 }
0757 }
0758 }
0759 }
0760 } else if (geoId.subdetId() == MuonSubdetId::CSC) {
0761
0762 CSCDetId did(geoId.rawId());
0763 int station = did.station();
0764 int ring = did.ring();
0765
0766
0767 TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp =
0768 hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
0769 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
0770 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
0771
0772 for (; hits_begin != hits_end; ++hits_begin) {
0773 muonCorrelatedHits.at(station - 1).push_back(*hits_begin);
0774 }
0775
0776 if ((station == 1 && ring == 3) && dtOverlapToCheck)
0777 cscOverlapToCheck = true;
0778
0779
0780 CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
0781 for (CSCRecHit2DCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
0782 if (!cscOverlapToCheck) {
0783 muonRecHits.at(station - 1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &*rechit));
0784 } else {
0785 tmpCSC1.push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &*rechit));
0786
0787
0788 MuonRecHitContainer temp = findPerpCluster(tmpCSC1);
0789 if (temp.empty())
0790 continue;
0791
0792 float center;
0793 if (temp.size() > 1) {
0794 center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp()) / 2.;
0795 } else {
0796 center = temp.front()->globalPosition().perp();
0797 }
0798 temp.clear();
0799
0800 if (center > 550.) {
0801 muonRecHits.at(2).insert(muonRecHits.at(2).end(), tmpCSC1.begin(), tmpCSC1.end());
0802 } else {
0803 muonRecHits.at(1).insert(muonRecHits.at(1).end(), tmpCSC1.begin(), tmpCSC1.end());
0804 }
0805 tmpCSC1.clear();
0806 }
0807 }
0808 } else if (geoId.subdetId() == MuonSubdetId::RPC) {
0809 LogTrace(category_) << "Wrong subdet id" << endl;
0810 }
0811 }
0812
0813
0814 for (int stat = 0; stat < 4; stat++) {
0815 theCorrelatedStationHits[stat] = muonCorrelatedHits.at(stat).size();
0816 theAllStationHits[stat] = muonRecHits[stat].size();
0817 }
0818 LogTrace(category_) << "Hits used by the segments, by station " << theCorrelatedStationHits.at(0) << " "
0819 << theCorrelatedStationHits.at(1) << " " << theCorrelatedStationHits.at(2) << " "
0820 << theCorrelatedStationHits.at(3) << endl;
0821
0822 LogTrace(category_) << "All DT 1D/CSC 2D hits, by station " << theAllStationHits.at(0) << " "
0823 << theAllStationHits.at(1) << " " << theAllStationHits.at(2) << " " << theAllStationHits.at(3)
0824 << endl;
0825
0826
0827 MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiBest;
0828 TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
0829
0830 for (int stat = 0; stat != 4; stat++) {
0831 const size_t nhit = muonRecHits[stat].size();
0832 if (nhit < 2)
0833 continue;
0834 muonRecHitsPhiBest.clear();
0835 muonRecHitsPhiBest.reserve(nhit);
0836
0837
0838
0839 stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi());
0840
0841
0842 std::vector<size_t> clUppers;
0843 for (size_t ihit = 0; ihit < nhit; ++ihit) {
0844 const size_t jhit = (ihit + 1) % nhit;
0845 const double phi1 = muonRecHits[stat].at(ihit)->globalPosition().barePhi();
0846 const double phi2 = muonRecHits[stat].at(jhit)->globalPosition().barePhi();
0847
0848 const double dphi = std::abs(deltaPhi(phi1, phi2));
0849 if (dphi >= 0.05)
0850 clUppers.push_back(ihit);
0851 }
0852
0853
0854 double dphimax = 0;
0855 if (clUppers.empty()) {
0856
0857 const double refPhi = muonRecHits[stat].at(0)->globalPosition().barePhi();
0858 double dphilo = 0, dphihi = 0;
0859 for (auto& hit : muonRecHits[stat]) {
0860 muonRecHitsPhiBest.push_back(hit);
0861 const double phi = hit->globalPosition().barePhi();
0862 dphilo = std::min(dphilo, deltaPhi(refPhi, phi));
0863 dphihi = std::max(dphihi, deltaPhi(refPhi, phi));
0864 }
0865 dphimax = std::abs(dphihi + dphilo);
0866 } else {
0867
0868
0869 size_t bestUpper = 0, bestLower = 0;
0870 for (auto icl = clUppers.begin(); icl != clUppers.end(); ++icl) {
0871
0872 const size_t upper = *icl;
0873
0874 const auto prevCl = (icl == clUppers.begin()) ? clUppers.end() : icl;
0875 const size_t lower = (*(prevCl - 1) + 1) % nhit;
0876
0877
0878 if (upper == lower)
0879 continue;
0880
0881 const double phi1 = muonRecHits[stat].at(upper)->globalPosition().barePhi();
0882 const double phi2 = muonRecHits[stat].at(lower)->globalPosition().barePhi();
0883
0884 const double dphi = std::abs(deltaPhi(phi1, phi2));
0885 if (dphimax < dphi) {
0886 dphimax = dphi;
0887 bestUpper = upper;
0888 bestLower = lower;
0889 }
0890 }
0891
0892 if (bestUpper > bestLower) {
0893 muonRecHitsPhiBest.reserve(bestUpper - bestLower + 1);
0894 std::copy(muonRecHits[stat].begin() + bestLower,
0895 muonRecHits[stat].begin() + bestUpper + 1,
0896 std::back_inserter(muonRecHitsPhiBest));
0897 } else if (bestUpper < bestLower) {
0898 muonRecHitsPhiBest.reserve(bestUpper + (nhit - bestLower) + 1);
0899 std::copy(muonRecHits[stat].begin(),
0900 muonRecHits[stat].begin() + bestUpper + 1,
0901 std::back_inserter(muonRecHitsPhiBest));
0902 std::copy(
0903 muonRecHits[stat].begin() + bestLower, muonRecHits[stat].end(), std::back_inserter(muonRecHitsPhiBest));
0904 }
0905 }
0906
0907
0908 if (!muonRecHitsPhiBest.empty()) {
0909 muonRecHits[stat] = muonRecHitsPhiBest;
0910 stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessAbsMag());
0911 GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
0912 theStationShowerTSize.at(stat) = refpoint.mag() * dphimax;
0913 }
0914
0915
0916 if (!muonCorrelatedHits.at(stat).empty()) {
0917 float dthetamax = 0;
0918 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iseed = muonCorrelatedHits.at(stat).begin();
0919 iseed != muonCorrelatedHits.at(stat).end();
0920 ++iseed) {
0921 if (!(*iseed)->isValid())
0922 continue;
0923 GlobalPoint refpoint = (*iseed)->globalPosition();
0924 muonRecHitsThetaTemp.clear();
0925 muonRecHitsThetaTemp = findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
0926 if (muonRecHitsThetaTemp.size() > 1) {
0927 float dtheta = fabs((float)muonRecHitsThetaTemp.back()->globalPosition().theta() -
0928 (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
0929 if (dtheta > dthetamax) {
0930 dthetamax = dtheta;
0931 muonRecHitsThetaBest = muonRecHitsThetaTemp;
0932 }
0933 }
0934 }
0935 }
0936
0937
0938 if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
0939 theStationShowerDeltaR.at(stat) = sqrt(pow(deltaPhi(muonRecHitsPhiBest.front()->globalPosition().barePhi(),
0940 muonRecHitsPhiBest.back()->globalPosition().barePhi()),
0941 2) +
0942 pow(muonRecHitsThetaBest.front()->globalPosition().theta() -
0943 muonRecHitsThetaBest.back()->globalPosition().theta(),
0944 2));
0945
0946 }
0947
0948 LogTrace(category_) << "deltaR around a track containing all the station hits, by station "
0949 << theStationShowerDeltaR.at(0) << " " << theStationShowerDeltaR.at(1) << " "
0950 << theStationShowerDeltaR.at(2) << " " << theStationShowerDeltaR.at(3) << endl;
0951
0952 LogTrace(category_) << "Transverse cluster size, by station " << theStationShowerTSize.at(0) << " "
0953 << theStationShowerTSize.at(1) << " " << theStationShowerTSize.at(2) << " "
0954 << theStationShowerTSize.at(3) << endl;
0955
0956 return;
0957 }