TICL/plugins/PFTICLProducer.cc

0001 // This producer converts a list of TICLCandidates to a list of PFCandidates.
0002
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008
0009 #include "DataFormats/Common/interface/View.h"
0010
0011 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0012 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0013
0014 #include "DataFormats/HGCalReco/interface/TICLCandidate.h"
0015
0016 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0017
0018 class PFTICLProducer : public edm::stream::EDProducer<> {
0019 public:
0020   PFTICLProducer(const edm::ParameterSet&);
0021   ~PFTICLProducer() override {}
0022
0023   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0024
0025   void produce(edm::Event&, const edm::EventSetup&) override;
0026
0027 private:
0028   // parameters
0029   const bool useMTDTiming_;
0030   const bool useTimingAverage_;
0031   const float timingQualityThreshold_;
0032   const bool energy_from_regression_;
0033   // inputs
0034   const edm::EDGetTokenT<edm::View<TICLCandidate>> ticl_candidates_;
0035   const edm::EDGetTokenT<edm::ValueMap<float>> srcTrackTime_, srcTrackTimeError_, srcTrackTimeQuality_;
0036   const edm::EDGetTokenT<reco::MuonCollection> muons_;
0037   // For PFMuonAlgo
0038   std::unique_ptr<PFMuonAlgo> pfmu_;
0039 };
0040
0041 DEFINE_FWK_MODULE(PFTICLProducer);
0042
0043 PFTICLProducer::PFTICLProducer(const edm::ParameterSet& conf)
0044     : useMTDTiming_(conf.getParameter<bool>("useMTDTiming")),
0045       useTimingAverage_(conf.getParameter<bool>("useTimingAverage")),
0046       timingQualityThreshold_(conf.getParameter<double>("timingQualityThreshold")),
0047       energy_from_regression_(conf.getParameter<bool>("energyFromRegression")),
0048       ticl_candidates_(consumes<edm::View<TICLCandidate>>(conf.getParameter<edm::InputTag>("ticlCandidateSrc"))),
0049       srcTrackTime_(consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeValueMap"))),
0050       srcTrackTimeError_(consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeErrorMap"))),
0051       srcTrackTimeQuality_(consumes<edm::ValueMap<float>>(conf.getParameter<edm::InputTag>("trackTimeQualityMap"))),
0052       muons_(consumes<reco::MuonCollection>(conf.getParameter<edm::InputTag>("muonSrc"))),
0053       pfmu_(std::make_unique<PFMuonAlgo>(conf.getParameterSet("pfMuonAlgoParameters"),
0054                                          false)) {  // postMuonCleaning = false
0055   produces<reco::PFCandidateCollection>();
0056 }
0057
0058 void PFTICLProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0059   edm::ParameterSetDescription desc;
0060   desc.add<edm::InputTag>("ticlCandidateSrc", edm::InputTag("ticlTrackstersMerge"));
0061   desc.add<edm::InputTag>("trackTimeValueMap", edm::InputTag("tofPID:t0"));
0062   desc.add<edm::InputTag>("trackTimeErrorMap", edm::InputTag("tofPID:sigmat0"));
0063   desc.add<edm::InputTag>("trackTimeQualityMap", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
0064   desc.add<bool>("energyFromRegression", true);
0065   desc.add<double>("timingQualityThreshold", 0.5);
0066   desc.add<bool>("useMTDTiming", true);
0067   desc.add<bool>("useTimingAverage", false);
0068   // For PFMuonAlgo
0069   desc.add<edm::InputTag>("muonSrc", edm::InputTag("muons1stStep"));
0070   edm::ParameterSetDescription psd_PFMuonAlgo;
0071   PFMuonAlgo::fillPSetDescription(psd_PFMuonAlgo);
0072   desc.add<edm::ParameterSetDescription>("pfMuonAlgoParameters", psd_PFMuonAlgo);
0073   //
0074   descriptions.add("pfTICLProducer", desc);
0075 }
0076
0077 void PFTICLProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0078   //get TICLCandidates
0079   edm::Handle<edm::View<TICLCandidate>> ticl_cand_h;
0080   evt.getByToken(ticl_candidates_, ticl_cand_h);
0081   const auto ticl_candidates = *ticl_cand_h;
0082   edm::Handle<edm::ValueMap<float>> trackTimeH, trackTimeErrH, trackTimeQualH;
0083   evt.getByToken(srcTrackTime_, trackTimeH);
0084   evt.getByToken(srcTrackTimeError_, trackTimeErrH);
0085   evt.getByToken(srcTrackTimeQuality_, trackTimeQualH);
0086   const auto muonH = evt.getHandle(muons_);
0087   const auto muons = *muonH;
0088
0089   auto candidates = std::make_unique<reco::PFCandidateCollection>();
0090
0091   for (const auto& ticl_cand : ticl_candidates) {
0092     const auto abs_pdg_id = std::abs(ticl_cand.pdgId());
0093     const auto charge = ticl_cand.charge();
0094     const auto& four_mom = ticl_cand.p4();
0095     float total_raw_energy = 0.f;
0096     float total_em_raw_energy = 0.f;
0097     for (const auto& t : ticl_cand.tracksters()) {
0098       total_raw_energy += t->raw_energy();
0099       total_em_raw_energy += t->raw_em_energy();
0100     }
0101     float ecal_energy_fraction = total_em_raw_energy / total_raw_energy;
0102     float ecal_energy = energy_from_regression_ ? ticl_cand.p4().energy() * ecal_energy_fraction
0103                                                 : ticl_cand.rawEnergy() * ecal_energy_fraction;
0104     float hcal_energy =
0105         energy_from_regression_ ? ticl_cand.p4().energy() - ecal_energy : ticl_cand.rawEnergy() - ecal_energy;
0106     // fix for floating point rounding could go slightly below 0
0107     hcal_energy = std::max(0.f, hcal_energy);
0108     reco::PFCandidate::ParticleType part_type;
0109     switch (abs_pdg_id) {
0110       case 11:
0111         part_type = reco::PFCandidate::e;
0112         break;
0113       case 13:
0114         part_type = reco::PFCandidate::mu;
0115         break;
0116       case 22:
0117         part_type = reco::PFCandidate::gamma;
0118         break;
0119       case 130:
0120         part_type = reco::PFCandidate::h0;
0121         break;
0122       case 211:
0123         part_type = reco::PFCandidate::h;
0124         break;
0125       // default also handles neutral pions (111) for the time being (not yet foreseen in PFCandidate)
0126       default:
0127         part_type = reco::PFCandidate::X;
0128     }
0129
0130     candidates->emplace_back(charge, four_mom, part_type);
0131
0132     auto& candidate = candidates->back();
0133     candidate.setEcalEnergy(ecal_energy, ecal_energy);
0134     candidate.setHcalEnergy(hcal_energy, hcal_energy);
0135     if (candidate.charge()) {  // otherwise PFCandidate throws
0136       // Construct edm::Ref from edm::Ptr. As of now, assumes type to be reco::Track. To be extended (either via
0137       // dynamic type checking or configuration) if additional track types are needed.
0138       reco::TrackRef trackref(ticl_cand.trackPtr().id(), int(ticl_cand.trackPtr().key()), &evt.productGetter());
0139       candidate.setTrackRef(trackref);
0140       // Utilize PFMuonAlgo
0141       const int muId = PFMuonAlgo::muAssocToTrack(trackref, muons);
0142       if (muId != -1) {
0143         const reco::MuonRef muonref = reco::MuonRef(muonH, muId);
0144         const bool allowLoose = (part_type == reco::PFCandidate::mu);
0145         // Redefine pfmuon candidate kinematics and add muonref
0146         pfmu_->reconstructMuon(candidate, muonref, allowLoose);
0147       }
0148     }
0149
0150     // HGCAL timing as default values
0151     auto time = ticl_cand.time();
0152     auto timeE = ticl_cand.timeError();
0153
0154     if (useMTDTiming_ and candidate.charge()) {
0155       // Ignore HGCAL timing until it will be TOF corrected
0156       time = -99.;
0157       timeE = -1.;
0158       // Check MTD timing availability
0159       const bool assocQuality = (*trackTimeQualH)[candidate.trackRef()] > timingQualityThreshold_;
0160       if (assocQuality) {
0161         const auto timeHGC = time;
0162         const auto timeEHGC = timeE;
0163         const auto timeMTD = (*trackTimeH)[candidate.trackRef()];
0164         const auto timeEMTD = (*trackTimeErrH)[candidate.trackRef()];
0165
0166         if (useTimingAverage_ && (timeEMTD > 0 && timeEHGC > 0)) {
0167           // Compute weighted average between HGCAL and MTD timing
0168           const auto invTimeESqHGC = pow(timeEHGC, -2);
0169           const auto invTimeESqMTD = pow(timeEMTD, -2);
0170           timeE = (invTimeESqHGC * invTimeESqMTD) / (invTimeESqHGC + invTimeESqMTD);
0171           time = (timeHGC * invTimeESqHGC + timeMTD * invTimeESqMTD) * timeE;
0172           timeE = sqrt(timeE);
0173         } else if (timeEMTD > 0) {  // Ignore HGCal timing until it will be TOF corrected
0174           time = timeMTD;
0175           timeE = timeEMTD;
0176         }
0177       }
0178     }
0179     candidate.setTime(time, timeE);
0180   }
0181
0182   evt.put(std::move(candidates));
0183 }