Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-21 22:40:04

0001 #include <string>
0002 
0003 #include "DataFormats/Candidate/interface/Candidate.h"
0004 #include "DataFormats/Common/interface/Association.h"
0005 #include "DataFormats/Common/interface/ValueMap.h"
0006 #include "DataFormats/Common/interface/View.h"
0007 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0008 #include "DataFormats/MuonReco/interface/Muon.h"
0009 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0010 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0011 #include "DataFormats/PatCandidates/interface/HcalDepthEnergyFractions.h"
0012 #include "DataFormats/PatCandidates/interface/Jet.h"
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/global/EDProducer.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024 
0025 /*#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0026 #include "MagneticField/Engine/interface/MagneticField.h"
0027 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0028 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0029 #include
0030 "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0032 */
0033 //#define CRAZYSORT
0034 
0035 namespace pat {
0036   /// conversion map from quality flags used in PV association and miniAOD one
0037   const static int qualityMap[8] = {1, 0, 1, 1, 4, 4, 5, 6};
0038 
0039   class PATPackedCandidateProducer : public edm::global::EDProducer<> {
0040   public:
0041     explicit PATPackedCandidateProducer(const edm::ParameterSet &);
0042     ~PATPackedCandidateProducer() override;
0043 
0044     void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0045 
0046     // sorting of cands to maximize the zlib compression
0047     static bool candsOrdering(pat::PackedCandidate const &i, pat::PackedCandidate const &j) {
0048       if (std::abs(i.charge()) == std::abs(j.charge())) {
0049         if (i.charge() != 0) {
0050           if (i.hasTrackDetails() and !j.hasTrackDetails())
0051             return true;
0052           if (!i.hasTrackDetails() and j.hasTrackDetails())
0053             return false;
0054           if (i.covarianceSchema() > j.covarianceSchema())
0055             return true;
0056           if (i.covarianceSchema() < j.covarianceSchema())
0057             return false;
0058         }
0059         if (i.vertexRef() == j.vertexRef())
0060           return i.eta() > j.eta();
0061         else
0062           return i.vertexRef().key() < j.vertexRef().key();
0063       }
0064       return std::abs(i.charge()) > std::abs(j.charge());
0065     }
0066 
0067     template <typename T>
0068     static std::vector<size_t> sort_indexes(const std::vector<T> &v) {
0069       std::vector<size_t> idx(v.size());
0070       for (size_t i = 0; i != idx.size(); ++i)
0071         idx[i] = i;
0072       std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return candsOrdering(v[i1], v[i2]); });
0073       return idx;
0074     }
0075 
0076   private:
0077     // if PuppiSrc && PuppiNoLepSrc are empty, usePuppi is false
0078     // otherwise assumes that if they are set, you wanted to use puppi and will
0079     // throw an exception if the puppis are not found
0080     const bool usePuppi_;
0081 
0082     const edm::EDGetTokenT<reco::PFCandidateCollection> Cands_;
0083     const edm::EDGetTokenT<reco::VertexCollection> PVs_;
0084     const edm::EDGetTokenT<edm::Association<reco::VertexCollection>> PVAsso_;
0085     const edm::EDGetTokenT<edm::ValueMap<int>> PVAssoQuality_;
0086     const edm::EDGetTokenT<reco::VertexCollection> PVOrigs_;
0087     const edm::EDGetTokenT<reco::TrackCollection> TKOrigs_;
0088     const edm::EDGetTokenT<edm::ValueMap<float>> PuppiWeight_;
0089     const edm::EDGetTokenT<edm::ValueMap<float>> PuppiWeightNoLep_;
0090     std::vector<edm::EDGetTokenT<edm::View<reco::Candidate>>> SVWhiteLists_;
0091     const bool storeChargedHadronIsolation_;
0092     const edm::EDGetTokenT<edm::ValueMap<bool>> ChargedHadronIsolation_;
0093 
0094     const double minPtForChargedHadronProperties_;
0095     const double minPtForTrackProperties_;
0096     const double minPtForLowQualityTrackProperties_;
0097     const int covarianceVersion_;
0098     const std::vector<int> covariancePackingSchemas_;
0099 
0100     const std::vector<int> pfCandidateTypesForHcalDepth_;
0101     const bool storeHcalDepthEndcapOnly_;
0102 
0103     const bool storeTiming_;
0104     const bool timeFromValueMap_;
0105     const edm::EDGetTokenT<edm::ValueMap<float>> t0Map_;
0106     const edm::EDGetTokenT<edm::ValueMap<float>> t0ErrMap_;
0107 
0108     // for debugging
0109     float calcDxy(float dx, float dy, float phi) const { return -dx * std::sin(phi) + dy * std::cos(phi); }
0110     float calcDz(reco::Candidate::Point p, reco::Candidate::Point v, const reco::Candidate &c) const {
0111       return p.Z() - v.Z() - ((p.X() - v.X()) * c.px() + (p.Y() - v.Y()) * c.py()) * c.pz() / (c.pt() * c.pt());
0112     }
0113   };
0114 }  // namespace pat
0115 
0116 pat::PATPackedCandidateProducer::PATPackedCandidateProducer(const edm::ParameterSet &iConfig)
0117     : usePuppi_(!iConfig.getParameter<edm::InputTag>("PuppiSrc").encode().empty() ||
0118                 !iConfig.getParameter<edm::InputTag>("PuppiNoLepSrc").encode().empty()),
0119       Cands_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("inputCollection"))),
0120       PVs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("inputVertices"))),
0121       PVAsso_(
0122           consumes<edm::Association<reco::VertexCollection>>(iConfig.getParameter<edm::InputTag>("vertexAssociator"))),
0123       PVAssoQuality_(consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("vertexAssociator"))),
0124       PVOrigs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("originalVertices"))),
0125       TKOrigs_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("originalTracks"))),
0126       PuppiWeight_(usePuppi_ ? consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("PuppiSrc"))
0127                              : edm::EDGetTokenT<edm::ValueMap<float>>()),
0128       PuppiWeightNoLep_(usePuppi_ ? consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("PuppiNoLepSrc"))
0129                                   : edm::EDGetTokenT<edm::ValueMap<float>>()),
0130       storeChargedHadronIsolation_(!iConfig.getParameter<edm::InputTag>("chargedHadronIsolation").encode().empty()),
0131       ChargedHadronIsolation_(
0132           consumes<edm::ValueMap<bool>>(iConfig.getParameter<edm::InputTag>("chargedHadronIsolation"))),
0133       minPtForChargedHadronProperties_(iConfig.getParameter<double>("minPtForChargedHadronProperties")),
0134       minPtForTrackProperties_(iConfig.getParameter<double>("minPtForTrackProperties")),
0135       minPtForLowQualityTrackProperties_(iConfig.getParameter<double>("minPtForLowQualityTrackProperties")),
0136       covarianceVersion_(iConfig.getParameter<int>("covarianceVersion")),
0137       covariancePackingSchemas_(iConfig.getParameter<std::vector<int>>("covariancePackingSchemas")),
0138       pfCandidateTypesForHcalDepth_(iConfig.getParameter<std::vector<int>>("pfCandidateTypesForHcalDepth")),
0139       storeHcalDepthEndcapOnly_(iConfig.getParameter<bool>("storeHcalDepthEndcapOnly")),
0140       storeTiming_(iConfig.getParameter<bool>("storeTiming")),
0141       timeFromValueMap_(!iConfig.getParameter<edm::InputTag>("timeMap").encode().empty() &&
0142                         !iConfig.getParameter<edm::InputTag>("timeMapErr").encode().empty()),
0143       t0Map_(timeFromValueMap_ ? consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeMap"))
0144                                : edm::EDGetTokenT<edm::ValueMap<float>>()),
0145       t0ErrMap_(timeFromValueMap_ ? consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeMapErr"))
0146                                   : edm::EDGetTokenT<edm::ValueMap<float>>()) {
0147   std::vector<edm::InputTag> sv_tags =
0148       iConfig.getParameter<std::vector<edm::InputTag>>("secondaryVerticesForWhiteList");
0149   for (const auto &itag : sv_tags) {
0150     SVWhiteLists_.push_back(consumes<edm::View<reco::Candidate>>(itag));
0151   }
0152 
0153   produces<std::vector<pat::PackedCandidate>>();
0154   produces<edm::Association<pat::PackedCandidateCollection>>();
0155   produces<edm::Association<reco::PFCandidateCollection>>();
0156 
0157   if (not pfCandidateTypesForHcalDepth_.empty())
0158     produces<edm::ValueMap<pat::HcalDepthEnergyFractions>>("hcalDepthEnergyFractions");
0159 }
0160 
0161 pat::PATPackedCandidateProducer::~PATPackedCandidateProducer() {}
0162 
0163 void pat::PATPackedCandidateProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0164   edm::Handle<reco::PFCandidateCollection> cands;
0165   iEvent.getByToken(Cands_, cands);
0166 
0167   edm::Handle<edm::ValueMap<float>> puppiWeight;
0168   edm::Handle<edm::ValueMap<float>> puppiWeightNoLep;
0169   if (usePuppi_) {
0170     iEvent.getByToken(PuppiWeight_, puppiWeight);
0171     iEvent.getByToken(PuppiWeightNoLep_, puppiWeightNoLep);
0172   }
0173 
0174   edm::Handle<reco::VertexCollection> PVOrigs;
0175   iEvent.getByToken(PVOrigs_, PVOrigs);
0176 
0177   edm::Handle<edm::Association<reco::VertexCollection>> assoHandle;
0178   iEvent.getByToken(PVAsso_, assoHandle);
0179   edm::Handle<edm::ValueMap<int>> assoQualityHandle;
0180   iEvent.getByToken(PVAssoQuality_, assoQualityHandle);
0181   const edm::Association<reco::VertexCollection> &associatedPV = *(assoHandle.product());
0182   const edm::ValueMap<int> &associationQuality = *(assoQualityHandle.product());
0183 
0184   edm::Handle<edm::ValueMap<bool>> chargedHadronIsolationHandle;
0185   if (storeChargedHadronIsolation_)
0186     iEvent.getByToken(ChargedHadronIsolation_, chargedHadronIsolationHandle);
0187 
0188   std::set<unsigned int> whiteList;
0189   std::set<reco::TrackRef> whiteListTk;
0190   for (auto itoken : SVWhiteLists_) {
0191     edm::Handle<edm::View<reco::Candidate>> svWhiteListHandle;
0192     iEvent.getByToken(itoken, svWhiteListHandle);
0193     const edm::View<reco::Candidate> &svWhiteList = *(svWhiteListHandle.product());
0194     for (unsigned int i = 0; i < svWhiteList.size(); i++) {
0195       // Whitelist via Ptrs
0196       for (unsigned int j = 0; j < svWhiteList[i].numberOfSourceCandidatePtrs(); j++) {
0197         const edm::Ptr<reco::Candidate> &c = svWhiteList[i].sourceCandidatePtr(j);
0198         if (c.id() == cands.id())
0199           whiteList.insert(c.key());
0200       }
0201       // Whitelist via RecoCharged
0202       for (auto dau = svWhiteList[i].begin(); dau != svWhiteList[i].end(); dau++) {
0203         const reco::RecoChargedCandidate *chCand = dynamic_cast<const reco::RecoChargedCandidate *>(&(*dau));
0204         if (chCand != nullptr) {
0205           whiteListTk.insert(chCand->track());
0206         }
0207       }
0208     }
0209   }
0210 
0211   edm::Handle<edm::ValueMap<float>> t0Map;
0212   edm::Handle<edm::ValueMap<float>> t0ErrMap;
0213   if (timeFromValueMap_) {
0214     iEvent.getByToken(t0Map_, t0Map);
0215     iEvent.getByToken(t0ErrMap_, t0ErrMap);
0216   }
0217 
0218   edm::Handle<reco::VertexCollection> PVs;
0219   iEvent.getByToken(PVs_, PVs);
0220   reco::VertexRef PV(PVs.id());
0221   reco::VertexRefProd PVRefProd(PVs);
0222   math::XYZPoint PVpos;
0223 
0224   std::vector<pat::HcalDepthEnergyFractions> hcalDepthEnergyFractions;
0225   hcalDepthEnergyFractions.reserve(cands->size());
0226   std::vector<pat::HcalDepthEnergyFractions> hcalDepthEnergyFractions_Ordered;
0227   hcalDepthEnergyFractions_Ordered.reserve(cands->size());
0228 
0229   edm::Handle<reco::TrackCollection> TKOrigs;
0230   iEvent.getByToken(TKOrigs_, TKOrigs);
0231   auto outPtrP = std::make_unique<std::vector<pat::PackedCandidate>>();
0232   std::vector<int> mapping(cands->size());
0233   std::vector<int> mappingReverse(cands->size());
0234   std::vector<int> mappingTk(TKOrigs->size(), -1);
0235 
0236   for (unsigned int ic = 0, nc = cands->size(); ic < nc; ++ic) {
0237     const reco::PFCandidate &cand = (*cands)[ic];
0238     const reco::Track *ctrack = nullptr;
0239     if ((abs(cand.pdgId()) == 11 || cand.pdgId() == 22) && cand.gsfTrackRef().isNonnull()) {
0240       ctrack = &*cand.gsfTrackRef();
0241     } else if (cand.trackRef().isNonnull()) {
0242       ctrack = &*cand.trackRef();
0243     }
0244     if (ctrack) {
0245       float dist = 1e99;
0246       int pvi = -1;
0247       for (size_t ii = 0; ii < PVs->size(); ii++) {
0248         float dz = std::abs(ctrack->dz(((*PVs)[ii]).position()));
0249         if (dz < dist) {
0250           pvi = ii;
0251           dist = dz;
0252         }
0253       }
0254       PV = reco::VertexRef(PVs, pvi);
0255       math::XYZPoint vtx = cand.vertex();
0256       pat::PackedCandidate::LostInnerHits lostHits = pat::PackedCandidate::noLostInnerHits;
0257       const reco::VertexRef &PVOrig = associatedPV[reco::CandidatePtr(cands, ic)];
0258       if (PVOrig.isNonnull())
0259         PV = reco::VertexRef(PVs,
0260                              PVOrig.key());  // WARNING: assume the PV slimmer is keeping same order
0261       int quality = associationQuality[reco::CandidatePtr(cands, ic)];
0262       //          if ((size_t)pvi!=PVOrig.key()) std::cout << "not closest in Z"
0263       //          << pvi << " " << PVOrig.key() << " " << cand.pt() << " " <<
0264       //          quality << std::endl; TrajectoryStateOnSurface tsos =
0265       //          extrapolator.extrapolate(trajectoryStateTransform::initialFreeState(*ctrack,&*magneticField),
0266       //          RecoVertex::convertPos(PV->position()));
0267       //   vtx = tsos.globalPosition();
0268       //          phiAtVtx = tsos.globalDirection().phi();
0269       vtx = ctrack->referencePoint();
0270       float ptTrk = ctrack->pt();
0271       float etaAtVtx = ctrack->eta();
0272       float phiAtVtx = ctrack->phi();
0273 
0274       int nlost = ctrack->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0275       if (nlost == 0) {
0276         if (ctrack->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
0277           lostHits = pat::PackedCandidate::validHitInFirstPixelBarrelLayer;
0278         }
0279       } else {
0280         lostHits = (nlost == 1 ? pat::PackedCandidate::oneLostInnerHit : pat::PackedCandidate::moreLostInnerHits);
0281       }
0282 
0283       outPtrP->push_back(
0284           pat::PackedCandidate(cand.polarP4(), vtx, ptTrk, etaAtVtx, phiAtVtx, cand.pdgId(), PVRefProd, PV.key()));
0285       outPtrP->back().setAssociationQuality(pat::PackedCandidate::PVAssociationQuality(qualityMap[quality]));
0286       outPtrP->back().setCovarianceVersion(covarianceVersion_);
0287       if (cand.trackRef().isNonnull() && PVOrig.isNonnull() && PVOrig->trackWeight(cand.trackRef()) > 0.5 &&
0288           quality == 7) {
0289         outPtrP->back().setAssociationQuality(pat::PackedCandidate::UsedInFitTight);
0290       }
0291       // properties of the best track
0292       outPtrP->back().setLostInnerHits(lostHits);
0293       if (outPtrP->back().pt() > minPtForTrackProperties_ || outPtrP->back().ptTrk() > minPtForTrackProperties_ ||
0294           whiteList.find(ic) != whiteList.end() ||
0295           (cand.trackRef().isNonnull() && whiteListTk.find(cand.trackRef()) != whiteListTk.end())) {
0296         outPtrP->back().setFirstHit(ctrack->hitPattern().getHitPattern(reco::HitPattern::TRACK_HITS, 0));
0297         if (abs(outPtrP->back().pdgId()) == 22) {
0298           outPtrP->back().setTrackProperties(*ctrack, covariancePackingSchemas_[4], covarianceVersion_);
0299         } else {
0300           if (ctrack->hitPattern().numberOfValidPixelHits() > 0) {
0301             outPtrP->back().setTrackProperties(*ctrack,
0302                                                covariancePackingSchemas_[0],
0303                                                covarianceVersion_);  // high quality
0304           } else {
0305             outPtrP->back().setTrackProperties(*ctrack, covariancePackingSchemas_[1], covarianceVersion_);
0306           }
0307         }
0308         // outPtrP->back().setTrackProperties(*ctrack,tsos.curvilinearError());
0309       } else {
0310         if (outPtrP->back().pt() > minPtForLowQualityTrackProperties_) {
0311           if (ctrack->hitPattern().numberOfValidPixelHits() > 0)
0312             outPtrP->back().setTrackProperties(*ctrack,
0313                                                covariancePackingSchemas_[2],
0314                                                covarianceVersion_);  // low quality, with pixels
0315           else
0316             outPtrP->back().setTrackProperties(*ctrack,
0317                                                covariancePackingSchemas_[3],
0318                                                covarianceVersion_);  // low quality, without pixels
0319         }
0320       }
0321 
0322       // these things are always for the CKF track
0323       outPtrP->back().setTrackHighPurity(cand.trackRef().isNonnull() &&
0324                                          cand.trackRef()->quality(reco::Track::highPurity));
0325       if (cand.muonRef().isNonnull()) {
0326         outPtrP->back().setMuonID(cand.muonRef()->isStandAloneMuon(), cand.muonRef()->isGlobalMuon());
0327       }
0328     } else {
0329       if (!PVs->empty()) {
0330         PV = reco::VertexRef(PVs, 0);
0331         PVpos = PV->position();
0332       }
0333 
0334       outPtrP->push_back(pat::PackedCandidate(
0335           cand.polarP4(), PVpos, cand.pt(), cand.eta(), cand.phi(), cand.pdgId(), PVRefProd, PV.key()));
0336       outPtrP->back().setAssociationQuality(
0337           pat::PackedCandidate::PVAssociationQuality(pat::PackedCandidate::UsedInFitTight));
0338     }
0339 
0340     // neutrals and isolated charged hadrons
0341 
0342     bool isIsolatedChargedHadron = false;
0343     if (storeChargedHadronIsolation_) {
0344       const edm::ValueMap<bool> &chargedHadronIsolation = *(chargedHadronIsolationHandle.product());
0345       isIsolatedChargedHadron =
0346           ((cand.pt() > minPtForChargedHadronProperties_) && (chargedHadronIsolation[reco::PFCandidateRef(cands, ic)]));
0347       outPtrP->back().setIsIsolatedChargedHadron(isIsolatedChargedHadron);
0348     }
0349 
0350     if (abs(cand.pdgId()) == 1 || abs(cand.pdgId()) == 130) {
0351       outPtrP->back().setHcalFraction(cand.hcalEnergy() / (cand.ecalEnergy() + cand.hcalEnergy()));
0352     } else if ((cand.charge() || abs(cand.pdgId()) == 22) && cand.pt() > 0.5) {
0353       outPtrP->back().setHcalFraction(cand.hcalEnergy() / (cand.ecalEnergy() + cand.hcalEnergy()));
0354       outPtrP->back().setCaloFraction((cand.hcalEnergy() + cand.ecalEnergy()) / cand.energy());
0355     } else {
0356       outPtrP->back().setHcalFraction(0);
0357       outPtrP->back().setCaloFraction(0);
0358     }
0359 
0360     if (isIsolatedChargedHadron) {
0361       outPtrP->back().setRawCaloFraction((cand.rawEcalEnergy() + cand.rawHcalEnergy()) / cand.energy());
0362       outPtrP->back().setRawHcalFraction(cand.rawHcalEnergy() / (cand.rawEcalEnergy() + cand.rawHcalEnergy()));
0363     } else {
0364       outPtrP->back().setRawCaloFraction(0);
0365       outPtrP->back().setRawHcalFraction(0);
0366     }
0367 
0368     std::vector<float> dummyVector;
0369     dummyVector.clear();
0370     pat::HcalDepthEnergyFractions hcalDepthEFrac(dummyVector);
0371 
0372     // storing HcalDepthEnergyFraction information
0373     if (std::find(pfCandidateTypesForHcalDepth_.begin(), pfCandidateTypesForHcalDepth_.end(), abs(cand.pdgId())) !=
0374         pfCandidateTypesForHcalDepth_.end()) {
0375       if (!storeHcalDepthEndcapOnly_ ||
0376           fabs(outPtrP->back().eta()) > 1.3) {  // storeHcalDepthEndcapOnly_==false -> store all eta of
0377                                                 // selected PF types, if true, only |eta|>1.3 of selected
0378                                                 // PF types will be stored
0379         std::vector<float> hcalDepthEnergyFractionTmp(cand.hcalDepthEnergyFractions().begin(),
0380                                                       cand.hcalDepthEnergyFractions().end());
0381         hcalDepthEFrac.reset(hcalDepthEnergyFractionTmp);
0382       }
0383     }
0384     hcalDepthEnergyFractions.push_back(hcalDepthEFrac);
0385 
0386     // specifically this is the PFLinker requirements to apply the e/gamma
0387     // regression
0388     if (cand.particleId() == reco::PFCandidate::e ||
0389         (cand.particleId() == reco::PFCandidate::gamma && cand.mva_nothing_gamma() > 0.)) {
0390       outPtrP->back().setGoodEgamma();
0391     }
0392 
0393     if (usePuppi_) {
0394       reco::PFCandidateRef pkref(cands, ic);
0395 
0396       float puppiWeightVal = (*puppiWeight)[pkref];
0397       float puppiWeightNoLepVal = (*puppiWeightNoLep)[pkref];
0398       outPtrP->back().setPuppiWeight(puppiWeightVal, puppiWeightNoLepVal);
0399     }
0400 
0401     if (storeTiming_) {
0402       if (timeFromValueMap_) {
0403         if (cand.trackRef().isNonnull()) {
0404           auto t0 = (*t0Map)[cand.trackRef()];
0405           auto t0Err = (*t0ErrMap)[cand.trackRef()];
0406           outPtrP->back().setTime(t0, t0Err);
0407         }
0408       } else {
0409         if (cand.isTimeValid()) {
0410           outPtrP->back().setTime(cand.time(), cand.timeError());
0411         }
0412       }
0413     }
0414 
0415     mapping[ic] = ic;  // trivial at the moment!
0416     if (cand.trackRef().isNonnull() && cand.trackRef().id() == TKOrigs.id()) {
0417       mappingTk[cand.trackRef().key()] = ic;
0418     }
0419   }
0420 
0421   auto outPtrPSorted = std::make_unique<std::vector<pat::PackedCandidate>>();
0422   std::vector<size_t> order = sort_indexes(*outPtrP);
0423   std::vector<size_t> reverseOrder(order.size());
0424   for (size_t i = 0, nc = cands->size(); i < nc; i++) {
0425     outPtrPSorted->push_back((*outPtrP)[order[i]]);
0426     reverseOrder[order[i]] = i;
0427     mappingReverse[order[i]] = i;
0428     hcalDepthEnergyFractions_Ordered.push_back(hcalDepthEnergyFractions[order[i]]);
0429   }
0430 
0431   // Fix track association for sorted candidates
0432   for (size_t i = 0, ntk = mappingTk.size(); i < ntk; i++) {
0433     if (mappingTk[i] >= 0)
0434       mappingTk[i] = reverseOrder[mappingTk[i]];
0435   }
0436 
0437   edm::OrphanHandle<pat::PackedCandidateCollection> oh = iEvent.put(std::move(outPtrPSorted));
0438 
0439   // now build the two maps
0440   auto pf2pc = std::make_unique<edm::Association<pat::PackedCandidateCollection>>(oh);
0441   auto pc2pf = std::make_unique<edm::Association<reco::PFCandidateCollection>>(cands);
0442   edm::Association<pat::PackedCandidateCollection>::Filler pf2pcFiller(*pf2pc);
0443   edm::Association<reco::PFCandidateCollection>::Filler pc2pfFiller(*pc2pf);
0444   pf2pcFiller.insert(cands, mappingReverse.begin(), mappingReverse.end());
0445   pc2pfFiller.insert(oh, order.begin(), order.end());
0446   // include also the mapping track -> packed PFCand
0447   pf2pcFiller.insert(TKOrigs, mappingTk.begin(), mappingTk.end());
0448 
0449   pf2pcFiller.fill();
0450   pc2pfFiller.fill();
0451   iEvent.put(std::move(pf2pc));
0452   iEvent.put(std::move(pc2pf));
0453 
0454   // HCAL depth energy fraction additions using ValueMap
0455   auto hcalDepthEnergyFractionsV = std::make_unique<edm::ValueMap<HcalDepthEnergyFractions>>();
0456   edm::ValueMap<HcalDepthEnergyFractions>::Filler fillerHcalDepthEnergyFractions(*hcalDepthEnergyFractionsV);
0457   fillerHcalDepthEnergyFractions.insert(
0458       cands, hcalDepthEnergyFractions_Ordered.begin(), hcalDepthEnergyFractions_Ordered.end());
0459   fillerHcalDepthEnergyFractions.fill();
0460 
0461   if (not pfCandidateTypesForHcalDepth_.empty())
0462     iEvent.put(std::move(hcalDepthEnergyFractionsV), "hcalDepthEnergyFractions");
0463 }
0464 
0465 using pat::PATPackedCandidateProducer;
0466 #include "FWCore/Framework/interface/MakerMacros.h"
0467 DEFINE_FWK_MODULE(PATPackedCandidateProducer);