PFSimProducer/plugins/SimPFProducer.cc

0001 // This producer assigns vertex times (with a specified resolution) to tracks.
0002 // The times are produced as valuemaps associated to tracks, so the track dataformat doesn't
0003 // need to be modified.
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   // parameters
0055   const double superClusterThreshold_, neutralEMThreshold_, neutralHADThreshold_;
0056   const bool useTiming_;
0057   const bool useTimingQuality_;
0058   const double timingQualityThreshold_;
0059
0060   // inputs
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   // tracking particle associators by order of preference
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   // simPFProducer
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   // if useTiming_
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 }  // namespace
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   //get associators
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   //get PFRecTrack
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   //get track collections
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   // get timing, if enabled
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   //get tracking particle collections
0210   edm::Handle<TrackingParticleCollection> TPCollectionH;
0211   evt.getByToken(trackingParticles_, TPCollectionH);
0212   //const TrackingParticleCollection&  TPCollection = *TPCollectionH;
0213
0214   // grab phony clustering information
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   // associate the reco tracks / gsf Tracks
0235   std::vector<reco::RecoToSimCollection> associatedTracks, associatedTracksGsf;
0236   for (const auto& associator : associators) {
0237     associatedTracks.emplace_back(associator->associateRecoToSim(TrackCollectionH, TPCollectionH));
0238     //associatedTracksGsf.emplace_back(associator->associateRecoToSim(GsfTrackCollectionH, TPCollectionH));
0239   }
0240
0241   // make blocks out of calo particles so we can have cluster references
0242   // likewise fill out superclusters
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;  // take seed pos as supercluster point
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   // list tracks so we can mark them as used and/or fight over them
0298   std::vector<bool> usedTrack(TrackCollection.size(), false),
0299       //usedGsfTrack(GsfTrackCollection.size(),false),
0300       usedSimCluster(SimClusters.size(), false);
0301
0302   auto candidates = std::make_unique<reco::PFCandidateCollection>();
0303   // in good particle flow fashion, start from the tracks and go out
0304   for (unsigned itk = 0; itk < TrackCollection.size(); ++itk) {
0305     auto tkRef = TrackCollection.refAt(itk);
0306     // skip tracks not selected by PF
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     // assured now that we are matched to a set of tracks
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       // check if track-mtd match is of sufficient quality
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     // bind to cluster if there is one and try to gather conversions, etc
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) {  // collect brems/conv. brems
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               //*TODO* cluster time is not reliable at the moment, so just keep time from the track if available
0385             }
0386           }
0387         }
0388       }
0389       // Now try to include also electrons that have been reconstructed using
0390       // the GraphCaloParticles. In particular, recover the cases in which the
0391       // tracking particle associated to the CaloParticle has not left any hits
0392       // in the calorimeters or, if it had, the cluster has been skipped due to
0393       // threshold requirements.
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     // remove tracks already used by muons
0411     if (MuonTrackToGeneralTrack.count(itk) || absPdgId == 13)
0412       candidates->pop_back();
0413   }
0414
0415   // now loop over the non-collected clusters in blocks
0416   // and turn them into neutral hadrons or photons
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 }