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 #include <numeric>
0019
0020 // user include files

0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/PluginManager/interface/PluginManager.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0026 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0029 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0030 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0031 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0032 #include "DataFormats/GeometrySurface/interface/Bounds.h"
0033 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0034 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0035 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0036 #include "MagneticField/Engine/interface/MagneticField.h"
0037 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0038 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0039 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0040 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0041 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0042 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0043 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0044 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0045 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0046 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0047 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0048 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
0049 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryBuilder.h"
0050 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0051 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0052 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0053
0054 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0055 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0056
0057 #include "DataFormats/MuonReco/interface/MuonShower.h"
0058 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0059 #include "DataFormats/MuonReco/interface/Muon.h"
0060 #include "DataFormats/TrackReco/interface/Track.h"
0061 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0062
0063 using namespace std;
0064 using namespace edm;
0065
0066 namespace {
0067   //

0068   // Sort function optimized for the case where the desired comparison

0069   // is extravagantly expensive.  To use, replace

0070   //

0071   // auto LessPhi = [](const auto& lhs, const auto& rhs) {

0072   //   return (lhs->globalPosition().barePhi() < rhs->globalPosition().barePhi())

0073   // }

0074   // stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi);

0075   //

0076   // with

0077   //

0078   // auto calcPhi = [](const auto& hit) { return hit->globalPosition().barePhi(); };

0079   // muonRecHits[stat] = sort_all_indexed(muonRecHits[stat], std::less(), calcPhi);

0080   //

0081   // This calculates the values to be compared once in O(n) time, instead of

0082   // O(n*log(n)) times in the comparison function.

0083   //

0084   template <typename RandomAccessSequence, typename Predicate, typename Transform>
0085   RandomAccessSequence sort_all_indexed(const RandomAccessSequence& s, Predicate p, Transform t) {
0086     std::vector<size_t> idx(s.size());
0087     std::iota(idx.begin(), idx.end(), 0);
0088
0089     // fill the cache of the values to be sorted

0090     using valueCacheType = std::invoke_result_t<decltype(t), typename RandomAccessSequence::value_type>;
0091     std::vector<valueCacheType> valcache(s.size());
0092     std::transform(s.begin(), s.end(), valcache.begin(), t);
0093
0094     // sort the indices of the value cache

0095     auto idxComp = [&valcache, p](auto i1, auto i2) { return p(valcache[i1], valcache[i2]); };
0096     std::stable_sort(idx.begin(), idx.end(), idxComp);
0097
0098     // fill the sorted output vector

0099     RandomAccessSequence r(s.size());
0100     for (size_t i = 0; i < s.size(); ++i) {
0101       r[i] = s[idx[i]];
0102     }
0103     return r;
0104   }
0105 }  // namespace

0106
0107 //

0108 // Constructor

0109 //

0110 MuonShowerInformationFiller::MuonShowerInformationFiller(const edm::ParameterSet& par, edm::ConsumesCollector& iC)
0111     : theService(nullptr),
0112       theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
0113       theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
0114       theCSCSegmentsLabel(par.getParameter<InputTag>("CSCSegmentLabel")),
0115       theDT4DRecSegmentLabel(par.getParameter<InputTag>("DT4DRecSegmentLabel")) {
0116   theDTRecHitToken = iC.consumes<DTRecHitCollection>(theDTRecHitLabel);
0117   theCSCRecHitToken = iC.consumes<CSCRecHit2DCollection>(theCSCRecHitLabel);
0118   theCSCSegmentsToken = iC.consumes<CSCSegmentCollection>(theCSCSegmentsLabel);
0119   theDT4DRecSegmentToken = iC.consumes<DTRecSegment4DCollection>(theDT4DRecSegmentLabel);
0120
0121   theTrackerToken = iC.esConsumes();
0122   theTrackingGeometryToken = iC.esConsumes();
0123   theFieldToken = iC.esConsumes();
0124   theCSCGeometryToken = iC.esConsumes();
0125   theDTGeometryToken = iC.esConsumes();
0126
0127   edm::ParameterSet serviceParameters = par.getParameter<edm::ParameterSet>("ServiceParameters");
0128   theService = new MuonServiceProxy(serviceParameters, edm::ConsumesCollector(iC));
0129
0130   auto trackerRecHitBuilderName = par.getParameter<string>("TrackerRecHitBuilder");
0131   theTrackerRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", trackerRecHitBuilderName));
0132   auto muonRecHitBuilderName = par.getParameter<string>("MuonRecHitBuilder");
0133   theMuonRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", muonRecHitBuilderName));
0134
0135   theCacheId_TRH = 0;
0136   theCacheId_MT = 0;
0137
0138   category_ = "MuonShowerInformationFiller";
0139
0140   for (int istat = 0; istat < 4; istat++) {
0141     theStationShowerDeltaR.push_back(0.);
0142     theStationShowerTSize.push_back(0.);
0143     theAllStationHits.push_back(0);
0144     theCorrelatedStationHits.push_back(0);
0145   }
0146 }
0147
0148 //

