MuonIdentification/src/MuonShowerInformationFiller.cc

0001 /**

0002  *  \package: MuonIdentification

0003  *  \class: MuonShowerInformationFiller

0004  *  Description: class for muon shower identification

0005  *

0006

0007  *

0008  *  \author: A. Svyatkovskiy, Purdue University

0009  *

0010  **/
0011
0012 #include "RecoMuon/MuonIdentification/interface/MuonShowerInformationFiller.h"
0013
0014 // system include files

0015 #include <memory>
0016 #include <algorithm>
0017 #include <iostream>
0018
0019 // user include files

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 // Constructor

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 // Destructor

0108 //

0109 MuonShowerInformationFiller::~MuonShowerInformationFiller() {
0110   if (theService)
0111     delete theService;
0112 }
0113
0114 //

0115 //Fill the MuonShower struct

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   // Update the services

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 // Set Event

0143 //

0144 void MuonShowerInformationFiller::setEvent(const edm::Event& event) {
0145   // get all the necesary products

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 // Set services

0161 //

0162 void MuonShowerInformationFiller::setServices(const EventSetup& setup) {
0163   // DetLayer Geometry

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   // Transient Rechit Builders

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 // Get hits owned by segments

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     // loop on segments 4D

0193     DTRecSegment4DCollection::id_iterator chamberIdIt;
0194     for (chamberIdIt = dtSegments->id_begin(); chamberIdIt != dtSegments->id_end(); ++chamberIdIt) {
0195       if (*chamberIdIt != chamberId)
0196         continue;
0197
0198       // Get the range for the corresponding ChamberId

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       // Get the range for the corresponding ChamberId

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       //unused      DetId thisID = (*ihit1)->geographicalId();

0250       //LocalPoint lp1dinsegHit = (*ihit1)->localPosition();

0251       GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();
0252
0253       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
0254            ihit2 != allhitscorrelated.end();
0255            ++ihit2) {
0256         //unused        DetId thisID2 = (*ihit2)->geographicalId();

0257         //LocalPoint lp1dinsegHit2 = (*ihit2)->localPosition();

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 // Find cluster

0273 //

0274
0275 TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::findThetaCluster(
0276     TransientTrackingRecHit::ConstRecHitContainer& muonRecHits, const GlobalPoint& refpoint) const {
0277   if (muonRecHits.empty())
0278     return muonRecHits;
0279
0280   //clustering step by theta

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 //Used to treat overlap region

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 // Get compatible dets

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     // crossing points of track with cylinder

0353     GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const BarrelDetLayer*>(*iLayer));
0354
0355     // check if point is inside the detector

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     // check if point is inside the detector

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 // Intersection point of track with barrel layer

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   // there are two crossing points, return the one that is in the same quadrant as point of extrapolation

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 // Intersection point of track with a forward layer

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 // GeomDets along the track in DT

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 // GeomDets along the track in CSC

0583 //

0584 vector<const GeomDet*> MuonShowerInformationFiller::cscPositionToDets(const GlobalPoint& gp) const {
0585   // determine the endcap side

0586   int endcap = 0;
0587   if (gp.z() > 0) {
0588     endcap = 1;
0589   } else {
0590     endcap = 2;
0591   }
0592
0593   // determine the csc station and range of rings

0594   int station = 5;
0595
0596   // check all rings in a station

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.;  //for all the rings except first rings for stations > 1

0610   float phistep2 = M_PI / 9.;
0611   float phigp = (float)gp.barePhi();
0612
0613   int ring = -1;
0614
0615   // determine the ring

0616   if (station == 1) {
0617     //FIX ME!!!

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     // we have 18 sectors in that case

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     // we have 36 sectors in that case

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   // check exceptional cases

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 //Set class members

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   // split 1D rechits by station

0702   vector<MuonRecHitContainer> muonRecHits(4);
0703
0704   // split rechits from segs by station

0705   vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);
0706
0707   // get vector of GeomDets compatible with a track

0708   vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);
0709
0710   // for special cases: CSC station 1

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     // get det id

0717     DetId geoId = (*igd)->geographicalId();
0718
0719     // skip tracker hits

0720     if (geoId.det() != DetId::Muon)
0721       continue;
0722
0723     // DT

0724     if (geoId.subdetId() == MuonSubdetId::DT) {
0725       // get station

0726       DTChamberId detid(geoId.rawId());
0727       int station = detid.station();
0728       int wheel = detid.wheel();
0729
0730       // get rechits from segments per station

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       //check overlap certain wheels and stations

0741       if (abs(wheel) == 2 && station != 4 && station != 1)
0742         dtOverlapToCheck = true;
0743
0744       // loop over all superlayers of a DT chamber

0745       for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1;
0746            ++isuperlayer) {
0747         // loop over all layers inside the superlayer

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       // get station

0762       CSCDetId did(geoId.rawId());
0763       int station = did.station();
0764       int ring = did.ring();
0765
0766       //get rechits from segments by station

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       // split 1D rechits by station

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           //sort by perp, then insert to appropriate container

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   }  //loop over GeomDets compatible with a track

0812
0813   // calculate number of all and correlated hits

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   //station shower sizes

0827   MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiBest;
0828   TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
0829   // send station hits to the clustering algorithm

0830   for (int stat = 0; stat != 4; stat++) {
0831     const size_t nhit = muonRecHits[stat].size();
0832     if (nhit < 2)
0833       continue;  // Require at least 2 hits

0834     muonRecHitsPhiBest.clear();
0835     muonRecHitsPhiBest.reserve(nhit);
0836
0837     // Cluster seeds by global position phi. Best cluster is chosen to give greatest dphi

0838     // Sort by phi (complexity = NLogN with enough memory, or = NLog^2N for insufficient mem)

0839     stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi());
0840
0841     // Search for gaps (complexity = N)

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     //station shower sizes

0854     double dphimax = 0;
0855     if (clUppers.empty()) {
0856       // No gaps - there is only one cluster. Take all of them

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       // Loop over gaps to find the one with maximum dphi(begin, end)

0868       // By construction, number of gap must be greater than 1.

0869       size_t bestUpper = 0, bestLower = 0;
0870       for (auto icl = clUppers.begin(); icl != clUppers.end(); ++icl) {
0871         // upper bound of this cluster

0872         const size_t upper = *icl;
0873         // lower bound is +1 of preceeding upper bound

0874         const auto prevCl = (icl == clUppers.begin()) ? clUppers.end() : icl;
0875         const size_t lower = (*(prevCl - 1) + 1) % nhit;
0876
0877         // At least two hit for a cluster

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     //fill showerTs

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     //for theta

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();  //starting from the one with smallest value of phi

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         }  //at least two hits

0934       }    //loop over seeds

0935     }      //not empty container2

0936
0937     //fill deltaRs

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   }  //loop over station

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 }