Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:02:21 

0001 /** \class AnalyticalTrackSelector
0002  *
0003  * selects a subset of a track collection, copying extra information on demand
0004  * 
0005  * \author Paolo Azzurri, Giovanni Petrucciani 
0006  *
0007  *
0008  *
0009  */
0010 
0011 #include <utility>
0012 #include <vector>
0013 #include <memory>
0014 #include <algorithm>
0015 #include <map>
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "DataFormats/VertexReco/interface/Vertex.h"
0025 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0026 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0027 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0028 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0029 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0030 
0031 #include "MultiTrackSelector.h"
0032 
0033 using namespace reco;
0034 
0035 class dso_hidden AnalyticalTrackSelector final : public MultiTrackSelector {
0036 private:
0037 public:
0038   /// constructor
0039   explicit AnalyticalTrackSelector(const edm::ParameterSet& cfg);
0040   /// destructor
0041   ~AnalyticalTrackSelector() override;
0042 
0043 private:
0044   typedef math::XYZPoint Point;
0045   /// process one event
0046   void run(edm::Event& evt, const edm::EventSetup& es) const override;
0047 
0048   /// copy only the tracks, not extras and rechits (for AOD)
0049   bool copyExtras_;
0050   /// copy also trajectories and trajectory->track associations
0051   bool copyTrajectories_;
0052   /// eta restrictions
0053   double minEta_;
0054   double maxEta_;
0055 
0056   edm::EDGetTokenT<std::vector<Trajectory>> srcTraj_;
0057   edm::EDGetTokenT<TrajTrackAssociationCollection> srcTass_;
0058 };
0059 
0060 #include "DataFormats/Common/interface/ValueMap.h"
0061 
0062 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0063 #include <Math/DistFunc.h>
0064 #include "TMath.h"
0065 
0066 AnalyticalTrackSelector::AnalyticalTrackSelector(const edm::ParameterSet& cfg) : MultiTrackSelector() {
0067   //Spoof the pset for each track selector!
0068   //Size is always 1!!!
0069   qualityToSet_.reserve(1);
0070   vtxNumber_.reserve(1);
0071   vertexCut_.reserve(1);
0072   res_par_.reserve(1);
0073   chi2n_par_.reserve(1);
0074   chi2n_no1Dmod_par_.reserve(1);
0075   d0_par1_.reserve(1);
0076   dz_par1_.reserve(1);
0077   d0_par2_.reserve(1);
0078   dz_par2_.reserve(1);
0079   applyAdaptedPVCuts_.reserve(1);
0080   max_d0_.reserve(1);
0081   max_z0_.reserve(1);
0082   nSigmaZ_.reserve(1);
0083   min_layers_.reserve(1);
0084   min_3Dlayers_.reserve(1);
0085   max_lostLayers_.reserve(1);
0086   min_hits_bypass_.reserve(1);
0087   applyAbsCutsIfNoPV_.reserve(1);
0088   max_d0NoPV_.reserve(1);
0089   max_z0NoPV_.reserve(1);
0090   preFilter_.reserve(1);
0091   max_relpterr_.reserve(1);
0092   min_nhits_.reserve(1);
0093   max_minMissHitOutOrIn_.reserve(1);
0094   max_lostHitFraction_.reserve(1);
0095   min_eta_.reserve(1);
0096   max_eta_.reserve(1);
0097   forest_.reserve(1);
0098   mvaType_.reserve(1);
0099   useMVA_.reserve(1);
0100 
0101   produces<edm::ValueMap<float>>("MVAVals");
0102   //foward compatibility
0103   produces<MVACollection>("MVAValues");
0104   useAnyMVA_ = false;
0105   forest_[0] = nullptr;
0106   if (cfg.exists("useAnyMVA"))
0107     useAnyMVA_ = cfg.getParameter<bool>("useAnyMVA");
0108 
0109   src_ = consumes<reco::TrackCollection>(cfg.getParameter<edm::InputTag>("src"));
0110   hSrc_ = consumes<TrackingRecHitCollection>(cfg.getParameter<edm::InputTag>("src"));
0111   beamspot_ = consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamspot"));
0112   useVertices_ = cfg.getParameter<bool>("useVertices");
0113   useVtxError_ = cfg.getParameter<bool>("useVtxError");
0114   if (useVertices_)
0115     vertices_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("vertices"));
0116   copyExtras_ = cfg.getUntrackedParameter<bool>("copyExtras", false);
0117   copyTrajectories_ = cfg.getUntrackedParameter<bool>("copyTrajectories", false);
0118   if (copyTrajectories_) {
0119     srcTraj_ = consumes<std::vector<Trajectory>>(cfg.getParameter<edm::InputTag>("src"));
0120     srcTass_ = consumes<TrajTrackAssociationCollection>(cfg.getParameter<edm::InputTag>("src"));
0121   }
0122 
0123   qualityToSet_.push_back(TrackBase::undefQuality);
0124   // parameters for vertex selection
0125   vtxNumber_.push_back(useVertices_ ? cfg.getParameter<int32_t>("vtxNumber") : 0);
0126   vertexCut_.push_back(useVertices_ ? cfg.getParameter<std::string>("vertexCut") : "");
0127   //  parameters for adapted optimal cuts on chi2 and primary vertex compatibility
0128   res_par_.push_back(cfg.getParameter<std::vector<double>>("res_par"));
0129   chi2n_par_.push_back(cfg.getParameter<double>("chi2n_par"));
0130   chi2n_no1Dmod_par_.push_back(cfg.getParameter<double>("chi2n_no1Dmod_par"));
0131   d0_par1_.push_back(cfg.getParameter<std::vector<double>>("d0_par1"));
0132   dz_par1_.push_back(cfg.getParameter<std::vector<double>>("dz_par1"));
0133   d0_par2_.push_back(cfg.getParameter<std::vector<double>>("d0_par2"));
0134   dz_par2_.push_back(cfg.getParameter<std::vector<double>>("dz_par2"));
0135 
0136   // Boolean indicating if adapted primary vertex compatibility cuts are to be applied.
0137   applyAdaptedPVCuts_.push_back(cfg.getParameter<bool>("applyAdaptedPVCuts"));
0138   // Impact parameter absolute cuts.
0139   max_d0_.push_back(cfg.getParameter<double>("max_d0"));
0140   max_z0_.push_back(cfg.getParameter<double>("max_z0"));
0141   nSigmaZ_.push_back(cfg.getParameter<double>("nSigmaZ"));
0142   // Cuts on numbers of layers with hits/3D hits/lost hits.
0143   min_layers_.push_back(cfg.getParameter<uint32_t>("minNumberLayers"));
0144   min_3Dlayers_.push_back(cfg.getParameter<uint32_t>("minNumber3DLayers"));
0145   max_lostLayers_.push_back(cfg.getParameter<uint32_t>("maxNumberLostLayers"));
0146   min_hits_bypass_.push_back(cfg.getParameter<uint32_t>("minHitsToBypassChecks"));
0147   max_relpterr_.push_back(cfg.getParameter<double>("max_relpterr"));
0148   min_nhits_.push_back(cfg.getParameter<uint32_t>("min_nhits"));
0149   max_minMissHitOutOrIn_.push_back(
0150       cfg.existsAs<int32_t>("max_minMissHitOutOrIn") ? cfg.getParameter<int32_t>("max_minMissHitOutOrIn") : 99);
0151   max_lostHitFraction_.push_back(
0152       cfg.existsAs<double>("max_lostHitFraction") ? cfg.getParameter<double>("max_lostHitFraction") : 1.0);
0153   min_eta_.push_back(cfg.getParameter<double>("min_eta"));
0154   max_eta_.push_back(cfg.getParameter<double>("max_eta"));
0155 
0156   // Flag to apply absolute cuts if no PV passes the selection
0157   applyAbsCutsIfNoPV_.push_back(cfg.getParameter<bool>("applyAbsCutsIfNoPV"));
0158   keepAllTracks_.push_back(cfg.exists("keepAllTracks") ? cfg.getParameter<bool>("keepAllTracks") : false);
0159 
0160   setQualityBit_.push_back(false);
0161   std::string qualityStr = cfg.getParameter<std::string>("qualityBit");
0162 
0163   if (d0_par1_[0].size() != 2 || dz_par1_[0].size() != 2 || d0_par2_[0].size() != 2 || dz_par2_[0].size() != 2) {
0164     edm::LogError("MisConfiguration") << "vector of size less then 2";
0165     throw;
0166   }
0167 
0168   if (cfg.exists("qualityBit")) {
0169     std::string qualityStr = cfg.getParameter<std::string>("qualityBit");
0170     if (!qualityStr.empty()) {
0171       setQualityBit_[0] = true;
0172       qualityToSet_[0] = TrackBase::qualityByName(cfg.getParameter<std::string>("qualityBit"));
0173     }
0174   }
0175 
0176   if (keepAllTracks_[0] && !setQualityBit_[0])
0177     throw cms::Exception("Configuration")
0178         << "If you set 'keepAllTracks' to true, you must specify which qualityBit to set.\n";
0179   if (setQualityBit_[0] && (qualityToSet_[0] == TrackBase::undefQuality))
0180     throw cms::Exception("Configuration")
0181         << "You can't set the quality bit " << cfg.getParameter<std::string>("qualityBit")
0182         << " as it is 'undefQuality' or unknown.\n";
0183   if (applyAbsCutsIfNoPV_[0]) {
0184     max_d0NoPV_.push_back(cfg.getParameter<double>("max_d0NoPV"));
0185     max_z0NoPV_.push_back(cfg.getParameter<double>("max_z0NoPV"));
0186   } else {  //dummy values
0187     max_d0NoPV_.push_back(0.);
0188     max_z0NoPV_.push_back(0.);
0189   }
0190 
0191   if (useAnyMVA_) {
0192     bool thisMVA = false;
0193     if (cfg.exists("useMVA"))
0194       thisMVA = cfg.getParameter<bool>("useMVA");
0195     useMVA_.push_back(thisMVA);
0196     if (thisMVA) {
0197       double minVal = -1;
0198       if (cfg.exists("minMVA"))
0199         minVal = cfg.getParameter<double>("minMVA");
0200       min_MVA_.push_back(minVal);
0201       mvaType_.push_back(cfg.exists("mvaType") ? cfg.getParameter<std::string>("mvaType") : "Detached");
0202       forestLabel_.push_back(cfg.exists("GBRForestLabel") ? cfg.getParameter<std::string>("GBRForestLabel")
0203                                                           : "MVASelectorIter0");
0204       useMVAonly_.push_back(cfg.exists("useMVAonly") ? cfg.getParameter<bool>("useMVAonly") : false);
0205     } else {
0206       min_MVA_.push_back(-9999.0);
0207       useMVAonly_.push_back(false);
0208       mvaType_.push_back("Detached");
0209       forestLabel_.push_back("MVASelectorIter0");
0210     }
0211   } else {
0212     useMVA_.push_back(false);
0213     useMVAonly_.push_back(false);
0214     min_MVA_.push_back(-9999.0);
0215     mvaType_.push_back("Detached");
0216     forestLabel_.push_back("MVASelectorIter0");
0217   }
0218 
0219   std::string alias(cfg.getParameter<std::string>("@module_label"));
0220   if (copyExtras_) {
0221     produces<reco::TrackExtraCollection>().setBranchAlias(alias + "TrackExtras");
0222     produces<TrackingRecHitCollection>().setBranchAlias(alias + "RecHits");
0223   }
0224   if (copyTrajectories_) {
0225     produces<std::vector<Trajectory>>().setBranchAlias(alias + "Trajectories");
0226     produces<TrajTrackAssociationCollection>().setBranchAlias(alias + "TrajectoryTrackAssociations");
0227   }
0228   // TrackCollection refers to TrackingRechit and TrackExtra
0229   // collections, need to declare its production after them to work
0230   // around a rare race condition in framework scheduling
0231   produces<reco::TrackCollection>().setBranchAlias(alias + "Tracks");
0232 }
0233 
0234 AnalyticalTrackSelector::~AnalyticalTrackSelector() {}
0235 
0236 void AnalyticalTrackSelector::run(edm::Event& evt, const edm::EventSetup& es) const {
0237   // storage....
0238   std::unique_ptr<reco::TrackCollection> selTracks_;
0239   std::unique_ptr<reco::TrackExtraCollection> selTrackExtras_;
0240   std::unique_ptr<TrackingRecHitCollection> selHits_;
0241   std::unique_ptr<std::vector<Trajectory>> selTrajs_;
0242   std::unique_ptr<std::vector<const Trajectory*>> selTrajPtrs_;
0243   std::unique_ptr<TrajTrackAssociationCollection> selTTAss_;
0244   reco::TrackRefProd rTracks_;
0245   reco::TrackExtraRefProd rTrackExtras_;
0246   TrackingRecHitRefProd rHits_;
0247   edm::RefProd<std::vector<Trajectory>> rTrajectories_;
0248   std::vector<reco::TrackRef> trackRefs_;
0249 
0250   using namespace std;
0251   using namespace edm;
0252   using namespace reco;
0253 
0254   Handle<TrackCollection> hSrcTrack;
0255   Handle<vector<Trajectory>> hTraj;
0256   Handle<vector<Trajectory>> hTrajP;
0257   Handle<TrajTrackAssociationCollection> hTTAss;
0258 
0259   // looking for the beam spot
0260   edm::Handle<reco::BeamSpot> hBsp;
0261   evt.getByToken(beamspot_, hBsp);
0262   reco::BeamSpot vertexBeamSpot;
0263   vertexBeamSpot = *hBsp;
0264 
0265   // Select good primary vertices for use in subsequent track selection
0266   const reco::VertexCollection dummyVtx;
0267   const reco::VertexCollection* vtxPtr = &dummyVtx;
0268   std::vector<Point> points;
0269   std::vector<float> vterr, vzerr;
0270   if (useVertices_) {
0271     vtxPtr = &evt.get(vertices_);
0272     selectVertices(0, *vtxPtr, points, vterr, vzerr);
0273     // Debug
0274     LogDebug("SelectVertex") << points.size() << " good pixel vertices";
0275   }
0276 
0277   // Get tracks
0278   evt.getByToken(src_, hSrcTrack);
0279   // get hits in track..
0280   Handle<TrackingRecHitCollection> hSrcHits;
0281   evt.getByToken(hSrc_, hSrcHits);
0282   const TrackingRecHitCollection& srcHits(*hSrcHits);
0283 
0284   selTracks_ = std::make_unique<TrackCollection>();
0285   rTracks_ = evt.getRefBeforePut<TrackCollection>();
0286   if (copyExtras_) {
0287     selTrackExtras_ = std::make_unique<TrackExtraCollection>();
0288     selHits_ = std::make_unique<TrackingRecHitCollection>();
0289     rHits_ = evt.getRefBeforePut<TrackingRecHitCollection>();
0290     rTrackExtras_ = evt.getRefBeforePut<TrackExtraCollection>();
0291   }
0292 
0293   if (copyTrajectories_)
0294     trackRefs_.resize(hSrcTrack->size());
0295 
0296   std::vector<float> mvaVals_(hSrcTrack->size(), -99.f);
0297   processMVA(evt, es, vertexBeamSpot, *vtxPtr, 0, mvaVals_, true);
0298 
0299   // Loop over tracks
0300   size_t current = 0;
0301   for (TrackCollection::const_iterator it = hSrcTrack->begin(), ed = hSrcTrack->end(); it != ed; ++it, ++current) {
0302     const Track& trk = *it;
0303     // Check if this track passes cuts
0304 
0305     LogTrace("TrackSelection") << "ready to check track with pt=" << trk.pt();
0306 
0307     float mvaVal = 0;
0308     if (useAnyMVA_)
0309       mvaVal = mvaVals_[current];
0310     bool ok = select(0, vertexBeamSpot, srcHits, trk, points, vterr, vzerr, mvaVal);
0311     if (!ok) {
0312       LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " NOT selected";
0313 
0314       if (copyTrajectories_)
0315         trackRefs_[current] = reco::TrackRef();
0316       if (!keepAllTracks_[0])
0317         continue;
0318     }
0319     LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " selected";
0320     selTracks_->push_back(Track(trk));  // clone and store
0321     if (ok && setQualityBit_[0]) {
0322       selTracks_->back().setQuality(qualityToSet_[0]);
0323       if (qualityToSet_[0] == TrackBase::tight) {
0324         selTracks_->back().setQuality(TrackBase::loose);
0325       } else if (qualityToSet_[0] == TrackBase::highPurity) {
0326         selTracks_->back().setQuality(TrackBase::loose);
0327         selTracks_->back().setQuality(TrackBase::tight);
0328       }
0329       if (!points.empty()) {
0330         if (qualityToSet_[0] == TrackBase::loose) {
0331           selTracks_->back().setQuality(TrackBase::looseSetWithPV);
0332         } else if (qualityToSet_[0] == TrackBase::highPurity) {
0333           selTracks_->back().setQuality(TrackBase::looseSetWithPV);
0334           selTracks_->back().setQuality(TrackBase::highPuritySetWithPV);
0335         }
0336       }
0337     }
0338     if (copyExtras_) {
0339       // TrackExtras
0340       selTrackExtras_->push_back(TrackExtra(trk.outerPosition(),
0341                                             trk.outerMomentum(),
0342                                             trk.outerOk(),
0343                                             trk.innerPosition(),
0344                                             trk.innerMomentum(),
0345                                             trk.innerOk(),
0346                                             trk.outerStateCovariance(),
0347                                             trk.outerDetId(),
0348                                             trk.innerStateCovariance(),
0349                                             trk.innerDetId(),
0350                                             trk.seedDirection(),
0351                                             trk.seedRef()));
0352       selTracks_->back().setExtra(TrackExtraRef(rTrackExtras_, selTrackExtras_->size() - 1));
0353       TrackExtra& tx = selTrackExtras_->back();
0354       tx.setResiduals(trk.residuals());
0355       // TrackingRecHits
0356       auto const firstHitIndex = selHits_->size();
0357       for (trackingRecHit_iterator hit = trk.recHitsBegin(); hit != trk.recHitsEnd(); ++hit) {
0358         selHits_->push_back((*hit)->clone());
0359       }
0360       tx.setHits(rHits_, firstHitIndex, selHits_->size() - firstHitIndex);
0361       tx.setTrajParams(trk.extra()->trajParams(), trk.extra()->chi2sX5());
0362     }
0363     if (copyTrajectories_) {
0364       trackRefs_[current] = TrackRef(rTracks_, selTracks_->size() - 1);
0365     }
0366   }
0367   if (copyTrajectories_) {
0368     Handle<vector<Trajectory>> hTraj;
0369     Handle<TrajTrackAssociationCollection> hTTAss;
0370     evt.getByToken(srcTass_, hTTAss);
0371     evt.getByToken(srcTraj_, hTraj);
0372     selTrajs_ = std::make_unique<std::vector<Trajectory>>();
0373     rTrajectories_ = evt.getRefBeforePut<vector<Trajectory>>();
0374     selTTAss_ = std::make_unique<TrajTrackAssociationCollection>(rTrajectories_, rTracks_);
0375     for (size_t i = 0, n = hTraj->size(); i < n; ++i) {
0376       Ref<vector<Trajectory>> trajRef(hTraj, i);
0377       TrajTrackAssociationCollection::const_iterator match = hTTAss->find(trajRef);
0378       if (match != hTTAss->end()) {
0379         const Ref<TrackCollection>& trkRef = match->val;
0380         short oldKey = static_cast<short>(trkRef.key());
0381         if (trackRefs_[oldKey].isNonnull()) {
0382           selTrajs_->push_back(Trajectory(*trajRef));
0383           selTTAss_->insert(Ref<vector<Trajectory>>(rTrajectories_, selTrajs_->size() - 1), trackRefs_[oldKey]);
0384         }
0385       }
0386     }
0387   }
0388 
0389   evt.put(std::move(selTracks_));
0390   if (copyExtras_) {
0391     evt.put(std::move(selTrackExtras_));
0392     evt.put(std::move(selHits_));
0393   }
0394   if (copyTrajectories_) {
0395     evt.put(std::move(selTrajs_));
0396     evt.put(std::move(selTTAss_));
0397   }
0398 }
0399 
0400 #include "FWCore/PluginManager/interface/ModuleDef.h"
0401 #include "FWCore/Framework/interface/MakerMacros.h"
0402 
0403 DEFINE_FWK_MODULE(AnalyticalTrackSelector);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!