Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-13 01:32:36 

0001 // -*- C++ -*-
0002 //
0003 // Package:    NtupleDump/TrackingNtuple
0004 // Class:      TrackingNtuple
0005 //
0006 /**\class TrackingNtuple TrackingNtuple.cc NtupleDump/TrackingNtuple/plugins/TrackingNtuple.cc
0007 
0008    Description: [one line class summary]
0009 
0010    Implementation:
0011    [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Giuseppe Cerati
0015 //         Created:  Tue, 25 Aug 2015 13:22:49 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025 
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/Utilities/interface/isFinite.h"
0031 #include "FWCore/ServiceRegistry/interface/Service.h"
0032 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0033 #include "CommonTools/Utils/interface/DynArray.h"
0034 #include "DataFormats/Provenance/interface/ProductID.h"
0035 #include "DataFormats/Common/interface/ContainerMask.h"
0036 
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/Utilities/interface/transform.h"
0039 
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0042 #include "DataFormats/TrackReco/interface/SeedStopInfo.h"
0043 
0044 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0045 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0046 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0047 
0048 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0049 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0050 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0051 #include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"
0052 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0053 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
0054 
0055 #include "MagneticField/Engine/interface/MagneticField.h"
0056 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0057 
0058 #include "DataFormats/SiPixelDetId/interface/PixelChannelIdentifier.h"
0059 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0060 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0061 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0062 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0063 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0064 #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h"
0065 
0066 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0067 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0068 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0069 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0070 
0071 #include "DataFormats/VertexReco/interface/Vertex.h"
0072 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0073 
0074 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0075 #include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
0076 #include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"
0077 #include "SimTracker/TrackAssociation/interface/ParametersDefinerForTP.h"
0078 #include "SimTracker/TrackAssociation/interface/TrackingParticleIP.h"
0079 #include "SimTracker/TrackAssociation/interface/trackAssociationChi2.h"
0080 #include "SimTracker/TrackAssociation/interface/trackHitsToClusterRefs.h"
0081 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0082 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0083 
0084 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0085 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0086 
0087 #include "SimTracker/TrackHistory/interface/HistoryBase.h"
0088 #include "HepPDT/ParticleID.hh"
0089 
0090 #include "Validation/RecoTrack/interface/trackFromSeedFitFailed.h"
0091 
0092 #include "RecoTracker/FinalTrackSelectors/interface/getBestVertex.h"
0093 
0094 #include <set>
0095 #include <map>
0096 #include <unordered_set>
0097 #include <unordered_map>
0098 #include <tuple>
0099 #include <utility>
0100 
0101 #include "TTree.h"
0102 
0103 /*
0104 todo: 
0105 add refitted hit position after track/seed fit
0106 add local angle, path length!
0107 */
0108 
0109 namespace {
0110   // This pattern is copied from QuickTrackAssociatorByHitsImpl. If
0111   // further needs arises, it wouldn't hurt to abstract it somehow.
0112   using TrackingParticleRefKeyToIndex = std::unordered_map<reco::RecoToSimCollection::index_type, size_t>;
0113   using TrackingVertexRefKeyToIndex = TrackingParticleRefKeyToIndex;
0114   using SimHitFullKey = std::pair<TrackPSimHitRef::key_type, edm::ProductID>;
0115   using SimHitRefKeyToIndex = std::map<SimHitFullKey, size_t>;
0116   using TrackingParticleRefKeyToCount = TrackingParticleRefKeyToIndex;
0117 
0118   std::string subdetstring(int subdet) {
0119     switch (subdet) {
0120       case StripSubdetector::TIB:
0121         return "- TIB";
0122       case StripSubdetector::TOB:
0123         return "- TOB";
0124       case StripSubdetector::TEC:
0125         return "- TEC";
0126       case StripSubdetector::TID:
0127         return "- TID";
0128       case PixelSubdetector::PixelBarrel:
0129         return "- PixBar";
0130       case PixelSubdetector::PixelEndcap:
0131         return "- PixFwd";
0132       default:
0133         return "UNKNOWN TRACKER HIT TYPE";
0134     }
0135   }
0136 
0137   struct ProductIDSetPrinter {
0138     ProductIDSetPrinter(const std::set<edm::ProductID>& set) : set_(set) {}
0139 
0140     void print(std::ostream& os) const {
0141       for (const auto& item : set_) {
0142         os << item << " ";
0143       }
0144     }
0145 
0146     const std::set<edm::ProductID>& set_;
0147   };
0148   std::ostream& operator<<(std::ostream& os, const ProductIDSetPrinter& o) {
0149     o.print(os);
0150     return os;
0151   }
0152   template <typename T>
0153   struct ProductIDMapPrinter {
0154     ProductIDMapPrinter(const std::map<edm::ProductID, T>& map) : map_(map) {}
0155 
0156     void print(std::ostream& os) const {
0157       for (const auto& item : map_) {
0158         os << item.first << " ";
0159       }
0160     }
0161 
0162     const std::map<edm::ProductID, T>& map_;
0163   };
0164   template <typename T>
0165   auto make_ProductIDMapPrinter(const std::map<edm::ProductID, T>& map) {
0166     return ProductIDMapPrinter<T>(map);
0167   }
0168   template <typename T>
0169   std::ostream& operator<<(std::ostream& os, const ProductIDMapPrinter<T>& o) {
0170     o.print(os);
0171     return os;
0172   }
0173 
0174   template <typename T>
0175   struct VectorPrinter {
0176     VectorPrinter(const std::vector<T>& vec) : vec_(vec) {}
0177 
0178     void print(std::ostream& os) const {
0179       for (const auto& item : vec_) {
0180         os << item << " ";
0181       }
0182     }
0183 
0184     const std::vector<T>& vec_;
0185   };
0186   template <typename T>
0187   auto make_VectorPrinter(const std::vector<T>& vec) {
0188     return VectorPrinter<T>(vec);
0189   }
0190   template <typename T>
0191   std::ostream& operator<<(std::ostream& os, const VectorPrinter<T>& o) {
0192     o.print(os);
0193     return os;
0194   }
0195 
0196   void checkProductID(const std::set<edm::ProductID>& set, const edm::ProductID& id, const char* name) {
0197     if (set.find(id) == set.end())
0198       throw cms::Exception("Configuration")
0199           << "Got " << name << " with a hit with ProductID " << id
0200           << " which does not match to the set of ProductID's for the hits: " << ProductIDSetPrinter(set)
0201           << ". Usually this is caused by a wrong hit collection in the configuration.";
0202   }
0203 
0204   template <typename SimLink, typename Func>
0205   void forEachMatchedSimLink(const edm::DetSet<SimLink>& digiSimLinks, uint32_t channel, Func func) {
0206     for (const auto& link : digiSimLinks) {
0207       if (link.channel() == channel) {
0208         func(link);
0209       }
0210     }
0211   }
0212 
0213   /// No-op function used in the trick of CombineDetId::impl2()
0214   template <typename... Args>
0215   void call_nop(Args&&... args) {}
0216 
0217   template <typename... Types>
0218   class CombineDetId {
0219   public:
0220     CombineDetId() {}
0221 
0222     /// Return the raw DetId, assumes that the first type is
0223     /// DetIdCommon that has operator[]()
0224     unsigned int operator[](size_t i) const { return std::get<0>(content_)[i]; }
0225 
0226     template <typename... Args>
0227     void book(Args&&... args) {
0228       impl([&](auto& vec) { vec.book(std::forward<Args>(args)...); });
0229     }
0230 
0231     template <typename... Args>
0232     void push_back(Args&&... args) {
0233       impl([&](auto& vec) { vec.push_back(std::forward<Args>(args)...); });
0234     }
0235 
0236     template <typename... Args>
0237     void resize(Args&&... args) {
0238       impl([&](auto& vec) { vec.resize(std::forward<Args>(args)...); });
0239     }
0240 
0241     template <typename... Args>
0242     void set(Args&&... args) {
0243       impl([&](auto& vec) { vec.set(std::forward<Args>(args)...); });
0244     }
0245 
0246     void clear() {
0247       impl([&](auto& vec) { vec.clear(); });
0248     }
0249 
0250   private:
0251     // Trick to not repeate std::index_sequence_for in each of the methods above
0252     template <typename F>
0253     void impl(F&& func) {
0254       impl2(std::index_sequence_for<Types...>{}, std::forward<F>(func));
0255     }
0256 
0257     // Trick to exploit parameter pack expansion in function call
0258     // arguments to call a member function for each tuple element
0259     // (with the same signature). The comma operator is needed to
0260     // return a value from the expression as an argument for the
0261     // call_nop.
0262     template <std::size_t... Is, typename F>
0263     void impl2(std::index_sequence<Is...>, F&& func) {
0264       call_nop((func(std::get<Is>(content_)), 0)...);
0265     }
0266 
0267     std::tuple<Types...> content_;
0268   };
0269 
0270   std::map<unsigned int, double> chargeFraction(const SiPixelCluster& cluster,
0271                                                 const DetId& detId,
0272                                                 const edm::DetSetVector<PixelDigiSimLink>& digiSimLink) {
0273     std::map<unsigned int, double> simTrackIdToAdc;
0274 
0275     auto idetset = digiSimLink.find(detId);
0276     if (idetset == digiSimLink.end())
0277       return simTrackIdToAdc;
0278 
0279     double adcSum = 0;
0280     PixelDigiSimLink found;
0281     for (int iPix = 0; iPix != cluster.size(); ++iPix) {
0282       const SiPixelCluster::Pixel& pixel = cluster.pixel(iPix);
0283       adcSum += pixel.adc;
0284       uint32_t channel = PixelChannelIdentifier::pixelToChannel(pixel.x, pixel.y);
0285       forEachMatchedSimLink(*idetset, channel, [&](const PixelDigiSimLink& simLink) {
0286         double& adc = simTrackIdToAdc[simLink.SimTrackId()];
0287         adc += pixel.adc * simLink.fraction();
0288       });
0289     }
0290 
0291     for (auto& pair : simTrackIdToAdc) {
0292       if (adcSum == 0.)
0293         pair.second = 0.;
0294       else
0295         pair.second /= adcSum;
0296     }
0297 
0298     return simTrackIdToAdc;
0299   }
0300 
0301   std::map<unsigned int, double> chargeFraction(const SiStripCluster& cluster,
0302                                                 const DetId& detId,
0303                                                 const edm::DetSetVector<StripDigiSimLink>& digiSimLink) {
0304     std::map<unsigned int, double> simTrackIdToAdc;
0305 
0306     auto idetset = digiSimLink.find(detId);
0307     if (idetset == digiSimLink.end())
0308       return simTrackIdToAdc;
0309 
0310     double adcSum = 0;
0311     StripDigiSimLink found;
0312     int first = cluster.firstStrip();
0313     for (size_t i = 0; i < cluster.amplitudes().size(); ++i) {
0314       adcSum += cluster.amplitudes()[i];
0315       forEachMatchedSimLink(*idetset, first + i, [&](const StripDigiSimLink& simLink) {
0316         double& adc = simTrackIdToAdc[simLink.SimTrackId()];
0317         adc += cluster.amplitudes()[i] * simLink.fraction();
0318       });
0319 
0320       for (const auto& pair : simTrackIdToAdc) {
0321         simTrackIdToAdc[pair.first] = (adcSum != 0. ? pair.second / adcSum : 0.);
0322       }
0323     }
0324     return simTrackIdToAdc;
0325   }
0326 
0327   std::map<unsigned int, double> chargeFraction(const Phase2TrackerCluster1D& cluster,
0328                                                 const DetId& detId,
0329                                                 const edm::DetSetVector<StripDigiSimLink>& digiSimLink) {
0330     std::map<unsigned int, double> simTrackIdToAdc;
0331     throw cms::Exception("LogicError") << "Not possible to use StripDigiSimLink with Phase2TrackerCluster1D! ";
0332     return simTrackIdToAdc;
0333   }
0334 
0335   //In the OT, there is no measurement of the charge, so no ADC value.
0336   //Only in the SSA chip (so in PSs) you have one "threshold" flag that tells you if the charge of at least one strip in the cluster exceeded 1.2 MIPs.
0337   std::map<unsigned int, double> chargeFraction(const Phase2TrackerCluster1D& cluster,
0338                                                 const DetId& detId,
0339                                                 const edm::DetSetVector<PixelDigiSimLink>& digiSimLink) {
0340     std::map<unsigned int, double> simTrackIdToAdc;
0341     return simTrackIdToAdc;
0342   }
0343 
0344   struct TrackTPMatch {
0345     int key = -1;
0346     int countClusters = 0;
0347   };
0348 
0349   template <typename HRange>
0350   TrackTPMatch findBestMatchingTrackingParticle(const HRange& hits,
0351                                                 const ClusterTPAssociation& clusterToTPMap,
0352                                                 const TrackingParticleRefKeyToIndex& tpKeyToIndex) {
0353     struct Count {
0354       int clusters = 0;
0355       size_t innermostHit = std::numeric_limits<size_t>::max();
0356     };
0357 
0358     std::vector<OmniClusterRef> clusters = track_associator::hitsToClusterRefs(hits.begin(), hits.end());
0359 
0360     std::unordered_map<int, Count> count;
0361     for (size_t iCluster = 0, end = clusters.size(); iCluster < end; ++iCluster) {
0362       const auto& clusterRef = clusters[iCluster];
0363 
0364       auto range = clusterToTPMap.equal_range(clusterRef);
0365       for (auto ip = range.first; ip != range.second; ++ip) {
0366         const auto tpKey = ip->second.key();
0367         if (tpKeyToIndex.find(tpKey) == tpKeyToIndex.end())  // filter out TPs not given as an input
0368           continue;
0369 
0370         auto& elem = count[tpKey];
0371         ++elem.clusters;
0372         elem.innermostHit = std::min(elem.innermostHit, iCluster);
0373       }
0374     }
0375 
0376     // In case there are many matches with the same number of clusters,
0377     // select the one with innermost hit
0378     TrackTPMatch best;
0379     int bestCount = 2;  // require >= 3 cluster for the best match
0380     size_t bestInnermostHit = std::numeric_limits<size_t>::max();
0381     for (auto& keyCount : count) {
0382       if (keyCount.second.clusters > bestCount ||
0383           (keyCount.second.clusters == bestCount && keyCount.second.innermostHit < bestInnermostHit)) {
0384         best.key = keyCount.first;
0385         best.countClusters = bestCount = keyCount.second.clusters;
0386         bestInnermostHit = keyCount.second.innermostHit;
0387       }
0388     }
0389 
0390     LogTrace("TrackingNtuple") << "findBestMatchingTrackingParticle key " << best.key;
0391 
0392     return best;
0393   }
0394 
0395   template <typename HRange>
0396   TrackTPMatch findMatchingTrackingParticleFromFirstHit(const HRange& hits,
0397                                                         const ClusterTPAssociation& clusterToTPMap,
0398                                                         const TrackingParticleRefKeyToIndex& tpKeyToIndex) {
0399     TrackTPMatch best;
0400 
0401     std::vector<OmniClusterRef> clusters = track_associator::hitsToClusterRefs(hits.begin(), hits.end());
0402     if (clusters.empty()) {
0403       return best;
0404     }
0405 
0406     auto operateCluster = [&](const auto& clusterRef, const auto& func) {
0407       auto range = clusterToTPMap.equal_range(clusterRef);
0408       for (auto ip = range.first; ip != range.second; ++ip) {
0409         const auto tpKey = ip->second.key();
0410         if (tpKeyToIndex.find(tpKey) == tpKeyToIndex.end())  // filter out TPs not given as an input
0411           continue;
0412         func(tpKey);
0413       }
0414     };
0415 
0416     std::vector<unsigned int>
0417         validTPs;  // first cluster can be associated to multiple TPs, use vector as set as this should be small
0418     auto iCluster = clusters.begin();
0419     operateCluster(*iCluster, [&](unsigned int tpKey) { validTPs.push_back(tpKey); });
0420     if (validTPs.empty()) {
0421       return best;
0422     }
0423     ++iCluster;
0424     ++best.countClusters;
0425 
0426     std::vector<bool> foundTPs(validTPs.size(), false);
0427     for (auto iEnd = clusters.end(); iCluster != iEnd; ++iCluster) {
0428       const auto& clusterRef = *iCluster;
0429 
0430       // find out to which first-cluster TPs this cluster is matched to
0431       operateCluster(clusterRef, [&](unsigned int tpKey) {
0432         auto found = std::find(cbegin(validTPs), cend(validTPs), tpKey);
0433         if (found != cend(validTPs)) {
0434           foundTPs[std::distance(cbegin(validTPs), found)] = true;
0435         }
0436       });
0437 
0438       // remove the non-found TPs
0439       auto iTP = validTPs.size();
0440       do {
0441         --iTP;
0442 
0443         if (!foundTPs[iTP]) {
0444           validTPs.erase(validTPs.begin() + iTP);
0445           foundTPs.erase(foundTPs.begin() + iTP);
0446         }
0447       } while (iTP > 0);
0448       if (!validTPs.empty()) {
0449         // for multiple TPs the "first one" is a bit arbitrary, but
0450         // I hope it is rare that a track would have many
0451         // consecutive hits matched to two TPs
0452         best.key = validTPs[0];
0453       } else {
0454         break;
0455       }
0456 
0457       std::fill(begin(foundTPs), end(foundTPs), false);
0458       ++best.countClusters;
0459     }
0460 
0461     // Reqquire >= 3 clusters for a match
0462     return best.countClusters >= 3 ? best : TrackTPMatch();
0463   }
0464 }  // namespace
0465 
0466 //
0467 // class declaration
0468 //
0469 
0470 class TrackingNtuple : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0471 public:
0472   explicit TrackingNtuple(const edm::ParameterSet&);
0473   ~TrackingNtuple() override;
0474 
0475   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0476 
0477 private:
0478   void analyze(const edm::Event&, const edm::EventSetup&) override;
0479 
0480   void clearVariables();
0481 
0482   enum class HitType { Pixel = 0, Strip = 1, Glued = 2, Invalid = 3, Phase2OT = 4, Unknown = 99 };
0483 
0484   // This gives the "best" classification of a reco hit
0485   // In case of reco hit mathing to multiple sim, smaller number is
0486   // considered better
0487   // To be kept in synch with class HitSimType in ntuple.py
0488   enum class HitSimType { Signal = 0, ITPileup = 1, OOTPileup = 2, Noise = 3, Unknown = 99 };
0489 
0490   using MVACollection = std::vector<float>;
0491   using QualityMaskCollection = std::vector<unsigned char>;
0492 
0493   using PixelMaskContainer = edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster>>;
0494   using StripMaskContainer = edm::ContainerMask<edmNew::DetSetVector<SiStripCluster>>;
0495 
0496   struct TPHitIndex {
0497     TPHitIndex(unsigned int tp = 0, unsigned int simHit = 0, float to = 0, unsigned int id = 0)
0498         : tpKey(tp), simHitIdx(simHit), tof(to), detId(id) {}
0499     unsigned int tpKey;
0500     unsigned int simHitIdx;
0501     float tof;
0502     unsigned int detId;
0503   };
0504   static bool tpHitIndexListLess(const TPHitIndex& i, const TPHitIndex& j) { return (i.tpKey < j.tpKey); }
0505   static bool tpHitIndexListLessSort(const TPHitIndex& i, const TPHitIndex& j) {
0506     if (i.tpKey == j.tpKey) {
0507       if (edm::isNotFinite(i.tof) && edm::isNotFinite(j.tof)) {
0508         return i.detId < j.detId;
0509       }
0510       return i.tof < j.tof;  // works as intended if either one is NaN
0511     }
0512     return i.tpKey < j.tpKey;
0513   }
0514 
0515   void fillBeamSpot(const reco::BeamSpot& bs);
0516   void fillPixelHits(const edm::Event& iEvent,
0517                      const TrackerGeometry& tracker,
0518                      const ClusterTPAssociation& clusterToTPMap,
0519                      const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0520                      const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
0521                      const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
0522                      const TrackerTopology& tTopo,
0523                      const SimHitRefKeyToIndex& simHitRefKeyToIndex,
0524                      std::set<edm::ProductID>& hitProductIds);
0525 
0526   void fillStripRphiStereoHits(const edm::Event& iEvent,
0527                                const TrackerGeometry& tracker,
0528                                const ClusterTPAssociation& clusterToTPMap,
0529                                const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0530                                const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
0531                                const edm::DetSetVector<StripDigiSimLink>& digiSimLink,
0532                                const TrackerTopology& tTopo,
0533                                const SimHitRefKeyToIndex& simHitRefKeyToIndex,
0534                                std::set<edm::ProductID>& hitProductIds);
0535 
0536   void fillStripMatchedHits(const edm::Event& iEvent,
0537                             const TrackerGeometry& tracker,
0538                             const TrackerTopology& tTopo,
0539                             std::vector<std::pair<int, int>>& monoStereoClusterList);
0540 
0541   size_t addStripMatchedHit(const SiStripMatchedRecHit2D& hit,
0542                             const TrackerGeometry& tracker,
0543                             const TrackerTopology& tTopo,
0544                             const std::vector<std::pair<uint64_t, StripMaskContainer const*>>& stripMasks,
0545                             std::vector<std::pair<int, int>>& monoStereoClusterList);
0546 
0547   void fillPhase2OTHits(const edm::Event& iEvent,
0548                         const ClusterTPAssociation& clusterToTPMap,
0549                         const TrackerGeometry& tracker,
0550                         const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0551                         const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
0552                         const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
0553                         const TrackerTopology& tTopo,
0554                         const SimHitRefKeyToIndex& simHitRefKeyToIndex,
0555                         std::set<edm::ProductID>& hitProductIds);
0556 
0557   void fillSeeds(const edm::Event& iEvent,
0558                  const TrackingParticleRefVector& tpCollection,
0559                  const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0560                  const reco::BeamSpot& bs,
0561                  const TrackerGeometry& tracker,
0562                  const reco::TrackToTrackingParticleAssociator& associatorByHits,
0563                  const ClusterTPAssociation& clusterToTPMap,
0564                  const MagneticField& theMF,
0565                  const TrackerTopology& tTopo,
0566                  std::vector<std::pair<int, int>>& monoStereoClusterList,
0567                  const std::set<edm::ProductID>& hitProductIds,
0568                  std::map<edm::ProductID, size_t>& seedToCollIndex);
0569 
0570   void fillTracks(const edm::RefToBaseVector<reco::Track>& tracks,
0571                   const TrackerGeometry& tracker,
0572                   const TrackingParticleRefVector& tpCollection,
0573                   const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0574                   const TrackingParticleRefKeyToCount& tpKeyToClusterCount,
0575                   const MagneticField& mf,
0576                   const reco::BeamSpot& bs,
0577                   const reco::VertexCollection& vertices,
0578                   const reco::TrackToTrackingParticleAssociator& associatorByHits,
0579                   const ClusterTPAssociation& clusterToTPMap,
0580                   const TrackerTopology& tTopo,
0581                   const std::set<edm::ProductID>& hitProductIds,
0582                   const std::map<edm::ProductID, size_t>& seedToCollIndex,
0583                   const std::vector<const MVACollection*>& mvaColls,
0584                   const std::vector<const QualityMaskCollection*>& qualColls);
0585 
0586   void fillCandidates(const edm::Handle<TrackCandidateCollection>& candsHandle,
0587                       int algo,
0588                       const TrackingParticleRefVector& tpCollection,
0589                       const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0590                       const TrackingParticleRefKeyToCount& tpKeyToClusterCount,
0591                       const MagneticField& mf,
0592                       const reco::BeamSpot& bs,
0593                       const reco::VertexCollection& vertices,
0594                       const reco::TrackToTrackingParticleAssociator& associatorByHits,
0595                       const ClusterTPAssociation& clusterToTPMap,
0596                       const TrackerGeometry& tracker,
0597                       const TrackerTopology& tTopo,
0598                       const std::set<edm::ProductID>& hitProductIds,
0599                       const std::map<edm::ProductID, size_t>& seedToCollIndex);
0600 
0601   void fillSimHits(const TrackerGeometry& tracker,
0602                    const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0603                    const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
0604                    const TrackerTopology& tTopo,
0605                    SimHitRefKeyToIndex& simHitRefKeyToIndex,
0606                    std::vector<TPHitIndex>& tpHitList);
0607 
0608   void fillTrackingParticles(const edm::Event& iEvent,
0609                              const edm::EventSetup& iSetup,
0610                              const edm::RefToBaseVector<reco::Track>& tracks,
0611                              const reco::BeamSpot& bs,
0612                              const TrackingParticleRefVector& tpCollection,
0613                              const TrackingVertexRefKeyToIndex& tvKeyToIndex,
0614                              const reco::TrackToTrackingParticleAssociator& associatorByHits,
0615                              const std::vector<TPHitIndex>& tpHitList,
0616                              const TrackingParticleRefKeyToCount& tpKeyToClusterCount);
0617 
0618   void fillTrackingParticlesForSeeds(const TrackingParticleRefVector& tpCollection,
0619                                      const reco::SimToRecoCollection& simRecColl,
0620                                      const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0621                                      const unsigned int seedOffset);
0622 
0623   void fillVertices(const reco::VertexCollection& vertices, const edm::RefToBaseVector<reco::Track>& tracks);
0624 
0625   void fillTrackingVertices(const TrackingVertexRefVector& trackingVertices,
0626                             const TrackingParticleRefKeyToIndex& tpKeyToIndex);
0627 
0628   struct SimHitData {
0629     std::vector<int> matchingSimHit;
0630     std::vector<float> chargeFraction;
0631     std::vector<float> xySignificance;
0632     std::vector<int> bunchCrossing;
0633     std::vector<int> event;
0634     HitSimType type = HitSimType::Unknown;
0635   };
0636 
0637   template <typename SimLink>
0638   SimHitData matchCluster(const OmniClusterRef& cluster,
0639                           DetId hitId,
0640                           int clusterKey,
0641                           const TrackingRecHit& hit,
0642                           const ClusterTPAssociation& clusterToTPMap,
0643                           const TrackingParticleRefKeyToIndex& tpKeyToIndex,
0644                           const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
0645                           const edm::DetSetVector<SimLink>& digiSimLinks,
0646                           const SimHitRefKeyToIndex& simHitRefKeyToIndex,
0647                           HitType hitType);
0648 
0649   // ----------member data ---------------------------
0650   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken_;
0651   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0652   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tGeomToken_;
0653 
0654   std::vector<edm::EDGetTokenT<edm::View<reco::Track>>> seedTokens_;
0655   std::vector<edm::EDGetTokenT<std::vector<SeedStopInfo>>> seedStopInfoTokens_;
0656 
0657   std::vector<edm::EDGetTokenT<TrackCandidateCollection>> candidateTokens_;
0658   edm::EDGetTokenT<edm::View<reco::Track>> trackToken_;
0659   std::vector<std::tuple<edm::EDGetTokenT<MVACollection>, edm::EDGetTokenT<QualityMaskCollection>>>
0660       mvaQualityCollectionTokens_;
0661   edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
0662   edm::EDGetTokenT<TrackingParticleRefVector> trackingParticleRefToken_;
0663   edm::EDGetTokenT<ClusterTPAssociation> clusterTPMapToken_;
0664   edm::EDGetTokenT<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitTPMapToken_;
0665   edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> trackAssociatorToken_;
0666   edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink>> pixelSimLinkToken_;
0667   edm::EDGetTokenT<edm::DetSetVector<StripDigiSimLink>> stripSimLinkToken_;
0668   edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink>> siphase2OTSimLinksToken_;
0669   bool includeStripHits_, includePhase2OTHits_;
0670   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0671   edm::EDGetTokenT<SiPixelRecHitCollection> pixelRecHitToken_;
0672   edm::EDGetTokenT<SiStripRecHit2DCollection> stripRphiRecHitToken_;
0673   edm::EDGetTokenT<SiStripRecHit2DCollection> stripStereoRecHitToken_;
0674   edm::EDGetTokenT<SiStripMatchedRecHit2DCollection> stripMatchedRecHitToken_;
0675   edm::EDGetTokenT<Phase2TrackerRecHit1DCollectionNew> phase2OTRecHitToken_;
0676   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0677   edm::EDGetTokenT<TrackingVertexCollection> trackingVertexToken_;
0678   edm::EDGetTokenT<edm::ValueMap<unsigned int>> tpNLayersToken_;
0679   edm::EDGetTokenT<edm::ValueMap<unsigned int>> tpNPixelLayersToken_;
0680   edm::EDGetTokenT<edm::ValueMap<unsigned int>> tpNStripStereoLayersToken_;
0681 
0682   std::vector<std::pair<unsigned int, edm::EDGetTokenT<PixelMaskContainer>>> pixelUseMaskTokens_;
0683   std::vector<std::pair<unsigned int, edm::EDGetTokenT<StripMaskContainer>>> stripUseMaskTokens_;
0684 
0685   std::string builderName_;
0686   const bool includeSeeds_;
0687   const bool includeTrackCandidates_;
0688   const bool addSeedCurvCov_;
0689   const bool includeAllHits_;
0690   const bool includeOnTrackHitData_;
0691   const bool includeMVA_;
0692   const bool includeTrackingParticles_;
0693   const bool includeOOT_;
0694   const bool keepEleSimHits_;
0695   const bool saveSimHitsP3_;
0696   const bool simHitBySignificance_;
0697 
0698   HistoryBase tracer_;
0699   ParametersDefinerForTP parametersDefiner_;
0700 
0701   TTree* t;
0702 
0703   // DetId branches
0704 #define BOOK(name) tree->Branch((prefix + "_" + #name).c_str(), &name);
0705   class DetIdCommon {
0706   public:
0707     DetIdCommon(){};
0708 
0709     unsigned int operator[](size_t i) const { return detId[i]; }
0710 
0711     void book(const std::string& prefix, TTree* tree) {
0712       BOOK(detId);
0713       BOOK(subdet);
0714       BOOK(layer);
0715       BOOK(side);
0716       BOOK(module);
0717       BOOK(moduleType);
0718     }
0719 
0720     void push_back(const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0721       detId.push_back(id.rawId());
0722       subdet.push_back(id.subdetId());
0723       layer.push_back(tTopo.layer(id));
0724       module.push_back(tTopo.module(id));
0725       moduleType.push_back(static_cast<int>(tracker.getDetectorType(id)));
0726 
0727       unsigned short s = 0;
0728       switch (id.subdetId()) {
0729         case StripSubdetector::TIB:
0730           s = tTopo.tibSide(id);
0731           break;
0732         case StripSubdetector::TOB:
0733           s = tTopo.tobSide(id);
0734           break;
0735         default:
0736           s = tTopo.side(id);
0737       }
0738       side.push_back(s);
0739     }
0740 
0741     void resize(size_t size) {
0742       detId.resize(size);
0743       subdet.resize(size);
0744       layer.resize(size);
0745       side.resize(size);
0746       module.resize(size);
0747       moduleType.resize(size);
0748     }
0749 
0750     void set(size_t index, const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0751       detId[index] = id.rawId();
0752       subdet[index] = id.subdetId();
0753       layer[index] = tTopo.layer(id);
0754       side[index] = tTopo.side(id);
0755       module[index] = tTopo.module(id);
0756       moduleType[index] = static_cast<int>(tracker.getDetectorType(id));
0757     }
0758 
0759     void clear() {
0760       detId.clear();
0761       subdet.clear();
0762       layer.clear();
0763       side.clear();
0764       module.clear();
0765       moduleType.clear();
0766     }
0767 
0768   private:
0769     std::vector<unsigned int> detId;
0770     std::vector<unsigned short> subdet;
0771     std::vector<unsigned short> layer;  // or disk/wheel
0772     std::vector<unsigned short> side;
0773     std::vector<unsigned short> module;
0774     std::vector<unsigned int> moduleType;
0775   };
0776 
0777   class DetIdPixelOnly {
0778   public:
0779     DetIdPixelOnly() {}
0780 
0781     void book(const std::string& prefix, TTree* tree) {
0782       BOOK(ladder);
0783       BOOK(blade);
0784       BOOK(panel);
0785     }
0786 
0787     void push_back(const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0788       const bool isBarrel = id.subdetId() == PixelSubdetector::PixelBarrel;
0789       ladder.push_back(isBarrel ? tTopo.pxbLadder(id) : 0);
0790       blade.push_back(isBarrel ? 0 : tTopo.pxfBlade(id));
0791       panel.push_back(isBarrel ? 0 : tTopo.pxfPanel(id));
0792     }
0793 
0794     void clear() {
0795       ladder.clear();
0796       blade.clear();
0797       panel.clear();
0798     }
0799 
0800   private:
0801     std::vector<unsigned short> ladder;
0802     std::vector<unsigned short> blade;
0803     std::vector<unsigned short> panel;
0804   };
0805 
0806   class DetIdOTCommon {
0807   public:
0808     DetIdOTCommon() {}
0809 
0810     void book(const std::string& prefix, TTree* tree) {
0811       BOOK(order);
0812       BOOK(ring);
0813       BOOK(rod);
0814     }
0815 
0816     void push_back(const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0817       const auto parsed = parse(tTopo, id);
0818       order.push_back(parsed.order);
0819       ring.push_back(parsed.ring);
0820       rod.push_back(parsed.rod);
0821     }
0822 
0823     void resize(size_t size) {
0824       order.resize(size);
0825       ring.resize(size);
0826       rod.resize(size);
0827     }
0828 
0829     void set(size_t index, const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0830       const auto parsed = parse(tTopo, id);
0831       order[index] = parsed.order;
0832       ring[index] = parsed.ring;
0833       rod[index] = parsed.rod;
0834     }
0835 
0836     void clear() {
0837       order.clear();
0838       ring.clear();
0839       rod.clear();
0840     }
0841 
0842   private:
0843     struct Parsed {
0844       // use int here instead of short to avoid compilation errors due
0845       // to narrowing conversion (less boilerplate than explicit static_casts)
0846       unsigned int order = 0;
0847       unsigned int ring = 0;
0848       unsigned int rod = 0;
0849     };
0850     Parsed parse(const TrackerTopology& tTopo, const DetId& id) const {
0851       switch (id.subdetId()) {
0852         case StripSubdetector::TIB:
0853           return Parsed{tTopo.tibOrder(id), 0, 0};
0854         case StripSubdetector::TID:
0855           return Parsed{tTopo.tidOrder(id), tTopo.tidRing(id), 0};
0856         case StripSubdetector::TOB:
0857           return Parsed{0, 0, tTopo.tobRod(id)};
0858         case StripSubdetector::TEC:
0859           return Parsed{tTopo.tecOrder(id), tTopo.tecRing(id), 0};
0860         default:
0861           return Parsed{};
0862       };
0863     }
0864 
0865     std::vector<unsigned short> order;
0866     std::vector<unsigned short> ring;
0867     std::vector<unsigned short> rod;
0868   };
0869 
0870   class DetIdStripOnly {
0871   public:
0872     DetIdStripOnly() {}
0873 
0874     void book(const std::string& prefix, TTree* tree) {
0875       BOOK(string);
0876       BOOK(petalNumber);
0877       BOOK(isStereo);
0878       BOOK(isRPhi);
0879       BOOK(isGlued);
0880     }
0881 
0882     void push_back(const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0883       const auto parsed = parse(tTopo, id);
0884       string.push_back(parsed.string);
0885       petalNumber.push_back(parsed.petalNumber);
0886       isStereo.push_back(tTopo.isStereo(id));
0887       isRPhi.push_back(tTopo.isRPhi(id));
0888       isGlued.push_back(parsed.glued);
0889     }
0890 
0891     void resize(size_t size) {
0892       string.resize(size);
0893       petalNumber.resize(size);
0894       isStereo.resize(size);
0895       isRPhi.resize(size);
0896       isGlued.resize(size);
0897     }
0898 
0899     void set(size_t index, const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0900       const auto parsed = parse(tTopo, id);
0901       string[index] = parsed.string;
0902       petalNumber[index] = parsed.petalNumber;
0903       isStereo[index] = tTopo.isStereo(id);
0904       isRPhi[index] = tTopo.isRPhi(id);
0905       isGlued[index] = parsed.glued;
0906     }
0907 
0908     void clear() {
0909       string.clear();
0910       isStereo.clear();
0911       isRPhi.clear();
0912       isGlued.clear();
0913       petalNumber.clear();
0914     }
0915 
0916   private:
0917     struct Parsed {
0918       // use int here instead of short to avoid compilation errors due
0919       // to narrowing conversion (less boilerplate than explicit static_casts)
0920       unsigned int string = 0;
0921       unsigned int petalNumber = 0;
0922       bool glued = false;
0923     };
0924     Parsed parse(const TrackerTopology& tTopo, const DetId& id) const {
0925       switch (id.subdetId()) {
0926         case StripSubdetector::TIB:
0927           return Parsed{tTopo.tibString(id), 0, tTopo.tibIsDoubleSide(id)};
0928         case StripSubdetector::TID:
0929           return Parsed{0, 0, tTopo.tidIsDoubleSide(id)};
0930         case StripSubdetector::TOB:
0931           return Parsed{0, 0, tTopo.tobIsDoubleSide(id)};
0932         case StripSubdetector::TEC:
0933           return Parsed{0, tTopo.tecPetalNumber(id), tTopo.tecIsDoubleSide(id)};
0934         default:
0935           return Parsed{};
0936       }
0937     }
0938 
0939     std::vector<unsigned short> string;
0940     std::vector<unsigned short> petalNumber;
0941     std::vector<unsigned short> isStereo;
0942     std::vector<unsigned short> isRPhi;
0943     std::vector<unsigned short> isGlued;
0944   };
0945 
0946   class DetIdPhase2OTOnly {
0947   public:
0948     DetIdPhase2OTOnly() {}
0949 
0950     void book(const std::string& prefix, TTree* tree) {
0951       BOOK(isLower);
0952       BOOK(isUpper);
0953       BOOK(isStack);
0954     }
0955 
0956     void push_back(const TrackerGeometry& tracker, const TrackerTopology& tTopo, const DetId& id) {
0957       isLower.push_back(tTopo.isLower(id));
0958       isUpper.push_back(tTopo.isUpper(id));
0959       isStack.push_back(tTopo.stack(id) ==
0960                         0);  // equivalent to *IsDoubleSide() but without the hardcoded layer+ring requirements
0961     }
0962 
0963     void clear() {
0964       isLower.clear();
0965       isUpper.clear();
0966       isStack.clear();
0967     }
0968 
0969   private:
0970     std::vector<unsigned short> isLower;
0971     std::vector<unsigned short> isUpper;
0972     std::vector<unsigned short> isStack;
0973   };
0974 #undef BOOK
0975 
0976   using DetIdPixel = CombineDetId<DetIdCommon, DetIdPixelOnly>;
0977   using DetIdStrip = CombineDetId<DetIdCommon, DetIdOTCommon, DetIdStripOnly>;
0978   using DetIdPhase2OT = CombineDetId<DetIdCommon, DetIdOTCommon, DetIdPhase2OTOnly>;
0979   using DetIdAll = CombineDetId<DetIdCommon, DetIdPixelOnly, DetIdOTCommon, DetIdStripOnly>;
0980   using DetIdAllPhase2 = CombineDetId<DetIdCommon, DetIdPixelOnly, DetIdOTCommon, DetIdPhase2OTOnly>;
0981 
0982   // event
0983   edm::RunNumber_t ev_run;
0984   edm::LuminosityBlockNumber_t ev_lumi;
0985   edm::EventNumber_t ev_event;
0986 
0987   ////////////////////
0988   // tracks
0989   // (first) index runs through tracks
0990   std::vector<float> trk_px;
0991   std::vector<float> trk_py;
0992   std::vector<float> trk_pz;
0993   std::vector<float> trk_pt;
0994   std::vector<float> trk_inner_px;
0995   std::vector<float> trk_inner_py;
0996   std::vector<float> trk_inner_pz;
0997   std::vector<float> trk_inner_pt;
0998   std::vector<float> trk_outer_px;
0999   std::vector<float> trk_outer_py;
1000   std::vector<float> trk_outer_pz;
1001   std::vector<float> trk_outer_pt;
1002   std::vector<float> trk_eta;
1003   std::vector<float> trk_lambda;
1004   std::vector<float> trk_cotTheta;
1005   std::vector<float> trk_phi;
1006   std::vector<float> trk_dxy;
1007   std::vector<float> trk_dz;
1008   std::vector<float> trk_dxyPV;
1009   std::vector<float> trk_dzPV;
1010   std::vector<float> trk_dxyClosestPV;
1011   std::vector<float> trk_dzClosestPV;
1012   std::vector<float> trk_ptErr;
1013   std::vector<float> trk_etaErr;
1014   std::vector<float> trk_lambdaErr;
1015   std::vector<float> trk_phiErr;
1016   std::vector<float> trk_dxyErr;
1017   std::vector<float> trk_dzErr;
1018   std::vector<float> trk_refpoint_x;
1019   std::vector<float> trk_refpoint_y;
1020   std::vector<float> trk_refpoint_z;
1021   std::vector<float> trk_nChi2;
1022   std::vector<float> trk_nChi2_1Dmod;
1023   std::vector<float> trk_ndof;
1024   std::vector<std::vector<float>> trk_mvas;
1025   std::vector<std::vector<unsigned short>> trk_qualityMasks;
1026   std::vector<int> trk_q;
1027   std::vector<unsigned int> trk_nValid;
1028   std::vector<unsigned int> trk_nLost;
1029   std::vector<unsigned int> trk_nInactive;
1030   std::vector<unsigned int> trk_nPixel;
1031   std::vector<unsigned int> trk_nStrip;
1032   std::vector<unsigned int> trk_nOuterLost;
1033   std::vector<unsigned int> trk_nInnerLost;
1034   std::vector<unsigned int> trk_nOuterInactive;
1035   std::vector<unsigned int> trk_nInnerInactive;
1036   std::vector<unsigned int> trk_nPixelLay;
1037   std::vector<unsigned int> trk_nStripLay;
1038   std::vector<unsigned int> trk_n3DLay;
1039   std::vector<unsigned int> trk_nLostLay;
1040   std::vector<unsigned int> trk_nCluster;
1041   std::vector<unsigned int> trk_algo;
1042   std::vector<unsigned int> trk_originalAlgo;
1043   std::vector<decltype(reco::TrackBase().algoMaskUL())> trk_algoMask;
1044   std::vector<unsigned short> trk_stopReason;
1045   std::vector<short> trk_isHP;
1046   std::vector<int> trk_seedIdx;
1047   std::vector<int> trk_vtxIdx;
1048   std::vector<short> trk_isTrue;
1049   std::vector<int> trk_bestSimTrkIdx;
1050   std::vector<float> trk_bestSimTrkShareFrac;
1051   std::vector<float> trk_bestSimTrkShareFracSimDenom;
1052   std::vector<float> trk_bestSimTrkShareFracSimClusterDenom;
1053   std::vector<float> trk_bestSimTrkNChi2;
1054   std::vector<int> trk_bestFromFirstHitSimTrkIdx;
1055   std::vector<float> trk_bestFromFirstHitSimTrkShareFrac;
1056   std::vector<float> trk_bestFromFirstHitSimTrkShareFracSimDenom;
1057   std::vector<float> trk_bestFromFirstHitSimTrkShareFracSimClusterDenom;
1058   std::vector<float> trk_bestFromFirstHitSimTrkNChi2;
1059   std::vector<std::vector<float>> trk_simTrkShareFrac;  // second index runs through matched TrackingParticles
1060   std::vector<std::vector<float>> trk_simTrkNChi2;      // second index runs through matched TrackingParticles
1061   std::vector<std::vector<int>> trk_simTrkIdx;          // second index runs through matched TrackingParticles
1062   std::vector<std::vector<int>> trk_hitIdx;             // second index runs through hits
1063   std::vector<std::vector<int>> trk_hitType;            // second index runs through hits
1064   ////////////////////
1065   // track candidates
1066   // (first) index runs through track candidates
1067   std::vector<short> tcand_pca_valid;  // pca: propagated params at BS
1068   std::vector<float> tcand_pca_px;
1069   std::vector<float> tcand_pca_py;
1070   std::vector<float> tcand_pca_pz;
1071   std::vector<float> tcand_pca_pt;
1072   std::vector<float> tcand_pca_eta;
1073   std::vector<float> tcand_pca_phi;
1074   std::vector<float> tcand_pca_dxy;
1075   std::vector<float> tcand_pca_dz;
1076   std::vector<float> tcand_pca_ptErr;
1077   std::vector<float> tcand_pca_etaErr;
1078   std::vector<float> tcand_pca_lambdaErr;
1079   std::vector<float> tcand_pca_phiErr;
1080   std::vector<float> tcand_pca_dxyErr;
1081   std::vector<float> tcand_pca_dzErr;
1082   std::vector<float> tcand_px;  // from PState
1083   std::vector<float> tcand_py;
1084   std::vector<float> tcand_pz;
1085   std::vector<float> tcand_pt;
1086   std::vector<float> tcand_x;
1087   std::vector<float> tcand_y;
1088   std::vector<float> tcand_z;
1089   std::vector<float> tcand_qbpErr;
1090   std::vector<float> tcand_lambdaErr;
1091   std::vector<float> tcand_phiErr;
1092   std::vector<float> tcand_xtErr;
1093   std::vector<float> tcand_ytErr;
1094   std::vector<float> tcand_ndof;  // sum(hit.dimention) - 5
1095   std::vector<int> tcand_q;
1096   std::vector<unsigned int> tcand_nValid;
1097   std::vector<unsigned int> tcand_nPixel;
1098   std::vector<unsigned int> tcand_nStrip;
1099   std::vector<unsigned int> tcand_nCluster;
1100   std::vector<unsigned int> tcand_algo;
1101   std::vector<unsigned short> tcand_stopReason;
1102   std::vector<int> tcand_seedIdx;
1103   std::vector<int> tcand_vtxIdx;
1104   std::vector<short> tcand_isTrue;
1105   std::vector<int> tcand_bestSimTrkIdx;
1106   std::vector<float> tcand_bestSimTrkShareFrac;
1107   std::vector<float> tcand_bestSimTrkShareFracSimDenom;
1108   std::vector<float> tcand_bestSimTrkShareFracSimClusterDenom;
1109   std::vector<float> tcand_bestSimTrkNChi2;
1110   std::vector<int> tcand_bestFromFirstHitSimTrkIdx;
1111   std::vector<float> tcand_bestFromFirstHitSimTrkShareFrac;
1112   std::vector<float> tcand_bestFromFirstHitSimTrkShareFracSimDenom;
1113   std::vector<float> tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom;
1114   std::vector<float> tcand_bestFromFirstHitSimTrkNChi2;
1115   std::vector<std::vector<float>> tcand_simTrkShareFrac;  // second index runs through matched TrackingParticles
1116   std::vector<std::vector<float>> tcand_simTrkNChi2;      // second index runs through matched TrackingParticles
1117   std::vector<std::vector<int>> tcand_simTrkIdx;          // second index runs through matched TrackingParticles
1118   std::vector<std::vector<int>> tcand_hitIdx;             // second index runs through hits
1119   std::vector<std::vector<int>> tcand_hitType;            // second index runs through hits
1120   ////////////////////
1121   // sim tracks
1122   // (first) index runs through TrackingParticles
1123   std::vector<int> sim_event;
1124   std::vector<int> sim_bunchCrossing;
1125   std::vector<int> sim_pdgId;
1126   std::vector<std::vector<int>> sim_genPdgIds;
1127   std::vector<int> sim_isFromBHadron;
1128   std::vector<float> sim_px;
1129   std::vector<float> sim_py;
1130   std::vector<float> sim_pz;
1131   std::vector<float> sim_pt;
1132   std::vector<float> sim_eta;
1133   std::vector<float> sim_phi;
1134   std::vector<float> sim_pca_pt;
1135   std::vector<float> sim_pca_eta;
1136   std::vector<float> sim_pca_lambda;
1137   std::vector<float> sim_pca_cotTheta;
1138   std::vector<float> sim_pca_phi;
1139   std::vector<float> sim_pca_dxy;
1140   std::vector<float> sim_pca_dz;
1141   std::vector<int> sim_q;
1142   // numbers of sim hits/layers
1143   std::vector<unsigned int> sim_nValid;
1144   std::vector<unsigned int> sim_nPixel;
1145   std::vector<unsigned int> sim_nStrip;
1146   std::vector<unsigned int> sim_nLay;
1147   std::vector<unsigned int> sim_nPixelLay;
1148   std::vector<unsigned int> sim_n3DLay;
1149   // number of sim hits as calculated in TrackingTruthAccumulator
1150   std::vector<unsigned int> sim_nTrackerHits;
1151   // number of clusters associated to TP
1152   std::vector<unsigned int> sim_nRecoClusters;
1153   // links to other objects
1154   std::vector<std::vector<int>> sim_trkIdx;          // second index runs through matched tracks
1155   std::vector<std::vector<float>> sim_trkShareFrac;  // second index runs through matched tracks
1156   std::vector<std::vector<int>> sim_seedIdx;         // second index runs through matched seeds
1157   std::vector<int> sim_parentVtxIdx;
1158   std::vector<std::vector<int>> sim_decayVtxIdx;  // second index runs through decay vertices
1159   std::vector<std::vector<int>> sim_simHitIdx;    // second index runs through SimHits
1160   ////////////////////
1161   // pixel hits
1162   // (first) index runs through hits
1163   std::vector<short> pix_isBarrel;
1164   DetIdPixel pix_detId;
1165   std::vector<std::vector<int>> pix_trkIdx;            // second index runs through tracks containing this hit
1166   std::vector<std::vector<float>> pix_onTrk_x;         // second index runs through tracks containing this hit
1167   std::vector<std::vector<float>> pix_onTrk_y;         // same for all pix_onTrk_*
1168   std::vector<std::vector<float>> pix_onTrk_z;         //
1169   std::vector<std::vector<float>> pix_onTrk_xx;        //
1170   std::vector<std::vector<float>> pix_onTrk_xy;        //
1171   std::vector<std::vector<float>> pix_onTrk_yy;        //
1172   std::vector<std::vector<float>> pix_onTrk_yz;        //
1173   std::vector<std::vector<float>> pix_onTrk_zz;        //
1174   std::vector<std::vector<float>> pix_onTrk_zx;        //
1175   std::vector<std::vector<int>> pix_tcandIdx;          // second index runs through candidates containing this hit
1176   std::vector<std::vector<int>> pix_seeIdx;            // second index runs through seeds containing this hit
1177   std::vector<std::vector<int>> pix_simHitIdx;         // second index runs through SimHits inducing this hit
1178   std::vector<std::vector<float>> pix_xySignificance;  // second index runs through SimHits inducing this hit
1179   std::vector<std::vector<float>> pix_chargeFraction;  // second index runs through SimHits inducing this hit
1180   std::vector<unsigned short> pix_simType;
1181   std::vector<float> pix_x;
1182   std::vector<float> pix_y;
1183   std::vector<float> pix_z;
1184   std::vector<float> pix_xx;
1185   std::vector<float> pix_xy;
1186   std::vector<float> pix_yy;
1187   std::vector<float> pix_yz;
1188   std::vector<float> pix_zz;
1189   std::vector<float> pix_zx;
1190   std::vector<float>
1191       pix_radL;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
1192   std::vector<float> pix_bbxi;
1193   std::vector<int> pix_clustSizeCol;
1194   std::vector<int> pix_clustSizeRow;
1195   std::vector<uint64_t> pix_usedMask;
1196   ////////////////////
1197   // strip hits
1198   // (first) index runs through hits
1199   std::vector<short> str_isBarrel;
1200   DetIdStrip str_detId;
1201   std::vector<std::vector<int>> str_trkIdx;            // second index runs through tracks containing this hit
1202   std::vector<std::vector<float>> str_onTrk_x;         // second index runs through tracks containing this hit
1203   std::vector<std::vector<float>> str_onTrk_y;         // same for all pix_onTrk_*
1204   std::vector<std::vector<float>> str_onTrk_z;         //
1205   std::vector<std::vector<float>> str_onTrk_xx;        //
1206   std::vector<std::vector<float>> str_onTrk_xy;        //
1207   std::vector<std::vector<float>> str_onTrk_yy;        //
1208   std::vector<std::vector<float>> str_onTrk_yz;        //
1209   std::vector<std::vector<float>> str_onTrk_zz;        //
1210   std::vector<std::vector<float>> str_onTrk_zx;        //
1211   std::vector<std::vector<int>> str_tcandIdx;          // second index runs through candidates containing this hit
1212   std::vector<std::vector<int>> str_seeIdx;            // second index runs through seeds containing this hit
1213   std::vector<std::vector<int>> str_simHitIdx;         // second index runs through SimHits inducing this hit
1214   std::vector<std::vector<float>> str_xySignificance;  // second index runs through SimHits inducing this hit
1215   std::vector<std::vector<float>> str_chargeFraction;  // second index runs through SimHits inducing this hit
1216   std::vector<unsigned short> str_simType;
1217   std::vector<float> str_x;
1218   std::vector<float> str_y;
1219   std::vector<float> str_z;
1220   std::vector<float> str_xx;
1221   std::vector<float> str_xy;
1222   std::vector<float> str_yy;
1223   std::vector<float> str_yz;
1224   std::vector<float> str_zz;
1225   std::vector<float> str_zx;
1226   std::vector<float>
1227       str_radL;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
1228   std::vector<float> str_bbxi;
1229   std::vector<float> str_chargePerCM;
1230   std::vector<int> str_clustSize;
1231   std::vector<uint64_t> str_usedMask;
1232   ////////////////////
1233   // strip matched hits
1234   // (first) index runs through hits
1235   std::vector<short> glu_isBarrel;
1236   DetIdStrip glu_detId;
1237   std::vector<int> glu_monoIdx;
1238   std::vector<int> glu_stereoIdx;
1239   std::vector<std::vector<int>> glu_seeIdx;  // second index runs through seeds containing this hit
1240   std::vector<float> glu_x;
1241   std::vector<float> glu_y;
1242   std::vector<float> glu_z;
1243   std::vector<float> glu_xx;
1244   std::vector<float> glu_xy;
1245   std::vector<float> glu_yy;
1246   std::vector<float> glu_yz;
1247   std::vector<float> glu_zz;
1248   std::vector<float> glu_zx;
1249   std::vector<float>
1250       glu_radL;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
1251   std::vector<float> glu_bbxi;
1252   std::vector<float> glu_chargePerCM;
1253   std::vector<int> glu_clustSizeMono;
1254   std::vector<int> glu_clustSizeStereo;
1255   std::vector<uint64_t> glu_usedMaskMono;
1256   std::vector<uint64_t> glu_usedMaskStereo;
1257   ////////////////////
1258   // phase2 Outer Tracker hits
1259   // (first) index runs through hits
1260   std::vector<short> ph2_isBarrel;
1261   DetIdPhase2OT ph2_detId;
1262   std::vector<std::vector<int>> ph2_trkIdx;            // second index runs through tracks containing this hit
1263   std::vector<std::vector<float>> ph2_onTrk_x;         // second index runs through tracks containing this hit
1264   std::vector<std::vector<float>> ph2_onTrk_y;         // same for all pix_onTrk_*
1265   std::vector<std::vector<float>> ph2_onTrk_z;         //
1266   std::vector<std::vector<float>> ph2_onTrk_xx;        //
1267   std::vector<std::vector<float>> ph2_onTrk_xy;        //
1268   std::vector<std::vector<float>> ph2_onTrk_yy;        //
1269   std::vector<std::vector<float>> ph2_onTrk_yz;        //
1270   std::vector<std::vector<float>> ph2_onTrk_zz;        //
1271   std::vector<std::vector<float>> ph2_onTrk_zx;        //
1272   std::vector<std::vector<int>> ph2_tcandIdx;          // second index runs through candidates containing this hit
1273   std::vector<std::vector<int>> ph2_seeIdx;            // second index runs through seeds containing this hit
1274   std::vector<std::vector<int>> ph2_simHitIdx;         // second index runs through SimHits inducing this hit
1275   std::vector<std::vector<float>> ph2_xySignificance;  // second index runs through SimHits inducing this hit
1276   //std::vector<std::vector<float>> ph2_chargeFraction; // Not supported at the moment for Phase2
1277   std::vector<unsigned short> ph2_simType;
1278   std::vector<float> ph2_x;
1279   std::vector<float> ph2_y;
1280   std::vector<float> ph2_z;
1281   std::vector<float> ph2_xx;
1282   std::vector<float> ph2_xy;
1283   std::vector<float> ph2_yy;
1284   std::vector<float> ph2_yz;
1285   std::vector<float> ph2_zz;
1286   std::vector<float> ph2_zx;
1287   std::vector<float>
1288       ph2_radL;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
1289   std::vector<float> ph2_bbxi;
1290 
1291   ////////////////////
1292   // invalid (missing/inactive/etc) hits
1293   // (first) index runs through hits
1294   std::vector<short> inv_isBarrel;
1295   DetIdAll inv_detId;
1296   DetIdAllPhase2 inv_detId_phase2;
1297   std::vector<unsigned short> inv_type;
1298   ////////////////////
1299   // sim hits
1300   // (first) index runs through hits
1301   DetIdAll simhit_detId;
1302   DetIdAllPhase2 simhit_detId_phase2;
1303   std::vector<float> simhit_x;
1304   std::vector<float> simhit_y;
1305   std::vector<float> simhit_z;
1306   std::vector<float> simhit_px;
1307   std::vector<float> simhit_py;
1308   std::vector<float> simhit_pz;
1309   std::vector<int> simhit_particle;
1310   std::vector<short> simhit_process;
1311   std::vector<float> simhit_eloss;
1312   std::vector<float> simhit_tof;
1313   //std::vector<unsigned int> simhit_simTrackId; // can be useful for debugging, but not much of general interest
1314   std::vector<int> simhit_simTrkIdx;
1315   std::vector<std::vector<int>> simhit_hitIdx;   // second index runs through induced reco hits
1316   std::vector<std::vector<int>> simhit_hitType;  // second index runs through induced reco hits
1317   ////////////////////
1318   // beam spot
1319   float bsp_x;
1320   float bsp_y;
1321   float bsp_z;
1322   float bsp_sigmax;
1323   float bsp_sigmay;
1324   float bsp_sigmaz;
1325   ////////////////////
1326   // seeds
1327   // (first) index runs through seeds
1328   std::vector<short> see_fitok;
1329   std::vector<float> see_px;
1330   std::vector<float> see_py;
1331   std::vector<float> see_pz;
1332   std::vector<float> see_pt;
1333   std::vector<float> see_eta;
1334   std::vector<float> see_phi;
1335   std::vector<float> see_dxy;
1336   std::vector<float> see_dz;
1337   std::vector<float> see_ptErr;
1338   std::vector<float> see_etaErr;
1339   std::vector<float> see_phiErr;
1340   std::vector<float> see_dxyErr;
1341   std::vector<float> see_dzErr;
1342   std::vector<float> see_chi2;
1343   std::vector<float> see_statePt;
1344   std::vector<float> see_stateTrajX;
1345   std::vector<float> see_stateTrajY;
1346   std::vector<float> see_stateTrajPx;
1347   std::vector<float> see_stateTrajPy;
1348   std::vector<float> see_stateTrajPz;
1349   std::vector<float> see_stateTrajGlbX;
1350   std::vector<float> see_stateTrajGlbY;
1351   std::vector<float> see_stateTrajGlbZ;
1352   std::vector<float> see_stateTrajGlbPx;
1353   std::vector<float> see_stateTrajGlbPy;
1354   std::vector<float> see_stateTrajGlbPz;
1355   std::vector<std::vector<float>> see_stateCurvCov;
1356   std::vector<int> see_q;
1357   std::vector<unsigned int> see_nValid;
1358   std::vector<unsigned int> see_nPixel;
1359   std::vector<unsigned int> see_nGlued;
1360   std::vector<unsigned int> see_nStrip;
1361   std::vector<unsigned int> see_nPhase2OT;
1362   std::vector<unsigned int> see_nCluster;
1363   std::vector<unsigned int> see_algo;
1364   std::vector<unsigned short> see_stopReason;
1365   std::vector<unsigned short> see_nCands;
1366   std::vector<int> see_trkIdx;
1367   std::vector<int> see_tcandIdx;
1368   std::vector<short> see_isTrue;
1369   std::vector<int> see_bestSimTrkIdx;
1370   std::vector<float> see_bestSimTrkShareFrac;
1371   std::vector<int> see_bestFromFirstHitSimTrkIdx;
1372   std::vector<float> see_bestFromFirstHitSimTrkShareFrac;
1373   std::vector<std::vector<float>> see_simTrkShareFrac;  // second index runs through matched TrackingParticles
1374   std::vector<std::vector<int>> see_simTrkIdx;          // second index runs through matched TrackingParticles
1375   std::vector<std::vector<int>> see_hitIdx;             // second index runs through hits
1376   std::vector<std::vector<int>> see_hitType;            // second index runs through hits
1377   //seed algo offset, index runs through iterations
1378   std::vector<unsigned int> see_offset;
1379 
1380   ////////////////////
1381   // Vertices
1382   // (first) index runs through vertices
1383   std::vector<float> vtx_x;
1384   std::vector<float> vtx_y;
1385   std::vector<float> vtx_z;
1386   std::vector<float> vtx_xErr;
1387   std::vector<float> vtx_yErr;
1388   std::vector<float> vtx_zErr;
1389   std::vector<float> vtx_ndof;
1390   std::vector<float> vtx_chi2;
1391   std::vector<short> vtx_fake;
1392   std::vector<short> vtx_valid;
1393   std::vector<std::vector<int>> vtx_trkIdx;  // second index runs through tracks used in the vertex fit
1394 
1395   ////////////////////
1396   // Tracking vertices
1397   // (first) index runs through TrackingVertices
1398   std::vector<int> simvtx_event;
1399   std::vector<int> simvtx_bunchCrossing;
1400   std::vector<unsigned int> simvtx_processType;  // only from first SimVertex of TrackingVertex
1401   std::vector<float> simvtx_x;
1402   std::vector<float> simvtx_y;
1403   std::vector<float> simvtx_z;
1404   std::vector<std::vector<int>> simvtx_sourceSimIdx;    // second index runs through source TrackingParticles
1405   std::vector<std::vector<int>> simvtx_daughterSimIdx;  // second index runs through daughter TrackingParticles
1406   std::vector<int> simpv_idx;
1407 };
1408 
1409 //
1410 // constructors and destructor
1411 //
1412 TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig)
1413     : mfToken_(esConsumes()),
1414       tTopoToken_(esConsumes()),
1415       tGeomToken_(esConsumes()),
1416       trackToken_(consumes<edm::View<reco::Track>>(iConfig.getUntrackedParameter<edm::InputTag>("tracks"))),
1417       clusterTPMapToken_(consumes<ClusterTPAssociation>(iConfig.getUntrackedParameter<edm::InputTag>("clusterTPMap"))),
1418       simHitTPMapToken_(consumes<SimHitTPAssociationProducer::SimHitTPAssociationList>(
1419           iConfig.getUntrackedParameter<edm::InputTag>("simHitTPMap"))),
1420       trackAssociatorToken_(consumes<reco::TrackToTrackingParticleAssociator>(
1421           iConfig.getUntrackedParameter<edm::InputTag>("trackAssociator"))),
1422       pixelSimLinkToken_(consumes<edm::DetSetVector<PixelDigiSimLink>>(
1423           iConfig.getUntrackedParameter<edm::InputTag>("pixelDigiSimLink"))),
1424       stripSimLinkToken_(consumes<edm::DetSetVector<StripDigiSimLink>>(
1425           iConfig.getUntrackedParameter<edm::InputTag>("stripDigiSimLink"))),
1426       siphase2OTSimLinksToken_(consumes<edm::DetSetVector<PixelDigiSimLink>>(
1427           iConfig.getUntrackedParameter<edm::InputTag>("phase2OTSimLink"))),
1428       includeStripHits_(!iConfig.getUntrackedParameter<edm::InputTag>("stripDigiSimLink").label().empty()),
1429       includePhase2OTHits_(!iConfig.getUntrackedParameter<edm::InputTag>("phase2OTSimLink").label().empty()),
1430       beamSpotToken_(consumes<reco::BeamSpot>(iConfig.getUntrackedParameter<edm::InputTag>("beamSpot"))),
1431       pixelRecHitToken_(
1432           consumes<SiPixelRecHitCollection>(iConfig.getUntrackedParameter<edm::InputTag>("pixelRecHits"))),
1433       stripRphiRecHitToken_(
1434           consumes<SiStripRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripRphiRecHits"))),
1435       stripStereoRecHitToken_(
1436           consumes<SiStripRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripStereoRecHits"))),
1437       stripMatchedRecHitToken_(consumes<SiStripMatchedRecHit2DCollection>(
1438           iConfig.getUntrackedParameter<edm::InputTag>("stripMatchedRecHits"))),
1439       phase2OTRecHitToken_(consumes<Phase2TrackerRecHit1DCollectionNew>(
1440           iConfig.getUntrackedParameter<edm::InputTag>("phase2OTRecHits"))),
1441       vertexToken_(consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertices"))),
1442       trackingVertexToken_(
1443           consumes<TrackingVertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trackingVertices"))),
1444       tpNLayersToken_(consumes<edm::ValueMap<unsigned int>>(
1445           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNlayers"))),
1446       tpNPixelLayersToken_(consumes<edm::ValueMap<unsigned int>>(
1447           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNpixellayers"))),
1448       tpNStripStereoLayersToken_(consumes<edm::ValueMap<unsigned int>>(
1449           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNstripstereolayers"))),
1450       includeSeeds_(iConfig.getUntrackedParameter<bool>("includeSeeds")),
1451       includeTrackCandidates_(iConfig.getUntrackedParameter<bool>("includeTrackCandidates")),
1452       addSeedCurvCov_(iConfig.getUntrackedParameter<bool>("addSeedCurvCov")),
1453       includeAllHits_(iConfig.getUntrackedParameter<bool>("includeAllHits")),
1454       includeOnTrackHitData_(iConfig.getUntrackedParameter<bool>("includeOnTrackHitData")),
1455       includeMVA_(iConfig.getUntrackedParameter<bool>("includeMVA")),
1456       includeTrackingParticles_(iConfig.getUntrackedParameter<bool>("includeTrackingParticles")),
1457       includeOOT_(iConfig.getUntrackedParameter<bool>("includeOOT")),
1458       keepEleSimHits_(iConfig.getUntrackedParameter<bool>("keepEleSimHits")),
1459       saveSimHitsP3_(iConfig.getUntrackedParameter<bool>("saveSimHitsP3")),
1460       simHitBySignificance_(iConfig.getUntrackedParameter<bool>("simHitBySignificance")),
1461       parametersDefiner_(iConfig.getUntrackedParameter<edm::InputTag>("beamSpot"), consumesCollector()) {
1462   if (includeSeeds_) {
1463     seedTokens_ =
1464         edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("seedTracks"),
1465                               [&](const edm::InputTag& tag) { return consumes<edm::View<reco::Track>>(tag); });
1466     seedStopInfoTokens_ =
1467         edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("trackCandidates"),
1468                               [&](const edm::InputTag& tag) { return consumes<std::vector<SeedStopInfo>>(tag); });
1469     if (seedTokens_.size() != seedStopInfoTokens_.size()) {
1470       throw cms::Exception("Configuration") << "Got " << seedTokens_.size() << " seed collections, but "
1471                                             << seedStopInfoTokens_.size() << " track candidate collections";
1472     }
1473   }
1474   if (includeTrackCandidates_)
1475     candidateTokens_ =
1476         edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("trackCandidates"),
1477                               [&](const edm::InputTag& tag) { return consumes<TrackCandidateCollection>(tag); });
1478 
1479   if (includeAllHits_) {
1480     if (includeStripHits_ && includePhase2OTHits_) {
1481       throw cms::Exception("Configuration")
1482           << "Both stripDigiSimLink and phase2OTSimLink are set, please set only either one (this information is used "
1483              "to infer if you're running phase0/1 or phase2 detector)";
1484     }
1485     if (!includeStripHits_ && !includePhase2OTHits_) {
1486       throw cms::Exception("Configuration")
1487           << "Neither stripDigiSimLink or phase2OTSimLink are set, please set either one.";
1488     }
1489 
1490     auto const& maskVPset = iConfig.getUntrackedParameterSetVector("clusterMasks");
1491     pixelUseMaskTokens_.reserve(maskVPset.size());
1492     stripUseMaskTokens_.reserve(maskVPset.size());
1493     for (auto const& mask : maskVPset) {
1494       auto index = mask.getUntrackedParameter<unsigned int>("index");
1495       assert(index < 64);
1496       pixelUseMaskTokens_.emplace_back(index,
1497                                        consumes<PixelMaskContainer>(mask.getUntrackedParameter<edm::InputTag>("src")));
1498       if (includeStripHits_)
1499         stripUseMaskTokens_.emplace_back(
1500             index, consumes<StripMaskContainer>(mask.getUntrackedParameter<edm::InputTag>("src")));
1501     }
1502   }
1503 
1504   const bool tpRef = iConfig.getUntrackedParameter<bool>("trackingParticlesRef");
1505   const auto tpTag = iConfig.getUntrackedParameter<edm::InputTag>("trackingParticles");
1506   if (tpRef) {
1507     trackingParticleRefToken_ = consumes<TrackingParticleRefVector>(tpTag);
1508   } else {
1509     trackingParticleToken_ = consumes<TrackingParticleCollection>(tpTag);
1510   }
1511 
1512   tracer_.depth(-2);  // as in SimTracker/TrackHistory/src/TrackClassifier.cc
1513 
1514   if (includeMVA_) {
1515     mvaQualityCollectionTokens_ = edm::vector_transform(
1516         iConfig.getUntrackedParameter<std::vector<std::string>>("trackMVAs"), [&](const std::string& tag) {
1517           return std::make_tuple(consumes<MVACollection>(edm::InputTag(tag, "MVAValues")),
1518                                  consumes<QualityMaskCollection>(edm::InputTag(tag, "QualityMasks")));
1519         });
1520   }
1521 
1522   usesResource(TFileService::kSharedResource);
1523   edm::Service<TFileService> fs;
1524   t = fs->make<TTree>("tree", "tree");
1525 
1526   t->Branch("event", &ev_event);
1527   t->Branch("lumi", &ev_lumi);
1528   t->Branch("run", &ev_run);
1529 
1530   //tracks
1531   t->Branch("trk_px", &trk_px);
1532   t->Branch("trk_py", &trk_py);
1533   t->Branch("trk_pz", &trk_pz);
1534   t->Branch("trk_pt", &trk_pt);
1535   t->Branch("trk_inner_px", &trk_inner_px);
1536   t->Branch("trk_inner_py", &trk_inner_py);
1537   t->Branch("trk_inner_pz", &trk_inner_pz);
1538   t->Branch("trk_inner_pt", &trk_inner_pt);
1539   t->Branch("trk_outer_px", &trk_outer_px);
1540   t->Branch("trk_outer_py", &trk_outer_py);
1541   t->Branch("trk_outer_pz", &trk_outer_pz);
1542   t->Branch("trk_outer_pt", &trk_outer_pt);
1543   t->Branch("trk_eta", &trk_eta);
1544   t->Branch("trk_lambda", &trk_lambda);
1545   t->Branch("trk_cotTheta", &trk_cotTheta);
1546   t->Branch("trk_phi", &trk_phi);
1547   t->Branch("trk_dxy", &trk_dxy);
1548   t->Branch("trk_dz", &trk_dz);
1549   t->Branch("trk_dxyPV", &trk_dxyPV);
1550   t->Branch("trk_dzPV", &trk_dzPV);
1551   t->Branch("trk_dxyClosestPV", &trk_dxyClosestPV);
1552   t->Branch("trk_dzClosestPV", &trk_dzClosestPV);
1553   t->Branch("trk_ptErr", &trk_ptErr);
1554   t->Branch("trk_etaErr", &trk_etaErr);
1555   t->Branch("trk_lambdaErr", &trk_lambdaErr);
1556   t->Branch("trk_phiErr", &trk_phiErr);
1557   t->Branch("trk_dxyErr", &trk_dxyErr);
1558   t->Branch("trk_dzErr", &trk_dzErr);
1559   t->Branch("trk_refpoint_x", &trk_refpoint_x);
1560   t->Branch("trk_refpoint_y", &trk_refpoint_y);
1561   t->Branch("trk_refpoint_z", &trk_refpoint_z);
1562   t->Branch("trk_nChi2", &trk_nChi2);
1563   t->Branch("trk_nChi2_1Dmod", &trk_nChi2_1Dmod);
1564   t->Branch("trk_ndof", &trk_ndof);
1565   if (includeMVA_) {
1566     trk_mvas.resize(mvaQualityCollectionTokens_.size());
1567     trk_qualityMasks.resize(mvaQualityCollectionTokens_.size());
1568     if (!trk_mvas.empty()) {
1569       t->Branch("trk_mva", &(trk_mvas[0]));
1570       t->Branch("trk_qualityMask", &(trk_qualityMasks[0]));
1571       for (size_t i = 1; i < trk_mvas.size(); ++i) {
1572         t->Branch(("trk_mva" + std::to_string(i + 1)).c_str(), &(trk_mvas[i]));
1573         t->Branch(("trk_qualityMask" + std::to_string(i + 1)).c_str(), &(trk_qualityMasks[i]));
1574       }
1575     }
1576   }
1577   t->Branch("trk_q", &trk_q);
1578   t->Branch("trk_nValid", &trk_nValid);
1579   t->Branch("trk_nLost", &trk_nLost);
1580   t->Branch("trk_nInactive", &trk_nInactive);
1581   t->Branch("trk_nPixel", &trk_nPixel);
1582   t->Branch("trk_nStrip", &trk_nStrip);
1583   t->Branch("trk_nOuterLost", &trk_nOuterLost);
1584   t->Branch("trk_nInnerLost", &trk_nInnerLost);
1585   t->Branch("trk_nOuterInactive", &trk_nOuterInactive);
1586   t->Branch("trk_nInnerInactive", &trk_nInnerInactive);
1587   t->Branch("trk_nPixelLay", &trk_nPixelLay);
1588   t->Branch("trk_nStripLay", &trk_nStripLay);
1589   t->Branch("trk_n3DLay", &trk_n3DLay);
1590   t->Branch("trk_nLostLay", &trk_nLostLay);
1591   t->Branch("trk_nCluster", &trk_nCluster);
1592   t->Branch("trk_algo", &trk_algo);
1593   t->Branch("trk_originalAlgo", &trk_originalAlgo);
1594   t->Branch("trk_algoMask", &trk_algoMask);
1595   t->Branch("trk_stopReason", &trk_stopReason);
1596   t->Branch("trk_isHP", &trk_isHP);
1597   if (includeSeeds_) {
1598     t->Branch("trk_seedIdx", &trk_seedIdx);
1599   }
1600   t->Branch("trk_vtxIdx", &trk_vtxIdx);
1601   if (includeTrackingParticles_) {
1602     t->Branch("trk_simTrkIdx", &trk_simTrkIdx);
1603     t->Branch("trk_simTrkShareFrac", &trk_simTrkShareFrac);
1604     t->Branch("trk_simTrkNChi2", &trk_simTrkNChi2);
1605     t->Branch("trk_bestSimTrkIdx", &trk_bestSimTrkIdx);
1606     t->Branch("trk_bestFromFirstHitSimTrkIdx", &trk_bestFromFirstHitSimTrkIdx);
1607   } else {
1608     t->Branch("trk_isTrue", &trk_isTrue);
1609   }
1610   t->Branch("trk_bestSimTrkShareFrac", &trk_bestSimTrkShareFrac);
1611   t->Branch("trk_bestSimTrkShareFracSimDenom", &trk_bestSimTrkShareFracSimDenom);
1612   t->Branch("trk_bestSimTrkShareFracSimClusterDenom", &trk_bestSimTrkShareFracSimClusterDenom);
1613   t->Branch("trk_bestSimTrkNChi2", &trk_bestSimTrkNChi2);
1614   t->Branch("trk_bestFromFirstHitSimTrkShareFrac", &trk_bestFromFirstHitSimTrkShareFrac);
1615   t->Branch("trk_bestFromFirstHitSimTrkShareFracSimDenom", &trk_bestFromFirstHitSimTrkShareFracSimDenom);
1616   t->Branch("trk_bestFromFirstHitSimTrkShareFracSimClusterDenom", &trk_bestFromFirstHitSimTrkShareFracSimClusterDenom);
1617   t->Branch("trk_bestFromFirstHitSimTrkNChi2", &trk_bestFromFirstHitSimTrkNChi2);
1618   if (includeAllHits_) {
1619     t->Branch("trk_hitIdx", &trk_hitIdx);
1620     t->Branch("trk_hitType", &trk_hitType);
1621   }
1622   if (includeTrackCandidates_) {
1623     t->Branch("tcand_pca_valid", &tcand_pca_valid);
1624     t->Branch("tcand_pca_px", &tcand_pca_px);
1625     t->Branch("tcand_pca_py", &tcand_pca_py);
1626     t->Branch("tcand_pca_pz", &tcand_pca_pz);
1627     t->Branch("tcand_pca_pt", &tcand_pca_pt);
1628     t->Branch("tcand_pca_eta", &tcand_pca_eta);
1629     t->Branch("tcand_pca_phi", &tcand_pca_phi);
1630     t->Branch("tcand_pca_dxy", &tcand_pca_dxy);
1631     t->Branch("tcand_pca_dz", &tcand_pca_dz);
1632     t->Branch("tcand_pca_ptErr", &tcand_pca_ptErr);
1633     t->Branch("tcand_pca_etaErr", &tcand_pca_etaErr);
1634     t->Branch("tcand_pca_lambdaErr", &tcand_pca_lambdaErr);
1635     t->Branch("tcand_pca_phiErr", &tcand_pca_phiErr);
1636     t->Branch("tcand_pca_dxyErr", &tcand_pca_dxyErr);
1637     t->Branch("tcand_pca_dzErr", &tcand_pca_dzErr);
1638     t->Branch("tcand_px", &tcand_px);
1639     t->Branch("tcand_py", &tcand_py);
1640     t->Branch("tcand_pz", &tcand_pz);
1641     t->Branch("tcand_pt", &tcand_pt);
1642     t->Branch("tcand_x", &tcand_x);
1643     t->Branch("tcand_y", &tcand_y);
1644     t->Branch("tcand_z", &tcand_z);
1645     t->Branch("tcand_qbpErr", &tcand_qbpErr);
1646     t->Branch("tcand_lambdaErr", &tcand_lambdaErr);
1647     t->Branch("tcand_phiErr", &tcand_phiErr);
1648     t->Branch("tcand_xtErr", &tcand_xtErr);
1649     t->Branch("tcand_ytErr", &tcand_ytErr);
1650     t->Branch("tcand_q", &tcand_q);
1651     t->Branch("tcand_ndof", &tcand_ndof);
1652     t->Branch("tcand_nValid", &tcand_nValid);
1653     t->Branch("tcand_nPixel", &tcand_nPixel);
1654     t->Branch("tcand_nStrip", &tcand_nStrip);
1655     t->Branch("tcand_nCluster", &tcand_nCluster);
1656     t->Branch("tcand_algo", &tcand_algo);
1657     t->Branch("tcand_stopReason", &tcand_stopReason);
1658     if (includeSeeds_) {
1659       t->Branch("tcand_seedIdx", &tcand_seedIdx);
1660     }
1661     t->Branch("tcand_vtxIdx", &tcand_vtxIdx);
1662     if (includeTrackingParticles_) {
1663       t->Branch("tcand_simTrkIdx", &tcand_simTrkIdx);
1664       t->Branch("tcand_simTrkShareFrac", &tcand_simTrkShareFrac);
1665       t->Branch("tcand_simTrkNChi2", &tcand_simTrkNChi2);
1666       t->Branch("tcand_bestSimTrkIdx", &tcand_bestSimTrkIdx);
1667       t->Branch("tcand_bestFromFirstHitSimTrkIdx", &tcand_bestFromFirstHitSimTrkIdx);
1668     } else {
1669       t->Branch("tcand_isTrue", &tcand_isTrue);
1670     }
1671     t->Branch("tcand_bestSimTrkShareFrac", &tcand_bestSimTrkShareFrac);
1672     t->Branch("tcand_bestSimTrkShareFracSimDenom", &tcand_bestSimTrkShareFracSimDenom);
1673     t->Branch("tcand_bestSimTrkShareFracSimClusterDenom", &tcand_bestSimTrkShareFracSimClusterDenom);
1674     t->Branch("tcand_bestSimTrkNChi2", &tcand_bestSimTrkNChi2);
1675     t->Branch("tcand_bestFromFirstHitSimTrkShareFrac", &tcand_bestFromFirstHitSimTrkShareFrac);
1676     t->Branch("tcand_bestFromFirstHitSimTrkShareFracSimDenom", &tcand_bestFromFirstHitSimTrkShareFracSimDenom);
1677     t->Branch("tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom",
1678               &tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom);
1679     t->Branch("tcand_bestFromFirstHitSimTrkNChi2", &tcand_bestFromFirstHitSimTrkNChi2);
1680     if (includeAllHits_) {
1681       t->Branch("tcand_hitIdx", &tcand_hitIdx);
1682       t->Branch("tcand_hitType", &tcand_hitType);
1683     }
1684   }
1685   if (includeTrackingParticles_) {
1686     //sim tracks
1687     t->Branch("sim_event", &sim_event);
1688     t->Branch("sim_bunchCrossing", &sim_bunchCrossing);
1689     t->Branch("sim_pdgId", &sim_pdgId);
1690     t->Branch("sim_genPdgIds", &sim_genPdgIds);
1691     t->Branch("sim_isFromBHadron", &sim_isFromBHadron);
1692     t->Branch("sim_px", &sim_px);
1693     t->Branch("sim_py", &sim_py);
1694     t->Branch("sim_pz", &sim_pz);
1695     t->Branch("sim_pt", &sim_pt);
1696     t->Branch("sim_eta", &sim_eta);
1697     t->Branch("sim_phi", &sim_phi);
1698     t->Branch("sim_pca_pt", &sim_pca_pt);
1699     t->Branch("sim_pca_eta", &sim_pca_eta);
1700     t->Branch("sim_pca_lambda", &sim_pca_lambda);
1701     t->Branch("sim_pca_cotTheta", &sim_pca_cotTheta);
1702     t->Branch("sim_pca_phi", &sim_pca_phi);
1703     t->Branch("sim_pca_dxy", &sim_pca_dxy);
1704     t->Branch("sim_pca_dz", &sim_pca_dz);
1705     t->Branch("sim_q", &sim_q);
1706     t->Branch("sim_nValid", &sim_nValid);
1707     t->Branch("sim_nPixel", &sim_nPixel);
1708     t->Branch("sim_nStrip", &sim_nStrip);
1709     t->Branch("sim_nLay", &sim_nLay);
1710     t->Branch("sim_nPixelLay", &sim_nPixelLay);
1711     t->Branch("sim_n3DLay", &sim_n3DLay);
1712     t->Branch("sim_nTrackerHits", &sim_nTrackerHits);
1713     t->Branch("sim_nRecoClusters", &sim_nRecoClusters);
1714     t->Branch("sim_trkIdx", &sim_trkIdx);
1715     t->Branch("sim_trkShareFrac", &sim_trkShareFrac);
1716     if (includeSeeds_) {
1717       t->Branch("sim_seedIdx", &sim_seedIdx);
1718     }
1719     t->Branch("sim_parentVtxIdx", &sim_parentVtxIdx);
1720     t->Branch("sim_decayVtxIdx", &sim_decayVtxIdx);
1721     if (includeAllHits_) {
1722       t->Branch("sim_simHitIdx", &sim_simHitIdx);
1723     }
1724   }
1725   if (includeAllHits_) {
1726     //pixels
1727     t->Branch("pix_isBarrel", &pix_isBarrel);
1728     pix_detId.book("pix", t);
1729     t->Branch("pix_trkIdx", &pix_trkIdx);
1730     if (includeOnTrackHitData_) {
1731       t->Branch("pix_onTrk_x", &pix_onTrk_x);
1732       t->Branch("pix_onTrk_y", &pix_onTrk_y);
1733       t->Branch("pix_onTrk_z", &pix_onTrk_z);
1734       t->Branch("pix_onTrk_xx", &pix_onTrk_xx);
1735       t->Branch("pix_onTrk_xy", &pix_onTrk_xy);
1736       t->Branch("pix_onTrk_yy", &pix_onTrk_yy);
1737       t->Branch("pix_onTrk_yz", &pix_onTrk_yz);
1738       t->Branch("pix_onTrk_zz", &pix_onTrk_zz);
1739       t->Branch("pix_onTrk_zx", &pix_onTrk_zx);
1740     }
1741     if (includeTrackCandidates_)
1742       t->Branch("pix_tcandIdx", &pix_tcandIdx);
1743     if (includeSeeds_) {
1744       t->Branch("pix_seeIdx", &pix_seeIdx);
1745     }
1746     if (includeTrackingParticles_) {
1747       t->Branch("pix_simHitIdx", &pix_simHitIdx);
1748       if (simHitBySignificance_) {
1749         t->Branch("pix_xySignificance", &pix_xySignificance);
1750       }
1751       t->Branch("pix_chargeFraction", &pix_chargeFraction);
1752       t->Branch("pix_simType", &pix_simType);
1753     }
1754     t->Branch("pix_x", &pix_x);
1755     t->Branch("pix_y", &pix_y);
1756     t->Branch("pix_z", &pix_z);
1757     t->Branch("pix_xx", &pix_xx);
1758     t->Branch("pix_xy", &pix_xy);
1759     t->Branch("pix_yy", &pix_yy);
1760     t->Branch("pix_yz", &pix_yz);
1761     t->Branch("pix_zz", &pix_zz);
1762     t->Branch("pix_zx", &pix_zx);
1763     t->Branch("pix_radL", &pix_radL);
1764     t->Branch("pix_bbxi", &pix_bbxi);
1765     t->Branch("pix_clustSizeCol", &pix_clustSizeCol);
1766     t->Branch("pix_clustSizeRow", &pix_clustSizeRow);
1767     t->Branch("pix_usedMask", &pix_usedMask);
1768     //strips
1769     if (includeStripHits_) {
1770       t->Branch("str_isBarrel", &str_isBarrel);
1771       str_detId.book("str", t);
1772       t->Branch("str_trkIdx", &str_trkIdx);
1773       if (includeOnTrackHitData_) {
1774         t->Branch("str_onTrk_x", &str_onTrk_x);
1775         t->Branch("str_onTrk_y", &str_onTrk_y);
1776         t->Branch("str_onTrk_z", &str_onTrk_z);
1777         t->Branch("str_onTrk_xx", &str_onTrk_xx);
1778         t->Branch("str_onTrk_xy", &str_onTrk_xy);
1779         t->Branch("str_onTrk_yy", &str_onTrk_yy);
1780         t->Branch("str_onTrk_yz", &str_onTrk_yz);
1781         t->Branch("str_onTrk_zz", &str_onTrk_zz);
1782         t->Branch("str_onTrk_zx", &str_onTrk_zx);
1783       }
1784       if (includeTrackCandidates_)
1785         t->Branch("str_tcandIdx", &str_tcandIdx);
1786       if (includeSeeds_) {
1787         t->Branch("str_seeIdx", &str_seeIdx);
1788       }
1789       if (includeTrackingParticles_) {
1790         t->Branch("str_simHitIdx", &str_simHitIdx);
1791         if (simHitBySignificance_) {
1792           t->Branch("str_xySignificance", &str_xySignificance);
1793         }
1794         t->Branch("str_chargeFraction", &str_chargeFraction);
1795         t->Branch("str_simType", &str_simType);
1796       }
1797       t->Branch("str_x", &str_x);
1798       t->Branch("str_y", &str_y);
1799       t->Branch("str_z", &str_z);
1800       t->Branch("str_xx", &str_xx);
1801       t->Branch("str_xy", &str_xy);
1802       t->Branch("str_yy", &str_yy);
1803       t->Branch("str_yz", &str_yz);
1804       t->Branch("str_zz", &str_zz);
1805       t->Branch("str_zx", &str_zx);
1806       t->Branch("str_radL", &str_radL);
1807       t->Branch("str_bbxi", &str_bbxi);
1808       t->Branch("str_chargePerCM", &str_chargePerCM);
1809       t->Branch("str_clustSize", &str_clustSize);
1810       t->Branch("str_usedMask", &str_usedMask);
1811       //matched hits
1812       t->Branch("glu_isBarrel", &glu_isBarrel);
1813       glu_detId.book("glu", t);
1814       t->Branch("glu_monoIdx", &glu_monoIdx);
1815       t->Branch("glu_stereoIdx", &glu_stereoIdx);
1816       if (includeSeeds_) {
1817         t->Branch("glu_seeIdx", &glu_seeIdx);
1818       }
1819       t->Branch("glu_x", &glu_x);
1820       t->Branch("glu_y", &glu_y);
1821       t->Branch("glu_z", &glu_z);
1822       t->Branch("glu_xx", &glu_xx);
1823       t->Branch("glu_xy", &glu_xy);
1824       t->Branch("glu_yy", &glu_yy);
1825       t->Branch("glu_yz", &glu_yz);
1826       t->Branch("glu_zz", &glu_zz);
1827       t->Branch("glu_zx", &glu_zx);
1828       t->Branch("glu_radL", &glu_radL);
1829       t->Branch("glu_bbxi", &glu_bbxi);
1830       t->Branch("glu_chargePerCM", &glu_chargePerCM);
1831       t->Branch("glu_clustSizeMono", &glu_clustSizeMono);
1832       t->Branch("glu_clustSizeStereo", &glu_clustSizeStereo);
1833       t->Branch("glu_usedMaskMono", &glu_usedMaskMono);
1834       t->Branch("glu_usedMaskStereo", &glu_usedMaskStereo);
1835     }
1836     //phase2 OT
1837     if (includePhase2OTHits_) {
1838       t->Branch("ph2_isBarrel", &ph2_isBarrel);
1839       ph2_detId.book("ph2", t);
1840       t->Branch("ph2_trkIdx", &ph2_trkIdx);
1841       if (includeOnTrackHitData_) {
1842         t->Branch("ph2_onTrk_x", &ph2_onTrk_x);
1843         t->Branch("ph2_onTrk_y", &ph2_onTrk_y);
1844         t->Branch("ph2_onTrk_z", &ph2_onTrk_z);
1845         t->Branch("ph2_onTrk_xx", &ph2_onTrk_xx);
1846         t->Branch("ph2_onTrk_xy", &ph2_onTrk_xy);
1847         t->Branch("ph2_onTrk_yy", &ph2_onTrk_yy);
1848         t->Branch("ph2_onTrk_yz", &ph2_onTrk_yz);
1849         t->Branch("ph2_onTrk_zz", &ph2_onTrk_zz);
1850         t->Branch("ph2_onTrk_zx", &ph2_onTrk_zx);
1851       }
1852       if (includeTrackCandidates_)
1853         t->Branch("ph2_tcandIdx", &ph2_tcandIdx);
1854       if (includeSeeds_) {
1855         t->Branch("ph2_seeIdx", &ph2_seeIdx);
1856       }
1857       if (includeTrackingParticles_) {
1858         t->Branch("ph2_simHitIdx", &ph2_simHitIdx);
1859         if (simHitBySignificance_) {
1860           t->Branch("ph2_xySignificance", &ph2_xySignificance);
1861         }
1862         t->Branch("ph2_simType", &ph2_simType);
1863       }
1864       t->Branch("ph2_x", &ph2_x);
1865       t->Branch("ph2_y", &ph2_y);
1866       t->Branch("ph2_z", &ph2_z);
1867       t->Branch("ph2_xx", &ph2_xx);
1868       t->Branch("ph2_xy", &ph2_xy);
1869       t->Branch("ph2_yy", &ph2_yy);
1870       t->Branch("ph2_yz", &ph2_yz);
1871       t->Branch("ph2_zz", &ph2_zz);
1872       t->Branch("ph2_zx", &ph2_zx);
1873       t->Branch("ph2_radL", &ph2_radL);
1874       t->Branch("ph2_bbxi", &ph2_bbxi);
1875       t->Branch("ph2_bbxi", &ph2_bbxi);
1876     }
1877     //invalid hits
1878     t->Branch("inv_isBarrel", &inv_isBarrel);
1879     if (includeStripHits_)
1880       inv_detId.book("inv", t);
1881     else
1882       inv_detId_phase2.book("inv", t);
1883     t->Branch("inv_type", &inv_type);
1884     //simhits
1885     if (includeTrackingParticles_) {
1886       if (includeStripHits_)
1887         simhit_detId.book("simhit", t);
1888       else
1889         simhit_detId_phase2.book("simhit", t);
1890       t->Branch("simhit_x", &simhit_x);
1891       t->Branch("simhit_y", &simhit_y);
1892       t->Branch("simhit_z", &simhit_z);
1893       if (saveSimHitsP3_) {
1894         t->Branch("simhit_px", &simhit_px);
1895         t->Branch("simhit_py", &simhit_py);
1896         t->Branch("simhit_pz", &simhit_pz);
1897       }
1898       t->Branch("simhit_particle", &simhit_particle);
1899       t->Branch("simhit_process", &simhit_process);
1900       t->Branch("simhit_eloss", &simhit_eloss);
1901       t->Branch("simhit_tof", &simhit_tof);
1902       t->Branch("simhit_simTrkIdx", &simhit_simTrkIdx);
1903       t->Branch("simhit_hitIdx", &simhit_hitIdx);
1904       t->Branch("simhit_hitType", &simhit_hitType);
1905     }
1906   }
1907   //beam spot
1908   t->Branch("bsp_x", &bsp_x, "bsp_x/F");
1909   t->Branch("bsp_y", &bsp_y, "bsp_y/F");
1910   t->Branch("bsp_z", &bsp_z, "bsp_z/F");
1911   t->Branch("bsp_sigmax", &bsp_sigmax, "bsp_sigmax/F");
1912   t->Branch("bsp_sigmay", &bsp_sigmay, "bsp_sigmay/F");
1913   t->Branch("bsp_sigmaz", &bsp_sigmaz, "bsp_sigmaz/F");
1914   if (includeSeeds_) {
1915     //seeds
1916     t->Branch("see_fitok", &see_fitok);
1917     t->Branch("see_px", &see_px);
1918     t->Branch("see_py", &see_py);
1919     t->Branch("see_pz", &see_pz);
1920     t->Branch("see_pt", &see_pt);
1921     t->Branch("see_eta", &see_eta);
1922     t->Branch("see_phi", &see_phi);
1923     t->Branch("see_dxy", &see_dxy);
1924     t->Branch("see_dz", &see_dz);
1925     t->Branch("see_ptErr", &see_ptErr);
1926     t->Branch("see_etaErr", &see_etaErr);
1927     t->Branch("see_phiErr", &see_phiErr);
1928     t->Branch("see_dxyErr", &see_dxyErr);
1929     t->Branch("see_dzErr", &see_dzErr);
1930     t->Branch("see_chi2", &see_chi2);
1931     t->Branch("see_statePt", &see_statePt);
1932     t->Branch("see_stateTrajX", &see_stateTrajX);
1933     t->Branch("see_stateTrajY", &see_stateTrajY);
1934     t->Branch("see_stateTrajPx", &see_stateTrajPx);
1935     t->Branch("see_stateTrajPy", &see_stateTrajPy);
1936     t->Branch("see_stateTrajPz", &see_stateTrajPz);
1937     t->Branch("see_stateTrajGlbX", &see_stateTrajGlbX);
1938     t->Branch("see_stateTrajGlbY", &see_stateTrajGlbY);
1939     t->Branch("see_stateTrajGlbZ", &see_stateTrajGlbZ);
1940     t->Branch("see_stateTrajGlbPx", &see_stateTrajGlbPx);
1941     t->Branch("see_stateTrajGlbPy", &see_stateTrajGlbPy);
1942     t->Branch("see_stateTrajGlbPz", &see_stateTrajGlbPz);
1943     if (addSeedCurvCov_) {
1944       t->Branch("see_stateCurvCov", &see_stateCurvCov);
1945     }
1946     t->Branch("see_q", &see_q);
1947     t->Branch("see_nValid", &see_nValid);
1948     t->Branch("see_nPixel", &see_nPixel);
1949     t->Branch("see_nGlued", &see_nGlued);
1950     t->Branch("see_nStrip", &see_nStrip);
1951     t->Branch("see_nPhase2OT", &see_nPhase2OT);
1952     t->Branch("see_nCluster", &see_nCluster);
1953     t->Branch("see_algo", &see_algo);
1954     t->Branch("see_stopReason", &see_stopReason);
1955     t->Branch("see_nCands", &see_nCands);
1956     t->Branch("see_trkIdx", &see_trkIdx);
1957     if (includeTrackCandidates_)
1958       t->Branch("see_tcandIdx", &see_tcandIdx);
1959     if (includeTrackingParticles_) {
1960       t->Branch("see_simTrkIdx", &see_simTrkIdx);
1961       t->Branch("see_simTrkShareFrac", &see_simTrkShareFrac);
1962       t->Branch("see_bestSimTrkIdx", &see_bestSimTrkIdx);
1963       t->Branch("see_bestFromFirstHitSimTrkIdx", &see_bestFromFirstHitSimTrkIdx);
1964     } else {
1965       t->Branch("see_isTrue", &see_isTrue);
1966     }
1967     t->Branch("see_bestSimTrkShareFrac", &see_bestSimTrkShareFrac);
1968     t->Branch("see_bestFromFirstHitSimTrkShareFrac", &see_bestFromFirstHitSimTrkShareFrac);
1969     if (includeAllHits_) {
1970       t->Branch("see_hitIdx", &see_hitIdx);
1971       t->Branch("see_hitType", &see_hitType);
1972     }
1973     //seed algo offset
1974     t->Branch("see_offset", &see_offset);
1975   }
1976 
1977   //vertices
1978   t->Branch("vtx_x", &vtx_x);
1979   t->Branch("vtx_y", &vtx_y);
1980   t->Branch("vtx_z", &vtx_z);
1981   t->Branch("vtx_xErr", &vtx_xErr);
1982   t->Branch("vtx_yErr", &vtx_yErr);
1983   t->Branch("vtx_zErr", &vtx_zErr);
1984   t->Branch("vtx_ndof", &vtx_ndof);
1985   t->Branch("vtx_chi2", &vtx_chi2);
1986   t->Branch("vtx_fake", &vtx_fake);
1987   t->Branch("vtx_valid", &vtx_valid);
1988   t->Branch("vtx_trkIdx", &vtx_trkIdx);
1989 
1990   // tracking vertices
1991   t->Branch("simvtx_event", &simvtx_event);
1992   t->Branch("simvtx_bunchCrossing", &simvtx_bunchCrossing);
1993   t->Branch("simvtx_processType", &simvtx_processType);
1994   t->Branch("simvtx_x", &simvtx_x);
1995   t->Branch("simvtx_y", &simvtx_y);
1996   t->Branch("simvtx_z", &simvtx_z);
1997   t->Branch("simvtx_sourceSimIdx", &simvtx_sourceSimIdx);
1998   t->Branch("simvtx_daughterSimIdx", &simvtx_daughterSimIdx);
1999 
2000   t->Branch("simpv_idx", &simpv_idx);
2001 
2002   //t->Branch("" , &);
2003 }
2004 
2005 TrackingNtuple::~TrackingNtuple() {
2006   // do anything here that needs to be done at desctruction time
2007   // (e.g. close files, deallocate resources etc.)
2008 }
2009 
2010 //
2011 // member functions
2012 //
2013 void TrackingNtuple::clearVariables() {
2014   ev_run = 0;
2015   ev_lumi = 0;
2016   ev_event = 0;
2017 
2018   //tracks
2019   trk_px.clear();
2020   trk_py.clear();
2021   trk_pz.clear();
2022   trk_pt.clear();
2023   trk_inner_px.clear();
2024   trk_inner_py.clear();
2025   trk_inner_pz.clear();
2026   trk_inner_pt.clear();
2027   trk_outer_px.clear();
2028   trk_outer_py.clear();
2029   trk_outer_pz.clear();
2030   trk_outer_pt.clear();
2031   trk_eta.clear();
2032   trk_lambda.clear();
2033   trk_cotTheta.clear();
2034   trk_phi.clear();
2035   trk_dxy.clear();
2036   trk_dz.clear();
2037   trk_dxyPV.clear();
2038   trk_dzPV.clear();
2039   trk_dxyClosestPV.clear();
2040   trk_dzClosestPV.clear();
2041   trk_ptErr.clear();
2042   trk_etaErr.clear();
2043   trk_lambdaErr.clear();
2044   trk_phiErr.clear();
2045   trk_dxyErr.clear();
2046   trk_dzErr.clear();
2047   trk_refpoint_x.clear();
2048   trk_refpoint_y.clear();
2049   trk_refpoint_z.clear();
2050   trk_nChi2.clear();
2051   trk_nChi2_1Dmod.clear();
2052   trk_ndof.clear();
2053   for (auto& mva : trk_mvas) {
2054     mva.clear();
2055   }
2056   for (auto& mask : trk_qualityMasks) {
2057     mask.clear();
2058   }
2059   trk_q.clear();
2060   trk_nValid.clear();
2061   trk_nLost.clear();
2062   trk_nInactive.clear();
2063   trk_nPixel.clear();
2064   trk_nStrip.clear();
2065   trk_nOuterLost.clear();
2066   trk_nInnerLost.clear();
2067   trk_nOuterInactive.clear();
2068   trk_nInnerInactive.clear();
2069   trk_nPixelLay.clear();
2070   trk_nStripLay.clear();
2071   trk_n3DLay.clear();
2072   trk_nLostLay.clear();
2073   trk_nCluster.clear();
2074   trk_algo.clear();
2075   trk_originalAlgo.clear();
2076   trk_algoMask.clear();
2077   trk_stopReason.clear();
2078   trk_isHP.clear();
2079   trk_seedIdx.clear();
2080   trk_vtxIdx.clear();
2081   trk_isTrue.clear();
2082   trk_bestSimTrkIdx.clear();
2083   trk_bestSimTrkShareFrac.clear();
2084   trk_bestSimTrkShareFracSimDenom.clear();
2085   trk_bestSimTrkShareFracSimClusterDenom.clear();
2086   trk_bestSimTrkNChi2.clear();
2087   trk_bestFromFirstHitSimTrkIdx.clear();
2088   trk_bestFromFirstHitSimTrkShareFrac.clear();
2089   trk_bestFromFirstHitSimTrkShareFracSimDenom.clear();
2090   trk_bestFromFirstHitSimTrkShareFracSimClusterDenom.clear();
2091   trk_bestFromFirstHitSimTrkNChi2.clear();
2092   trk_simTrkIdx.clear();
2093   trk_simTrkShareFrac.clear();
2094   trk_simTrkNChi2.clear();
2095   trk_hitIdx.clear();
2096   trk_hitType.clear();
2097   //track candidates
2098   tcand_pca_valid.clear();
2099   tcand_pca_px.clear();
2100   tcand_pca_py.clear();
2101   tcand_pca_pz.clear();
2102   tcand_pca_pt.clear();
2103   tcand_pca_eta.clear();
2104   tcand_pca_phi.clear();
2105   tcand_pca_dxy.clear();
2106   tcand_pca_dz.clear();
2107   tcand_pca_ptErr.clear();
2108   tcand_pca_etaErr.clear();
2109   tcand_pca_lambdaErr.clear();
2110   tcand_pca_phiErr.clear();
2111   tcand_pca_dxyErr.clear();
2112   tcand_pca_dzErr.clear();
2113   tcand_px.clear();
2114   tcand_py.clear();
2115   tcand_pz.clear();
2116   tcand_pt.clear();
2117   tcand_x.clear();
2118   tcand_y.clear();
2119   tcand_z.clear();
2120   tcand_qbpErr.clear();
2121   tcand_lambdaErr.clear();
2122   tcand_phiErr.clear();
2123   tcand_xtErr.clear();
2124   tcand_ytErr.clear();
2125   tcand_ndof.clear();
2126   tcand_q.clear();
2127   tcand_nValid.clear();
2128   tcand_nPixel.clear();
2129   tcand_nStrip.clear();
2130   tcand_nCluster.clear();
2131   tcand_algo.clear();
2132   tcand_stopReason.clear();
2133   tcand_seedIdx.clear();
2134   tcand_vtxIdx.clear();
2135   tcand_isTrue.clear();
2136   tcand_bestSimTrkIdx.clear();
2137   tcand_bestSimTrkShareFrac.clear();
2138   tcand_bestSimTrkShareFracSimDenom.clear();
2139   tcand_bestSimTrkShareFracSimClusterDenom.clear();
2140   tcand_bestSimTrkNChi2.clear();
2141   tcand_bestFromFirstHitSimTrkIdx.clear();
2142   tcand_bestFromFirstHitSimTrkShareFrac.clear();
2143   tcand_bestFromFirstHitSimTrkShareFracSimDenom.clear();
2144   tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom.clear();
2145   tcand_bestFromFirstHitSimTrkNChi2.clear();
2146   tcand_simTrkShareFrac.clear();
2147   tcand_simTrkNChi2.clear();
2148   tcand_simTrkIdx.clear();
2149   tcand_hitIdx.clear();
2150   tcand_hitType.clear();
2151   //sim tracks
2152   sim_event.clear();
2153   sim_bunchCrossing.clear();
2154   sim_pdgId.clear();
2155   sim_genPdgIds.clear();
2156   sim_isFromBHadron.clear();
2157   sim_px.clear();
2158   sim_py.clear();
2159   sim_pz.clear();
2160   sim_pt.clear();
2161   sim_eta.clear();
2162   sim_phi.clear();
2163   sim_pca_pt.clear();
2164   sim_pca_eta.clear();
2165   sim_pca_lambda.clear();
2166   sim_pca_cotTheta.clear();
2167   sim_pca_phi.clear();
2168   sim_pca_dxy.clear();
2169   sim_pca_dz.clear();
2170   sim_q.clear();
2171   sim_nValid.clear();
2172   sim_nPixel.clear();
2173   sim_nStrip.clear();
2174   sim_nLay.clear();
2175   sim_nPixelLay.clear();
2176   sim_n3DLay.clear();
2177   sim_nTrackerHits.clear();
2178   sim_nRecoClusters.clear();
2179   sim_trkIdx.clear();
2180   sim_seedIdx.clear();
2181   sim_trkShareFrac.clear();
2182   sim_parentVtxIdx.clear();
2183   sim_decayVtxIdx.clear();
2184   sim_simHitIdx.clear();
2185   //pixels
2186   pix_isBarrel.clear();
2187   pix_detId.clear();
2188   pix_trkIdx.clear();
2189   pix_onTrk_x.clear();
2190   pix_onTrk_y.clear();
2191   pix_onTrk_z.clear();
2192   pix_onTrk_xx.clear();
2193   pix_onTrk_xy.clear();
2194   pix_onTrk_yy.clear();
2195   pix_onTrk_yz.clear();
2196   pix_onTrk_zz.clear();
2197   pix_onTrk_zx.clear();
2198   pix_tcandIdx.clear();
2199   pix_seeIdx.clear();
2200   pix_simHitIdx.clear();
2201   pix_xySignificance.clear();
2202   pix_chargeFraction.clear();
2203   pix_simType.clear();
2204   pix_x.clear();
2205   pix_y.clear();
2206   pix_z.clear();
2207   pix_xx.clear();
2208   pix_xy.clear();
2209   pix_yy.clear();
2210   pix_yz.clear();
2211   pix_zz.clear();
2212   pix_zx.clear();
2213   pix_radL.clear();
2214   pix_bbxi.clear();
2215   pix_clustSizeCol.clear();
2216   pix_clustSizeRow.clear();
2217   pix_usedMask.clear();
2218   //strips
2219   str_isBarrel.clear();
2220   str_detId.clear();
2221   str_trkIdx.clear();
2222   str_onTrk_x.clear();
2223   str_onTrk_y.clear();
2224   str_onTrk_z.clear();
2225   str_onTrk_xx.clear();
2226   str_onTrk_xy.clear();
2227   str_onTrk_yy.clear();
2228   str_onTrk_yz.clear();
2229   str_onTrk_zz.clear();
2230   str_onTrk_zx.clear();
2231   str_tcandIdx.clear();
2232   str_seeIdx.clear();
2233   str_simHitIdx.clear();
2234   str_xySignificance.clear();
2235   str_chargeFraction.clear();
2236   str_simType.clear();
2237   str_x.clear();
2238   str_y.clear();
2239   str_z.clear();
2240   str_xx.clear();
2241   str_xy.clear();
2242   str_yy.clear();
2243   str_yz.clear();
2244   str_zz.clear();
2245   str_zx.clear();
2246   str_radL.clear();
2247   str_bbxi.clear();
2248   str_chargePerCM.clear();
2249   str_clustSize.clear();
2250   str_usedMask.clear();
2251   //matched hits
2252   glu_isBarrel.clear();
2253   glu_detId.clear();
2254   glu_monoIdx.clear();
2255   glu_stereoIdx.clear();
2256   glu_seeIdx.clear();
2257   glu_x.clear();
2258   glu_y.clear();
2259   glu_z.clear();
2260   glu_xx.clear();
2261   glu_xy.clear();
2262   glu_yy.clear();
2263   glu_yz.clear();
2264   glu_zz.clear();
2265   glu_zx.clear();
2266   glu_radL.clear();
2267   glu_bbxi.clear();
2268   glu_chargePerCM.clear();
2269   glu_clustSizeMono.clear();
2270   glu_clustSizeStereo.clear();
2271   glu_usedMaskMono.clear();
2272   glu_usedMaskStereo.clear();
2273   //phase2 OT
2274   ph2_isBarrel.clear();
2275   ph2_detId.clear();
2276   ph2_trkIdx.clear();
2277   ph2_onTrk_x.clear();
2278   ph2_onTrk_y.clear();
2279   ph2_onTrk_z.clear();
2280   ph2_onTrk_xx.clear();
2281   ph2_onTrk_xy.clear();
2282   ph2_onTrk_yy.clear();
2283   ph2_onTrk_yz.clear();
2284   ph2_onTrk_zz.clear();
2285   ph2_onTrk_zx.clear();
2286   ph2_tcandIdx.clear();
2287   ph2_seeIdx.clear();
2288   ph2_xySignificance.clear();
2289   ph2_simHitIdx.clear();
2290   ph2_simType.clear();
2291   ph2_x.clear();
2292   ph2_y.clear();
2293   ph2_z.clear();
2294   ph2_xx.clear();
2295   ph2_xy.clear();
2296   ph2_yy.clear();
2297   ph2_yz.clear();
2298   ph2_zz.clear();
2299   ph2_zx.clear();
2300   ph2_radL.clear();
2301   ph2_bbxi.clear();
2302   //invalid hits
2303   inv_isBarrel.clear();
2304   inv_detId.clear();
2305   inv_detId_phase2.clear();
2306   inv_type.clear();
2307   // simhits
2308   simhit_detId.clear();
2309   simhit_detId_phase2.clear();
2310   simhit_x.clear();
2311   simhit_y.clear();
2312   simhit_z.clear();
2313   simhit_px.clear();
2314   simhit_py.clear();
2315   simhit_pz.clear();
2316   simhit_particle.clear();
2317   simhit_process.clear();
2318   simhit_eloss.clear();
2319   simhit_tof.clear();
2320   //simhit_simTrackId.clear();
2321   simhit_simTrkIdx.clear();
2322   simhit_hitIdx.clear();
2323   simhit_hitType.clear();
2324   //beamspot
2325   bsp_x = -9999.;
2326   bsp_y = -9999.;
2327   bsp_z = -9999.;
2328   bsp_sigmax = -9999.;
2329   bsp_sigmay = -9999.;
2330   bsp_sigmaz = -9999.;
2331   //seeds
2332   see_fitok.clear();
2333   see_px.clear();
2334   see_py.clear();
2335   see_pz.clear();
2336   see_pt.clear();
2337   see_eta.clear();
2338   see_phi.clear();
2339   see_dxy.clear();
2340   see_dz.clear();
2341   see_ptErr.clear();
2342   see_etaErr.clear();
2343   see_phiErr.clear();
2344   see_dxyErr.clear();
2345   see_dzErr.clear();
2346   see_chi2.clear();
2347   see_statePt.clear();
2348   see_stateTrajX.clear();
2349   see_stateTrajY.clear();
2350   see_stateTrajPx.clear();
2351   see_stateTrajPy.clear();
2352   see_stateTrajPz.clear();
2353   see_stateTrajGlbX.clear();
2354   see_stateTrajGlbY.clear();
2355   see_stateTrajGlbZ.clear();
2356   see_stateTrajGlbPx.clear();
2357   see_stateTrajGlbPy.clear();
2358   see_stateTrajGlbPz.clear();
2359   see_stateCurvCov.clear();
2360   see_q.clear();
2361   see_nValid.clear();
2362   see_nPixel.clear();
2363   see_nGlued.clear();
2364   see_nStrip.clear();
2365   see_nPhase2OT.clear();
2366   see_nCluster.clear();
2367   see_algo.clear();
2368   see_stopReason.clear();
2369   see_nCands.clear();
2370   see_trkIdx.clear();
2371   see_tcandIdx.clear();
2372   see_isTrue.clear();
2373   see_bestSimTrkIdx.clear();
2374   see_bestSimTrkShareFrac.clear();
2375   see_bestFromFirstHitSimTrkIdx.clear();
2376   see_bestFromFirstHitSimTrkShareFrac.clear();
2377   see_simTrkIdx.clear();
2378   see_simTrkShareFrac.clear();
2379   see_hitIdx.clear();
2380   see_hitType.clear();
2381   //seed algo offset
2382   see_offset.clear();
2383 
2384   // vertices
2385   vtx_x.clear();
2386   vtx_y.clear();
2387   vtx_z.clear();
2388   vtx_xErr.clear();
2389   vtx_yErr.clear();
2390   vtx_zErr.clear();
2391   vtx_ndof.clear();
2392   vtx_chi2.clear();
2393   vtx_fake.clear();
2394   vtx_valid.clear();
2395   vtx_trkIdx.clear();
2396 
2397   // Tracking vertices
2398   simvtx_event.clear();
2399   simvtx_bunchCrossing.clear();
2400   simvtx_processType.clear();
2401   simvtx_x.clear();
2402   simvtx_y.clear();
2403   simvtx_z.clear();
2404   simvtx_sourceSimIdx.clear();
2405   simvtx_daughterSimIdx.clear();
2406   simpv_idx.clear();
2407 }
2408 
2409 // ------------ method called for each event  ------------
2410 void TrackingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
2411   using namespace edm;
2412   using namespace reco;
2413   using namespace std;
2414 
2415   const auto& mf = iSetup.getData(mfToken_);
2416   const TrackerTopology& tTopo = iSetup.getData(tTopoToken_);
2417   const TrackerGeometry& tracker = iSetup.getData(tGeomToken_);
2418 
2419   edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator;
2420   iEvent.getByToken(trackAssociatorToken_, theAssociator);
2421   const reco::TrackToTrackingParticleAssociator& associatorByHits = *theAssociator;
2422 
2423   LogDebug("TrackingNtuple") << "Analyzing new event";
2424 
2425   //initialize tree variables
2426   clearVariables();
2427 
2428   // FIXME: we really need to move to edm::View for reading the
2429   // TrackingParticles... Unfortunately it has non-trivial
2430   // consequences on the associator/association interfaces etc.
2431   TrackingParticleRefVector tmpTP;
2432   const TrackingParticleRefVector* tmpTPptr = nullptr;
2433   edm::Handle<TrackingParticleCollection> TPCollectionH;
2434   edm::Handle<TrackingParticleRefVector> TPCollectionHRefVector;
2435 
2436   if (!trackingParticleToken_.isUninitialized()) {
2437     iEvent.getByToken(trackingParticleToken_, TPCollectionH);
2438     for (size_t i = 0, size = TPCollectionH->size(); i < size; ++i) {
2439       tmpTP.push_back(TrackingParticleRef(TPCollectionH, i));
2440     }
2441     tmpTPptr = &tmpTP;
2442   } else {
2443     iEvent.getByToken(trackingParticleRefToken_, TPCollectionHRefVector);
2444     tmpTPptr = TPCollectionHRefVector.product();
2445   }
2446   const TrackingParticleRefVector& tpCollection = *tmpTPptr;
2447 
2448   // Fill mapping from Ref::key() to index
2449   TrackingParticleRefKeyToIndex tpKeyToIndex;
2450   for (size_t i = 0; i < tpCollection.size(); ++i) {
2451     tpKeyToIndex[tpCollection[i].key()] = i;
2452   }
2453 
2454   // tracking vertices
2455   edm::Handle<TrackingVertexCollection> htv;
2456   iEvent.getByToken(trackingVertexToken_, htv);
2457   const TrackingVertexCollection& tvs = *htv;
2458 
2459   // Fill mapping from Ref::key() to index
2460   TrackingVertexRefVector tvRefs;
2461   TrackingVertexRefKeyToIndex tvKeyToIndex;
2462   for (size_t i = 0; i < tvs.size(); ++i) {
2463     const TrackingVertex& v = tvs[i];
2464     if (!includeOOT_ && v.eventId().bunchCrossing() != 0)  // Ignore OOTPU
2465       continue;
2466     tvKeyToIndex[i] = tvRefs.size();
2467     tvRefs.push_back(TrackingVertexRef(htv, i));
2468   }
2469 
2470   //get association maps, etc.
2471   Handle<ClusterTPAssociation> pCluster2TPListH;
2472   iEvent.getByToken(clusterTPMapToken_, pCluster2TPListH);
2473   const ClusterTPAssociation& clusterToTPMap = *pCluster2TPListH;
2474   edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
2475   iEvent.getByToken(simHitTPMapToken_, simHitsTPAssoc);
2476 
2477   // TP -> cluster count
2478   TrackingParticleRefKeyToCount tpKeyToClusterCount;
2479   for (const auto& clusterTP : clusterToTPMap) {
2480     tpKeyToClusterCount[clusterTP.second.key()] += 1;
2481   }
2482 
2483   // SimHit key -> index mapping
2484   SimHitRefKeyToIndex simHitRefKeyToIndex;
2485 
2486   //make a list to link TrackingParticles to its simhits
2487   std::vector<TPHitIndex> tpHitList;
2488 
2489   // Count the number of reco cluster per TP
2490 
2491   std::set<edm::ProductID> hitProductIds;
2492   std::map<edm::ProductID, size_t> seedCollToOffset;
2493 
2494   ev_run = iEvent.id().run();
2495   ev_lumi = iEvent.id().luminosityBlock();
2496   ev_event = iEvent.id().event();
2497 
2498   // Digi->Sim links for pixels and strips
2499   edm::Handle<edm::DetSetVector<PixelDigiSimLink>> pixelDigiSimLinksHandle;
2500   iEvent.getByToken(pixelSimLinkToken_, pixelDigiSimLinksHandle);
2501   const auto& pixelDigiSimLinks = *pixelDigiSimLinksHandle;
2502 
2503   edm::Handle<edm::DetSetVector<StripDigiSimLink>> stripDigiSimLinksHandle;
2504   iEvent.getByToken(stripSimLinkToken_, stripDigiSimLinksHandle);
2505 
2506   // Phase2 OT DigiSimLink
2507   edm::Handle<edm::DetSetVector<PixelDigiSimLink>> siphase2OTSimLinksHandle;
2508   iEvent.getByToken(siphase2OTSimLinksToken_, siphase2OTSimLinksHandle);
2509 
2510   //beamspot
2511   Handle<reco::BeamSpot> recoBeamSpotHandle;
2512   iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
2513   BeamSpot const& bs = *recoBeamSpotHandle;
2514   fillBeamSpot(bs);
2515 
2516   //prapare list to link matched hits to collection
2517   vector<pair<int, int>> monoStereoClusterList;
2518   if (includeAllHits_) {
2519     // simhits, only if TPs are saved as well
2520     if (includeTrackingParticles_) {
2521       fillSimHits(tracker, tpKeyToIndex, *simHitsTPAssoc, tTopo, simHitRefKeyToIndex, tpHitList);
2522     }
2523 
2524     //pixel hits
2525     fillPixelHits(iEvent,
2526                   tracker,
2527                   clusterToTPMap,
2528                   tpKeyToIndex,
2529                   *simHitsTPAssoc,
2530                   pixelDigiSimLinks,
2531                   tTopo,
2532                   simHitRefKeyToIndex,
2533                   hitProductIds);
2534 
2535     //strip hits
2536     if (includeStripHits_) {
2537       LogDebug("TrackingNtuple") << "foundStripSimLink";
2538       const auto& stripDigiSimLinks = *stripDigiSimLinksHandle;
2539       fillStripRphiStereoHits(iEvent,
2540                               tracker,
2541                               clusterToTPMap,
2542                               tpKeyToIndex,
2543                               *simHitsTPAssoc,
2544                               stripDigiSimLinks,
2545                               tTopo,
2546                               simHitRefKeyToIndex,
2547                               hitProductIds);
2548 
2549       //matched hits
2550       fillStripMatchedHits(iEvent, tracker, tTopo, monoStereoClusterList);
2551     }
2552 
2553     if (includePhase2OTHits_) {
2554       LogDebug("TrackingNtuple") << "foundPhase2OTSimLinks";
2555       const auto& phase2OTSimLinks = *siphase2OTSimLinksHandle;
2556       fillPhase2OTHits(iEvent,
2557                        clusterToTPMap,
2558                        tracker,
2559                        tpKeyToIndex,
2560                        *simHitsTPAssoc,
2561                        phase2OTSimLinks,
2562                        tTopo,
2563                        simHitRefKeyToIndex,
2564                        hitProductIds);
2565     }
2566   }
2567 
2568   //seeds
2569   if (includeSeeds_) {
2570     fillSeeds(iEvent,
2571               tpCollection,
2572               tpKeyToIndex,
2573               bs,
2574               tracker,
2575               associatorByHits,
2576               clusterToTPMap,
2577               mf,
2578               tTopo,
2579               monoStereoClusterList,
2580               hitProductIds,
2581               seedCollToOffset);
2582   }
2583 
2584   //tracks
2585   edm::Handle<edm::View<reco::Track>> tracksHandle;
2586   iEvent.getByToken(trackToken_, tracksHandle);
2587   const edm::View<reco::Track>& tracks = *tracksHandle;
2588   // The associator interfaces really need to be fixed...
2589   edm::RefToBaseVector<reco::Track> trackRefs;
2590   for (edm::View<Track>::size_type i = 0; i < tracks.size(); ++i) {
2591     trackRefs.push_back(tracks.refAt(i));
2592   }
2593   std::vector<const MVACollection*> mvaColls;
2594   std::vector<const QualityMaskCollection*> qualColls;
2595   if (includeMVA_) {
2596     edm::Handle<MVACollection> hmva;
2597     edm::Handle<QualityMaskCollection> hqual;
2598 
2599     for (const auto& tokenTpl : mvaQualityCollectionTokens_) {
2600       iEvent.getByToken(std::get<0>(tokenTpl), hmva);
2601       iEvent.getByToken(std::get<1>(tokenTpl), hqual);
2602 
2603       mvaColls.push_back(hmva.product());
2604       qualColls.push_back(hqual.product());
2605       if (mvaColls.back()->size() != tracks.size()) {
2606         throw cms::Exception("Configuration")
2607             << "Inconsistency in track collection and MVA sizes. Track collection has " << tracks.size()
2608             << " tracks, whereas the MVA " << (mvaColls.size() - 1) << " has " << mvaColls.back()->size()
2609             << " entries. Double-check your configuration.";
2610       }
2611       if (qualColls.back()->size() != tracks.size()) {
2612         throw cms::Exception("Configuration")
2613             << "Inconsistency in track collection and quality mask sizes. Track collection has " << tracks.size()
2614             << " tracks, whereas the quality mask " << (qualColls.size() - 1) << " has " << qualColls.back()->size()
2615             << " entries. Double-check your configuration.";
2616       }
2617     }
2618   }
2619 
2620   edm::Handle<reco::VertexCollection> vertices;
2621   iEvent.getByToken(vertexToken_, vertices);
2622 
2623   fillTracks(trackRefs,
2624              tracker,
2625              tpCollection,
2626              tpKeyToIndex,
2627              tpKeyToClusterCount,
2628              mf,
2629              bs,
2630              *vertices,
2631              associatorByHits,
2632              clusterToTPMap,
2633              tTopo,
2634              hitProductIds,
2635              seedCollToOffset,
2636              mvaColls,
2637              qualColls);
2638 
2639   if (includeTrackCandidates_) {
2640     //candidates
2641     for (auto const& token : candidateTokens_) {
2642       auto const tCandsHandle = iEvent.getHandle(token);
2643 
2644       edm::EDConsumerBase::Labels labels;
2645       labelsForToken(token, labels);
2646       TString label = labels.module;
2647       label.ReplaceAll("TrackCandidates", "");
2648       label.ReplaceAll("muonSeeded", "muonSeededStep");
2649       int algo = reco::TrackBase::algoByName(label.Data());
2650 
2651       fillCandidates(tCandsHandle,
2652                      algo,
2653                      tpCollection,
2654                      tpKeyToIndex,
2655                      tpKeyToClusterCount,
2656                      mf,
2657                      bs,
2658                      *vertices,
2659                      associatorByHits,
2660                      clusterToTPMap,
2661                      tracker,
2662                      tTopo,
2663                      hitProductIds,
2664                      seedCollToOffset);
2665     }
2666   }
2667 
2668   //tracking particles
2669   //sort association maps with simHits
2670   std::sort(tpHitList.begin(), tpHitList.end(), tpHitIndexListLessSort);
2671   fillTrackingParticles(
2672       iEvent, iSetup, trackRefs, bs, tpCollection, tvKeyToIndex, associatorByHits, tpHitList, tpKeyToClusterCount);
2673 
2674   // vertices
2675   fillVertices(*vertices, trackRefs);
2676 
2677   // tracking vertices
2678   fillTrackingVertices(tvRefs, tpKeyToIndex);
2679 
2680   t->Fill();
2681 }
2682 
2683 void TrackingNtuple::fillBeamSpot(const reco::BeamSpot& bs) {
2684   bsp_x = bs.x0();
2685   bsp_y = bs.y0();
2686   bsp_z = bs.x0();
2687   bsp_sigmax = bs.BeamWidthX();
2688   bsp_sigmay = bs.BeamWidthY();
2689   bsp_sigmaz = bs.sigmaZ();
2690 }
2691 
2692 namespace {
2693   template <typename SimLink>
2694   struct GetCluster;
2695   template <>
2696   struct GetCluster<PixelDigiSimLink> {
2697     static const SiPixelCluster& call(const OmniClusterRef& cluster) { return cluster.pixelCluster(); }
2698   };
2699   template <>
2700   struct GetCluster<StripDigiSimLink> {
2701     static const SiStripCluster& call(const OmniClusterRef& cluster) { return cluster.stripCluster(); }
2702   };
2703 }  // namespace
2704 
2705 template <typename SimLink>
2706 TrackingNtuple::SimHitData TrackingNtuple::matchCluster(
2707     const OmniClusterRef& cluster,
2708     DetId hitId,
2709     int clusterKey,
2710     const TrackingRecHit& hit,
2711     const ClusterTPAssociation& clusterToTPMap,
2712     const TrackingParticleRefKeyToIndex& tpKeyToIndex,
2713     const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
2714     const edm::DetSetVector<SimLink>& digiSimLinks,
2715     const SimHitRefKeyToIndex& simHitRefKeyToIndex,
2716     HitType hitType) {
2717   SimHitData ret;
2718 
2719   std::map<unsigned int, double> simTrackIdToChargeFraction;
2720   if (hitType == HitType::Phase2OT)
2721     simTrackIdToChargeFraction = chargeFraction(cluster.phase2OTCluster(), hitId, digiSimLinks);
2722   else
2723     simTrackIdToChargeFraction = chargeFraction(GetCluster<SimLink>::call(cluster), hitId, digiSimLinks);
2724 
2725   float h_x = 0, h_y = 0;
2726   float h_xx = 0, h_xy = 0, h_yy = 0;
2727   if (simHitBySignificance_) {
2728     h_x = hit.localPosition().x();
2729     h_y = hit.localPosition().y();
2730     h_xx = hit.localPositionError().xx();
2731     h_xy = hit.localPositionError().xy();
2732     h_yy = hit.localPositionError().yy();
2733   }
2734 
2735   ret.type = HitSimType::Noise;
2736   auto range = clusterToTPMap.equal_range(cluster);
2737   if (range.first != range.second) {
2738     for (auto ip = range.first; ip != range.second; ++ip) {
2739       const TrackingParticleRef& trackingParticle = ip->second;
2740 
2741       // Find out if the cluster is from signal/ITPU/OOTPU
2742       const auto event = trackingParticle->eventId().event();
2743       const auto bx = trackingParticle->eventId().bunchCrossing();
2744       HitSimType type = HitSimType::OOTPileup;
2745       if (bx == 0) {
2746         type = (event == 0 ? HitSimType::Signal : HitSimType::ITPileup);
2747       } else
2748         type = HitSimType::OOTPileup;
2749       ret.type = static_cast<HitSimType>(std::min(static_cast<int>(ret.type), static_cast<int>(type)));
2750 
2751       // Limit to only input TrackingParticles (usually signal+ITPU)
2752       auto tpIndex = tpKeyToIndex.find(trackingParticle.key());
2753       if (tpIndex == tpKeyToIndex.end())
2754         continue;
2755 
2756       //now get the corresponding sim hit
2757       std::pair<TrackingParticleRef, TrackPSimHitRef> simHitTPpairWithDummyTP(trackingParticle, TrackPSimHitRef());
2758       //SimHit is dummy: for simHitTPAssociationListGreater sorting only the TP is needed
2759       auto range = std::equal_range(simHitsTPAssoc.begin(),
2760                                     simHitsTPAssoc.end(),
2761                                     simHitTPpairWithDummyTP,
2762                                     SimHitTPAssociationProducer::simHitTPAssociationListGreater);
2763       bool foundSimHit = false;
2764       bool foundElectron = false;
2765       int foundElectrons = 0;
2766       int foundNonElectrons = 0;
2767       for (auto ip = range.first; ip != range.second; ++ip) {
2768         TrackPSimHitRef TPhit = ip->second;
2769         DetId dId = DetId(TPhit->detUnitId());
2770         if (dId.rawId() == hitId.rawId()) {
2771           // skip electron SimHits for non-electron TPs also here
2772           if (std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
2773             foundElectrons++;
2774           } else {
2775             foundNonElectrons++;
2776           }
2777         }
2778       }
2779 
2780       float minSignificance = 1e12;
2781       if (simHitBySignificance_) {  //save the best matching hit
2782 
2783         int simHitKey = -1;
2784         edm::ProductID simHitID;
2785         for (auto ip = range.first; ip != range.second; ++ip) {
2786           TrackPSimHitRef TPhit = ip->second;
2787           DetId dId = DetId(TPhit->detUnitId());
2788           if (dId.rawId() == hitId.rawId()) {
2789             // skip electron SimHits for non-electron TPs also here
2790             if (std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
2791               foundElectron = true;
2792               if (!keepEleSimHits_)
2793                 continue;
2794             }
2795 
2796             float sx = TPhit->localPosition().x();
2797             float sy = TPhit->localPosition().y();
2798             float dx = sx - h_x;
2799             float dy = sy - h_y;
2800             float sig = (dx * dx * h_yy - 2 * dx * dy * h_xy + dy * dy * h_xx) / (h_xx * h_yy - h_xy * h_xy);
2801 
2802             if (sig < minSignificance) {
2803               minSignificance = sig;
2804               foundSimHit = true;
2805               simHitKey = TPhit.key();
2806               simHitID = TPhit.id();
2807             }
2808           }
2809         }  //loop over matching hits
2810 
2811         auto simHitIndex = simHitRefKeyToIndex.at(std::make_pair(simHitKey, simHitID));
2812         ret.matchingSimHit.push_back(simHitIndex);
2813 
2814         double chargeFraction = 0.;
2815         for (const SimTrack& simtrk : trackingParticle->g4Tracks()) {
2816           auto found = simTrackIdToChargeFraction.find(simtrk.trackId());
2817           if (found != simTrackIdToChargeFraction.end()) {
2818             chargeFraction += found->second;
2819           }
2820         }
2821         ret.xySignificance.push_back(minSignificance);
2822         ret.chargeFraction.push_back(chargeFraction);
2823 
2824         // only for debug prints
2825         ret.bunchCrossing.push_back(bx);
2826         ret.event.push_back(event);
2827 
2828         simhit_hitIdx[simHitIndex].push_back(clusterKey);
2829         simhit_hitType[simHitIndex].push_back(static_cast<int>(hitType));
2830 
2831       } else {  //save all matching hits
2832         for (auto ip = range.first; ip != range.second; ++ip) {
2833           TrackPSimHitRef TPhit = ip->second;
2834           DetId dId = DetId(TPhit->detUnitId());
2835           if (dId.rawId() == hitId.rawId()) {
2836             // skip electron SimHits for non-electron TPs also here
2837             if (std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
2838               foundElectron = true;
2839               if (!keepEleSimHits_)
2840                 continue;
2841               if (foundNonElectrons > 0)
2842                 continue;  //prioritize: skip electrons if non-electrons are present
2843             }
2844 
2845             foundSimHit = true;
2846             auto simHitKey = TPhit.key();
2847             auto simHitID = TPhit.id();
2848 
2849             auto simHitIndex = simHitRefKeyToIndex.at(std::make_pair(simHitKey, simHitID));
2850             ret.matchingSimHit.push_back(simHitIndex);
2851 
2852             double chargeFraction = 0.;
2853             for (const SimTrack& simtrk : trackingParticle->g4Tracks()) {
2854               auto found = simTrackIdToChargeFraction.find(simtrk.trackId());
2855               if (found != simTrackIdToChargeFraction.end()) {
2856                 chargeFraction += found->second;
2857               }
2858             }
2859             ret.xySignificance.push_back(minSignificance);
2860             ret.chargeFraction.push_back(chargeFraction);
2861 
2862             // only for debug prints
2863             ret.bunchCrossing.push_back(bx);
2864             ret.event.push_back(event);
2865 
2866             simhit_hitIdx[simHitIndex].push_back(clusterKey);
2867             simhit_hitType[simHitIndex].push_back(static_cast<int>(hitType));
2868           }
2869         }
2870       }  //if/else simHitBySignificance_
2871       if (!foundSimHit) {
2872         // In case we didn't find a simhit because of filtered-out
2873         // electron SimHit, just ignore the missing SimHit.
2874         if (foundElectron && !keepEleSimHits_)
2875           continue;
2876 
2877         auto ex = cms::Exception("LogicError")
2878                   << "Did not find SimHit for reco hit DetId " << hitId.rawId() << " for TP " << trackingParticle.key()
2879                   << " bx:event " << bx << ":" << event << " PDGid " << trackingParticle->pdgId() << " q "
2880                   << trackingParticle->charge() << " p4 " << trackingParticle->p4() << " nG4 "
2881                   << trackingParticle->g4Tracks().size() << ".\nFound SimHits from detectors ";
2882         for (auto ip = range.first; ip != range.second; ++ip) {
2883           TrackPSimHitRef TPhit = ip->second;
2884           DetId dId = DetId(TPhit->detUnitId());
2885           ex << dId.rawId() << " ";
2886         }
2887         if (trackingParticle->eventId().event() != 0) {
2888           ex << "\nSince this is a TrackingParticle from pileup, check that you're running the pileup mixing in "
2889                 "playback mode.";
2890         }
2891         throw ex;
2892       }
2893     }
2894   }
2895 
2896   return ret;
2897 }
2898 
2899 void TrackingNtuple::fillSimHits(const TrackerGeometry& tracker,
2900                                  const TrackingParticleRefKeyToIndex& tpKeyToIndex,
2901                                  const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
2902                                  const TrackerTopology& tTopo,
2903                                  SimHitRefKeyToIndex& simHitRefKeyToIndex,
2904                                  std::vector<TPHitIndex>& tpHitList) {
2905   for (const auto& assoc : simHitsTPAssoc) {
2906     auto tpKey = assoc.first.key();
2907 
2908     // SimHitTPAssociationList can contain more TrackingParticles than
2909     // what are given to this EDAnalyzer, so we can filter those out here.
2910     auto found = tpKeyToIndex.find(tpKey);
2911     if (found == tpKeyToIndex.end())
2912       continue;
2913     const auto tpIndex = found->second;
2914 
2915     // skip non-tracker simhits (mostly muons)
2916     const auto& simhit = *(assoc.second);
2917     auto detId = DetId(simhit.detUnitId());
2918     if (detId.det() != DetId::Tracker)
2919       continue;
2920 
2921     // Skip electron SimHits for non-electron TrackingParticles to
2922     // filter out delta rays. The delta ray hits just confuse. If we
2923     // need them later, let's add them as a separate "collection" of
2924     // hits of a TP
2925     const TrackingParticle& tp = *(assoc.first);
2926     if (!keepEleSimHits_ && std::abs(simhit.particleType()) == 11 && std::abs(tp.pdgId()) != 11)
2927       continue;
2928 
2929     auto simHitKey = std::make_pair(assoc.second.key(), assoc.second.id());
2930 
2931     if (simHitRefKeyToIndex.find(simHitKey) != simHitRefKeyToIndex.end()) {
2932       for (const auto& assoc2 : simHitsTPAssoc) {
2933         if (std::make_pair(assoc2.second.key(), assoc2.second.id()) == simHitKey) {
2934 #ifdef EDM_ML_DEBUG
2935           auto range1 = std::equal_range(simHitsTPAssoc.begin(),
2936                                          simHitsTPAssoc.end(),
2937                                          std::make_pair(assoc.first, TrackPSimHitRef()),
2938                                          SimHitTPAssociationProducer::simHitTPAssociationListGreater);
2939           auto range2 = std::equal_range(simHitsTPAssoc.begin(),
2940                                          simHitsTPAssoc.end(),
2941                                          std::make_pair(assoc2.first, TrackPSimHitRef()),
2942                                          SimHitTPAssociationProducer::simHitTPAssociationListGreater);
2943 
2944           LogTrace("TrackingNtuple") << "Earlier TP " << assoc2.first.key() << " SimTrack Ids";
2945           for (const auto& simTrack : assoc2.first->g4Tracks()) {
2946             edm::LogPrint("TrackingNtuple") << " SimTrack " << simTrack.trackId() << " BX:event "
2947                                             << simTrack.eventId().bunchCrossing() << ":" << simTrack.eventId().event();
2948           }
2949           for (auto iHit = range2.first; iHit != range2.second; ++iHit) {
2950             LogTrace("TrackingNtuple") << " SimHit " << iHit->second.key() << " " << iHit->second.id() << " tof "
2951                                        << iHit->second->tof() << " trackId " << iHit->second->trackId() << " BX:event "
2952                                        << iHit->second->eventId().bunchCrossing() << ":"
2953                                        << iHit->second->eventId().event();
2954           }
2955           LogTrace("TrackingNtuple") << "Current TP " << assoc.first.key() << " SimTrack Ids";
2956           for (const auto& simTrack : assoc.first->g4Tracks()) {
2957             edm::LogPrint("TrackingNtuple") << " SimTrack " << simTrack.trackId() << " BX:event "
2958                                             << simTrack.eventId().bunchCrossing() << ":" << simTrack.eventId().event();
2959           }
2960           for (auto iHit = range1.first; iHit != range1.second; ++iHit) {
2961             LogTrace("TrackingNtuple") << " SimHit " << iHit->second.key() << " " << iHit->second.id() << " tof "
2962                                        << iHit->second->tof() << " trackId " << iHit->second->trackId() << " BX:event "
2963                                        << iHit->second->eventId().bunchCrossing() << ":"
2964                                        << iHit->second->eventId().event();
2965           }
2966 #endif
2967 
2968           throw cms::Exception("LogicError")
2969               << "Got second time the SimHit " << simHitKey.first << " of " << simHitKey.second
2970               << ", first time with TrackingParticle " << assoc2.first.key() << ", now with " << tpKey;
2971         }
2972       }
2973       throw cms::Exception("LogicError") << "Got second time the SimHit " << simHitKey.first << " of "
2974                                          << simHitKey.second << ", now with TrackingParticle " << tpKey
2975                                          << ", but I didn't find the first occurrance!";
2976     }
2977 
2978     auto det = tracker.idToDetUnit(detId);
2979     if (!det)
2980       throw cms::Exception("LogicError") << "Did not find a det unit for DetId " << simhit.detUnitId()
2981                                          << " from tracker geometry";
2982 
2983     const auto pos = det->surface().toGlobal(simhit.localPosition());
2984     const float tof = simhit.timeOfFlight();
2985 
2986     const auto simHitIndex = simhit_x.size();
2987     simHitRefKeyToIndex[simHitKey] = simHitIndex;
2988 
2989     if (includeStripHits_)
2990       simhit_detId.push_back(tracker, tTopo, detId);
2991     else
2992       simhit_detId_phase2.push_back(tracker, tTopo, detId);
2993     simhit_x.push_back(pos.x());
2994     simhit_y.push_back(pos.y());
2995     simhit_z.push_back(pos.z());
2996     if (saveSimHitsP3_) {
2997       const auto mom = det->surface().toGlobal(simhit.momentumAtEntry());
2998       simhit_px.push_back(mom.x());
2999       simhit_py.push_back(mom.y());
3000       simhit_pz.push_back(mom.z());
3001     }
3002     simhit_particle.push_back(simhit.particleType());
3003     simhit_process.push_back(simhit.processType());
3004     simhit_eloss.push_back(simhit.energyLoss());
3005     simhit_tof.push_back(tof);
3006     //simhit_simTrackId.push_back(simhit.trackId());
3007 
3008     simhit_simTrkIdx.push_back(tpIndex);
3009 
3010     simhit_hitIdx.emplace_back();   // filled in matchCluster
3011     simhit_hitType.emplace_back();  // filled in matchCluster
3012 
3013     tpHitList.emplace_back(tpKey, simHitIndex, tof, simhit.detUnitId());
3014   }
3015 }
3016 
3017 void TrackingNtuple::fillPixelHits(const edm::Event& iEvent,
3018                                    const TrackerGeometry& tracker,
3019                                    const ClusterTPAssociation& clusterToTPMap,
3020                                    const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3021                                    const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
3022                                    const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
3023                                    const TrackerTopology& tTopo,
3024                                    const SimHitRefKeyToIndex& simHitRefKeyToIndex,
3025                                    std::set<edm::ProductID>& hitProductIds) {
3026   std::vector<std::pair<uint64_t, PixelMaskContainer const*>> pixelMasks;
3027   pixelMasks.reserve(pixelUseMaskTokens_.size());
3028   for (const auto& itoken : pixelUseMaskTokens_) {
3029     edm::Handle<PixelMaskContainer> aH;
3030     iEvent.getByToken(itoken.second, aH);
3031     pixelMasks.emplace_back(1 << itoken.first, aH.product());
3032   }
3033   auto pixUsedMask = [&pixelMasks](size_t key) {
3034     uint64_t mask = 0;
3035     for (auto const& m : pixelMasks) {
3036       if (m.second->mask(key))
3037         mask |= m.first;
3038     }
3039     return mask;
3040   };
3041 
3042   edm::Handle<SiPixelRecHitCollection> pixelHits;
3043   iEvent.getByToken(pixelRecHitToken_, pixelHits);
3044   for (auto it = pixelHits->begin(); it != pixelHits->end(); it++) {
3045     const DetId hitId = it->detId();
3046     for (auto hit = it->begin(); hit != it->end(); hit++) {
3047       hitProductIds.insert(hit->cluster().id());
3048 
3049       const int key = hit->cluster().key();
3050       const int lay = tTopo.layer(hitId);
3051 
3052       pix_isBarrel.push_back(hitId.subdetId() == 1);
3053       pix_detId.push_back(tracker, tTopo, hitId);
3054       pix_trkIdx.emplace_back();  // filled in fillTracks
3055       pix_onTrk_x.emplace_back();
3056       pix_onTrk_y.emplace_back();
3057       pix_onTrk_z.emplace_back();
3058       pix_onTrk_xx.emplace_back();
3059       pix_onTrk_xy.emplace_back();
3060       pix_onTrk_yy.emplace_back();
3061       pix_onTrk_yz.emplace_back();
3062       pix_onTrk_zz.emplace_back();
3063       pix_onTrk_zx.emplace_back();
3064       pix_tcandIdx.emplace_back();  // filled in fillCandidates
3065       pix_seeIdx.emplace_back();    // filled in fillSeeds
3066       pix_x.push_back(hit->globalPosition().x());
3067       pix_y.push_back(hit->globalPosition().y());
3068       pix_z.push_back(hit->globalPosition().z());
3069       pix_xx.push_back(hit->globalPositionError().cxx());
3070       pix_xy.push_back(hit->globalPositionError().cyx());
3071       pix_yy.push_back(hit->globalPositionError().cyy());
3072       pix_yz.push_back(hit->globalPositionError().czy());
3073       pix_zz.push_back(hit->globalPositionError().czz());
3074       pix_zx.push_back(hit->globalPositionError().czx());
3075       pix_radL.push_back(hit->surface()->mediumProperties().radLen());
3076       pix_bbxi.push_back(hit->surface()->mediumProperties().xi());
3077       pix_clustSizeCol.push_back(hit->cluster()->sizeY());
3078       pix_clustSizeRow.push_back(hit->cluster()->sizeX());
3079       pix_usedMask.push_back(pixUsedMask(hit->firstClusterRef().key()));
3080 
3081       LogTrace("TrackingNtuple") << "pixHit cluster=" << key << " subdId=" << hitId.subdetId() << " lay=" << lay
3082                                  << " rawId=" << hitId.rawId() << " pos =" << hit->globalPosition();
3083       if (includeTrackingParticles_) {
3084         SimHitData simHitData = matchCluster(hit->firstClusterRef(),
3085                                              hitId,
3086                                              key,
3087                                              *hit,
3088                                              clusterToTPMap,
3089                                              tpKeyToIndex,
3090                                              simHitsTPAssoc,
3091                                              digiSimLink,
3092                                              simHitRefKeyToIndex,
3093                                              HitType::Pixel);
3094         pix_simHitIdx.push_back(simHitData.matchingSimHit);
3095         pix_simType.push_back(static_cast<int>(simHitData.type));
3096         pix_xySignificance.push_back(simHitData.xySignificance);
3097         pix_chargeFraction.push_back(simHitData.chargeFraction);
3098         LogTrace("TrackingNtuple") << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
3099         if (!simHitData.matchingSimHit.empty()) {
3100           const auto simHitIdx = simHitData.matchingSimHit[0];
3101           LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx << " simHitPos="
3102                                      << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3103                                      << " energyLoss=" << simhit_eloss[simHitIdx]
3104                                      << " particleType=" << simhit_particle[simHitIdx]
3105                                      << " processType=" << simhit_process[simHitIdx]
3106                                      << " bunchCrossing=" << simHitData.bunchCrossing[0]
3107                                      << " event=" << simHitData.event[0];
3108         }
3109       }
3110     }
3111   }
3112 }
3113 
3114 void TrackingNtuple::fillStripRphiStereoHits(const edm::Event& iEvent,
3115                                              const TrackerGeometry& tracker,
3116                                              const ClusterTPAssociation& clusterToTPMap,
3117                                              const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3118                                              const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
3119                                              const edm::DetSetVector<StripDigiSimLink>& digiSimLink,
3120                                              const TrackerTopology& tTopo,
3121                                              const SimHitRefKeyToIndex& simHitRefKeyToIndex,
3122                                              std::set<edm::ProductID>& hitProductIds) {
3123   std::vector<std::pair<uint64_t, StripMaskContainer const*>> stripMasks;
3124   stripMasks.reserve(stripUseMaskTokens_.size());
3125   for (const auto& itoken : stripUseMaskTokens_) {
3126     edm::Handle<StripMaskContainer> aH;
3127     iEvent.getByToken(itoken.second, aH);
3128     stripMasks.emplace_back(1 << itoken.first, aH.product());
3129   }
3130   auto strUsedMask = [&stripMasks](size_t key) {
3131     uint64_t mask = 0;
3132     for (auto const& m : stripMasks) {
3133       if (m.second->mask(key))
3134         mask |= m.first;
3135     }
3136     return mask;
3137   };
3138 
3139   //index strip hit branches by cluster index
3140   edm::Handle<SiStripRecHit2DCollection> rphiHits;
3141   iEvent.getByToken(stripRphiRecHitToken_, rphiHits);
3142   edm::Handle<SiStripRecHit2DCollection> stereoHits;
3143   iEvent.getByToken(stripStereoRecHitToken_, stereoHits);
3144   int totalStripHits = rphiHits->dataSize() + stereoHits->dataSize();
3145   str_isBarrel.resize(totalStripHits);
3146   str_detId.resize(totalStripHits);
3147   str_trkIdx.resize(totalStripHits);  // filled in fillTracks
3148   str_onTrk_x.resize(totalStripHits);
3149   str_onTrk_y.resize(totalStripHits);
3150   str_onTrk_z.resize(totalStripHits);
3151   str_onTrk_xx.resize(totalStripHits);
3152   str_onTrk_xy.resize(totalStripHits);
3153   str_onTrk_yy.resize(totalStripHits);
3154   str_onTrk_yz.resize(totalStripHits);
3155   str_onTrk_zz.resize(totalStripHits);
3156   str_onTrk_zx.resize(totalStripHits);
3157   str_tcandIdx.resize(totalStripHits);  // filled in fillCandidates
3158   str_seeIdx.resize(totalStripHits);    // filled in fillSeeds
3159   str_simHitIdx.resize(totalStripHits);
3160   str_simType.resize(totalStripHits);
3161   str_chargeFraction.resize(totalStripHits);
3162   str_x.resize(totalStripHits);
3163   str_y.resize(totalStripHits);
3164   str_z.resize(totalStripHits);
3165   str_xx.resize(totalStripHits);
3166   str_xy.resize(totalStripHits);
3167   str_yy.resize(totalStripHits);
3168   str_yz.resize(totalStripHits);
3169   str_zz.resize(totalStripHits);
3170   str_zx.resize(totalStripHits);
3171   str_xySignificance.resize(totalStripHits);
3172   str_chargeFraction.resize(totalStripHits);
3173   str_radL.resize(totalStripHits);
3174   str_bbxi.resize(totalStripHits);
3175   str_chargePerCM.resize(totalStripHits);
3176   str_clustSize.resize(totalStripHits);
3177   str_usedMask.resize(totalStripHits);
3178 
3179   auto fill = [&](const SiStripRecHit2DCollection& hits, const char* name) {
3180     for (const auto& detset : hits) {
3181       const DetId hitId = detset.detId();
3182       for (const auto& hit : detset) {
3183         hitProductIds.insert(hit.cluster().id());
3184 
3185         const int key = hit.cluster().key();
3186         const int lay = tTopo.layer(hitId);
3187         str_isBarrel[key] = (hitId.subdetId() == StripSubdetector::TIB || hitId.subdetId() == StripSubdetector::TOB);
3188         str_detId.set(key, tracker, tTopo, hitId);
3189         str_x[key] = hit.globalPosition().x();
3190         str_y[key] = hit.globalPosition().y();
3191         str_z[key] = hit.globalPosition().z();
3192         str_xx[key] = hit.globalPositionError().cxx();
3193         str_xy[key] = hit.globalPositionError().cyx();
3194         str_yy[key] = hit.globalPositionError().cyy();
3195         str_yz[key] = hit.globalPositionError().czy();
3196         str_zz[key] = hit.globalPositionError().czz();
3197         str_zx[key] = hit.globalPositionError().czx();
3198         str_radL[key] = hit.surface()->mediumProperties().radLen();
3199         str_bbxi[key] = hit.surface()->mediumProperties().xi();
3200         str_chargePerCM[key] = siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster());
3201         str_clustSize[key] = hit.cluster()->amplitudes().size();
3202         str_usedMask[key] = strUsedMask(key);
3203         LogTrace("TrackingNtuple") << name << " cluster=" << key << " subdId=" << hitId.subdetId() << " lay=" << lay
3204                                    << " rawId=" << hitId.rawId() << " pos =" << hit.globalPosition();
3205 
3206         if (includeTrackingParticles_) {
3207           SimHitData simHitData = matchCluster(hit.firstClusterRef(),
3208                                                hitId,
3209                                                key,
3210                                                hit,
3211                                                clusterToTPMap,
3212                                                tpKeyToIndex,
3213                                                simHitsTPAssoc,
3214                                                digiSimLink,
3215                                                simHitRefKeyToIndex,
3216                                                HitType::Strip);
3217           str_simHitIdx[key] = simHitData.matchingSimHit;
3218           str_simType[key] = static_cast<int>(simHitData.type);
3219           str_xySignificance[key] = simHitData.xySignificance;
3220           str_chargeFraction[key] = simHitData.chargeFraction;
3221           LogTrace("TrackingNtuple") << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
3222           if (!simHitData.matchingSimHit.empty()) {
3223             const auto simHitIdx = simHitData.matchingSimHit[0];
3224             LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx << " simHitPos="
3225                                        << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3226                                        << " simHitPos="
3227                                        << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3228                                        << " energyLoss=" << simhit_eloss[simHitIdx]
3229                                        << " particleType=" << simhit_particle[simHitIdx]
3230                                        << " processType=" << simhit_process[simHitIdx]
3231                                        << " bunchCrossing=" << simHitData.bunchCrossing[0]
3232                                        << " event=" << simHitData.event[0];
3233           }
3234         }
3235       }
3236     }
3237   };
3238 
3239   fill(*rphiHits, "stripRPhiHit");
3240   fill(*stereoHits, "stripStereoHit");
3241 }
3242 
3243 size_t TrackingNtuple::addStripMatchedHit(const SiStripMatchedRecHit2D& hit,
3244                                           const TrackerGeometry& tracker,
3245                                           const TrackerTopology& tTopo,
3246                                           const std::vector<std::pair<uint64_t, StripMaskContainer const*>>& stripMasks,
3247                                           std::vector<std::pair<int, int>>& monoStereoClusterList) {
3248   auto strUsedMask = [&stripMasks](size_t key) {
3249     uint64_t mask = 0;
3250     for (auto const& m : stripMasks) {
3251       if (m.second->mask(key))
3252         mask |= m.first;
3253     }
3254     return mask;
3255   };
3256 
3257   const auto hitId = hit.geographicalId();
3258   const int lay = tTopo.layer(hitId);
3259   monoStereoClusterList.emplace_back(hit.monoHit().cluster().key(), hit.stereoHit().cluster().key());
3260   glu_isBarrel.push_back((hitId.subdetId() == StripSubdetector::TIB || hitId.subdetId() == StripSubdetector::TOB));
3261   glu_detId.push_back(tracker, tTopo, hitId);
3262   glu_monoIdx.push_back(hit.monoHit().cluster().key());
3263   glu_stereoIdx.push_back(hit.stereoHit().cluster().key());
3264   glu_seeIdx.emplace_back();  // filled in fillSeeds
3265   glu_x.push_back(hit.globalPosition().x());
3266   glu_y.push_back(hit.globalPosition().y());
3267   glu_z.push_back(hit.globalPosition().z());
3268   glu_xx.push_back(hit.globalPositionError().cxx());
3269   glu_xy.push_back(hit.globalPositionError().cyx());
3270   glu_yy.push_back(hit.globalPositionError().cyy());
3271   glu_yz.push_back(hit.globalPositionError().czy());
3272   glu_zz.push_back(hit.globalPositionError().czz());
3273   glu_zx.push_back(hit.globalPositionError().czx());
3274   glu_radL.push_back(hit.surface()->mediumProperties().radLen());
3275   glu_bbxi.push_back(hit.surface()->mediumProperties().xi());
3276   glu_chargePerCM.push_back(siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()));
3277   glu_clustSizeMono.push_back(hit.monoHit().cluster()->amplitudes().size());
3278   glu_clustSizeStereo.push_back(hit.stereoHit().cluster()->amplitudes().size());
3279   glu_usedMaskMono.push_back(strUsedMask(hit.monoHit().cluster().key()));
3280   glu_usedMaskStereo.push_back(strUsedMask(hit.stereoHit().cluster().key()));
3281   LogTrace("TrackingNtuple") << "stripMatchedHit"
3282                              << " cluster0=" << hit.stereoHit().cluster().key()
3283                              << " cluster1=" << hit.monoHit().cluster().key() << " subdId=" << hitId.subdetId()
3284                              << " lay=" << lay << " rawId=" << hitId.rawId() << " pos =" << hit.globalPosition();
3285   return glu_isBarrel.size() - 1;
3286 }
3287 
3288 void TrackingNtuple::fillStripMatchedHits(const edm::Event& iEvent,
3289                                           const TrackerGeometry& tracker,
3290                                           const TrackerTopology& tTopo,
3291                                           std::vector<std::pair<int, int>>& monoStereoClusterList) {
3292   std::vector<std::pair<uint64_t, StripMaskContainer const*>> stripMasks;
3293   stripMasks.reserve(stripUseMaskTokens_.size());
3294   for (const auto& itoken : stripUseMaskTokens_) {
3295     edm::Handle<StripMaskContainer> aH;
3296     iEvent.getByToken(itoken.second, aH);
3297     stripMasks.emplace_back(1 << itoken.first, aH.product());
3298   }
3299 
3300   edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
3301   iEvent.getByToken(stripMatchedRecHitToken_, matchedHits);
3302   for (auto it = matchedHits->begin(); it != matchedHits->end(); it++) {
3303     for (auto hit = it->begin(); hit != it->end(); hit++) {
3304       addStripMatchedHit(*hit, tracker, tTopo, stripMasks, monoStereoClusterList);
3305     }
3306   }
3307 }
3308 
3309 void TrackingNtuple::fillPhase2OTHits(const edm::Event& iEvent,
3310                                       const ClusterTPAssociation& clusterToTPMap,
3311                                       const TrackerGeometry& tracker,
3312                                       const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3313                                       const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
3314                                       const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
3315                                       const TrackerTopology& tTopo,
3316                                       const SimHitRefKeyToIndex& simHitRefKeyToIndex,
3317                                       std::set<edm::ProductID>& hitProductIds) {
3318   edm::Handle<Phase2TrackerRecHit1DCollectionNew> phase2OTHits;
3319   iEvent.getByToken(phase2OTRecHitToken_, phase2OTHits);
3320   for (auto it = phase2OTHits->begin(); it != phase2OTHits->end(); it++) {
3321     const DetId hitId = it->detId();
3322     for (auto hit = it->begin(); hit != it->end(); hit++) {
3323       hitProductIds.insert(hit->cluster().id());
3324 
3325       const int key = hit->cluster().key();
3326       const int lay = tTopo.layer(hitId);
3327 
3328       ph2_isBarrel.push_back(hitId.subdetId() == 1);
3329       ph2_detId.push_back(tracker, tTopo, hitId);
3330       ph2_trkIdx.emplace_back();  // filled in fillTracks
3331       ph2_onTrk_x.emplace_back();
3332       ph2_onTrk_y.emplace_back();
3333       ph2_onTrk_z.emplace_back();
3334       ph2_onTrk_xx.emplace_back();
3335       ph2_onTrk_xy.emplace_back();
3336       ph2_onTrk_yy.emplace_back();
3337       ph2_onTrk_yz.emplace_back();
3338       ph2_onTrk_zz.emplace_back();
3339       ph2_onTrk_zx.emplace_back();
3340       ph2_tcandIdx.emplace_back();  // filled in fillCandidates
3341       ph2_seeIdx.emplace_back();    // filled in fillSeeds
3342       ph2_x.push_back(hit->globalPosition().x());
3343       ph2_y.push_back(hit->globalPosition().y());
3344       ph2_z.push_back(hit->globalPosition().z());
3345       ph2_xx.push_back(hit->globalPositionError().cxx());
3346       ph2_xy.push_back(hit->globalPositionError().cyx());
3347       ph2_yy.push_back(hit->globalPositionError().cyy());
3348       ph2_yz.push_back(hit->globalPositionError().czy());
3349       ph2_zz.push_back(hit->globalPositionError().czz());
3350       ph2_zx.push_back(hit->globalPositionError().czx());
3351       ph2_radL.push_back(hit->surface()->mediumProperties().radLen());
3352       ph2_bbxi.push_back(hit->surface()->mediumProperties().xi());
3353 
3354       LogTrace("TrackingNtuple") << "phase2 OT cluster=" << key << " subdId=" << hitId.subdetId() << " lay=" << lay
3355                                  << " rawId=" << hitId.rawId() << " pos =" << hit->globalPosition();
3356 
3357       if (includeTrackingParticles_) {
3358         SimHitData simHitData = matchCluster(hit->firstClusterRef(),
3359                                              hitId,
3360                                              key,
3361                                              *hit,
3362                                              clusterToTPMap,
3363                                              tpKeyToIndex,
3364                                              simHitsTPAssoc,
3365                                              digiSimLink,
3366                                              simHitRefKeyToIndex,
3367                                              HitType::Phase2OT);
3368         ph2_xySignificance.push_back(simHitData.xySignificance);
3369         ph2_simHitIdx.push_back(simHitData.matchingSimHit);
3370         ph2_simType.push_back(static_cast<int>(simHitData.type));
3371         LogTrace("TrackingNtuple") << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
3372         if (!simHitData.matchingSimHit.empty()) {
3373           const auto simHitIdx = simHitData.matchingSimHit[0];
3374           LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx << " simHitPos="
3375                                      << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3376                                      << " energyLoss=" << simhit_eloss[simHitIdx]
3377                                      << " particleType=" << simhit_particle[simHitIdx]
3378                                      << " processType=" << simhit_process[simHitIdx]
3379                                      << " bunchCrossing=" << simHitData.bunchCrossing[0]
3380                                      << " event=" << simHitData.event[0];
3381         }
3382       }
3383     }
3384   }
3385 }
3386 
3387 void TrackingNtuple::fillSeeds(const edm::Event& iEvent,
3388                                const TrackingParticleRefVector& tpCollection,
3389                                const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3390                                const reco::BeamSpot& bs,
3391                                const TrackerGeometry& tracker,
3392                                const reco::TrackToTrackingParticleAssociator& associatorByHits,
3393                                const ClusterTPAssociation& clusterToTPMap,
3394                                const MagneticField& theMF,
3395                                const TrackerTopology& tTopo,
3396                                std::vector<std::pair<int, int>>& monoStereoClusterList,
3397                                const std::set<edm::ProductID>& hitProductIds,
3398                                std::map<edm::ProductID, size_t>& seedCollToOffset) {
3399   TSCBLBuilderNoMaterial tscblBuilder;
3400   for (size_t iColl = 0; iColl < seedTokens_.size(); ++iColl) {
3401     const auto& seedToken = seedTokens_[iColl];
3402 
3403     edm::Handle<edm::View<reco::Track>> seedTracksHandle;
3404     iEvent.getByToken(seedToken, seedTracksHandle);
3405     const auto& seedTracks = *seedTracksHandle;
3406 
3407     if (seedTracks.empty())
3408       continue;
3409 
3410     edm::EDConsumerBase::Labels labels;
3411     labelsForToken(seedToken, labels);
3412 
3413     const auto& seedStopInfoToken = seedStopInfoTokens_[iColl];
3414     edm::Handle<std::vector<SeedStopInfo>> seedStopInfoHandle;
3415     iEvent.getByToken(seedStopInfoToken, seedStopInfoHandle);
3416     const auto& seedStopInfos = *seedStopInfoHandle;
3417     if (seedTracks.size() != seedStopInfos.size()) {
3418       edm::EDConsumerBase::Labels labels2;
3419       labelsForToken(seedStopInfoToken, labels2);
3420 
3421       throw cms::Exception("LogicError") << "Got " << seedTracks.size() << " seeds, but " << seedStopInfos.size()
3422                                          << " seed stopping infos for collections " << labels.module << ", "
3423                                          << labels2.module;
3424     }
3425 
3426     std::vector<std::pair<uint64_t, StripMaskContainer const*>> stripMasks;
3427     stripMasks.reserve(stripUseMaskTokens_.size());
3428     for (const auto& itoken : stripUseMaskTokens_) {
3429       edm::Handle<StripMaskContainer> aH;
3430       iEvent.getByToken(itoken.second, aH);
3431       stripMasks.emplace_back(1 << itoken.first, aH.product());
3432     }
3433 
3434     // The associator interfaces really need to be fixed...
3435     edm::RefToBaseVector<reco::Track> seedTrackRefs;
3436     for (edm::View<reco::Track>::size_type i = 0; i < seedTracks.size(); ++i) {
3437       seedTrackRefs.push_back(seedTracks.refAt(i));
3438     }
3439     reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(seedTrackRefs, tpCollection);
3440     reco::SimToRecoCollection simRecColl = associatorByHits.associateSimToReco(seedTrackRefs, tpCollection);
3441 
3442     TString label = labels.module;
3443     //format label to match algoName
3444     label.ReplaceAll("seedTracks", "");
3445     label.ReplaceAll("Seeds", "");
3446     label.ReplaceAll("muonSeeded", "muonSeededStep");
3447     //for HLT seeds
3448     label.ReplaceAll("FromPixelTracks", "");
3449     label.ReplaceAll("PFLowPixel", "");
3450     label.ReplaceAll("hltDoubletRecovery", "pixelPairStep");  //random choice
3451     int algo = reco::TrackBase::algoByName(label.Data());
3452 
3453     edm::ProductID id = seedTracks[0].seedRef().id();
3454     const auto offset = see_fitok.size();
3455     auto inserted = seedCollToOffset.emplace(id, offset);
3456     if (!inserted.second)
3457       throw cms::Exception("Configuration")
3458           << "Trying to add seeds with ProductID " << id << " for a second time from collection " << labels.module
3459           << ", seed algo " << label << ". Typically this is caused by a configuration problem.";
3460     see_offset.push_back(offset);
3461 
3462     LogTrace("TrackingNtuple") << "NEW SEED LABEL: " << label << " size: " << seedTracks.size() << " algo=" << algo
3463                                << " ProductID " << id;
3464 
3465     for (size_t iSeed = 0; iSeed < seedTrackRefs.size(); ++iSeed) {
3466       const auto& seedTrackRef = seedTrackRefs[iSeed];
3467       const auto& seedTrack = *seedTrackRef;
3468       const auto& seedRef = seedTrack.seedRef();
3469       const auto& seed = *seedRef;
3470 
3471       const auto seedStopInfo = seedStopInfos[iSeed];
3472 
3473       if (seedRef.id() != id)
3474         throw cms::Exception("LogicError")
3475             << "All tracks in 'TracksFromSeeds' collection should point to seeds in the same collection. Now the "
3476                "element 0 had ProductID "
3477             << id << " while the element " << seedTrackRef.key() << " had " << seedTrackRef.id()
3478             << ". The source collection is " << labels.module << ".";
3479 
3480       std::vector<int> tpIdx;
3481       std::vector<float> sharedFraction;
3482       auto foundTPs = recSimColl.find(seedTrackRef);
3483       if (foundTPs != recSimColl.end()) {
3484         for (const auto& tpQuality : foundTPs->val) {
3485           tpIdx.push_back(tpKeyToIndex.at(tpQuality.first.key()));
3486           sharedFraction.push_back(tpQuality.second);
3487         }
3488       }
3489 
3490       // Search for a best-matching TrackingParticle for a seed
3491       const int nHits = seedTrack.numberOfValidHits();
3492       const auto clusters = track_associator::hitsToClusterRefs(
3493           seedTrack.recHitsBegin(),
3494           seedTrack.recHitsEnd());  // TODO: this function is called 3 times per track, try to reduce
3495       const int nClusters = clusters.size();
3496       const auto bestKeyCount = findBestMatchingTrackingParticle(seedTrack.recHits(), clusterToTPMap, tpKeyToIndex);
3497       const float bestShareFrac =
3498           nClusters > 0 ? static_cast<float>(bestKeyCount.countClusters) / static_cast<float>(nClusters) : 0;
3499       // Another way starting from the first hit of the seed
3500       const auto bestFirstHitKeyCount =
3501           findMatchingTrackingParticleFromFirstHit(seedTrack.recHits(), clusterToTPMap, tpKeyToIndex);
3502       const float bestFirstHitShareFrac =
3503           nClusters > 0 ? static_cast<float>(bestFirstHitKeyCount.countClusters) / static_cast<float>(nClusters) : 0;
3504 
3505       const bool seedFitOk = !trackFromSeedFitFailed(seedTrack);
3506       const int charge = seedTrack.charge();
3507       const float pt = seedFitOk ? seedTrack.pt() : 0;
3508       const float eta = seedFitOk ? seedTrack.eta() : 0;
3509       const float phi = seedFitOk ? seedTrack.phi() : 0;
3510 
3511       const auto seedIndex = see_fitok.size();
3512 
3513       see_fitok.push_back(seedFitOk);
3514 
3515       see_px.push_back(seedFitOk ? seedTrack.px() : 0);
3516       see_py.push_back(seedFitOk ? seedTrack.py() : 0);
3517       see_pz.push_back(seedFitOk ? seedTrack.pz() : 0);
3518       see_pt.push_back(pt);
3519       see_eta.push_back(eta);
3520       see_phi.push_back(phi);
3521       see_q.push_back(charge);
3522       see_nValid.push_back(nHits);
3523 
3524       see_dxy.push_back(seedFitOk ? seedTrack.dxy(bs.position()) : 0);
3525       see_dz.push_back(seedFitOk ? seedTrack.dz(bs.position()) : 0);
3526       see_ptErr.push_back(seedFitOk ? seedTrack.ptError() : 0);
3527       see_etaErr.push_back(seedFitOk ? seedTrack.etaError() : 0);
3528       see_phiErr.push_back(seedFitOk ? seedTrack.phiError() : 0);
3529       see_dxyErr.push_back(seedFitOk ? seedTrack.dxyError() : 0);
3530       see_dzErr.push_back(seedFitOk ? seedTrack.dzError() : 0);
3531       see_algo.push_back(algo);
3532       see_stopReason.push_back(seedStopInfo.stopReasonUC());
3533       see_nCands.push_back(seedStopInfo.candidatesPerSeed());
3534 
3535       const auto& state = seedTrack.seedRef()->startingState();
3536       const auto& pos = state.parameters().position();
3537       const auto& mom = state.parameters().momentum();
3538       see_statePt.push_back(state.pt());
3539       see_stateTrajX.push_back(pos.x());
3540       see_stateTrajY.push_back(pos.y());
3541       see_stateTrajPx.push_back(mom.x());
3542       see_stateTrajPy.push_back(mom.y());
3543       see_stateTrajPz.push_back(mom.z());
3544 
3545       ///the following is useful for analysis in global coords at seed hit surface
3546       const TrackingRecHit* lastRecHit = &*(seed.recHits().end() - 1);
3547       TrajectoryStateOnSurface tsos =
3548           trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), &theMF);
3549       auto const& stateGlobal = tsos.globalParameters();
3550       see_stateTrajGlbX.push_back(stateGlobal.position().x());
3551       see_stateTrajGlbY.push_back(stateGlobal.position().y());
3552       see_stateTrajGlbZ.push_back(stateGlobal.position().z());
3553       see_stateTrajGlbPx.push_back(stateGlobal.momentum().x());
3554       see_stateTrajGlbPy.push_back(stateGlobal.momentum().y());
3555       see_stateTrajGlbPz.push_back(stateGlobal.momentum().z());
3556       if (addSeedCurvCov_) {
3557         auto const& stateCcov = tsos.curvilinearError().matrix();
3558         std::vector<float> cov(15);
3559         auto covP = cov.begin();
3560         for (auto const val : stateCcov)
3561           *(covP++) = val;  //row-major
3562         see_stateCurvCov.push_back(std::move(cov));
3563       }
3564 
3565       see_trkIdx.push_back(-1);    // to be set correctly in fillTracks
3566       see_tcandIdx.push_back(-1);  // to be set correctly in fillCandidates
3567       if (includeTrackingParticles_) {
3568         see_simTrkIdx.push_back(tpIdx);
3569         see_simTrkShareFrac.push_back(sharedFraction);
3570         see_bestSimTrkIdx.push_back(bestKeyCount.key >= 0 ? tpKeyToIndex.at(bestKeyCount.key) : -1);
3571         see_bestFromFirstHitSimTrkIdx.push_back(
3572             bestFirstHitKeyCount.key >= 0 ? tpKeyToIndex.at(bestFirstHitKeyCount.key) : -1);
3573       } else {
3574         see_isTrue.push_back(!tpIdx.empty());
3575       }
3576       see_bestSimTrkShareFrac.push_back(bestShareFrac);
3577       see_bestFromFirstHitSimTrkShareFrac.push_back(bestFirstHitShareFrac);
3578 
3579       /// Hmm, the following could make sense instead of plain failing if propagation to beam line fails
3580       /*
3581       const TrackingRecHit* lastRecHit = &*(seed.recHits().second-1);
3582       TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( itSeed->startingState(), lastRecHit->surface(), &theMF);
3583       float pt  = state.globalParameters().momentum().perp();
3584       float eta = state.globalParameters().momentum().eta();
3585       float phi = state.globalParameters().momentum().phi();
3586       see_px      .push_back( state.globalParameters().momentum().x() );
3587       see_py      .push_back( state.globalParameters().momentum().y() );
3588       see_pz      .push_back( state.globalParameters().momentum().z() );
3589       */
3590 
3591       std::vector<int> hitIdx;
3592       std::vector<int> hitType;
3593 
3594       for (auto const& hit : seed.recHits()) {
3595         int subid = hit.geographicalId().subdetId();
3596         if (subid == (int)PixelSubdetector::PixelBarrel || subid == (int)PixelSubdetector::PixelEndcap) {
3597           const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&hit);
3598           const auto& clusterRef = bhit->firstClusterRef();
3599           const auto clusterKey = clusterRef.cluster_pixel().key();
3600           if (includeAllHits_) {
3601             checkProductID(hitProductIds, clusterRef.id(), "seed");
3602             pix_seeIdx[clusterKey].push_back(seedIndex);
3603           }
3604           hitIdx.push_back(clusterKey);
3605           hitType.push_back(static_cast<int>(HitType::Pixel));
3606         } else if (subid == (int)StripSubdetector::TOB || subid == (int)StripSubdetector::TID ||
3607                    subid == (int)StripSubdetector::TIB || subid == (int)StripSubdetector::TEC) {
3608           if (trackerHitRTTI::isMatched(hit)) {
3609             const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(&hit);
3610             if (includeAllHits_) {
3611               checkProductID(hitProductIds, matchedHit->monoClusterRef().id(), "seed");
3612               checkProductID(hitProductIds, matchedHit->stereoClusterRef().id(), "seed");
3613             }
3614             int monoIdx = matchedHit->monoClusterRef().key();
3615             int stereoIdx = matchedHit->stereoClusterRef().key();
3616 
3617             std::vector<std::pair<int, int>>::iterator pos =
3618                 find(monoStereoClusterList.begin(), monoStereoClusterList.end(), std::make_pair(monoIdx, stereoIdx));
3619             size_t gluedIndex = -1;
3620             if (pos != monoStereoClusterList.end()) {
3621               gluedIndex = std::distance(monoStereoClusterList.begin(), pos);
3622             } else {
3623               // We can encounter glued hits not in the input
3624               // SiStripMatchedRecHit2DCollection, e.g. via muon
3625               // outside-in seeds (or anything taking hits from
3626               // MeasurementTrackerEvent). So let's add them here.
3627               gluedIndex = addStripMatchedHit(*matchedHit, tracker, tTopo, stripMasks, monoStereoClusterList);
3628             }
3629 
3630             if (includeAllHits_)
3631               glu_seeIdx[gluedIndex].push_back(seedIndex);
3632             hitIdx.push_back(gluedIndex);
3633             hitType.push_back(static_cast<int>(HitType::Glued));
3634           } else {
3635             const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&hit);
3636             const auto& clusterRef = bhit->firstClusterRef();
3637             unsigned int clusterKey;
3638             if (clusterRef.isPhase2()) {
3639               clusterKey = clusterRef.cluster_phase2OT().key();
3640             } else {
3641               clusterKey = clusterRef.cluster_strip().key();
3642             }
3643 
3644             if (includeAllHits_) {
3645               checkProductID(hitProductIds, clusterRef.id(), "seed");
3646               if (clusterRef.isPhase2()) {
3647                 ph2_seeIdx[clusterKey].push_back(seedIndex);
3648               } else {
3649                 str_seeIdx[clusterKey].push_back(seedIndex);
3650               }
3651             }
3652 
3653             hitIdx.push_back(clusterKey);
3654             if (clusterRef.isPhase2()) {
3655               hitType.push_back(static_cast<int>(HitType::Phase2OT));
3656             } else {
3657               hitType.push_back(static_cast<int>(HitType::Strip));
3658             }
3659           }
3660         } else {
3661           LogTrace("TrackingNtuple") << " not pixel and not Strip detector";
3662         }
3663       }
3664       see_hitIdx.push_back(hitIdx);
3665       see_hitType.push_back(hitType);
3666       see_nPixel.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Pixel)));
3667       see_nGlued.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Glued)));
3668       see_nStrip.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Strip)));
3669       see_nPhase2OT.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Phase2OT)));
3670       see_nCluster.push_back(nClusters);
3671       //the part below is not strictly needed
3672       float chi2 = -1;
3673       if (nHits == 2) {
3674         auto const recHit0 = seed.recHits().begin();
3675         auto const recHit1 = seed.recHits().begin() + 1;
3676         std::vector<GlobalPoint> gp(2);
3677         std::vector<GlobalError> ge(2);
3678         gp[0] = recHit0->globalPosition();
3679         ge[0] = recHit0->globalPositionError();
3680         gp[1] = recHit1->globalPosition();
3681         ge[1] = recHit1->globalPositionError();
3682         LogTrace("TrackingNtuple")
3683             << "seed " << seedTrackRef.key() << " pt=" << pt << " eta=" << eta << " phi=" << phi << " q=" << charge
3684             << " - PAIR - ids: " << recHit0->geographicalId().rawId() << " " << recHit1->geographicalId().rawId()
3685             << " hitpos: " << gp[0] << " " << gp[1] << " trans0: "
3686             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[0]->globalPosition()
3687                                                     : GlobalPoint(0, 0, 0))
3688             << " "
3689             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[1]->globalPosition()
3690                                                     : GlobalPoint(0, 0, 0))
3691             << " trans1: "
3692             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[0]->globalPosition()
3693                                                     : GlobalPoint(0, 0, 0))
3694             << " "
3695             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[1]->globalPosition()
3696                                                     : GlobalPoint(0, 0, 0))
3697             << " eta,phi: " << gp[0].eta() << "," << gp[0].phi();
3698       } else if (nHits == 3) {
3699         auto const recHit0 = seed.recHits().begin();
3700         auto const recHit1 = seed.recHits().begin() + 1;
3701         auto const recHit2 = seed.recHits().begin() + 2;
3702         declareDynArray(GlobalPoint, 4, gp);
3703         declareDynArray(GlobalError, 4, ge);
3704         declareDynArray(bool, 4, bl);
3705         gp[0] = recHit0->globalPosition();
3706         ge[0] = recHit0->globalPositionError();
3707         int subid0 = recHit0->geographicalId().subdetId();
3708         bl[0] = (subid0 == StripSubdetector::TIB || subid0 == StripSubdetector::TOB ||
3709                  subid0 == (int)PixelSubdetector::PixelBarrel);
3710         gp[1] = recHit1->globalPosition();
3711         ge[1] = recHit1->globalPositionError();
3712         int subid1 = recHit1->geographicalId().subdetId();
3713         bl[1] = (subid1 == StripSubdetector::TIB || subid1 == StripSubdetector::TOB ||
3714                  subid1 == (int)PixelSubdetector::PixelBarrel);
3715         gp[2] = recHit2->globalPosition();
3716         ge[2] = recHit2->globalPositionError();
3717         int subid2 = recHit2->geographicalId().subdetId();
3718         bl[2] = (subid2 == StripSubdetector::TIB || subid2 == StripSubdetector::TOB ||
3719                  subid2 == (int)PixelSubdetector::PixelBarrel);
3720         RZLine rzLine(gp, ge, bl);
3721         float seed_chi2 = rzLine.chi2();
3722         //float seed_pt = state.globalParameters().momentum().perp();
3723         float seed_pt = pt;
3724         LogTrace("TrackingNtuple")
3725             << "seed " << seedTrackRef.key() << " pt=" << pt << " eta=" << eta << " phi=" << phi << " q=" << charge
3726             << " - TRIPLET - ids: " << recHit0->geographicalId().rawId() << " " << recHit1->geographicalId().rawId()
3727             << " " << recHit2->geographicalId().rawId() << " hitpos: " << gp[0] << " " << gp[1] << " " << gp[2]
3728             << " trans0: "
3729             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[0]->globalPosition()
3730                                                     : GlobalPoint(0, 0, 0))
3731             << " "
3732             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[1]->globalPosition()
3733                                                     : GlobalPoint(0, 0, 0))
3734             << " trans1: "
3735             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[0]->globalPosition()
3736                                                     : GlobalPoint(0, 0, 0))
3737             << " "
3738             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[1]->globalPosition()
3739                                                     : GlobalPoint(0, 0, 0))
3740             << " trans2: "
3741             << (recHit2->transientHits().size() > 1 ? recHit2->transientHits()[0]->globalPosition()
3742                                                     : GlobalPoint(0, 0, 0))
3743             << " "
3744             << (recHit2->transientHits().size() > 1 ? recHit2->transientHits()[1]->globalPosition()
3745                                                     : GlobalPoint(0, 0, 0))
3746             << " local: "
3747             << recHit2->localPosition()
3748             //<< " tsos pos, mom: " << state.globalPosition()<<" "<<state.globalMomentum()
3749             << " eta,phi: " << gp[0].eta() << "," << gp[0].phi() << " pt,chi2: " << seed_pt << "," << seed_chi2;
3750         chi2 = seed_chi2;
3751       }
3752       see_chi2.push_back(chi2);
3753     }
3754 
3755     fillTrackingParticlesForSeeds(tpCollection, simRecColl, tpKeyToIndex, offset);
3756   }
3757 }
3758 
3759 void TrackingNtuple::fillTracks(const edm::RefToBaseVector<reco::Track>& tracks,
3760                                 const TrackerGeometry& tracker,
3761                                 const TrackingParticleRefVector& tpCollection,
3762                                 const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3763                                 const TrackingParticleRefKeyToCount& tpKeyToClusterCount,
3764                                 const MagneticField& mf,
3765                                 const reco::BeamSpot& bs,
3766                                 const reco::VertexCollection& vertices,
3767                                 const reco::TrackToTrackingParticleAssociator& associatorByHits,
3768                                 const ClusterTPAssociation& clusterToTPMap,
3769                                 const TrackerTopology& tTopo,
3770                                 const std::set<edm::ProductID>& hitProductIds,
3771                                 const std::map<edm::ProductID, size_t>& seedCollToOffset,
3772                                 const std::vector<const MVACollection*>& mvaColls,
3773                                 const std::vector<const QualityMaskCollection*>& qualColls) {
3774   reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(tracks, tpCollection);
3775   edm::EDConsumerBase::Labels labels;
3776   labelsForToken(trackToken_, labels);
3777   LogTrace("TrackingNtuple") << "NEW TRACK LABEL: " << labels.module;
3778 
3779   auto pvPosition = vertices[0].position();
3780 
3781   for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) {
3782     const auto& itTrack = tracks[iTrack];
3783     int charge = itTrack->charge();
3784     float pt = itTrack->pt();
3785     float eta = itTrack->eta();
3786     const double lambda = itTrack->lambda();
3787     float chi2 = itTrack->normalizedChi2();
3788     float ndof = itTrack->ndof();
3789     float phi = itTrack->phi();
3790     int nHits = itTrack->numberOfValidHits();
3791     const reco::HitPattern& hp = itTrack->hitPattern();
3792 
3793     const auto& tkParam = itTrack->parameters();
3794     auto tkCov = itTrack->covariance();
3795     tkCov.Invert();
3796 
3797     // Standard track-TP matching
3798     int nSimHits = 0;
3799     bool isSimMatched = false;
3800     std::vector<int> tpIdx;
3801     std::vector<float> sharedFraction;
3802     std::vector<float> tpChi2;
3803     auto foundTPs = recSimColl.find(itTrack);
3804     if (foundTPs != recSimColl.end()) {
3805       if (!foundTPs->val.empty()) {
3806         nSimHits = foundTPs->val[0].first->numberOfTrackerHits();
3807         isSimMatched = true;
3808       }
3809       for (const auto& tpQuality : foundTPs->val) {
3810         tpIdx.push_back(tpKeyToIndex.at(tpQuality.first.key()));
3811         sharedFraction.push_back(tpQuality.second);
3812         tpChi2.push_back(track_associator::trackAssociationChi2(tkParam, tkCov, *(tpCollection[tpIdx.back()]), mf, bs));
3813       }
3814     }
3815 
3816     // Search for a best-matching TrackingParticle for a track
3817     const auto clusters = track_associator::hitsToClusterRefs(
3818         itTrack->recHitsBegin(),
3819         itTrack->recHitsEnd());  // TODO: this function is called 3 times per track, try to reduce
3820     const int nClusters = clusters.size();
3821 
3822     const auto bestKeyCount = findBestMatchingTrackingParticle(itTrack->recHits(), clusterToTPMap, tpKeyToIndex);
3823     const float bestShareFrac = static_cast<float>(bestKeyCount.countClusters) / static_cast<float>(nClusters);
3824     float bestShareFracSimDenom = 0;
3825     float bestShareFracSimClusterDenom = 0;
3826     float bestChi2 = -1;
3827     if (bestKeyCount.key >= 0) {
3828       bestShareFracSimDenom =
3829           static_cast<float>(bestKeyCount.countClusters) /
3830           static_cast<float>(tpCollection[tpKeyToIndex.at(bestKeyCount.key)]->numberOfTrackerHits());
3831       bestShareFracSimClusterDenom =
3832           static_cast<float>(bestKeyCount.countClusters) / static_cast<float>(tpKeyToClusterCount.at(bestKeyCount.key));
3833       bestChi2 = track_associator::trackAssociationChi2(
3834           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestKeyCount.key)]), mf, bs);
3835     }
3836     // Another way starting from the first hit of the track
3837     const auto bestFirstHitKeyCount =
3838         findMatchingTrackingParticleFromFirstHit(itTrack->recHits(), clusterToTPMap, tpKeyToIndex);
3839     const float bestFirstHitShareFrac =
3840         static_cast<float>(bestFirstHitKeyCount.countClusters) / static_cast<float>(nClusters);
3841     float bestFirstHitShareFracSimDenom = 0;
3842     float bestFirstHitShareFracSimClusterDenom = 0;
3843     float bestFirstHitChi2 = -1;
3844     if (bestFirstHitKeyCount.key >= 0) {
3845       bestFirstHitShareFracSimDenom =
3846           static_cast<float>(bestFirstHitKeyCount.countClusters) /
3847           static_cast<float>(tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]->numberOfTrackerHits());
3848       bestFirstHitShareFracSimClusterDenom = static_cast<float>(bestFirstHitKeyCount.countClusters) /
3849                                              static_cast<float>(tpKeyToClusterCount.at(bestFirstHitKeyCount.key));
3850       bestFirstHitChi2 = track_associator::trackAssociationChi2(
3851           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]), mf, bs);
3852     }
3853 
3854     float chi2_1Dmod = chi2;
3855     int count1dhits = 0;
3856     for (auto iHit = itTrack->recHitsBegin(), iEnd = itTrack->recHitsEnd(); iHit != iEnd; ++iHit) {
3857       const TrackingRecHit& hit = **iHit;
3858       if (hit.isValid() && typeid(hit) == typeid(SiStripRecHit1D))
3859         ++count1dhits;
3860     }
3861     if (count1dhits > 0) {
3862       chi2_1Dmod = (chi2 + count1dhits) / (ndof + count1dhits);
3863     }
3864 
3865     Point bestPV = getBestVertex(*itTrack, vertices);
3866 
3867     trk_px.push_back(itTrack->px());
3868     trk_py.push_back(itTrack->py());
3869     trk_pz.push_back(itTrack->pz());
3870     trk_pt.push_back(pt);
3871     trk_inner_px.push_back(itTrack->innerMomentum().x());
3872     trk_inner_py.push_back(itTrack->innerMomentum().y());
3873     trk_inner_pz.push_back(itTrack->innerMomentum().z());
3874     trk_inner_pt.push_back(itTrack->innerMomentum().rho());
3875     trk_outer_px.push_back(itTrack->outerMomentum().x());
3876     trk_outer_py.push_back(itTrack->outerMomentum().y());
3877     trk_outer_pz.push_back(itTrack->outerMomentum().z());
3878     trk_outer_pt.push_back(itTrack->outerMomentum().rho());
3879     trk_eta.push_back(eta);
3880     trk_lambda.push_back(lambda);
3881     trk_cotTheta.push_back(1 / tan(M_PI * 0.5 - lambda));
3882     trk_phi.push_back(phi);
3883     trk_dxy.push_back(itTrack->dxy(bs.position()));
3884     trk_dz.push_back(itTrack->dz(bs.position()));
3885     trk_dxyPV.push_back(itTrack->dxy(pvPosition));
3886     trk_dzPV.push_back(itTrack->dz(pvPosition));
3887     trk_dxyClosestPV.push_back(itTrack->dxy(bestPV));
3888     trk_dzClosestPV.push_back(itTrack->dz(bestPV));
3889     trk_ptErr.push_back(itTrack->ptError());
3890     trk_etaErr.push_back(itTrack->etaError());
3891     trk_lambdaErr.push_back(itTrack->lambdaError());
3892     trk_phiErr.push_back(itTrack->phiError());
3893     trk_dxyErr.push_back(itTrack->dxyError());
3894     trk_dzErr.push_back(itTrack->dzError());
3895     trk_refpoint_x.push_back(itTrack->vx());
3896     trk_refpoint_y.push_back(itTrack->vy());
3897     trk_refpoint_z.push_back(itTrack->vz());
3898     trk_nChi2.push_back(chi2);
3899     trk_nChi2_1Dmod.push_back(chi2_1Dmod);
3900     trk_ndof.push_back(ndof);
3901     trk_q.push_back(charge);
3902     trk_nValid.push_back(hp.numberOfValidHits());
3903     trk_nLost.push_back(hp.numberOfLostHits(reco::HitPattern::TRACK_HITS));
3904     trk_nInactive.push_back(hp.trackerLayersTotallyOffOrBad(reco::HitPattern::TRACK_HITS));
3905     trk_nPixel.push_back(hp.numberOfValidPixelHits());
3906     trk_nStrip.push_back(hp.numberOfValidStripHits());
3907     trk_nOuterLost.push_back(hp.numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS));
3908     trk_nInnerLost.push_back(hp.numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS));
3909     trk_nOuterInactive.push_back(hp.trackerLayersTotallyOffOrBad(reco::HitPattern::MISSING_OUTER_HITS));
3910     trk_nInnerInactive.push_back(hp.trackerLayersTotallyOffOrBad(reco::HitPattern::MISSING_INNER_HITS));
3911     trk_nPixelLay.push_back(hp.pixelLayersWithMeasurement());
3912     trk_nStripLay.push_back(hp.stripLayersWithMeasurement());
3913     trk_n3DLay.push_back(hp.numberOfValidStripLayersWithMonoAndStereo() + hp.pixelLayersWithMeasurement());
3914     trk_nLostLay.push_back(hp.trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS));
3915     trk_nCluster.push_back(nClusters);
3916     trk_algo.push_back(itTrack->algo());
3917     trk_originalAlgo.push_back(itTrack->originalAlgo());
3918     trk_algoMask.push_back(itTrack->algoMaskUL());
3919     trk_stopReason.push_back(itTrack->stopReason());
3920     trk_isHP.push_back(itTrack->quality(reco::TrackBase::highPurity));
3921     if (includeMVA_) {
3922       for (size_t i = 0; i < trk_mvas.size(); ++i) {
3923         trk_mvas[i].push_back((*(mvaColls[i]))[iTrack]);
3924         trk_qualityMasks[i].push_back((*(qualColls[i]))[iTrack]);
3925       }
3926     }
3927     if (includeSeeds_) {
3928       auto offset = seedCollToOffset.find(itTrack->seedRef().id());
3929       if (offset == seedCollToOffset.end()) {
3930         throw cms::Exception("Configuration")
3931             << "Track algo '" << reco::TrackBase::algoName(itTrack->algo()) << "' originalAlgo '"
3932             << reco::TrackBase::algoName(itTrack->originalAlgo()) << "' refers to seed collection "
3933             << itTrack->seedRef().id()
3934             << ", but that seed collection is not given as an input. The following collections were given as an input "
3935             << make_ProductIDMapPrinter(seedCollToOffset);
3936       }
3937 
3938       const auto seedIndex = offset->second + itTrack->seedRef().key();
3939       trk_seedIdx.push_back(seedIndex);
3940       if (see_trkIdx[seedIndex] != -1) {
3941         throw cms::Exception("LogicError") << "Track index has already been set for seed " << seedIndex << " to "
3942                                            << see_trkIdx[seedIndex] << "; was trying to set it to " << iTrack;
3943       }
3944       see_trkIdx[seedIndex] = iTrack;
3945     }
3946     trk_vtxIdx.push_back(-1);  // to be set correctly in fillVertices
3947     if (includeTrackingParticles_) {
3948       trk_simTrkIdx.push_back(tpIdx);
3949       trk_simTrkShareFrac.push_back(sharedFraction);
3950       trk_simTrkNChi2.push_back(tpChi2);
3951       trk_bestSimTrkIdx.push_back(bestKeyCount.key >= 0 ? tpKeyToIndex.at(bestKeyCount.key) : -1);
3952       trk_bestFromFirstHitSimTrkIdx.push_back(bestFirstHitKeyCount.key >= 0 ? tpKeyToIndex.at(bestFirstHitKeyCount.key)
3953                                                                             : -1);
3954     } else {
3955       trk_isTrue.push_back(!tpIdx.empty());
3956     }
3957     trk_bestSimTrkShareFrac.push_back(bestShareFrac);
3958     trk_bestSimTrkShareFracSimDenom.push_back(bestShareFracSimDenom);
3959     trk_bestSimTrkShareFracSimClusterDenom.push_back(bestShareFracSimClusterDenom);
3960     trk_bestSimTrkNChi2.push_back(bestChi2);
3961     trk_bestFromFirstHitSimTrkShareFrac.push_back(bestFirstHitShareFrac);
3962     trk_bestFromFirstHitSimTrkShareFracSimDenom.push_back(bestFirstHitShareFracSimDenom);
3963     trk_bestFromFirstHitSimTrkShareFracSimClusterDenom.push_back(bestFirstHitShareFracSimClusterDenom);
3964     trk_bestFromFirstHitSimTrkNChi2.push_back(bestFirstHitChi2);
3965 
3966     LogTrace("TrackingNtuple") << "Track #" << itTrack.key() << " with q=" << charge << ", pT=" << pt
3967                                << " GeV, eta: " << eta << ", phi: " << phi << ", chi2=" << chi2 << ", Nhits=" << nHits
3968                                << ", algo=" << itTrack->algoName(itTrack->algo()).c_str()
3969                                << " hp=" << itTrack->quality(reco::TrackBase::highPurity)
3970                                << " seed#=" << itTrack->seedRef().key() << " simMatch=" << isSimMatched
3971                                << " nSimHits=" << nSimHits
3972                                << " sharedFraction=" << (sharedFraction.empty() ? -1 : sharedFraction[0])
3973                                << " tpIdx=" << (tpIdx.empty() ? -1 : tpIdx[0]);
3974     std::vector<int> hitIdx;
3975     std::vector<int> hitType;
3976 
3977     for (auto i = itTrack->recHitsBegin(); i != itTrack->recHitsEnd(); i++) {
3978       auto hit = *i;
3979       DetId hitId = hit->geographicalId();
3980       LogTrace("TrackingNtuple") << "hit #" << std::distance(itTrack->recHitsBegin(), i)
3981                                  << " subdet=" << hitId.subdetId();
3982       if (hitId.det() != DetId::Tracker)
3983         continue;
3984 
3985       LogTrace("TrackingNtuple") << " " << subdetstring(hitId.subdetId()) << " " << tTopo.layer(hitId);
3986 
3987       if (hit->isValid()) {
3988         //ugly... but works
3989         const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*hit);
3990         const auto& clusterRef = bhit->firstClusterRef();
3991         unsigned int clusterKey;
3992         if (clusterRef.isPixel()) {
3993           clusterKey = clusterRef.cluster_pixel().key();
3994         } else if (clusterRef.isPhase2()) {
3995           clusterKey = clusterRef.cluster_phase2OT().key();
3996         } else {
3997           clusterKey = clusterRef.cluster_strip().key();
3998         }
3999 
4000         LogTrace("TrackingNtuple") << " id: " << hitId.rawId() << " - globalPos =" << hit->globalPosition()
4001                                    << " cluster=" << clusterKey << " clusterRef ID=" << clusterRef.id()
4002                                    << " eta,phi: " << hit->globalPosition().eta() << "," << hit->globalPosition().phi();
4003         if (includeAllHits_) {
4004           checkProductID(hitProductIds, clusterRef.id(), "track");
4005           if (clusterRef.isPixel()) {
4006             pix_trkIdx[clusterKey].push_back(iTrack);
4007             pix_onTrk_x[clusterKey].push_back(hit->globalPosition().x());
4008             pix_onTrk_y[clusterKey].push_back(hit->globalPosition().y());
4009             pix_onTrk_z[clusterKey].push_back(hit->globalPosition().z());
4010             pix_onTrk_xx[clusterKey].push_back(hit->globalPositionError().cxx());
4011             pix_onTrk_xy[clusterKey].push_back(hit->globalPositionError().cyx());
4012             pix_onTrk_yy[clusterKey].push_back(hit->globalPositionError().cyy());
4013             pix_onTrk_yz[clusterKey].push_back(hit->globalPositionError().czy());
4014             pix_onTrk_zz[clusterKey].push_back(hit->globalPositionError().czz());
4015             pix_onTrk_zx[clusterKey].push_back(hit->globalPositionError().czx());
4016           } else if (clusterRef.isPhase2()) {
4017             ph2_trkIdx[clusterKey].push_back(iTrack);
4018             ph2_onTrk_x[clusterKey].push_back(hit->globalPosition().x());
4019             ph2_onTrk_y[clusterKey].push_back(hit->globalPosition().y());
4020             ph2_onTrk_z[clusterKey].push_back(hit->globalPosition().z());
4021             ph2_onTrk_xx[clusterKey].push_back(hit->globalPositionError().cxx());
4022             ph2_onTrk_xy[clusterKey].push_back(hit->globalPositionError().cyx());
4023             ph2_onTrk_yy[clusterKey].push_back(hit->globalPositionError().cyy());
4024             ph2_onTrk_yz[clusterKey].push_back(hit->globalPositionError().czy());
4025             ph2_onTrk_zz[clusterKey].push_back(hit->globalPositionError().czz());
4026             ph2_onTrk_zx[clusterKey].push_back(hit->globalPositionError().czx());
4027           } else {
4028             str_trkIdx[clusterKey].push_back(iTrack);
4029             str_onTrk_x[clusterKey].push_back(hit->globalPosition().x());
4030             str_onTrk_y[clusterKey].push_back(hit->globalPosition().y());
4031             str_onTrk_z[clusterKey].push_back(hit->globalPosition().z());
4032             str_onTrk_xx[clusterKey].push_back(hit->globalPositionError().cxx());
4033             str_onTrk_xy[clusterKey].push_back(hit->globalPositionError().cyx());
4034             str_onTrk_yy[clusterKey].push_back(hit->globalPositionError().cyy());
4035             str_onTrk_yz[clusterKey].push_back(hit->globalPositionError().czy());
4036             str_onTrk_zz[clusterKey].push_back(hit->globalPositionError().czz());
4037             str_onTrk_zx[clusterKey].push_back(hit->globalPositionError().czx());
4038           }
4039         }
4040 
4041         hitIdx.push_back(clusterKey);
4042         if (clusterRef.isPixel()) {
4043           hitType.push_back(static_cast<int>(HitType::Pixel));
4044         } else if (clusterRef.isPhase2()) {
4045           hitType.push_back(static_cast<int>(HitType::Phase2OT));
4046         } else {
4047           hitType.push_back(static_cast<int>(HitType::Strip));
4048         }
4049       } else {
4050         LogTrace("TrackingNtuple") << " - invalid hit";
4051 
4052         hitIdx.push_back(inv_isBarrel.size());
4053         hitType.push_back(static_cast<int>(HitType::Invalid));
4054 
4055         inv_isBarrel.push_back(hitId.subdetId() == 1);
4056         if (includeStripHits_)
4057           inv_detId.push_back(tracker, tTopo, hitId);
4058         else
4059           inv_detId_phase2.push_back(tracker, tTopo, hitId);
4060         inv_type.push_back(hit->getType());
4061       }
4062     }
4063 
4064     trk_hitIdx.push_back(hitIdx);
4065     trk_hitType.push_back(hitType);
4066   }
4067 }
4068 
4069 void TrackingNtuple::fillCandidates(const edm::Handle<TrackCandidateCollection>& candsHandle,
4070                                     int algo,
4071                                     const TrackingParticleRefVector& tpCollection,
4072                                     const TrackingParticleRefKeyToIndex& tpKeyToIndex,
4073                                     const TrackingParticleRefKeyToCount& tpKeyToClusterCount,
4074                                     const MagneticField& mf,
4075                                     const reco::BeamSpot& bs,
4076                                     const reco::VertexCollection& vertices,
4077                                     const reco::TrackToTrackingParticleAssociator& associatorByHits,
4078                                     const ClusterTPAssociation& clusterToTPMap,
4079                                     const TrackerGeometry& tracker,
4080                                     const TrackerTopology& tTopo,
4081                                     const std::set<edm::ProductID>& hitProductIds,
4082                                     const std::map<edm::ProductID, size_t>& seedCollToOffset) {
4083   reco::RecoToSimCollectionTCandidate recSimColl = associatorByHits.associateRecoToSim(candsHandle, tpCollection);
4084 
4085   TSCBLBuilderNoMaterial tscblBuilder;
4086 
4087   auto const& cands = *candsHandle;
4088   for (size_t iCand = 0; iCand < cands.size(); ++iCand) {
4089     const auto& aCand = cands[iCand];
4090     const edm::Ref<TrackCandidateCollection> aCandRef(candsHandle, iCand);
4091 
4092     //get parameters and errors from the candidate state
4093     auto const& pState = aCand.trajectoryStateOnDet();
4094     TrajectoryStateOnSurface state =
4095         trajectoryStateTransform::transientState(pState, &(tracker.idToDet(pState.detId())->surface()), &mf);
4096     TrajectoryStateClosestToBeamLine tbStateAtPCA = tscblBuilder(*state.freeState(), bs);
4097     if (!tbStateAtPCA.isValid()) {
4098       edm::LogVerbatim("TrackBuilding") << "TrajectoryStateClosestToBeamLine not valid";
4099     }
4100 
4101     auto const& stateAtPCA = tbStateAtPCA.trackStateAtPCA();
4102     auto v0 = stateAtPCA.position();
4103     auto p = stateAtPCA.momentum();
4104     math::XYZPoint pos(v0.x(), v0.y(), v0.z());
4105     math::XYZVector mom(p.x(), p.y(), p.z());
4106 
4107     //pseduo track for access to easy methods
4108     static const reco::Track::CovarianceMatrix dummyCov = AlgebraicMatrixID();
4109     reco::Track trk(
4110         0, 0, pos, mom, stateAtPCA.charge(), tbStateAtPCA.isValid() ? stateAtPCA.curvilinearError().matrix() : dummyCov);
4111 
4112     const auto& tkParam = trk.parameters();
4113     auto tkCov = trk.covariance();
4114     tkCov.Invert();
4115 
4116     // Standard track-TP matching
4117     int nSimHits = 0;
4118     bool isSimMatched = false;
4119     std::vector<int> tpIdx;
4120     std::vector<float> sharedFraction;
4121     std::vector<float> tpChi2;
4122     auto foundTPs = recSimColl.find(aCandRef);
4123     if (foundTPs != recSimColl.end()) {
4124       if (!foundTPs->val.empty()) {
4125         nSimHits = foundTPs->val[0].first->numberOfTrackerHits();
4126         isSimMatched = true;
4127       }
4128       for (const auto& tpQuality : foundTPs->val) {
4129         tpIdx.push_back(tpKeyToIndex.at(tpQuality.first.key()));
4130         sharedFraction.push_back(tpQuality.second);
4131         tpChi2.push_back(track_associator::trackAssociationChi2(tkParam, tkCov, *(tpCollection[tpIdx.back()]), mf, bs));
4132       }
4133     }
4134 
4135     // Search for a best-matching TrackingParticle for a track
4136     const auto clusters = track_associator::hitsToClusterRefs(aCand.recHits().begin(), aCand.recHits().end());
4137     const int nClusters = clusters.size();
4138 
4139     const auto bestKeyCount = findBestMatchingTrackingParticle(aCand.recHits(), clusterToTPMap, tpKeyToIndex);
4140     const float bestCountF = bestKeyCount.countClusters;
4141     const float bestShareFrac = bestCountF / nClusters;
4142     float bestShareFracSimDenom = 0;
4143     float bestShareFracSimClusterDenom = 0;
4144     float bestChi2 = -1;
4145     if (bestKeyCount.key >= 0) {
4146       bestShareFracSimDenom = bestCountF / tpCollection[tpKeyToIndex.at(bestKeyCount.key)]->numberOfTrackerHits();
4147       bestShareFracSimClusterDenom = bestCountF / tpKeyToClusterCount.at(bestKeyCount.key);
4148       bestChi2 = track_associator::trackAssociationChi2(
4149           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestKeyCount.key)]), mf, bs);
4150     }
4151     // Another way starting from the first hit of the track
4152     const auto bestFirstHitKeyCount =
4153         findMatchingTrackingParticleFromFirstHit(aCand.recHits(), clusterToTPMap, tpKeyToIndex);
4154     const float bestFirstCountF = bestFirstHitKeyCount.countClusters;
4155     const float bestFirstHitShareFrac = bestFirstCountF / nClusters;
4156     float bestFirstHitShareFracSimDenom = 0;
4157     float bestFirstHitShareFracSimClusterDenom = 0;
4158     float bestFirstHitChi2 = -1;
4159     if (bestFirstHitKeyCount.key >= 0) {
4160       bestFirstHitShareFracSimDenom =
4161           bestFirstCountF / tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]->numberOfTrackerHits();
4162       bestFirstHitShareFracSimClusterDenom = bestFirstCountF / tpKeyToClusterCount.at(bestFirstHitKeyCount.key);
4163       bestFirstHitChi2 = track_associator::trackAssociationChi2(
4164           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]), mf, bs);
4165     }
4166 
4167     auto iglobCand = tcand_pca_valid.size();  //global cand index
4168     tcand_pca_valid.push_back(tbStateAtPCA.isValid());
4169     tcand_pca_px.push_back(trk.px());
4170     tcand_pca_py.push_back(trk.py());
4171     tcand_pca_pz.push_back(trk.pz());
4172     tcand_pca_pt.push_back(trk.pt());
4173     tcand_pca_eta.push_back(trk.eta());
4174     tcand_pca_phi.push_back(trk.phi());
4175     tcand_pca_dxy.push_back(trk.dxy());
4176     tcand_pca_dz.push_back(trk.dz());
4177     tcand_pca_ptErr.push_back(trk.ptError());
4178     tcand_pca_etaErr.push_back(trk.etaError());
4179     tcand_pca_lambdaErr.push_back(trk.lambdaError());
4180     tcand_pca_phiErr.push_back(trk.phiError());
4181     tcand_pca_dxyErr.push_back(trk.dxyError());
4182     tcand_pca_dzErr.push_back(trk.dzError());
4183     tcand_px.push_back(state.globalMomentum().x());
4184     tcand_py.push_back(state.globalMomentum().y());
4185     tcand_pz.push_back(state.globalMomentum().z());
4186     tcand_pt.push_back(state.globalMomentum().perp());
4187     tcand_x.push_back(state.globalPosition().x());
4188     tcand_y.push_back(state.globalPosition().y());
4189     tcand_z.push_back(state.globalPosition().z());
4190 
4191     auto const& pStateCov = state.curvilinearError().matrix();
4192     tcand_qbpErr.push_back(sqrt(pStateCov(0, 0)));
4193     tcand_lambdaErr.push_back(sqrt(pStateCov(1, 1)));
4194     tcand_phiErr.push_back(sqrt(pStateCov(2, 2)));
4195     tcand_xtErr.push_back(sqrt(pStateCov(3, 3)));
4196     tcand_ytErr.push_back(sqrt(pStateCov(4, 4)));
4197 
4198     int ndof = -5;
4199     int nValid = 0;
4200     int nPixel = 0;
4201     int nStrip = 0;
4202     for (auto const& hit : aCand.recHits()) {
4203       if (hit.isValid()) {
4204         ndof += hit.dimension();
4205         nValid++;
4206 
4207         auto const subdet = hit.geographicalId().subdetId();
4208         if (subdet == PixelSubdetector::PixelBarrel || subdet == PixelSubdetector::PixelEndcap)
4209           nPixel++;
4210         else
4211           nStrip++;
4212       }
4213     }
4214     tcand_ndof.push_back(ndof);
4215     tcand_q.push_back(trk.charge());
4216     tcand_nValid.push_back(nValid);
4217     tcand_nPixel.push_back(nPixel);
4218     tcand_nStrip.push_back(nStrip);
4219     tcand_nCluster.push_back(nClusters);
4220     tcand_algo.push_back(algo);
4221     tcand_stopReason.push_back(aCand.stopReason());
4222     if (includeSeeds_) {
4223       auto offset = seedCollToOffset.find(aCand.seedRef().id());
4224       if (offset == seedCollToOffset.end()) {
4225         throw cms::Exception("Configuration")
4226             << "Track candidate refers to seed collection " << aCand.seedRef().id()
4227             << ", but that seed collection is not given as an input. The following collections were given as an input "
4228             << make_ProductIDMapPrinter(seedCollToOffset);
4229       }
4230 
4231       const auto seedIndex = offset->second + aCand.seedRef().key();
4232       tcand_seedIdx.push_back(seedIndex);
4233       if (see_tcandIdx[seedIndex] != -1) {
4234         throw cms::Exception("LogicError")
4235             << "Track cand index has already been set for seed " << seedIndex << " to " << see_tcandIdx[seedIndex]
4236             << "; was trying to set it to " << iglobCand << " current " << iCand;
4237       }
4238       see_tcandIdx[seedIndex] = iglobCand;
4239     }
4240     tcand_vtxIdx.push_back(-1);  // to be set correctly in fillVertices
4241     if (includeTrackingParticles_) {
4242       tcand_simTrkIdx.push_back(tpIdx);
4243       tcand_simTrkShareFrac.push_back(sharedFraction);
4244       tcand_simTrkNChi2.push_back(tpChi2);
4245       tcand_bestSimTrkIdx.push_back(bestKeyCount.key >= 0 ? tpKeyToIndex.at(bestKeyCount.key) : -1);
4246       tcand_bestFromFirstHitSimTrkIdx.push_back(
4247           bestFirstHitKeyCount.key >= 0 ? tpKeyToIndex.at(bestFirstHitKeyCount.key) : -1);
4248     } else {
4249       tcand_isTrue.push_back(!tpIdx.empty());
4250     }
4251     tcand_bestSimTrkShareFrac.push_back(bestShareFrac);
4252     tcand_bestSimTrkShareFracSimDenom.push_back(bestShareFracSimDenom);
4253     tcand_bestSimTrkShareFracSimClusterDenom.push_back(bestShareFracSimClusterDenom);
4254     tcand_bestSimTrkNChi2.push_back(bestChi2);
4255     tcand_bestFromFirstHitSimTrkShareFrac.push_back(bestFirstHitShareFrac);
4256     tcand_bestFromFirstHitSimTrkShareFracSimDenom.push_back(bestFirstHitShareFracSimDenom);
4257     tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom.push_back(bestFirstHitShareFracSimClusterDenom);
4258     tcand_bestFromFirstHitSimTrkNChi2.push_back(bestFirstHitChi2);
4259 
4260     LogTrace("TrackingNtuple") << "Track cand #" << iCand << " glob " << iglobCand << " with q=" << trk.charge()
4261                                << ", pT=" << trk.pt() << " GeV, eta: " << trk.eta() << ", phi: " << trk.phi()
4262                                << ", nValid=" << nValid << " seed#=" << aCand.seedRef().key()
4263                                << " simMatch=" << isSimMatched << " nSimHits=" << nSimHits
4264                                << " sharedFraction=" << (sharedFraction.empty() ? -1 : sharedFraction[0])
4265                                << " tpIdx=" << (tpIdx.empty() ? -1 : tpIdx[0]);
4266     std::vector<int> hitIdx;
4267     std::vector<int> hitType;
4268 
4269     for (auto i = aCand.recHits().begin(); i != aCand.recHits().end(); i++) {
4270       const TrackingRecHit* hit = &*i;
4271       DetId hitId = hit->geographicalId();
4272       LogTrace("TrackingNtuple") << "hit #" << std::distance(aCand.recHits().begin(), i)
4273                                  << " subdet=" << hitId.subdetId();
4274       if (hitId.det() != DetId::Tracker)
4275         continue;
4276 
4277       LogTrace("TrackingNtuple") << " " << subdetstring(hitId.subdetId()) << " " << tTopo.layer(hitId);
4278 
4279       if (hit->isValid()) {
4280         //ugly... but works
4281         const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*hit);
4282         const auto& clusterRef = bhit->firstClusterRef();
4283         unsigned int clusterKey;
4284         if (clusterRef.isPixel()) {
4285           clusterKey = clusterRef.cluster_pixel().key();
4286         } else if (clusterRef.isPhase2()) {
4287           clusterKey = clusterRef.cluster_phase2OT().key();
4288         } else {
4289           clusterKey = clusterRef.cluster_strip().key();
4290         }
4291 
4292         LogTrace("TrackingNtuple") << " id: " << hitId.rawId() << " - globalPos =" << hit->globalPosition()
4293                                    << " cluster=" << clusterKey << " clusterRef ID=" << clusterRef.id()
4294                                    << " eta,phi: " << hit->globalPosition().eta() << "," << hit->globalPosition().phi();
4295         if (includeAllHits_) {
4296           checkProductID(hitProductIds, clusterRef.id(), "track");
4297           if (clusterRef.isPixel()) {
4298             pix_tcandIdx[clusterKey].push_back(iglobCand);
4299           } else if (clusterRef.isPhase2()) {
4300             ph2_tcandIdx[clusterKey].push_back(iglobCand);
4301           } else {
4302             str_tcandIdx[clusterKey].push_back(iglobCand);
4303           }
4304         }
4305 
4306         hitIdx.push_back(clusterKey);
4307         if (clusterRef.isPixel()) {
4308           hitType.push_back(static_cast<int>(HitType::Pixel));
4309         } else if (clusterRef.isPhase2()) {
4310           hitType.push_back(static_cast<int>(HitType::Phase2OT));
4311         } else {
4312           hitType.push_back(static_cast<int>(HitType::Strip));
4313         }
4314       }
4315     }
4316 
4317     tcand_hitIdx.push_back(hitIdx);
4318     tcand_hitType.push_back(hitType);
4319   }
4320 }
4321 
4322 void TrackingNtuple::fillTrackingParticles(const edm::Event& iEvent,
4323                                            const edm::EventSetup& iSetup,
4324                                            const edm::RefToBaseVector<reco::Track>& tracks,
4325                                            const reco::BeamSpot& bs,
4326                                            const TrackingParticleRefVector& tpCollection,
4327                                            const TrackingVertexRefKeyToIndex& tvKeyToIndex,
4328                                            const reco::TrackToTrackingParticleAssociator& associatorByHits,
4329                                            const std::vector<TPHitIndex>& tpHitList,
4330                                            const TrackingParticleRefKeyToCount& tpKeyToClusterCount) {
4331   // Number of 3D layers for TPs
4332   edm::Handle<edm::ValueMap<unsigned int>> tpNLayersH;
4333   iEvent.getByToken(tpNLayersToken_, tpNLayersH);
4334   const auto& nLayers_tPCeff = *tpNLayersH;
4335 
4336   iEvent.getByToken(tpNPixelLayersToken_, tpNLayersH);
4337   const auto& nPixelLayers_tPCeff = *tpNLayersH;
4338 
4339   iEvent.getByToken(tpNStripStereoLayersToken_, tpNLayersH);
4340   const auto& nStripMonoAndStereoLayers_tPCeff = *tpNLayersH;
4341 
4342   reco::SimToRecoCollection simRecColl = associatorByHits.associateSimToReco(tracks, tpCollection);
4343 
4344   for (const TrackingParticleRef& tp : tpCollection) {
4345     LogTrace("TrackingNtuple") << "tracking particle pt=" << tp->pt() << " eta=" << tp->eta() << " phi=" << tp->phi();
4346     bool isRecoMatched = false;
4347     std::vector<int> tkIdx;
4348     std::vector<float> sharedFraction;
4349     auto foundTracks = simRecColl.find(tp);
4350     if (foundTracks != simRecColl.end()) {
4351       isRecoMatched = true;
4352       for (const auto& trackQuality : foundTracks->val) {
4353         sharedFraction.push_back(trackQuality.second);
4354         tkIdx.push_back(trackQuality.first.key());
4355       }
4356     }
4357 
4358     sim_genPdgIds.emplace_back();
4359     for (const auto& genRef : tp->genParticles()) {
4360       if (genRef.isNonnull())
4361         sim_genPdgIds.back().push_back(genRef->pdgId());
4362     }
4363 
4364     bool isFromBHadron = false;
4365     // Logic is similar to SimTracker/TrackHistory
4366     if (tracer_.evaluate(tp)) {  // ignore TP if history can not be traced
4367       // following is from TrackClassifier::processesAtGenerator()
4368       HistoryBase::RecoGenParticleTrail const& recoGenParticleTrail = tracer_.recoGenParticleTrail();
4369       for (const auto& particle : recoGenParticleTrail) {
4370         HepPDT::ParticleID particleID(particle->pdgId());
4371         if (particleID.hasBottom()) {
4372           isFromBHadron = true;
4373           break;
4374         }
4375       }
4376     }
4377 
4378     LogTrace("TrackingNtuple") << "matched to tracks = " << make_VectorPrinter(tkIdx)
4379                                << " isRecoMatched=" << isRecoMatched;
4380     sim_event.push_back(tp->eventId().event());
4381     sim_bunchCrossing.push_back(tp->eventId().bunchCrossing());
4382     sim_pdgId.push_back(tp->pdgId());
4383     sim_isFromBHadron.push_back(isFromBHadron);
4384     sim_px.push_back(tp->px());
4385     sim_py.push_back(tp->py());
4386     sim_pz.push_back(tp->pz());
4387     sim_pt.push_back(tp->pt());
4388     sim_eta.push_back(tp->eta());
4389     sim_phi.push_back(tp->phi());
4390     sim_q.push_back(tp->charge());
4391     sim_trkIdx.push_back(tkIdx);
4392     sim_trkShareFrac.push_back(sharedFraction);
4393     sim_parentVtxIdx.push_back(tvKeyToIndex.at(tp->parentVertex().key()));
4394     std::vector<int> decayIdx;
4395     for (const auto& v : tp->decayVertices())
4396       decayIdx.push_back(tvKeyToIndex.at(v.key()));
4397     sim_decayVtxIdx.push_back(decayIdx);
4398 
4399     //Calcualte the impact parameters w.r.t. PCA
4400     TrackingParticle::Vector momentum = parametersDefiner_.momentum(iEvent, iSetup, tp);
4401     TrackingParticle::Point vertex = parametersDefiner_.vertex(iEvent, iSetup, tp);
4402     auto dxySim = TrackingParticleIP::dxy(vertex, momentum, bs.position());
4403     auto dzSim = TrackingParticleIP::dz(vertex, momentum, bs.position());
4404     const double lambdaSim = M_PI / 2 - momentum.theta();
4405     sim_pca_pt.push_back(std::sqrt(momentum.perp2()));
4406     sim_pca_eta.push_back(momentum.Eta());
4407     sim_pca_lambda.push_back(lambdaSim);
4408     sim_pca_cotTheta.push_back(1 / tan(M_PI * 0.5 - lambdaSim));
4409     sim_pca_phi.push_back(momentum.phi());
4410     sim_pca_dxy.push_back(dxySim);
4411     sim_pca_dz.push_back(dzSim);
4412 
4413     std::vector<int> hitIdx;
4414     int nPixel = 0, nStrip = 0;
4415     auto rangeHit = std::equal_range(tpHitList.begin(), tpHitList.end(), TPHitIndex(tp.key()), tpHitIndexListLess);
4416     for (auto ip = rangeHit.first; ip != rangeHit.second; ++ip) {
4417       auto type = HitType::Unknown;
4418       if (!simhit_hitType[ip->simHitIdx].empty())
4419         type = static_cast<HitType>(simhit_hitType[ip->simHitIdx][0]);
4420       LogTrace("TrackingNtuple") << "simhit=" << ip->simHitIdx << " type=" << static_cast<int>(type);
4421       hitIdx.push_back(ip->simHitIdx);
4422       const auto detid = DetId(includeStripHits_ ? simhit_detId[ip->simHitIdx] : simhit_detId_phase2[ip->simHitIdx]);
4423       if (detid.det() != DetId::Tracker) {
4424         throw cms::Exception("LogicError") << "Encountered SimHit for TP " << tp.key() << " with DetId "
4425                                            << detid.rawId() << " whose det() is not Tracker but " << detid.det();
4426       }
4427       const auto subdet = detid.subdetId();
4428       switch (subdet) {
4429         case PixelSubdetector::PixelBarrel:
4430         case PixelSubdetector::PixelEndcap:
4431           ++nPixel;
4432           break;
4433         case StripSubdetector::TIB:
4434         case StripSubdetector::TID:
4435         case StripSubdetector::TOB:
4436         case StripSubdetector::TEC:
4437           ++nStrip;
4438           break;
4439         default:
4440           throw cms::Exception("LogicError") << "Encountered SimHit for TP " << tp.key() << " with DetId "
4441                                              << detid.rawId() << " whose subdet is not recognized, is " << subdet;
4442       };
4443     }
4444     sim_nValid.push_back(hitIdx.size());
4445     sim_nPixel.push_back(nPixel);
4446     sim_nStrip.push_back(nStrip);
4447 
4448     const auto nSimLayers = nLayers_tPCeff[tp];
4449     const auto nSimPixelLayers = nPixelLayers_tPCeff[tp];
4450     const auto nSimStripMonoAndStereoLayers = nStripMonoAndStereoLayers_tPCeff[tp];
4451     sim_nLay.push_back(nSimLayers);
4452     sim_nPixelLay.push_back(nSimPixelLayers);
4453     sim_n3DLay.push_back(nSimPixelLayers + nSimStripMonoAndStereoLayers);
4454 
4455     sim_nTrackerHits.push_back(tp->numberOfTrackerHits());
4456     auto found = tpKeyToClusterCount.find(tp.key());
4457     sim_nRecoClusters.push_back(found != cend(tpKeyToClusterCount) ? found->second : 0);
4458 
4459     sim_simHitIdx.push_back(hitIdx);
4460   }
4461 }
4462 
4463 // called from fillSeeds
4464 void TrackingNtuple::fillTrackingParticlesForSeeds(const TrackingParticleRefVector& tpCollection,
4465                                                    const reco::SimToRecoCollection& simRecColl,
4466                                                    const TrackingParticleRefKeyToIndex& tpKeyToIndex,
4467                                                    const unsigned int seedOffset) {
4468   if (sim_seedIdx.empty())  // first call
4469     sim_seedIdx.resize(tpCollection.size());
4470 
4471   for (const auto& keyVal : simRecColl) {
4472     const auto& tpRef = keyVal.key;
4473     auto found = tpKeyToIndex.find(tpRef.key());
4474     if (found == tpKeyToIndex.end())
4475       throw cms::Exception("Assert") << __FILE__ << ":" << __LINE__ << " fillTrackingParticlesForSeeds: tpRef.key() "
4476                                      << tpRef.key() << " not found from tpKeyToIndex. tpKeyToIndex size "
4477                                      << tpKeyToIndex.size();
4478     const auto tpIndex = found->second;
4479     for (const auto& pair : keyVal.val) {
4480       const auto& seedRef = pair.first->seedRef();
4481       sim_seedIdx[tpIndex].push_back(seedOffset + seedRef.key());
4482     }
4483   }
4484 }
4485 
4486 void TrackingNtuple::fillVertices(const reco::VertexCollection& vertices,
4487                                   const edm::RefToBaseVector<reco::Track>& tracks) {
4488   for (size_t iVertex = 0, size = vertices.size(); iVertex < size; ++iVertex) {
4489     const reco::Vertex& vertex = vertices[iVertex];
4490     vtx_x.push_back(vertex.x());
4491     vtx_y.push_back(vertex.y());
4492     vtx_z.push_back(vertex.z());
4493     vtx_xErr.push_back(vertex.xError());
4494     vtx_yErr.push_back(vertex.yError());
4495     vtx_zErr.push_back(vertex.zError());
4496     vtx_chi2.push_back(vertex.chi2());
4497     vtx_ndof.push_back(vertex.ndof());
4498     vtx_fake.push_back(vertex.isFake());
4499     vtx_valid.push_back(vertex.isValid());
4500 
4501     std::vector<int> trkIdx;
4502     for (auto iTrack = vertex.tracks_begin(); iTrack != vertex.tracks_end(); ++iTrack) {
4503       // Ignore link if vertex was fit from a track collection different from the input
4504       if (iTrack->id() != tracks.id())
4505         continue;
4506 
4507       trkIdx.push_back(iTrack->key());
4508 
4509       if (trk_vtxIdx[iTrack->key()] != -1) {
4510         throw cms::Exception("LogicError") << "Vertex index has already been set for track " << iTrack->key() << " to "
4511                                            << trk_vtxIdx[iTrack->key()] << "; was trying to set it to " << iVertex;
4512       }
4513       trk_vtxIdx[iTrack->key()] = iVertex;
4514     }
4515     vtx_trkIdx.push_back(trkIdx);
4516   }
4517 }
4518 
4519 void TrackingNtuple::fillTrackingVertices(const TrackingVertexRefVector& trackingVertices,
4520                                           const TrackingParticleRefKeyToIndex& tpKeyToIndex) {
4521   int current_event = -1;
4522   for (const auto& ref : trackingVertices) {
4523     const TrackingVertex v = *ref;
4524     if (v.eventId().event() != current_event) {
4525       // next PV
4526       current_event = v.eventId().event();
4527       simpv_idx.push_back(simvtx_x.size());
4528     }
4529 
4530     unsigned int processType = std::numeric_limits<unsigned int>::max();
4531     if (!v.g4Vertices().empty()) {
4532       processType = v.g4Vertices()[0].processType();
4533     }
4534 
4535     simvtx_event.push_back(v.eventId().event());
4536     simvtx_bunchCrossing.push_back(v.eventId().bunchCrossing());
4537     simvtx_processType.push_back(processType);
4538