Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:55 

0001 /**
0002  *  Class: GlobalMuonRefitter
0003  *
0004  *  Description:
0005  *
0006  *
0007  *
0008  *  Authors :
0009  *  P. Traczyk, SINS Warsaw
0010  *
0011  *  \modified by C. Calabria, INFN & Universita Bari
0012  *  \modified by D. Nash, Northeastern University
0013  *  \modified by C. Caputo, UCLouvain
0014  *
0015  **/
0016 
0017 #include "RecoMuon/GlobalTrackingTools/interface/GlobalMuonRefitter.h"
0018 
0019 //---------------
0020 // C++ Headers --
0021 //---------------
0022 
0023 #include <iostream>
0024 #include <iomanip>
0025 #include <algorithm>
0026 
0027 //-------------------------------
0028 // Collaborating Class Headers --
0029 //-------------------------------
0030 
0031 #include "FWCore/Framework/interface/Event.h"
0032 
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0035 #include "TrackingTools/TrackFitters/interface/RecHitLessByDet.h"
0036 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0037 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0038 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0039 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0040 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0041 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0042 
0043 #include "DataFormats/DetId/interface/DetId.h"
0044 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0045 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0046 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0047 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0048 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0049 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0050 #include "Geometry/DTGeometry/interface/DTLayer.h"
0051 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
0052 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0053 #include <DataFormats/GEMRecHit/interface/GEMRecHit.h>
0054 #include <DataFormats/GEMRecHit/interface/ME0Segment.h>
0055 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0056 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0057 
0058 #include "DataFormats/TrackReco/interface/Track.h"
0059 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0060 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0061 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0062 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0063 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0064 #include "RecoMuon/TrackingTools/interface/MuonCandidate.h"
0065 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0066 #include "RecoMuon/GlobalTrackingTools/interface/DynamicTruncation.h"
0067 
0068 using namespace std;
0069 using namespace edm;
0070 
0071 //----------------
0072 // Constructors --
0073 //----------------
0074 
0075 GlobalMuonRefitter::GlobalMuonRefitter(const edm::ParameterSet& par,
0076                                        const MuonServiceProxy* service,
0077                                        edm::ConsumesCollector& iC)
0078     : theCosmicFlag(par.getParameter<bool>("PropDirForCosmics")),
0079       theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
0080       theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
0081       theGEMRecHitLabel(par.getParameter<InputTag>("GEMRecHitLabel")),
0082       theME0RecHitLabel(par.getParameter<InputTag>("ME0RecHitLabel")),
0083       theService(service),
0084       theDynamicTruncationConfig(iC) {
0085   theCategory = par.getUntrackedParameter<string>("Category", "Muon|RecoMuon|GlobalMuon|GlobalMuonRefitter");
0086 
0087   theHitThreshold = par.getParameter<int>("HitThreshold");
0088   theDTChi2Cut = par.getParameter<double>("Chi2CutDT");
0089   theCSCChi2Cut = par.getParameter<double>("Chi2CutCSC");
0090   theRPCChi2Cut = par.getParameter<double>("Chi2CutRPC");
0091   theGEMChi2Cut = par.getParameter<double>("Chi2CutGEM");
0092   theME0Chi2Cut = par.getParameter<double>("Chi2CutME0");
0093 
0094   // Refit direction
0095   string refitDirectionName = par.getParameter<string>("RefitDirection");
0096 
0097   if (refitDirectionName == "insideOut")
0098     theRefitDirection = insideOut;
0099   else if (refitDirectionName == "outsideIn")
0100     theRefitDirection = outsideIn;
0101   else
0102     throw cms::Exception("TrackTransformer constructor")
0103         << "Wrong refit direction chosen in TrackTransformer ParameterSet"
0104         << "\n"
0105         << "Possible choices are:"
0106         << "\n"
0107         << "RefitDirection = insideOut or RefitDirection = outsideIn";
0108 
0109   theFitterToken = iC.esConsumes(edm::ESInputTag("", par.getParameter<string>("Fitter")));
0110   thePropagatorName = par.getParameter<string>("Propagator");
0111 
0112   theSkipStation = par.getParameter<int>("SkipStation");
0113   theTrackerSkipSystem = par.getParameter<int>("TrackerSkipSystem");
0114   theTrackerSkipSection = par.getParameter<int>("TrackerSkipSection");  //layer, wheel, or disk depending on the system
0115 
0116   theTrackerRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", par.getParameter<string>("TrackerRecHitBuilder")));
0117   theMuonRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", par.getParameter<string>("MuonRecHitBuilder")));
0118 
0119   theRPCInTheFit = par.getParameter<bool>("RefitRPCHits");
0120 
0121   theDYTthrs = par.getParameter<std::vector<int> >("DYTthrs");
0122   theDYTselector = par.getParameter<int>("DYTselector");
0123   theDYTupdator = par.getParameter<bool>("DYTupdator");
0124   theDYTuseAPE = par.getParameter<bool>("DYTuseAPE");
0125   theDYTParThrsMode = par.getParameter<bool>("DYTuseThrsParametrization");
0126   if (theDYTParThrsMode)
0127     theDYTthrsParameters = par.getParameter<edm::ParameterSet>("DYTthrsParameters");
0128   dytInfo = new reco::DYTInfo();
0129 
0130   if (par.existsAs<double>("RescaleErrorFactor")) {
0131     theRescaleErrorFactor = par.getParameter<double>("RescaleErrorFactor");
0132     edm::LogWarning("GlobalMuonRefitter") << "using error rescale factor " << theRescaleErrorFactor;
0133   } else
0134     theRescaleErrorFactor = 1000.;
0135 
0136   theCacheId_TRH = 0;
0137   theDTRecHitToken = iC.consumes<DTRecHitCollection>(theDTRecHitLabel);
0138   theCSCRecHitToken = iC.consumes<CSCRecHit2DCollection>(theCSCRecHitLabel);
0139   theGEMRecHitToken = iC.consumes<GEMRecHitCollection>(theGEMRecHitLabel);
0140   theME0RecHitToken = iC.consumes<ME0SegmentCollection>(theME0RecHitLabel);
0141   CSCSegmentsToken = iC.consumes<CSCSegmentCollection>(InputTag("cscSegments"));
0142   all4DSegmentsToken = iC.consumes<DTRecSegment4DCollection>(InputTag("dt4DSegments"));
0143 }
0144 
0145 //--------------
0146 // Destructor --
0147 //--------------
0148 
0149 GlobalMuonRefitter::~GlobalMuonRefitter() { delete dytInfo; }
0150 
0151 //
0152 // set Event
0153 //
0154 void GlobalMuonRefitter::setEvent(const edm::Event& event) {
0155   event.getByToken(theDTRecHitToken, theDTRecHits);
0156   event.getByToken(theCSCRecHitToken, theCSCRecHits);
0157   event.getByToken(theGEMRecHitToken, theGEMRecHits);
0158   event.getByToken(theME0RecHitToken, theME0RecHits);
0159   event.getByToken(CSCSegmentsToken, CSCSegments);
0160   event.getByToken(all4DSegmentsToken, all4DSegments);
0161 }
0162 
0163 void GlobalMuonRefitter::setServices(const EventSetup& setup) {
0164   theFitter = setup.getData(theFitterToken).clone();
0165 
0166   // Transient Rechit Builders
0167   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
0168   if (newCacheId_TRH != theCacheId_TRH) {
0169     LogDebug(theCategory) << "TransientRecHitRecord changed!";
0170     theTrackerRecHitBuilder = &setup.getData(theTrackerRecHitBuilderToken);
0171     theMuonRecHitBuilder = &setup.getData(theMuonRecHitBuilderToken);
0172     hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder)->cloner();
0173   }
0174   theFitter->setHitCloner(&hitCloner);
0175   theEventSetup = &setup;
0176 }
0177 
0178 //
0179 // build a combined tracker-muon trajectory
0180 //
0181 vector<Trajectory> GlobalMuonRefitter::refit(const reco::Track& globalTrack,
0182                                              const int theMuonHitsOption,
0183                                              const TrackerTopology* tTopo) const {
0184   LogTrace(theCategory) << " *** GlobalMuonRefitter *** option " << theMuonHitsOption << endl;
0185 
0186   ConstRecHitContainer allRecHitsTemp;  // all muon rechits temp
0187 
0188   reco::TransientTrack track(globalTrack, &*(theService->magneticField()), theService->trackingGeometry());
0189 
0190   auto tkbuilder = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder);
0191 
0192   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit)
0193     if ((*hit)->isValid()) {
0194       if ((*hit)->geographicalId().det() == DetId::Tracker)
0195         allRecHitsTemp.push_back((**hit).cloneForFit(*tkbuilder->geometry()->idToDet((**hit).geographicalId())));
0196       else if ((*hit)->geographicalId().det() == DetId::Muon) {
0197         if ((*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit) {
0198           LogTrace(theCategory) << "RPC Rec Hit discarged";
0199           continue;
0200         }
0201         allRecHitsTemp.push_back(theMuonRecHitBuilder->build(&**hit));
0202       }
0203     }
0204   vector<Trajectory> refitted = refit(globalTrack, track, allRecHitsTemp, theMuonHitsOption, tTopo);
0205   return refitted;
0206 }
0207 
0208 //
0209 // build a combined tracker-muon trajectory
0210 //
0211 vector<Trajectory> GlobalMuonRefitter::refit(const reco::Track& globalTrack,
0212                                              const reco::TransientTrack track,
0213                                              const TransientTrackingRecHit::ConstRecHitContainer& allRecHitsTemp,
0214                                              const int theMuonHitsOption,
0215                                              const TrackerTopology* tTopo) const {
0216   // MuonHitsOption: 0 - tracker only
0217   //                 1 - include all muon hits
0218   //                 2 - include only first muon hit(s)
0219   //                 3 - include only selected muon hits
0220   //                 4 - redo pattern recognition with dynamic truncation
0221 
0222   vector<int> stationHits(4, 0);
0223   map<DetId, int> hitMap;
0224 
0225   ConstRecHitContainer allRecHits;       // all muon rechits
0226   ConstRecHitContainer fmsRecHits;       // only first muon rechits
0227   ConstRecHitContainer selectedRecHits;  // selected muon rechits
0228   ConstRecHitContainer DYTRecHits;       // rec hits from dynamic truncation algorithm
0229 
0230   LogTrace(theCategory) << " *** GlobalMuonRefitter *** option " << theMuonHitsOption << endl;
0231   LogTrace(theCategory) << " Track momentum before refit: " << globalTrack.pt() << endl;
0232   LogTrace(theCategory) << " Hits size before : " << allRecHitsTemp.size() << endl;
0233 
0234   allRecHits = getRidOfSelectStationHits(allRecHitsTemp, tTopo);
0235   //    printHits(allRecHits);
0236   LogTrace(theCategory) << " Hits size: " << allRecHits.size() << endl;
0237 
0238   vector<Trajectory> outputTraj;
0239 
0240   if ((theMuonHitsOption == 1) || (theMuonHitsOption == 3) || (theMuonHitsOption == 4)) {
0241     // refit the full track with all muon hits
0242     vector<Trajectory> globalTraj = transform(globalTrack, track, allRecHits);
0243 
0244     if (globalTraj.empty()) {
0245       LogTrace(theCategory) << "No trajectory from the TrackTransformer!" << endl;
0246       return vector<Trajectory>();
0247     }
0248 
0249     LogTrace(theCategory) << " Initial trajectory state: "
0250                           << globalTraj.front().lastMeasurement().updatedState().freeState()->parameters() << endl;
0251 
0252     if (theMuonHitsOption == 1)
0253       outputTraj.push_back(globalTraj.front());
0254 
0255     if (theMuonHitsOption == 3) {
0256       checkMuonHits(globalTrack, allRecHits, hitMap);
0257       selectedRecHits = selectMuonHits(globalTraj.front(), hitMap);
0258       LogTrace(theCategory) << " Selected hits size: " << selectedRecHits.size() << endl;
0259       outputTraj = transform(globalTrack, track, selectedRecHits);
0260     }
0261 
0262     if (theMuonHitsOption == 4) {
0263       //
0264       // DYT 2.0
0265       //
0266       DynamicTruncation dytRefit(theDynamicTruncationConfig, *theEventSetup, *theService);
0267       dytRefit.setProd(all4DSegments, CSCSegments);
0268       dytRefit.setSelector(theDYTselector);
0269       dytRefit.setThr(theDYTthrs);
0270       dytRefit.setUpdateState(theDYTupdator);
0271       dytRefit.setUseAPE(theDYTuseAPE);
0272       if (theDYTParThrsMode) {
0273         dytRefit.setParThrsMode(theDYTParThrsMode);
0274         dytRefit.setThrsMap(theDYTthrsParameters);
0275         dytRefit.setRecoP(globalTrack.p());
0276         dytRefit.setRecoEta(globalTrack.eta());
0277       }
0278       DYTRecHits = dytRefit.filter(globalTraj.front());
0279       dytInfo->CopyFrom(dytRefit.getDYTInfo());
0280       if ((DYTRecHits.size() > 1) &&
0281           (DYTRecHits.front()->globalPosition().mag() > DYTRecHits.back()->globalPosition().mag()))
0282         stable_sort(DYTRecHits.begin(), DYTRecHits.end(), RecHitLessByDet(alongMomentum));
0283       outputTraj = transform(globalTrack, track, DYTRecHits);
0284     }
0285 
0286   } else if (theMuonHitsOption == 2) {
0287     getFirstHits(globalTrack, allRecHits, fmsRecHits);
0288     outputTraj = transform(globalTrack, track, fmsRecHits);
0289   }
0290 
0291   if (!outputTraj.empty()) {
0292     LogTrace(theCategory) << "Refitted pt: "
0293                           << outputTraj.front().firstMeasurement().updatedState().globalParameters().momentum().perp()
0294                           << endl;
0295     return outputTraj;
0296   } else {
0297     LogTrace(theCategory) << "No refitted Tracks... " << endl;
0298     return vector<Trajectory>();
0299   }
0300 }
0301 
0302 //
0303 //
0304 //
0305 void GlobalMuonRefitter::checkMuonHits(const reco::Track& muon,
0306                                        ConstRecHitContainer& all,
0307                                        map<DetId, int>& hitMap) const {
0308   LogTrace(theCategory) << " GlobalMuonRefitter::checkMuonHits " << endl;
0309 
0310   float coneSize = 20.0;
0311 
0312   // loop through all muon hits and calculate the maximum # of hits in each chamber
0313   for (ConstRecHitContainer::const_iterator imrh = all.begin(); imrh != all.end(); imrh++) {
0314     if ((*imrh != nullptr) && !(*imrh)->isValid())
0315       continue;
0316 
0317     int detRecHits = 0;
0318     MuonRecHitContainer dRecHits;
0319 
0320     DetId id = (*imrh)->geographicalId();
0321     DetId chamberId;
0322 
0323     // Skip tracker hits
0324     if (id.det() != DetId::Muon)
0325       continue;
0326 
0327     if (id.subdetId() == MuonSubdetId::DT) {
0328       DTChamberId did(id.rawId());
0329       chamberId = did;
0330 
0331       if ((*imrh)->dimension() > 1) {
0332         std::vector<const TrackingRecHit*> hits2d = (*imrh)->recHits();
0333         for (std::vector<const TrackingRecHit*>::const_iterator hit2d = hits2d.begin(); hit2d != hits2d.end();
0334              hit2d++) {
0335           if ((*hit2d)->dimension() > 1) {
0336             std::vector<const TrackingRecHit*> hits1d = (*hit2d)->recHits();
0337             for (std::vector<const TrackingRecHit*>::const_iterator hit1d = hits1d.begin(); hit1d != hits1d.end();
0338                  hit1d++) {
0339               DetId id1 = (*hit1d)->geographicalId();
0340               DTLayerId lid(id1.rawId());
0341               // Get the 1d DT RechHits from this layer
0342               DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
0343               int layerHits = 0;
0344               for (DTRecHitCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0345                 double rhitDistance = fabs(ir->localPosition().x() - (**hit1d).localPosition().x());
0346                 if (rhitDistance < coneSize)
0347                   layerHits++;
0348                 LogTrace(theCategory) << "       " << (ir)->localPosition() << "  " << (**hit1d).localPosition()
0349                                       << " Distance: " << rhitDistance << " recHits: " << layerHits
0350                                       << "  SL: " << lid.superLayer() << endl;
0351               }
0352               if (layerHits > detRecHits)
0353                 detRecHits = layerHits;
0354             }
0355           } else {
0356             DTLayerId lid(id.rawId());
0357             // Get the 1d DT RechHits from this layer
0358             DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
0359             for (DTRecHitCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0360               double rhitDistance = fabs(ir->localPosition().x() - (**imrh).localPosition().x());
0361               if (rhitDistance < coneSize)
0362                 detRecHits++;
0363               LogTrace(theCategory) << "       " << (ir)->localPosition() << "  " << (**imrh).localPosition()
0364                                     << " Distance: " << rhitDistance << " recHits: " << detRecHits << endl;
0365             }
0366           }
0367         }
0368 
0369       } else {
0370         DTLayerId lid(id.rawId());
0371 
0372         // Get the 1d DT RechHits from this layer
0373         DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
0374 
0375         for (DTRecHitCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0376           double rhitDistance = fabs(ir->localPosition().x() - (**imrh).localPosition().x());
0377           if (rhitDistance < coneSize)
0378             detRecHits++;
0379           LogTrace(theCategory) << "       " << (ir)->localPosition() << "  " << (**imrh).localPosition()
0380                                 << " Distance: " << rhitDistance << " recHits: " << detRecHits << endl;
0381         }
0382       }
0383     }  // end of if DT
0384     else if (id.subdetId() == MuonSubdetId::CSC) {
0385       CSCDetId did(id.rawId());
0386       chamberId = did.chamberId();
0387 
0388       if ((*imrh)->recHits().size() > 1) {
0389         std::vector<const TrackingRecHit*> hits2d = (*imrh)->recHits();
0390         for (std::vector<const TrackingRecHit*>::const_iterator hit2d = hits2d.begin(); hit2d != hits2d.end();
0391              hit2d++) {
0392           DetId id1 = (*hit2d)->geographicalId();
0393           CSCDetId lid(id1.rawId());
0394 
0395           // Get the CSC Rechits from this layer
0396           CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(lid);
0397           int layerHits = 0;
0398 
0399           for (CSCRecHit2DCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0400             double rhitDistance = (ir->localPosition() - (**hit2d).localPosition()).mag();
0401             if (rhitDistance < coneSize)
0402               layerHits++;
0403             LogTrace(theCategory) << ir->localPosition() << "  " << (**hit2d).localPosition()
0404                                   << " Distance: " << rhitDistance << " recHits: " << layerHits << endl;
0405           }
0406           if (layerHits > detRecHits)
0407             detRecHits = layerHits;
0408         }
0409       } else {
0410         // Get the CSC Rechits from this layer
0411         CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
0412 
0413         for (CSCRecHit2DCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0414           double rhitDistance = (ir->localPosition() - (**imrh).localPosition()).mag();
0415           if (rhitDistance < coneSize)
0416             detRecHits++;
0417           LogTrace(theCategory) << ir->localPosition() << "  " << (**imrh).localPosition()
0418                                 << " Distance: " << rhitDistance << " recHits: " << detRecHits << endl;
0419         }
0420       }
0421     }  //end of CSC if
0422     else if (id.subdetId() == MuonSubdetId::GEM) {
0423       GEMDetId did(id.rawId());
0424       chamberId = did.chamberId();
0425 
0426       if ((*imrh)->recHits().size() > 1) {
0427         std::vector<const TrackingRecHit*> hits2d = (*imrh)->recHits();
0428         for (std::vector<const TrackingRecHit*>::const_iterator hit2d = hits2d.begin(); hit2d != hits2d.end();
0429              hit2d++) {
0430           DetId id1 = (*hit2d)->geographicalId();
0431           GEMDetId lid(id1.rawId());
0432 
0433           // Get the GEM Rechits from this layer
0434           GEMRecHitCollection::range dRecHits = theGEMRecHits->get(lid);
0435           int layerHits = 0;
0436 
0437           for (GEMRecHitCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0438             double rhitDistance = (ir->localPosition() - (**hit2d).localPosition()).mag();
0439             if (rhitDistance < coneSize)
0440               layerHits++;
0441             LogTrace(theCategory) << ir->localPosition() << "  " << (**hit2d).localPosition()
0442                                   << " Distance: " << rhitDistance << " recHits: " << layerHits << endl;
0443           }
0444           if (layerHits > detRecHits)
0445             detRecHits = layerHits;
0446         }
0447       } else {
0448         // Get the GEM Rechits from this layer
0449         GEMRecHitCollection::range dRecHits = theGEMRecHits->get(did);
0450 
0451         for (GEMRecHitCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0452           double rhitDistance = (ir->localPosition() - (**imrh).localPosition()).mag();
0453           if (rhitDistance < coneSize)
0454             detRecHits++;
0455           LogTrace(theCategory) << ir->localPosition() << "  " << (**imrh).localPosition()
0456                                 << " Distance: " << rhitDistance << " recHits: " << detRecHits << endl;
0457         }
0458       }
0459     }  //end of GEM if
0460     else if (id.subdetId() == MuonSubdetId::ME0) {
0461       ME0DetId did(id.rawId());
0462       chamberId = did.chamberId();
0463 
0464       if ((*imrh)->recHits().size() > 1) {
0465         std::vector<const TrackingRecHit*> hits2d = (*imrh)->recHits();
0466         for (std::vector<const TrackingRecHit*>::const_iterator hit2d = hits2d.begin(); hit2d != hits2d.end();
0467              hit2d++) {
0468           DetId id1 = (*hit2d)->geographicalId();
0469           ME0DetId lid(id1.rawId());
0470 
0471           // Get the ME0 Rechits from this layer
0472           ME0SegmentCollection::range dRecHits = theME0RecHits->get(lid);
0473           int layerHits = 0;
0474 
0475           for (ME0SegmentCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0476             double rhitDistance = (ir->localPosition() - (**hit2d).localPosition()).mag();
0477             if (rhitDistance < coneSize)
0478               layerHits++;
0479             LogTrace(theCategory) << ir->localPosition() << "  " << (**hit2d).localPosition()
0480                                   << " Distance: " << rhitDistance << " recHits: " << layerHits << endl;
0481           }
0482           if (layerHits > detRecHits)
0483             detRecHits = layerHits;
0484         }
0485       } else {
0486         // Get the ME0 Rechits from this layer
0487         ME0SegmentCollection::range dRecHits = theME0RecHits->get(did);
0488 
0489         for (ME0SegmentCollection::const_iterator ir = dRecHits.first; ir != dRecHits.second; ir++) {
0490           double rhitDistance = (ir->localPosition() - (**imrh).localPosition()).mag();
0491           if (rhitDistance < coneSize)
0492             detRecHits++;
0493           LogTrace(theCategory) << ir->localPosition() << "  " << (**imrh).localPosition()
0494                                 << " Distance: " << rhitDistance << " recHits: " << detRecHits << endl;
0495         }
0496       }
0497     }  //end of ME0 if
0498     else {
0499       if (id.subdetId() != MuonSubdetId::RPC)
0500         LogError(theCategory) << " Wrong Hit Type ";
0501       continue;
0502     }
0503 
0504     map<DetId, int>::iterator imap = hitMap.find(chamberId);
0505     if (imap != hitMap.end()) {
0506       if (detRecHits > imap->second)
0507         imap->second = detRecHits;
0508     } else
0509       hitMap[chamberId] = detRecHits;
0510 
0511   }  // end of loop over muon rechits
0512 
0513   for (map<DetId, int>::iterator imap = hitMap.begin(); imap != hitMap.end(); imap++)
0514     LogTrace(theCategory) << " Station " << imap->first.rawId() << ": " << imap->second << endl;
0515 
0516   LogTrace(theCategory) << "CheckMuonHits: " << all.size();
0517 
0518   // check order of muon measurements
0519   if ((all.size() > 1) && (all.front()->globalPosition().mag() > all.back()->globalPosition().mag())) {
0520     LogTrace(theCategory) << "reverse order: ";
0521     stable_sort(all.begin(), all.end(), RecHitLessByDet(alongMomentum));
0522   }
0523 }
0524 
0525 //
0526 // Get the hits from the first muon station (containing hits)
0527 //
0528 void GlobalMuonRefitter::getFirstHits(const reco::Track& muon,
0529                                       ConstRecHitContainer& all,
0530                                       ConstRecHitContainer& first) const {
0531   LogTrace(theCategory) << " GlobalMuonRefitter::getFirstHits\nall rechits length:" << all.size() << endl;
0532   first.clear();
0533 
0534   int station_to_keep = 999;
0535   vector<int> stations;
0536   for (ConstRecHitContainer::const_iterator ihit = all.begin(); ihit != all.end(); ++ihit) {
0537     int station = 0;
0538     bool use_it = true;
0539     DetId id = (*ihit)->geographicalId();
0540     unsigned raw_id = id.rawId();
0541     if (!(*ihit)->isValid())
0542       station = -1;
0543     else {
0544       if (id.det() == DetId::Muon) {
0545         switch (id.subdetId()) {
0546           case MuonSubdetId::DT:
0547             station = DTChamberId(raw_id).station();
0548             break;
0549           case MuonSubdetId::CSC:
0550             station = CSCDetId(raw_id).station();
0551             break;
0552           case MuonSubdetId::GEM:
0553             station = GEMDetId(raw_id).station();
0554             break;
0555           case MuonSubdetId::ME0:
0556             station = ME0DetId(raw_id).station();
0557             break;
0558           case MuonSubdetId::RPC:
0559             station = RPCDetId(raw_id).station();
0560             use_it = false;
0561             break;
0562         }
0563       }
0564     }
0565 
0566     if (use_it && station > 0 && station < station_to_keep)
0567       station_to_keep = station;
0568     stations.push_back(station);
0569     LogTrace(theCategory) << "rawId: " << raw_id << " station = " << station << " station_to_keep is now "
0570                           << station_to_keep;
0571   }
0572 
0573   if (station_to_keep <= 0 || station_to_keep > 4 || stations.size() != all.size())
0574     LogInfo(theCategory) << "failed to getFirstHits (all muon hits are outliers/bad ?)! station_to_keep = "
0575                          << station_to_keep << " stations.size " << stations.size() << " all.size " << all.size();
0576 
0577   for (unsigned i = 0; i < stations.size(); ++i)
0578     if (stations[i] >= 0 && stations[i] <= station_to_keep)
0579       first.push_back(all[i]);
0580 
0581   return;
0582 }
0583 
0584 //
0585 // select muon hits compatible with trajectory;
0586 // check hits in chambers with showers
0587 //
0588 GlobalMuonRefitter::ConstRecHitContainer GlobalMuonRefitter::selectMuonHits(const Trajectory& traj,
0589                                                                             const map<DetId, int>& hitMap) const {
0590   ConstRecHitContainer muonRecHits;
0591   const double globalChi2Cut = 200.0;
0592 
0593   vector<TrajectoryMeasurement> muonMeasurements = traj.measurements();
0594 
0595   // loop through all muon hits and skip hits with bad chi2 in chambers with high occupancy
0596   for (std::vector<TrajectoryMeasurement>::const_iterator im = muonMeasurements.begin(); im != muonMeasurements.end();
0597        im++) {
0598     if (!(*im).recHit()->isValid())
0599       continue;
0600     if ((*im).recHit()->det()->geographicalId().det() != DetId::Muon) {
0601       //      if ( ( chi2ndf < globalChi2Cut ) )
0602       muonRecHits.push_back((*im).recHit());
0603       continue;
0604     }
0605     const MuonTransientTrackingRecHit* immrh = dynamic_cast<const MuonTransientTrackingRecHit*>((*im).recHit().get());
0606 
0607     DetId id = immrh->geographicalId();
0608     DetId chamberId;
0609     int threshold = 0;
0610     double chi2Cut = 0.0;
0611 
0612     // get station of hit if it is in DT
0613     if ((*immrh).isDT()) {
0614       DTChamberId did(id.rawId());
0615       chamberId = did;
0616       threshold = theHitThreshold;
0617       chi2Cut = theDTChi2Cut;
0618     }
0619     // get station of hit if it is in CSC
0620     else if ((*immrh).isCSC()) {
0621       CSCDetId did(id.rawId());
0622       chamberId = did.chamberId();
0623       threshold = theHitThreshold;
0624       chi2Cut = theCSCChi2Cut;
0625     }
0626     // get station of hit if it is in GEM
0627     else if ((*immrh).isGEM()) {
0628       GEMDetId did(id.rawId());
0629       chamberId = did.chamberId();
0630       threshold = theHitThreshold;
0631       chi2Cut = theGEMChi2Cut;
0632     }
0633     // get station of hit if it is in ME0
0634     else if ((*immrh).isME0()) {
0635       ME0DetId did(id.rawId());
0636       chamberId = did.chamberId();
0637       threshold = theHitThreshold;
0638       chi2Cut = theME0Chi2Cut;
0639     }
0640     // get station of hit if it is in RPC
0641     else if ((*immrh).isRPC()) {
0642       RPCDetId rpcid(id.rawId());
0643       chamberId = rpcid;
0644       threshold = theHitThreshold;
0645       chi2Cut = theRPCChi2Cut;
0646     } else
0647       continue;
0648 
0649     double chi2ndf = (*im).estimate() / (*im).recHit()->dimension();
0650 
0651     bool keep = true;
0652     map<DetId, int>::const_iterator imap = hitMap.find(chamberId);
0653     if (imap != hitMap.end())
0654       if (imap->second > threshold)
0655         keep = false;
0656 
0657     if ((keep || (chi2ndf < chi2Cut)) && (chi2ndf < globalChi2Cut)) {
0658       muonRecHits.push_back((*im).recHit());
0659     } else {
0660       LogTrace(theCategory) << "Skip hit: " << id.rawId() << " chi2=" << chi2ndf << " ( threshold: " << chi2Cut
0661                             << ") Det: " << imap->second << endl;
0662     }
0663   }
0664 
0665   // check order of rechits
0666   reverse(muonRecHits.begin(), muonRecHits.end());
0667   return muonRecHits;
0668 }
0669 
0670 //
0671 // print RecHits
0672 //
0673 void GlobalMuonRefitter::printHits(const ConstRecHitContainer& hits) const {
0674   LogTrace(theCategory) << "Used RecHits: " << hits.size();
0675   for (ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
0676     if (!(*ir)->isValid()) {
0677       LogTrace(theCategory) << "invalid RecHit";
0678       continue;
0679     }
0680 
0681     const GlobalPoint& pos = (*ir)->globalPosition();
0682 
0683     LogTrace(theCategory) << "r = " << sqrt(pos.x() * pos.x() + pos.y() * pos.y()) << "  z = " << pos.z()
0684                           << "  dimension = " << (*ir)->dimension()
0685                           << "  det = " << (*ir)->det()->geographicalId().det()
0686                           << "  subdet = " << (*ir)->det()->subDetector()
0687                           << "  raw id = " << (*ir)->det()->geographicalId().rawId();
0688   }
0689 }
0690 
0691 //
0692 // add Trajectory* to TrackCand if not already present
0693 //
0694 GlobalMuonRefitter::RefitDirection GlobalMuonRefitter::checkRecHitsOrdering(
0695     const TransientTrackingRecHit::ConstRecHitContainer& recHits) const {
0696   if (!recHits.empty()) {
0697     ConstRecHitContainer::const_iterator frontHit = recHits.begin();
0698     ConstRecHitContainer::const_iterator backHit = recHits.end() - 1;
0699     while (!(*frontHit)->isValid() && frontHit != backHit) {
0700       frontHit++;
0701     }
0702     while (!(*backHit)->isValid() && backHit != frontHit) {
0703       backHit--;
0704     }
0705 
0706     double rFirst = (*frontHit)->globalPosition().mag();
0707     double rLast = (*backHit)->globalPosition().mag();
0708 
0709     if (rFirst < rLast)
0710       return insideOut;
0711     else if (rFirst > rLast)
0712       return outsideIn;
0713     else {
0714       LogError(theCategory) << "Impossible determine the rechits order" << endl;
0715       return undetermined;
0716     }
0717   } else {
0718     LogError(theCategory) << "Impossible determine the rechits order" << endl;
0719     return undetermined;
0720   }
0721 }
0722 
0723 //
0724 // Convert Tracks into Trajectories with a given set of hits
0725 //
0726 vector<Trajectory> GlobalMuonRefitter::transform(
0727     const reco::Track& newTrack,
0728     const reco::TransientTrack track,
0729     const TransientTrackingRecHit::ConstRecHitContainer& urecHitsForReFit) const {
0730   TransientTrackingRecHit::ConstRecHitContainer recHitsForReFit = urecHitsForReFit;
0731   LogTrace(theCategory) << "GlobalMuonRefitter::transform: " << recHitsForReFit.size() << " hits:";
0732   printHits(recHitsForReFit);
0733 
0734   if (recHitsForReFit.size() < 2)
0735     return vector<Trajectory>();
0736 
0737   // Check the order of the rechits
0738   RefitDirection recHitsOrder = checkRecHitsOrdering(recHitsForReFit);
0739 
0740   LogTrace(theCategory) << "checkRecHitsOrdering() returned " << recHitsOrder << ", theRefitDirection is "
0741                         << theRefitDirection << " (insideOut == " << insideOut << ", outsideIn == " << outsideIn << ")";
0742 
0743   // Reverse the order in the case of inconsistency between the fit direction and the rechit order
0744   if (theRefitDirection != recHitsOrder)
0745     reverse(recHitsForReFit.begin(), recHitsForReFit.end());
0746 
0747   // Even though we checked the rechits' ordering above, we may have
0748   // already flipped them elsewhere (getFirstHits() is such a
0749   // culprit). Use the global positions of the states and the desired
0750   // refit direction to find the starting TSOS.
0751   TrajectoryStateOnSurface firstTSOS, lastTSOS;
0752   unsigned int innerId;  //UNUSED: outerId;
0753   bool order_swapped = track.outermostMeasurementState().globalPosition().mag() <
0754                        track.innermostMeasurementState().globalPosition().mag();
0755   bool inner_is_first;
0756   LogTrace(theCategory) << "order swapped? " << order_swapped;
0757 
0758   // Fill the starting state, depending on the ordering above.
0759   if ((theRefitDirection == insideOut && !order_swapped) || (theRefitDirection == outsideIn && order_swapped)) {
0760     innerId = newTrack.innerDetId();
0761     //UNUSED:    outerId   = newTrack.outerDetId();
0762     firstTSOS = track.innermostMeasurementState();
0763     lastTSOS = track.outermostMeasurementState();
0764     inner_is_first = true;
0765   } else {
0766     innerId = newTrack.outerDetId();
0767     //UNUSED:    outerId   = newTrack.innerDetId();
0768     firstTSOS = track.outermostMeasurementState();
0769     lastTSOS = track.innermostMeasurementState();
0770     inner_is_first = false;
0771   }
0772 
0773   LogTrace(theCategory) << "firstTSOS: inner_is_first? " << inner_is_first << " globalPosition is "
0774                         << firstTSOS.globalPosition() << " innerId is " << innerId;
0775 
0776   if (!firstTSOS.isValid()) {
0777     LogWarning(theCategory) << "Error wrong initial state!" << endl;
0778     return vector<Trajectory>();
0779   }
0780 
0781   firstTSOS.rescaleError(theRescaleErrorFactor);
0782 
0783   // This is the only way to get a TrajectorySeed with settable propagation direction
0784   PTrajectoryStateOnDet garbage1;
0785   edm::OwnVector<TrackingRecHit> garbage2;
0786 
0787   // These lines cause the code to ignore completely what was set
0788   // above, and force propDir for tracks from collisions!
0789   //  if(propDir == alongMomentum && theRefitDirection == outsideIn)  propDir=oppositeToMomentum;
0790   //  if(propDir == oppositeToMomentum && theRefitDirection == insideOut) propDir=alongMomentum;
0791 
0792   const TrajectoryStateOnSurface& tsosForDir = inner_is_first ? lastTSOS : firstTSOS;
0793   PropagationDirection propDir =
0794       (tsosForDir.globalPosition().basicVector().dot(tsosForDir.globalMomentum().basicVector()) > 0)
0795           ? alongMomentum
0796           : oppositeToMomentum;
0797   LogTrace(theCategory) << "propDir based on firstTSOS x dot p is " << propDir << " (alongMomentum == " << alongMomentum
0798                         << ", oppositeToMomentum == " << oppositeToMomentum << ")";
0799 
0800   // Additional propagation diretcion determination logic for cosmic muons
0801   if (theCosmicFlag) {
0802     PropagationDirection propDir_first =
0803         (firstTSOS.globalPosition().basicVector().dot(firstTSOS.globalMomentum().basicVector()) > 0)
0804             ? alongMomentum
0805             : oppositeToMomentum;
0806     PropagationDirection propDir_last =
0807         (lastTSOS.globalPosition().basicVector().dot(lastTSOS.globalMomentum().basicVector()) > 0) ? alongMomentum
0808                                                                                                    : oppositeToMomentum;
0809     LogTrace(theCategory) << "propDir_first " << propDir_first << ", propdir_last " << propDir_last << " : they "
0810                           << (propDir_first == propDir_last ? "agree" : "disagree");
0811 
0812     int y_count = 0;
0813     for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator it = recHitsForReFit.begin();
0814          it != recHitsForReFit.end();
0815          ++it) {
0816       if ((*it)->globalPosition().y() > 0)
0817         ++y_count;
0818       else
0819         --y_count;
0820     }
0821 
0822     PropagationDirection propDir_ycount = alongMomentum;
0823     if (y_count > 0) {
0824       if (theRefitDirection == insideOut)
0825         propDir_ycount = oppositeToMomentum;
0826       else if (theRefitDirection == outsideIn)
0827         propDir_ycount = alongMomentum;
0828     } else {
0829       if (theRefitDirection == insideOut)
0830         propDir_ycount = alongMomentum;
0831       else if (theRefitDirection == outsideIn)
0832         propDir_ycount = oppositeToMomentum;
0833     }
0834 
0835     LogTrace(theCategory) << "y_count = " << y_count << "; based on geometrically-outermost TSOS, propDir is "
0836                           << propDir << ": " << (propDir == propDir_ycount ? "agrees" : "disagrees")
0837                           << " with ycount determination";
0838 
0839     if (propDir_first != propDir_last) {
0840       LogTrace(theCategory) << "since first/last disagreed, using y_count propDir";
0841       propDir = propDir_ycount;
0842     }
0843   }
0844 
0845   TrajectorySeed seed(garbage1, garbage2, propDir);
0846 
0847   if (recHitsForReFit.front()->geographicalId() != DetId(innerId)) {
0848     LogDebug(theCategory) << "Propagation occured" << endl;
0849     LogTrace(theCategory) << "propagating firstTSOS at " << firstTSOS.globalPosition()
0850                           << " to first rechit with surface pos "
0851                           << recHitsForReFit.front()->det()->surface().toGlobal(LocalPoint(0, 0, 0));
0852     firstTSOS =
0853         theService->propagator(thePropagatorName)->propagate(firstTSOS, recHitsForReFit.front()->det()->surface());
0854     if (!firstTSOS.isValid()) {
0855       LogDebug(theCategory) << "Propagation error!" << endl;
0856       return vector<Trajectory>();
0857     }
0858   }
0859 
0860   LogDebug(theCategory) << " GlobalMuonRefitter : theFitter " << propDir << endl;
0861   LogDebug(theCategory) << "                      First TSOS: " << firstTSOS.globalPosition()
0862                         << "  p=" << firstTSOS.globalMomentum() << " = " << firstTSOS.globalMomentum().mag() << endl;
0863 
0864   LogDebug(theCategory) << "                      Starting seed: "
0865                         << " nHits= " << seed.nHits() << " tsos: " << seed.startingState().parameters().position()
0866                         << "  p=" << seed.startingState().parameters().momentum() << endl;
0867 
0868   LogDebug(theCategory) << "                      RecHits: " << recHitsForReFit.size() << endl;
0869 
0870   vector<Trajectory> trajectories = theFitter->fit(seed, recHitsForReFit, firstTSOS);
0871 
0872   if (trajectories.empty()) {
0873     LogDebug(theCategory) << "No Track refitted!" << endl;
0874     return vector<Trajectory>();
0875   }
0876   return trajectories;
0877 }
0878 
0879 //
0880 // Remove Selected Station Rec Hits
0881 //
0882 GlobalMuonRefitter::ConstRecHitContainer GlobalMuonRefitter::getRidOfSelectStationHits(
0883     const ConstRecHitContainer& hits, const TrackerTopology* tTopo) const {
0884   ConstRecHitContainer results;
0885   ConstRecHitContainer::const_iterator it = hits.begin();
0886   for (; it != hits.end(); it++) {
0887     DetId id = (*it)->geographicalId();
0888 
0889     //Check that this is a Muon hit that we're toying with -- else pass on this because the hacker is a moron / not careful
0890 
0891     if (id.det() == DetId::Tracker && theTrackerSkipSystem > 0) {
0892       int layer = -999;
0893       int disk = -999;
0894       int wheel = -999;
0895       if (id.subdetId() == theTrackerSkipSystem) {
0896         //                              continue;  //caveat that just removes the whole system from refitting
0897 
0898         if (theTrackerSkipSystem == PXB) {
0899           layer = tTopo->pxbLayer(id);
0900         }
0901         if (theTrackerSkipSystem == TIB) {
0902           layer = tTopo->tibLayer(id);
0903         }
0904 
0905         if (theTrackerSkipSystem == TOB) {
0906           layer = tTopo->tobLayer(id);
0907         }
0908         if (theTrackerSkipSystem == PXF) {
0909           disk = tTopo->pxfDisk(id);
0910         }
0911         if (theTrackerSkipSystem == TID) {
0912           wheel = tTopo->tidWheel(id);
0913         }
0914         if (theTrackerSkipSystem == TEC) {
0915           wheel = tTopo->tecWheel(id);
0916         }
0917         if (theTrackerSkipSection >= 0 && layer == theTrackerSkipSection)
0918           continue;
0919         if (theTrackerSkipSection >= 0 && disk == theTrackerSkipSection)
0920           continue;
0921         if (theTrackerSkipSection >= 0 && wheel == theTrackerSkipSection)
0922           continue;
0923       }
0924     }
0925 
0926     if (id.det() == DetId::Muon && theSkipStation) {
0927       int station = -999;
0928       //UNUSED:      int wheel = -999;
0929       if (id.subdetId() == MuonSubdetId::DT) {
0930         DTChamberId did(id.rawId());
0931         station = did.station();
0932         //UNUSED:   wheel = did.wheel();
0933       } else if (id.subdetId() == MuonSubdetId::CSC) {
0934         CSCDetId did(id.rawId());
0935         station = did.station();
0936       } else if (id.subdetId() == MuonSubdetId::GEM) {
0937         GEMDetId did(id.rawId());
0938         station = did.station();
0939       } else if (id.subdetId() == MuonSubdetId::ME0) {
0940         ME0DetId did(id.rawId());
0941         station = did.station();
0942       } else if (id.subdetId() == MuonSubdetId::RPC) {
0943         RPCDetId rpcid(id.rawId());
0944         station = rpcid.station();
0945       }
0946       if (station == theSkipStation)
0947         continue;
0948     }
0949     results.push_back(*it);
0950   }
0951   return results;
0952 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!