0149 // Destructor

0150 //

0151 MuonShowerInformationFiller::~MuonShowerInformationFiller() {
0152   if (theService)
0153     delete theService;
0154 }
0155
0156 //

0157 //Fill the MuonShower struct

0158 //

0159 reco::MuonShower MuonShowerInformationFiller::fillShowerInformation(const reco::Muon& muon,
0160                                                                     const edm::Event& iEvent,
0161                                                                     const edm::EventSetup& iSetup) {
0162   reco::MuonShower returnShower;
0163
0164   // Update the services

0165   theService->update(iSetup);
0166   setEvent(iEvent);
0167   setServices(theService->eventSetup());
0168
0169   fillHitsByStation(muon);
0170   std::vector<int> nStationHits = theAllStationHits;
0171   std::vector<int> nStationCorrelatedHits = theCorrelatedStationHits;
0172   std::vector<float> stationShowerSizeT = theStationShowerTSize;
0173   std::vector<float> stationShowerDeltaR = theStationShowerDeltaR;
0174
0175   returnShower.nStationHits = nStationHits;
0176   returnShower.nStationCorrelatedHits = nStationCorrelatedHits;
0177   returnShower.stationShowerSizeT = stationShowerSizeT;
0178   returnShower.stationShowerDeltaR = stationShowerDeltaR;
0179
0180   return returnShower;
0181 }
0182
0183 //

0184 // Set Event

0185 //

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

0188   event.getByToken(theDTRecHitToken, theDTRecHits);
0189   event.getByToken(theCSCRecHitToken, theCSCRecHits);
0190   event.getByToken(theCSCSegmentsToken, theCSCSegments);
0191   event.getByToken(theDT4DRecSegmentToken, theDT4DRecSegments);
0192
0193   for (int istat = 0; istat < 4; istat++) {
0194     theStationShowerDeltaR.at(istat) = 0.;
0195     theStationShowerTSize.at(istat) = 0.;
0196     theAllStationHits.at(istat) = 0;
0197     theCorrelatedStationHits.at(istat) = 0;
0198   }
0199 }
0200
0201 //

0202 // Set services

0203 //

0204 void MuonShowerInformationFiller::setServices(const EventSetup& setup) {
0205   // DetLayer Geometry

0206   theTrackingGeometry = setup.getHandle(theTrackingGeometryToken);
0207   theField = setup.getHandle(theFieldToken);
0208   theTracker = setup.getHandle(theTrackerToken);
0209   theCSCGeometry = setup.getHandle(theCSCGeometryToken);
0210   theDTGeometry = setup.getHandle(theDTGeometryToken);
0211
0212   // Transient Rechit Builders

0213   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
0214   if (newCacheId_TRH != theCacheId_TRH) {
0215     theTrackerRecHitBuilder = setup.getHandle(theTrackerRecHitBuilderToken);
0216     theMuonRecHitBuilder = setup.getHandle(theMuonRecHitBuilderToken);
0217   }
0218 }
0219
0220 //

0221 // Get hits owned by segments

0222 //

