File indexing completed on 2024-06-22 02:24:08
0001
0002
0003
0004
0005 #include "CLHEP/Units/SystemOfUnits.h"
0006 #include "DataFormats/Common/interface/ValueMap.h"
0007 #include "DataFormats/Common/interface/View.h"
0008 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0009 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0010 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0011 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0012 #include "DataFormats/MuonReco/interface/Muon.h"
0013 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0014 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0015 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0016 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/global/EDProducer.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0028 #include "FWCore/Utilities/interface/isFinite.h"
0029 #include "FWCore/Utilities/interface/transform.h"
0030 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0031 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0032 #include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
0033 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0034 #include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
0035 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0036 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0037 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0038 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0039 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0040 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0041
0042 #include <unordered_map>
0043 #include <memory>
0044
0045 class SimPFProducer : public edm::global::EDProducer<> {
0046 public:
0047 SimPFProducer(const edm::ParameterSet&);
0048
0049 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0052
0053 private:
0054
0055 const double superClusterThreshold_, neutralEMThreshold_, neutralHADThreshold_;
0056 const bool useTiming_;
0057 const bool useTimingQuality_;
0058 const double timingQualityThreshold_;
0059
0060
0061 const edm::EDGetTokenT<edm::View<reco::PFRecTrack>> pfRecTracks_;
0062 const edm::EDGetTokenT<edm::View<reco::Track>> tracks_;
0063 const edm::EDGetTokenT<edm::View<reco::Track>> gsfTracks_;
0064 const edm::EDGetTokenT<reco::MuonCollection> muons_;
0065 const edm::EDGetTokenT<edm::ValueMap<float>> srcTrackTime_, srcTrackTimeError_, srcTrackTimeQuality_;
0066 const edm::EDGetTokenT<edm::ValueMap<float>> srcGsfTrackTime_, srcGsfTrackTimeError_, srcGsfTrackTimeQuality_;
0067 const edm::EDGetTokenT<TrackingParticleCollection> trackingParticles_;
0068 const edm::EDGetTokenT<SimClusterCollection> simClustersTruth_;
0069 const edm::EDGetTokenT<CaloParticleCollection> caloParticles_;
0070 const edm::EDGetTokenT<std::vector<reco::PFCluster>> simClusters_;
0071
0072 const std::vector<edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator>> associators_;
0073 };
0074
0075 #include "FWCore/Framework/interface/MakerMacros.h"
0076 DEFINE_FWK_MODULE(SimPFProducer);
0077
0078 void SimPFProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0079
0080 edm::ParameterSetDescription desc;
0081 desc.add<edm::InputTag>("simClustersSrc", {"particleFlowClusterHGCalFromSimCl"});
0082 desc.add<edm::InputTag>("trackSrc", {"generalTracks"});
0083 desc.add<std::vector<edm::InputTag>>("associators",
0084 {
0085 {"quickTrackAssociatorByHits"},
0086 });
0087 desc.add<edm::InputTag>("pfRecTrackSrc", {"hgcalTrackCollection", "TracksInHGCal"});
0088 desc.add<edm::InputTag>("trackingParticleSrc", {"mix", "MergedTrackTruth"});
0089 desc.add<double>("neutralEMThreshold", 0.25);
0090 desc.add<edm::InputTag>("caloParticlesSrc", {"mix", "MergedCaloTruth"});
0091 desc.add<double>("superClusterThreshold", 4.0);
0092 desc.add<edm::InputTag>("simClusterTruthSrc", {"mix", "MergedCaloTruth"});
0093 desc.add<edm::InputTag>("muonSrc", {"muons1stStep"});
0094 desc.add<double>("neutralHADThreshold", 0.25);
0095 desc.add<edm::InputTag>("gsfTrackSrc", {"electronGsfTracks"});
0096
0097
0098 desc.add<bool>("useTiming", false);
0099 desc.add<bool>("useTimingQuality", false);
0100 desc.add<edm::InputTag>("trackTimeValueMap", {});
0101 desc.add<edm::InputTag>("trackTimeErrorMap", {});
0102 desc.add<edm::InputTag>("trackTimeQualityMap", {});
0103 desc.add<double>("timingQualityThreshold", 0);
0104 desc.add<edm::InputTag>("gsfTrackTimeValueMap", {});
0105 desc.add<edm::InputTag>("gsfTrackTimeErrorMap", {});
0106 desc.add<edm::InputTag>("gsfTrackTimeQualityMap", {});
0107
0108 descriptions.addWithDefaultLabel(desc);
0109 }
0110
0111 namespace {
0112
0113 template <typename T>
0114 uint64_t hashSimInfo(const T& simTruth, size_t i = 0) {
0115 uint64_t evtid = simTruth.eventId().rawId();
0116 uint64_t trackid = simTruth.g4Tracks()[i].trackId();
0117 return ((evtid << 3) + 23401923) ^ trackid;
0118 };
0119 }
0120
0121 SimPFProducer::SimPFProducer(const edm::ParameterSet& conf)
0122 : superClusterThreshold_(conf.getParameter<double>("superClusterThreshold")),
0123 neutralEMThreshold_(conf.getParameter<double>("neutralEMThreshold")),
0124 neutralHADThreshold_(conf.getParameter<double>("neutralHADThreshold")),
0125 useTiming_(conf.getParameter<bool>("useTiming")),
0126 useTimingQuality_(conf.getParameter<bool>("useTimingQuality")),
0127 timingQualityThreshold_(useTimingQuality_ ? conf.getParameter<double>("timingQualityThreshold") : -99.),
0128 pfRecTracks_(consumes<edm::View<reco::PFRecTrack>>(conf.getParameter<edm::InputTag>("pfRecTrackSrc"))),
0129 tracks_(consumes<edm::View<reco::Track>>(conf.getParameter<edm::InputTag>("trackSrc"))),
0130 gsfTracks_(consumes<edm::View<reco::Track>>(conf.getParameter<edm::InputTag>("gsfTrackSrc"))),
0131 muons_(consumes<reco::MuonCollection>(conf.getParameter<edm::InputTag>("muonSrc"))),
0132 srcTrackTime_(useTiming_ ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeValueMap"))
0133 : edm::EDGetTokenT<edm::ValueMap<float>>()),
0134 srcTrackTimeError_(useTiming_
0135 ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeErrorMap"))
0136 : edm::EDGetTokenT<edm::ValueMap<float>>()),
0137 srcTrackTimeQuality_(useTimingQuality_
0138 ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeQualityMap"))
0139 : edm::EDGetTokenT<edm::ValueMap<float>>()),
0140 srcGsfTrackTime_(useTiming_
0141 ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("gsfTrackTimeValueMap"))
0142 : edm::EDGetTokenT<edm::ValueMap<float>>()),
0143 srcGsfTrackTimeError_(
0144 useTiming_ ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("gsfTrackTimeErrorMap"))
0145 : edm::EDGetTokenT<edm::ValueMap<float>>()),
0146 srcGsfTrackTimeQuality_(
0147 useTimingQuality_ ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("gsfTrackTimeQualityMap"))
0148 : edm::EDGetTokenT<edm::ValueMap<float>>()),
0149 trackingParticles_(consumes<TrackingParticleCollection>(conf.getParameter<edm::InputTag>("trackingParticleSrc"))),
0150 simClustersTruth_(consumes<SimClusterCollection>(conf.getParameter<edm::InputTag>("simClusterTruthSrc"))),
0151 caloParticles_(consumes<CaloParticleCollection>(conf.getParameter<edm::InputTag>("caloParticlesSrc"))),
0152 simClusters_(consumes<std::vector<reco::PFCluster>>(conf.getParameter<edm::InputTag>("simClustersSrc"))),
0153 associators_(edm::vector_transform(
0154 conf.getParameter<std::vector<edm::InputTag>>("associators"),
0155 [this](const edm::InputTag& tag) { return this->consumes<reco::TrackToTrackingParticleAssociator>(tag); })) {
0156 produces<reco::PFBlockCollection>();
0157 produces<reco::SuperClusterCollection>("perfect");
0158 produces<reco::PFCandidateCollection>();
0159 }
0160
0161 void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const {
0162
0163 std::vector<edm::Handle<reco::TrackToTrackingParticleAssociator>> associators;
0164 for (const auto& token : associators_) {
0165 associators.emplace_back();
0166 auto& back = associators.back();
0167 evt.getByToken(token, back);
0168 }
0169
0170
0171 edm::Handle<edm::View<reco::PFRecTrack>> PFTrackCollectionH;
0172 evt.getByToken(pfRecTracks_, PFTrackCollectionH);
0173 const edm::View<reco::PFRecTrack> PFTrackCollection = *PFTrackCollectionH;
0174 std::unordered_set<unsigned> PFTrackToGeneralTrack;
0175 for (unsigned i = 0; i < PFTrackCollection.size(); ++i) {
0176 const auto ptr = PFTrackCollection.ptrAt(i);
0177 PFTrackToGeneralTrack.insert(ptr->trackRef().key());
0178 }
0179
0180
0181 edm::Handle<edm::View<reco::Track>> TrackCollectionH;
0182 evt.getByToken(tracks_, TrackCollectionH);
0183 const edm::View<reco::Track>& TrackCollection = *TrackCollectionH;
0184
0185 edm::Handle<reco::MuonCollection> muons;
0186 evt.getByToken(muons_, muons);
0187 std::unordered_set<unsigned> MuonTrackToGeneralTrack;
0188 for (auto const& mu : *muons.product()) {
0189 reco::TrackRef muTrkRef = mu.track();
0190 if (muTrkRef.isNonnull())
0191 MuonTrackToGeneralTrack.insert(muTrkRef.key());
0192 }
0193
0194
0195 edm::Handle<edm::ValueMap<float>> trackTimeH, trackTimeErrH, trackTimeQualH, gsfTrackTimeH, gsfTrackTimeErrH,
0196 gsfTrackTimeQualH;
0197 if (useTiming_) {
0198 evt.getByToken(srcTrackTime_, trackTimeH);
0199 evt.getByToken(srcTrackTimeError_, trackTimeErrH);
0200 evt.getByToken(srcGsfTrackTime_, gsfTrackTimeH);
0201 evt.getByToken(srcGsfTrackTimeError_, gsfTrackTimeErrH);
0202
0203 if (useTimingQuality_) {
0204 evt.getByToken(srcTrackTimeQuality_, trackTimeQualH);
0205 evt.getByToken(srcGsfTrackTimeQuality_, gsfTrackTimeQualH);
0206 }
0207 }
0208
0209
0210 edm::Handle<TrackingParticleCollection> TPCollectionH;
0211 evt.getByToken(trackingParticles_, TPCollectionH);
0212
0213
0214
0215 edm::Handle<SimClusterCollection> SimClustersTruthH;
0216 evt.getByToken(simClustersTruth_, SimClustersTruthH);
0217 const SimClusterCollection& SimClustersTruth = *SimClustersTruthH;
0218
0219 edm::Handle<CaloParticleCollection> CaloParticlesH;
0220 evt.getByToken(caloParticles_, CaloParticlesH);
0221 const CaloParticleCollection& CaloParticles = *CaloParticlesH;
0222
0223 edm::Handle<std::vector<reco::PFCluster>> SimClustersH;
0224 evt.getByToken(simClusters_, SimClustersH);
0225 const std::vector<reco::PFCluster>& SimClusters = *SimClustersH;
0226
0227 std::unordered_map<uint64_t, size_t> hashToSimCluster;
0228
0229 for (unsigned i = 0; i < SimClustersTruth.size(); ++i) {
0230 const auto& simTruth = SimClustersTruth[i];
0231 hashToSimCluster[hashSimInfo(simTruth)] = i;
0232 }
0233
0234
0235 std::vector<reco::RecoToSimCollection> associatedTracks, associatedTracksGsf;
0236 for (const auto& associator : associators) {
0237 associatedTracks.emplace_back(associator->associateRecoToSim(TrackCollectionH, TPCollectionH));
0238
0239 }
0240
0241
0242
0243 auto superclusters = std::make_unique<reco::SuperClusterCollection>();
0244 auto blocks = std::make_unique<reco::PFBlockCollection>();
0245 std::unordered_map<size_t, size_t> simCluster2Block;
0246 std::unordered_map<size_t, size_t> simCluster2BlockIndex;
0247 std::unordered_multimap<size_t, size_t> caloParticle2SimCluster;
0248 std::vector<int> caloParticle2SuperCluster;
0249 for (unsigned icp = 0; icp < CaloParticles.size(); ++icp) {
0250 blocks->emplace_back();
0251 auto& block = blocks->back();
0252 const auto& simclusters = CaloParticles[icp].simClusters();
0253 double pttot = 0.0;
0254 double etot = 0.0;
0255 std::vector<size_t> good_simclusters;
0256 for (unsigned isc = 0; isc < simclusters.size(); ++isc) {
0257 auto simc = simclusters[isc];
0258 auto pdgId = std::abs(simc->pdgId());
0259 edm::Ref<std::vector<reco::PFCluster>> clusterRef(SimClustersH, simc.key());
0260 if (((pdgId == 22 || pdgId == 11) && clusterRef->energy() > neutralEMThreshold_) ||
0261 clusterRef->energy() > neutralHADThreshold_) {
0262 good_simclusters.push_back(isc);
0263 etot += clusterRef->energy();
0264 pttot += clusterRef->pt();
0265 auto bec = std::make_unique<reco::PFBlockElementCluster>(clusterRef, reco::PFBlockElement::HGCAL);
0266 block.addElement(bec.get());
0267 simCluster2Block[simc.key()] = icp;
0268 simCluster2BlockIndex[simc.key()] = bec->index();
0269 caloParticle2SimCluster.emplace(CaloParticles[icp].g4Tracks()[0].trackId(), simc.key());
0270 }
0271 }
0272
0273 auto pdgId = std::abs(CaloParticles[icp].pdgId());
0274
0275 caloParticle2SuperCluster.push_back(-1);
0276 if ((pdgId == 22 || pdgId == 11) && pttot > superClusterThreshold_) {
0277 caloParticle2SuperCluster[icp] = superclusters->size();
0278
0279 math::XYZPoint seedpos;
0280 reco::CaloClusterPtr seed;
0281 reco::CaloClusterPtrVector clusters;
0282 for (auto idx : good_simclusters) {
0283 edm::Ptr<reco::PFCluster> ptr(SimClustersH, simclusters[idx].key());
0284 clusters.push_back(ptr);
0285 if (seed.isNull() || seed->energy() < ptr->energy()) {
0286 seed = ptr;
0287 seedpos = ptr->position();
0288 }
0289 }
0290 superclusters->emplace_back(etot, seedpos, seed, clusters);
0291 }
0292 }
0293
0294 auto blocksHandle = evt.put(std::move(blocks));
0295 auto superClustersHandle = evt.put(std::move(superclusters), "perfect");
0296
0297
0298 std::vector<bool> usedTrack(TrackCollection.size(), false),
0299
0300 usedSimCluster(SimClusters.size(), false);
0301
0302 auto candidates = std::make_unique<reco::PFCandidateCollection>();
0303
0304 for (unsigned itk = 0; itk < TrackCollection.size(); ++itk) {
0305 auto tkRef = TrackCollection.refAt(itk);
0306
0307 if (PFTrackToGeneralTrack.count(itk) == 0)
0308 continue;
0309 reco::RecoToSimCollection::const_iterator assoc_tps = associatedTracks.back().end();
0310 for (const auto& association : associatedTracks) {
0311 assoc_tps = association.find(tkRef);
0312 if (assoc_tps != association.end())
0313 break;
0314 }
0315 if (assoc_tps == associatedTracks.back().end())
0316 continue;
0317
0318 const auto& matches = assoc_tps->val;
0319
0320 const auto absPdgId = std::abs(matches[0].first->pdgId());
0321 const auto charge = tkRef->charge();
0322 const auto three_mom = tkRef->momentum();
0323 constexpr double mpion2 = 0.13957 * 0.13957;
0324 double energy = std::sqrt(three_mom.mag2() + mpion2);
0325 math::XYZTLorentzVector trk_p4(three_mom.x(), three_mom.y(), three_mom.z(), energy);
0326
0327 reco::PFCandidate::ParticleType part_type;
0328
0329 switch (absPdgId) {
0330 case 11:
0331 part_type = reco::PFCandidate::e;
0332 break;
0333 case 13:
0334 part_type = reco::PFCandidate::mu;
0335 break;
0336 default:
0337 part_type = reco::PFCandidate::h;
0338 }
0339
0340 candidates->emplace_back(charge, trk_p4, part_type);
0341 auto& candidate = candidates->back();
0342
0343 candidate.setTrackRef(tkRef.castTo<reco::TrackRef>());
0344
0345 if (useTiming_) {
0346
0347 const bool assocQuality = useTimingQuality_ ? (*trackTimeQualH)[tkRef] > timingQualityThreshold_ : true;
0348 if (assocQuality) {
0349 candidate.setTime((*trackTimeH)[tkRef], (*trackTimeErrH)[tkRef]);
0350 } else {
0351 candidate.setTime(0., -1.);
0352 }
0353 }
0354
0355
0356 for (const auto& match : matches) {
0357 uint64_t hash = hashSimInfo(*(match.first));
0358 if (hashToSimCluster.count(hash)) {
0359 auto simcHash = hashToSimCluster[hash];
0360
0361 if (!usedSimCluster[simcHash]) {
0362 if (simCluster2Block.count(simcHash) && simCluster2BlockIndex.count(simcHash)) {
0363 size_t block = simCluster2Block.find(simcHash)->second;
0364 size_t blockIdx = simCluster2BlockIndex.find(simcHash)->second;
0365 edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle, block);
0366 candidate.addElementInBlock(blockRef, blockIdx);
0367 usedSimCluster[simcHash] = true;
0368 }
0369 }
0370 if (absPdgId == 11) {
0371 if (simCluster2Block.count(simcHash)) {
0372 auto block_index = simCluster2Block.find(simcHash)->second;
0373 auto supercluster_index = caloParticle2SuperCluster[block_index];
0374 if (supercluster_index != -1) {
0375 edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle, block_index);
0376 for (const auto& elem : blockRef->elements()) {
0377 const auto& ref = elem.clusterRef();
0378 if (!usedSimCluster[ref.key()]) {
0379 candidate.addElementInBlock(blockRef, elem.index());
0380 usedSimCluster[ref.key()] = true;
0381 }
0382 }
0383
0384
0385 }
0386 }
0387 }
0388 }
0389
0390
0391
0392
0393
0394 if (caloParticle2SimCluster.count(match.first->g4Tracks()[0].trackId())) {
0395 auto range = caloParticle2SimCluster.equal_range(match.first->g4Tracks()[0].trackId());
0396 for (auto it = range.first; it != range.second; ++it) {
0397 if (!usedSimCluster[it->second]) {
0398 usedSimCluster[it->second] = true;
0399 if (simCluster2Block.find(it->second) != simCluster2Block.end()) {
0400 size_t block = simCluster2Block.find(it->second)->second;
0401 size_t blockIdx = simCluster2BlockIndex.find(it->second)->second;
0402 edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle, block);
0403 candidate.addElementInBlock(blockRef, blockIdx);
0404 }
0405 }
0406 }
0407 }
0408 }
0409 usedTrack[tkRef.key()] = true;
0410
0411 if (MuonTrackToGeneralTrack.count(itk) || absPdgId == 13)
0412 candidates->pop_back();
0413 }
0414
0415
0416
0417 const auto& theblocks = *blocksHandle;
0418 for (unsigned ibl = 0; ibl < theblocks.size(); ++ibl) {
0419 reco::PFBlockRef blref(blocksHandle, ibl);
0420 const auto& elements = theblocks[ibl].elements();
0421 for (const auto& elem : elements) {
0422 const auto& ref = elem.clusterRef();
0423 const auto& simtruth = SimClustersTruth[ref.key()];
0424 reco::PFCandidate::ParticleType part_type;
0425 if (!usedSimCluster[ref.key()]) {
0426 auto absPdgId = std::abs(simtruth.pdgId());
0427 switch (absPdgId) {
0428 case 11:
0429 case 22:
0430 part_type = reco::PFCandidate::gamma;
0431 break;
0432 default:
0433 part_type = reco::PFCandidate::h0;
0434 }
0435 const auto three_mom = (ref->position() - math::XYZPoint(0, 0, 0)).unit() * ref->correctedEnergy();
0436 math::XYZTLorentzVector clu_p4(three_mom.x(), three_mom.y(), three_mom.z(), ref->correctedEnergy());
0437 candidates->emplace_back(0, clu_p4, part_type);
0438 auto& candidate = candidates->back();
0439 candidate.addElementInBlock(blref, elem.index());
0440 }
0441 }
0442 }
0443
0444 evt.put(std::move(candidates));
0445 }