RecoTrack/plugins/TrackingNtuple.cc

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