0223 TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::hitsFromSegments(
0224     const GeomDet* geomDet,
0225     edm::Handle<DTRecSegment4DCollection> dtSegments,
0226     edm::Handle<CSCSegmentCollection> cscSegments) const {
0227   MuonTransientTrackingRecHit::MuonRecHitContainer segments;
0228
0229   DetId geoId = geomDet->geographicalId();
0230
0231   if (geoId.subdetId() == MuonSubdetId::DT) {
0232     DTChamberId chamberId(geoId.rawId());
0233
0234     // loop on segments 4D

0235     DTRecSegment4DCollection::id_iterator chamberIdIt;
0236     for (chamberIdIt = dtSegments->id_begin(); chamberIdIt != dtSegments->id_end(); ++chamberIdIt) {
0237       if (*chamberIdIt != chamberId)
0238         continue;
0239
0240       // Get the range for the corresponding ChamberId

0241       DTRecSegment4DCollection::range range = dtSegments->get((*chamberIdIt));
0242
0243       for (DTRecSegment4DCollection::const_iterator iseg = range.first; iseg != range.second; ++iseg) {
0244         if (iseg->dimension() != 4)
0245           continue;
0246         segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*iseg));
0247       }
0248     }
0249   } else if (geoId.subdetId() == MuonSubdetId::CSC) {
0250     CSCDetId did(geoId.rawId());
0251
0252     for (CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin(); chamberId != cscSegments->id_end();
0253          ++chamberId) {
0254       if ((*chamberId).chamber() != did.chamber())
0255         continue;
0256
0257       // Get the range for the corresponding ChamberId

0258       CSCSegmentCollection::range range = cscSegments->get((*chamberId));
0259
0260       for (CSCSegmentCollection::const_iterator iseg = range.first; iseg != range.second; ++iseg) {
0261         if (iseg->dimension() != 3)
0262           continue;
0263         segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*iseg));
0264       }
0265     }
0266   } else {
0267     LogTrace(category_) << "Segments are not built in RPCs" << endl;
0268   }
0269
0270   TransientTrackingRecHit::ConstRecHitContainer allhitscorrelated;
0271
0272   if (segments.empty())
0273     return allhitscorrelated;
0274
0275   TransientTrackingRecHit::ConstRecHitPointer muonRecHit(segments.front());
0276   allhitscorrelated = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, 2);
0277
0278   if (segments.size() == 1)
0279     return allhitscorrelated;
0280
0281   for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
0282        iseg != segments.end();
0283        ++iseg) {
0284     TransientTrackingRecHit::ConstRecHitPointer muonRecHit((*iseg));
0285     TransientTrackingRecHit::ConstRecHitContainer hits1 =
0286         MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, 2);
0287
0288     for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin(); ihit1 != hits1.end();
0289          ++ihit1) {
0290       bool usedbefore = false;
0291       //unused      DetId thisID = (*ihit1)->geographicalId();

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

0293       GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();
0294
0295       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
0296            ihit2 != allhitscorrelated.end();
0297            ++ihit2) {
0298         //unused        DetId thisID2 = (*ihit2)->geographicalId();

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

0300         GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();
0301
0302         if ((gp1dinsegHit2 - gp1dinsegHit).mag() < 1.0)
0303           usedbefore = true;
0304       }
0305       if (!usedbefore)
0306         allhitscorrelated.push_back(*ihit1);
0307     }
0308   }
0309
0310   return allhitscorrelated;
0311 }
0312
0313 //

0314 // Find cluster

0315 //

0316
0317 TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::findThetaCluster(
0318     const TransientTrackingRecHit::ConstRecHitContainer& muonRecHitsIn, const GlobalPoint& refpoint) const {
0319   if (muonRecHitsIn.empty())
0320     return muonRecHitsIn;
0321
0322   //clustering step by theta

0323   float step = 0.05;
0324   TransientTrackingRecHit::ConstRecHitContainer result;
0325
0326   auto muonRecHitsTmp = sort_all_indexed(muonRecHitsIn, std::less(), AbsDThetaTransform(refpoint));
0327
0328   for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHitsTmp.begin();
0329        ihit != muonRecHitsTmp.end() - 1;
0330        ++ihit) {
0331     if (fabs((*(ihit + 1))->globalPosition().theta() - (*ihit)->globalPosition().theta()) < step) {
0332       result.push_back(*ihit);
0333     } else {
0334       break;
0335     }
0336   }
0337
0338   return result;
0339 }
0340
0341 //

0342 //Used to treat overlap region

0343 //

0344 MuonTransientTrackingRecHit::MuonRecHitContainer MuonShowerInformationFiller::findPerpCluster(
0345     const MuonTransientTrackingRecHit::MuonRecHitContainer& muonRecHitsIn) const {
0346   if (muonRecHitsIn.empty())
0347     return muonRecHitsIn;
0348
0349   auto muonRecHitsTmp = sort_all_indexed(muonRecHitsIn, std::less(), PerpTransform());
0350
0351   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator seedhit =
0352       min_element(muonRecHitsTmp.begin(), muonRecHitsTmp.end(), LessPerp());
0353
0354   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
0355   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;
0356
0357   float step = 0.1;
0358   while (ihigh != muonRecHitsTmp.end() - 1 &&
0359          (fabs((*(ihigh + 1))->globalPosition().perp() - (*ihigh)->globalPosition().perp()) < step)) {
0360     ihigh++;
0361   }
0362   while (ilow != muonRecHitsTmp.begin() &&
0363          (fabs((*ilow)->globalPosition().perp() - (*(ilow - 1))->globalPosition().perp()) < step)) {
0364     ilow--;
0365   }
0366
0367   MuonTransientTrackingRecHit::MuonRecHitContainer result(ilow, ihigh);
0368
0369   return result;
0370 }
0371
0372 //

0373 // Get compatible dets

0374 //

