MuonSeedGenerator/src/MuonCSCSeedFromRecHits.cc

0001 #include "RecoMuon/MuonSeedGenerator/src/MuonCSCSeedFromRecHits.h"
0002 #include "RecoMuon/MuonSeedGenerator/src/MuonSeedPtExtractor.h"
0003 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0004 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0005 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0006 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0007 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0008 #include "DataFormats/Math/interface/deltaPhi.h"
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "Geometry/CSCGeometry/interface/CSCChamberSpecs.h"
0012 #include "Geometry/CSCGeometry/interface/CSCChamber.h"
0013 #include <iomanip>
0014
0015 MuonCSCSeedFromRecHits::MuonCSCSeedFromRecHits() : MuonSeedFromRecHits() {}
0016
0017 TrajectorySeed MuonCSCSeedFromRecHits::seed() const {
0018   TrajectorySeed result;
0019   if (theRhits.size() == 1) {
0020     //return createSeed(100., 100., theRhits[0]);
0021     makeDefaultSeed(result);
0022     return result;
0023   }
0024   //@@ doesn't handle overlap between ME11 and ME12 correctly
0025   // sort by station
0026   MuonRecHitContainer station1Hits, station2Hits, station3Hits, station4Hits;
0027   for (MuonRecHitContainer::const_iterator iter = theRhits.begin(), end = theRhits.end(); iter != end; ++iter) {
0028     int station = CSCDetId((*iter)->geographicalId().rawId()).station();
0029     if ((*iter)->isME0()) {
0030       station = 1;  //ME0DetId((*iter)->geographicalId().rawId()).station();
0031     }
0032
0033     if (station == 1) {
0034       station1Hits.push_back(*iter);
0035     } else if (station == 2) {
0036       station2Hits.push_back(*iter);
0037     } else if (station == 3) {
0038       station3Hits.push_back(*iter);
0039     } else if (station == 4) {
0040       station4Hits.push_back(*iter);
0041     }
0042   }
0043
0044   //std::cout << "Station hits " << station1Hits.size() << " "
0045   //                            << station2Hits.size() << " "
0046   //                             << station3Hits.size() << std::endl;
0047
0048   // see whether station 2 or station 3 is better
0049   MuonRecHitContainer* betterSecondHits = &station2Hits;
0050   MuonRecHitContainer* notAsGoodSecondHits = &station3Hits;
0051   if (!station2Hits.empty() && !station3Hits.empty()) {
0052     // swap if station 3 has better quailty
0053     if (segmentQuality(station3Hits[0]) < segmentQuality(station2Hits[0])) {
0054       betterSecondHits = &station3Hits;
0055       notAsGoodSecondHits = &station2Hits;
0056     }
0057   }
0058
0059   // now try to make pairs
0060   if (makeSeed(station1Hits, *betterSecondHits, result)) {
0061     return result;
0062   }
0063   if (makeSeed(station1Hits, *notAsGoodSecondHits, result)) {
0064     return result;
0065   }
0066   if (makeSeed(station2Hits, station3Hits, result)) {
0067     return result;
0068   }
0069   if (makeSeed(station1Hits, station4Hits, result)) {
0070     return result;
0071   }
0072
0073   // no luck
0074   makeDefaultSeed(result);
0075   return result;
0076 }
0077
0078 bool MuonCSCSeedFromRecHits::makeSeed(const MuonRecHitContainer& hits1,
0079                                       const MuonRecHitContainer& hits2,
0080                                       TrajectorySeed& seed) const {
0081   for (MuonRecHitContainer::const_iterator itr1 = hits1.begin(), end1 = hits1.end(); itr1 != end1; ++itr1) {
0082     CSCDetId cscId1((*itr1)->geographicalId().rawId());
0083     //int type1 = CSCChamberSpecs::whatChamberType(cscId1.station(), cscId1.ring());
0084
0085     for (MuonRecHitContainer::const_iterator itr2 = hits2.begin(), end2 = hits2.end(); itr2 != end2; ++itr2) {
0086       CSCDetId cscId2((*itr2)->geographicalId().rawId());
0087       //int type2 = CSCChamberSpecs::whatChamberType(cscId2.station(), cscId2.ring());
0088
0089       // take the first pair that comes along.  Probably want to rank them later
0090       std::vector<double> pts = thePtExtractor->pT_extract(*itr1, *itr2);
0091
0092       double pt = pts[0];
0093       double sigmapt = pts[1];
0094       double minpt = 3.;
0095
0096       // if too small, probably an error.  Keep trying.
0097       if (fabs(pt) > minpt) {
0098         double maxpt = 2000.;
0099         if (pt > maxpt) {
0100           pt = maxpt;
0101           sigmapt = maxpt;
0102         }
0103         if (pt < -maxpt) {
0104           pt = -maxpt;
0105           sigmapt = maxpt;
0106         }
0107
0108         // get the position and direction from the higher-quality segment
0109         ConstMuonRecHitPointer bestSeg = bestEndcapHit(theRhits);
0110         seed = createSeed(pt, sigmapt, bestSeg);
0111
0112         //std::cout << "FITTED TIMESPT " << pt << " dphi " << dphi << " eta " << eta << std::endl;
0113         return true;
0114       }
0115     }
0116   }
0117
0118   // guess it didn't find one
0119   //std::cout << "NOTHING FOUND ! " << std::endl;
0120   return false;
0121 }
0122
0123 //typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer;
0124 int MuonCSCSeedFromRecHits::segmentQuality(ConstMuonRecHitPointer segment) const {
0125   int Nchi2 = 0;
0126   int quality = 0;
0127   int nhits = segment->recHits().size();
0128   if (segment->chi2() / (nhits * 2. - 4.) > 3.)
0129     Nchi2 = 1;
0130   if (segment->chi2() / (nhits * 2. - 4.) > 9.)
0131     Nchi2 = 2;
0132
0133   if (nhits > 4)
0134     quality = 1 + Nchi2;
0135   if (nhits == 4)
0136     quality = 3 + Nchi2;
0137   if (nhits == 3)
0138     quality = 5 + Nchi2;
0139
0140   float dPhiGloDir = fabs(deltaPhi(segment->globalPosition().barePhi(), segment->globalDirection().barePhi()));
0141
0142   if (dPhiGloDir > .2)
0143     ++quality;
0144   // add a penalty for being ME1A if the chamber is ganged
0145   if (segment->isCSC() and CSCDetId(segment->geographicalId()).ring() == 4) {
0146     const auto chamber = dynamic_cast<const CSCChamber*>(segment->det());
0147     if (chamber->specs()->gangedStrips())
0148       ++quality;
0149   }
0150
0151   return quality;
0152 }
0153
0154 MuonCSCSeedFromRecHits::ConstMuonRecHitPointer MuonCSCSeedFromRecHits::bestEndcapHit(
0155     const MuonRecHitContainer& endcapHits) const {
0156   MuonRecHitPointer me1 = nullptr, meit = nullptr;
0157   float dPhiGloDir = .0;          //  +v
0158   float bestdPhiGloDir = M_PI;    //  +v
0159   int quality1 = 0, quality = 0;  //  +v  I= 5,6-p. / II= 4p.  / III= 3p.
0160
0161   for (MuonRecHitContainer::const_iterator iter = endcapHits.begin(); iter != endcapHits.end(); iter++) {
0162     if (!(*iter)->isCSC() && !(*iter)->isME0())
0163       continue;
0164
0165     // tmp compar. Glob-Dir for the same tr-segm:
0166
0167     meit = *iter;
0168
0169     quality = segmentQuality(meit);
0170
0171     dPhiGloDir = fabs(deltaPhi(meit->globalPosition().barePhi(), meit->globalDirection().barePhi()));
0172
0173     if (!me1) {
0174       me1 = meit;
0175       quality1 = quality;
0176       bestdPhiGloDir = dPhiGloDir;
0177     }
0178
0179     if (me1) {
0180       if (!me1->isValid()) {
0181         me1 = meit;
0182         quality1 = quality;
0183         bestdPhiGloDir = dPhiGloDir;
0184       }
0185
0186       if (me1->isValid() && quality < quality1) {
0187         me1 = meit;
0188         quality1 = quality;
0189         bestdPhiGloDir = dPhiGloDir;
0190       }
0191
0192       if (me1->isValid() && bestdPhiGloDir > .03) {
0193         if (dPhiGloDir < bestdPhiGloDir - .01 && quality == quality1) {
0194           me1 = meit;
0195           quality1 = quality;
0196           bestdPhiGloDir = dPhiGloDir;
0197         }
0198       }
0199     }
0200
0201   }  //  iter
0202
0203   return me1;
0204 }
0205
0206 void MuonCSCSeedFromRecHits::makeDefaultSeed(TrajectorySeed& seed) const {
0207   //Search ME1  ...
0208   ConstMuonRecHitPointer me1 = bestEndcapHit(theRhits);
0209   //  bool good=false;
0210
0211   if (me1 && me1->isValid()) {
0212     //revert if a LogTrace or smth is necessary    good =
0213     createDefaultEndcapSeed(me1, seed);
0214   }
0215 }
0216
0217 bool MuonCSCSeedFromRecHits::createDefaultEndcapSeed(ConstMuonRecHitPointer last, TrajectorySeed& seed) const {
0218   //float momentum = computeDefaultPt(last);
0219   std::vector<double> momentum = thePtExtractor->pT_extract(last, last);
0220   seed = createSeed(momentum[0], momentum[1], last);
0221   return true;
0222 }
0223
0224 void MuonCSCSeedFromRecHits::analyze() const {
0225   for (MuonRecHitContainer::const_iterator iter = theRhits.begin(); iter != theRhits.end(); iter++) {
0226     if (!(*iter)->isCSC())
0227       continue;
0228     for (MuonRecHitContainer::const_iterator iter2 = iter + 1; iter2 != theRhits.end(); ++iter2) {
0229       if (!(*iter2)->isCSC())
0230         continue;
0231
0232       CSCDetId cscId1((*iter)->geographicalId().rawId());
0233       CSCDetId cscId2((*iter2)->geographicalId().rawId());
0234       double dphi = deltaPhi((**iter).globalPosition().barePhi(), (**iter2).globalPosition().barePhi());
0235
0236       int type1 = cscId1.iChamberType();
0237       int type2 = cscId2.iChamberType();
0238
0239       // say the lower station first
0240       if (type1 < type2) {
0241         std::cout << "HITPAIRA," << type1 << type2 << "," << dphi << "," << (**iter).globalPosition().eta()
0242                   << std::endl;
0243       }
0244       if (type2 < type1) {
0245         std::cout << "HITPAIRB," << type2 << type1 << "," << -dphi << "," << (**iter2).globalPosition().eta()
0246                   << std::endl;
0247       }
0248       if (type1 == type2) {
0249         std::cout << "HITPAIRSAMESTATION," << type1 << cscId1.ring() << cscId2.ring() << "," << dphi << ","
0250                   << (**iter).globalPosition().eta() << std::endl;
0251       }
0252     }
0253   }
0254 }