Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-02-28 23:34:30 

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 aSeedRef = a->seedRef();
0290       if (!aSeedRef) {
0291         LogDebug("DuplicateTrackMerger") << "aSeedRef is null!";
0292         return false;
0293       }
0294       const auto aSeed = a->seedRef().key();
0295       const auto bSeedRef = b->seedRef();
0296       if (!bSeedRef) {
0297         LogDebug("DuplicateTrackMerger") << "bSeedRef is null!";
0298         return false;
0299       }
0300       const auto bSeed = b->seedRef().key();
0301       return ((ev == 6903 && ((aOriAlgo == 23 && aSeed == 695 && bOriAlgo == 5 && bSeed == 652) ||
0302                               (aOriAlgo == 23 && aSeed == 400 && bOriAlgo == 7 && bSeed == 156) ||
0303                               (aOriAlgo == 4 && aSeed == 914 && bOriAlgo == 22 && bSeed == 503) ||
0304                               (aOriAlgo == 5 && aSeed == 809 && bOriAlgo == 4 && bSeed == 1030) ||
0305                               (aOriAlgo == 23 && aSeed == 749 && bOriAlgo == 5 && bSeed == 659) ||
0306                               (aOriAlgo == 4 && aSeed == 1053 && bOriAlgo == 23 && bSeed == 1035) ||
0307                               (aOriAlgo == 4 && aSeed == 810 && bOriAlgo == 5 && bSeed == 666) ||
0308                               (aOriAlgo == 4 && aSeed == 974 && bOriAlgo == 5 && bSeed == 778))) ||
0309               (ev == 6904 && ((aOriAlgo == 23 && aSeed == 526 && bOriAlgo == 5 && bSeed == 307) ||
0310                               (aOriAlgo == 4 && aSeed == 559 && bOriAlgo == 22 && bSeed == 472))) ||
0311               (ev == 6902 && ((aOriAlgo == 4 && aSeed == 750 && bOriAlgo == 22 && bSeed == 340) ||
0312                               (aOriAlgo == 4 && aSeed == 906 && bOriAlgo == 5 && bSeed == 609) ||
0313                               (aOriAlgo == 4 && aSeed == 724 && bOriAlgo == 5 && bSeed == 528) ||
0314                               (aOriAlgo == 4 && aSeed == 943 && bOriAlgo == 23 && bSeed == 739) ||
0315                               (aOriAlgo == 8 && aSeed == 2 && bOriAlgo == 9 && bSeed == 2282) ||
0316                               (aOriAlgo == 23 && aSeed == 827 && bOriAlgo == 5 && bSeed == 656) ||
0317                               (aOriAlgo == 22 && aSeed == 667 && bOriAlgo == 7 && bSeed == 516))));
0318     };
0319 #endif
0320 
0321     // cache few "heavy to compute quantities
0322     int nTracks = 0;
0323     declareDynArray(const reco::Track *, tracks.size(), selTracks);
0324     declareDynArray(unsigned int, tracks.size(), nValidHits);
0325     declareDynArray(unsigned int, tracks.size(), oriIndex);
0326     for (auto i = 0U; i < tracks.size(); i++) {
0327       const reco::Track *rt1 = &tracks[i];
0328       if (rt1->innerMomentum().perp2() < minpT2_)
0329         continue;
0330       selTracks[nTracks] = rt1;
0331       nValidHits[nTracks] = rt1->numberOfValidHits();  // yes it is extremely heavy!
0332       oriIndex[nTracks] = i;
0333       ++nTracks;
0334     }
0335 
0336     for (int i = 0; i < nTracks; i++) {
0337       const reco::Track *rt1 = selTracks[i];
0338       for (int j = i + 1; j < nTracks; j++) {
0339         const reco::Track *rt2 = selTracks[j];
0340 
0341 #ifdef EDM_ML_DEBUG
0342         debug_ = false;
0343         if (test(rt1, rt2) || test(rt2, rt1)) {
0344           debug_ = true;
0345           LogTrace("DuplicateTrackMerger")
0346               << "Track1 " << i << " originalAlgo " << rt1->originalAlgo() << " seed " << rt1->seedRef().key() << " pT "
0347               << std::sqrt(rt1->innerMomentum().perp2()) << " charge " << rt1->charge() << " outerPosition2 "
0348               << rt1->outerPosition().perp2() << "\n"
0349               << "Track2 " << j << " originalAlgo " << rt2->originalAlgo() << " seed " << rt2->seedRef().key() << " pT "
0350               << std::sqrt(rt2->innerMomentum().perp2()) << " charge " << rt2->charge() << " outerPosition2 "
0351               << rt2->outerPosition().perp2();
0352         }
0353 #endif
0354 
0355         if (rt1->charge() != rt2->charge())
0356           continue;
0357         auto cosT = (*rt1).momentum().Dot((*rt2).momentum());  // not normalized!
0358         IfLogTrace(debug_, "DuplicateTrackMerger") << " cosT " << cosT;
0359         if (cosT < 0.)
0360           continue;
0361         cosT /= std::sqrt((*rt1).momentum().Mag2() * (*rt2).momentum().Mag2());
0362 
0363         const reco::Track *t1, *t2;
0364         unsigned int nhv1, nhv2;
0365         if (rt1->outerPosition().perp2() < rt2->outerPosition().perp2()) {
0366           t1 = rt1;
0367           nhv1 = nValidHits[i];
0368           t2 = rt2;
0369           nhv2 = nValidHits[j];
0370         } else {
0371           t1 = rt2;
0372           nhv1 = nValidHits[j];
0373           t2 = rt1;
0374           nhv2 = nValidHits[i];
0375         }
0376         auto deltaR3d2 = (t1->outerPosition() - t2->innerPosition()).mag2();
0377 
0378         if (t1->outerPosition().perp2() > t2->innerPosition().perp2())
0379           deltaR3d2 *= -1.0;
0380         IfLogTrace(debug_, "DuplicateTrackMerger")
0381             << " deltaR3d2 " << deltaR3d2 << " t1.outerPos2 " << t1->outerPosition().perp2() << " t2.innerPos2 "
0382             << t2->innerPosition().perp2();
0383 
0384         bool compatible = false;
0385         DuplicateTrackType duplType;
0386         if (deltaR3d2 >= minDeltaR3d2_) {
0387           compatible = checkForDisjointTracks(t1, t2, tscpBuilder);
0388           duplType = DuplicateTrackType::Disjoint;
0389         } else {
0390           compatible = checkForOverlappingTracks(t1, t2, nhv1, nhv2, cosT);
0391           duplType = DuplicateTrackType::Overlapping;
0392         }
0393         if (!compatible)
0394           continue;
0395 
0396         IfLogTrace(debug_, "DuplicateTrackMerger") << " marking as duplicates" << oriIndex[i] << ',' << oriIndex[j];
0397         out_duplicateCandidates->push_back(merger_.merge(*t1, *t2, duplType));
0398         out_candidateMap->emplace_back(oriIndex[i], oriIndex[j]);
0399 
0400 #ifdef VI_STAT
0401         ++stat.nCand;
0402         //    auto cosT = float((*t1).momentum().unit().Dot((*t2).momentum().unit()));
0403         if (cosT > 0)
0404           update_minimum(stat.maxCos, float(cosT));
0405         else
0406           ++stat.nLoop0;
0407 #endif
0408       }
0409     }
0410     iEvent.put(std::move(out_duplicateCandidates), "candidates");
0411     iEvent.put(std::move(out_candidateMap), "candidateMap");
0412   }
0413 
0414   bool DuplicateTrackMerger::checkForDisjointTracks(const reco::Track *t1,
0415                                                     const reco::Track *t2,
0416                                                     TSCPBuilderNoMaterial &tscpBuilder) const {
0417     IfLogTrace(debug_, "DuplicateTrackMerger") << " Checking for disjoint duplicates";
0418 
0419     FreeTrajectoryState fts1 = trajectoryStateTransform::outerFreeState(*t1, &*magfield_, false);
0420     FreeTrajectoryState fts2 = trajectoryStateTransform::innerFreeState(*t2, &*magfield_, false);
0421     GlobalPoint avgPoint((t1->outerPosition().x() + t2->innerPosition().x()) * 0.5,
0422                          (t1->outerPosition().y() + t2->innerPosition().y()) * 0.5,
0423                          (t1->outerPosition().z() + t2->innerPosition().z()) * 0.5);
0424     TrajectoryStateClosestToPoint TSCP1 = tscpBuilder(fts1, avgPoint);
0425     TrajectoryStateClosestToPoint TSCP2 = tscpBuilder(fts2, avgPoint);
0426     IfLogTrace(debug_, "DuplicateTrackMerger")
0427         << " TSCP1.isValid " << TSCP1.isValid() << " TSCP2.isValid " << TSCP2.isValid();
0428     if (!TSCP1.isValid())
0429       return false;
0430     if (!TSCP2.isValid())
0431       return false;
0432 
0433     const FreeTrajectoryState &ftsn1 = TSCP1.theState();
0434     const FreeTrajectoryState &ftsn2 = TSCP2.theState();
0435 
0436     IfLogTrace(debug_, "DuplicateTrackMerger") << " DCA2 " << (ftsn2.position() - ftsn1.position()).mag2();
0437     if ((ftsn2.position() - ftsn1.position()).mag2() > maxDCA2_)
0438       return false;
0439 
0440     auto qoverp1 = ftsn1.signedInverseMomentum();
0441     auto qoverp2 = ftsn2.signedInverseMomentum();
0442     float tmva_dqoverp_ = qoverp1 - qoverp2;
0443     IfLogTrace(debug_, "DuplicateTrackMerger") << " dqoverp " << tmva_dqoverp_;
0444     if (std::abs(tmva_dqoverp_) > maxDQoP_)
0445       return false;
0446 
0447     //auto pp = [&](TrajectoryStateClosestToPoint const & ts) { std::cout << ' ' << ts.perigeeParameters().vector()[0] << '/'  << ts.perigeeError().transverseCurvatureError();};
0448     //if(qoverp1*qoverp2 <0) { std::cout << "charge different " << qoverp1 <<',' << qoverp2; pp(TSCP1); pp(TSCP2); std::cout << std::endl;}
0449 
0450     // lambda = pi/2 - theta....  so l1-l2 == t2-t1
0451     float tmva_dlambda_ = ftsn2.momentum().theta() - ftsn1.momentum().theta();
0452     IfLogTrace(debug_, "DuplicateTrackMerger") << " dlambda " << tmva_dlambda_;
0453     if (std::abs(tmva_dlambda_) > maxDLambda_)
0454       return false;
0455 
0456     auto phi1 = ftsn1.momentum().phi();
0457     auto phi2 = ftsn2.momentum().phi();
0458     float tmva_dphi_ = phi1 - phi2;
0459     if (std::abs(tmva_dphi_) > float(M_PI))
0460       tmva_dphi_ = 2.f * float(M_PI) - std::abs(tmva_dphi_);
0461     IfLogTrace(debug_, "DuplicateTrackMerger") << " dphi " << tmva_dphi_;
0462     if (std::abs(tmva_dphi_) > maxDPhi_)
0463       return false;
0464 
0465     auto dxy1 =
0466         (-ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) / TSCP1.pt();
0467     auto dxy2 =
0468         (-ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) / TSCP2.pt();
0469     float tmva_ddxy_ = dxy1 - dxy2;
0470     IfLogTrace(debug_, "DuplicateTrackMerger") << " ddxy " << tmva_ddxy_;
0471     if (std::abs(tmva_ddxy_) > maxDdxy_)
0472       return false;
0473 
0474     auto dsz1 = ftsn1.position().z() * TSCP1.pt() / TSCP1.momentum().mag() -
0475                 (ftsn1.position().x() * ftsn1.momentum().y() + ftsn1.position().y() * ftsn1.momentum().x()) /
0476                     TSCP1.pt() * ftsn1.momentum().z() / ftsn1.momentum().mag();
0477     auto dsz2 = ftsn2.position().z() * TSCP2.pt() / TSCP2.momentum().mag() -
0478                 (ftsn2.position().x() * ftsn2.momentum().y() + ftsn2.position().y() * ftsn2.momentum().x()) /
0479                     TSCP2.pt() * ftsn2.momentum().z() / ftsn2.momentum().mag();
0480     float tmva_ddsz_ = dsz1 - dsz2;
0481     IfLogTrace(debug_, "DuplicateTrackMerger") << " ddsz " << tmva_ddsz_;
0482     if (std::abs(tmva_ddsz_) > maxDdsz_)
0483       return false;
0484 
0485     float tmva_d3dr_ = avgPoint.perp();
0486     float tmva_d3dz_ = avgPoint.z();
0487     float tmva_outer_nMissingInner_ = t2->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0488     float tmva_inner_nMissingOuter_ = t1->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
0489 
0490     float gbrVals_[9];
0491     gbrVals_[0] = tmva_ddsz_;
0492     gbrVals_[1] = tmva_ddxy_;
0493     gbrVals_[2] = tmva_dphi_;
0494     gbrVals_[3] = tmva_dlambda_;
0495     gbrVals_[4] = tmva_dqoverp_;
0496     gbrVals_[5] = tmva_d3dr_;
0497     gbrVals_[6] = tmva_d3dz_;
0498     gbrVals_[7] = tmva_outer_nMissingInner_;
0499     gbrVals_[8] = tmva_inner_nMissingOuter_;
0500 
0501     auto mvaBDTG = forest_->GetClassifier(gbrVals_);
0502     IfLogTrace(debug_, "DuplicateTrackMerger") << " mvaBDTG " << mvaBDTG;
0503     if (mvaBDTG < minBDTG_)
0504       return false;
0505 
0506     //  std::cout << "to merge " << mvaBDTG << ' ' << std::copysign(std::sqrt(std::abs(deltaR3d2)),deltaR3d2) << ' ' << tmva_dphi_ << ' ' << TSCP1.pt() <<'/'<<TSCP2.pt() << std::endl;
0507     return true;
0508   }
0509 
0510   bool DuplicateTrackMerger::checkForOverlappingTracks(
0511       const reco::Track *t1, const reco::Track *t2, unsigned int nvh1, unsigned int nvh2, double cosT) const {
0512     // ensure t1 is the shorter track
0513     if (nvh2 < nvh1) {
0514       std::swap(t1, t2);
0515       std::swap(nvh1, nvh2);
0516     }
0517 
0518     IfLogTrace(debug_, "DuplicateTrackMerger")
0519         << " Checking for overlapping duplicates, cosT " << cosT << " t1 hits " << nvh1;
0520     if (cosT < overlapCheckMinCosT_)
0521       return false;
0522     if (nvh1 > overlapCheckMaxHits_)
0523       return false;
0524 
0525     // find the hit on the longer track on layer of the first hit of the shorter track
0526     auto findHitOnT2 = [&](const TrackingRecHit *hit1) {
0527       const auto hitDet = hit1->geographicalId().det();
0528       const auto hitSubdet = hit1->geographicalId().subdetId();
0529       const auto hitLayer = ttopo_->layer(hit1->geographicalId());
0530       return std::find_if(t2->recHitsBegin(), t2->recHitsEnd(), [&](const TrackingRecHit *hit2) {
0531         const auto &detId = hit2->geographicalId();
0532         return (detId.det() == hitDet && detId.subdetId() == hitSubdet && ttopo_->layer(detId) == hitLayer);
0533       });
0534     };
0535 
0536     auto t1HitIter = t1->recHitsBegin();
0537     if (!(*t1HitIter)->isValid()) {
0538       IfLogTrace(debug_, "DuplicateTrackMerger") << " first t1 hit invalid";
0539       return false;
0540     }
0541     auto t2HitIter = findHitOnT2(*t1HitIter);
0542     if (t2HitIter == t2->recHitsEnd()) {
0543       // if first hit not found, try with second
0544       // if that fails, then reject
0545       ++t1HitIter;
0546       assert(t1HitIter != t1->recHitsEnd());
0547 
0548       if (!(*t1HitIter)->isValid()) {
0549         IfLogTrace(debug_, "DuplicateTrackMerger") << " second t1 hit invalid";
0550         return false;
0551       }
0552       t2HitIter = findHitOnT2(*t1HitIter);
0553       if (t2HitIter == t2->recHitsEnd())
0554         return false;
0555     }
0556     IfLogTrace(debug_, "DuplicateTrackMerger")
0557         << " starting overlap check from t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0558         << std::distance(t2->recHitsBegin(), t2HitIter);
0559 
0560     auto prevSubdet = (*t1HitIter)->geographicalId().subdetId();
0561     const TrajectoryStateOnSurface tsosInner = trajectoryStateTransform::innerStateOnSurface(*t2, *geom_, magfield_);
0562 
0563     ++t1HitIter;
0564     ++t2HitIter;
0565     unsigned int missedLayers = 0;
0566     while (t1HitIter != t1->recHitsEnd() && t2HitIter != t2->recHitsEnd()) {
0567       // in case of invalid hits, reject immediately
0568       if ((*t1HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t1HitIter) ||
0569           (*t2HitIter)->getType() != TrackingRecHit::valid || trackerHitRTTI::isUndef(**t2HitIter)) {
0570         IfLogTrace(debug_, "DuplicateTrackMerger")
0571             << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0572             << std::distance(t2->recHitsBegin(), t2HitIter) << " either is invalid, types t1 "
0573             << (*t1HitIter)->getType() << " t2 " << (*t2HitIter)->getType();
0574         return false;
0575       }
0576 
0577       const auto &t1DetId = (*t1HitIter)->geographicalId();
0578       const auto &t2DetId = (*t2HitIter)->geographicalId();
0579 
0580       const auto t1Det = t1DetId.det();
0581       const auto t2Det = t2DetId.det();
0582       if (t1Det != DetId::Tracker || t2Det != DetId::Tracker) {
0583         IfLogTrace(debug_, "DuplicateTrackMerger") << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter)
0584                                                    << " t2 hit " << std::distance(t2->recHitsBegin(), t2HitIter)
0585                                                    << " either not from Tracker, dets t1 " << t1Det << " t2 " << t2Det;
0586         return false;
0587       }
0588 
0589       const auto t1Subdet = t1DetId.subdetId();
0590       const auto t1Layer = ttopo_->layer(t1DetId);
0591 
0592       // reject if hits have the same DetId but are different
0593       if (t1DetId == t2DetId) {
0594         if (!(*t1HitIter)->sharesInput(*t2HitIter, TrackingRecHit::all)) {
0595           IfLogTrace(debug_, "DuplicateTrackMerger")
0596               << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0597               << std::distance(t2->recHitsBegin(), t2HitIter) << " same DetId (" << t1DetId.rawId()
0598               << ") but do not share all input";
0599           return false;
0600         }
0601       } else {
0602         const auto t2Subdet = t2DetId.subdetId();
0603         const auto t2Layer = ttopo_->layer(t2DetId);
0604 
0605         // reject if hits are on different layers
0606         if (t1Subdet != t2Subdet || t1Layer != t2Layer) {
0607           bool recovered = false;
0608           // but try to recover first by checking if either one has skipped over a layer
0609           if (t1Subdet == prevSubdet && t2Subdet != prevSubdet) {
0610             // t1 has a layer t2 doesn't
0611             ++t1HitIter;
0612             recovered = true;
0613           } else if (t1Subdet != prevSubdet && t2Subdet == prevSubdet) {
0614             // t2 has a layer t1 doesn't
0615             ++t2HitIter;
0616             recovered = true;
0617           } else if (t1Subdet == t2Subdet) {
0618             prevSubdet = t1Subdet;
0619             // same subdet, so layer must be different
0620             if (t2Layer > t1Layer) {
0621               // t1 has a layer t2 doesn't
0622               ++t1HitIter;
0623               recovered = true;
0624             } else if (t1Layer > t2Layer) {
0625               // t2 has a layer t1 doesn't
0626               ++t2HitIter;
0627               recovered = true;
0628             }
0629           }
0630           if (recovered) {
0631             ++missedLayers;
0632             if (missedLayers > overlapCheckMaxMissingLayers_) {
0633               IfLogTrace(debug_, "DuplicateTrackMerger") << " max number of missed layers exceeded";
0634               return false;
0635             }
0636             continue;
0637           }
0638 
0639           IfLogTrace(debug_, "DuplicateTrackMerger")
0640               << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0641               << std::distance(t2->recHitsBegin(), t2HitIter) << " are on different layers (subdet, layer) t1 "
0642               << t1Subdet << "," << t1Layer << " t2 " << t2Subdet << "," << t2Layer;
0643           return false;
0644         }
0645         // reject if same layer (but not same hit) in non-pixel detector
0646         else if (t1Subdet != PixelSubdetector::PixelBarrel && t1Subdet != PixelSubdetector::PixelEndcap) {
0647           IfLogTrace(debug_, "DuplicateTrackMerger")
0648               << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0649               << std::distance(t2->recHitsBegin(), t2HitIter) << " are on same layer, but in non-pixel detector (det "
0650               << t1Det << " subdet " << t1Subdet << " layer " << t1Layer << ")";
0651           return false;
0652         }
0653       }
0654 
0655       // Propagate longer track to the shorter track hit surface, check compatibility
0656       TrajectoryStateOnSurface tsosPropagated = propagator_->propagate(tsosInner, (*t1HitIter)->det()->surface());
0657       if (!tsosPropagated.isValid()) {  // reject immediately if TSOS is not valid
0658         IfLogTrace(debug_, "DuplicateTrackMerger")
0659             << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0660             << std::distance(t2->recHitsBegin(), t2HitIter) << " TSOS not valid";
0661         return false;
0662       }
0663       auto passChi2Pair = chi2Estimator_->estimate(tsosPropagated, **t1HitIter);
0664       if (!passChi2Pair.first) {
0665         IfLogTrace(debug_, "DuplicateTrackMerger")
0666             << " t1 hit " << std::distance(t1->recHitsBegin(), t1HitIter) << " t2 hit "
0667             << std::distance(t2->recHitsBegin(), t2HitIter) << " hit chi2 compatibility failed with chi2 "
0668             << passChi2Pair.second;
0669         return false;
0670       }
0671 
0672       prevSubdet = t1Subdet;
0673       ++t1HitIter;
0674       ++t2HitIter;
0675     }
0676     if (t1HitIter != t1->recHitsEnd()) {
0677       IfLogTrace(debug_, "DuplicateTrackMerger") << " hits on t2 ended before hits on t1";
0678       return false;
0679     }
0680 
0681     IfLogTrace(debug_, "DuplicateTrackMerger") << " all hits on t2 are on layers whre t1 has also a hit";
0682     return true;
0683   }
0684 }  // namespace
0685 
0686 #include "FWCore/PluginManager/interface/ModuleDef.h"
0687 #include "FWCore/Framework/interface/MakerMacros.h"
0688 
0689 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!