0375 vector<const GeomDet*> MuonShowerInformationFiller::getCompatibleDets(const reco::Track& track) const {
0376   vector<const GeomDet*> total;
0377   total.reserve(1000);
0378
0379   LogTrace(category_) << "Consider a track " << track.p() << " eta: " << track.eta() << " phi " << track.phi() << endl;
0380
0381   TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(
0382       track, *theService->trackingGeometry(), &*theService->magneticField());
0383   TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(
0384       track, *theService->trackingGeometry(), &*theService->magneticField());
0385
0386   GlobalPoint innerPos = innerTsos.globalPosition();
0387   GlobalPoint outerPos = outerTsos.globalPosition();
0388
0389   vector<GlobalPoint> allCrossingPoints;
0390
0391   const vector<const DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();
0392
0393   for (auto iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {
0394     // crossing points of track with cylinder

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

0398     if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500) &&
0399         (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0)))
0400       allCrossingPoints.push_back(xPoint);
0401   }
0402
0403   allCrossingPoints = sort_all_indexed(allCrossingPoints, std::less(), MagTransform(innerPos));
0404
0405   vector<const GeomDet*> tempDT;
0406
0407   for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
0408     tempDT = dtPositionToDets(*ipos);
0409     vector<const GeomDet*>::const_iterator begin = tempDT.begin();
0410     vector<const GeomDet*>::const_iterator end = tempDT.end();
0411
0412     for (; begin != end; ++begin) {
0413       total.push_back(*begin);
0414     }
0415   }
0416   allCrossingPoints.clear();
0417
0418   const vector<const DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
0419   for (auto iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {
0420     GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const ForwardDetLayer*>(*iLayer));
0421
0422     // check if point is inside the detector

0423     if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500.0) &&
0424         (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0)))
0425       allCrossingPoints.push_back(xPoint);
0426   }
0427   allCrossingPoints = sort_all_indexed(allCrossingPoints, std::less(), MagTransform(innerPos));
0428
0429   vector<const GeomDet*> tempCSC;
0430   for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
0431     tempCSC = cscPositionToDets(*ipos);
0432     vector<const GeomDet*>::const_iterator begin = tempCSC.begin();
0433     vector<const GeomDet*>::const_iterator end = tempCSC.end();
0434
0435     for (; begin != end; ++begin) {
0436       total.push_back(*begin);
0437     }
0438   }
0439
0440   return total;
0441 }
0442
0443 //

0444 // Intersection point of track with barrel layer

0445 //

0446 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0447                                                        const GlobalPoint& p2,
0448                                                        const BarrelDetLayer* dl) const {
0449   const BoundCylinder& bc = dl->specificSurface();
0450   return crossingPoint(p1, p2, bc);
0451 }
0452
0453 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0454                                                        const GlobalPoint& p2,
0455                                                        const Cylinder& cyl) const {
0456   float radius = cyl.radius();
0457
0458   GlobalVector dp = p1 - p2;
0459   float slope = dp.x() / dp.y();
0460   float a = p1.x() - slope * p1.y();
0461
0462   float n2 = (1 + slope * slope);
0463   float n1 = 2 * a * slope;
0464   float n0 = a * a - radius * radius;
0465
0466   float y1 = 9999;
0467   float y2 = 9999;
0468   if (n1 * n1 - 4 * n2 * n0 > 0) {
0469     y1 = (-n1 + sqrt(n1 * n1 - 4 * n2 * n0)) / (2 * n2);
0470     y2 = (-n1 - sqrt(n1 * n1 - 4 * n2 * n0)) / (2 * n2);
0471   }
0472
0473   float x1 = p1.x() + slope * (y1 - p1.y());
0474   float x2 = p1.x() + slope * (y2 - p1.y());
0475
0476   float slopeZ = dp.z() / dp.y();
0477
0478   float z1 = p1.z() + slopeZ * (y1 - p1.y());
0479   float z2 = p1.z() + slopeZ * (y2 - p1.y());
0480
0481   // there are two crossing points, return the one that is in the same quadrant as point of extrapolation

0482   if ((p2.x() * x1 > 0) && (y1 * p2.y() > 0) && (z1 * p2.z() > 0)) {
0483     return GlobalPoint(x1, y1, z1);
0484   } else {
0485     return GlobalPoint(x2, y2, z2);
0486   }
0487 }
0488
0489 //

0490 // Intersection point of track with a forward layer

0491 //

0492 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0493                                                        const GlobalPoint& p2,
0494                                                        const ForwardDetLayer* dl) const {
0495   const BoundDisk& bc = dl->specificSurface();
0496   return crossingPoint(p1, p2, bc);
0497 }
0498
0499 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
0500                                                        const GlobalPoint& p2,
0501                                                        const BoundDisk& disk) const {
0502   float diskZ = disk.position().z();
0503   int endcap = diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
0504   diskZ = diskZ + endcap * dynamic_cast<const SimpleDiskBounds&>(disk.bounds()).thickness() / 2.;
0505
0506   GlobalVector dp = p1 - p2;
0507
0508   float slopeZ = dp.z() / dp.y();
0509   float y1 = diskZ / slopeZ;
0510
0511   float slopeX = dp.z() / dp.x();
0512   float x1 = diskZ / slopeX;
0513
0514   float z1 = diskZ;
0515
0516   if (p2.z() * z1 > 0) {
0517     return GlobalPoint(x1, y1, z1);
0518   } else {
0519     return GlobalPoint(0, 0, 0);
0520   }
0521 }
0522
0523 //

