MuonSeedGenerator/src/MuonOverlapSeedFromRecHits.cc

0001 #include "RecoMuon/MuonSeedGenerator/src/MuonOverlapSeedFromRecHits.h"
0002 #include "RecoMuon/MuonSeedGenerator/src/MuonDTSeedFromRecHits.h"
0003 #include "RecoMuon/MuonSeedGenerator/src/MuonCSCSeedFromRecHits.h"
0004 #include "RecoMuon/MuonSeedGenerator/src/MuonSeedPtExtractor.h"
0005 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0006 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0007 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0008 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "Geometry/CSCGeometry/interface/CSCChamberSpecs.h"
0011 #include <iomanip>
0012
0013 MuonOverlapSeedFromRecHits::MuonOverlapSeedFromRecHits() : MuonSeedFromRecHits() {}
0014
0015 std::vector<TrajectorySeed> MuonOverlapSeedFromRecHits::seeds() const {
0016   std::vector<TrajectorySeed> result;
0017   //@@ doesn't handle overlap between ME11 and ME12 correctly
0018   // sort by station
0019   MuonRecHitContainer barrelHits, endcapHits;
0020   for (MuonRecHitContainer::const_iterator iter = theRhits.begin(), end = theRhits.end(); iter != end; ++iter) {
0021     if ((*iter)->isDT()) {
0022       DTChamberId dtId((**iter).geographicalId().rawId());
0023       // try not doing seeds that start in DT station 2, if there'as a good single segment seed
0024       if (dtId.station() == 1 || (dtId.station() == 2 && (*iter)->dimension() == 2)) {
0025         barrelHits.push_back(*iter);
0026       }
0027     } else {
0028       endcapHits.push_back(*iter);
0029     }
0030   }
0031
0032   ConstMuonRecHitPointer bestSegment = bestHit(barrelHits, endcapHits);
0033   for (MuonRecHitContainer::const_iterator barrelHitItr = barrelHits.begin(), lastBarrelHit = barrelHits.end();
0034        barrelHitItr != lastBarrelHit;
0035        ++barrelHitItr) {
0036     for (MuonRecHitContainer::const_iterator endcapHitItr = endcapHits.begin(), lastEndcapHit = endcapHits.end();
0037          endcapHitItr != lastEndcapHit;
0038          ++endcapHitItr) {
0039       TrajectorySeed seed;
0040       bool good = makeSeed(*barrelHitItr, *endcapHitItr, bestSegment, seed);
0041       if (good)
0042         result.push_back(seed);
0043       // try just one seed
0044       return result;
0045     }
0046   }
0047
0048   //std::cout << "Overlap hits " << barrelHits.size() << " "
0049   //                             << endcapHits.size() << std::endl;
0050
0051   return result;
0052 }
0053
0054 bool MuonOverlapSeedFromRecHits::makeSeed(MuonTransientTrackingRecHit::ConstMuonRecHitPointer barrelHit,
0055                                           MuonTransientTrackingRecHit::ConstMuonRecHitPointer endcapHit,
0056                                           MuonTransientTrackingRecHit::ConstMuonRecHitPointer bestSegment,
0057                                           TrajectorySeed& result) const {
0058   std::vector<double> pts = thePtExtractor->pT_extract(barrelHit, endcapHit);
0059   double minpt = 3.;
0060   double pt = pts[0];
0061   double sigmapt = pts[1];
0062   // if too small, probably an error.  Keep trying.
0063   if (pt != 0) {
0064     if (fabs(pt) > minpt) {
0065       double maxpt = 2000.;
0066       if (pt > maxpt) {
0067         pt = maxpt;
0068         sigmapt = maxpt;
0069       }
0070       if (pt < -maxpt) {
0071         pt = -maxpt;
0072         sigmapt = maxpt;
0073       }
0074     }
0075
0076     result = createSeed(pt, sigmapt, bestSegment);
0077     //std::cout << "OVERLAPFITTED PT " << pt << " dphi " << dphi << " eta " << eta << std::endl;
0078     return true;
0079   }
0080   return false;
0081 }
0082
0083 MuonTransientTrackingRecHit::ConstMuonRecHitPointer MuonOverlapSeedFromRecHits::bestHit(
0084     const MuonTransientTrackingRecHit::MuonRecHitContainer& barrelHits,
0085     const MuonTransientTrackingRecHit::MuonRecHitContainer& endcapHits) const {
0086   MuonDTSeedFromRecHits dtSeeder;
0087   MuonCSCSeedFromRecHits cscSeeder;
0088
0089   ConstMuonRecHitPointer result;
0090   if (barrelHits.size() > endcapHits.size()) {
0091     result = dtSeeder.bestBarrelHit(barrelHits);
0092     if (result->dimension() == 2)
0093       result = cscSeeder.bestEndcapHit(endcapHits);
0094   } else {
0095     result = cscSeeder.bestEndcapHit(endcapHits);
0096   }
0097   return result;
0098 }