TrackerSeedGenerator/src/TrackerSeedCleaner.cc

0001 /*
0002  * \class TrackerSeedCleaner
0003  *  Reference class for seeds cleaning
0004  *  Seeds Cleaner based on sharedHits cleaning, direction cleaning and pt cleaning
0005     \author A. Grelli -  Purdue University, Pavia University
0006  */
0007
0008 #include "RecoMuon/TrackerSeedGenerator/interface/TrackerSeedCleaner.h"
0009
0010 //---------------
0011 // C++ Headers --
0012 //---------------
0013 #include <vector>
0014
0015 //-------------------------------
0016 // Collaborating Class Headers --
0017 //-------------------------------
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020
0021 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0022 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0023 #include "DataFormats/Math/interface/deltaPhi.h"
0024
0025 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0026 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0027 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0028 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0029 #include "RecoTracker/TkTrackingRegions/interface/TkTrackingRegionsMargin.h"
0030 #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h"
0031
0032 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0033
0034 using namespace std;
0035 using namespace edm;
0036
0037 //
0038 // inizialization
0039 //
0040 void TrackerSeedCleaner::init(const MuonServiceProxy* service) { theProxyService = service; }
0041
0042 //
0043 //
0044 //
0045 void TrackerSeedCleaner::setEvent(const edm::Event& event) { event.getByToken(beamspotToken_, bsHandle_); }
0046
0047 //
0048 // clean seeds
0049 //
0050 void TrackerSeedCleaner::clean(const reco::TrackRef& muR,
0051                                const RectangularEtaPhiTrackingRegion& region,
0052                                tkSeeds& seeds) {
0053   // call the shared input cleaner
0054   if (cleanBySharedHits) {
0055     seeds = nonRedundantSeeds(seeds);
0056   }
0057
0058   theTTRHBuilder = theProxyService->eventSetup().getHandle(theTTRHBuilderToken);
0059
0060   LogDebug("TrackerSeedCleaner") << seeds.size() << " trajectory seeds to the events before cleaning" << endl;
0061
0062   //check the validity otherwise vertexing
0063   const reco::BeamSpot& bs = *bsHandle_;
0064   /*reco track and seeds as arguments. Seeds eta and phi are checked and
0065    based on deviation from L2 eta and phi seed is accepted or not*/
0066
0067   std::vector<TrajectorySeed> result;
0068
0069   TSCBLBuilderNoMaterial tscblBuilder;
0070   // PerigeeConversions tspConverter;
0071   for (TrajectorySeedCollection::iterator seed = seeds.begin(); seed < seeds.end(); ++seed) {
0072     if (seed->nHits() < 2)
0073       continue;
0074     //get parameters and errors from the seed state
0075     TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(&*(seed->recHits().end() - 1));
0076     TrajectoryStateOnSurface state = trajectoryStateTransform::transientState(
0077         seed->startingState(), recHit->surface(), theProxyService->magneticField().product());
0078
0079     TrajectoryStateClosestToBeamLine tsAtClosestApproachSeed =
0080         tscblBuilder(*state.freeState(), bs);  //as in TrackProducerAlgorithms
0081     if (!tsAtClosestApproachSeed.isValid())
0082       continue;
0083     GlobalPoint vSeed1 = tsAtClosestApproachSeed.trackStateAtPCA().position();
0084     GlobalVector pSeed = tsAtClosestApproachSeed.trackStateAtPCA().momentum();
0085     GlobalPoint vSeed(vSeed1.x() - bs.x0(), vSeed1.y() - bs.y0(), vSeed1.z() - bs.z0());
0086
0087     //eta,phi info from seeds
0088     double etaSeed = state.globalMomentum().eta();
0089     double phiSeed = pSeed.phi();
0090
0091     //if the limits are too stringent rescale limits
0092     typedef PixelRecoRange<float> Range;
0093     typedef TkTrackingRegionsMargin<float> Margin;
0094
0095     Range etaRange = region.etaRange();
0096     double etaLimit = (fabs(fabs(etaRange.max()) - fabs(etaRange.mean())) < 0.1)
0097                           ? 0.1
0098                           : fabs(fabs(etaRange.max()) - fabs(etaRange.mean()));
0099
0100     Margin phiMargin = region.phiMargin();
0101     double phiLimit = (phiMargin.right() < 0.1) ? 0.1 : phiMargin.right();
0102
0103     double ptSeed = pSeed.perp();
0104     double ptMin = (region.ptMin() > 3.5) ? 3.5 : region.ptMin();
0105     // Clean
0106     bool inEtaRange = etaSeed >= (etaRange.mean() - etaLimit) && etaSeed <= (etaRange.mean() + etaLimit);
0107     bool inPhiRange = (fabs(deltaPhi(phiSeed, double(region.direction().phi()))) < phiLimit);
0108     // pt cleaner
0109     bool inPtRange = ptSeed >= ptMin && ptSeed <= 2 * (muR->pt());
0110
0111     // save efficiency don't clean triplets with pt cleaner
0112     if (seed->nHits() == 3)
0113       inPtRange = true;
0114
0115     // use pt and angle cleaners
0116     if (inPtRange && usePt_Cleaner && !useDirection_Cleaner) {
0117       result.push_back(*seed);
0118       LogDebug("TrackerSeedCleaner") << " Keeping the seed : this seed passed pt selection";
0119     }
0120
0121     // use only angle default option
0122     if (inEtaRange && inPhiRange && !usePt_Cleaner && useDirection_Cleaner) {
0123       result.push_back(*seed);
0124       LogDebug("TrackerSeedCleaner") << " Keeping the seed : this seed passed direction selection";
0125     }
0126
0127     // use all the cleaners
0128     if (inEtaRange && inPhiRange && inPtRange && usePt_Cleaner && useDirection_Cleaner) {
0129       result.push_back(*seed);
0130       LogDebug("TrackerSeedCleaner") << " Keeping the seed : this seed passed direction and pt selection";
0131     }
0132
0133     LogDebug("TrackerSeedCleaner") << " eta for current seed " << etaSeed << "\n"
0134                                    << " phi for current seed " << phiSeed << "\n"
0135                                    << " eta for L2 track  " << muR->eta() << "\n"
0136                                    << " phi for L2 track  " << muR->phi() << "\n";
0137   }
0138
0139   //the new seeds collection
0140   if (!result.empty() && (useDirection_Cleaner || usePt_Cleaner))
0141     seeds.swap(result);
0142
0143   LogDebug("TrackerSeedCleaner") << seeds.size() << " trajectory seeds to the events after cleaning" << endl;
0144
0145   return;
0146 }
0147
0148 TrackerSeedCleaner::tkSeeds TrackerSeedCleaner::nonRedundantSeeds(tkSeeds const& seeds) const {
0149   std::vector<uint> idxTriplets{}, idxNonTriplets{};
0150   idxTriplets.reserve(seeds.size());
0151   idxNonTriplets.reserve(seeds.size());
0152
0153   for (uint i1 = 0; i1 < seeds.size(); ++i1) {
0154     auto const& s1 = seeds[i1];
0155     if (s1.nHits() == 3)
0156       idxTriplets.emplace_back(i1);
0157     else
0158       idxNonTriplets.emplace_back(i1);
0159   }
0160
0161   if (idxTriplets.empty()) {
0162     return seeds;
0163   }
0164
0165   std::vector<bool> keepSeedFlags(seeds.size(), true);
0166   for (uint j1 = 0; j1 < idxNonTriplets.size(); ++j1) {
0167     auto const i1 = idxNonTriplets[j1];
0168     auto const& seed = seeds[i1];
0169     keepSeedFlags[i1] = seedIsNotRedundant(seeds, seed, idxTriplets);
0170   }
0171
0172   tkSeeds result{};
0173   result.reserve(seeds.size());
0174
0175   for (uint i1 = 0; i1 < seeds.size(); ++i1) {
0176     if (keepSeedFlags[i1]) {
0177       result.emplace_back(seeds[i1]);
0178     }
0179   }
0180
0181   return result;
0182 }
0183
0184 bool TrackerSeedCleaner::seedIsNotRedundant(tkSeeds const& seeds,
0185                                             TrajectorySeed const& s1,
0186                                             std::vector<uint> const& otherIdxs) const {
0187   auto const& rh1s = s1.recHits();
0188   for (uint j2 = 0; j2 < otherIdxs.size(); ++j2) {
0189     auto const& s2 = seeds[otherIdxs[j2]];
0190     // number of shared hits
0191     uint shared = 0;
0192     for (auto const& h2 : s2.recHits()) {
0193       for (auto const& h1 : rh1s) {
0194         if (h2.sharesInput(&h1, TrackingRecHit::all)) {
0195           ++shared;
0196         }
0197       }
0198     }
0199     if (shared == 2) {
0200       return false;
0201     }
0202   }
0203
0204   return true;
0205 }