0524 // GeomDets along the track in DT

0525 //

0526 vector<const GeomDet*> MuonShowerInformationFiller::dtPositionToDets(const GlobalPoint& gp) const {
0527   int minwheel = -3;
0528   int maxwheel = -3;
0529   if (gp.z() < -680.0) {
0530     minwheel = -3;
0531     maxwheel = -3;
0532   } else if (gp.z() < -396.0) {
0533     minwheel = -2;
0534     maxwheel = -1;
0535   } else if (gp.z() < -126.8) {
0536     minwheel = -2;
0537     maxwheel = 0;
0538   } else if (gp.z() < 126.8) {
0539     minwheel = -1;
0540     maxwheel = 1;
0541   } else if (gp.z() < 396.0) {
0542     minwheel = 0;
0543     maxwheel = 2;
0544   } else if (gp.z() < 680.0) {
0545     minwheel = 1;
0546     maxwheel = 2;
0547   } else {
0548     minwheel = 3;
0549     maxwheel = 3;
0550   }
0551
0552   int station = 5;
0553   if (gp.perp() > 680.0 && gp.perp() < 755.0)
0554     station = 4;
0555   else if (gp.perp() > 580.0)
0556     station = 3;
0557   else if (gp.perp() > 480.0)
0558     station = 2;
0559   else if (gp.perp() > 380.0)
0560     station = 1;
0561   else
0562     station = 0;
0563
0564   vector<int> sectors;
0565
0566   float phistep = M_PI / 6;
0567
0568   float phigp = (float)gp.barePhi();
0569
0570   if (fabs(deltaPhi(phigp, 0 * phistep)) < phistep)
0571     sectors.push_back(1);
0572   if (fabs(deltaPhi(phigp, phistep)) < phistep)
0573     sectors.push_back(2);
0574   if (fabs(deltaPhi(phigp, 2 * phistep)) < phistep)
0575     sectors.push_back(3);
0576   if (fabs(deltaPhi(phigp, 3 * phistep)) < phistep) {
0577     sectors.push_back(4);
0578     if (station == 4)
0579       sectors.push_back(13);
0580   }
0581   if (fabs(deltaPhi(phigp, 4 * phistep)) < phistep)
0582     sectors.push_back(5);
0583   if (fabs(deltaPhi(phigp, 5 * phistep)) < phistep)
0584     sectors.push_back(6);
0585   if (fabs(deltaPhi(phigp, 6 * phistep)) < phistep)
0586     sectors.push_back(7);
0587   if (fabs(deltaPhi(phigp, 7 * phistep)) < phistep)
0588     sectors.push_back(8);
0589   if (fabs(deltaPhi(phigp, 8 * phistep)) < phistep)
0590     sectors.push_back(9);
0591   if (fabs(deltaPhi(phigp, 9 * phistep)) < phistep) {
0592     sectors.push_back(10);
0593     if (station == 4)
0594       sectors.push_back(14);
0595   }
0596   if (fabs(deltaPhi(phigp, 10 * phistep)) < phistep)
0597     sectors.push_back(11);
0598   if (fabs(deltaPhi(phigp, 11 * phistep)) < phistep)
0599     sectors.push_back(12);
0600
0601   LogTrace(category_) << "DT position to dets" << endl;
0602   LogTrace(category_) << "number of sectors to consider: " << sectors.size() << endl;
0603   LogTrace(category_) << "station: " << station << " wheels: " << minwheel << " " << maxwheel << endl;
0604
0605   vector<const GeomDet*> result;
0606   if (station > 4 || station < 1)
0607     return result;
0608   if (minwheel > 2 || maxwheel < -2)
0609     return result;
0610
0611   for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
0612     for (int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
0613       DTChamberId chamberid(iwheel, station, (*isector));
0614       result.push_back(theService->trackingGeometry()->idToDet(chamberid));
0615     }
0616   }
0617
0618   LogTrace(category_) << "number of GeomDets for this track: " << result.size() << endl;
0619
0620   return result;
0621 }
0622
0623 //

0624 // GeomDets along the track in CSC

0625 //

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

0628   int endcap = 0;
0629   if (gp.z() > 0) {
0630     endcap = 1;
0631   } else {
0632     endcap = 2;
0633   }
0634
0635   // determine the csc station and range of rings

0636   int station = 5;
0637
0638   // check all rings in a station

0639   if (fabs(gp.z()) > 1000. && fabs(gp.z()) < 1055.0) {
0640     station = 4;
0641   } else if (fabs(gp.z()) > 910.0 && fabs(gp.z()) < 965.0) {
0642     station = 3;
0643   } else if (fabs(gp.z()) > 800.0 && fabs(gp.z()) < 860.0) {
0644     station = 2;
0645   } else if (fabs(gp.z()) > 570.0 && fabs(gp.z()) < 730.0) {
0646     station = 1;
0647   }
0648
0649   vector<int> sectors;
0650
0651   float phistep1 = M_PI / 18.;  //for all the rings except first rings for stations > 1

