Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:04

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/src/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     int idx = 0;
0297     for (auto const& r1 : seeds[i].recHits()) {
0298       idx++;
0299       const GeomDet* gdet = theService->trackingGeometry()->idToDet(r1.geographicalId());
0300       DetId geoId = gdet->geographicalId();
0301 
0302       if (geoId.subdetId() == MuonSubdetId::DT) {
0303         DTChamberId DT_Id(r1.geographicalId());
0304         //std::cout<<" ID:"<<DT_Id <<" pos:"<< r1->localPosition()  <<std::endl;
0305         if (DT_Id.station() != 1)
0306           continue;
0307         withFirstLayer = true;
0308       }
0309       if (geoId.subdetId() == MuonSubdetId::CSC) {
0310         idx++;
0311         CSCDetId CSC_Id = CSCDetId(r1.geographicalId());
0312         //std::cout<<" ID:"<<CSC_Id <<" pos:"<< r1->localPosition()  <<std::endl;
0313         if (CSC_Id.station() != 1)
0314           continue;
0315         withFirstLayer = true;
0316       }
0317     }
0318     bool goodseed = (longSeed && withFirstLayer) ? true : false;
0319 
0320     if (goodseed == good)
0321       theCandidate.push_back(seeds[i]);
0322   }
0323   return theCandidate;
0324 }
0325 
0326 std::vector<SeedContainer> MuonSeedCleaner::GroupSeeds(std::vector<TrajectorySeed>& seeds) {
0327   std::vector<SeedContainer> seedCollection;
0328   seedCollection.clear();
0329   std::vector<TrajectorySeed> theGroup;
0330   std::vector<bool> usedSeed(seeds.size(), false);
0331 
0332   // categorize seeds by comparing overlapping segments or a certian eta-phi cone
0333   for (size_t i = 0; i < seeds.size(); i++) {
0334     if (usedSeed[i])
0335       continue;
0336     theGroup.push_back(seeds[i]);
0337     usedSeed[i] = true;
0338 
0339     GlobalPoint pos1 = SeedPosition(seeds[i]);
0340 
0341     for (size_t j = i + 1; j < seeds.size(); j++) {
0342       // 1.1 seeds with overlaaping segments will be grouped together
0343       unsigned int overlapping = OverlapSegments(seeds[i], seeds[j]);
0344       if (!usedSeed[j] && overlapping > 0) {
0345         // reject the identical seeds
0346         if (seeds[i].nHits() == overlapping && seeds[j].nHits() == overlapping) {
0347           usedSeed[j] = true;
0348           continue;
0349         }
0350         theGroup.push_back(seeds[j]);
0351         usedSeed[j] = true;
0352       }
0353       if (usedSeed[j])
0354         continue;
0355 
0356       // 1.2 seeds in a certain cone are grouped together
0357       GlobalPoint pos2 = SeedPosition(seeds[j]);
0358       double dh = pos1.eta() - pos2.eta();
0359       double df = pos1.phi() - pos2.phi();
0360       double dR = sqrt((dh * dh) + (df * df));
0361 
0362       if (dR > 0.3 && seeds[j].nHits() == 1)
0363         continue;
0364       if (dR > 0.2 && seeds[j].nHits() > 1)
0365         continue;
0366       theGroup.push_back(seeds[j]);
0367       usedSeed[j] = true;
0368     }
0369     sort(theGroup.begin(), theGroup.end(), lengthSorting);
0370     seedCollection.push_back(theGroup);
0371     //std::cout<<" group "<<seedCollection.size() <<" w/"<< theGroup.size() <<" seeds"<<std::endl;
0372     theGroup.clear();
0373   }
0374   return seedCollection;
0375 }
0376 
0377 unsigned int MuonSeedCleaner::OverlapSegments(const TrajectorySeed& seed1, const TrajectorySeed& seed2) {
0378   unsigned int overlapping = 0;
0379   for (auto const& r1 : seed1.recHits()) {
0380     DetId id1 = r1.geographicalId();
0381     const GeomDet* gdet1 = theService->trackingGeometry()->idToDet(id1);
0382     GlobalPoint gp1 = gdet1->toGlobal(r1.localPosition());
0383 
0384     for (auto const& r2 : seed2.recHits()) {
0385       DetId id2 = r2.geographicalId();
0386       if (id1 != id2)
0387         continue;
0388 
0389       const GeomDet* gdet2 = theService->trackingGeometry()->idToDet(id2);
0390       GlobalPoint gp2 = gdet2->toGlobal(r2.localPosition());
0391 
0392       double dx = gp1.x() - gp2.x();
0393       double dy = gp1.y() - gp2.y();
0394       double dz = gp1.z() - gp2.z();
0395       double dL = sqrt(dx * dx + dy * dy + dz * dz);
0396 
0397       if (dL < 1.)
0398         overlapping++;
0399     }
0400   }
0401   return overlapping;
0402 }
0403 
0404 double MuonSeedCleaner::SeedChi2(const TrajectorySeed& seed) {
0405   double theChi2 = 0.;
0406   for (auto const& r1 : seed.recHits()) {
0407     //std::cout<<"    segmet : "<<it <<std::endl;
0408     theChi2 += NChi2OfSegment(r1);
0409   }
0410   theChi2 = theChi2 / seed.nHits();
0411 
0412   //std::cout<<" final Length :"<<NSegs<<std::endl;
0413   return theChi2;
0414 }
0415 
0416 int MuonSeedCleaner::SeedLength(const TrajectorySeed& seed) {
0417   int theHits = 0;
0418   for (auto const& recHit : seed.recHits()) {
0419     //std::cout<<"    segmet : "<<it <<std::endl;
0420     theHits += NRecHitsFromSegment(recHit);
0421   }
0422 
0423   //std::cout<<" final Length :"<<NSegs<<std::endl;
0424   return theHits;
0425 }
0426 
0427 GlobalPoint MuonSeedCleaner::SeedPosition(const TrajectorySeed& seed) {
0428   PTrajectoryStateOnDet pTSOD = seed.startingState();
0429   DetId SeedDetId(pTSOD.detId());
0430   const GeomDet* geoDet = theService->trackingGeometry()->idToDet(SeedDetId);
0431   TrajectoryStateOnSurface SeedTSOS =
0432       trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
0433   GlobalPoint pos = SeedTSOS.globalPosition();
0434 
0435   return pos;
0436 }
0437 
0438 GlobalVector MuonSeedCleaner::SeedMomentum(const TrajectorySeed& seed) {
0439   PTrajectoryStateOnDet pTSOD = seed.startingState();
0440   DetId SeedDetId(pTSOD.detId());
0441   const GeomDet* geoDet = theService->trackingGeometry()->idToDet(SeedDetId);
0442   TrajectoryStateOnSurface SeedTSOS =
0443       trajectoryStateTransform::transientState(pTSOD, &(geoDet->surface()), &*theService->magneticField());
0444   GlobalVector mom = SeedTSOS.globalMomentum();
0445 
0446   return mom;
0447 }
0448 
0449 int MuonSeedCleaner::NRecHitsFromSegment(const TrackingRecHit& rhit) {
0450   int NRechits = 0;
0451   const GeomDet* gdet = theService->trackingGeometry()->idToDet(rhit.geographicalId());
0452   MuonTransientTrackingRecHit::MuonRecHitPointer theSeg =
0453       MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone());
0454 
0455   DetId geoId = gdet->geographicalId();
0456   if (geoId.subdetId() == MuonSubdetId::DT) {
0457     DTChamberId DT_Id(rhit.geographicalId());
0458     std::vector<TrackingRecHit*> DThits = theSeg->recHits();
0459     int dt1DHits = 0;
0460     for (size_t j = 0; j < DThits.size(); j++) {
0461       dt1DHits += (DThits[j]->recHits()).size();
0462     }
0463     NRechits = dt1DHits;
0464   }
0465 
0466   if (geoId.subdetId() == MuonSubdetId::CSC) {
0467     NRechits = (theSeg->recHits()).size();
0468   }
0469   return NRechits;
0470 }
0471 
0472 int MuonSeedCleaner::NRecHitsFromSegment(MuonTransientTrackingRecHit* rhit) {
0473   int NRechits = 0;
0474   DetId geoId = rhit->geographicalId();
0475   if (geoId.subdetId() == MuonSubdetId::DT) {
0476     DTChamberId DT_Id(geoId);
0477     std::vector<TrackingRecHit*> DThits = rhit->recHits();
0478     int dt1DHits = 0;
0479     for (size_t j = 0; j < DThits.size(); j++) {
0480       dt1DHits += (DThits[j]->recHits()).size();
0481     }
0482     NRechits = dt1DHits;
0483     //std::cout<<" D_rh("<< dt1DHits  <<") " ;
0484   }
0485   if (geoId.subdetId() == MuonSubdetId::CSC) {
0486     NRechits = (rhit->recHits()).size();
0487     //std::cout<<" C_rh("<<(rhit->recHits()).size() <<") " ;
0488   }
0489   return NRechits;
0490 }
0491 
0492 double MuonSeedCleaner::NChi2OfSegment(const TrackingRecHit& rhit) {
0493   double NChi2 = 999999.;
0494   const GeomDet* gdet = theService->trackingGeometry()->idToDet(rhit.geographicalId());
0495   MuonTransientTrackingRecHit::MuonRecHitPointer theSeg =
0496       MuonTransientTrackingRecHit::specificBuild(gdet, rhit.clone());
0497 
0498   double dof = static_cast<double>(theSeg->degreesOfFreedom());
0499   NChi2 = theSeg->chi2() / dof;
0500   //std::cout<<" Chi2 = "<< NChi2  <<" |" ;
0501 
0502   return NChi2;
0503 }