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