0652   float phistep2 = M_PI / 9.;
0653   float phigp = (float)gp.barePhi();
0654
0655   int ring = -1;
0656
0657   // determine the ring

0658   if (station == 1) {
0659     //FIX ME!!!

0660     if (gp.perp() > 100 && gp.perp() < 270)
0661       ring = 1;
0662     else if (gp.perp() > 270 && gp.perp() < 450)
0663       ring = 2;
0664     else if (gp.perp() > 450 && gp.perp() < 695)
0665       ring = 3;
0666     else if (gp.perp() > 100 && gp.perp() < 270)
0667       ring = 4;
0668
0669   } else if (station == 2) {
0670     if (gp.perp() > 140 && gp.perp() < 350)
0671       ring = 1;
0672     else if (gp.perp() > 350 && gp.perp() < 700)
0673       ring = 2;
0674
0675   } else if (station == 3) {
0676     if (gp.perp() > 160 && gp.perp() < 350)
0677       ring = 1;
0678     else if (gp.perp() > 350 && gp.perp() < 700)
0679       ring = 2;
0680
0681   } else if (station == 4) {
0682     if (gp.perp() > 175 && gp.perp() < 350)
0683       ring = 1;
0684     else if (gp.perp() > 350 && gp.perp() < 700)
0685       ring = 2;
0686   }
0687
0688   if (station > 1 && ring == 1) {
0689     // we have 18 sectors in that case

0690     for (int i = 0; i < 18; i++) {
0691       if (fabs(deltaPhi(phigp, i * phistep2)) < phistep2)
0692         sectors.push_back(i + 1);
0693     }
0694
0695   } else {
0696     // we have 36 sectors in that case

0697     for (int i = 0; i < 36; i++) {
0698       if (fabs(deltaPhi(phigp, i * phistep1)) < phistep1)
0699         sectors.push_back(i + 1);
0700     }
0701   }
0702
0703   LogTrace(category_) << "CSC position to dets" << endl;
0704   LogTrace(category_) << "ring: " << ring << endl;
0705   LogTrace(category_) << "endcap: " << endcap << endl;
0706   LogTrace(category_) << "station: " << station << endl;
0707   LogTrace(category_) << "CSC number of sectors to consider: " << sectors.size() << endl;
0708
0709   // check exceptional cases

0710   vector<const GeomDet*> result;
0711   if (station > 4 || station < 1)
0712     return result;
0713   if (endcap == 0)
0714     return result;
0715   if (ring == -1)
0716     return result;
0717
0718   int minlayer = 1;
0719   int maxlayer = 6;
0720
0721   for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
0722     for (int ilayer = minlayer; ilayer != maxlayer + 1; ++ilayer) {
0723       CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
0724       result.push_back(theService->trackingGeometry()->idToDet(cscid));
0725     }
0726   }
0727
0728   return result;
0729 }
0730
0731 //

0732 //Set class members

0733 //

0734 void MuonShowerInformationFiller::fillHitsByStation(const reco::Muon& muon) {
0735   reco::TrackRef track;
0736   if (muon.isGlobalMuon())
0737     track = muon.globalTrack();
0738   else if (muon.isStandAloneMuon())
0739     track = muon.outerTrack();
0740   else
0741     return;
0742
0743   // split 1D rechits by station

0744   vector<MuonRecHitContainer> muonRecHits(4);
0745
0746   // split rechits from segs by station

0747   vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);
0748
0749   // get vector of GeomDets compatible with a track

0750   vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);
0751
0752   // for special cases: CSC station 1

