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