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.addOptional<edm::InputTag>("trackTimeValueMap");
0099   desc.addOptional<edm::InputTag>("trackTimeErrorMap");
0100   desc.addOptional<edm::InputTag>("trackTimeQualityMap");
0101   desc.addOptional<double>("timingQualityThreshold");
0102   desc.addOptional<edm::InputTag>("gsfTrackTimeValueMap");
0103   desc.addOptional<edm::InputTag>("gsfTrackTimeErrorMap");
0104   desc.addOptional<edm::InputTag>("gsfTrackTimeQualityMap");
0105
0106   descriptions.add("simPFProducer", desc);
0107 }
0108
0109 namespace {
0110
0111   template <typename T>
0112   uint64_t hashSimInfo(const T& simTruth, size_t i = 0) {
0113     uint64_t evtid = simTruth.eventId().rawId();
0114     uint64_t trackid = simTruth.g4Tracks()[i].trackId();
0115     return ((evtid << 3) + 23401923) ^ trackid;
0116   };
0117 }  // namespace
0118
0119 SimPFProducer::SimPFProducer(const edm::ParameterSet& conf)
0120     : superClusterThreshold_(conf.getParameter<double>("superClusterThreshold")),
0121       neutralEMThreshold_(conf.getParameter<double>("neutralEMThreshold")),
0122       neutralHADThreshold_(conf.getParameter<double>("neutralHADThreshold")),
0123       useTiming_(conf.existsAs<edm::InputTag>("trackTimeValueMap")),
0124       useTimingQuality_(conf.existsAs<edm::InputTag>("trackTimeQualityMap")),
0125       timingQualityThreshold_(useTimingQuality_ ? conf.getParameter<double>("timingQualityThreshold") : -99.),
0126       pfRecTracks_(consumes<edm::View<reco::PFRecTrack>>(conf.getParameter<edm::InputTag>("pfRecTrackSrc"))),
0127       tracks_(consumes<edm::View<reco::Track>>(conf.getParameter<edm::InputTag>("trackSrc"))),
0128       gsfTracks_(consumes<edm::View<reco::Track>>(conf.getParameter<edm::InputTag>("gsfTrackSrc"))),
0129       muons_(consumes<reco::MuonCollection>(conf.getParameter<edm::InputTag>("muonSrc"))),
0130       srcTrackTime_(useTiming_ ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeValueMap"))
0131                                : edm::EDGetTokenT<edm::ValueMap<float>>()),
0132       srcTrackTimeError_(useTiming_
0133                              ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeErrorMap"))
0134                              : edm::EDGetTokenT<edm::ValueMap<float>>()),
0135       srcTrackTimeQuality_(useTimingQuality_
0136                                ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeQualityMap"))
0137                                : edm::EDGetTokenT<edm::ValueMap<float>>()),
0138       srcGsfTrackTime_(useTiming_
0139                            ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("gsfTrackTimeValueMap"))
0140                            : edm::EDGetTokenT<edm::ValueMap<float>>()),
0141       srcGsfTrackTimeError_(
0142           useTiming_ ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("gsfTrackTimeErrorMap"))
0143                      : edm::EDGetTokenT<edm::ValueMap<float>>()),
0144       srcGsfTrackTimeQuality_(
0145           useTimingQuality_ ? consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("gsfTrackTimeQualityMap"))
0146                             : edm::EDGetTokenT<edm::ValueMap<float>>()),
0147       trackingParticles_(consumes<TrackingParticleCollection>(conf.getParameter<edm::InputTag>("trackingParticleSrc"))),
0148       simClustersTruth_(consumes<SimClusterCollection>(conf.getParameter<edm::InputTag>("simClusterTruthSrc"))),
0149       caloParticles_(consumes<CaloParticleCollection>(conf.getParameter<edm::InputTag>("caloParticlesSrc"))),
0150       simClusters_(consumes<std::vector<reco::PFCluster>>(conf.getParameter<edm::InputTag>("simClustersSrc"))),
0151       associators_(edm::vector_transform(
0152           conf.getParameter<std::vector<edm::InputTag>>("associators"),
0153           [this](const edm::InputTag& tag) { return this->consumes<reco::TrackToTrackingParticleAssociator>(tag); })) {
0154   produces<reco::PFBlockCollection>();
0155   produces<reco::SuperClusterCollection>("perfect");
0156   produces<reco::PFCandidateCollection>();
0157 }
0158
0159 void SimPFProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const {
0160   //get associators
0161   std::vector<edm::Handle<reco::TrackToTrackingParticleAssociator>> associators;
0162   for (const auto& token : associators_) {
0163     associators.emplace_back();
0164     auto& back = associators.back();
0165     evt.getByToken(token, back);
0166   }
0167
0168   //get PFRecTrack
0169   edm::Handle<edm::View<reco::PFRecTrack>> PFTrackCollectionH;
0170   evt.getByToken(pfRecTracks_, PFTrackCollectionH);
0171   const edm::View<reco::PFRecTrack> PFTrackCollection = *PFTrackCollectionH;
0172   std::unordered_set<unsigned> PFTrackToGeneralTrack;
0173   for (unsigned i = 0; i < PFTrackCollection.size(); ++i) {
0174     const auto ptr = PFTrackCollection.ptrAt(i);
0175     PFTrackToGeneralTrack.insert(ptr->trackRef().key());
0176   }
0177
0178   //get track collections
0179   edm::Handle<edm::View<reco::Track>> TrackCollectionH;
0180   evt.getByToken(tracks_, TrackCollectionH);
0181   const edm::View<reco::Track>& TrackCollection = *TrackCollectionH;
0182
0183   edm::Handle<reco::MuonCollection> muons;
0184   evt.getByToken(muons_, muons);
0185   std::unordered_set<unsigned> MuonTrackToGeneralTrack;
0186   for (auto const& mu : *muons.product()) {
0187     reco::TrackRef muTrkRef = mu.track();
0188     if (muTrkRef.isNonnull())
0189       MuonTrackToGeneralTrack.insert(muTrkRef.key());
0190   }
0191
0192   // get timing, if enabled
0193   edm::Handle<edm::ValueMap<float>> trackTimeH, trackTimeErrH, trackTimeQualH, gsfTrackTimeH, gsfTrackTimeErrH,
0194       gsfTrackTimeQualH;
0195   if (useTiming_) {
0196     evt.getByToken(srcTrackTime_, trackTimeH);
0197     evt.getByToken(srcTrackTimeError_, trackTimeErrH);
0198     evt.getByToken(srcGsfTrackTime_, gsfTrackTimeH);
0199     evt.getByToken(srcGsfTrackTimeError_, gsfTrackTimeErrH);
0200
0201     if (useTimingQuality_) {
0202       evt.getByToken(srcTrackTimeQuality_, trackTimeQualH);
0203       evt.getByToken(srcGsfTrackTimeQuality_, gsfTrackTimeQualH);
0204     }
0205   }
0206
0207   //get tracking particle collections
0208   edm::Handle<TrackingParticleCollection> TPCollectionH;
0209   evt.getByToken(trackingParticles_, TPCollectionH);
0210   //const TrackingParticleCollection&  TPCollection = *TPCollectionH;
0211
0212   // grab phony clustering information
0213   edm::Handle<SimClusterCollection> SimClustersTruthH;
0214   evt.getByToken(simClustersTruth_, SimClustersTruthH);
0215   const SimClusterCollection& SimClustersTruth = *SimClustersTruthH;
0216
0217   edm::Handle<CaloParticleCollection> CaloParticlesH;
0218   evt.getByToken(caloParticles_, CaloParticlesH);
0219   const CaloParticleCollection& CaloParticles = *CaloParticlesH;
0220
0221   edm::Handle<std::vector<reco::PFCluster>> SimClustersH;
0222   evt.getByToken(simClusters_, SimClustersH);
0223   const std::vector<reco::PFCluster>& SimClusters = *SimClustersH;
0224
0225   std::unordered_map<uint64_t, size_t> hashToSimCluster;
0226
0227   for (unsigned i = 0; i < SimClustersTruth.size(); ++i) {
0228     const auto& simTruth = SimClustersTruth[i];
0229     hashToSimCluster[hashSimInfo(simTruth)] = i;
0230   }
0231
0232   // associate the reco tracks / gsf Tracks
0233   std::vector<reco::RecoToSimCollection> associatedTracks, associatedTracksGsf;
0234   for (const auto& associator : associators) {
0235     associatedTracks.emplace_back(associator->associateRecoToSim(TrackCollectionH, TPCollectionH));
0236     //associatedTracksGsf.emplace_back(associator->associateRecoToSim(GsfTrackCollectionH, TPCollectionH));
0237   }
0238
0239   // make blocks out of calo particles so we can have cluster references
0240   // likewise fill out superclusters
0241   auto superclusters = std::make_unique<reco::SuperClusterCollection>();
0242   auto blocks = std::make_unique<reco::PFBlockCollection>();
0243   std::unordered_map<size_t, size_t> simCluster2Block;
0244   std::unordered_map<size_t, size_t> simCluster2BlockIndex;
0245   std::unordered_multimap<size_t, size_t> caloParticle2SimCluster;
0246   std::vector<int> caloParticle2SuperCluster;
0247   for (unsigned icp = 0; icp < CaloParticles.size(); ++icp) {
0248     blocks->emplace_back();
0249     auto& block = blocks->back();
0250     const auto& simclusters = CaloParticles[icp].simClusters();
0251     double pttot = 0.0;
0252     double etot = 0.0;
0253     std::vector<size_t> good_simclusters;
0254     for (unsigned isc = 0; isc < simclusters.size(); ++isc) {
0255       auto simc = simclusters[isc];
0256       auto pdgId = std::abs(simc->pdgId());
0257       edm::Ref<std::vector<reco::PFCluster>> clusterRef(SimClustersH, simc.key());
0258       if (((pdgId == 22 || pdgId == 11) && clusterRef->energy() > neutralEMThreshold_) ||
0259           clusterRef->energy() > neutralHADThreshold_) {
0260         good_simclusters.push_back(isc);
0261         etot += clusterRef->energy();
0262         pttot += clusterRef->pt();
0263         auto bec = std::make_unique<reco::PFBlockElementCluster>(clusterRef, reco::PFBlockElement::HGCAL);
0264         block.addElement(bec.get());
0265         simCluster2Block[simc.key()] = icp;
0266         simCluster2BlockIndex[simc.key()] = bec->index();
0267         caloParticle2SimCluster.emplace(CaloParticles[icp].g4Tracks()[0].trackId(), simc.key());
0268       }
0269     }
0270
0271     auto pdgId = std::abs(CaloParticles[icp].pdgId());
0272
0273     caloParticle2SuperCluster.push_back(-1);
0274     if ((pdgId == 22 || pdgId == 11) && pttot > superClusterThreshold_) {
0275       caloParticle2SuperCluster[icp] = superclusters->size();
0276
0277       math::XYZPoint seedpos;  // take seed pos as supercluster point
0278       reco::CaloClusterPtr seed;
0279       reco::CaloClusterPtrVector clusters;
0280       for (auto idx : good_simclusters) {
0281         edm::Ptr<reco::PFCluster> ptr(SimClustersH, simclusters[idx].key());
0282         clusters.push_back(ptr);
0283         if (seed.isNull() || seed->energy() < ptr->energy()) {
0284           seed = ptr;
0285           seedpos = ptr->position();
0286         }
0287       }
0288       superclusters->emplace_back(etot, seedpos, seed, clusters);
0289     }
0290   }
0291
0292   auto blocksHandle = evt.put(std::move(blocks));
0293   auto superClustersHandle = evt.put(std::move(superclusters), "perfect");
0294
0295   // list tracks so we can mark them as used and/or fight over them
0296   std::vector<bool> usedTrack(TrackCollection.size(), false),
0297       //usedGsfTrack(GsfTrackCollection.size(),false),
0298       usedSimCluster(SimClusters.size(), false);
0299
0300   auto candidates = std::make_unique<reco::PFCandidateCollection>();
0301   // in good particle flow fashion, start from the tracks and go out
0302   for (unsigned itk = 0; itk < TrackCollection.size(); ++itk) {
0303     auto tkRef = TrackCollection.refAt(itk);
0304     // skip tracks not selected by PF
0305     if (PFTrackToGeneralTrack.count(itk) == 0)
0306       continue;
0307     reco::RecoToSimCollection::const_iterator assoc_tps = associatedTracks.back().end();
0308     for (const auto& association : associatedTracks) {
0309       assoc_tps = association.find(tkRef);
0310       if (assoc_tps != association.end())
0311         break;
0312     }
0313     if (assoc_tps == associatedTracks.back().end())
0314       continue;
0315     // assured now that we are matched to a set of tracks
0316     const auto& matches = assoc_tps->val;
0317
0318     const auto absPdgId = std::abs(matches[0].first->pdgId());
0319     const auto charge = tkRef->charge();
0320     const auto three_mom = tkRef->momentum();
0321     constexpr double mpion2 = 0.13957 * 0.13957;
0322     double energy = std::sqrt(three_mom.mag2() + mpion2);
0323     math::XYZTLorentzVector trk_p4(three_mom.x(), three_mom.y(), three_mom.z(), energy);
0324
0325     reco::PFCandidate::ParticleType part_type;
0326
0327     switch (absPdgId) {
0328       case 11:
0329         part_type = reco::PFCandidate::e;
0330         break;
0331       case 13:
0332         part_type = reco::PFCandidate::mu;
0333         break;
0334       default:
0335         part_type = reco::PFCandidate::h;
0336     }
0337
0338     candidates->emplace_back(charge, trk_p4, part_type);
0339     auto& candidate = candidates->back();
0340
0341     candidate.setTrackRef(tkRef.castTo<reco::TrackRef>());
0342
0343     if (useTiming_) {
0344       // check if track-mtd match is of sufficient quality
0345       const bool assocQuality = useTimingQuality_ ? (*trackTimeQualH)[tkRef] > timingQualityThreshold_ : true;
0346       if (assocQuality) {
0347         candidate.setTime((*trackTimeH)[tkRef], (*trackTimeErrH)[tkRef]);
0348       } else {
0349         candidate.setTime(0., -1.);
0350       }
0351     }
0352
0353     // bind to cluster if there is one and try to gather conversions, etc
0354     for (const auto& match : matches) {
0355       uint64_t hash = hashSimInfo(*(match.first));
0356       if (hashToSimCluster.count(hash)) {
0357         auto simcHash = hashToSimCluster[hash];
0358
0359         if (!usedSimCluster[simcHash]) {
0360           if (simCluster2Block.count(simcHash) && simCluster2BlockIndex.count(simcHash)) {
0361             size_t block = simCluster2Block.find(simcHash)->second;
0362             size_t blockIdx = simCluster2BlockIndex.find(simcHash)->second;
0363             edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle, block);
0364             candidate.addElementInBlock(blockRef, blockIdx);
0365             usedSimCluster[simcHash] = true;
0366           }
0367         }
0368         if (absPdgId == 11) {  // collect brems/conv. brems
0369           if (simCluster2Block.count(simcHash)) {
0370             auto block_index = simCluster2Block.find(simcHash)->second;
0371             auto supercluster_index = caloParticle2SuperCluster[block_index];
0372             if (supercluster_index != -1) {
0373               edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle, block_index);
0374               for (const auto& elem : blockRef->elements()) {
0375                 const auto& ref = elem.clusterRef();
0376                 if (!usedSimCluster[ref.key()]) {
0377                   candidate.addElementInBlock(blockRef, elem.index());
0378                   usedSimCluster[ref.key()] = true;
0379                 }
0380               }
0381
0382               //*TODO* cluster time is not reliable at the moment, so just keep time from the track if available
0383             }
0384           }
0385         }
0386       }
0387       // Now try to include also electrons that have been reconstructed using
0388       // the GraphCaloParticles. In particular, recover the cases in which the
0389       // tracking particle associated to the CaloParticle has not left any hits
0390       // in the calorimeters or, if it had, the cluster has been skipped due to
0391       // threshold requirements.
0392       if (caloParticle2SimCluster.count(match.first->g4Tracks()[0].trackId())) {
0393         auto range = caloParticle2SimCluster.equal_range(match.first->g4Tracks()[0].trackId());
0394         for (auto it = range.first; it != range.second; ++it) {
0395           if (!usedSimCluster[it->second]) {
0396             usedSimCluster[it->second] = true;
0397             if (simCluster2Block.find(it->second) != simCluster2Block.end()) {
0398               size_t block = simCluster2Block.find(it->second)->second;
0399               size_t blockIdx = simCluster2BlockIndex.find(it->second)->second;
0400               edm::Ref<reco::PFBlockCollection> blockRef(blocksHandle, block);
0401               candidate.addElementInBlock(blockRef, blockIdx);
0402             }
0403           }
0404         }
0405       }
0406     }
0407     usedTrack[tkRef.key()] = true;
0408     // remove tracks already used by muons
0409     if (MuonTrackToGeneralTrack.count(itk) || absPdgId == 13)
0410       candidates->pop_back();
0411   }
0412
0413   // now loop over the non-collected clusters in blocks
0414   // and turn them into neutral hadrons or photons
0415   const auto& theblocks = *blocksHandle;
0416   for (unsigned ibl = 0; ibl < theblocks.size(); ++ibl) {
0417     reco::PFBlockRef blref(blocksHandle, ibl);
0418     const auto& elements = theblocks[ibl].elements();
0419     for (const auto& elem : elements) {
0420       const auto& ref = elem.clusterRef();
0421       const auto& simtruth = SimClustersTruth[ref.key()];
0422       reco::PFCandidate::ParticleType part_type;
0423       if (!usedSimCluster[ref.key()]) {
0424         auto absPdgId = std::abs(simtruth.pdgId());
0425         switch (absPdgId) {
0426           case 11:
0427           case 22:
0428             part_type = reco::PFCandidate::gamma;
0429             break;
0430           default:
0431             part_type = reco::PFCandidate::h0;
0432         }
0433         const auto three_mom = (ref->position() - math::XYZPoint(0, 0, 0)).unit() * ref->correctedEnergy();
0434         math::XYZTLorentzVector clu_p4(three_mom.x(), three_mom.y(), three_mom.z(), ref->correctedEnergy());
0435         candidates->emplace_back(0, clu_p4, part_type);
0436         auto& candidate = candidates->back();
0437         candidate.addElementInBlock(blref, elem.index());
0438       }
0439     }
0440   }
0441
0442   evt.put(std::move(candidates));
0443 }