MuonSeedGenerator/src/MuonSeedCleaner.cc

0001 /**
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Shih-Chuan Kao, Dominique Fortin - UCR
0005  */
0006
0007 #include <RecoMuon/MuonSeedGenerator/interface/MuonSeedCleaner.h>
0008
0009 // Data Formats
0010 #include <DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h>
0011 #include <DataFormats/TrajectorySeed/interface/TrajectorySeed.h>
0012 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
0013 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0014 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
0015 #include <DataFormats/CSCRecHit/interface/CSCSegment.h>
0016 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
0017 #include <DataFormats/DTRecHit/interface/DTRecSegment4D.h>
0018
0019 // Geometry
0020 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0021 #include <TrackingTools/DetLayers/interface/DetLayer.h>
0022 #include <RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h>
0023 #include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>
0024 #include <RecoMuon/Records/interface/MuonRecoGeometryRecord.h>
0025
0026 // muon service
0027 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
0028 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
0029
0030 // Framework
0031 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include <FWCore/Framework/interface/ESHandle.h>
0034 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0035 #include <DataFormats/Common/interface/Handle.h>
0036
0037 // C++
0038 #include <vector>
0039 #include <deque>
0040 #include <utility>
0041
0042 //typedef std::pair<double, TrajectorySeed> seedpr ;
0043 //static bool ptDecreasing(const seedpr s1, const seedpr s2) { return ( s1.first > s2.first ); }
0044 static bool lengthSorting(const TrajectorySeed& s1, const TrajectorySeed& s2) { return (s1.nHits() > s2.nHits()); }
0045
0046 /*
0047  * Constructor
0048  */
0049 MuonSeedCleaner::MuonSeedCleaner(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC) {
0050   // Local Debug flag
0051   debug = pset.getParameter<bool>("DebugMuonSeed");
0052
0053   // muon service
0054   edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
0055   theService = new MuonServiceProxy(serviceParameters, std::move(iC));
0056 }
0057
0058 /*
0059  * Destructor
0060  */
0061 MuonSeedCleaner::~MuonSeedCleaner() {
0062   if (theService)
0063     delete theService;
0064 }
0065
0066 /***********************************
0067  *
0068  *  Seed Cleaner
0069  *
0070  ***********************************/
0071
0072 std::vector<TrajectorySeed> MuonSeedCleaner::seedCleaner(const edm::EventSetup& eventSetup,
0073                                                          std::vector<TrajectorySeed>& seeds) {
0074   theService->update(eventSetup);
0075
0076   std::vector<TrajectorySeed> FinalSeeds;
0077
0078   // group the seeds
0079   std::vector<SeedContainer> theCollection = GroupSeeds(seeds);
0080
0081   // ckeck each group and pick the good one
0082   for (size_t i = 0; i < theCollection.size(); i++) {
0083     // separate seeds w/ more than 1 segments and w/ 1st layer segment information
0084     SeedContainer goodSeeds = SeedCandidates(theCollection[i], true);
0085     SeedContainer otherSeeds = SeedCandidates(theCollection[i], false);
0086     if (MomentumFilter(goodSeeds)) {
0087       //std::cout<<" == type1 "<<std::endl;
0088       TrajectorySeed bestSeed = Chi2LengthSelection(goodSeeds);
0089       FinalSeeds.push_back(bestSeed);
0090
0091       GlobalPoint seedgp = SeedPosition(bestSeed);
0092       double eta = fabs(seedgp.eta());
0093       if (goodSeeds.size() > 2 && eta > 1.5) {
0094         TrajectorySeed anotherSeed = MoreRecHits(goodSeeds);
0095         FinalSeeds.push_back(anotherSeed);
0096       }
0097     } else if (MomentumFilter(otherSeeds)) {
0098       //std::cout<<" == type2 "<<std::endl;
0099       TrajectorySeed bestSeed = MoreRecHits(otherSeeds);
0100       FinalSeeds.push_back(bestSeed);
0101
0102       GlobalPoint seedgp = SeedPosition(bestSeed);
0103       double eta = fabs(seedgp.eta());
0104       if (otherSeeds.size() > 2 && eta > 1.5) {
0105         TrajectorySeed anotherSeed = LeanHighMomentum(otherSeeds);
0106         FinalSeeds.push_back(anotherSeed);
0107       }
0108     } else {
0109       //std::cout<<" == type3 "<<std::endl;
0110       TrajectorySeed bestSeed = LeanHighMomentum(theCollection[i]);
0111       FinalSeeds.push_back(bestSeed);
0112
0113       GlobalPoint seedgp = SeedPosition(bestSeed);
0114       double eta = fabs(seedgp.eta());
0115       if (theCollection.size() > 2 && eta > 1.5) {
0116         TrajectorySeed anotherSeed = BiggerCone(theCollection[i]);
0117         FinalSeeds.push_back(anotherSeed);
0118       }
0119     }
0120   }
0121   return FinalSeeds;
0122 }
0123
0124 TrajectorySeed MuonSeedCleaner::Chi2LengthSelection(std::vector<TrajectorySeed>& seeds) {
0125   if (seeds.size() == 1)
0126     return seeds[0];
0127
0128   int winner = 0;
0129   int moreHits = 0;
0130   double bestChi2 = 99999.;
0131   for (size_t i = 0; i < seeds.size(); i++) {
0132     // 1. fill out the Nchi2 of segments of the seed
0133     //GlobalVector mom = SeedMomentum( seeds[i] ); // temporary use for debugging
0134     //double pt = sqrt( (mom.x()*mom.x()) + (mom.y()*mom.y()) );
0135     //std::cout<<" > SEED"<<i<<"  pt:"<<pt<< std::endl;
0136
0137     double theChi2 = SeedChi2(seeds[i]);
0138     double dChi2 = fabs(1. - (theChi2 / bestChi2));
0139     int theHits = seeds[i].nHits();
0140     int dHits = theHits - moreHits;
0141     //std::cout<<" -----  "<<std::endl;
0142
0143     // 2. better chi2
0144     if (theChi2 < bestChi2 && dChi2 > 0.05) {
0145       winner = static_cast<int>(i);
0146       bestChi2 = theChi2;
0147       moreHits = theHits;
0148     }
0149     // 3. if chi2 is not much better, pick more rechits one
0150     if (theChi2 >= bestChi2 && dChi2 < 0.05 && dHits > 0) {
0151       winner = static_cast<int>(i);
0152       bestChi2 = theChi2;
0153       moreHits = theHits;
0154     }
0155   }
0156   //std::cout<<" Winner is "<< winner <<std::endl;
0157   TrajectorySeed theSeed = seeds[winner];
0158   seeds.erase(seeds.begin() + winner);
0159   return theSeed;
0160 }
0161
0162 TrajectorySeed MuonSeedCleaner::BiggerCone(std::vector<TrajectorySeed>& seeds) {
0163   if (seeds.size() == 1)
0164     return seeds[0];
0165
0166   float biggerProjErr = 9999.;
0167   int winner = 0;
0168   AlgebraicSymMatrix mat(5, 0);
0169   for (size_t i = 0; i < seeds.size(); i++) {
0170     auto r1 = seeds[i].recHits().begin();
0171     mat = r1->parametersError().similarityT(r1->projectionMatrix());
0172
0173     int NRecHits = NRecHitsFromSegment(*r1);
0174
0175     float ddx = mat[1][1];
0176     float ddy = mat[2][2];
0177     float dxx = mat[3][3];
0178     float dyy = mat[4][4];
0179     float projectErr = sqrt((ddx * 10000.) + (ddy * 10000.) + dxx + dyy);
0180
0181     if (NRecHits < 5)
0182       continue;
0183     if (projectErr < biggerProjErr)
0184       continue;
0185
0186     winner = static_cast<int>(i);
0187     biggerProjErr = projectErr;
0188   }
0189   TrajectorySeed theSeed = seeds[winner];
0190   seeds.erase(seeds.begin() + winner);
0191   return theSeed;
0192 }
0193
0194 TrajectorySeed MuonSeedCleaner::LeanHighMomentum(std::vector<TrajectorySeed>& seeds) {
0195   if (seeds.size() == 1)
0196     return seeds[0];
0197
0198   double highestPt = 0.;
0199   int winner = 0;
0200   for (size_t i = 0; i < seeds.size(); i++) {
0201     GlobalVector mom = SeedMomentum(seeds[i]);
0202     double pt = sqrt((mom.x() * mom.x()) + (mom.y() * mom.y()));
0203     if (pt > highestPt) {
0204       winner = static_cast<int>(i);
0205       highestPt = pt;
0206     }
0207   }
0208   TrajectorySeed theSeed = seeds[winner];
0209   seeds.erase(seeds.begin() + winner);
0210   return theSeed;
0211 }
0212
0213 TrajectorySeed MuonSeedCleaner::MoreRecHits(std::vector<TrajectorySeed>& seeds) {
0214   if (seeds.size() == 1)
0215     return seeds[0];
0216
0217   int winner = 0;
0218   int moreHits = 0;
0219   double betterChi2 = 99999.;
0220   for (size_t i = 0; i < seeds.size(); i++) {
0221     int theHits = 0;
0222     for (auto const& r1 : seeds[i].recHits()) {
0223       theHits += NRecHitsFromSegment(r1);
0224     }
0225
0226     double theChi2 = SeedChi2(seeds[i]);
0227
0228     if (theHits == moreHits && theChi2 < betterChi2) {
0229       betterChi2 = theChi2;
0230       winner = static_cast<int>(i);
0231     }
0232     if (theHits > moreHits) {
0233       moreHits = theHits;
0234       betterChi2 = theChi2;
0235       winner = static_cast<int>(i);
0236     }
0237   }
0238   TrajectorySeed theSeed = seeds[winner];
0239   seeds.erase(seeds.begin() + winner);
0240   return theSeed;
0241 }
0242
0243 SeedContainer MuonSeedCleaner::LengthFilter(std::vector<TrajectorySeed>& seeds) {
0244   SeedContainer longSeeds;
0245   int NSegs = 0;
0246   for (size_t i = 0; i < seeds.size(); i++) {
0247     int theLength = static_cast<int>(seeds[i].nHits());
0248     if (theLength > NSegs) {
0249       NSegs = theLength;
0250       longSeeds.clear();
0251       longSeeds.push_back(seeds[i]);
0252     } else if (theLength == NSegs) {
0253       longSeeds.push_back(seeds[i]);
0254     } else {
0255       continue;
0256     }
0257   }
0258   //std::cout<<" final Length :"<<NSegs<<std::endl;
0259
0260   return longSeeds;
0261 }
0262
0263 bool MuonSeedCleaner::MomentumFilter(std::vector<TrajectorySeed>& seeds) {
0264   bool findgoodMomentum = false;
0265   SeedContainer goodMomentumSeeds = seeds;
0266   seeds.clear();
0267   for (size_t i = 0; i < goodMomentumSeeds.size(); i++) {
0268     GlobalVector mom = SeedMomentum(goodMomentumSeeds[i]);
0269     double pt = sqrt((mom.x() * mom.x()) + (mom.y() * mom.y()));
0270     if (pt < 6. || pt > 2000.)
0271       continue;
0272     //if ( pt < 6. ) continue;
0273     //std::cout<<" passed momentum :"<< pt <<std::endl;
0274     seeds.push_back(goodMomentumSeeds[i]);
0275     findgoodMomentum = true;
0276   }
0277   if (seeds.empty())
0278     seeds = goodMomentumSeeds;
0279
0280   return findgoodMomentum;
0281 }
0282
0283 SeedContainer MuonSeedCleaner::SeedCandidates(std::vector<TrajectorySeed>& seeds, bool good) {
0284   SeedContainer theCandidate;
0285   theCandidate.clear();
0286
0287   bool longSeed = false;
0288   bool withFirstLayer = false;
0289
0290   //std::cout<<"***** Seed Classification *****"<< seeds.size() <<std::endl;
0291   for (size_t i = 0; i < seeds.size(); i++) {
0292     if (seeds[i].nHits() > 1)
0293       longSeed = true;
0294     //std::cout<<"  Seed: "<<i<<" w/"<<seeds[i].nHits()<<" segs "<<std::endl;
0295     // looking for 1st layer segment
0296     for (auto const& r1 : seeds[i].recHits()) {
0297       const GeomDet* gdet = theService->trackingGeometry()->idToDet(r1.geographicalId());
0298       DetId geoId = gdet->geographicalId();
0299
0300       if (geoId.subdetId() == MuonSubdetId::DT) {
0301         DTChamberId DT_Id(r1.geographicalId());
0302         //std::cout<<" ID:"<<DT_Id <<" pos:"<< r1->localPosition()  <<std::endl;
0303         if (DT_Id.station() != 1)
0304           continue;
0305         withFirstLayer = true;
0306       }
0307       if (geoId.subdetId() == MuonSubdetId::CSC) {
0308         CSCDetId CSC_Id = CSCDetId(r1.geographicalId());
0309         //std::cout<<" ID:"<<CSC_Id <<" pos:"<< r1->localPosition()  <<std::endl;
0310         if (CSC_Id.station() != 1)
0311           continue;
0312         withFirstLayer = true;
0313       }
0314     }
0315     bool goodseed = (longSeed && withFirstLayer) ? true : false;
0316
0317     if (goodseed == good)
0318       theCandidate.push_back(seeds[i]);
0319   }
0320   return theCandidate;
0321 }
0322
0323 std::vector<SeedContainer> MuonSeedCleaner::GroupSeeds(std::vector<TrajectorySeed>& seeds) {
0324   std::vector<SeedContainer> seedCollection;
0325   seedCollection.clear();
0326   std::vector<TrajectorySeed> theGroup;
0327   std::vector<bool> usedSeed(seeds.size(), false);
0328
0329   // categorize seeds by comparing overlapping segments or a certian eta-phi cone
0330   for (size_t i = 0; i < seeds.size(); i++) {
0331     if (usedSeed[i])
0332       continue;
0333     theGroup.push_back(seeds[i]);
0334     usedSeed[i] = true;
0335
0336     GlobalPoint pos1 = SeedPosition(seeds[i]);
0337
0338     for (size_t j = i + 1; j < seeds.size(); j++) {
0339       // 1.1 seeds with overlaaping segments will be grouped together
0340       unsigned int overlapping = OverlapSegments(seeds[i], seeds[j]);
0341       if (!usedSeed[j] && overlapping > 0) {
0342         // reject the identical seeds
0343         if (seeds[i].nHits() == overlapping && seeds[j].nHits() == overlapping) {
0344           usedSeed[j] = true;
0345           continue;
0346         }
0347         theGroup.push_back(seeds[j]);
0348         usedSeed[j] = true;
0349       }
0350       if (usedSeed[j])
0351         continue;
0352
0353       // 1.2 seeds in a certain cone are grouped together
0354       GlobalPoint pos2 = SeedPosition(seeds[j]);
0355       double dh = pos1.eta() - pos2.eta();
0356       double df = pos1.phi() - pos2.phi();
0357       double dR = sqrt((dh * dh) + (df * df));
0358
0359       if (dR > 0.3 && seeds[j].nHits() == 1)
0360         continue;
0361       if (dR > 0.2 && seeds[j].nHits() > 1)
0362         continue;
0363       theGroup.push_back(seeds[j]);
0364       usedSeed[j] = true;
0365     }
0366     sort(theGroup.begin(), theGroup.end(), lengthSorting);
0367     seedCollection.push_back(theGroup);
0368     //std::cout<<" group "<<seedCollection.size() <<" w/"<< theGroup.size() <<" seeds"<<std::endl;
0369     theGroup.clear();
0370   }
0371   return seedCollection;
0372 }
0373
0374 unsigned int MuonSeedCleaner::OverlapSegments(const TrajectorySeed& seed1, const TrajectorySeed& seed2) {
0375   unsigned int overlapping = 0;
0376   for (auto const& r1 : seed1.recHits()) {
0377     DetId id1 = r1.geographicalId();
0378     const GeomDet* gdet1 = theService->trackingGeometry()->idToDet(id1);
0379     GlobalPoint gp1 = gdet1->toGlobal(r1.localPosition());
0380
0381     for (auto const& r2 : seed2.recHits()) {
0382       DetId id2 = r2.geographicalId();
0383       if (id1 != id2)
0384         continue;
0385
0386       const GeomDet* gdet2 = theService->trackingGeometry()->idToDet(id2);
0387       GlobalPoint gp2 = gdet2->toGlobal(r2.localPosition());
0388
0389       double dx = gp1.x() - gp2.x();
0390       double dy = gp1.y() - gp2.y();
0391       double dz = gp1.z() - gp2.z();
0392       double dL = sqrt(dx * dx + dy * dy + dz * dz);
0393
0394       if (dL < 1.)
0395         overlapping++;
0396     }
0397   }
0398   return overlapping;
0399 }
0400
0401 double MuonSeedCleaner::SeedChi2(const TrajectorySeed& seed) {
0402   double theChi2 = 0.;
0403   for (auto const& r1 : seed.recHits()) {
0404     //std::cout<<"    segmet : "<<it <<std::endl;
0405     theChi2 += NChi2OfSegment(r1);
0406   }
0407   theChi2 = theChi2 / seed.nHits();
0408
0409   //std::cout<<" final Length :"<<NSegs<<std::endl;
0410   return theChi2;
0411 }
0412
0413 int MuonSeedCleaner::SeedLength(const TrajectorySeed& seed) {
0414   int theHits = 0;
0415   for (auto const& recHit : seed.recHits()) {
0416     //std::cout<<"    segmet : "<<it <<std::endl;
0417     theHits += NRecHitsFromSegment(recHit);
0418   }
0419
0420   //std::cout<<" final Length :"<<NSegs<<std::endl;
0421   return theHits;
0422 }
0423
0424 GlobalPoint MuonSeedCleaner::SeedPosition(const TrajectorySeed& seed) {
0425   PTrajectoryStateOnDet pTSOD = seed.startingState();
0426   DetId SeedDetId(pTSOD.detId());
0427   const GeomDet* geoDet = theService->trackingGeometry()->idToDet(SeedDetId);
0428   TrajectoryStateOnSurface SeedTSOS =
0429       trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
0430   GlobalPoint pos = SeedTSOS.globalPosition();
0431
0432   return pos;
0433 }
0434
0435 GlobalVector MuonSeedCleaner::SeedMomentum(const TrajectorySeed& seed) {
0436   PTrajectoryStateOnDet pTSOD = seed.startingState();
0437   DetId SeedDetId(pTSOD.detId());
0438   const GeomDet* geoDet = theService->trackingGeometry()->idToDet(SeedDetId);
0439   TrajectoryStateOnSurface SeedTSOS =
0440       trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
0441   GlobalVector mom = SeedTSOS.globalMomentum();
0442
0443   return mom;
0444 }
0445
0446 int MuonSeedCleaner::NRecHitsFromSegment(const TrackingRecHit& rhit) {
0447   int NRechits = 0;
0448   const GeomDet* gdet = theService->trackingGeometry()->idToDet(rhit.geographicalId());
0449   MuonTransientTrackingRecHit::MuonRecHitPointer theSeg =
0450       MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone());
0451
0452   DetId geoId = gdet->geographicalId();
0453   if (geoId.subdetId() == MuonSubdetId::DT) {
0454     DTChamberId DT_Id(rhit.geographicalId());
0455     std::vector<TrackingRecHit*> DThits = theSeg->recHits();
0456     int dt1DHits = 0;
0457     for (size_t j = 0; j < DThits.size(); j++) {
0458       dt1DHits += (DThits[j]->recHits()).size();
0459     }
0460     NRechits = dt1DHits;
0461   }
0462
0463   if (geoId.subdetId() == MuonSubdetId::CSC) {
0464     NRechits = (theSeg->recHits()).size();
0465   }
0466   return NRechits;
0467 }
0468
0469 int MuonSeedCleaner::NRecHitsFromSegment(MuonTransientTrackingRecHit* rhit) {
0470   int NRechits = 0;
0471   DetId geoId = rhit->geographicalId();
0472   if (geoId.subdetId() == MuonSubdetId::DT) {
0473     DTChamberId DT_Id(geoId);
0474     std::vector<TrackingRecHit*> DThits = rhit->recHits();
0475     int dt1DHits = 0;
0476     for (size_t j = 0; j < DThits.size(); j++) {
0477       dt1DHits += (DThits[j]->recHits()).size();
0478     }
0479     NRechits = dt1DHits;
0480     //std::cout<<" D_rh("<< dt1DHits  <<") " ;
0481   }
0482   if (geoId.subdetId() == MuonSubdetId::CSC) {
0483     NRechits = (rhit->recHits()).size();
0484     //std::cout<<" C_rh("<<(rhit->recHits()).size() <<") " ;
0485   }
0486   return NRechits;
0487 }
0488
0489 double MuonSeedCleaner::NChi2OfSegment(const TrackingRecHit& rhit) {
0490   double NChi2 = 999999.;
0491   const GeomDet* gdet = theService->trackingGeometry()->idToDet(rhit.geographicalId());
0492   MuonTransientTrackingRecHit::MuonRecHitPointer theSeg =
0493       MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone());
0494
0495   double dof = static_cast<double>(theSeg->degreesOfFreedom());
0496   NChi2 = theSeg->chi2() / dof;
0497   //std::cout<<" Chi2 = "<< NChi2  <<" |" ;
0498
0499   return NChi2;
0500 }