0753   MuonRecHitContainer tmpCSC1;
0754   bool dtOverlapToCheck = false;
0755   bool cscOverlapToCheck = false;
0756
0757   for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++) {
0758     // get det id

0759     DetId geoId = (*igd)->geographicalId();
0760
0761     // skip tracker hits

0762     if (geoId.det() != DetId::Muon)
0763       continue;
0764
0765     // DT

0766     if (geoId.subdetId() == MuonSubdetId::DT) {
0767       // get station

0768       DTChamberId detid(geoId.rawId());
0769       int station = detid.station();
0770       int wheel = detid.wheel();
0771
0772       // get rechits from segments per station

0773       TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp =
0774           hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
0775       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
0776       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
0777
0778       for (; hits_begin != hits_end; ++hits_begin) {
0779         muonCorrelatedHits.at(station - 1).push_back(*hits_begin);
0780       }
0781
0782       //check overlap certain wheels and stations

0783       if (abs(wheel) == 2 && station != 4 && station != 1)
0784         dtOverlapToCheck = true;
0785
0786       // loop over all superlayers of a DT chamber

0787       for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1;
0788            ++isuperlayer) {
0789         // loop over all layers inside the superlayer

0790         for (int ilayer = DTChamberId::minLayerId; ilayer != DTChamberId::maxLayerId + 1; ++ilayer) {
0791           DTLayerId lid(detid, isuperlayer, ilayer);
0792           DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
0793           for (DTRecHitCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
0794             vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
0795             for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end();
0796                  ++irechit) {
0797               muonRecHits.at(station - 1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &**irechit));
0798             }
0799           }
0800         }
0801       }
0802     } else if (geoId.subdetId() == MuonSubdetId::CSC) {
0803       // get station

0804       CSCDetId did(geoId.rawId());
0805       int station = did.station();
0806       int ring = did.ring();
0807
0808       //get rechits from segments by station

0809       TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp =
0810           hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
0811       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
0812       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
0813
0814       for (; hits_begin != hits_end; ++hits_begin) {
0815         muonCorrelatedHits.at(station - 1).push_back(*hits_begin);
0816       }
0817
0818       if ((station == 1 && ring == 3) && dtOverlapToCheck)
0819         cscOverlapToCheck = true;
0820
0821       // split 1D rechits by station

0822       CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
0823       for (CSCRecHit2DCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
0824         if (!cscOverlapToCheck) {
0825           muonRecHits.at(station - 1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &*rechit));
0826         } else {
0827           tmpCSC1.push_back(MuonTransientTrackingRecHit::specificBuild((&**igd), &*rechit));
0828
0829           //sort by perp, then insert to appropriate container

0830           MuonRecHitContainer temp = findPerpCluster(tmpCSC1);
0831           if (temp.empty())
0832             continue;
0833
0834           float center;
0835           if (temp.size() > 1) {
0836             center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp()) / 2.;
0837           } else {
0838             center = temp.front()->globalPosition().perp();
0839           }
0840           temp.clear();
0841
0842           if (center > 550.) {
0843             muonRecHits.at(2).insert(muonRecHits.at(2).end(), tmpCSC1.begin(), tmpCSC1.end());
0844           } else {
0845             muonRecHits.at(1).insert(muonRecHits.at(1).end(), tmpCSC1.begin(), tmpCSC1.end());
0846           }
0847           tmpCSC1.clear();
0848         }
0849       }
0850     } else if (geoId.subdetId() == MuonSubdetId::RPC) {
0851       LogTrace(category_) << "Wrong subdet id" << endl;
0852     }
0853   }  //loop over GeomDets compatible with a track

0854
0855   // calculate number of all and correlated hits

0856   for (int stat = 0; stat < 4; stat++) {
0857     theCorrelatedStationHits[stat] = muonCorrelatedHits.at(stat).size();
0858     theAllStationHits[stat] = muonRecHits[stat].size();
0859   }
0860   LogTrace(category_) << "Hits used by the segments, by station " << theCorrelatedStationHits.at(0) << " "
0861                       << theCorrelatedStationHits.at(1) << " " << theCorrelatedStationHits.at(2) << " "
0862                       << theCorrelatedStationHits.at(3) << endl;
0863
0864   LogTrace(category_) << "All DT 1D/CSC 2D  hits, by station " << theAllStationHits.at(0) << " "
0865                       << theAllStationHits.at(1) << " " << theAllStationHits.at(2) << " " << theAllStationHits.at(3)
0866                       << endl;
0867
0868   //station shower sizes

0869   MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiBest;
0870   TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
0871   // send station hits to the clustering algorithm

