Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:22:16 

0001 /** \class DuplicateTrackMerger
0002  * 
0003  * selects pairs of tracks that should be single tracks
0004  *
0005  * \author Matthew Walker
0006  */
0007 
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Utilities/interface/ESGetToken.h"
0015 
0016 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 
0019 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0020 
0021 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0022 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0023 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0024 #include "TrackMerger.h"
0025 
0026 #include "CommonTools/Utils/interface/DynArray.h"
0027 #include <vector>
0028 #include <algorithm>
0029 #include <string>
0030 #include <iostream>
0031 #include <atomic>
0032 
0033 #include "CondFormats/GBRForest/interface/GBRForest.h"
0034 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
0035 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0036 
0037 using namespace reco;
0038 namespace {
0039   class DuplicateTrackMerger final : public edm::stream::EDProducer<> {
0040   public:
0041     /// constructor
0042     explicit DuplicateTrackMerger(const edm::ParameterSet &iPara);
0043     /// destructor
0044     ~DuplicateTrackMerger() override;
0045 
0046     using CandidateToDuplicate = std::vector<std::pair<int, int>>;
0047 
0048     static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0049 
0050   private:
0051     /// produce one event
0052     void produce(edm::Event &, const edm::EventSetup &) override;
0053 
0054     bool checkForDisjointTracks(const reco::Track *t1, const reco::Track *t2, TSCPBuilderNoMaterial &tscpBuilder) const;
0055     bool checkForOverlappingTracks(
0056         const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const;
0057 
0058   private:
0059     /// MVA discriminator
0060     const GBRForest *forest_;
0061 
0062     /// MVA weights file
0063     std::string dbFileName_;
0064     bool useForestFromDB_;
0065     std::string forestLabel_;
0066 
0067     std::string propagatorName_;
0068     std::string chi2EstimatorName_;
0069 
0070     /// track input collection
0071     edm::EDGetTokenT<reco::TrackCollection> trackSource_;
0072     /// minDeltaR3d cut value
0073     double minDeltaR3d2_;
0074     /// minBDTG cut value
0075     double minBDTG_;
0076     ///min pT cut value
0077     double minpT2_;
0078     ///min p cut value
0079     double minP_;
0080     ///max distance between two tracks at closest approach
0081     float maxDCA2_;
0082     ///max difference in phi between two tracks
0083     float maxDPhi_;
0084     ///max difference in Lambda between two tracks
0085     float maxDLambda_;
0086     ///max difference in transverse impact parameter between two tracks
0087     float maxDdxy_;
0088     ///max difference in longitudinal impact parameter between two tracks
0089     float maxDdsz_;
0090     ///max difference in q/p between two tracks
0091     float maxDQoP_;
0092     /// max number of hits for shorter track for the overlap check
0093     unsigned int overlapCheckMaxHits_;
0094     /// max number of missing layers for the overlap check
0095     unsigned int overlapCheckMaxMissingLayers_;
0096     /// min cosT for the overlap check
0097     double overlapCheckMinCosT_;
0098 
0099     const MagneticField *magfield_;
0100     const TrackerTopology *ttopo_;
0101     const TrackerGeometry *geom_;
0102     const Propagator *propagator_;
0103     const Chi2MeasurementEstimatorBase *chi2Estimator_;
0104 
0105     edm::ESGetToken<GBRForest, GBRWrapperRcd> forestToken_;
0106     edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0107     edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
0108     edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geometryToken_;
0109     edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
0110     edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> estimatorToken_;
0111 
0112     ///Merger
0113     TrackMerger merger_;
0114 
0115 #ifdef EDM_ML_DEBUG
0116     bool debug_;
0117 #endif
0118   };
0119 }  // namespace
0120 
0121 #include "FWCore/Framework/interface/Event.h"
0122 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0123 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0124 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0125 #include "FWCore/Framework/interface/EventSetup.h"
0126 #include "FWCore/Framework/interface/ESHandle.h"
0127 #include "TFile.h"
0128 
0129 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0130 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0131 
0132 #include "DuplicateTrackType.h"
0133 
0134 namespace {
0135 
0136   void DuplicateTrackMerger::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0137     edm::ParameterSetDescription desc;
0138     desc.add<edm::InputTag>("source", edm::InputTag());
0139     desc.add<double>("minDeltaR3d", -4.0);
0140     desc.add<double>("minBDTG", -0.1);
0141     desc.add<double>("minpT", 0.2);
0142     desc.add<double>("minP", 0.4);
0143     desc.add<double>("maxDCA", 30.0);
0144     desc.add<double>("maxDPhi", 0.30);
0145     desc.add<double>("maxDLambda", 0.30);
0146     desc.add<double>("maxDdsz", 10.0);
0147     desc.add<double>("maxDdxy", 10.0);
0148     desc.add<double>("maxDQoP", 0.25);
0149     desc.add<unsigned>("overlapCheckMaxHits", 4);
0150     desc.add<unsigned>("overlapCheckMaxMissingLayers", 1);
0151     desc.add<double>("overlapCheckMinCosT", 0.99);
0152     desc.add<std::string>("forestLabel", "MVADuplicate");
0153     desc.add<std::string>("GBRForestFileName", "");
0154     desc.add<bool>("useInnermostState", true);
0155     desc.add<std::string>("ttrhBuilderName", "WithAngleAndTemplate");
0156     desc.add<std::string>("propagatorName", "PropagatorWithMaterial");
0157     desc.add<std::string>("chi2EstimatorName", "DuplicateTrackMergerChi2Est");
0158     descriptions.add("DuplicateTrackMerger", desc);
0159   }
0160 
0161   DuplicateTrackMerger::DuplicateTrackMerger(const edm::ParameterSet &iPara)
0162       : forest_(nullptr), merger_(iPara, consumesCollector()) {
0163     trackSource_ = consumes<reco::TrackCollection>(iPara.getParameter<edm::InputTag>("source"));
0164     minDeltaR3d2_ = iPara.getParameter<double>("minDeltaR3d");
0165     minDeltaR3d2_ *= std::abs(minDeltaR3d2_);
0166     minBDTG_ = iPara.getParameter<double>("minBDTG");
0167     minpT2_ = iPara.getParameter<double>("minpT");
0168     minpT2_ *= minpT2_;
0169     minP_ = iPara.getParameter<double>("minP");
0170     maxDCA2_ = iPara.getParameter<double>("maxDCA");
0171     maxDCA2_ *= maxDCA2_;
0172     maxDPhi_ = iPara.getParameter<double>("maxDPhi");
0173     maxDLambda_ = iPara.getParameter<double>("maxDLambda");
0174     maxDdsz_ = iPara.getParameter<double>("maxDdsz");
0175     maxDdxy_ = iPara.getParameter<double>("maxDdxy");
0176     maxDQoP_ = iPara.getParameter<double>("maxDQoP");
0177     overlapCheckMaxHits_ = iPara.getParameter<unsigned>("overlapCheckMaxHits");
0178     overlapCheckMaxMissingLayers_ = iPara.getParameter<unsigned>("overlapCheckMaxMissingLayers");
0179     overlapCheckMinCosT_ = iPara.getParameter<double>("overlapCheckMinCosT");
0180 
0181     produces<std::vector<TrackCandidate>>("candidates");
0182     produces<CandidateToDuplicate>("candidateMap");
0183 
0184     forestLabel_ = iPara.getParameter<std::string>("forestLabel");
0185 
0186     dbFileName_ = iPara.getParameter<std::string>("GBRForestFileName");
0187     useForestFromDB_ = dbFileName_.empty();
0188 
0189     propagatorName_ = iPara.getParameter<std::string>("propagatorName");
0190     chi2EstimatorName_ = iPara.getParameter<std::string>("chi2EstimatorName");
0191     if (useForestFromDB_) {
0192       forestToken_ = esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", forestLabel_));
0193     }
0194     magFieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0195     trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0196     geometryToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0197     propagatorToken_ = esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", propagatorName_));
0198     estimatorToken_ =
0199         esConsumes<Chi2MeasurementEstimatorBase, TrackingComponentsRecord>(edm::ESInputTag("", chi2EstimatorName_));
0200 
0201     /*
0202   tmvaReader_ = new TMVA::Reader("!Color:Silent");
0203   tmvaReader_->AddVariable("ddsz",&tmva_ddsz_);
0204   tmvaReader_->AddVariable("ddxy",&tmva_ddxy_);
0205   tmvaReader_->AddVariable("dphi",&tmva_dphi_);
0206   tmvaReader_->AddVariable("dlambda",&tmva_dlambda_);
0207   tmvaReader_->AddVariable("dqoverp",&tmva_dqoverp_);
0208   tmvaReader_->AddVariable("delta3d_r",&tmva_d3dr_);
0209   tmvaReader_->AddVariable("delta3d_z",&tmva_d3dz_);
0210   tmvaReader_->AddVariable("outer_nMissingInner",&tmva_outer_nMissingInner_);
0211   tmvaReader_->AddVariable("inner_nMissingOuter",&tmva_inner_nMissingOuter_);
0212   tmvaReader_->BookMVA("BDTG",mvaFilePath);
0213   */
0214   }
0215 
0216   DuplicateTrackMerger::~DuplicateTrackMerger() {
0217     if (!useForestFromDB_)
0218       delete forest_;
0219   }
0220 
0221 #ifdef VI_STAT
0222   struct Stat {
0223     Stat() : maxCos(1.1), nCand(0), nLoop0(0) {}
0224     ~Stat() { std::cout << "Stats " << nCand << ' ' << nLoop0 << ' ' << maxCos << std::endl; }
0225     std::atomic<float> maxCos;
0226     std::atomic<int> nCand, nLoop0;
0227   };
0228   Stat stat;
0229 #endif
0230 
0231   template <typename T>
0232   void update_maximum(std::atomic<T> &maximum_value, T const &value) noexcept {
0233     T prev_value = maximum_value;
0234     while (prev_value < value && !maximum_value.compare_exchange_weak(prev_value, value))
0235       ;
0236   }
0237 
0238   template <typename T>
0239   void update_minimum(std::atomic<T> &minimum_value, T const &value) noexcept {
0240     T prev_value = minimum_value;
0241     while (prev_value > value && !minimum_value.compare_exchange_weak(prev_value, value))
0242       ;
0243   }
0244 
0245   void DuplicateTrackMerger::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0246     merger_.init(iSetup);
0247 
0248     if (!forest_) {
0249       if (useForestFromDB_) {
0250         edm::ESHandle<GBRForest> forestHandle = iSetup.getHandle(forestToken_);
0251         forest_ = forestHandle.product();
0252       } else {
0253         TFile gbrfile(dbFileName_.c_str());
0254         forest_ = dynamic_cast<const GBRForest *>(gbrfile.Get(forestLabel_.c_str()));
0255       }
0256     }
0257 
0258     //edm::Handle<edm::View<reco::Track> >handle;
0259     edm::Handle<reco::TrackCollection> handle;
0260     iEvent.getByToken(trackSource_, handle);
0261     auto const &tracks = *handle;
0262 
0263     edm::ESHandle<MagneticField> hmagfield = iSetup.getHandle(magFieldToken_);
0264     magfield_ = hmagfield.product();
0265 
0266     edm::ESHandle<TrackerTopology> httopo = iSetup.getHandle(trackerTopoToken_);
0267     ttopo_ = httopo.product();
0268 
0269     edm::ESHandle<TrackerGeometry> hgeom = iSetup.getHandle(geometryToken_);
0270     geom_ = hgeom.product();
0271 
0272     edm::ESHandle<Propagator> hpropagator = iSetup.getHandle(propagatorToken_);
0273     propagator_ = hpropagator.product();
0274 
0275     edm::ESHandle<Chi2MeasurementEstimatorBase> hestimator = iSetup.getHandle(estimatorToken_);
0276     chi2Estimator_ = hestimator.product();
0277 
0278     TSCPBuilderNoMaterial tscpBuilder;
0279     auto out_duplicateCandidates = std::make_unique<std::vector<TrackCandidate>>();
0280 
0281     auto out_candidateMap = std::make_unique<CandidateToDuplicate>();
0282     LogDebug("DuplicateTrackMerger") << "Number of tracks to be checked for merging: " << tracks.size();
0283 
0284 #ifdef EDM_ML_DEBUG
0285     auto test = [&](const reco::Track *a, const reco::Track *b) {
0286       const auto ev = iEvent.id().event();
0287       const auto aOriAlgo = a->originalAlgo();
0288       const auto bOriAlgo = b->originalAlgo();
0289       const auto aSeed = a->seedRef().key();
0290       const auto bSeed = b->seedRef().key();
0291       return ((ev == 6903 && ((aOriAlgo == 23 && aSeed == 695 && bOriAlgo == 5 && bSeed == 652) ||
0292                               (aOriAlgo == 23 && aSeed == 400 && bOriAlgo == 7 && bSeed == 156) ||
0293                               (aOriAlgo == 4 && aSeed == 914 && bOriAlgo == 22 && bSeed == 503) ||
0294                               (aOriAlgo == 5 && aSeed == 809 && bOriAlgo == 4 && bSeed == 1030) ||
0295                               (aOriAlgo == 23 && aSeed == 749 && bOriAlgo == 5 && bSeed == 659) ||
0296                               (aOriAlgo == 4 && aSeed == 1053 && bOriAlgo == 23 && bSeed == 1035) ||
0297                               (aOriAlgo == 4 && aSeed == 810 && bOriAlgo == 5 && bSeed == 666) ||
0298                               (aOriAlgo == 4 && aSeed == 974 && bOriAlgo == 5 && bSeed == 778))) ||
0299               (ev == 6904 && ((aOriAlgo == 23 && aSeed == 526 && bOriAlgo == 5 && bSeed == 307) ||
0300                               (aOriAlgo == 4 && aSeed == 559 && bOriAlgo == 22 && bSeed == 472))) ||
0301               (ev == 6902 && ((aOriAlgo == 4 && aSeed == 750 && bOriAlgo == 22 && bSeed == 340) ||
0302                               (aOriAlgo == 4 && aSeed == 906 && bOriAlgo == 5 && bSeed == 609) ||
0303                               (aOriAlgo == 4 && aSeed == 724 && bOriAlgo == 5 && bSeed == 528) ||
0304                               (aOriAlgo == 4 && aSeed == 943 && bOriAlgo == 23 && bSeed == 739) ||
0305                               (aOriAlgo == 8 && aSeed == 2 && bOriAlgo == 9 && bSeed == 2282) ||
0306                               (aOriAlgo == 23 && aSeed == 827 && bOriAlgo == 5 && bSeed == 656) ||
0307                               (aOriAlgo == 22 && aSeed == 667 && bOriAlgo == 7 && bSeed == 516))));
0308     };
0309 #endif
0310 
0311     // cache few "heavy to compute quantities
0312     int nTracks = 0;
0313     declareDynArray(const reco::Track *, tracks.size(), selTracks);
0314     declareDynArray(unsigned int, tracks.size(), nValidHits);
0315     declareDynArray(unsigned int, tracks.size(), oriIndex);
0316     for (auto i = 0U; i < tracks.size(); i++) {
0317       const reco::Track *rt1 = &tracks[i];
0318       if (rt1->innerMomentum().perp2() < minpT2_)
0319         continue;
0320       selTracks[nTracks] = rt1;
0321       nValidHits[nTracks] = rt1->numberOfValidHits();  // yes it is extremely heavy!
0322       oriIndex[nTracks] = i;
0323       ++nTracks;
0324     }
0325 
0326     for (int i = 0; i < nTracks; i++) {
0327       const reco::Track *rt1 = selTracks[i];
0328       for (int j = i + 1; j < nTracks; j++) {
0329         const reco::Track *rt2 = selTracks[j];
0330 
0331 #ifdef EDM_ML_DEBUG
0332         debug_ = false;
0333         if (test(rt1, rt2) || test(rt2, rt1)) {
0334           debug_ = true;
0335           LogTrace("DuplicateTrackMerger")
0336               << "Track1 " << i << " originalAlgo " << rt1->originalAlgo() << " seed " << rt1->seedRef().key() << " pT "
0337               << std::sqrt(rt1->innerMomentum().perp2()) << " charge " << rt1->charge() << " outerPosition2 "
0338               << rt1->outerPosition().perp2() << "\n"
0339               << "Track2 " << j << " originalAlgo " << rt2->originalAlgo() << " seed " << rt2->seedRef().key() << " pT "
0340               << std::sqrt(rt2->innerMomentum().perp2()) << " charge " << rt2->charge() << " outerPosition2 "
0341               << rt2->outerPosition().perp2();
0342         }
0343 #endif
0344 
0345         if (rt1->charge() != rt2->charge())
0346           continue;
0347         auto cosT = (*rt1).momentum().Dot((*rt2).momentum());  // not normalized!
0348         IfLogTrace(debug_, "DuplicateTrackMerger") << " cosT " << cosT;
0349         if (cosT < 0.)
0350           continue;
0351         cosT /= std::sqrt((*rt1).momentum().Mag2() * (*rt2).momentum().Mag2());
0352 
0353         const reco::Track *t1, *t2;
0354         unsigned int nhv1, nhv2;
0355         if (rt1->outerPosition().perp2() < rt2->outerPosition().perp2()) {
0356           t1 = rt1;
0357           nhv1 = nValidHits[i];
0358           t2 = rt2;
0359           nhv2 = nValidHits[j];
0360         } else {
0361           t1 = rt2;
0362           nhv1 = nValidHits[j];
0363           t2 = rt1;
0364           nhv2 = nValidHits[i];
0365         }
0366         auto deltaR3d2 = (t1->outerPosition() - t2->innerPosition()).mag2();
0367 
0368         if (t1->outerPosition().perp2() > t2->innerPosition().perp2())
0369           deltaR3d2 *= -1.0;
0370         IfLogTrace(debug_, "DuplicateTrackMerger")
0371             << " deltaR3d2 " << deltaR3d2 << " t1.outerPos2 " << t1->outerPosition().perp2() << " t2.innerPos2 "
0372             << t2->innerPosition().perp2();
0373 
0374         bool compatible = false;
0375         DuplicateTrackType duplType;
0376         if (deltaR3d2 >= minDeltaR3d2_) {
0377           compatible = checkForDisjointTracks(t1, t2, tscpBuilder);
0378           duplType = DuplicateTrackType::Disjoint;
0379         } else {
0380           compatible = checkForOverlappingTracks(t1, t2, nhv1, nhv2, cosT);
0381           duplType = DuplicateTrackType::Overlapping;
0382         }
0383         if (!compatible)
0384           continue;
0385 
0386         IfLogTrace(debug_, "DuplicateTrackMerger") << " marking as duplicates" << oriIndex[i] << ',' << oriIndex[j];
0387         out_duplicateCandidates->push_back(merger_.merge(*t1, *t2, duplType));
0388         out_candidateMap->emplace_back(oriIndex[i], oriIndex[j]);
0389 
0390 #ifdef VI_STAT
0391         ++stat.nCand;
0392         //    auto cosT = float((*t1).momentum().unit().Dot((*t2).momentum().unit()));
0393         if (cosT > 0)
0394           update_minimum(stat.maxCos, float(cosT));
0395         else
0396           ++stat.nLoop0;
0397 #endif
0398       }
0399     }
0400     iEvent.put(std::move(out_duplicateCandidates), "candidates");
0401     iEvent.put(std::move(out_candidateMap), "candidateMap");
0402   }
0403 
0404   bool DuplicateTrackMerger::checkForDisjointTracks(const reco::Track *t1,
0405                                                     const reco::Track *t2,
0406                                                     TSCPBuilderNoMaterial &tscpBuilder) const {
0407     IfLogTrace(debug_, "DuplicateTrackMerger") << " Checking for disjoint duplicates";
0408 
0409     FreeTrajectoryState fts1 = trajectoryStateTransform::outerFreeState(*t1, &*magfield_, false);
0410     FreeTrajectoryState fts2 = trajectoryStateTransform::innerFreeState(*t2, &*magfield_, false);
0411     GlobalPoint avgPoint((t1->outerPosition().x() + t2->innerPosition().x()) * 0.5,
0412                          (t1->outerPosition().y() + t2->innerPosition().y()) * 0.5,
0413                          (t1->outerPosition().z() + t2->innerPosition().z()) * 0.5);
0414     TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
0415     TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
0416     IfLogTrace(debug_, "DuplicateTrackMerger")
0417         << " TSCP1.isValid " << TSCP1.isValid() << " TSCP2.isValid " << TSCP2.isValid();
0418     if (!TSCP1.isValid())
0419       return false;
0420     if (!TSCP2.isValid())
0421       return false;
0422 
0423     const FreeTrajectoryState &ftsn1 = TSCP1.theState();
0424     const FreeTrajectoryState &ftsn2 = TSCP2.theState();
0425 
0426     IfLogTrace(debug_, "DuplicateTrackMerger") << " DCA2 " << (ftsn2.position() - ftsn1.position()).mag2();
0427     if ((ftsn2.position() - ftsn1.position()).mag2() > maxDCA2_)
0428       return false;
0429 
0430     auto qoverp1 = ftsn1.signedInverseMomentum();
0431     auto qoverp2 = ftsn2.signedInverseMomentum();
0432     float tmva_dqoverp_ = qoverp1 - qoverp2;
0433     IfLogTrace(debug_, "DuplicateTrackMerger") << " dqoverp " << tmva_dqoverp_;
0434     if (std::abs(tmva_dqoverp_) > maxDQoP_)
0435       return false;
0436 
0437     //auto pp = [&](TrajectoryStateClosestToPoint const & ts) { std::cout << ' ' << ts.perigeeParameters().vector()[0] << '/'  << ts.perigeeError().transverseCurvatureError();};
0438     //if(qoverp1*qoverp2 <0) { std::cout << "charge different " << qoverp1 <<',' << qoverp2; pp(TSCP1); pp(TSCP2); std::cout << std::endl;}
0439 
0440     // lambda = pi/2 - theta....  so l1-l2 == t2-t1
0441     float tmva_dlambda_ = ftsn2.momentum().theta() - ftsn1.momentum().theta();
0442     IfLogTrace(debug_, "DuplicateTrackMerger") << " dlambda " << tmva_dlambda_;
0443     if (std::abs(tmva_dlambda_) > maxDLambda_)
0444       return false;
0445 
0446     auto phi1 = ftsn1.momentum().phi();
0447     auto phi2 = ftsn2.momentum().phi();
0448     float tmva_dphi_ = phi1 - phi2;
0449     if (std::abs(tmva_dphi_) > float(M_PI))
0450       tmva_dphi_ = 2.f * float(M_PI) - std::abs(tmva_dphi_);
0451     IfLogTrace(debug_, "DuplicateTrackMerger") << " dphi " << tmva_dphi_;
0452     if (std::abs(tmva_dphi_) > maxDPhi_)
0453       return false;
0454 
0455     auto dxy1 =
0456         (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) / TSCP1.pt();
0457     auto dxy2 =
0458         (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) / TSCP2.pt();
0459     float tmva_ddxy_ = dxy1 - dxy2;
0460     IfLogTrace(debug_, "DuplicateTrackMerger") << " ddxy " << tmva_ddxy_;
0461     if (std::abs(tmva_ddxy_) > maxDdxy_)
0462       return false;
0463 
0464     auto dsz1 = ftsn1.position().z() * TSCP1.pt() / TSCP1.momentum().mag() -
0465                 (ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) /
0466                     TSCP1.pt() * ftsn1.momentum().z() / ftsn1.momentum().mag();
0467     auto dsz2 = ftsn2.position().z() * TSCP2.pt() / TSCP2.momentum().mag() -
0468                 (ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) /
0469                     TSCP2.pt() * ftsn2.momentum().z() / ftsn2.momentum().mag();
0470     float tmva_ddsz_ = dsz1 - dsz2;
0471     IfLogTrace(debug_, "DuplicateTrackMerger") << " ddsz " << tmva_ddsz_;
0472     if (std::abs(tmva_ddsz_) > maxDdsz_)
0473       return false;
0474 
0475     float tmva_d3dr_ = avgPoint.perp();
0476     float tmva_d3dz_ = avgPoint.z();
0477     float tmva_outer_nMissingInner_ = t2->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0478     float tmva_inner_nMissingOuter_ = t1->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
0479 
0480     float gbrVals_[9];
0481     gbrVals_[0] = tmva_ddsz_;
0482     gbrVals_[1] = tmva_ddxy_;
0483     gbrVals_[2] = tmva_dphi_;
0484     gbrVals_[3] = tmva_dlambda_;
0485     gbrVals_[4] = tmva_dqoverp_;
0486     gbrVals_[5] = tmva_d3dr_;
0487     gbrVals_[6] = tmva_d3dz_;
0488     gbrVals_[7] = tmva_outer_nMissingInner_;
0489     gbrVals_[8] = tmva_inner_nMissingOuter_;
0490 
0491     auto mvaBDTG = forest_->GetClassifier(gbrVals_);
0492     IfLogTrace(debug_, "DuplicateTrackMerger") << " mvaBDTG " << mvaBDTG;
0493     if (mvaBDTG < minBDTG_)
0494       return false;
0495 
0496     //  std::cout << "to merge " << mvaBDTG << ' ' << std::copysign(std::sqrt(std::abs(deltaR3d2)),deltaR3d2) << ' ' << tmva_dphi_ << ' ' << TSCP1.pt() <<'/'<<TSCP2.pt() << std::endl;
0497     return true;
0498   }
0499 
0500   bool DuplicateTrackMerger::checkForOverlappingTracks(
0501       const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const {
0502     // ensure t1 is the shorter track
0503     if (nvh2 < nvh1) {
0504       std::swap(t1, t2);
0505       std::swap(nvh1, nvh2);
0506     }
0507 
0508     IfLogTrace(debug_, "DuplicateTrackMerger")
0509         << " Checking for overlapping duplicates, cosT " << cosT << " t1 hits " << nvh1;
0510     if (cosT < overlapCheckMinCosT_)
0511       return false;
0512     if (nvh1 > overlapCheckMaxHits_)
0513       return false;
0514 
0515     // find the hit on the longer track on layer of the first hit of the shorter track
0516     auto findHitOnT2 = [&](const TrackingRecHit *hit1) {
0517       const auto hitDet = hit1->geographicalId().det();
0518       const auto hitSubdet = hit1->geographicalId().subdetId();
0519       const auto hitLayer = ttopo_->layer(hit1->geographicalId());
0520       return std::find_if(t2->recHitsBegin(), t2->recHitsEnd(), [&](const TrackingRecHit *hit2) {
0521         const auto &detId = hit2->geographicalId();
0522         return (detId.det() == hitDet && detId.subdetId() == hitSubdet && ttopo_->layer(detId) == hitLayer);
0523       });
0524     };
0525 
0526     auto t1HitIter = t1->recHitsBegin();
0527     if (!(*t1HitIter)->isValid()) {
0528       IfLogTrace(debug_, "DuplicateTrackMerger") << " first t1 hit invalid";
0529       return false;
0530     }
0531     auto t2HitIter = findHitOnT2(*t1HitIter);
0532     if (t2HitIter == t2->recHitsEnd()) {
0533       // if first hit not found, try with second
0534       // if that fails, then reject
0535       ++t1HitIter;
0536       assert(t1HitIter != t1->recHitsEnd());
0537 
0538       if (!(*t1HitIter)->isValid()) {
0539         IfLogTrace(debug_, "DuplicateTrackMerger") << " second t1 hit invalid";
0540         return false;
0541       }
0542       t2HitIter = findHitOnT2(*t1HitIter);
0543       if (t2HitIter == t2->recHitsEnd())
0544         return false;
0545     }
0546     IfLogTrace(debug_, "DuplicateTrackMerger")
0547         << " starting overlap check from t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0548         << std::distance(t2->recHitsBegin(), t2HitIter);
0549 
0550     auto prevSubdet = (*t1HitIter)->geographicalId().subdetId();
0551     const TrajectoryStateOnSurface tsosInner = trajectoryStateTransform::innerStateOnSurface(*t2, *geom_, magfield_);
0552 
0553     ++t1HitIter;
0554     ++t2HitIter;
0555     unsigned int missedLayers = 0;
0556     while (t1HitIter != t1->recHitsEnd() && t2HitIter != t2->recHitsEnd()) {
0557       // in case of invalid hits, reject immediately
0558       if ((*t1HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t1HitIter) ||
0559           (*t2HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t2HitIter)) {
0560         IfLogTrace(debug_, "DuplicateTrackMerger")
0561             << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0562             << std::distance(t2->recHitsBegin(), t2HitIter) << " either is invalid, types t1 "
0563             << (*t1HitIter)->getType() << " t2 " << (*t2HitIter)->getType();
0564         return false;
0565       }
0566 
0567       const auto &t1DetId = (*t1HitIter)->geographicalId();
0568       const auto &t2DetId = (*t2HitIter)->geographicalId();
0569 
0570       const auto t1Det = t1DetId.det();
0571       const auto t2Det = t2DetId.det();
0572       if (t1Det != DetId::Tracker || t2Det != DetId::Tracker) {
0573         IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
0574                                                    << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
0575                                                    << " either not from Tracker, dets t1 " << t1Det << " t2 " << t2Det;
0576         return false;
0577       }
0578 
0579       const auto t1Subdet = t1DetId.subdetId();
0580       const auto t1Layer = ttopo_->layer(t1DetId);
0581 
0582       // reject if hits have the same DetId but are different
0583       if (t1DetId == t2DetId) {
0584         if (!(*t1HitIter)->sharesInput(*t2HitIter, TrackingRecHit::all)) {
0585           IfLogTrace(debug_, "DuplicateTrackMerger")
0586               << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0587               << std::distance(t2->recHitsBegin(), t2HitIter) << " same DetId (" << t1DetId.rawId()
0588               << ") but do not share all input";
0589           return false;
0590         }
0591       } else {
0592         const auto t2Subdet = t2DetId.subdetId();
0593         const auto t2Layer = ttopo_->layer(t2DetId);
0594 
0595         // reject if hits are on different layers
0596         if (t1Subdet != t2Subdet || t1Layer != t2Layer) {
0597           bool recovered = false;
0598           // but try to recover first by checking if either one has skipped over a layer
0599           if (t1Subdet == prevSubdet && t2Subdet != prevSubdet) {
0600             // t1 has a layer t2 doesn't
0601             ++t1HitIter;
0602             recovered = true;
0603           } else if (t1Subdet != prevSubdet && t2Subdet == prevSubdet) {
0604             // t2 has a layer t1 doesn't
0605             ++t2HitIter;
0606             recovered = true;
0607           } else if (t1Subdet == t2Subdet) {
0608             prevSubdet = t1Subdet;
0609             // same subdet, so layer must be different
0610             if (t2Layer > t1Layer) {
0611               // t1 has a layer t2 doesn't
0612               ++t1HitIter;
0613               recovered = true;
0614             } else if (t1Layer > t2Layer) {
0615               // t2 has a layer t1 doesn't
0616               ++t2HitIter;
0617               recovered = true;
0618             }
0619           }
0620           if (recovered) {
0621             ++missedLayers;
0622             if (missedLayers > overlapCheckMaxMissingLayers_) {
0623               IfLogTrace(debug_, "DuplicateTrackMerger") << " max number of missed layers exceeded";
0624               return false;
0625             }
0626             continue;
0627           }
0628 
0629           IfLogTrace(debug_, "DuplicateTrackMerger")
0630               << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0631               << std::distance(t2->recHitsBegin(), t2HitIter) << " are on different layers (subdet, layer) t1 "
0632               << t1Subdet << "," << t1Layer << " t2 " << t2Subdet << "," << t2Layer;
0633           return false;
0634         }
0635         // reject if same layer (but not same hit) in non-pixel detector
0636         else if (t1Subdet != PixelSubdetector::PixelBarrel && t1Subdet != PixelSubdetector::PixelEndcap) {
0637           IfLogTrace(debug_, "DuplicateTrackMerger")
0638               << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0639               << std::distance(t2->recHitsBegin(), t2HitIter) << " are on same layer, but in non-pixel detector (det "
0640               << t1Det << " subdet " << t1Subdet << " layer " << t1Layer << ")";
0641           return false;
0642         }
0643       }
0644 
0645       // Propagate longer track to the shorter track hit surface, check compatibility
0646       TrajectoryStateOnSurface tsosPropagated = propagator_->propagate(tsosInner, (*t1HitIter)->det()->surface());
0647       if (!tsosPropagated.isValid()) {  // reject immediately if TSOS is not valid
0648         IfLogTrace(debug_, "DuplicateTrackMerger")
0649             << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0650             << std::distance(t2->recHitsBegin(), t2HitIter) << " TSOS not valid";
0651         return false;
0652       }
0653       auto passChi2Pair = chi2Estimator_->estimate(tsosPropagated, **t1HitIter);
0654       if (!passChi2Pair.first) {
0655         IfLogTrace(debug_, "DuplicateTrackMerger")
0656             << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0657             << std::distance(t2->recHitsBegin(), t2HitIter) << " hit chi2 compatibility failed with chi2 "
0658             << passChi2Pair.second;
0659         return false;
0660       }
0661 
0662       prevSubdet = t1Subdet;
0663       ++t1HitIter;
0664       ++t2HitIter;
0665     }
0666     if (t1HitIter != t1->recHitsEnd()) {
0667       IfLogTrace(debug_, "DuplicateTrackMerger") << " hits on t2 ended before hits on t1";
0668       return false;
0669     }
0670 
0671     IfLogTrace(debug_, "DuplicateTrackMerger") << " all hits on t2 are on layers whre t1 has also a hit";
0672     return true;
0673   }
0674 }  // namespace
0675 
0676 #include "FWCore/PluginManager/interface/ModuleDef.h"
0677 #include "FWCore/Framework/interface/MakerMacros.h"
0678 
0679 DEFINE_FWK_MODULE(DuplicateTrackMerger);
[ 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!