Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:24 

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 "Validation/RecoTrack/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 
1296   ////////////////////
1297   // invalid (missing/inactive/etc) hits
1298   // (first) index runs through hits
1299   std::vector<short> inv_isBarrel;
1300   DetIdAll inv_detId;
1301   DetIdAllPhase2 inv_detId_phase2;
1302   std::vector<unsigned short> inv_type;
1303   ////////////////////
1304   // sim hits
1305   // (first) index runs through hits
1306   DetIdAll simhit_detId;
1307   DetIdAllPhase2 simhit_detId_phase2;
1308   std::vector<float> simhit_x;
1309   std::vector<float> simhit_y;
1310   std::vector<float> simhit_z;
1311   std::vector<float> simhit_px;
1312   std::vector<float> simhit_py;
1313   std::vector<float> simhit_pz;
1314   std::vector<int> simhit_particle;
1315   std::vector<short> simhit_process;
1316   std::vector<float> simhit_eloss;
1317   std::vector<float> simhit_tof;
1318   //std::vector<unsigned int> simhit_simTrackId; // can be useful for debugging, but not much of general interest
1319   std::vector<int> simhit_simTrkIdx;
1320   std::vector<std::vector<int>> simhit_hitIdx;   // second index runs through induced reco hits
1321   std::vector<std::vector<int>> simhit_hitType;  // second index runs through induced reco hits
1322   ////////////////////
1323   // beam spot
1324   float bsp_x;
1325   float bsp_y;
1326   float bsp_z;
1327   float bsp_sigmax;
1328   float bsp_sigmay;
1329   float bsp_sigmaz;
1330   ////////////////////
1331   // seeds
1332   // (first) index runs through seeds
1333   std::vector<short> see_fitok;
1334   std::vector<float> see_px;
1335   std::vector<float> see_py;
1336   std::vector<float> see_pz;
1337   std::vector<float> see_pt;
1338   std::vector<float> see_eta;
1339   std::vector<float> see_phi;
1340   std::vector<float> see_dxy;
1341   std::vector<float> see_dz;
1342   std::vector<float> see_ptErr;
1343   std::vector<float> see_etaErr;
1344   std::vector<float> see_phiErr;
1345   std::vector<float> see_dxyErr;
1346   std::vector<float> see_dzErr;
1347   std::vector<float> see_chi2;
1348   std::vector<float> see_statePt;
1349   std::vector<float> see_stateTrajX;
1350   std::vector<float> see_stateTrajY;
1351   std::vector<float> see_stateTrajPx;
1352   std::vector<float> see_stateTrajPy;
1353   std::vector<float> see_stateTrajPz;
1354   std::vector<float> see_stateTrajGlbX;
1355   std::vector<float> see_stateTrajGlbY;
1356   std::vector<float> see_stateTrajGlbZ;
1357   std::vector<float> see_stateTrajGlbPx;
1358   std::vector<float> see_stateTrajGlbPy;
1359   std::vector<float> see_stateTrajGlbPz;
1360   std::vector<std::vector<float>> see_stateCurvCov;
1361   std::vector<int> see_q;
1362   std::vector<unsigned int> see_nValid;
1363   std::vector<unsigned int> see_nPixel;
1364   std::vector<unsigned int> see_nGlued;
1365   std::vector<unsigned int> see_nStrip;
1366   std::vector<unsigned int> see_nPhase2OT;
1367   std::vector<unsigned int> see_nCluster;
1368   std::vector<unsigned int> see_algo;
1369   std::vector<unsigned short> see_stopReason;
1370   std::vector<unsigned short> see_nCands;
1371   std::vector<int> see_trkIdx;
1372   std::vector<int> see_tcandIdx;
1373   std::vector<short> see_isTrue;
1374   std::vector<int> see_bestSimTrkIdx;
1375   std::vector<float> see_bestSimTrkShareFrac;
1376   std::vector<int> see_bestFromFirstHitSimTrkIdx;
1377   std::vector<float> see_bestFromFirstHitSimTrkShareFrac;
1378   std::vector<std::vector<float>> see_simTrkShareFrac;  // second index runs through matched TrackingParticles
1379   std::vector<std::vector<int>> see_simTrkIdx;          // second index runs through matched TrackingParticles
1380   std::vector<std::vector<int>> see_hitIdx;             // second index runs through hits
1381   std::vector<std::vector<int>> see_hitType;            // second index runs through hits
1382   //seed algo offset, index runs through iterations
1383   std::vector<unsigned int> see_offset;
1384 
1385   ////////////////////
1386   // Vertices
1387   // (first) index runs through vertices
1388   std::vector<float> vtx_x;
1389   std::vector<float> vtx_y;
1390   std::vector<float> vtx_z;
1391   std::vector<float> vtx_xErr;
1392   std::vector<float> vtx_yErr;
1393   std::vector<float> vtx_zErr;
1394   std::vector<float> vtx_ndof;
1395   std::vector<float> vtx_chi2;
1396   std::vector<short> vtx_fake;
1397   std::vector<short> vtx_valid;
1398   std::vector<std::vector<int>> vtx_trkIdx;  // second index runs through tracks used in the vertex fit
1399 
1400   ////////////////////
1401   // Tracking vertices
1402   // (first) index runs through TrackingVertices
1403   std::vector<int> simvtx_event;
1404   std::vector<int> simvtx_bunchCrossing;
1405   std::vector<unsigned int> simvtx_processType;  // only from first SimVertex of TrackingVertex
1406   std::vector<float> simvtx_x;
1407   std::vector<float> simvtx_y;
1408   std::vector<float> simvtx_z;
1409   std::vector<std::vector<int>> simvtx_sourceSimIdx;    // second index runs through source TrackingParticles
1410   std::vector<std::vector<int>> simvtx_daughterSimIdx;  // second index runs through daughter TrackingParticles
1411   std::vector<int> simpv_idx;
1412 };
1413 
1414 //
1415 // constructors and destructor
1416 //
1417 TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig)
1418     : mfToken_(esConsumes()),
1419       tTopoToken_(esConsumes()),
1420       tGeomToken_(esConsumes()),
1421       trackToken_(consumes<edm::View<reco::Track>>(iConfig.getUntrackedParameter<edm::InputTag>("tracks"))),
1422       clusterTPMapToken_(consumes<ClusterTPAssociation>(iConfig.getUntrackedParameter<edm::InputTag>("clusterTPMap"))),
1423       simHitTPMapToken_(consumes<SimHitTPAssociationProducer::SimHitTPAssociationList>(
1424           iConfig.getUntrackedParameter<edm::InputTag>("simHitTPMap"))),
1425       trackAssociatorToken_(consumes<reco::TrackToTrackingParticleAssociator>(
1426           iConfig.getUntrackedParameter<edm::InputTag>("trackAssociator"))),
1427       pixelSimLinkToken_(consumes<edm::DetSetVector<PixelDigiSimLink>>(
1428           iConfig.getUntrackedParameter<edm::InputTag>("pixelDigiSimLink"))),
1429       stripSimLinkToken_(consumes<edm::DetSetVector<StripDigiSimLink>>(
1430           iConfig.getUntrackedParameter<edm::InputTag>("stripDigiSimLink"))),
1431       siphase2OTSimLinksToken_(consumes<edm::DetSetVector<PixelDigiSimLink>>(
1432           iConfig.getUntrackedParameter<edm::InputTag>("phase2OTSimLink"))),
1433       includeStripHits_(!iConfig.getUntrackedParameter<edm::InputTag>("stripDigiSimLink").label().empty()),
1434       includePhase2OTHits_(!iConfig.getUntrackedParameter<edm::InputTag>("phase2OTSimLink").label().empty()),
1435       beamSpotToken_(consumes<reco::BeamSpot>(iConfig.getUntrackedParameter<edm::InputTag>("beamSpot"))),
1436       pixelRecHitToken_(
1437           consumes<SiPixelRecHitCollection>(iConfig.getUntrackedParameter<edm::InputTag>("pixelRecHits"))),
1438       stripRphiRecHitToken_(
1439           consumes<SiStripRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripRphiRecHits"))),
1440       stripStereoRecHitToken_(
1441           consumes<SiStripRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripStereoRecHits"))),
1442       stripMatchedRecHitToken_(consumes<SiStripMatchedRecHit2DCollection>(
1443           iConfig.getUntrackedParameter<edm::InputTag>("stripMatchedRecHits"))),
1444       phase2OTRecHitToken_(consumes<Phase2TrackerRecHit1DCollectionNew>(
1445           iConfig.getUntrackedParameter<edm::InputTag>("phase2OTRecHits"))),
1446       vertexToken_(consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertices"))),
1447       trackingVertexToken_(
1448           consumes<TrackingVertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trackingVertices"))),
1449       tpNLayersToken_(consumes<edm::ValueMap<unsigned int>>(
1450           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNlayers"))),
1451       tpNPixelLayersToken_(consumes<edm::ValueMap<unsigned int>>(
1452           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNpixellayers"))),
1453       tpNStripStereoLayersToken_(consumes<edm::ValueMap<unsigned int>>(
1454           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNstripstereolayers"))),
1455       includeSeeds_(iConfig.getUntrackedParameter<bool>("includeSeeds")),
1456       includeTrackCandidates_(iConfig.getUntrackedParameter<bool>("includeTrackCandidates")),
1457       addSeedCurvCov_(iConfig.getUntrackedParameter<bool>("addSeedCurvCov")),
1458       includeAllHits_(iConfig.getUntrackedParameter<bool>("includeAllHits")),
1459       includeOnTrackHitData_(iConfig.getUntrackedParameter<bool>("includeOnTrackHitData")),
1460       includeMVA_(iConfig.getUntrackedParameter<bool>("includeMVA")),
1461       includeTrackingParticles_(iConfig.getUntrackedParameter<bool>("includeTrackingParticles")),
1462       includeOOT_(iConfig.getUntrackedParameter<bool>("includeOOT")),
1463       keepEleSimHits_(iConfig.getUntrackedParameter<bool>("keepEleSimHits")),
1464       saveSimHitsP3_(iConfig.getUntrackedParameter<bool>("saveSimHitsP3")),
1465       simHitBySignificance_(iConfig.getUntrackedParameter<bool>("simHitBySignificance")),
1466       parametersDefiner_(iConfig.getUntrackedParameter<edm::InputTag>("beamSpot"), consumesCollector()) {
1467   if (includeSeeds_) {
1468     seedTokens_ =
1469         edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("seedTracks"),
1470                               [&](const edm::InputTag& tag) { return consumes<edm::View<reco::Track>>(tag); });
1471     seedStopInfoTokens_ =
1472         edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("trackCandidates"),
1473                               [&](const edm::InputTag& tag) { return consumes<std::vector<SeedStopInfo>>(tag); });
1474     if (seedTokens_.size() != seedStopInfoTokens_.size()) {
1475       throw cms::Exception("Configuration") << "Got " << seedTokens_.size() << " seed collections, but "
1476                                             << seedStopInfoTokens_.size() << " track candidate collections";
1477     }
1478   }
1479   if (includeTrackCandidates_)
1480     candidateTokens_ =
1481         edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("trackCandidates"),
1482                               [&](const edm::InputTag& tag) { return consumes<TrackCandidateCollection>(tag); });
1483 
1484   if (includeAllHits_) {
1485     if (includeStripHits_ && includePhase2OTHits_) {
1486       throw cms::Exception("Configuration")
1487           << "Both stripDigiSimLink and phase2OTSimLink are set, please set only either one (this information is used "
1488              "to infer if you're running phase0/1 or phase2 detector)";
1489     }
1490     if (!includeStripHits_ && !includePhase2OTHits_) {
1491       throw cms::Exception("Configuration")
1492           << "Neither stripDigiSimLink or phase2OTSimLink are set, please set either one.";
1493     }
1494 
1495     auto const& maskVPset = iConfig.getUntrackedParameterSetVector("clusterMasks");
1496     pixelUseMaskTokens_.reserve(maskVPset.size());
1497     stripUseMaskTokens_.reserve(maskVPset.size());
1498     ph2OTUseMaskTokens_.reserve(maskVPset.size());
1499     for (auto const& mask : maskVPset) {
1500       auto index = mask.getUntrackedParameter<unsigned int>("index");
1501       assert(index < 64);
1502       pixelUseMaskTokens_.emplace_back(index,
1503                                        consumes<PixelMaskContainer>(mask.getUntrackedParameter<edm::InputTag>("src")));
1504       if (includeStripHits_)
1505         stripUseMaskTokens_.emplace_back(
1506             index, consumes<StripMaskContainer>(mask.getUntrackedParameter<edm::InputTag>("src")));
1507       if (includePhase2OTHits_)
1508         ph2OTUseMaskTokens_.emplace_back(
1509             index, consumes<Phase2OTMaskContainer>(mask.getUntrackedParameter<edm::InputTag>("src")));
1510     }
1511   }
1512 
1513   const bool tpRef = iConfig.getUntrackedParameter<bool>("trackingParticlesRef");
1514   const auto tpTag = iConfig.getUntrackedParameter<edm::InputTag>("trackingParticles");
1515   if (tpRef) {
1516     trackingParticleRefToken_ = consumes<TrackingParticleRefVector>(tpTag);
1517   } else {
1518     trackingParticleToken_ = consumes<TrackingParticleCollection>(tpTag);
1519   }
1520 
1521   tracer_.depth(-2);  // as in SimTracker/TrackHistory/src/TrackClassifier.cc
1522 
1523   if (includeMVA_) {
1524     mvaQualityCollectionTokens_ = edm::vector_transform(
1525         iConfig.getUntrackedParameter<std::vector<std::string>>("trackMVAs"), [&](const std::string& tag) {
1526           return std::make_tuple(consumes<MVACollection>(edm::InputTag(tag, "MVAValues")),
1527                                  consumes<QualityMaskCollection>(edm::InputTag(tag, "QualityMasks")));
1528         });
1529   }
1530 
1531   usesResource(TFileService::kSharedResource);
1532   edm::Service<TFileService> fs;
1533   t = fs->make<TTree>("tree", "tree");
1534 
1535   t->Branch("event", &ev_event);
1536   t->Branch("lumi", &ev_lumi);
1537   t->Branch("run", &ev_run);
1538 
1539   //tracks
1540   t->Branch("trk_px", &trk_px);
1541   t->Branch("trk_py", &trk_py);
1542   t->Branch("trk_pz", &trk_pz);
1543   t->Branch("trk_pt", &trk_pt);
1544   t->Branch("trk_inner_px", &trk_inner_px);
1545   t->Branch("trk_inner_py", &trk_inner_py);
1546   t->Branch("trk_inner_pz", &trk_inner_pz);
1547   t->Branch("trk_inner_pt", &trk_inner_pt);
1548   t->Branch("trk_outer_px", &trk_outer_px);
1549   t->Branch("trk_outer_py", &trk_outer_py);
1550   t->Branch("trk_outer_pz", &trk_outer_pz);
1551   t->Branch("trk_outer_pt", &trk_outer_pt);
1552   t->Branch("trk_eta", &trk_eta);
1553   t->Branch("trk_lambda", &trk_lambda);
1554   t->Branch("trk_cotTheta", &trk_cotTheta);
1555   t->Branch("trk_phi", &trk_phi);
1556   t->Branch("trk_dxy", &trk_dxy);
1557   t->Branch("trk_dz", &trk_dz);
1558   t->Branch("trk_dxyPV", &trk_dxyPV);
1559   t->Branch("trk_dzPV", &trk_dzPV);
1560   t->Branch("trk_dxyClosestPV", &trk_dxyClosestPV);
1561   t->Branch("trk_dzClosestPV", &trk_dzClosestPV);
1562   t->Branch("trk_ptErr", &trk_ptErr);
1563   t->Branch("trk_etaErr", &trk_etaErr);
1564   t->Branch("trk_lambdaErr", &trk_lambdaErr);
1565   t->Branch("trk_phiErr", &trk_phiErr);
1566   t->Branch("trk_dxyErr", &trk_dxyErr);
1567   t->Branch("trk_dzErr", &trk_dzErr);
1568   t->Branch("trk_refpoint_x", &trk_refpoint_x);
1569   t->Branch("trk_refpoint_y", &trk_refpoint_y);
1570   t->Branch("trk_refpoint_z", &trk_refpoint_z);
1571   t->Branch("trk_nChi2", &trk_nChi2);
1572   t->Branch("trk_nChi2_1Dmod", &trk_nChi2_1Dmod);
1573   t->Branch("trk_ndof", &trk_ndof);
1574   if (includeMVA_) {
1575     trk_mvas.resize(mvaQualityCollectionTokens_.size());
1576     trk_qualityMasks.resize(mvaQualityCollectionTokens_.size());
1577     if (!trk_mvas.empty()) {
1578       t->Branch("trk_mva", &(trk_mvas[0]));
1579       t->Branch("trk_qualityMask", &(trk_qualityMasks[0]));
1580       for (size_t i = 1; i < trk_mvas.size(); ++i) {
1581         t->Branch(("trk_mva" + std::to_string(i + 1)).c_str(), &(trk_mvas[i]));
1582         t->Branch(("trk_qualityMask" + std::to_string(i + 1)).c_str(), &(trk_qualityMasks[i]));
1583       }
1584     }
1585   }
1586   t->Branch("trk_q", &trk_q);
1587   t->Branch("trk_nValid", &trk_nValid);
1588   t->Branch("trk_nLost", &trk_nLost);
1589   t->Branch("trk_nInactive", &trk_nInactive);
1590   t->Branch("trk_nPixel", &trk_nPixel);
1591   t->Branch("trk_nStrip", &trk_nStrip);
1592   t->Branch("trk_nOuterLost", &trk_nOuterLost);
1593   t->Branch("trk_nInnerLost", &trk_nInnerLost);
1594   t->Branch("trk_nOuterInactive", &trk_nOuterInactive);
1595   t->Branch("trk_nInnerInactive", &trk_nInnerInactive);
1596   t->Branch("trk_nPixelLay", &trk_nPixelLay);
1597   t->Branch("trk_nStripLay", &trk_nStripLay);
1598   t->Branch("trk_n3DLay", &trk_n3DLay);
1599   t->Branch("trk_nLostLay", &trk_nLostLay);
1600   t->Branch("trk_nCluster", &trk_nCluster);
1601   t->Branch("trk_algo", &trk_algo);
1602   t->Branch("trk_originalAlgo", &trk_originalAlgo);
1603   t->Branch("trk_algoMask", &trk_algoMask);
1604   t->Branch("trk_stopReason", &trk_stopReason);
1605   t->Branch("trk_isHP", &trk_isHP);
1606   if (includeSeeds_) {
1607     t->Branch("trk_seedIdx", &trk_seedIdx);
1608   }
1609   t->Branch("trk_vtxIdx", &trk_vtxIdx);
1610   if (includeTrackingParticles_) {
1611     t->Branch("trk_simTrkIdx", &trk_simTrkIdx);
1612     t->Branch("trk_simTrkShareFrac", &trk_simTrkShareFrac);
1613     t->Branch("trk_simTrkNChi2", &trk_simTrkNChi2);
1614     t->Branch("trk_bestSimTrkIdx", &trk_bestSimTrkIdx);
1615     t->Branch("trk_bestFromFirstHitSimTrkIdx", &trk_bestFromFirstHitSimTrkIdx);
1616   } else {
1617     t->Branch("trk_isTrue", &trk_isTrue);
1618   }
1619   t->Branch("trk_bestSimTrkShareFrac", &trk_bestSimTrkShareFrac);
1620   t->Branch("trk_bestSimTrkShareFracSimDenom", &trk_bestSimTrkShareFracSimDenom);
1621   t->Branch("trk_bestSimTrkShareFracSimClusterDenom", &trk_bestSimTrkShareFracSimClusterDenom);
1622   t->Branch("trk_bestSimTrkNChi2", &trk_bestSimTrkNChi2);
1623   t->Branch("trk_bestFromFirstHitSimTrkShareFrac", &trk_bestFromFirstHitSimTrkShareFrac);
1624   t->Branch("trk_bestFromFirstHitSimTrkShareFracSimDenom", &trk_bestFromFirstHitSimTrkShareFracSimDenom);
1625   t->Branch("trk_bestFromFirstHitSimTrkShareFracSimClusterDenom", &trk_bestFromFirstHitSimTrkShareFracSimClusterDenom);
1626   t->Branch("trk_bestFromFirstHitSimTrkNChi2", &trk_bestFromFirstHitSimTrkNChi2);
1627   if (includeAllHits_) {
1628     t->Branch("trk_hitIdx", &trk_hitIdx);
1629     t->Branch("trk_hitType", &trk_hitType);
1630   }
1631   if (includeTrackCandidates_) {
1632     t->Branch("tcand_pca_valid", &tcand_pca_valid);
1633     t->Branch("tcand_pca_px", &tcand_pca_px);
1634     t->Branch("tcand_pca_py", &tcand_pca_py);
1635     t->Branch("tcand_pca_pz", &tcand_pca_pz);
1636     t->Branch("tcand_pca_pt", &tcand_pca_pt);
1637     t->Branch("tcand_pca_eta", &tcand_pca_eta);
1638     t->Branch("tcand_pca_phi", &tcand_pca_phi);
1639     t->Branch("tcand_pca_dxy", &tcand_pca_dxy);
1640     t->Branch("tcand_pca_dz", &tcand_pca_dz);
1641     t->Branch("tcand_pca_ptErr", &tcand_pca_ptErr);
1642     t->Branch("tcand_pca_etaErr", &tcand_pca_etaErr);
1643     t->Branch("tcand_pca_lambdaErr", &tcand_pca_lambdaErr);
1644     t->Branch("tcand_pca_phiErr", &tcand_pca_phiErr);
1645     t->Branch("tcand_pca_dxyErr", &tcand_pca_dxyErr);
1646     t->Branch("tcand_pca_dzErr", &tcand_pca_dzErr);
1647     t->Branch("tcand_px", &tcand_px);
1648     t->Branch("tcand_py", &tcand_py);
1649     t->Branch("tcand_pz", &tcand_pz);
1650     t->Branch("tcand_pt", &tcand_pt);
1651     t->Branch("tcand_x", &tcand_x);
1652     t->Branch("tcand_y", &tcand_y);
1653     t->Branch("tcand_z", &tcand_z);
1654     t->Branch("tcand_qbpErr", &tcand_qbpErr);
1655     t->Branch("tcand_lambdaErr", &tcand_lambdaErr);
1656     t->Branch("tcand_phiErr", &tcand_phiErr);
1657     t->Branch("tcand_xtErr", &tcand_xtErr);
1658     t->Branch("tcand_ytErr", &tcand_ytErr);
1659     t->Branch("tcand_q", &tcand_q);
1660     t->Branch("tcand_ndof", &tcand_ndof);
1661     t->Branch("tcand_nValid", &tcand_nValid);
1662     t->Branch("tcand_nPixel", &tcand_nPixel);
1663     t->Branch("tcand_nStrip", &tcand_nStrip);
1664     t->Branch("tcand_nCluster", &tcand_nCluster);
1665     t->Branch("tcand_algo", &tcand_algo);
1666     t->Branch("tcand_stopReason", &tcand_stopReason);
1667     if (includeSeeds_) {
1668       t->Branch("tcand_seedIdx", &tcand_seedIdx);
1669     }
1670     t->Branch("tcand_vtxIdx", &tcand_vtxIdx);
1671     if (includeTrackingParticles_) {
1672       t->Branch("tcand_simTrkIdx", &tcand_simTrkIdx);
1673       t->Branch("tcand_simTrkShareFrac", &tcand_simTrkShareFrac);
1674       t->Branch("tcand_simTrkNChi2", &tcand_simTrkNChi2);
1675       t->Branch("tcand_bestSimTrkIdx", &tcand_bestSimTrkIdx);
1676       t->Branch("tcand_bestFromFirstHitSimTrkIdx", &tcand_bestFromFirstHitSimTrkIdx);
1677     } else {
1678       t->Branch("tcand_isTrue", &tcand_isTrue);
1679     }
1680     t->Branch("tcand_bestSimTrkShareFrac", &tcand_bestSimTrkShareFrac);
1681     t->Branch("tcand_bestSimTrkShareFracSimDenom", &tcand_bestSimTrkShareFracSimDenom);
1682     t->Branch("tcand_bestSimTrkShareFracSimClusterDenom", &tcand_bestSimTrkShareFracSimClusterDenom);
1683     t->Branch("tcand_bestSimTrkNChi2", &tcand_bestSimTrkNChi2);
1684     t->Branch("tcand_bestFromFirstHitSimTrkShareFrac", &tcand_bestFromFirstHitSimTrkShareFrac);
1685     t->Branch("tcand_bestFromFirstHitSimTrkShareFracSimDenom", &tcand_bestFromFirstHitSimTrkShareFracSimDenom);
1686     t->Branch("tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom",
1687               &tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom);
1688     t->Branch("tcand_bestFromFirstHitSimTrkNChi2", &tcand_bestFromFirstHitSimTrkNChi2);
1689     if (includeAllHits_) {
1690       t->Branch("tcand_hitIdx", &tcand_hitIdx);
1691       t->Branch("tcand_hitType", &tcand_hitType);
1692     }
1693   }
1694   if (includeTrackingParticles_) {
1695     //sim tracks
1696     t->Branch("sim_event", &sim_event);
1697     t->Branch("sim_bunchCrossing", &sim_bunchCrossing);
1698     t->Branch("sim_pdgId", &sim_pdgId);
1699     t->Branch("sim_genPdgIds", &sim_genPdgIds);
1700     t->Branch("sim_isFromBHadron", &sim_isFromBHadron);
1701     t->Branch("sim_px", &sim_px);
1702     t->Branch("sim_py", &sim_py);
1703     t->Branch("sim_pz", &sim_pz);
1704     t->Branch("sim_pt", &sim_pt);
1705     t->Branch("sim_eta", &sim_eta);
1706     t->Branch("sim_phi", &sim_phi);
1707     t->Branch("sim_pca_pt", &sim_pca_pt);
1708     t->Branch("sim_pca_eta", &sim_pca_eta);
1709     t->Branch("sim_pca_lambda", &sim_pca_lambda);
1710     t->Branch("sim_pca_cotTheta", &sim_pca_cotTheta);
1711     t->Branch("sim_pca_phi", &sim_pca_phi);
1712     t->Branch("sim_pca_dxy", &sim_pca_dxy);
1713     t->Branch("sim_pca_dz", &sim_pca_dz);
1714     t->Branch("sim_q", &sim_q);
1715     t->Branch("sim_nValid", &sim_nValid);
1716     t->Branch("sim_nPixel", &sim_nPixel);
1717     t->Branch("sim_nStrip", &sim_nStrip);
1718     t->Branch("sim_nLay", &sim_nLay);
1719     t->Branch("sim_nPixelLay", &sim_nPixelLay);
1720     t->Branch("sim_n3DLay", &sim_n3DLay);
1721     t->Branch("sim_nTrackerHits", &sim_nTrackerHits);
1722     t->Branch("sim_nRecoClusters", &sim_nRecoClusters);
1723     t->Branch("sim_trkIdx", &sim_trkIdx);
1724     t->Branch("sim_trkShareFrac", &sim_trkShareFrac);
1725     if (includeSeeds_) {
1726       t->Branch("sim_seedIdx", &sim_seedIdx);
1727     }
1728     t->Branch("sim_parentVtxIdx", &sim_parentVtxIdx);
1729     t->Branch("sim_decayVtxIdx", &sim_decayVtxIdx);
1730     if (includeAllHits_) {
1731       t->Branch("sim_simHitIdx", &sim_simHitIdx);
1732     }
1733   }
1734   if (includeAllHits_) {
1735     //pixels
1736     t->Branch("pix_isBarrel", &pix_isBarrel);
1737     pix_detId.book("pix", t);
1738     t->Branch("pix_trkIdx", &pix_trkIdx);
1739     if (includeOnTrackHitData_) {
1740       t->Branch("pix_onTrk_x", &pix_onTrk_x);
1741       t->Branch("pix_onTrk_y", &pix_onTrk_y);
1742       t->Branch("pix_onTrk_z", &pix_onTrk_z);
1743       t->Branch("pix_onTrk_xx", &pix_onTrk_xx);
1744       t->Branch("pix_onTrk_xy", &pix_onTrk_xy);
1745       t->Branch("pix_onTrk_yy", &pix_onTrk_yy);
1746       t->Branch("pix_onTrk_yz", &pix_onTrk_yz);
1747       t->Branch("pix_onTrk_zz", &pix_onTrk_zz);
1748       t->Branch("pix_onTrk_zx", &pix_onTrk_zx);
1749     }
1750     if (includeTrackCandidates_)
1751       t->Branch("pix_tcandIdx", &pix_tcandIdx);
1752     if (includeSeeds_) {
1753       t->Branch("pix_seeIdx", &pix_seeIdx);
1754     }
1755     if (includeTrackingParticles_) {
1756       t->Branch("pix_simHitIdx", &pix_simHitIdx);
1757       if (simHitBySignificance_) {
1758         t->Branch("pix_xySignificance", &pix_xySignificance);
1759       }
1760       t->Branch("pix_chargeFraction", &pix_chargeFraction);
1761       t->Branch("pix_simType", &pix_simType);
1762     }
1763     t->Branch("pix_x", &pix_x);
1764     t->Branch("pix_y", &pix_y);
1765     t->Branch("pix_z", &pix_z);
1766     t->Branch("pix_xx", &pix_xx);
1767     t->Branch("pix_xy", &pix_xy);
1768     t->Branch("pix_yy", &pix_yy);
1769     t->Branch("pix_yz", &pix_yz);
1770     t->Branch("pix_zz", &pix_zz);
1771     t->Branch("pix_zx", &pix_zx);
1772     t->Branch("pix_radL", &pix_radL);
1773     t->Branch("pix_bbxi", &pix_bbxi);
1774     t->Branch("pix_clustSizeCol", &pix_clustSizeCol);
1775     t->Branch("pix_clustSizeRow", &pix_clustSizeRow);
1776     t->Branch("pix_usedMask", &pix_usedMask);
1777     //strips
1778     if (includeStripHits_) {
1779       t->Branch("str_isBarrel", &str_isBarrel);
1780       str_detId.book("str", t);
1781       t->Branch("str_trkIdx", &str_trkIdx);
1782       if (includeOnTrackHitData_) {
1783         t->Branch("str_onTrk_x", &str_onTrk_x);
1784         t->Branch("str_onTrk_y", &str_onTrk_y);
1785         t->Branch("str_onTrk_z", &str_onTrk_z);
1786         t->Branch("str_onTrk_xx", &str_onTrk_xx);
1787         t->Branch("str_onTrk_xy", &str_onTrk_xy);
1788         t->Branch("str_onTrk_yy", &str_onTrk_yy);
1789         t->Branch("str_onTrk_yz", &str_onTrk_yz);
1790         t->Branch("str_onTrk_zz", &str_onTrk_zz);
1791         t->Branch("str_onTrk_zx", &str_onTrk_zx);
1792       }
1793       if (includeTrackCandidates_)
1794         t->Branch("str_tcandIdx", &str_tcandIdx);
1795       if (includeSeeds_) {
1796         t->Branch("str_seeIdx", &str_seeIdx);
1797       }
1798       if (includeTrackingParticles_) {
1799         t->Branch("str_simHitIdx", &str_simHitIdx);
1800         if (simHitBySignificance_) {
1801           t->Branch("str_xySignificance", &str_xySignificance);
1802         }
1803         t->Branch("str_chargeFraction", &str_chargeFraction);
1804         t->Branch("str_simType", &str_simType);
1805       }
1806       t->Branch("str_x", &str_x);
1807       t->Branch("str_y", &str_y);
1808       t->Branch("str_z", &str_z);
1809       t->Branch("str_xx", &str_xx);
1810       t->Branch("str_xy", &str_xy);
1811       t->Branch("str_yy", &str_yy);
1812       t->Branch("str_yz", &str_yz);
1813       t->Branch("str_zz", &str_zz);
1814       t->Branch("str_zx", &str_zx);
1815       t->Branch("str_radL", &str_radL);
1816       t->Branch("str_bbxi", &str_bbxi);
1817       t->Branch("str_chargePerCM", &str_chargePerCM);
1818       t->Branch("str_clustSize", &str_clustSize);
1819       t->Branch("str_usedMask", &str_usedMask);
1820       //matched hits
1821       t->Branch("glu_isBarrel", &glu_isBarrel);
1822       glu_detId.book("glu", t);
1823       t->Branch("glu_monoIdx", &glu_monoIdx);
1824       t->Branch("glu_stereoIdx", &glu_stereoIdx);
1825       if (includeSeeds_) {
1826         t->Branch("glu_seeIdx", &glu_seeIdx);
1827       }
1828       t->Branch("glu_x", &glu_x);
1829       t->Branch("glu_y", &glu_y);
1830       t->Branch("glu_z", &glu_z);
1831       t->Branch("glu_xx", &glu_xx);
1832       t->Branch("glu_xy", &glu_xy);
1833       t->Branch("glu_yy", &glu_yy);
1834       t->Branch("glu_yz", &glu_yz);
1835       t->Branch("glu_zz", &glu_zz);
1836       t->Branch("glu_zx", &glu_zx);
1837       t->Branch("glu_radL", &glu_radL);
1838       t->Branch("glu_bbxi", &glu_bbxi);
1839       t->Branch("glu_chargePerCM", &glu_chargePerCM);
1840       t->Branch("glu_clustSizeMono", &glu_clustSizeMono);
1841       t->Branch("glu_clustSizeStereo", &glu_clustSizeStereo);
1842       t->Branch("glu_usedMaskMono", &glu_usedMaskMono);
1843       t->Branch("glu_usedMaskStereo", &glu_usedMaskStereo);
1844     }
1845     //phase2 OT
1846     if (includePhase2OTHits_) {
1847       t->Branch("ph2_isBarrel", &ph2_isBarrel);
1848       ph2_detId.book("ph2", t);
1849       t->Branch("ph2_trkIdx", &ph2_trkIdx);
1850       if (includeOnTrackHitData_) {
1851         t->Branch("ph2_onTrk_x", &ph2_onTrk_x);
1852         t->Branch("ph2_onTrk_y", &ph2_onTrk_y);
1853         t->Branch("ph2_onTrk_z", &ph2_onTrk_z);
1854         t->Branch("ph2_onTrk_xx", &ph2_onTrk_xx);
1855         t->Branch("ph2_onTrk_xy", &ph2_onTrk_xy);
1856         t->Branch("ph2_onTrk_yy", &ph2_onTrk_yy);
1857         t->Branch("ph2_onTrk_yz", &ph2_onTrk_yz);
1858         t->Branch("ph2_onTrk_zz", &ph2_onTrk_zz);
1859         t->Branch("ph2_onTrk_zx", &ph2_onTrk_zx);
1860       }
1861       if (includeTrackCandidates_)
1862         t->Branch("ph2_tcandIdx", &ph2_tcandIdx);
1863       if (includeSeeds_) {
1864         t->Branch("ph2_seeIdx", &ph2_seeIdx);
1865       }
1866       if (includeTrackingParticles_) {
1867         t->Branch("ph2_simHitIdx", &ph2_simHitIdx);
1868         if (simHitBySignificance_) {
1869           t->Branch("ph2_xySignificance", &ph2_xySignificance);
1870         }
1871         t->Branch("ph2_simType", &ph2_simType);
1872       }
1873       t->Branch("ph2_x", &ph2_x);
1874       t->Branch("ph2_y", &ph2_y);
1875       t->Branch("ph2_z", &ph2_z);
1876       t->Branch("ph2_xx", &ph2_xx);
1877       t->Branch("ph2_xy", &ph2_xy);
1878       t->Branch("ph2_yy", &ph2_yy);
1879       t->Branch("ph2_yz", &ph2_yz);
1880       t->Branch("ph2_zz", &ph2_zz);
1881       t->Branch("ph2_zx", &ph2_zx);
1882       t->Branch("ph2_radL", &ph2_radL);
1883       t->Branch("ph2_bbxi", &ph2_bbxi);
1884       t->Branch("ph2_usedMask", &ph2_usedMask);
1885     }
1886     //invalid hits
1887     t->Branch("inv_isBarrel", &inv_isBarrel);
1888     if (includeStripHits_)
1889       inv_detId.book("inv", t);
1890     else
1891       inv_detId_phase2.book("inv", t);
1892     t->Branch("inv_type", &inv_type);
1893     //simhits
1894     if (includeTrackingParticles_) {
1895       if (includeStripHits_)
1896         simhit_detId.book("simhit", t);
1897       else
1898         simhit_detId_phase2.book("simhit", t);
1899       t->Branch("simhit_x", &simhit_x);
1900       t->Branch("simhit_y", &simhit_y);
1901       t->Branch("simhit_z", &simhit_z);
1902       if (saveSimHitsP3_) {
1903         t->Branch("simhit_px", &simhit_px);
1904         t->Branch("simhit_py", &simhit_py);
1905         t->Branch("simhit_pz", &simhit_pz);
1906       }
1907       t->Branch("simhit_particle", &simhit_particle);
1908       t->Branch("simhit_process", &simhit_process);
1909       t->Branch("simhit_eloss", &simhit_eloss);
1910       t->Branch("simhit_tof", &simhit_tof);
1911       t->Branch("simhit_simTrkIdx", &simhit_simTrkIdx);
1912       t->Branch("simhit_hitIdx", &simhit_hitIdx);
1913       t->Branch("simhit_hitType", &simhit_hitType);
1914     }
1915   }
1916   //beam spot
1917   t->Branch("bsp_x", &bsp_x, "bsp_x/F");
1918   t->Branch("bsp_y", &bsp_y, "bsp_y/F");
1919   t->Branch("bsp_z", &bsp_z, "bsp_z/F");
1920   t->Branch("bsp_sigmax", &bsp_sigmax, "bsp_sigmax/F");
1921   t->Branch("bsp_sigmay", &bsp_sigmay, "bsp_sigmay/F");
1922   t->Branch("bsp_sigmaz", &bsp_sigmaz, "bsp_sigmaz/F");
1923   if (includeSeeds_) {
1924     //seeds
1925     t->Branch("see_fitok", &see_fitok);
1926     t->Branch("see_px", &see_px);
1927     t->Branch("see_py", &see_py);
1928     t->Branch("see_pz", &see_pz);
1929     t->Branch("see_pt", &see_pt);
1930     t->Branch("see_eta", &see_eta);
1931     t->Branch("see_phi", &see_phi);
1932     t->Branch("see_dxy", &see_dxy);
1933     t->Branch("see_dz", &see_dz);
1934     t->Branch("see_ptErr", &see_ptErr);
1935     t->Branch("see_etaErr", &see_etaErr);
1936     t->Branch("see_phiErr", &see_phiErr);
1937     t->Branch("see_dxyErr", &see_dxyErr);
1938     t->Branch("see_dzErr", &see_dzErr);
1939     t->Branch("see_chi2", &see_chi2);
1940     t->Branch("see_statePt", &see_statePt);
1941     t->Branch("see_stateTrajX", &see_stateTrajX);
1942     t->Branch("see_stateTrajY", &see_stateTrajY);
1943     t->Branch("see_stateTrajPx", &see_stateTrajPx);
1944     t->Branch("see_stateTrajPy", &see_stateTrajPy);
1945     t->Branch("see_stateTrajPz", &see_stateTrajPz);
1946     t->Branch("see_stateTrajGlbX", &see_stateTrajGlbX);
1947     t->Branch("see_stateTrajGlbY", &see_stateTrajGlbY);
1948     t->Branch("see_stateTrajGlbZ", &see_stateTrajGlbZ);
1949     t->Branch("see_stateTrajGlbPx", &see_stateTrajGlbPx);
1950     t->Branch("see_stateTrajGlbPy", &see_stateTrajGlbPy);
1951     t->Branch("see_stateTrajGlbPz", &see_stateTrajGlbPz);
1952     if (addSeedCurvCov_) {
1953       t->Branch("see_stateCurvCov", &see_stateCurvCov);
1954     }
1955     t->Branch("see_q", &see_q);
1956     t->Branch("see_nValid", &see_nValid);
1957     t->Branch("see_nPixel", &see_nPixel);
1958     t->Branch("see_nGlued", &see_nGlued);
1959     t->Branch("see_nStrip", &see_nStrip);
1960     t->Branch("see_nPhase2OT", &see_nPhase2OT);
1961     t->Branch("see_nCluster", &see_nCluster);
1962     t->Branch("see_algo", &see_algo);
1963     t->Branch("see_stopReason", &see_stopReason);
1964     t->Branch("see_nCands", &see_nCands);
1965     t->Branch("see_trkIdx", &see_trkIdx);
1966     if (includeTrackCandidates_)
1967       t->Branch("see_tcandIdx", &see_tcandIdx);
1968     if (includeTrackingParticles_) {
1969       t->Branch("see_simTrkIdx", &see_simTrkIdx);
1970       t->Branch("see_simTrkShareFrac", &see_simTrkShareFrac);
1971       t->Branch("see_bestSimTrkIdx", &see_bestSimTrkIdx);
1972       t->Branch("see_bestFromFirstHitSimTrkIdx", &see_bestFromFirstHitSimTrkIdx);
1973     } else {
1974       t->Branch("see_isTrue", &see_isTrue);
1975     }
1976     t->Branch("see_bestSimTrkShareFrac", &see_bestSimTrkShareFrac);
1977     t->Branch("see_bestFromFirstHitSimTrkShareFrac", &see_bestFromFirstHitSimTrkShareFrac);
1978     if (includeAllHits_) {
1979       t->Branch("see_hitIdx", &see_hitIdx);
1980       t->Branch("see_hitType", &see_hitType);
1981     }
1982     //seed algo offset
1983     t->Branch("see_offset", &see_offset);
1984   }
1985 
1986   //vertices
1987   t->Branch("vtx_x", &vtx_x);
1988   t->Branch("vtx_y", &vtx_y);
1989   t->Branch("vtx_z", &vtx_z);
1990   t->Branch("vtx_xErr", &vtx_xErr);
1991   t->Branch("vtx_yErr", &vtx_yErr);
1992   t->Branch("vtx_zErr", &vtx_zErr);
1993   t->Branch("vtx_ndof", &vtx_ndof);
1994   t->Branch("vtx_chi2", &vtx_chi2);
1995   t->Branch("vtx_fake", &vtx_fake);
1996   t->Branch("vtx_valid", &vtx_valid);
1997   t->Branch("vtx_trkIdx", &vtx_trkIdx);
1998 
1999   // tracking vertices
2000   t->Branch("simvtx_event", &simvtx_event);
2001   t->Branch("simvtx_bunchCrossing", &simvtx_bunchCrossing);
2002   t->Branch("simvtx_processType", &simvtx_processType);
2003   t->Branch("simvtx_x", &simvtx_x);
2004   t->Branch("simvtx_y", &simvtx_y);
2005   t->Branch("simvtx_z", &simvtx_z);
2006   t->Branch("simvtx_sourceSimIdx", &simvtx_sourceSimIdx);
2007   t->Branch("simvtx_daughterSimIdx", &simvtx_daughterSimIdx);
2008 
2009   t->Branch("simpv_idx", &simpv_idx);
2010 
2011   //t->Branch("" , &);
2012 }
2013 
2014 TrackingNtuple::~TrackingNtuple() {
2015   // do anything here that needs to be done at desctruction time
2016   // (e.g. close files, deallocate resources etc.)
2017 }
2018 
2019 //
2020 // member functions
2021 //
2022 void TrackingNtuple::clearVariables() {
2023   ev_run = 0;
2024   ev_lumi = 0;
2025   ev_event = 0;
2026 
2027   //tracks
2028   trk_px.clear();
2029   trk_py.clear();
2030   trk_pz.clear();
2031   trk_pt.clear();
2032   trk_inner_px.clear();
2033   trk_inner_py.clear();
2034   trk_inner_pz.clear();
2035   trk_inner_pt.clear();
2036   trk_outer_px.clear();
2037   trk_outer_py.clear();
2038   trk_outer_pz.clear();
2039   trk_outer_pt.clear();
2040   trk_eta.clear();
2041   trk_lambda.clear();
2042   trk_cotTheta.clear();
2043   trk_phi.clear();
2044   trk_dxy.clear();
2045   trk_dz.clear();
2046   trk_dxyPV.clear();
2047   trk_dzPV.clear();
2048   trk_dxyClosestPV.clear();
2049   trk_dzClosestPV.clear();
2050   trk_ptErr.clear();
2051   trk_etaErr.clear();
2052   trk_lambdaErr.clear();
2053   trk_phiErr.clear();
2054   trk_dxyErr.clear();
2055   trk_dzErr.clear();
2056   trk_refpoint_x.clear();
2057   trk_refpoint_y.clear();
2058   trk_refpoint_z.clear();
2059   trk_nChi2.clear();
2060   trk_nChi2_1Dmod.clear();
2061   trk_ndof.clear();
2062   for (auto& mva : trk_mvas) {
2063     mva.clear();
2064   }
2065   for (auto& mask : trk_qualityMasks) {
2066     mask.clear();
2067   }
2068   trk_q.clear();
2069   trk_nValid.clear();
2070   trk_nLost.clear();
2071   trk_nInactive.clear();
2072   trk_nPixel.clear();
2073   trk_nStrip.clear();
2074   trk_nOuterLost.clear();
2075   trk_nInnerLost.clear();
2076   trk_nOuterInactive.clear();
2077   trk_nInnerInactive.clear();
2078   trk_nPixelLay.clear();
2079   trk_nStripLay.clear();
2080   trk_n3DLay.clear();
2081   trk_nLostLay.clear();
2082   trk_nCluster.clear();
2083   trk_algo.clear();
2084   trk_originalAlgo.clear();
2085   trk_algoMask.clear();
2086   trk_stopReason.clear();
2087   trk_isHP.clear();
2088   trk_seedIdx.clear();
2089   trk_vtxIdx.clear();
2090   trk_isTrue.clear();
2091   trk_bestSimTrkIdx.clear();
2092   trk_bestSimTrkShareFrac.clear();
2093   trk_bestSimTrkShareFracSimDenom.clear();
2094   trk_bestSimTrkShareFracSimClusterDenom.clear();
2095   trk_bestSimTrkNChi2.clear();
2096   trk_bestFromFirstHitSimTrkIdx.clear();
2097   trk_bestFromFirstHitSimTrkShareFrac.clear();
2098   trk_bestFromFirstHitSimTrkShareFracSimDenom.clear();
2099   trk_bestFromFirstHitSimTrkShareFracSimClusterDenom.clear();
2100   trk_bestFromFirstHitSimTrkNChi2.clear();
2101   trk_simTrkIdx.clear();
2102   trk_simTrkShareFrac.clear();
2103   trk_simTrkNChi2.clear();
2104   trk_hitIdx.clear();
2105   trk_hitType.clear();
2106   //track candidates
2107   tcand_pca_valid.clear();
2108   tcand_pca_px.clear();
2109   tcand_pca_py.clear();
2110   tcand_pca_pz.clear();
2111   tcand_pca_pt.clear();
2112   tcand_pca_eta.clear();
2113   tcand_pca_phi.clear();
2114   tcand_pca_dxy.clear();
2115   tcand_pca_dz.clear();
2116   tcand_pca_ptErr.clear();
2117   tcand_pca_etaErr.clear();
2118   tcand_pca_lambdaErr.clear();
2119   tcand_pca_phiErr.clear();
2120   tcand_pca_dxyErr.clear();
2121   tcand_pca_dzErr.clear();
2122   tcand_px.clear();
2123   tcand_py.clear();
2124   tcand_pz.clear();
2125   tcand_pt.clear();
2126   tcand_x.clear();
2127   tcand_y.clear();
2128   tcand_z.clear();
2129   tcand_qbpErr.clear();
2130   tcand_lambdaErr.clear();
2131   tcand_phiErr.clear();
2132   tcand_xtErr.clear();
2133   tcand_ytErr.clear();
2134   tcand_ndof.clear();
2135   tcand_q.clear();
2136   tcand_nValid.clear();
2137   tcand_nPixel.clear();
2138   tcand_nStrip.clear();
2139   tcand_nCluster.clear();
2140   tcand_algo.clear();
2141   tcand_stopReason.clear();
2142   tcand_seedIdx.clear();
2143   tcand_vtxIdx.clear();
2144   tcand_isTrue.clear();
2145   tcand_bestSimTrkIdx.clear();
2146   tcand_bestSimTrkShareFrac.clear();
2147   tcand_bestSimTrkShareFracSimDenom.clear();
2148   tcand_bestSimTrkShareFracSimClusterDenom.clear();
2149   tcand_bestSimTrkNChi2.clear();
2150   tcand_bestFromFirstHitSimTrkIdx.clear();
2151   tcand_bestFromFirstHitSimTrkShareFrac.clear();
2152   tcand_bestFromFirstHitSimTrkShareFracSimDenom.clear();
2153   tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom.clear();
2154   tcand_bestFromFirstHitSimTrkNChi2.clear();
2155   tcand_simTrkShareFrac.clear();
2156   tcand_simTrkNChi2.clear();
2157   tcand_simTrkIdx.clear();
2158   tcand_hitIdx.clear();
2159   tcand_hitType.clear();
2160   //sim tracks
2161   sim_event.clear();
2162   sim_bunchCrossing.clear();
2163   sim_pdgId.clear();
2164   sim_genPdgIds.clear();
2165   sim_isFromBHadron.clear();
2166   sim_px.clear();
2167   sim_py.clear();
2168   sim_pz.clear();
2169   sim_pt.clear();
2170   sim_eta.clear();
2171   sim_phi.clear();
2172   sim_pca_pt.clear();
2173   sim_pca_eta.clear();
2174   sim_pca_lambda.clear();
2175   sim_pca_cotTheta.clear();
2176   sim_pca_phi.clear();
2177   sim_pca_dxy.clear();
2178   sim_pca_dz.clear();
2179   sim_q.clear();
2180   sim_nValid.clear();
2181   sim_nPixel.clear();
2182   sim_nStrip.clear();
2183   sim_nLay.clear();
2184   sim_nPixelLay.clear();
2185   sim_n3DLay.clear();
2186   sim_nTrackerHits.clear();
2187   sim_nRecoClusters.clear();
2188   sim_trkIdx.clear();
2189   sim_seedIdx.clear();
2190   sim_trkShareFrac.clear();
2191   sim_parentVtxIdx.clear();
2192   sim_decayVtxIdx.clear();
2193   sim_simHitIdx.clear();
2194   //pixels
2195   pix_isBarrel.clear();
2196   pix_detId.clear();
2197   pix_trkIdx.clear();
2198   pix_onTrk_x.clear();
2199   pix_onTrk_y.clear();
2200   pix_onTrk_z.clear();
2201   pix_onTrk_xx.clear();
2202   pix_onTrk_xy.clear();
2203   pix_onTrk_yy.clear();
2204   pix_onTrk_yz.clear();
2205   pix_onTrk_zz.clear();
2206   pix_onTrk_zx.clear();
2207   pix_tcandIdx.clear();
2208   pix_seeIdx.clear();
2209   pix_simHitIdx.clear();
2210   pix_xySignificance.clear();
2211   pix_chargeFraction.clear();
2212   pix_simType.clear();
2213   pix_x.clear();
2214   pix_y.clear();
2215   pix_z.clear();
2216   pix_xx.clear();
2217   pix_xy.clear();
2218   pix_yy.clear();
2219   pix_yz.clear();
2220   pix_zz.clear();
2221   pix_zx.clear();
2222   pix_radL.clear();
2223   pix_bbxi.clear();
2224   pix_clustSizeCol.clear();
2225   pix_clustSizeRow.clear();
2226   pix_usedMask.clear();
2227   //strips
2228   str_isBarrel.clear();
2229   str_detId.clear();
2230   str_trkIdx.clear();
2231   str_onTrk_x.clear();
2232   str_onTrk_y.clear();
2233   str_onTrk_z.clear();
2234   str_onTrk_xx.clear();
2235   str_onTrk_xy.clear();
2236   str_onTrk_yy.clear();
2237   str_onTrk_yz.clear();
2238   str_onTrk_zz.clear();
2239   str_onTrk_zx.clear();
2240   str_tcandIdx.clear();
2241   str_seeIdx.clear();
2242   str_simHitIdx.clear();
2243   str_xySignificance.clear();
2244   str_chargeFraction.clear();
2245   str_simType.clear();
2246   str_x.clear();
2247   str_y.clear();
2248   str_z.clear();
2249   str_xx.clear();
2250   str_xy.clear();
2251   str_yy.clear();
2252   str_yz.clear();
2253   str_zz.clear();
2254   str_zx.clear();
2255   str_radL.clear();
2256   str_bbxi.clear();
2257   str_chargePerCM.clear();
2258   str_clustSize.clear();
2259   str_usedMask.clear();
2260   //matched hits
2261   glu_isBarrel.clear();
2262   glu_detId.clear();
2263   glu_monoIdx.clear();
2264   glu_stereoIdx.clear();
2265   glu_seeIdx.clear();
2266   glu_x.clear();
2267   glu_y.clear();
2268   glu_z.clear();
2269   glu_xx.clear();
2270   glu_xy.clear();
2271   glu_yy.clear();
2272   glu_yz.clear();
2273   glu_zz.clear();
2274   glu_zx.clear();
2275   glu_radL.clear();
2276   glu_bbxi.clear();
2277   glu_chargePerCM.clear();
2278   glu_clustSizeMono.clear();
2279   glu_clustSizeStereo.clear();
2280   glu_usedMaskMono.clear();
2281   glu_usedMaskStereo.clear();
2282   //phase2 OT
2283   ph2_isBarrel.clear();
2284   ph2_detId.clear();
2285   ph2_trkIdx.clear();
2286   ph2_onTrk_x.clear();
2287   ph2_onTrk_y.clear();
2288   ph2_onTrk_z.clear();
2289   ph2_onTrk_xx.clear();
2290   ph2_onTrk_xy.clear();
2291   ph2_onTrk_yy.clear();
2292   ph2_onTrk_yz.clear();
2293   ph2_onTrk_zz.clear();
2294   ph2_onTrk_zx.clear();
2295   ph2_tcandIdx.clear();
2296   ph2_seeIdx.clear();
2297   ph2_xySignificance.clear();
2298   ph2_simHitIdx.clear();
2299   ph2_simType.clear();
2300   ph2_x.clear();
2301   ph2_y.clear();
2302   ph2_z.clear();
2303   ph2_xx.clear();
2304   ph2_xy.clear();
2305   ph2_yy.clear();
2306   ph2_yz.clear();
2307   ph2_zz.clear();
2308   ph2_zx.clear();
2309   ph2_radL.clear();
2310   ph2_bbxi.clear();
2311   ph2_usedMask.clear();
2312   //invalid hits
2313   inv_isBarrel.clear();
2314   inv_detId.clear();
2315   inv_detId_phase2.clear();
2316   inv_type.clear();
2317   // simhits
2318   simhit_detId.clear();
2319   simhit_detId_phase2.clear();
2320   simhit_x.clear();
2321   simhit_y.clear();
2322   simhit_z.clear();
2323   simhit_px.clear();
2324   simhit_py.clear();
2325   simhit_pz.clear();
2326   simhit_particle.clear();
2327   simhit_process.clear();
2328   simhit_eloss.clear();
2329   simhit_tof.clear();
2330   //simhit_simTrackId.clear();
2331   simhit_simTrkIdx.clear();
2332   simhit_hitIdx.clear();
2333   simhit_hitType.clear();
2334   //beamspot
2335   bsp_x = -9999.;
2336   bsp_y = -9999.;
2337   bsp_z = -9999.;
2338   bsp_sigmax = -9999.;
2339   bsp_sigmay = -9999.;
2340   bsp_sigmaz = -9999.;
2341   //seeds
2342   see_fitok.clear();
2343   see_px.clear();
2344   see_py.clear();
2345   see_pz.clear();
2346   see_pt.clear();
2347   see_eta.clear();
2348   see_phi.clear();
2349   see_dxy.clear();
2350   see_dz.clear();
2351   see_ptErr.clear();
2352   see_etaErr.clear();
2353   see_phiErr.clear();
2354   see_dxyErr.clear();
2355   see_dzErr.clear();
2356   see_chi2.clear();
2357   see_statePt.clear();
2358   see_stateTrajX.clear();
2359   see_stateTrajY.clear();
2360   see_stateTrajPx.clear();
2361   see_stateTrajPy.clear();
2362   see_stateTrajPz.clear();
2363   see_stateTrajGlbX.clear();
2364   see_stateTrajGlbY.clear();
2365   see_stateTrajGlbZ.clear();
2366   see_stateTrajGlbPx.clear();
2367   see_stateTrajGlbPy.clear();
2368   see_stateTrajGlbPz.clear();
2369   see_stateCurvCov.clear();
2370   see_q.clear();
2371   see_nValid.clear();
2372   see_nPixel.clear();
2373   see_nGlued.clear();
2374   see_nStrip.clear();
2375   see_nPhase2OT.clear();
2376   see_nCluster.clear();
2377   see_algo.clear();
2378   see_stopReason.clear();
2379   see_nCands.clear();
2380   see_trkIdx.clear();
2381   see_tcandIdx.clear();
2382   see_isTrue.clear();
2383   see_bestSimTrkIdx.clear();
2384   see_bestSimTrkShareFrac.clear();
2385   see_bestFromFirstHitSimTrkIdx.clear();
2386   see_bestFromFirstHitSimTrkShareFrac.clear();
2387   see_simTrkIdx.clear();
2388   see_simTrkShareFrac.clear();
2389   see_hitIdx.clear();
2390   see_hitType.clear();
2391   //seed algo offset
2392   see_offset.clear();
2393 
2394   // vertices
2395   vtx_x.clear();
2396   vtx_y.clear();
2397   vtx_z.clear();
2398   vtx_xErr.clear();
2399   vtx_yErr.clear();
2400   vtx_zErr.clear();
2401   vtx_ndof.clear();
2402   vtx_chi2.clear();
2403   vtx_fake.clear();
2404   vtx_valid.clear();
2405   vtx_trkIdx.clear();
2406 
2407   // Tracking vertices
2408   simvtx_event.clear();
2409   simvtx_bunchCrossing.clear();
2410   simvtx_processType.clear();
2411   simvtx_x.clear();
2412   simvtx_y.clear();
2413   simvtx_z.clear();
2414   simvtx_sourceSimIdx.clear();
2415   simvtx_daughterSimIdx.clear();
2416   simpv_idx.clear();
2417 }
2418 
2419 // ------------ method called for each event  ------------
2420 void TrackingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
2421   using namespace edm;
2422   using namespace reco;
2423   using namespace std;
2424 
2425   const auto& mf = iSetup.getData(mfToken_);
2426   const TrackerTopology& tTopo = iSetup.getData(tTopoToken_);
2427   const TrackerGeometry& tracker = iSetup.getData(tGeomToken_);
2428 
2429   edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator;
2430   iEvent.getByToken(trackAssociatorToken_, theAssociator);
2431   const reco::TrackToTrackingParticleAssociator& associatorByHits = *theAssociator;
2432 
2433   LogDebug("TrackingNtuple") << "Analyzing new event";
2434 
2435   //initialize tree variables
2436   clearVariables();
2437 
2438   // FIXME: we really need to move to edm::View for reading the
2439   // TrackingParticles... Unfortunately it has non-trivial
2440   // consequences on the associator/association interfaces etc.
2441   TrackingParticleRefVector tmpTP;
2442   const TrackingParticleRefVector* tmpTPptr = nullptr;
2443   edm::Handle<TrackingParticleCollection> TPCollectionH;
2444   edm::Handle<TrackingParticleRefVector> TPCollectionHRefVector;
2445 
2446   if (!trackingParticleToken_.isUninitialized()) {
2447     iEvent.getByToken(trackingParticleToken_, TPCollectionH);
2448     for (size_t i = 0, size = TPCollectionH->size(); i < size; ++i) {
2449       tmpTP.push_back(TrackingParticleRef(TPCollectionH, i));
2450     }
2451     tmpTPptr = &tmpTP;
2452   } else {
2453     iEvent.getByToken(trackingParticleRefToken_, TPCollectionHRefVector);
2454     tmpTPptr = TPCollectionHRefVector.product();
2455   }
2456   const TrackingParticleRefVector& tpCollection = *tmpTPptr;
2457 
2458   // Fill mapping from Ref::key() to index
2459   TrackingParticleRefKeyToIndex tpKeyToIndex;
2460   for (size_t i = 0; i < tpCollection.size(); ++i) {
2461     tpKeyToIndex[tpCollection[i].key()] = i;
2462   }
2463 
2464   // tracking vertices
2465   edm::Handle<TrackingVertexCollection> htv;
2466   iEvent.getByToken(trackingVertexToken_, htv);
2467   const TrackingVertexCollection& tvs = *htv;
2468 
2469   // Fill mapping from Ref::key() to index
2470   TrackingVertexRefVector tvRefs;
2471   TrackingVertexRefKeyToIndex tvKeyToIndex;
2472   for (size_t i = 0; i < tvs.size(); ++i) {
2473     const TrackingVertex& v = tvs[i];
2474     if (!includeOOT_ && v.eventId().bunchCrossing() != 0)  // Ignore OOTPU
2475       continue;
2476     tvKeyToIndex[i] = tvRefs.size();
2477     tvRefs.push_back(TrackingVertexRef(htv, i));
2478   }
2479 
2480   //get association maps, etc.
2481   Handle<ClusterTPAssociation> pCluster2TPListH;
2482   iEvent.getByToken(clusterTPMapToken_, pCluster2TPListH);
2483   const ClusterTPAssociation& clusterToTPMap = *pCluster2TPListH;
2484   edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
2485   iEvent.getByToken(simHitTPMapToken_, simHitsTPAssoc);
2486 
2487   // TP -> cluster count
2488   TrackingParticleRefKeyToCount tpKeyToClusterCount;
2489   for (const auto& clusterTP : clusterToTPMap) {
2490     tpKeyToClusterCount[clusterTP.second.key()] += 1;
2491   }
2492 
2493   // SimHit key -> index mapping
2494   SimHitRefKeyToIndex simHitRefKeyToIndex;
2495 
2496   //make a list to link TrackingParticles to its simhits
2497   std::vector<TPHitIndex> tpHitList;
2498 
2499   // Count the number of reco cluster per TP
2500 
2501   std::set<edm::ProductID> hitProductIds;
2502   std::map<edm::ProductID, size_t> seedCollToOffset;
2503 
2504   ev_run = iEvent.id().run();
2505   ev_lumi = iEvent.id().luminosityBlock();
2506   ev_event = iEvent.id().event();
2507 
2508   // Digi->Sim links for pixels and strips
2509   edm::Handle<edm::DetSetVector<PixelDigiSimLink>> pixelDigiSimLinksHandle;
2510   iEvent.getByToken(pixelSimLinkToken_, pixelDigiSimLinksHandle);
2511   const auto& pixelDigiSimLinks = *pixelDigiSimLinksHandle;
2512 
2513   edm::Handle<edm::DetSetVector<StripDigiSimLink>> stripDigiSimLinksHandle;
2514   iEvent.getByToken(stripSimLinkToken_, stripDigiSimLinksHandle);
2515 
2516   // Phase2 OT DigiSimLink
2517   edm::Handle<edm::DetSetVector<PixelDigiSimLink>> siphase2OTSimLinksHandle;
2518   iEvent.getByToken(siphase2OTSimLinksToken_, siphase2OTSimLinksHandle);
2519 
2520   //beamspot
2521   Handle<reco::BeamSpot> recoBeamSpotHandle;
2522   iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
2523   BeamSpot const& bs = *recoBeamSpotHandle;
2524   fillBeamSpot(bs);
2525 
2526   //prapare list to link matched hits to collection
2527   vector<pair<int, int>> monoStereoClusterList;
2528   if (includeAllHits_) {
2529     // simhits, only if TPs are saved as well
2530     if (includeTrackingParticles_) {
2531       fillSimHits(tracker, tpKeyToIndex, *simHitsTPAssoc, tTopo, simHitRefKeyToIndex, tpHitList);
2532     }
2533 
2534     //pixel hits
2535     fillPixelHits(iEvent,
2536                   tracker,
2537                   clusterToTPMap,
2538                   tpKeyToIndex,
2539                   *simHitsTPAssoc,
2540                   pixelDigiSimLinks,
2541                   tTopo,
2542                   simHitRefKeyToIndex,
2543                   hitProductIds);
2544 
2545     //strip hits
2546     if (includeStripHits_) {
2547       LogDebug("TrackingNtuple") << "foundStripSimLink";
2548       const auto& stripDigiSimLinks = *stripDigiSimLinksHandle;
2549       fillStripRphiStereoHits(iEvent,
2550                               tracker,
2551                               clusterToTPMap,
2552                               tpKeyToIndex,
2553                               *simHitsTPAssoc,
2554                               stripDigiSimLinks,
2555                               tTopo,
2556                               simHitRefKeyToIndex,
2557                               hitProductIds);
2558 
2559       //matched hits
2560       fillStripMatchedHits(iEvent, tracker, tTopo, monoStereoClusterList);
2561     }
2562 
2563     if (includePhase2OTHits_) {
2564       LogDebug("TrackingNtuple") << "foundPhase2OTSimLinks";
2565       const auto& phase2OTSimLinks = *siphase2OTSimLinksHandle;
2566       fillPhase2OTHits(iEvent,
2567                        clusterToTPMap,
2568                        tracker,
2569                        tpKeyToIndex,
2570                        *simHitsTPAssoc,
2571                        phase2OTSimLinks,
2572                        tTopo,
2573                        simHitRefKeyToIndex,
2574                        hitProductIds);
2575     }
2576   }
2577 
2578   //seeds
2579   if (includeSeeds_) {
2580     fillSeeds(iEvent,
2581               tpCollection,
2582               tpKeyToIndex,
2583               bs,
2584               tracker,
2585               associatorByHits,
2586               clusterToTPMap,
2587               mf,
2588               tTopo,
2589               monoStereoClusterList,
2590               hitProductIds,
2591               seedCollToOffset);
2592   }
2593 
2594   //tracks
2595   edm::Handle<edm::View<reco::Track>> tracksHandle;
2596   iEvent.getByToken(trackToken_, tracksHandle);
2597   const edm::View<reco::Track>& tracks = *tracksHandle;
2598   // The associator interfaces really need to be fixed...
2599   edm::RefToBaseVector<reco::Track> trackRefs;
2600   for (edm::View<Track>::size_type i = 0; i < tracks.size(); ++i) {
2601     trackRefs.push_back(tracks.refAt(i));
2602   }
2603   std::vector<const MVACollection*> mvaColls;
2604   std::vector<const QualityMaskCollection*> qualColls;
2605   if (includeMVA_) {
2606     edm::Handle<MVACollection> hmva;
2607     edm::Handle<QualityMaskCollection> hqual;
2608 
2609     for (const auto& tokenTpl : mvaQualityCollectionTokens_) {
2610       iEvent.getByToken(std::get<0>(tokenTpl), hmva);
2611       iEvent.getByToken(std::get<1>(tokenTpl), hqual);
2612 
2613       mvaColls.push_back(hmva.product());
2614       qualColls.push_back(hqual.product());
2615       if (mvaColls.back()->size() != tracks.size()) {
2616         throw cms::Exception("Configuration")
2617             << "Inconsistency in track collection and MVA sizes. Track collection has " << tracks.size()
2618             << " tracks, whereas the MVA " << (mvaColls.size() - 1) << " has " << mvaColls.back()->size()
2619             << " entries. Double-check your configuration.";
2620       }
2621       if (qualColls.back()->size() != tracks.size()) {
2622         throw cms::Exception("Configuration")
2623             << "Inconsistency in track collection and quality mask sizes. Track collection has " << tracks.size()
2624             << " tracks, whereas the quality mask " << (qualColls.size() - 1) << " has " << qualColls.back()->size()
2625             << " entries. Double-check your configuration.";
2626       }
2627     }
2628   }
2629 
2630   edm::Handle<reco::VertexCollection> vertices;
2631   iEvent.getByToken(vertexToken_, vertices);
2632 
2633   fillTracks(trackRefs,
2634              tracker,
2635              tpCollection,
2636              tpKeyToIndex,
2637              tpKeyToClusterCount,
2638              mf,
2639              bs,
2640              *vertices,
2641              associatorByHits,
2642              clusterToTPMap,
2643              tTopo,
2644              hitProductIds,
2645              seedCollToOffset,
2646              mvaColls,
2647              qualColls);
2648 
2649   if (includeTrackCandidates_) {
2650     //candidates
2651     for (auto const& token : candidateTokens_) {
2652       auto const tCandsHandle = iEvent.getHandle(token);
2653 
2654       edm::EDConsumerBase::Labels labels;
2655       labelsForToken(token, labels);
2656       TString label = labels.module;
2657       label.ReplaceAll("TrackCandidates", "");
2658       label.ReplaceAll("muonSeeded", "muonSeededStep");
2659       int algo = reco::TrackBase::algoByName(label.Data());
2660 
2661       fillCandidates(tCandsHandle,
2662                      algo,
2663                      tpCollection,
2664                      tpKeyToIndex,
2665                      tpKeyToClusterCount,
2666                      mf,
2667                      bs,
2668                      *vertices,
2669                      associatorByHits,
2670                      clusterToTPMap,
2671                      tracker,
2672                      tTopo,
2673                      hitProductIds,
2674                      seedCollToOffset);
2675     }
2676   }
2677 
2678   //tracking particles
2679   //sort association maps with simHits
2680   std::sort(tpHitList.begin(), tpHitList.end(), tpHitIndexListLessSort);
2681   fillTrackingParticles(
2682       iEvent, iSetup, trackRefs, bs, tpCollection, tvKeyToIndex, associatorByHits, tpHitList, tpKeyToClusterCount);
2683 
2684   // vertices
2685   fillVertices(*vertices, trackRefs);
2686 
2687   // tracking vertices
2688   fillTrackingVertices(tvRefs, tpKeyToIndex);
2689 
2690   t->Fill();
2691 }
2692 
2693 void TrackingNtuple::fillBeamSpot(const reco::BeamSpot& bs) {
2694   bsp_x = bs.x0();
2695   bsp_y = bs.y0();
2696   bsp_z = bs.x0();
2697   bsp_sigmax = bs.BeamWidthX();
2698   bsp_sigmay = bs.BeamWidthY();
2699   bsp_sigmaz = bs.sigmaZ();
2700 }
2701 
2702 namespace {
2703   template <typename SimLink>
2704   struct GetCluster;
2705   template <>
2706   struct GetCluster<PixelDigiSimLink> {
2707     static const SiPixelCluster& call(const OmniClusterRef& cluster) { return cluster.pixelCluster(); }
2708   };
2709   template <>
2710   struct GetCluster<StripDigiSimLink> {
2711     static const SiStripCluster& call(const OmniClusterRef& cluster) { return cluster.stripCluster(); }
2712   };
2713 }  // namespace
2714 
2715 template <typename SimLink>
2716 TrackingNtuple::SimHitData TrackingNtuple::matchCluster(
2717     const OmniClusterRef& cluster,
2718     DetId hitId,
2719     int clusterKey,
2720     const TrackingRecHit& hit,
2721     const ClusterTPAssociation& clusterToTPMap,
2722     const TrackingParticleRefKeyToIndex& tpKeyToIndex,
2723     const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
2724     const edm::DetSetVector<SimLink>& digiSimLinks,
2725     const SimHitRefKeyToIndex& simHitRefKeyToIndex,
2726     HitType hitType) {
2727   SimHitData ret;
2728 
2729   std::map<unsigned int, double> simTrackIdToChargeFraction;
2730   if (hitType == HitType::Phase2OT)
2731     simTrackIdToChargeFraction = chargeFraction(cluster.phase2OTCluster(), hitId, digiSimLinks);
2732   else
2733     simTrackIdToChargeFraction = chargeFraction(GetCluster<SimLink>::call(cluster), hitId, digiSimLinks);
2734 
2735   float h_x = 0, h_y = 0;
2736   float h_xx = 0, h_xy = 0, h_yy = 0;
2737   if (simHitBySignificance_) {
2738     h_x = hit.localPosition().x();
2739     h_y = hit.localPosition().y();
2740     h_xx = hit.localPositionError().xx();
2741     h_xy = hit.localPositionError().xy();
2742     h_yy = hit.localPositionError().yy();
2743   }
2744 
2745   ret.type = HitSimType::Noise;
2746   auto range = clusterToTPMap.equal_range(cluster);
2747   if (range.first != range.second) {
2748     for (auto ip = range.first; ip != range.second; ++ip) {
2749       const TrackingParticleRef& trackingParticle = ip->second;
2750 
2751       // Find out if the cluster is from signal/ITPU/OOTPU
2752       const auto event = trackingParticle->eventId().event();
2753       const auto bx = trackingParticle->eventId().bunchCrossing();
2754       HitSimType type = HitSimType::OOTPileup;
2755       if (bx == 0) {
2756         type = (event == 0 ? HitSimType::Signal : HitSimType::ITPileup);
2757       } else
2758         type = HitSimType::OOTPileup;
2759       ret.type = static_cast<HitSimType>(std::min(static_cast<int>(ret.type), static_cast<int>(type)));
2760 
2761       // Limit to only input TrackingParticles (usually signal+ITPU)
2762       auto tpIndex = tpKeyToIndex.find(trackingParticle.key());
2763       if (tpIndex == tpKeyToIndex.end())
2764         continue;
2765 
2766       //now get the corresponding sim hit
2767       std::pair<TrackingParticleRef, TrackPSimHitRef> simHitTPpairWithDummyTP(trackingParticle, TrackPSimHitRef());
2768       //SimHit is dummy: for simHitTPAssociationListGreater sorting only the TP is needed
2769       auto range = std::equal_range(simHitsTPAssoc.begin(),
2770                                     simHitsTPAssoc.end(),
2771                                     simHitTPpairWithDummyTP,
2772                                     SimHitTPAssociationProducer::simHitTPAssociationListGreater);
2773       bool foundSimHit = false;
2774       bool foundElectron = false;
2775       int foundNonElectrons = 0;
2776       for (auto ip = range.first; ip != range.second; ++ip) {
2777         TrackPSimHitRef TPhit = ip->second;
2778         DetId dId = DetId(TPhit->detUnitId());
2779         if (dId.rawId() == hitId.rawId()) {
2780           // skip electron SimHits for non-electron TPs also here
2781           if (std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
2782             continue;  //electrons found
2783           } else {
2784             foundNonElectrons++;
2785           }
2786         }
2787       }
2788 
2789       float minSignificance = 1e12;
2790       if (simHitBySignificance_) {  //save the best matching hit
2791 
2792         int simHitKey = -1;
2793         edm::ProductID simHitID;
2794         for (auto ip = range.first; ip != range.second; ++ip) {
2795           TrackPSimHitRef TPhit = ip->second;
2796           DetId dId = DetId(TPhit->detUnitId());
2797           if (dId.rawId() == hitId.rawId()) {
2798             // skip electron SimHits for non-electron TPs also here
2799             if (std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
2800               foundElectron = true;
2801               if (!keepEleSimHits_)
2802                 continue;
2803             }
2804 
2805             float sx = TPhit->localPosition().x();
2806             float sy = TPhit->localPosition().y();
2807             float dx = sx - h_x;
2808             float dy = sy - h_y;
2809             float sig = (dx * dx * h_yy - 2 * dx * dy * h_xy + dy * dy * h_xx) / (h_xx * h_yy - h_xy * h_xy);
2810 
2811             if (sig < minSignificance) {
2812               minSignificance = sig;
2813               foundSimHit = true;
2814               simHitKey = TPhit.key();
2815               simHitID = TPhit.id();
2816             }
2817           }
2818         }  //loop over matching hits
2819 
2820         auto simHitIndex = simHitRefKeyToIndex.at(std::make_pair(simHitKey, simHitID));
2821         ret.matchingSimHit.push_back(simHitIndex);
2822 
2823         double chargeFraction = 0.;
2824         for (const SimTrack& simtrk : trackingParticle->g4Tracks()) {
2825           auto found = simTrackIdToChargeFraction.find(simtrk.trackId());
2826           if (found != simTrackIdToChargeFraction.end()) {
2827             chargeFraction += found->second;
2828           }
2829         }
2830         ret.xySignificance.push_back(minSignificance);
2831         ret.chargeFraction.push_back(chargeFraction);
2832 
2833         // only for debug prints
2834         ret.bunchCrossing.push_back(bx);
2835         ret.event.push_back(event);
2836 
2837         simhit_hitIdx[simHitIndex].push_back(clusterKey);
2838         simhit_hitType[simHitIndex].push_back(static_cast<int>(hitType));
2839 
2840       } else {  //save all matching hits
2841         for (auto ip = range.first; ip != range.second; ++ip) {
2842           TrackPSimHitRef TPhit = ip->second;
2843           DetId dId = DetId(TPhit->detUnitId());
2844           if (dId.rawId() == hitId.rawId()) {
2845             // skip electron SimHits for non-electron TPs also here
2846             if (std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
2847               foundElectron = true;
2848               if (!keepEleSimHits_)
2849                 continue;
2850               if (foundNonElectrons > 0)
2851                 continue;  //prioritize: skip electrons if non-electrons are present
2852             }
2853 
2854             foundSimHit = true;
2855             auto simHitKey = TPhit.key();
2856             auto simHitID = TPhit.id();
2857 
2858             auto simHitIndex = simHitRefKeyToIndex.at(std::make_pair(simHitKey, simHitID));
2859             ret.matchingSimHit.push_back(simHitIndex);
2860 
2861             double chargeFraction = 0.;
2862             for (const SimTrack& simtrk : trackingParticle->g4Tracks()) {
2863               auto found = simTrackIdToChargeFraction.find(simtrk.trackId());
2864               if (found != simTrackIdToChargeFraction.end()) {
2865                 chargeFraction += found->second;
2866               }
2867             }
2868             ret.xySignificance.push_back(minSignificance);
2869             ret.chargeFraction.push_back(chargeFraction);
2870 
2871             // only for debug prints
2872             ret.bunchCrossing.push_back(bx);
2873             ret.event.push_back(event);
2874 
2875             simhit_hitIdx[simHitIndex].push_back(clusterKey);
2876             simhit_hitType[simHitIndex].push_back(static_cast<int>(hitType));
2877           }
2878         }
2879       }  //if/else simHitBySignificance_
2880       if (!foundSimHit) {
2881         // In case we didn't find a simhit because of filtered-out
2882         // electron SimHit, just ignore the missing SimHit.
2883         if (foundElectron && !keepEleSimHits_)
2884           continue;
2885 
2886         auto ex = cms::Exception("LogicError")
2887                   << "Did not find SimHit for reco hit DetId " << hitId.rawId() << " for TP " << trackingParticle.key()
2888                   << " bx:event " << bx << ":" << event << " PDGid " << trackingParticle->pdgId() << " q "
2889                   << trackingParticle->charge() << " p4 " << trackingParticle->p4() << " nG4 "
2890                   << trackingParticle->g4Tracks().size() << ".\nFound SimHits from detectors ";
2891         for (auto ip = range.first; ip != range.second; ++ip) {
2892           TrackPSimHitRef TPhit = ip->second;
2893           DetId dId = DetId(TPhit->detUnitId());
2894           ex << dId.rawId() << " ";
2895         }
2896         if (trackingParticle->eventId().event() != 0) {
2897           ex << "\nSince this is a TrackingParticle from pileup, check that you're running the pileup mixing in "
2898                 "playback mode.";
2899         }
2900         throw ex;
2901       }
2902     }
2903   }
2904 
2905   return ret;
2906 }
2907 
2908 void TrackingNtuple::fillSimHits(const TrackerGeometry& tracker,
2909                                  const TrackingParticleRefKeyToIndex& tpKeyToIndex,
2910                                  const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
2911                                  const TrackerTopology& tTopo,
2912                                  SimHitRefKeyToIndex& simHitRefKeyToIndex,
2913                                  std::vector<TPHitIndex>& tpHitList) {
2914   for (const auto& assoc : simHitsTPAssoc) {
2915     auto tpKey = assoc.first.key();
2916 
2917     // SimHitTPAssociationList can contain more TrackingParticles than
2918     // what are given to this EDAnalyzer, so we can filter those out here.
2919     auto found = tpKeyToIndex.find(tpKey);
2920     if (found == tpKeyToIndex.end())
2921       continue;
2922     const auto tpIndex = found->second;
2923 
2924     // skip non-tracker simhits (mostly muons)
2925     const auto& simhit = *(assoc.second);
2926     auto detId = DetId(simhit.detUnitId());
2927     if (detId.det() != DetId::Tracker)
2928       continue;
2929 
2930     // Skip electron SimHits for non-electron TrackingParticles to
2931     // filter out delta rays. The delta ray hits just confuse. If we
2932     // need them later, let's add them as a separate "collection" of
2933     // hits of a TP
2934     const TrackingParticle& tp = *(assoc.first);
2935     if (!keepEleSimHits_ && std::abs(simhit.particleType()) == 11 && std::abs(tp.pdgId()) != 11)
2936       continue;
2937 
2938     auto simHitKey = std::make_pair(assoc.second.key(), assoc.second.id());
2939 
2940     if (simHitRefKeyToIndex.find(simHitKey) != simHitRefKeyToIndex.end()) {
2941       for (const auto& assoc2 : simHitsTPAssoc) {
2942         if (std::make_pair(assoc2.second.key(), assoc2.second.id()) == simHitKey) {
2943 #ifdef EDM_ML_DEBUG
2944           auto range1 = std::equal_range(simHitsTPAssoc.begin(),
2945                                          simHitsTPAssoc.end(),
2946                                          std::make_pair(assoc.first, TrackPSimHitRef()),
2947                                          SimHitTPAssociationProducer::simHitTPAssociationListGreater);
2948           auto range2 = std::equal_range(simHitsTPAssoc.begin(),
2949                                          simHitsTPAssoc.end(),
2950                                          std::make_pair(assoc2.first, TrackPSimHitRef()),
2951                                          SimHitTPAssociationProducer::simHitTPAssociationListGreater);
2952 
2953           LogTrace("TrackingNtuple") << "Earlier TP " << assoc2.first.key() << " SimTrack Ids";
2954           for (const auto& simTrack : assoc2.first->g4Tracks()) {
2955             edm::LogPrint("TrackingNtuple") << " SimTrack " << simTrack.trackId() << " BX:event "
2956                                             << simTrack.eventId().bunchCrossing() << ":" << simTrack.eventId().event();
2957           }
2958           for (auto iHit = range2.first; iHit != range2.second; ++iHit) {
2959             LogTrace("TrackingNtuple") << " SimHit " << iHit->second.key() << " " << iHit->second.id() << " tof "
2960                                        << iHit->second->tof() << " trackId " << iHit->second->trackId() << " BX:event "
2961                                        << iHit->second->eventId().bunchCrossing() << ":"
2962                                        << iHit->second->eventId().event();
2963           }
2964           LogTrace("TrackingNtuple") << "Current TP " << assoc.first.key() << " SimTrack Ids";
2965           for (const auto& simTrack : assoc.first->g4Tracks()) {
2966             edm::LogPrint("TrackingNtuple") << " SimTrack " << simTrack.trackId() << " BX:event "
2967                                             << simTrack.eventId().bunchCrossing() << ":" << simTrack.eventId().event();
2968           }
2969           for (auto iHit = range1.first; iHit != range1.second; ++iHit) {
2970             LogTrace("TrackingNtuple") << " SimHit " << iHit->second.key() << " " << iHit->second.id() << " tof "
2971                                        << iHit->second->tof() << " trackId " << iHit->second->trackId() << " BX:event "
2972                                        << iHit->second->eventId().bunchCrossing() << ":"
2973                                        << iHit->second->eventId().event();
2974           }
2975 #endif
2976 
2977           throw cms::Exception("LogicError")
2978               << "Got second time the SimHit " << simHitKey.first << " of " << simHitKey.second
2979               << ", first time with TrackingParticle " << assoc2.first.key() << ", now with " << tpKey;
2980         }
2981       }
2982       throw cms::Exception("LogicError") << "Got second time the SimHit " << simHitKey.first << " of "
2983                                          << simHitKey.second << ", now with TrackingParticle " << tpKey
2984                                          << ", but I didn't find the first occurrance!";
2985     }
2986 
2987     auto det = tracker.idToDetUnit(detId);
2988     if (!det)
2989       throw cms::Exception("LogicError") << "Did not find a det unit for DetId " << simhit.detUnitId()
2990                                          << " from tracker geometry";
2991 
2992     const auto pos = det->surface().toGlobal(simhit.localPosition());
2993     const float tof = simhit.timeOfFlight();
2994 
2995     const auto simHitIndex = simhit_x.size();
2996     simHitRefKeyToIndex[simHitKey] = simHitIndex;
2997 
2998     if (includeStripHits_)
2999       simhit_detId.push_back(tracker, tTopo, detId);
3000     else
3001       simhit_detId_phase2.push_back(tracker, tTopo, detId);
3002     simhit_x.push_back(pos.x());
3003     simhit_y.push_back(pos.y());
3004     simhit_z.push_back(pos.z());
3005     if (saveSimHitsP3_) {
3006       const auto mom = det->surface().toGlobal(simhit.momentumAtEntry());
3007       simhit_px.push_back(mom.x());
3008       simhit_py.push_back(mom.y());
3009       simhit_pz.push_back(mom.z());
3010     }
3011     simhit_particle.push_back(simhit.particleType());
3012     simhit_process.push_back(simhit.processType());
3013     simhit_eloss.push_back(simhit.energyLoss());
3014     simhit_tof.push_back(tof);
3015     //simhit_simTrackId.push_back(simhit.trackId());
3016 
3017     simhit_simTrkIdx.push_back(tpIndex);
3018 
3019     simhit_hitIdx.emplace_back();   // filled in matchCluster
3020     simhit_hitType.emplace_back();  // filled in matchCluster
3021 
3022     tpHitList.emplace_back(tpKey, simHitIndex, tof, simhit.detUnitId());
3023   }
3024 }
3025 
3026 void TrackingNtuple::fillPixelHits(const edm::Event& iEvent,
3027                                    const TrackerGeometry& tracker,
3028                                    const ClusterTPAssociation& clusterToTPMap,
3029                                    const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3030                                    const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
3031                                    const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
3032                                    const TrackerTopology& tTopo,
3033                                    const SimHitRefKeyToIndex& simHitRefKeyToIndex,
3034                                    std::set<edm::ProductID>& hitProductIds) {
3035   std::vector<std::pair<uint64_t, PixelMaskContainer const*>> pixelMasks;
3036   pixelMasks.reserve(pixelUseMaskTokens_.size());
3037   for (const auto& itoken : pixelUseMaskTokens_) {
3038     edm::Handle<PixelMaskContainer> aH;
3039     iEvent.getByToken(itoken.second, aH);
3040     pixelMasks.emplace_back(1 << itoken.first, aH.product());
3041   }
3042   auto pixUsedMask = [&pixelMasks](size_t key) {
3043     uint64_t mask = 0;
3044     for (auto const& m : pixelMasks) {
3045       if (m.second->mask(key))
3046         mask |= m.first;
3047     }
3048     return mask;
3049   };
3050 
3051   edm::Handle<SiPixelRecHitCollection> pixelHits;
3052   iEvent.getByToken(pixelRecHitToken_, pixelHits);
3053   for (auto it = pixelHits->begin(); it != pixelHits->end(); it++) {
3054     const DetId hitId = it->detId();
3055     for (auto hit = it->begin(); hit != it->end(); hit++) {
3056       hitProductIds.insert(hit->cluster().id());
3057 
3058       const int key = hit->cluster().key();
3059       const int lay = tTopo.layer(hitId);
3060 
3061       pix_isBarrel.push_back(hitId.subdetId() == 1);
3062       pix_detId.push_back(tracker, tTopo, hitId);
3063       pix_trkIdx.emplace_back();  // filled in fillTracks
3064       pix_onTrk_x.emplace_back();
3065       pix_onTrk_y.emplace_back();
3066       pix_onTrk_z.emplace_back();
3067       pix_onTrk_xx.emplace_back();
3068       pix_onTrk_xy.emplace_back();
3069       pix_onTrk_yy.emplace_back();
3070       pix_onTrk_yz.emplace_back();
3071       pix_onTrk_zz.emplace_back();
3072       pix_onTrk_zx.emplace_back();
3073       pix_tcandIdx.emplace_back();  // filled in fillCandidates
3074       pix_seeIdx.emplace_back();    // filled in fillSeeds
3075       pix_x.push_back(hit->globalPosition().x());
3076       pix_y.push_back(hit->globalPosition().y());
3077       pix_z.push_back(hit->globalPosition().z());
3078       pix_xx.push_back(hit->globalPositionError().cxx());
3079       pix_xy.push_back(hit->globalPositionError().cyx());
3080       pix_yy.push_back(hit->globalPositionError().cyy());
3081       pix_yz.push_back(hit->globalPositionError().czy());
3082       pix_zz.push_back(hit->globalPositionError().czz());
3083       pix_zx.push_back(hit->globalPositionError().czx());
3084       pix_radL.push_back(hit->surface()->mediumProperties().radLen());
3085       pix_bbxi.push_back(hit->surface()->mediumProperties().xi());
3086       pix_clustSizeCol.push_back(hit->cluster()->sizeY());
3087       pix_clustSizeRow.push_back(hit->cluster()->sizeX());
3088       pix_usedMask.push_back(pixUsedMask(hit->firstClusterRef().key()));
3089 
3090       LogTrace("TrackingNtuple") << "pixHit cluster=" << key << " subdId=" << hitId.subdetId() << " lay=" << lay
3091                                  << " rawId=" << hitId.rawId() << " pos =" << hit->globalPosition();
3092       if (includeTrackingParticles_) {
3093         SimHitData simHitData = matchCluster(hit->firstClusterRef(),
3094                                              hitId,
3095                                              key,
3096                                              *hit,
3097                                              clusterToTPMap,
3098                                              tpKeyToIndex,
3099                                              simHitsTPAssoc,
3100                                              digiSimLink,
3101                                              simHitRefKeyToIndex,
3102                                              HitType::Pixel);
3103         pix_simHitIdx.push_back(simHitData.matchingSimHit);
3104         pix_simType.push_back(static_cast<int>(simHitData.type));
3105         pix_xySignificance.push_back(simHitData.xySignificance);
3106         pix_chargeFraction.push_back(simHitData.chargeFraction);
3107         LogTrace("TrackingNtuple") << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
3108         if (!simHitData.matchingSimHit.empty()) {
3109           const auto simHitIdx = simHitData.matchingSimHit[0];
3110           LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx << " simHitPos="
3111                                      << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3112                                      << " energyLoss=" << simhit_eloss[simHitIdx]
3113                                      << " particleType=" << simhit_particle[simHitIdx]
3114                                      << " processType=" << simhit_process[simHitIdx]
3115                                      << " bunchCrossing=" << simHitData.bunchCrossing[0]
3116                                      << " event=" << simHitData.event[0];
3117         }
3118       }
3119     }
3120   }
3121 }
3122 
3123 void TrackingNtuple::fillStripRphiStereoHits(const edm::Event& iEvent,
3124                                              const TrackerGeometry& tracker,
3125                                              const ClusterTPAssociation& clusterToTPMap,
3126                                              const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3127                                              const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
3128                                              const edm::DetSetVector<StripDigiSimLink>& digiSimLink,
3129                                              const TrackerTopology& tTopo,
3130                                              const SimHitRefKeyToIndex& simHitRefKeyToIndex,
3131                                              std::set<edm::ProductID>& hitProductIds) {
3132   std::vector<std::pair<uint64_t, StripMaskContainer const*>> stripMasks;
3133   stripMasks.reserve(stripUseMaskTokens_.size());
3134   for (const auto& itoken : stripUseMaskTokens_) {
3135     edm::Handle<StripMaskContainer> aH;
3136     iEvent.getByToken(itoken.second, aH);
3137     stripMasks.emplace_back(1 << itoken.first, aH.product());
3138   }
3139   auto strUsedMask = [&stripMasks](size_t key) {
3140     uint64_t mask = 0;
3141     for (auto const& m : stripMasks) {
3142       if (m.second->mask(key))
3143         mask |= m.first;
3144     }
3145     return mask;
3146   };
3147 
3148   //index strip hit branches by cluster index
3149   edm::Handle<SiStripRecHit2DCollection> rphiHits;
3150   iEvent.getByToken(stripRphiRecHitToken_, rphiHits);
3151   edm::Handle<SiStripRecHit2DCollection> stereoHits;
3152   iEvent.getByToken(stripStereoRecHitToken_, stereoHits);
3153   int totalStripHits = rphiHits->dataSize() + stereoHits->dataSize();
3154   str_isBarrel.resize(totalStripHits);
3155   str_detId.resize(totalStripHits);
3156   str_trkIdx.resize(totalStripHits);  // filled in fillTracks
3157   str_onTrk_x.resize(totalStripHits);
3158   str_onTrk_y.resize(totalStripHits);
3159   str_onTrk_z.resize(totalStripHits);
3160   str_onTrk_xx.resize(totalStripHits);
3161   str_onTrk_xy.resize(totalStripHits);
3162   str_onTrk_yy.resize(totalStripHits);
3163   str_onTrk_yz.resize(totalStripHits);
3164   str_onTrk_zz.resize(totalStripHits);
3165   str_onTrk_zx.resize(totalStripHits);
3166   str_tcandIdx.resize(totalStripHits);  // filled in fillCandidates
3167   str_seeIdx.resize(totalStripHits);    // filled in fillSeeds
3168   str_simHitIdx.resize(totalStripHits);
3169   str_simType.resize(totalStripHits);
3170   str_chargeFraction.resize(totalStripHits);
3171   str_x.resize(totalStripHits);
3172   str_y.resize(totalStripHits);
3173   str_z.resize(totalStripHits);
3174   str_xx.resize(totalStripHits);
3175   str_xy.resize(totalStripHits);
3176   str_yy.resize(totalStripHits);
3177   str_yz.resize(totalStripHits);
3178   str_zz.resize(totalStripHits);
3179   str_zx.resize(totalStripHits);
3180   str_xySignificance.resize(totalStripHits);
3181   str_chargeFraction.resize(totalStripHits);
3182   str_radL.resize(totalStripHits);
3183   str_bbxi.resize(totalStripHits);
3184   str_chargePerCM.resize(totalStripHits);
3185   str_clustSize.resize(totalStripHits);
3186   str_usedMask.resize(totalStripHits);
3187 
3188   auto fill = [&](const SiStripRecHit2DCollection& hits, const char* name) {
3189     for (const auto& detset : hits) {
3190       const DetId hitId = detset.detId();
3191       for (const auto& hit : detset) {
3192         hitProductIds.insert(hit.cluster().id());
3193 
3194         const int key = hit.cluster().key();
3195         const int lay = tTopo.layer(hitId);
3196         str_isBarrel[key] = (hitId.subdetId() == StripSubdetector::TIB || hitId.subdetId() == StripSubdetector::TOB);
3197         str_detId.set(key, tracker, tTopo, hitId);
3198         str_x[key] = hit.globalPosition().x();
3199         str_y[key] = hit.globalPosition().y();
3200         str_z[key] = hit.globalPosition().z();
3201         str_xx[key] = hit.globalPositionError().cxx();
3202         str_xy[key] = hit.globalPositionError().cyx();
3203         str_yy[key] = hit.globalPositionError().cyy();
3204         str_yz[key] = hit.globalPositionError().czy();
3205         str_zz[key] = hit.globalPositionError().czz();
3206         str_zx[key] = hit.globalPositionError().czx();
3207         str_radL[key] = hit.surface()->mediumProperties().radLen();
3208         str_bbxi[key] = hit.surface()->mediumProperties().xi();
3209         str_chargePerCM[key] = siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster());
3210         str_clustSize[key] = hit.cluster()->amplitudes().size();
3211         str_usedMask[key] = strUsedMask(key);
3212         LogTrace("TrackingNtuple") << name << " cluster=" << key << " subdId=" << hitId.subdetId() << " lay=" << lay
3213                                    << " rawId=" << hitId.rawId() << " pos =" << hit.globalPosition();
3214 
3215         if (includeTrackingParticles_) {
3216           SimHitData simHitData = matchCluster(hit.firstClusterRef(),
3217                                                hitId,
3218                                                key,
3219                                                hit,
3220                                                clusterToTPMap,
3221                                                tpKeyToIndex,
3222                                                simHitsTPAssoc,
3223                                                digiSimLink,
3224                                                simHitRefKeyToIndex,
3225                                                HitType::Strip);
3226           str_simHitIdx[key] = simHitData.matchingSimHit;
3227           str_simType[key] = static_cast<int>(simHitData.type);
3228           str_xySignificance[key] = simHitData.xySignificance;
3229           str_chargeFraction[key] = simHitData.chargeFraction;
3230           LogTrace("TrackingNtuple") << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
3231           if (!simHitData.matchingSimHit.empty()) {
3232             const auto simHitIdx = simHitData.matchingSimHit[0];
3233             LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx << " simHitPos="
3234                                        << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3235                                        << " simHitPos="
3236                                        << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3237                                        << " energyLoss=" << simhit_eloss[simHitIdx]
3238                                        << " particleType=" << simhit_particle[simHitIdx]
3239                                        << " processType=" << simhit_process[simHitIdx]
3240                                        << " bunchCrossing=" << simHitData.bunchCrossing[0]
3241                                        << " event=" << simHitData.event[0];
3242           }
3243         }
3244       }
3245     }
3246   };
3247 
3248   fill(*rphiHits, "stripRPhiHit");
3249   fill(*stereoHits, "stripStereoHit");
3250 }
3251 
3252 size_t TrackingNtuple::addStripMatchedHit(const SiStripMatchedRecHit2D& hit,
3253                                           const TrackerGeometry& tracker,
3254                                           const TrackerTopology& tTopo,
3255                                           const std::vector<std::pair<uint64_t, StripMaskContainer const*>>& stripMasks,
3256                                           std::vector<std::pair<int, int>>& monoStereoClusterList) {
3257   auto strUsedMask = [&stripMasks](size_t key) {
3258     uint64_t mask = 0;
3259     for (auto const& m : stripMasks) {
3260       if (m.second->mask(key))
3261         mask |= m.first;
3262     }
3263     return mask;
3264   };
3265 
3266   const auto hitId = hit.geographicalId();
3267   const int lay = tTopo.layer(hitId);
3268   monoStereoClusterList.emplace_back(hit.monoHit().cluster().key(), hit.stereoHit().cluster().key());
3269   glu_isBarrel.push_back((hitId.subdetId() == StripSubdetector::TIB || hitId.subdetId() == StripSubdetector::TOB));
3270   glu_detId.push_back(tracker, tTopo, hitId);
3271   glu_monoIdx.push_back(hit.monoHit().cluster().key());
3272   glu_stereoIdx.push_back(hit.stereoHit().cluster().key());
3273   glu_seeIdx.emplace_back();  // filled in fillSeeds
3274   glu_x.push_back(hit.globalPosition().x());
3275   glu_y.push_back(hit.globalPosition().y());
3276   glu_z.push_back(hit.globalPosition().z());
3277   glu_xx.push_back(hit.globalPositionError().cxx());
3278   glu_xy.push_back(hit.globalPositionError().cyx());
3279   glu_yy.push_back(hit.globalPositionError().cyy());
3280   glu_yz.push_back(hit.globalPositionError().czy());
3281   glu_zz.push_back(hit.globalPositionError().czz());
3282   glu_zx.push_back(hit.globalPositionError().czx());
3283   glu_radL.push_back(hit.surface()->mediumProperties().radLen());
3284   glu_bbxi.push_back(hit.surface()->mediumProperties().xi());
3285   glu_chargePerCM.push_back(siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()));
3286   glu_clustSizeMono.push_back(hit.monoHit().cluster()->amplitudes().size());
3287   glu_clustSizeStereo.push_back(hit.stereoHit().cluster()->amplitudes().size());
3288   glu_usedMaskMono.push_back(strUsedMask(hit.monoHit().cluster().key()));
3289   glu_usedMaskStereo.push_back(strUsedMask(hit.stereoHit().cluster().key()));
3290   LogTrace("TrackingNtuple") << "stripMatchedHit"
3291                              << " cluster0=" << hit.stereoHit().cluster().key()
3292                              << " cluster1=" << hit.monoHit().cluster().key() << " subdId=" << hitId.subdetId()
3293                              << " lay=" << lay << " rawId=" << hitId.rawId() << " pos =" << hit.globalPosition();
3294   return glu_isBarrel.size() - 1;
3295 }
3296 
3297 void TrackingNtuple::fillStripMatchedHits(const edm::Event& iEvent,
3298                                           const TrackerGeometry& tracker,
3299                                           const TrackerTopology& tTopo,
3300                                           std::vector<std::pair<int, int>>& monoStereoClusterList) {
3301   std::vector<std::pair<uint64_t, StripMaskContainer const*>> stripMasks;
3302   stripMasks.reserve(stripUseMaskTokens_.size());
3303   for (const auto& itoken : stripUseMaskTokens_) {
3304     edm::Handle<StripMaskContainer> aH;
3305     iEvent.getByToken(itoken.second, aH);
3306     stripMasks.emplace_back(1 << itoken.first, aH.product());
3307   }
3308 
3309   edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
3310   iEvent.getByToken(stripMatchedRecHitToken_, matchedHits);
3311   for (auto it = matchedHits->begin(); it != matchedHits->end(); it++) {
3312     for (auto hit = it->begin(); hit != it->end(); hit++) {
3313       addStripMatchedHit(*hit, tracker, tTopo, stripMasks, monoStereoClusterList);
3314     }
3315   }
3316 }
3317 
3318 void TrackingNtuple::fillPhase2OTHits(const edm::Event& iEvent,
3319                                       const ClusterTPAssociation& clusterToTPMap,
3320                                       const TrackerGeometry& tracker,
3321                                       const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3322                                       const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
3323                                       const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
3324                                       const TrackerTopology& tTopo,
3325                                       const SimHitRefKeyToIndex& simHitRefKeyToIndex,
3326                                       std::set<edm::ProductID>& hitProductIds) {
3327   std::vector<std::pair<uint64_t, Phase2OTMaskContainer const*>> phase2OTMasks;
3328   phase2OTMasks.reserve(ph2OTUseMaskTokens_.size());
3329   for (const auto& itoken : ph2OTUseMaskTokens_) {
3330     edm::Handle<Phase2OTMaskContainer> aH;
3331     iEvent.getByToken(itoken.second, aH);
3332     phase2OTMasks.emplace_back(1 << itoken.first, aH.product());
3333   }
3334   auto ph2OTUsedMask = [&phase2OTMasks](size_t key) {
3335     uint64_t mask = 0;
3336     for (auto const& m : phase2OTMasks) {
3337       if (m.second->mask(key))
3338         mask |= m.first;
3339     }
3340     return mask;
3341   };
3342 
3343   edm::Handle<Phase2TrackerRecHit1DCollectionNew> phase2OTHits;
3344   iEvent.getByToken(phase2OTRecHitToken_, phase2OTHits);
3345   for (auto it = phase2OTHits->begin(); it != phase2OTHits->end(); it++) {
3346     const DetId hitId = it->detId();
3347     for (auto hit = it->begin(); hit != it->end(); hit++) {
3348       hitProductIds.insert(hit->cluster().id());
3349 
3350       const int key = hit->cluster().key();
3351       const int lay = tTopo.layer(hitId);
3352 
3353       ph2_isBarrel.push_back(hitId.subdetId() == 1);
3354       ph2_detId.push_back(tracker, tTopo, hitId);
3355       ph2_trkIdx.emplace_back();  // filled in fillTracks
3356       ph2_onTrk_x.emplace_back();
3357       ph2_onTrk_y.emplace_back();
3358       ph2_onTrk_z.emplace_back();
3359       ph2_onTrk_xx.emplace_back();
3360       ph2_onTrk_xy.emplace_back();
3361       ph2_onTrk_yy.emplace_back();
3362       ph2_onTrk_yz.emplace_back();
3363       ph2_onTrk_zz.emplace_back();
3364       ph2_onTrk_zx.emplace_back();
3365       ph2_tcandIdx.emplace_back();  // filled in fillCandidates
3366       ph2_seeIdx.emplace_back();    // filled in fillSeeds
3367       ph2_x.push_back(hit->globalPosition().x());
3368       ph2_y.push_back(hit->globalPosition().y());
3369       ph2_z.push_back(hit->globalPosition().z());
3370       ph2_xx.push_back(hit->globalPositionError().cxx());
3371       ph2_xy.push_back(hit->globalPositionError().cyx());
3372       ph2_yy.push_back(hit->globalPositionError().cyy());
3373       ph2_yz.push_back(hit->globalPositionError().czy());
3374       ph2_zz.push_back(hit->globalPositionError().czz());
3375       ph2_zx.push_back(hit->globalPositionError().czx());
3376       ph2_radL.push_back(hit->surface()->mediumProperties().radLen());
3377       ph2_bbxi.push_back(hit->surface()->mediumProperties().xi());
3378       ph2_usedMask.push_back(ph2OTUsedMask(hit->firstClusterRef().key()));
3379 
3380       LogTrace("TrackingNtuple") << "phase2 OT cluster=" << key << " subdId=" << hitId.subdetId() << " lay=" << lay
3381                                  << " rawId=" << hitId.rawId() << " pos =" << hit->globalPosition();
3382 
3383       if (includeTrackingParticles_) {
3384         SimHitData simHitData = matchCluster(hit->firstClusterRef(),
3385                                              hitId,
3386                                              key,
3387                                              *hit,
3388                                              clusterToTPMap,
3389                                              tpKeyToIndex,
3390                                              simHitsTPAssoc,
3391                                              digiSimLink,
3392                                              simHitRefKeyToIndex,
3393                                              HitType::Phase2OT);
3394         ph2_xySignificance.push_back(simHitData.xySignificance);
3395         ph2_simHitIdx.push_back(simHitData.matchingSimHit);
3396         ph2_simType.push_back(static_cast<int>(simHitData.type));
3397         LogTrace("TrackingNtuple") << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
3398         if (!simHitData.matchingSimHit.empty()) {
3399           const auto simHitIdx = simHitData.matchingSimHit[0];
3400           LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx << " simHitPos="
3401                                      << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
3402                                      << " energyLoss=" << simhit_eloss[simHitIdx]
3403                                      << " particleType=" << simhit_particle[simHitIdx]
3404                                      << " processType=" << simhit_process[simHitIdx]
3405                                      << " bunchCrossing=" << simHitData.bunchCrossing[0]
3406                                      << " event=" << simHitData.event[0];
3407         }
3408       }
3409     }
3410   }
3411 }
3412 
3413 void TrackingNtuple::fillSeeds(const edm::Event& iEvent,
3414                                const TrackingParticleRefVector& tpCollection,
3415                                const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3416                                const reco::BeamSpot& bs,
3417                                const TrackerGeometry& tracker,
3418                                const reco::TrackToTrackingParticleAssociator& associatorByHits,
3419                                const ClusterTPAssociation& clusterToTPMap,
3420                                const MagneticField& theMF,
3421                                const TrackerTopology& tTopo,
3422                                std::vector<std::pair<int, int>>& monoStereoClusterList,
3423                                const std::set<edm::ProductID>& hitProductIds,
3424                                std::map<edm::ProductID, size_t>& seedCollToOffset) {
3425   TSCBLBuilderNoMaterial tscblBuilder;
3426   for (size_t iColl = 0; iColl < seedTokens_.size(); ++iColl) {
3427     const auto& seedToken = seedTokens_[iColl];
3428 
3429     edm::Handle<edm::View<reco::Track>> seedTracksHandle;
3430     iEvent.getByToken(seedToken, seedTracksHandle);
3431     const auto& seedTracks = *seedTracksHandle;
3432 
3433     if (seedTracks.empty())
3434       continue;
3435 
3436     edm::EDConsumerBase::Labels labels;
3437     labelsForToken(seedToken, labels);
3438 
3439     const auto& seedStopInfoToken = seedStopInfoTokens_[iColl];
3440     edm::Handle<std::vector<SeedStopInfo>> seedStopInfoHandle;
3441     iEvent.getByToken(seedStopInfoToken, seedStopInfoHandle);
3442     const auto& seedStopInfos = *seedStopInfoHandle;
3443     if (seedTracks.size() != seedStopInfos.size()) {
3444       edm::EDConsumerBase::Labels labels2;
3445       labelsForToken(seedStopInfoToken, labels2);
3446 
3447       throw cms::Exception("LogicError") << "Got " << seedTracks.size() << " seeds, but " << seedStopInfos.size()
3448                                          << " seed stopping infos for collections " << labels.module << ", "
3449                                          << labels2.module;
3450     }
3451 
3452     std::vector<std::pair<uint64_t, StripMaskContainer const*>> stripMasks;
3453     stripMasks.reserve(stripUseMaskTokens_.size());
3454     for (const auto& itoken : stripUseMaskTokens_) {
3455       edm::Handle<StripMaskContainer> aH;
3456       iEvent.getByToken(itoken.second, aH);
3457       stripMasks.emplace_back(1 << itoken.first, aH.product());
3458     }
3459 
3460     // The associator interfaces really need to be fixed...
3461     edm::RefToBaseVector<reco::Track> seedTrackRefs;
3462     for (edm::View<reco::Track>::size_type i = 0; i < seedTracks.size(); ++i) {
3463       seedTrackRefs.push_back(seedTracks.refAt(i));
3464     }
3465     reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(seedTrackRefs, tpCollection);
3466     reco::SimToRecoCollection simRecColl = associatorByHits.associateSimToReco(seedTrackRefs, tpCollection);
3467 
3468     TString label = labels.module;
3469     //format label to match algoName
3470     label.ReplaceAll("seedTracks", "");
3471     label.ReplaceAll("Seeds", "");
3472     label.ReplaceAll("muonSeeded", "muonSeededStep");
3473     //for HLT seeds
3474     label.ReplaceAll("FromPixelTracks", "");
3475     label.ReplaceAll("PFLowPixel", "");
3476     label.ReplaceAll("hltDoubletRecovery", "pixelPairStep");  //random choice
3477     int algo = reco::TrackBase::algoByName(label.Data());
3478 
3479     edm::ProductID id = seedTracks[0].seedRef().id();
3480     const auto offset = see_fitok.size();
3481     auto inserted = seedCollToOffset.emplace(id, offset);
3482     if (!inserted.second)
3483       throw cms::Exception("Configuration")
3484           << "Trying to add seeds with ProductID " << id << " for a second time from collection " << labels.module
3485           << ", seed algo " << label << ". Typically this is caused by a configuration problem.";
3486     see_offset.push_back(offset);
3487 
3488     LogTrace("TrackingNtuple") << "NEW SEED LABEL: " << label << " size: " << seedTracks.size() << " algo=" << algo
3489                                << " ProductID " << id;
3490 
3491     for (size_t iSeed = 0; iSeed < seedTrackRefs.size(); ++iSeed) {
3492       const auto& seedTrackRef = seedTrackRefs[iSeed];
3493       const auto& seedTrack = *seedTrackRef;
3494       const auto& seedRef = seedTrack.seedRef();
3495       const auto& seed = *seedRef;
3496 
3497       const auto seedStopInfo = seedStopInfos[iSeed];
3498 
3499       if (seedRef.id() != id)
3500         throw cms::Exception("LogicError")
3501             << "All tracks in 'TracksFromSeeds' collection should point to seeds in the same collection. Now the "
3502                "element 0 had ProductID "
3503             << id << " while the element " << seedTrackRef.key() << " had " << seedTrackRef.id()
3504             << ". The source collection is " << labels.module << ".";
3505 
3506       std::vector<int> tpIdx;
3507       std::vector<float> sharedFraction;
3508       auto foundTPs = recSimColl.find(seedTrackRef);
3509       if (foundTPs != recSimColl.end()) {
3510         for (const auto& tpQuality : foundTPs->val) {
3511           tpIdx.push_back(tpKeyToIndex.at(tpQuality.first.key()));
3512           sharedFraction.push_back(tpQuality.second);
3513         }
3514       }
3515 
3516       // Search for a best-matching TrackingParticle for a seed
3517       const int nHits = seedTrack.numberOfValidHits();
3518       const auto clusters = track_associator::hitsToClusterRefs(
3519           seedTrack.recHitsBegin(),
3520           seedTrack.recHitsEnd());  // TODO: this function is called 3 times per track, try to reduce
3521       const int nClusters = clusters.size();
3522       const auto bestKeyCount = findBestMatchingTrackingParticle(seedTrack.recHits(), clusterToTPMap, tpKeyToIndex);
3523       const float bestShareFrac =
3524           nClusters > 0 ? static_cast<float>(bestKeyCount.countClusters) / static_cast<float>(nClusters) : 0;
3525       // Another way starting from the first hit of the seed
3526       const auto bestFirstHitKeyCount =
3527           findMatchingTrackingParticleFromFirstHit(seedTrack.recHits(), clusterToTPMap, tpKeyToIndex);
3528       const float bestFirstHitShareFrac =
3529           nClusters > 0 ? static_cast<float>(bestFirstHitKeyCount.countClusters) / static_cast<float>(nClusters) : 0;
3530 
3531       const bool seedFitOk = !trackFromSeedFitFailed(seedTrack);
3532       const int charge = seedTrack.charge();
3533       const float pt = seedFitOk ? seedTrack.pt() : 0;
3534       const float eta = seedFitOk ? seedTrack.eta() : 0;
3535       const float phi = seedFitOk ? seedTrack.phi() : 0;
3536 
3537       const auto seedIndex = see_fitok.size();
3538 
3539       see_fitok.push_back(seedFitOk);
3540 
3541       see_px.push_back(seedFitOk ? seedTrack.px() : 0);
3542       see_py.push_back(seedFitOk ? seedTrack.py() : 0);
3543       see_pz.push_back(seedFitOk ? seedTrack.pz() : 0);
3544       see_pt.push_back(pt);
3545       see_eta.push_back(eta);
3546       see_phi.push_back(phi);
3547       see_q.push_back(charge);
3548       see_nValid.push_back(nHits);
3549 
3550       see_dxy.push_back(seedFitOk ? seedTrack.dxy(bs.position()) : 0);
3551       see_dz.push_back(seedFitOk ? seedTrack.dz(bs.position()) : 0);
3552       see_ptErr.push_back(seedFitOk ? seedTrack.ptError() : 0);
3553       see_etaErr.push_back(seedFitOk ? seedTrack.etaError() : 0);
3554       see_phiErr.push_back(seedFitOk ? seedTrack.phiError() : 0);
3555       see_dxyErr.push_back(seedFitOk ? seedTrack.dxyError() : 0);
3556       see_dzErr.push_back(seedFitOk ? seedTrack.dzError() : 0);
3557       see_algo.push_back(algo);
3558       see_stopReason.push_back(seedStopInfo.stopReasonUC());
3559       see_nCands.push_back(seedStopInfo.candidatesPerSeed());
3560 
3561       const auto& state = seedTrack.seedRef()->startingState();
3562       const auto& pos = state.parameters().position();
3563       const auto& mom = state.parameters().momentum();
3564       see_statePt.push_back(state.pt());
3565       see_stateTrajX.push_back(pos.x());
3566       see_stateTrajY.push_back(pos.y());
3567       see_stateTrajPx.push_back(mom.x());
3568       see_stateTrajPy.push_back(mom.y());
3569       see_stateTrajPz.push_back(mom.z());
3570 
3571       ///the following is useful for analysis in global coords at seed hit surface
3572       const TrackingRecHit* lastRecHit = &*(seed.recHits().end() - 1);
3573       TrajectoryStateOnSurface tsos =
3574           trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), &theMF);
3575       auto const& stateGlobal = tsos.globalParameters();
3576       see_stateTrajGlbX.push_back(stateGlobal.position().x());
3577       see_stateTrajGlbY.push_back(stateGlobal.position().y());
3578       see_stateTrajGlbZ.push_back(stateGlobal.position().z());
3579       see_stateTrajGlbPx.push_back(stateGlobal.momentum().x());
3580       see_stateTrajGlbPy.push_back(stateGlobal.momentum().y());
3581       see_stateTrajGlbPz.push_back(stateGlobal.momentum().z());
3582       if (addSeedCurvCov_) {
3583         auto const& stateCcov = tsos.curvilinearError().matrix();
3584         std::vector<float> cov(15);
3585         auto covP = cov.begin();
3586         for (auto const val : stateCcov)
3587           *(covP++) = val;  //row-major
3588         see_stateCurvCov.push_back(std::move(cov));
3589       }
3590 
3591       see_trkIdx.push_back(-1);    // to be set correctly in fillTracks
3592       see_tcandIdx.push_back(-1);  // to be set correctly in fillCandidates
3593       if (includeTrackingParticles_) {
3594         see_simTrkIdx.push_back(tpIdx);
3595         see_simTrkShareFrac.push_back(sharedFraction);
3596         see_bestSimTrkIdx.push_back(bestKeyCount.key >= 0 ? tpKeyToIndex.at(bestKeyCount.key) : -1);
3597         see_bestFromFirstHitSimTrkIdx.push_back(
3598             bestFirstHitKeyCount.key >= 0 ? tpKeyToIndex.at(bestFirstHitKeyCount.key) : -1);
3599       } else {
3600         see_isTrue.push_back(!tpIdx.empty());
3601       }
3602       see_bestSimTrkShareFrac.push_back(bestShareFrac);
3603       see_bestFromFirstHitSimTrkShareFrac.push_back(bestFirstHitShareFrac);
3604 
3605       /// Hmm, the following could make sense instead of plain failing if propagation to beam line fails
3606       /*
3607       const TrackingRecHit* lastRecHit = &*(seed.recHits().second-1);
3608       TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( itSeed->startingState(), lastRecHit->surface(), &theMF);
3609       float pt  = state.globalParameters().momentum().perp();
3610       float eta = state.globalParameters().momentum().eta();
3611       float phi = state.globalParameters().momentum().phi();
3612       see_px      .push_back( state.globalParameters().momentum().x() );
3613       see_py      .push_back( state.globalParameters().momentum().y() );
3614       see_pz      .push_back( state.globalParameters().momentum().z() );
3615       */
3616 
3617       std::vector<int> hitIdx;
3618       std::vector<int> hitType;
3619 
3620       for (auto const& hit : seed.recHits()) {
3621         int subid = hit.geographicalId().subdetId();
3622         if (subid == (int)PixelSubdetector::PixelBarrel || subid == (int)PixelSubdetector::PixelEndcap) {
3623           const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&hit);
3624           const auto& clusterRef = bhit->firstClusterRef();
3625           const auto clusterKey = clusterRef.cluster_pixel().key();
3626           if (includeAllHits_) {
3627             checkProductID(hitProductIds, clusterRef.id(), "seed");
3628             pix_seeIdx[clusterKey].push_back(seedIndex);
3629           }
3630           hitIdx.push_back(clusterKey);
3631           hitType.push_back(static_cast<int>(HitType::Pixel));
3632         } else if (subid == (int)StripSubdetector::TOB || subid == (int)StripSubdetector::TID ||
3633                    subid == (int)StripSubdetector::TIB || subid == (int)StripSubdetector::TEC) {
3634           if (trackerHitRTTI::isMatched(hit)) {
3635             const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(&hit);
3636             if (includeAllHits_) {
3637               checkProductID(hitProductIds, matchedHit->monoClusterRef().id(), "seed");
3638               checkProductID(hitProductIds, matchedHit->stereoClusterRef().id(), "seed");
3639             }
3640             int monoIdx = matchedHit->monoClusterRef().key();
3641             int stereoIdx = matchedHit->stereoClusterRef().key();
3642 
3643             std::vector<std::pair<int, int>>::iterator pos =
3644                 find(monoStereoClusterList.begin(), monoStereoClusterList.end(), std::make_pair(monoIdx, stereoIdx));
3645             size_t gluedIndex = -1;
3646             if (pos != monoStereoClusterList.end()) {
3647               gluedIndex = std::distance(monoStereoClusterList.begin(), pos);
3648             } else {
3649               // We can encounter glued hits not in the input
3650               // SiStripMatchedRecHit2DCollection, e.g. via muon
3651               // outside-in seeds (or anything taking hits from
3652               // MeasurementTrackerEvent). So let's add them here.
3653               gluedIndex = addStripMatchedHit(*matchedHit, tracker, tTopo, stripMasks, monoStereoClusterList);
3654             }
3655 
3656             if (includeAllHits_)
3657               glu_seeIdx[gluedIndex].push_back(seedIndex);
3658             hitIdx.push_back(gluedIndex);
3659             hitType.push_back(static_cast<int>(HitType::Glued));
3660           } else {
3661             const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&hit);
3662             const auto& clusterRef = bhit->firstClusterRef();
3663             unsigned int clusterKey;
3664             if (clusterRef.isPhase2()) {
3665               clusterKey = clusterRef.cluster_phase2OT().key();
3666             } else {
3667               clusterKey = clusterRef.cluster_strip().key();
3668             }
3669 
3670             if (includeAllHits_) {
3671               checkProductID(hitProductIds, clusterRef.id(), "seed");
3672               if (clusterRef.isPhase2()) {
3673                 ph2_seeIdx[clusterKey].push_back(seedIndex);
3674               } else {
3675                 str_seeIdx[clusterKey].push_back(seedIndex);
3676               }
3677             }
3678 
3679             hitIdx.push_back(clusterKey);
3680             if (clusterRef.isPhase2()) {
3681               hitType.push_back(static_cast<int>(HitType::Phase2OT));
3682             } else {
3683               hitType.push_back(static_cast<int>(HitType::Strip));
3684             }
3685           }
3686         } else {
3687           LogTrace("TrackingNtuple") << " not pixel and not Strip detector";
3688         }
3689       }
3690       see_hitIdx.push_back(hitIdx);
3691       see_hitType.push_back(hitType);
3692       see_nPixel.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Pixel)));
3693       see_nGlued.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Glued)));
3694       see_nStrip.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Strip)));
3695       see_nPhase2OT.push_back(std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Phase2OT)));
3696       see_nCluster.push_back(nClusters);
3697       //the part below is not strictly needed
3698       float chi2 = -1;
3699       if (nHits == 2) {
3700         auto const recHit0 = seed.recHits().begin();
3701         auto const recHit1 = seed.recHits().begin() + 1;
3702         std::vector<GlobalPoint> gp(2);
3703         std::vector<GlobalError> ge(2);
3704         gp[0] = recHit0->globalPosition();
3705         ge[0] = recHit0->globalPositionError();
3706         gp[1] = recHit1->globalPosition();
3707         ge[1] = recHit1->globalPositionError();
3708         LogTrace("TrackingNtuple")
3709             << "seed " << seedTrackRef.key() << " pt=" << pt << " eta=" << eta << " phi=" << phi << " q=" << charge
3710             << " - PAIR - ids: " << recHit0->geographicalId().rawId() << " " << recHit1->geographicalId().rawId()
3711             << " hitpos: " << gp[0] << " " << gp[1] << " trans0: "
3712             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[0]->globalPosition()
3713                                                     : GlobalPoint(0, 0, 0))
3714             << " "
3715             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[1]->globalPosition()
3716                                                     : GlobalPoint(0, 0, 0))
3717             << " trans1: "
3718             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[0]->globalPosition()
3719                                                     : GlobalPoint(0, 0, 0))
3720             << " "
3721             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[1]->globalPosition()
3722                                                     : GlobalPoint(0, 0, 0))
3723             << " eta,phi: " << gp[0].eta() << "," << gp[0].phi();
3724       } else if (nHits == 3) {
3725         auto const recHit0 = seed.recHits().begin();
3726         auto const recHit1 = seed.recHits().begin() + 1;
3727         auto const recHit2 = seed.recHits().begin() + 2;
3728         declareDynArray(GlobalPoint, 4, gp);
3729         declareDynArray(GlobalError, 4, ge);
3730         declareDynArray(bool, 4, bl);
3731         gp[0] = recHit0->globalPosition();
3732         ge[0] = recHit0->globalPositionError();
3733         int subid0 = recHit0->geographicalId().subdetId();
3734         bl[0] = (subid0 == StripSubdetector::TIB || subid0 == StripSubdetector::TOB ||
3735                  subid0 == (int)PixelSubdetector::PixelBarrel);
3736         gp[1] = recHit1->globalPosition();
3737         ge[1] = recHit1->globalPositionError();
3738         int subid1 = recHit1->geographicalId().subdetId();
3739         bl[1] = (subid1 == StripSubdetector::TIB || subid1 == StripSubdetector::TOB ||
3740                  subid1 == (int)PixelSubdetector::PixelBarrel);
3741         gp[2] = recHit2->globalPosition();
3742         ge[2] = recHit2->globalPositionError();
3743         int subid2 = recHit2->geographicalId().subdetId();
3744         bl[2] = (subid2 == StripSubdetector::TIB || subid2 == StripSubdetector::TOB ||
3745                  subid2 == (int)PixelSubdetector::PixelBarrel);
3746         RZLine rzLine(gp, ge, bl);
3747         float seed_chi2 = rzLine.chi2();
3748         //float seed_pt = state.globalParameters().momentum().perp();
3749         float seed_pt = pt;
3750         LogTrace("TrackingNtuple")
3751             << "seed " << seedTrackRef.key() << " pt=" << pt << " eta=" << eta << " phi=" << phi << " q=" << charge
3752             << " - TRIPLET - ids: " << recHit0->geographicalId().rawId() << " " << recHit1->geographicalId().rawId()
3753             << " " << recHit2->geographicalId().rawId() << " hitpos: " << gp[0] << " " << gp[1] << " " << gp[2]
3754             << " trans0: "
3755             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[0]->globalPosition()
3756                                                     : GlobalPoint(0, 0, 0))
3757             << " "
3758             << (recHit0->transientHits().size() > 1 ? recHit0->transientHits()[1]->globalPosition()
3759                                                     : GlobalPoint(0, 0, 0))
3760             << " trans1: "
3761             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[0]->globalPosition()
3762                                                     : GlobalPoint(0, 0, 0))
3763             << " "
3764             << (recHit1->transientHits().size() > 1 ? recHit1->transientHits()[1]->globalPosition()
3765                                                     : GlobalPoint(0, 0, 0))
3766             << " trans2: "
3767             << (recHit2->transientHits().size() > 1 ? recHit2->transientHits()[0]->globalPosition()
3768                                                     : GlobalPoint(0, 0, 0))
3769             << " "
3770             << (recHit2->transientHits().size() > 1 ? recHit2->transientHits()[1]->globalPosition()
3771                                                     : GlobalPoint(0, 0, 0))
3772             << " local: "
3773             << recHit2->localPosition()
3774             //<< " tsos pos, mom: " << state.globalPosition()<<" "<<state.globalMomentum()
3775             << " eta,phi: " << gp[0].eta() << "," << gp[0].phi() << " pt,chi2: " << seed_pt << "," << seed_chi2;
3776         chi2 = seed_chi2;
3777       }
3778       see_chi2.push_back(chi2);
3779     }
3780 
3781     fillTrackingParticlesForSeeds(tpCollection, simRecColl, tpKeyToIndex, offset);
3782   }
3783 }
3784 
3785 void TrackingNtuple::fillTracks(const edm::RefToBaseVector<reco::Track>& tracks,
3786                                 const TrackerGeometry& tracker,
3787                                 const TrackingParticleRefVector& tpCollection,
3788                                 const TrackingParticleRefKeyToIndex& tpKeyToIndex,
3789                                 const TrackingParticleRefKeyToCount& tpKeyToClusterCount,
3790                                 const MagneticField& mf,
3791                                 const reco::BeamSpot& bs,
3792                                 const reco::VertexCollection& vertices,
3793                                 const reco::TrackToTrackingParticleAssociator& associatorByHits,
3794                                 const ClusterTPAssociation& clusterToTPMap,
3795                                 const TrackerTopology& tTopo,
3796                                 const std::set<edm::ProductID>& hitProductIds,
3797                                 const std::map<edm::ProductID, size_t>& seedCollToOffset,
3798                                 const std::vector<const MVACollection*>& mvaColls,
3799                                 const std::vector<const QualityMaskCollection*>& qualColls) {
3800   reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(tracks, tpCollection);
3801   edm::EDConsumerBase::Labels labels;
3802   labelsForToken(trackToken_, labels);
3803   LogTrace("TrackingNtuple") << "NEW TRACK LABEL: " << labels.module;
3804 
3805   auto pvPosition = vertices[0].position();
3806 
3807   for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) {
3808     const auto& itTrack = tracks[iTrack];
3809     int charge = itTrack->charge();
3810     float pt = itTrack->pt();
3811     float eta = itTrack->eta();
3812     const double lambda = itTrack->lambda();
3813     float chi2 = itTrack->normalizedChi2();
3814     float ndof = itTrack->ndof();
3815     float phi = itTrack->phi();
3816     int nHits = itTrack->numberOfValidHits();
3817     const reco::HitPattern& hp = itTrack->hitPattern();
3818 
3819     const auto& tkParam = itTrack->parameters();
3820     auto tkCov = itTrack->covariance();
3821     tkCov.Invert();
3822 
3823     // Standard track-TP matching
3824     int nSimHits = 0;
3825     bool isSimMatched = false;
3826     std::vector<int> tpIdx;
3827     std::vector<float> sharedFraction;
3828     std::vector<float> tpChi2;
3829     auto foundTPs = recSimColl.find(itTrack);
3830     if (foundTPs != recSimColl.end()) {
3831       if (!foundTPs->val.empty()) {
3832         nSimHits = foundTPs->val[0].first->numberOfTrackerHits();
3833         isSimMatched = true;
3834       }
3835       for (const auto& tpQuality : foundTPs->val) {
3836         tpIdx.push_back(tpKeyToIndex.at(tpQuality.first.key()));
3837         sharedFraction.push_back(tpQuality.second);
3838         tpChi2.push_back(track_associator::trackAssociationChi2(tkParam, tkCov, *(tpCollection[tpIdx.back()]), mf, bs));
3839       }
3840     }
3841 
3842     // Search for a best-matching TrackingParticle for a track
3843     const auto clusters = track_associator::hitsToClusterRefs(
3844         itTrack->recHitsBegin(),
3845         itTrack->recHitsEnd());  // TODO: this function is called 3 times per track, try to reduce
3846     const int nClusters = clusters.size();
3847 
3848     const auto bestKeyCount = findBestMatchingTrackingParticle(itTrack->recHits(), clusterToTPMap, tpKeyToIndex);
3849     const float bestShareFrac = static_cast<float>(bestKeyCount.countClusters) / static_cast<float>(nClusters);
3850     float bestShareFracSimDenom = 0;
3851     float bestShareFracSimClusterDenom = 0;
3852     float bestChi2 = -1;
3853     if (bestKeyCount.key >= 0) {
3854       bestShareFracSimDenom =
3855           static_cast<float>(bestKeyCount.countClusters) /
3856           static_cast<float>(tpCollection[tpKeyToIndex.at(bestKeyCount.key)]->numberOfTrackerHits());
3857       bestShareFracSimClusterDenom =
3858           static_cast<float>(bestKeyCount.countClusters) / static_cast<float>(tpKeyToClusterCount.at(bestKeyCount.key));
3859       bestChi2 = track_associator::trackAssociationChi2(
3860           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestKeyCount.key)]), mf, bs);
3861     }
3862     // Another way starting from the first hit of the track
3863     const auto bestFirstHitKeyCount =
3864         findMatchingTrackingParticleFromFirstHit(itTrack->recHits(), clusterToTPMap, tpKeyToIndex);
3865     const float bestFirstHitShareFrac =
3866         static_cast<float>(bestFirstHitKeyCount.countClusters) / static_cast<float>(nClusters);
3867     float bestFirstHitShareFracSimDenom = 0;
3868     float bestFirstHitShareFracSimClusterDenom = 0;
3869     float bestFirstHitChi2 = -1;
3870     if (bestFirstHitKeyCount.key >= 0) {
3871       bestFirstHitShareFracSimDenom =
3872           static_cast<float>(bestFirstHitKeyCount.countClusters) /
3873           static_cast<float>(tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]->numberOfTrackerHits());
3874       bestFirstHitShareFracSimClusterDenom = static_cast<float>(bestFirstHitKeyCount.countClusters) /
3875                                              static_cast<float>(tpKeyToClusterCount.at(bestFirstHitKeyCount.key));
3876       bestFirstHitChi2 = track_associator::trackAssociationChi2(
3877           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]), mf, bs);
3878     }
3879 
3880     float chi2_1Dmod = chi2;
3881     int count1dhits = 0;
3882     for (auto iHit = itTrack->recHitsBegin(), iEnd = itTrack->recHitsEnd(); iHit != iEnd; ++iHit) {
3883       const TrackingRecHit& hit = **iHit;
3884       if (hit.isValid() && typeid(hit) == typeid(SiStripRecHit1D))
3885         ++count1dhits;
3886     }
3887     if (count1dhits > 0) {
3888       chi2_1Dmod = (chi2 + count1dhits) / (ndof + count1dhits);
3889     }
3890 
3891     Point bestPV = getBestVertex(*itTrack, vertices);
3892 
3893     trk_px.push_back(itTrack->px());
3894     trk_py.push_back(itTrack->py());
3895     trk_pz.push_back(itTrack->pz());
3896     trk_pt.push_back(pt);
3897     trk_inner_px.push_back(itTrack->innerMomentum().x());
3898     trk_inner_py.push_back(itTrack->innerMomentum().y());
3899     trk_inner_pz.push_back(itTrack->innerMomentum().z());
3900     trk_inner_pt.push_back(itTrack->innerMomentum().rho());
3901     trk_outer_px.push_back(itTrack->outerMomentum().x());
3902     trk_outer_py.push_back(itTrack->outerMomentum().y());
3903     trk_outer_pz.push_back(itTrack->outerMomentum().z());
3904     trk_outer_pt.push_back(itTrack->outerMomentum().rho());
3905     trk_eta.push_back(eta);
3906     trk_lambda.push_back(lambda);
3907     trk_cotTheta.push_back(1 / tan(M_PI * 0.5 - lambda));
3908     trk_phi.push_back(phi);
3909     trk_dxy.push_back(itTrack->dxy(bs.position()));
3910     trk_dz.push_back(itTrack->dz(bs.position()));
3911     trk_dxyPV.push_back(itTrack->dxy(pvPosition));
3912     trk_dzPV.push_back(itTrack->dz(pvPosition));
3913     trk_dxyClosestPV.push_back(itTrack->dxy(bestPV));
3914     trk_dzClosestPV.push_back(itTrack->dz(bestPV));
3915     trk_ptErr.push_back(itTrack->ptError());
3916     trk_etaErr.push_back(itTrack->etaError());
3917     trk_lambdaErr.push_back(itTrack->lambdaError());
3918     trk_phiErr.push_back(itTrack->phiError());
3919     trk_dxyErr.push_back(itTrack->dxyError());
3920     trk_dzErr.push_back(itTrack->dzError());
3921     trk_refpoint_x.push_back(itTrack->vx());
3922     trk_refpoint_y.push_back(itTrack->vy());
3923     trk_refpoint_z.push_back(itTrack->vz());
3924     trk_nChi2.push_back(chi2);
3925     trk_nChi2_1Dmod.push_back(chi2_1Dmod);
3926     trk_ndof.push_back(ndof);
3927     trk_q.push_back(charge);
3928     trk_nValid.push_back(hp.numberOfValidHits());
3929     trk_nLost.push_back(hp.numberOfLostHits(reco::HitPattern::TRACK_HITS));
3930     trk_nInactive.push_back(hp.trackerLayersTotallyOffOrBad(reco::HitPattern::TRACK_HITS));
3931     trk_nPixel.push_back(hp.numberOfValidPixelHits());
3932     trk_nStrip.push_back(hp.numberOfValidStripHits());
3933     trk_nOuterLost.push_back(hp.numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS));
3934     trk_nInnerLost.push_back(hp.numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS));
3935     trk_nOuterInactive.push_back(hp.trackerLayersTotallyOffOrBad(reco::HitPattern::MISSING_OUTER_HITS));
3936     trk_nInnerInactive.push_back(hp.trackerLayersTotallyOffOrBad(reco::HitPattern::MISSING_INNER_HITS));
3937     trk_nPixelLay.push_back(hp.pixelLayersWithMeasurement());
3938     trk_nStripLay.push_back(hp.stripLayersWithMeasurement());
3939     trk_n3DLay.push_back(hp.numberOfValidStripLayersWithMonoAndStereo() + hp.pixelLayersWithMeasurement());
3940     trk_nLostLay.push_back(hp.trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS));
3941     trk_nCluster.push_back(nClusters);
3942     trk_algo.push_back(itTrack->algo());
3943     trk_originalAlgo.push_back(itTrack->originalAlgo());
3944     trk_algoMask.push_back(itTrack->algoMaskUL());
3945     trk_stopReason.push_back(itTrack->stopReason());
3946     trk_isHP.push_back(itTrack->quality(reco::TrackBase::highPurity));
3947     if (includeMVA_) {
3948       for (size_t i = 0; i < trk_mvas.size(); ++i) {
3949         trk_mvas[i].push_back((*(mvaColls[i]))[iTrack]);
3950         trk_qualityMasks[i].push_back((*(qualColls[i]))[iTrack]);
3951       }
3952     }
3953     if (includeSeeds_) {
3954       auto offset = seedCollToOffset.find(itTrack->seedRef().id());
3955       if (offset == seedCollToOffset.end()) {
3956         throw cms::Exception("Configuration")
3957             << "Track algo '" << reco::TrackBase::algoName(itTrack->algo()) << "' originalAlgo '"
3958             << reco::TrackBase::algoName(itTrack->originalAlgo()) << "' refers to seed collection "
3959             << itTrack->seedRef().id()
3960             << ", but that seed collection is not given as an input. The following collections were given as an input "
3961             << make_ProductIDMapPrinter(seedCollToOffset);
3962       }
3963 
3964       const auto seedIndex = offset->second + itTrack->seedRef().key();
3965       trk_seedIdx.push_back(seedIndex);
3966       if (see_trkIdx[seedIndex] != -1) {
3967         throw cms::Exception("LogicError") << "Track index has already been set for seed " << seedIndex << " to "
3968                                            << see_trkIdx[seedIndex] << "; was trying to set it to " << iTrack;
3969       }
3970       see_trkIdx[seedIndex] = iTrack;
3971     }
3972     trk_vtxIdx.push_back(-1);  // to be set correctly in fillVertices
3973     if (includeTrackingParticles_) {
3974       trk_simTrkIdx.push_back(tpIdx);
3975       trk_simTrkShareFrac.push_back(sharedFraction);
3976       trk_simTrkNChi2.push_back(tpChi2);
3977       trk_bestSimTrkIdx.push_back(bestKeyCount.key >= 0 ? tpKeyToIndex.at(bestKeyCount.key) : -1);
3978       trk_bestFromFirstHitSimTrkIdx.push_back(bestFirstHitKeyCount.key >= 0 ? tpKeyToIndex.at(bestFirstHitKeyCount.key)
3979                                                                             : -1);
3980     } else {
3981       trk_isTrue.push_back(!tpIdx.empty());
3982     }
3983     trk_bestSimTrkShareFrac.push_back(bestShareFrac);
3984     trk_bestSimTrkShareFracSimDenom.push_back(bestShareFracSimDenom);
3985     trk_bestSimTrkShareFracSimClusterDenom.push_back(bestShareFracSimClusterDenom);
3986     trk_bestSimTrkNChi2.push_back(bestChi2);
3987     trk_bestFromFirstHitSimTrkShareFrac.push_back(bestFirstHitShareFrac);
3988     trk_bestFromFirstHitSimTrkShareFracSimDenom.push_back(bestFirstHitShareFracSimDenom);
3989     trk_bestFromFirstHitSimTrkShareFracSimClusterDenom.push_back(bestFirstHitShareFracSimClusterDenom);
3990     trk_bestFromFirstHitSimTrkNChi2.push_back(bestFirstHitChi2);
3991 
3992     LogTrace("TrackingNtuple") << "Track #" << itTrack.key() << " with q=" << charge << ", pT=" << pt
3993                                << " GeV, eta: " << eta << ", phi: " << phi << ", chi2=" << chi2 << ", Nhits=" << nHits
3994                                << ", algo=" << itTrack->algoName(itTrack->algo()).c_str()
3995                                << " hp=" << itTrack->quality(reco::TrackBase::highPurity)
3996                                << " seed#=" << itTrack->seedRef().key() << " simMatch=" << isSimMatched
3997                                << " nSimHits=" << nSimHits
3998                                << " sharedFraction=" << (sharedFraction.empty() ? -1 : sharedFraction[0])
3999                                << " tpIdx=" << (tpIdx.empty() ? -1 : tpIdx[0]);
4000     std::vector<int> hitIdx;
4001     std::vector<int> hitType;
4002 
4003     for (auto i = itTrack->recHitsBegin(); i != itTrack->recHitsEnd(); i++) {
4004       auto hit = *i;
4005       DetId hitId = hit->geographicalId();
4006       LogTrace("TrackingNtuple") << "hit #" << std::distance(itTrack->recHitsBegin(), i)
4007                                  << " subdet=" << hitId.subdetId();
4008       if (hitId.det() != DetId::Tracker)
4009         continue;
4010 
4011       LogTrace("TrackingNtuple") << " " << subdetstring(hitId.subdetId()) << " " << tTopo.layer(hitId);
4012 
4013       if (hit->isValid()) {
4014         //ugly... but works
4015         const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*hit);
4016         const auto& clusterRef = bhit->firstClusterRef();
4017         unsigned int clusterKey;
4018         if (clusterRef.isPixel()) {
4019           clusterKey = clusterRef.cluster_pixel().key();
4020         } else if (clusterRef.isPhase2()) {
4021           clusterKey = clusterRef.cluster_phase2OT().key();
4022         } else {
4023           clusterKey = clusterRef.cluster_strip().key();
4024         }
4025 
4026         LogTrace("TrackingNtuple") << " id: " << hitId.rawId() << " - globalPos =" << hit->globalPosition()
4027                                    << " cluster=" << clusterKey << " clusterRef ID=" << clusterRef.id()
4028                                    << " eta,phi: " << hit->globalPosition().eta() << "," << hit->globalPosition().phi();
4029         if (includeAllHits_) {
4030           checkProductID(hitProductIds, clusterRef.id(), "track");
4031           if (clusterRef.isPixel()) {
4032             pix_trkIdx[clusterKey].push_back(iTrack);
4033             pix_onTrk_x[clusterKey].push_back(hit->globalPosition().x());
4034             pix_onTrk_y[clusterKey].push_back(hit->globalPosition().y());
4035             pix_onTrk_z[clusterKey].push_back(hit->globalPosition().z());
4036             pix_onTrk_xx[clusterKey].push_back(hit->globalPositionError().cxx());
4037             pix_onTrk_xy[clusterKey].push_back(hit->globalPositionError().cyx());
4038             pix_onTrk_yy[clusterKey].push_back(hit->globalPositionError().cyy());
4039             pix_onTrk_yz[clusterKey].push_back(hit->globalPositionError().czy());
4040             pix_onTrk_zz[clusterKey].push_back(hit->globalPositionError().czz());
4041             pix_onTrk_zx[clusterKey].push_back(hit->globalPositionError().czx());
4042           } else if (clusterRef.isPhase2()) {
4043             ph2_trkIdx[clusterKey].push_back(iTrack);
4044             ph2_onTrk_x[clusterKey].push_back(hit->globalPosition().x());
4045             ph2_onTrk_y[clusterKey].push_back(hit->globalPosition().y());
4046             ph2_onTrk_z[clusterKey].push_back(hit->globalPosition().z());
4047             ph2_onTrk_xx[clusterKey].push_back(hit->globalPositionError().cxx());
4048             ph2_onTrk_xy[clusterKey].push_back(hit->globalPositionError().cyx());
4049             ph2_onTrk_yy[clusterKey].push_back(hit->globalPositionError().cyy());
4050             ph2_onTrk_yz[clusterKey].push_back(hit->globalPositionError().czy());
4051             ph2_onTrk_zz[clusterKey].push_back(hit->globalPositionError().czz());
4052             ph2_onTrk_zx[clusterKey].push_back(hit->globalPositionError().czx());
4053           } else {
4054             str_trkIdx[clusterKey].push_back(iTrack);
4055             str_onTrk_x[clusterKey].push_back(hit->globalPosition().x());
4056             str_onTrk_y[clusterKey].push_back(hit->globalPosition().y());
4057             str_onTrk_z[clusterKey].push_back(hit->globalPosition().z());
4058             str_onTrk_xx[clusterKey].push_back(hit->globalPositionError().cxx());
4059             str_onTrk_xy[clusterKey].push_back(hit->globalPositionError().cyx());
4060             str_onTrk_yy[clusterKey].push_back(hit->globalPositionError().cyy());
4061             str_onTrk_yz[clusterKey].push_back(hit->globalPositionError().czy());
4062             str_onTrk_zz[clusterKey].push_back(hit->globalPositionError().czz());
4063             str_onTrk_zx[clusterKey].push_back(hit->globalPositionError().czx());
4064           }
4065         }
4066 
4067         hitIdx.push_back(clusterKey);
4068         if (clusterRef.isPixel()) {
4069           hitType.push_back(static_cast<int>(HitType::Pixel));
4070         } else if (clusterRef.isPhase2()) {
4071           hitType.push_back(static_cast<int>(HitType::Phase2OT));
4072         } else {
4073           hitType.push_back(static_cast<int>(HitType::Strip));
4074         }
4075       } else {
4076         LogTrace("TrackingNtuple") << " - invalid hit";
4077 
4078         hitIdx.push_back(inv_isBarrel.size());
4079         hitType.push_back(static_cast<int>(HitType::Invalid));
4080 
4081         inv_isBarrel.push_back(hitId.subdetId() == 1);
4082         if (includeStripHits_)
4083           inv_detId.push_back(tracker, tTopo, hitId);
4084         else
4085           inv_detId_phase2.push_back(tracker, tTopo, hitId);
4086         inv_type.push_back(hit->getType());
4087       }
4088     }
4089 
4090     trk_hitIdx.push_back(hitIdx);
4091     trk_hitType.push_back(hitType);
4092   }
4093 }
4094 
4095 void TrackingNtuple::fillCandidates(const edm::Handle<TrackCandidateCollection>& candsHandle,
4096                                     int algo,
4097                                     const TrackingParticleRefVector& tpCollection,
4098                                     const TrackingParticleRefKeyToIndex& tpKeyToIndex,
4099                                     const TrackingParticleRefKeyToCount& tpKeyToClusterCount,
4100                                     const MagneticField& mf,
4101                                     const reco::BeamSpot& bs,
4102                                     const reco::VertexCollection& vertices,
4103                                     const reco::TrackToTrackingParticleAssociator& associatorByHits,
4104                                     const ClusterTPAssociation& clusterToTPMap,
4105                                     const TrackerGeometry& tracker,
4106                                     const TrackerTopology& tTopo,
4107                                     const std::set<edm::ProductID>& hitProductIds,
4108                                     const std::map<edm::ProductID, size_t>& seedCollToOffset) {
4109   reco::RecoToSimCollectionTCandidate recSimColl = associatorByHits.associateRecoToSim(candsHandle, tpCollection);
4110 
4111   TSCBLBuilderNoMaterial tscblBuilder;
4112 
4113   auto const& cands = *candsHandle;
4114   for (size_t iCand = 0; iCand < cands.size(); ++iCand) {
4115     const auto& aCand = cands[iCand];
4116     const edm::Ref<TrackCandidateCollection> aCandRef(candsHandle, iCand);
4117 
4118     //get parameters and errors from the candidate state
4119     auto const& pState = aCand.trajectoryStateOnDet();
4120     TrajectoryStateOnSurface state =
4121         trajectoryStateTransform::transientState(pState, &(tracker.idToDet(pState.detId())->surface()), &mf);
4122     TrajectoryStateClosestToBeamLine tbStateAtPCA = tscblBuilder(*state.freeState(), bs);
4123     if (!tbStateAtPCA.isValid()) {
4124       edm::LogVerbatim("TrackBuilding") << "TrajectoryStateClosestToBeamLine not valid";
4125     }
4126 
4127     auto const& stateAtPCA = tbStateAtPCA.trackStateAtPCA();
4128     auto v0 = stateAtPCA.position();
4129     auto p = stateAtPCA.momentum();
4130     math::XYZPoint pos(v0.x(), v0.y(), v0.z());
4131     math::XYZVector mom(p.x(), p.y(), p.z());
4132 
4133     //pseduo track for access to easy methods
4134     static const reco::Track::CovarianceMatrix dummyCov = AlgebraicMatrixID();
4135     reco::Track trk(
4136         0, 0, pos, mom, stateAtPCA.charge(), tbStateAtPCA.isValid() ? stateAtPCA.curvilinearError().matrix() : dummyCov);
4137 
4138     const auto& tkParam = trk.parameters();
4139     auto tkCov = trk.covariance();
4140     tkCov.Invert();
4141 
4142     // Standard track-TP matching
4143     int nSimHits = 0;
4144     bool isSimMatched = false;
4145     std::vector<int> tpIdx;
4146     std::vector<float> sharedFraction;
4147     std::vector<float> tpChi2;
4148     auto foundTPs = recSimColl.find(aCandRef);
4149     if (foundTPs != recSimColl.end()) {
4150       if (!foundTPs->val.empty()) {
4151         nSimHits = foundTPs->val[0].first->numberOfTrackerHits();
4152         isSimMatched = true;
4153       }
4154       for (const auto& tpQuality : foundTPs->val) {
4155         tpIdx.push_back(tpKeyToIndex.at(tpQuality.first.key()));
4156         sharedFraction.push_back(tpQuality.second);
4157         tpChi2.push_back(track_associator::trackAssociationChi2(tkParam, tkCov, *(tpCollection[tpIdx.back()]), mf, bs));
4158       }
4159     }
4160 
4161     // Search for a best-matching TrackingParticle for a track
4162     const auto clusters = track_associator::hitsToClusterRefs(aCand.recHits().begin(), aCand.recHits().end());
4163     const int nClusters = clusters.size();
4164 
4165     const auto bestKeyCount = findBestMatchingTrackingParticle(aCand.recHits(), clusterToTPMap, tpKeyToIndex);
4166     const float bestCountF = bestKeyCount.countClusters;
4167     const float bestShareFrac = bestCountF / nClusters;
4168     float bestShareFracSimDenom = 0;
4169     float bestShareFracSimClusterDenom = 0;
4170     float bestChi2 = -1;
4171     if (bestKeyCount.key >= 0) {
4172       bestShareFracSimDenom = bestCountF / tpCollection[tpKeyToIndex.at(bestKeyCount.key)]->numberOfTrackerHits();
4173       bestShareFracSimClusterDenom = bestCountF / tpKeyToClusterCount.at(bestKeyCount.key);
4174       bestChi2 = track_associator::trackAssociationChi2(
4175           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestKeyCount.key)]), mf, bs);
4176     }
4177     // Another way starting from the first hit of the track
4178     const auto bestFirstHitKeyCount =
4179         findMatchingTrackingParticleFromFirstHit(aCand.recHits(), clusterToTPMap, tpKeyToIndex);
4180     const float bestFirstCountF = bestFirstHitKeyCount.countClusters;
4181     const float bestFirstHitShareFrac = bestFirstCountF / nClusters;
4182     float bestFirstHitShareFracSimDenom = 0;
4183     float bestFirstHitShareFracSimClusterDenom = 0;
4184     float bestFirstHitChi2 = -1;
4185     if (bestFirstHitKeyCount.key >= 0) {
4186       bestFirstHitShareFracSimDenom =
4187           bestFirstCountF / tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]->numberOfTrackerHits();
4188       bestFirstHitShareFracSimClusterDenom = bestFirstCountF / tpKeyToClusterCount.at(bestFirstHitKeyCount.key);
4189       bestFirstHitChi2 = track_associator::trackAssociationChi2(
4190           tkParam, tkCov, *(tpCollection[tpKeyToIndex.at(bestFirstHitKeyCount.key)]), mf, bs);
4191     }
4192 
4193     auto iglobCand = tcand_pca_valid.size();  //global cand index
4194     tcand_pca_valid.push_back(tbStateAtPCA.isValid());
4195     tcand_pca_px.push_back(trk.px());
4196     tcand_pca_py.push_back(trk.py());
4197     tcand_pca_pz.push_back(trk.pz());
4198     tcand_pca_pt.push_back(trk.pt());
4199     tcand_pca_eta.push_back(trk.eta());
4200     tcand_pca_phi.push_back(trk.phi());
4201     tcand_pca_dxy.push_back(trk.dxy());
4202     tcand_pca_dz.push_back(trk.dz());
4203     tcand_pca_ptErr.push_back(trk.ptError());
4204     tcand_pca_etaErr.push_back(trk.etaError());
4205     tcand_pca_lambdaErr.push_back(trk.lambdaError());
4206     tcand_pca_phiErr.push_back(trk.phiError());
4207     tcand_pca_dxyErr.push_back(trk.dxyError());
4208     tcand_pca_dzErr.push_back(trk.dzError());
4209     tcand_px.push_back(state.globalMomentum().x());
4210     tcand_py.push_back(state.globalMomentum().y());
4211     tcand_pz.push_back(state.globalMomentum().z());
4212     tcand_pt.push_back(state.globalMomentum().perp());
4213     tcand_x.push_back(state.globalPosition().x());
4214     tcand_y.push_back(state.globalPosition().y());
4215     tcand_z.push_back(state.globalPosition().z());
4216 
4217     auto const& pStateCov = state.curvilinearError().matrix();
4218     tcand_qbpErr.push_back(sqrt(pStateCov(0, 0)));
4219     tcand_lambdaErr.push_back(sqrt(pStateCov(1, 1)));
4220     tcand_phiErr.push_back(sqrt(pStateCov(2, 2)));
4221     tcand_xtErr.push_back(sqrt(pStateCov(3, 3)));
4222     tcand_ytErr.push_back(sqrt(pStateCov(4, 4)));
4223 
4224     int ndof = -5;
4225     int nValid = 0;
4226     int nPixel = 0;
4227     int nStrip = 0;
4228     for (auto const& hit : aCand.recHits()) {
4229       if (hit.isValid()) {
4230         ndof += hit.dimension();
4231         nValid++;
4232 
4233         auto const subdet = hit.geographicalId().subdetId();
4234         if (subdet == PixelSubdetector::PixelBarrel || subdet == PixelSubdetector::PixelEndcap)
4235           nPixel++;
4236         else
4237           nStrip++;
4238       }
4239     }
4240     tcand_ndof.push_back(ndof);
4241     tcand_q.push_back(trk.charge());
4242     tcand_nValid.push_back(nValid);
4243     tcand_nPixel.push_back(nPixel);
4244     tcand_nStrip.push_back(nStrip);
4245     tcand_nCluster.push_back(nClusters);
4246     tcand_algo.push_back(algo);
4247     tcand_stopReason.push_back(aCand.stopReason());
4248     if (includeSeeds_) {
4249       auto offset = seedCollToOffset.find(aCand.seedRef().id());
4250       if (offset == seedCollToOffset.end()) {
4251         throw cms::Exception("Configuration")
4252             << "Track candidate refers to seed collection " << aCand.seedRef().id()
4253             << ", but that seed collection is not given as an input. The following collections were given as an input "
4254             << make_ProductIDMapPrinter(seedCollToOffset);
4255       }
4256 
4257       const auto seedIndex = offset->second + aCand.seedRef().key();
4258       tcand_seedIdx.push_back(seedIndex);
4259       if (see_tcandIdx[seedIndex] != -1) {
4260         throw cms::Exception("LogicError")
4261             << "Track cand index has already been set for seed " << seedIndex << " to " << see_tcandIdx[seedIndex]
4262             << "; was trying to set it to " << iglobCand << " current " << iCand;
4263       }
4264       see_tcandIdx[seedIndex] = iglobCand;
4265     }
4266     tcand_vtxIdx.push_back(-1);  // to be set correctly in fillVertices
4267     if (includeTrackingParticles_) {
4268       tcand_simTrkIdx.push_back(tpIdx);
4269       tcand_simTrkShareFrac.push_back(sharedFraction);
4270       tcand_simTrkNChi2.push_back(tpChi2);
4271       tcand_bestSimTrkIdx.push_back(bestKeyCount.key >= 0 ? tpKeyToIndex.at(bestKeyCount.key) : -1);
4272       tcand_bestFromFirstHitSimTrkIdx.push_back(
4273           bestFirstHitKeyCount.key >= 0 ? tpKeyToIndex.at(bestFirstHitKeyCount.key) : -1);
4274     } else {
4275       tcand_isTrue.push_back(!tpIdx.empty());
4276     }
4277     tcand_bestSimTrkShareFrac.push_back(bestShareFrac);
4278     tcand_bestSimTrkShareFracSimDenom.push_back(bestShareFracSimDenom);
4279     tcand_bestSimTrkShareFracSimClusterDenom.push_back(bestShareFracSimClusterDenom);
4280     tcand_bestSimTrkNChi2.push_back(bestChi2);
4281     tcand_bestFromFirstHitSimTrkShareFrac.push_back(bestFirstHitShareFrac);
4282     tcand_bestFromFirstHitSimTrkShareFracSimDenom.push_back(bestFirstHitShareFracSimDenom);
4283     tcand_bestFromFirstHitSimTrkShareFracSimClusterDenom.push_back(bestFirstHitShareFracSimClusterDenom);
4284     tcand_bestFromFirstHitSimTrkNChi2.push_back(bestFirstHitChi2);
4285 
4286     LogTrace("TrackingNtuple") << "Track cand #" << iCand << " glob " << iglobCand << " with q=" << trk.charge()
4287                                << ", pT=" << trk.pt() << " GeV, eta: " << trk.eta() << ", phi: " << trk.phi()
4288                                << ", nValid=" << nValid << " seed#=" << aCand.seedRef().key()
4289                                << " simMatch=" << isSimMatched << " nSimHits=" << nSimHits
4290                                << " sharedFraction=" << (sharedFraction.empty() ? -1 : sharedFraction[0])
4291                                << " tpIdx=" << (tpIdx.empty() ? -1 : tpIdx[0]);
4292     std::vector<int> hitIdx;
4293     std::vector<int> hitType;
4294 
4295     for (auto i = aCand.recHits().begin(); i != aCand.recHits().end(); i++) {
4296       const TrackingRecHit* hit = &*i;
4297       DetId hitId = hit->geographicalId();
4298       LogTrace("TrackingNtuple") << "hit #" << std::distance(aCand.recHits().begin(), i)
4299                                  << " subdet=" << hitId.subdetId();
4300       if (hitId.det() != DetId::Tracker)
4301         continue;
4302 
4303       LogTrace("TrackingNtuple") << " " << subdetstring(hitId.subdetId()) << " " << tTopo.layer(hitId);
4304 
4305       if (hit->isValid()) {
4306         //ugly... but works
4307         const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*hit);
4308         const auto& clusterRef = bhit->firstClusterRef();
4309         unsigned int clusterKey;
4310         if (clusterRef.isPixel()) {
4311           clusterKey = clusterRef.cluster_pixel().key();
4312         } else if (clusterRef.isPhase2()) {
4313           clusterKey = clusterRef.cluster_phase2OT().key();
4314         } else {
4315           clusterKey = clusterRef.cluster_strip().key();
4316         }
4317 
4318         LogTrace("TrackingNtuple") << " id: " << hitId.rawId() << " - globalPos =" << hit->globalPosition()
4319                                    << " cluster=" << clusterKey << " clusterRef ID=" << clusterRef.id()
4320                                    << " eta,phi: " << hit->globalPosition().eta() << "," << hit->globalPosition().phi();
4321         if (includeAllHits_) {
4322           checkProductID(hitProductIds, clusterRef.id(), "track");
4323           if (clusterRef.isPixel()) {
4324             pix_tcandIdx[clusterKey].push_back(iglobCand);
4325           } else if (clusterRef.isPhase2()) {
4326             ph2_tcandIdx[clusterKey].push_back(iglobCand);
4327           } else {
4328             str_tcandIdx[clusterKey].push_back(iglobCand);
4329           }
4330         }
4331 
4332         hitIdx.push_back(clusterKey);
4333         if (clusterRef.isPixel()) {
4334           hitType.push_back(static_cast<int>(HitType::Pixel));
4335         } else if (clusterRef.isPhase2()) {
4336           hitType.push_back(static_cast<int>(HitType::Phase2OT));
4337         } else {
4338           hitType.push_back(static_cast<int>(HitType::Strip));
4339         }
4340       }
4341     }
4342 
4343     tcand_hitIdx.push_back(hitIdx);
4344     tcand_hitType.push_back(hitType);
4345   }
4346 }
4347 
4348 void TrackingNtuple::fillTrackingParticles(const edm::Event& iEvent,
4349                                            const edm::EventSetup& iSetup,
4350                                            const edm::RefToBaseVector<reco::Track>& tracks,
4351                                            const reco::BeamSpot& bs,
4352                                            const TrackingParticleRefVector& tpCollection,
4353                                            const TrackingVertexRefKeyToIndex& tvKeyToIndex,
4354                                            const reco::TrackToTrackingParticleAssociator& associatorByHits,
4355                                            const std::vector<TPHitIndex>& tpHitList,
4356                                            const TrackingParticleRefKeyToCount& tpKeyToClusterCount) {
4357   // Number of 3D layers for TPs
4358   edm::Handle<edm::ValueMap<unsigned int>> tpNLayersH;
4359   iEvent.getByToken(tpNLayersToken_, tpNLayersH);
4360   const auto& nLayers_tPCeff = *tpNLayersH;
4361 
4362   iEvent.getByToken(tpNPixelLayersToken_, tpNLayersH);
4363   const auto& nPixelLayers_tPCeff = *tpNLayersH;
4364 
4365   iEvent.getByToken(tpNStripStereoLayersToken_, tpNLayersH);
4366   const auto& nStripMonoAndStereoLayers_tPCeff = *tpNLayersH;
4367 
4368   reco::SimToRecoCollection simRecColl = associatorByHits.associateSimToReco(tracks, tpCollection);
4369 
4370   for (const TrackingParticleRef& tp : tpCollection) {
4371     LogTrace("TrackingNtuple") << "tracking particle pt=" << tp->pt() << " eta=" << tp->eta() << " phi=" << tp->phi();
4372     bool isRecoMatched = false;
4373     std::vector<int> tkIdx;
4374     std::vector<float> sharedFraction;
4375     auto foundTracks = simRecColl.find(tp);
4376     if (foundTracks != simRecColl.end()) {
4377       isRecoMatched = true;
4378       for (const auto& trackQuality : foundTracks->val) {
4379         sharedFraction.push_back(trackQuality.second);
4380         tkIdx.push_back(trackQuality.first.key());
4381       }
4382     }
4383 
4384     sim_genPdgIds.emplace_back();
4385     for (const auto& genRef : tp->genParticles()) {
4386       if (genRef.isNonnull())
4387         sim_genPdgIds.back().push_back(genRef->pdgId());
4388     }
4389 
4390     bool isFromBHadron = false;
4391     // Logic is similar to SimTracker/TrackHistory
4392     if (tracer_.evaluate(tp)) {  // ignore TP if history can not be traced
4393       // following is from TrackClassifier::processesAtGenerator()
4394       HistoryBase::RecoGenParticleTrail const& recoGenParticleTrail = tracer_.recoGenParticleTrail();
4395       for (const auto& particle : recoGenParticleTrail) {
4396         HepPDT::ParticleID particleID(particle->pdgId());
4397         if (particleID.hasBottom()) {
4398           isFromBHadron = true;
4399           break;
4400         }
4401       }
4402     }
4403 
4404     LogTrace("TrackingNtuple") << "matched to tracks = " << make_VectorPrinter(tkIdx)
4405                                << " isRecoMatched=" << isRecoMatched;
4406     sim_event.push_back(tp->eventId().event());
4407     sim_bunchCrossing.push_back(tp->eventId().bunchCrossing());
4408     sim_pdgId.push_back(tp->pdgId());
4409     sim_isFromBHadron.push_back(isFromBHadron);
4410     sim_px.push_back(tp->px());
4411     sim_py.push_back(tp->py());
4412     sim_pz.push_back(tp->pz());
4413     sim_pt.push_back(tp->pt());
4414     sim_eta.push_back(tp->eta());
4415     sim_phi.push_back(tp->phi());
4416     sim_q.push_back(tp->charge());
4417     sim_trkIdx.push_back(tkIdx);
4418     sim_trkShareFrac.push_back(sharedFraction);
4419     sim_parentVtxIdx.push_back(tvKeyToIndex.at(tp->parentVertex().key()));
4420     std::vector<int> decayIdx;
4421     for (const auto& v : tp->decayVertices())
4422       decayIdx.push_back(tvKeyToIndex.at(v.key()));
4423     sim_decayVtxIdx.push_back(decayIdx);
4424 
4425     //Calcualte the impact parameters w.r.t. PCA
4426     TrackingParticle::Vector momentum = parametersDefiner_.momentum(iEvent, iSetup, tp);
4427     TrackingParticle::Point vertex = parametersDefiner_.vertex(iEvent, iSetup, tp);
4428     auto dxySim = TrackingParticleIP::dxy(vertex, momentum, bs.position());
4429     auto dzSim = TrackingParticleIP::dz(vertex, momentum, bs.position());
4430     const double lambdaSim = M_PI / 2 - momentum.theta();
4431     sim_pca_pt.push_back(std::sqrt(momentum.perp2()));
4432     sim_pca_eta.push_back(momentum.Eta());
4433     sim_pca_lambda.push_back(lambdaSim);
4434     sim_pca_cotTheta.push_back(1 / tan(M_PI * 0.5 - lambdaSim));
4435     sim_pca_phi.push_back(momentum.phi());
4436     sim_pca_dxy.push_back(dxySim);
4437     sim_pca_dz.push_back(dzSim);
4438 
4439     std::vector<int> hitIdx;
4440     int nPixel = 0, nStrip = 0;
4441     auto rangeHit = std::equal_range(tpHitList.begin(), tpHitList.end(), TPHitIndex(tp.key()), tpHitIndexListLess);
4442     for (auto ip = rangeHit.first; ip != rangeHit.second; ++ip) {
4443       auto type = HitType::Unknown;
4444       if (!simhit_hitType[ip->simHitIdx].empty())
4445         type = static_cast<HitType>(simhit_hitType[ip->simHitIdx][0]);
4446       LogTrace("TrackingNtuple") << "simhit=" << ip->simHitIdx << " type=" << static_cast<int>(type);
4447       hitIdx.push_back(ip->simHitIdx);
4448       const auto detid = DetId(includeStripHits_ ? simhit_detId[ip->simHitIdx] : simhit_detId_phase2[ip->simHitIdx]);
4449       if (detid.det() != DetId::Tracker) {
4450         throw cms::Exception("LogicError") << "Encountered SimHit for TP " << tp.key() << " with DetId "
4451                                            << detid.rawId() << " whose det() is not Tracker but " << detid.det();
4452       }
4453       const auto subdet = detid.subdetId();
4454       switch (subdet) {
4455         case PixelSubdetector::PixelBarrel:
4456         case PixelSubdetector::PixelEndcap:
4457           ++nPixel;
4458           break;
4459         case StripSubdetector::TIB:
4460         case StripSubdetector::TID:
4461         case StripSubdetector::TOB:
4462         case StripSubdetector::TEC:
4463           ++nStrip;
4464           break;
4465         default:
4466           throw cms::Exception("LogicError") << "Encountered SimHit for TP " << tp.key() << " with DetId "
4467                                              << detid.rawId() << " whose subdet is not recognized, is " << subdet;
4468       };
4469     }
4470     sim_nValid.push_back(hitIdx.size());
4471     sim_nPixel.push_back(nPixel);
4472     sim_nStrip.push_back(nStrip);
4473 
4474     const auto nSimLayers = nLayers_tPCeff[tp];
4475     const auto nSimPixelLayers = nPixelLayers_tPCeff[tp];
4476     const auto nSimStripMonoAndStereoLayers = nStripMonoAndStereoLayers_tPCeff[tp];
4477     sim_nLay.push_back(nSimLayers);
4478     sim_nPixelLay.push_back(nSimPixelLayers);
4479     sim_n3DLay.push_back(nSimPixelLayers + nSimStripMonoAndStereoLayers);
4480 
4481     sim_nTrackerHits.push_back(tp->numberOfTrackerHits());
4482     auto found = tpKeyToClusterCount.find(tp.key());
4483     sim_nRecoClusters.push_back(found != cend(tpKeyToClusterCount) ? found->second : 0);
4484 
4485     sim_simHitIdx.push_back(hitIdx);
4486   }
4487 }
4488 
4489 // called from fillSeeds
4490 void TrackingNtuple::fillTrackingParticlesForSeeds(const TrackingParticleRefVector& tpCollection,
4491                                                    const reco::SimToRecoCollection& simRecColl,
4492                                                    const TrackingParticleRefKeyToIndex& tpKeyToIndex,
4493                                                    const unsigned int seedOffset) {
4494   if (sim_seedIdx.empty())  // first call
4495     sim_seedIdx.resize(tpCollection.size());
4496 
4497   for (const auto& keyVal : simRecColl) {
4498     const auto& tpRef = keyVal.key;
4499     auto found = tpKeyToIndex.find(tpRef.key());
4500     if (found == tpKeyToIndex.end())
4501       throw cms::Exception("Assert") << __FILE__ << ":" << __LINE__ << " fillTrackingParticlesForSeeds: tpRef.key() "
4502                                      << tpRef.key() << " not found from tpKeyToIndex. tpKeyToIndex size "
4503                                      << tpKeyToIndex.size();
4504     const auto tpIndex = found->second;
4505     for (const auto& pair : keyVal.val) {
4506       const auto& seedRef = pair.first->seedRef();
4507       sim_seedIdx[tpIndex].push_back(seedOffset + seedRef.key());
4508     }
4509   }
4510 }
4511 
4512 void TrackingNtuple::fillVertices(const reco::VertexCollection& vertices,
4513                                   const edm::RefToBaseVector<reco::Track>& tracks) {
4514   for (size_t iVertex = 0, size = vertices.size(); iVertex < size; ++iVertex) {
4515     const reco::Vertex& vertex = vertices[iVertex];
4516     vtx_x.push_back(vertex.x());
4517     vtx_y.push_back(vertex.y());
4518     vtx_z.push_back(vertex.z());
4519     vtx_xErr.push_back(vertex.xError());
4520     vtx_yErr.push_back(vertex.yError());
4521     vtx_zErr.push_back(vertex.zError());
4522     vtx_chi2.push_back(vertex.chi2());
4523     vtx_ndof.push_back(vertex.ndof());
4524     vtx_fake.push_back(vertex.isFake());
4525     vtx_valid.push_back(vertex.isValid());
4526 
4527     std::vector<int> trkIdx;
4528     for (auto iTrack = vertex.tracks_begin(); iTrack != vertex.tracks_end(); ++iTrack) {
4529       // Ignore link if vertex was fit from a track collection different from the input
4530       if (iTrack->id() != tracks.id())
4531         continue;
4532 
4533       trkIdx.push_back(iTrack->key());
4534 
4535       if (trk_vtxIdx[iTrack->key()] != -1) {
4536         throw cms::Exception("LogicError") << "Vertex index has already been set for track " << iTrack->key() << " to "
4537                                            << trk_vtxIdx[iTrack->key()] << "; was trying to set it to " << iVertex;
4538       }
4539       trk_vtxIdx[iTrack->key()] = iVertex;
4540     }
4541     vtx_trkIdx.push_back(trkIdx);
4542   }
4543 }
4544 
4545 void TrackingNtuple::fillTrackingVertices(const TrackingVertexRefVector& trackingVertices,
4546                                           const TrackingParticleRefKeyToIndex& tpKeyToIndex) {
4547   int current_event = -1;
4548   for (const auto& ref : trackingVertices) {
4549     const TrackingVertex v = *ref;
4550     if (v.eventId().event() != current_event) {
4551       // next PV
4552       current_event = v.eventId().event();
4553       simpv_idx.push_back(simvtx_x.size());
4554     }
4555 
4556     unsigned int processType = std::numeric_limits<unsigned int>::max();
4557     if (!v.g4Vertices().empty()) {
4558       processType = v.g4Vertices()[0].processType();
4559     }
4560 
4561     simvtx_event.push_back(v.eventId().event());
4562     simvtx_bunchCrossing.push_back(v.eventId().bunchCrossing());
4563     simvtx_processType.push_back(processType);
4564 
4565     simvtx_x.push_back(v.position().x());
4566     simvtx_y.push_back(v.position().y());
4567     simvtx_z.push_back(v.position().z());
4568 
4569     auto fill = [&](const TrackingParticleRefVector& tps, std::vector<int>& idx) {
4570       for (const auto& tpRef : tps) {
4571         auto found = tpKeyToIndex.find(tpRef.key());
4572         if (found != tpKeyToIndex.end()) {
4573           idx.push_back(found->second);
4574         }
4575       }
4576     };
4577 
4578     std::vector<int> sourceIdx;
4579     std::vector<int> daughterIdx;
4580     fill(v.sourceTracks(), sourceIdx);
4581     fill(v.daughterTracks(), daughterIdx);
4582 
4583     simvtx_sourceSimIdx.push_back(sourceIdx);
4584     simvtx_daughterSimIdx.push_back(daughterIdx);
4585   }
4586 }
4587 
4588 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
4589 void TrackingNtuple::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
4590   //The following says we do not know what parameters are allowed so do no validation
4591   // Please change this to state exactly what you do use, even if it is no parameters
4592   edm::ParameterSetDescription desc;
4593   desc.addUntracked<std::vector<edm::InputTag>>(
4594       "seedTracks",
4595       std::vector<edm::InputTag>{edm::InputTag("seedTracksinitialStepSeeds"),
4596                                  edm::InputTag("seedTracksdetachedTripletStepSeeds"),
4597                                  edm::InputTag("seedTrackspixelPairStepSeeds"),
4598                                  edm::InputTag("seedTrackslowPtTripletStepSeeds"),
4599                                  edm::InputTag("seedTracksmixedTripletStepSeeds"),
4600                                  edm::InputTag("seedTrackspixelLessStepSeeds"),
4601                                  edm::InputTag("seedTrackstobTecStepSeeds"),
4602                                  edm::InputTag("seedTracksjetCoreRegionalStepSeeds"),
4603                                  edm::InputTag("seedTracksmuonSeededSeedsInOut"),
4604                                  edm::InputTag("seedTracksmuonSeededSeedsOutIn")});
4605   desc.addUntracked<std::vector<edm::InputTag>>(
4606       "trackCandidates",
4607       std::vector<edm::InputTag>{edm::InputTag("initialStepTrackCandidates"),
4608                                  edm::InputTag("detachedTripletStepTrackCandidates"),
4609                                  edm::InputTag("pixelPairStepTrackCandidates"),
4610                                  edm::InputTag("lowPtTripletStepTrackCandidates"),
4611                                  edm::InputTag("mixedTripletStepTrackCandidates"),
4612                                  edm::InputTag("pixelLessStepTrackCandidates"),
4613                                  edm::InputTag("tobTecStepTrackCandidates"),
4614                                  edm::InputTag("jetCoreRegionalStepTrackCandidates"),
4615                                  edm::InputTag("muonSeededTrackCandidatesInOut"),
4616                                  edm::InputTag("muonSeededTrackCandidatesOutIn")});
4617   desc.addUntracked<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
4618   desc.addUntracked<std::vector<std::string>>("trackMVAs", std::vector<std::string>{{"generalTracks"}});
4619 
4620   edm::ParameterSetDescription cMaskDesc;
4621   cMaskDesc.addUntracked<unsigned int>("index");
4622   cMaskDesc.addUntracked<edm::InputTag>("src");
4623   std::vector<edm::ParameterSet> cMasks;
4624   auto addMask = [&cMasks](reco::Track::TrackAlgorithm algo) {
4625     edm::ParameterSet ps;
4626     ps.addUntrackedParameter<unsigned int>("index", static_cast<unsigned int>(algo));
4627     ps.addUntrackedParameter<edm::InputTag>("src", {reco::Track::algoName(algo) + "Clusters"});
4628     cMasks.push_back(ps);
4629   };
4630   addMask(reco::Track::detachedQuadStep);
4631   addMask(reco::Track::highPtTripletStep);
4632   addMask(reco::Track::detachedTripletStep);
4633   addMask(reco::Track::lowPtQuadStep);
4634   addMask(reco::Track::lowPtTripletStep);
4635   addMask(reco::Track::mixedTripletStep);
4636   addMask(reco::Track::pixelLessStep);
4637   addMask(reco::Track::pixelPairStep);
4638   addMask(reco::Track::tobTecStep);
4639   desc.addVPSetUntracked("clusterMasks", cMaskDesc, cMasks);
4640 
4641   desc.addUntracked<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
4642   desc.addUntracked<bool>("trackingParticlesRef", false);
4643   desc.addUntracked<edm::InputTag>("clusterTPMap", edm::InputTag("tpClusterProducer"));
4644   desc.addUntracked<edm::InputTag>("simHitTPMap", edm::InputTag("simHitTPAssocProducer"));
4645   desc.addUntracked<edm::InputTag>("trackAssociator", edm::InputTag("quickTrackAssociatorByHits"));
4646   desc.addUntracked<edm::InputTag>("pixelDigiSimLink", edm::InputTag("simSiPixelDigis"));
4647   desc.addUntracked<edm::InputTag>("stripDigiSimLink", edm::InputTag("simSiStripDigis"));
4648   desc.addUntracked<edm::InputTag>("phase2OTSimLink", edm::InputTag(""));
4649   desc.addUntracked<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
4650   desc.addUntracked<edm::InputTag>("pixelRecHits", edm::InputTag("siPixelRecHits"));
4651   desc.addUntracked<edm::InputTag>("stripRphiRecHits", edm::InputTag("siStripMatchedRecHits", "rphiRecHit"));
4652   desc.addUntracked<edm::InputTag>("stripStereoRecHits", edm::InputTag("siStripMatchedRecHits", "stereoRecHit"));
4653   desc.addUntracked<edm::InputTag>("stripMatchedRecHits", edm::InputTag("siStripMatchedRecHits", "matchedRecHit"));
4654   desc.addUntracked<edm::InputTag>("phase2OTRecHits", edm::InputTag("siPhase2RecHits"));
4655   desc.addUntracked<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
4656   desc.addUntracked<edm::InputTag>("trackingVertices", edm::InputTag("mix", "MergedTrackTruth"));
4657   desc.addUntracked<edm::InputTag>("trackingParticleNlayers",
4658                                    edm::InputTag("trackingParticleNumberOfLayersProducer", "trackerLayers"));
4659   desc.addUntracked<edm::InputTag>("trackingParticleNpixellayers",
4660                                    edm::InputTag("trackingParticleNumberOfLayersProducer", "pixelLayers"));
4661   desc.addUntracked<edm::InputTag>("trackingParticleNstripstereolayers",
4662                                    edm::InputTag("trackingParticleNumberOfLayersProducer", "stripStereoLayers"));
4663   desc.addUntracked<std::string>("TTRHBuilder", "WithTrackAngle")
4664       ->setComment("currently not used: keep for possible future use");
4665   desc.addUntracked<bool>("includeSeeds", false);
4666   desc.addUntracked<bool>("includeTrackCandidates", false);
4667   desc.addUntracked<bool>("addSeedCurvCov", false);
4668   desc.addUntracked<bool>("includeAllHits", false);
4669   desc.addUntracked<bool>("includeOnTrackHitData", false);
4670   desc.addUntracked<bool>("includeMVA", true);
4671   desc.addUntracked<bool>("includeTrackingParticles", true);
4672   desc.addUntracked<bool>("includeOOT", false);
4673   desc.addUntracked<bool>("keepEleSimHits", false);
4674   desc.addUntracked<bool>("saveSimHitsP3", false);
4675   desc.addUntracked<bool>("simHitBySignificance", false);
4676   descriptions.add("trackingNtuple", desc);
4677 }
4678 
4679 //define this as a plug-in
4680 DEFINE_FWK_MODULE(TrackingNtuple);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!