0872   for (int stat = 0; stat != 4; stat++) {
0873     const size_t nhit = muonRecHits[stat].size();
0874     if (nhit < 2)
0875       continue;  // Require at least 2 hits

0876     muonRecHitsPhiBest.clear();
0877     muonRecHitsPhiBest.reserve(nhit);
0878
0879     // Cluster seeds by global position phi. Best cluster is chosen to give greatest dphi

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

0881     muonRecHits[stat] = sort_all_indexed(muonRecHits[stat], std::less(), PhiTransform());
0882
0883     // Search for gaps (complexity = N)

0884     std::vector<size_t> clUppers;
0885     for (size_t ihit = 0; ihit < nhit; ++ihit) {
0886       const size_t jhit = (ihit + 1) % nhit;
0887       const double phi1 = muonRecHits[stat].at(ihit)->globalPosition().barePhi();
0888       const double phi2 = muonRecHits[stat].at(jhit)->globalPosition().barePhi();
0889
0890       const double dphi = std::abs(deltaPhi(phi1, phi2));
0891       if (dphi >= 0.05)
0892         clUppers.push_back(ihit);
0893     }
0894
0895     //station shower sizes

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

0899       const double refPhi = muonRecHits[stat].at(0)->globalPosition().barePhi();
0900       double dphilo = 0, dphihi = 0;
0901       for (auto& hit : muonRecHits[stat]) {
0902         muonRecHitsPhiBest.push_back(hit);
0903         const double phi = hit->globalPosition().barePhi();
0904         dphilo = std::min(dphilo, deltaPhi(refPhi, phi));
0905         dphihi = std::max(dphihi, deltaPhi(refPhi, phi));
0906       }
0907       dphimax = std::abs(dphihi + dphilo);
0908     } else {
0909       // Loop over gaps to find the one with maximum dphi(begin, end)

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

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

0914         const size_t upper = *icl;
0915         // lower bound is +1 of preceeding upper bound

0916         const auto prevCl = (icl == clUppers.begin()) ? clUppers.end() : icl;
0917         const size_t lower = (*(prevCl - 1) + 1) % nhit;
0918
0919         // At least two hit for a cluster

0920         if (upper == lower)
0921           continue;
0922
0923         const double phi1 = muonRecHits[stat].at(upper)->globalPosition().barePhi();
0924         const double phi2 = muonRecHits[stat].at(lower)->globalPosition().barePhi();
0925
0926         const double dphi = std::abs(deltaPhi(phi1, phi2));
0927         if (dphimax < dphi) {
0928           dphimax = dphi;
0929           bestUpper = upper;
0930           bestLower = lower;
0931         }
0932       }
0933
0934       if (bestUpper > bestLower) {
0935         muonRecHitsPhiBest.reserve(bestUpper - bestLower + 1);
0936         std::copy(muonRecHits[stat].begin() + bestLower,
0937                   muonRecHits[stat].begin() + bestUpper + 1,
0938                   std::back_inserter(muonRecHitsPhiBest));
0939       } else if (bestUpper < bestLower) {
0940         muonRecHitsPhiBest.reserve(bestUpper + (nhit - bestLower) + 1);
0941         std::copy(muonRecHits[stat].begin(),
0942                   muonRecHits[stat].begin() + bestUpper + 1,
0943                   std::back_inserter(muonRecHitsPhiBest));
0944         std::copy(
0945             muonRecHits[stat].begin() + bestLower, muonRecHits[stat].end(), std::back_inserter(muonRecHitsPhiBest));
0946       }
0947     }
0948
0949     //fill showerTs

0950     if (!muonRecHitsPhiBest.empty()) {
0951       muonRecHits[stat] = sort_all_indexed(muonRecHitsPhiBest, std::less(), AbsMagTransform());
0952       GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
0953       theStationShowerTSize.at(stat) = refpoint.mag() * dphimax;
0954     }
0955
0956     //for theta

0957     if (!muonCorrelatedHits.at(stat).empty()) {
0958       float dthetamax = 0;
0959       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iseed = muonCorrelatedHits.at(stat).begin();
0960            iseed != muonCorrelatedHits.at(stat).end();
0961            ++iseed) {
0962         if (!(*iseed)->isValid())
0963           continue;
0964         GlobalPoint refpoint = (*iseed)->globalPosition();  //starting from the one with smallest value of phi

0965         muonRecHitsThetaTemp.clear();
0966         muonRecHitsThetaTemp = findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
0967         if (muonRecHitsThetaTemp.size() > 1) {
0968           float dtheta = fabs((float)muonRecHitsThetaTemp.back()->globalPosition().theta() -
0969                               (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
0970           if (dtheta > dthetamax) {
0971             dthetamax = dtheta;
0972             muonRecHitsThetaBest = muonRecHitsThetaTemp;
0973           }
0974         }  //at least two hits

0975       }    //loop over seeds

0976     }      //not empty container2

0977
0978     //fill deltaRs

0979     if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
0980       theStationShowerDeltaR.at(stat) = sqrt(pow(deltaPhi(muonRecHitsPhiBest.front()->globalPosition().barePhi(),
0981                                                           muonRecHitsPhiBest.back()->globalPosition().barePhi()),
0982                                                  2) +
0983                                              pow(muonRecHitsThetaBest.front()->globalPosition().theta() -
0984                                                      muonRecHitsThetaBest.back()->globalPosition().theta(),
0985                                                  2));
0986
0987   }  //loop over station

0988
0989   LogTrace(category_) << "deltaR around a track containing all the station hits, by station "
0990                       << theStationShowerDeltaR.at(0) << " " << theStationShowerDeltaR.at(1) << " "
0991                       << theStationShowerDeltaR.at(2) << " " << theStationShowerDeltaR.at(3) << endl;
0992
0993   LogTrace(category_) << "Transverse cluster size, by station " << theStationShowerTSize.at(0) << " "
0994                       << theStationShowerTSize.at(1) << " " << theStationShowerTSize.at(2) << " "
0995                       << theStationShowerTSize.at(3) << endl;
0996
0997   return;
0998 }