diff --git a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc index 39c248356a0d8..5f64b43b10ccc 100644 --- a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc +++ b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc @@ -141,6 +141,7 @@ SAMuon Phase2L1TGMTSAMuonProducer::Convertl1tMuon(const l1t::Muon& mu, const int int bstart = 0; wordtype word(0); + bstart = wordconcat(word, bstart, pt > 0, 1); bstart = wordconcat(word, bstart, pt, BITSGTPT); bstart = wordconcat(word, bstart, phi, BITSGTPHI); bstart = wordconcat(word, bstart, eta, BITSGTETA); diff --git a/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml b/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml index 28ac6a5e6e71c..d3bea1e062d62 100644 --- a/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml +++ b/L1Trigger/Phase2L1ParticleFlow/BuildFile.xml @@ -1,15 +1,14 @@ + - - - + diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/DiscretePFInputs.h b/L1Trigger/Phase2L1ParticleFlow/interface/DiscretePFInputs.h deleted file mode 100644 index 67f55afc2026c..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/DiscretePFInputs.h +++ /dev/null @@ -1,338 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_H -#define L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_H - -#if defined(__GXX_EXPERIMENTAL_CXX0X__) or defined(CMSSW) -#include -#include -#define L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE -#else -#include -#endif - -namespace l1t { - class PFTrack; - class PFCluster; - class PFCandidate; - class Muon; -} // namespace l1t - -// the serialization may be hidden if needed -#include -#include - -namespace l1tpf_impl { - - struct CaloCluster { - int16_t hwPt; - int16_t hwEmPt; - int16_t hwPtErr; - int16_t hwEta; - int16_t hwPhi; - uint16_t hwFlags; - bool isEM, used; - const l1t::PFCluster *src; - - // sorting - bool operator<(const CaloCluster &other) const { return hwPt > other.hwPt; } - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - static constexpr float PT_SCALE = 4.0; // quantize in units of 0.25 GeV (can be changed) - static constexpr float ETAPHI_FACTOR = 4; // size of an ecal crystal in phi in integer units (our choice) - static constexpr float ETAPHI_SCALE = - ETAPHI_FACTOR * - (180. / M_PI); // M_PI/180 is the size of an ECal crystal; we make a grid that is 4 times that size - static constexpr int16_t PHI_WRAP = 360 * ETAPHI_FACTOR; // what is 3.14 in integer - - static int16_t ptToInt16(float pt) { // avoid overflows - return std::min(round(pt * CaloCluster::PT_SCALE), std::numeric_limits::max()); - } - - // filling from floating point - void fill(float pt, - float emPt, - float ptErr, - float eta, - float phi, - bool em, - unsigned int flags, - const l1t::PFCluster *source = nullptr) { - hwPt = CaloCluster::ptToInt16(pt); - hwEmPt = CaloCluster::ptToInt16(emPt); - hwPtErr = CaloCluster::ptToInt16(ptErr); - hwEta = round(eta * CaloCluster::ETAPHI_SCALE); - hwPhi = int16_t(round(phi * CaloCluster::ETAPHI_SCALE)) % CaloCluster::PHI_WRAP; - isEM = em; - used = false; - hwFlags = flags; - src = source; - } - - float floatPt() const { return float(hwPt) / CaloCluster::PT_SCALE; } - float floatEmPt() const { return float(hwEmPt) / CaloCluster::PT_SCALE; } - float floatPtErr() const { return float(hwPtErr) / CaloCluster::PT_SCALE; } - static float minFloatPt() { return float(1.0) / CaloCluster::PT_SCALE; } - float floatEta() const { return float(hwEta) / CaloCluster::ETAPHI_SCALE; } - float floatPhi() const { return float(hwPhi) / CaloCluster::ETAPHI_SCALE; } - void setFloatPt(float pt) { hwPt = round(pt * CaloCluster::PT_SCALE); } - void setFloatEmPt(float emPt) { hwEmPt = round(emPt * CaloCluster::PT_SCALE); } -#endif - }; - - // https://twiki.cern.ch/twiki/bin/view/CMS/L1TriggerPhase2InterfaceSpecifications - struct InputTrack { - uint16_t hwInvpt; - int32_t hwVtxEta; - int32_t hwVtxPhi; - bool hwCharge; - int16_t hwZ0; - uint16_t hwChi2, hwStubs; - uint16_t hwFlags; - const l1t::PFTrack *src; - - enum QualityFlags { PFLOOSE = 1, PFTIGHT = 2, TKEG = 4 }; - bool quality(QualityFlags q) const { return hwFlags & q; } - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - static constexpr float INVPT_SCALE = 2E4; // 1%/pt @ 100 GeV is 2 bits - static constexpr float VTX_PHI_SCALE = 1 / 1.6E-3; // 5 micro rad is 2 bits - static constexpr float VTX_ETA_SCALE = 1 / 1E-4; // no idea, but assume it's somewhat worse than phi - static constexpr float Z0_SCALE = 20; // 1mm is 2 bits - static constexpr int32_t VTX_ETA_1p3 = 1.3 * InputTrack::VTX_ETA_SCALE; - - // filling from floating point - void fillInput( - float pt, float eta, float phi, int charge, float dz, unsigned int flags, const l1t::PFTrack *source = nullptr) { - hwInvpt = std::min(round(1 / pt * InputTrack::INVPT_SCALE), std::numeric_limits::max()); - hwVtxEta = round(eta * InputTrack::VTX_ETA_SCALE); - hwVtxPhi = round(phi * InputTrack::VTX_PHI_SCALE); - hwCharge = (charge > 0); - hwZ0 = round(dz * InputTrack::Z0_SCALE); - hwFlags = flags; - src = source; - } - - float floatVtxPt() const { return 1 / (float(hwInvpt) / InputTrack::INVPT_SCALE); } - float floatVtxEta() const { return float(hwVtxEta) / InputTrack::VTX_ETA_SCALE; } - float floatVtxPhi() const { return float(hwVtxPhi) / InputTrack::VTX_PHI_SCALE; } - float floatDZ() const { return float(hwZ0) / InputTrack::Z0_SCALE; } - int intCharge() const { return hwCharge ? +1 : -1; } -#endif - }; - - struct PropagatedTrack : public InputTrack { - int16_t hwPt; - int16_t hwPtErr; - int16_t hwCaloPtErr; - int16_t hwEta; // at calo - int16_t hwPhi; // at calo - bool muonLink; - bool used; // note: this flag is not used in the default PF, but is used in alternative algos - bool fromPV; - - // sorting - bool operator<(const PropagatedTrack &other) const { return hwPt > other.hwPt; } - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - void fillPropagated( - float pt, float ptErr, float caloPtErr, float caloEta, float caloPhi, unsigned int quality, bool isMuon) { - hwPt = CaloCluster::ptToInt16(pt); - hwPtErr = CaloCluster::ptToInt16(ptErr); - hwCaloPtErr = CaloCluster::ptToInt16(caloPtErr); - // saturation protection - if (hwPt == std::numeric_limits::max()) { - hwCaloPtErr = hwPt / 4; - } - hwEta = round(caloEta * CaloCluster::ETAPHI_SCALE); - hwPhi = int16_t(round(caloPhi * CaloCluster::ETAPHI_SCALE)) % CaloCluster::PHI_WRAP; - muonLink = isMuon; - used = false; - } - - float floatPt() const { return float(hwPt) / CaloCluster::PT_SCALE; } - float floatPtErr() const { return float(hwPtErr) / CaloCluster::PT_SCALE; } - float floatCaloPtErr() const { return float(hwCaloPtErr) / CaloCluster::PT_SCALE; } - float floatEta() const { return float(hwEta) / CaloCluster::ETAPHI_SCALE; } - float floatPhi() const { return float(hwPhi) / CaloCluster::ETAPHI_SCALE; } -#endif - }; - - struct Muon { - int16_t hwPt; - int16_t hwEta; // at calo - int16_t hwPhi; // at calo - uint16_t hwFlags; - bool hwCharge; - const l1t::Muon *src; - - // sorting - bool operator<(const Muon &other) const { return hwPt > other.hwPt; } - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - void fill(float pt, float eta, float phi, int charge, unsigned int flags, const l1t::Muon *source = nullptr) { - // we assume we use the same discrete ieta, iphi grid for all particles - hwPt = round(pt * CaloCluster::PT_SCALE); - hwEta = round(eta * CaloCluster::ETAPHI_SCALE); - hwPhi = int16_t(round(phi * CaloCluster::ETAPHI_SCALE)) % CaloCluster::PHI_WRAP; - hwCharge = (charge > 0); - hwFlags = flags; - src = source; - } - float floatPt() const { return float(hwPt) / CaloCluster::PT_SCALE; } - float floatEta() const { return float(hwEta) / CaloCluster::ETAPHI_SCALE; } - float floatPhi() const { return float(hwPhi) / CaloCluster::ETAPHI_SCALE; } - int intCharge() const { return hwCharge ? +1 : -1; } -#endif - }; - - struct PFParticle { - int16_t hwPt; - int16_t hwEta; // at calo face - int16_t hwPhi; - uint8_t hwId; // CH=0, EL=1, NH=2, GAMMA=3, MU=4 - int16_t hwVtxEta; // propagate back to Vtx for charged particles (if useful?) - int16_t hwVtxPhi; - uint16_t hwFlags; - CaloCluster cluster; - PropagatedTrack track; - bool chargedPV; - uint16_t hwPuppiWeight; - uint16_t hwStatus; // for debugging - const l1t::Muon *muonsrc; - const l1t::PFCandidate *src; - - // sorting - bool operator<(const PFParticle &other) const { return hwPt > other.hwPt; } - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - static constexpr float PUPPI_SCALE = 100; - - float floatPt() const { return float(hwPt) / CaloCluster::PT_SCALE; } - float floatEta() const { return float(hwEta) / CaloCluster::ETAPHI_SCALE; } - float floatPhi() const { return float(hwPhi) / CaloCluster::ETAPHI_SCALE; } - float floatVtxEta() const { - return (track.hwPt > 0 ? track.floatVtxEta() : float(hwVtxEta) / CaloCluster::ETAPHI_SCALE); - } - float floatVtxPhi() const { - return (track.hwPt > 0 ? track.floatVtxPhi() : float(hwVtxPhi) / CaloCluster::ETAPHI_SCALE); - } - float floatDZ() const { return float(track.hwZ0) / InputTrack::Z0_SCALE; } - float floatPuppiW() const { return float(hwPuppiWeight) / PUPPI_SCALE; } - int intCharge() const { return (track.hwPt > 0 ? track.intCharge() : 0); } - void setPuppiW(float w) { hwPuppiWeight = std::round(w * PUPPI_SCALE); } - void setFloatPt(float pt) { hwPt = round(pt * CaloCluster::PT_SCALE); } -#endif - }; - - struct EGParticle { - int16_t hwPt; - int16_t hwEta; // at calo face - int16_t hwPhi; - uint16_t hwQual; - - // FIXME: an index would also do... - CaloCluster cluster; - - // sorting - bool operator<(const EGParticle &other) const { return hwPt > other.hwPt; } - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - void setFloatPt(float pt) { hwPt = round(pt * CaloCluster::PT_SCALE); } - float floatPt() const { return float(hwPt) / CaloCluster::PT_SCALE; } - float floatEta() const { return float(hwEta) / CaloCluster::ETAPHI_SCALE; } - float floatPhi() const { return float(hwPhi) / CaloCluster::ETAPHI_SCALE; } -#endif - }; - - struct EGIso { - // FIXME: eventually only one iso will be saved - uint16_t hwIso; - uint16_t hwPFIso; - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - static constexpr float ISO_SCALE = 100; - void setIso(float iso, uint16_t &hwIso) { hwIso = round(iso * EGIso::ISO_SCALE); } - void setIso(float iso) { setIso(iso, hwIso); } - void setPFIso(float iso) { setIso(iso, hwPFIso); } - - float getFloatIso(uint16_t hwIso) const { return float(hwIso) / EGIso::ISO_SCALE; } - float floatIso() const { return getFloatIso(hwIso); } - float floatPFIso() const { return getFloatIso(hwPFIso); } -#endif - }; - - struct EGIsoPV : public EGIso { - uint16_t hwIsoPV; - uint16_t hwPFIsoPV; - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - void setIsoPV(float iso) { setIso(iso, hwIsoPV); } - void setPFIsoPV(float iso) { setIso(iso, hwPFIsoPV); } - - float floatIsoPV() const { return getFloatIso(hwIsoPV); } - float floatPFIsoPV() const { return getFloatIso(hwPFIsoPV); } -#endif - }; - - struct EGIsoEleParticle : public EGParticle, public EGIso { - // track parameters for electrons - int16_t hwVtxEta; - int16_t hwVtxPhi; - int16_t hwZ0; - bool hwCharge; - PropagatedTrack track; - -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - - float floatVtxEta() const { return float(hwVtxEta) / InputTrack::VTX_ETA_SCALE; } - float floatVtxPhi() const { return float(hwVtxPhi) / InputTrack::VTX_PHI_SCALE; } - float floatDZ() const { return float(track.hwZ0) / InputTrack::Z0_SCALE; } - int intCharge() const { return hwCharge ? +1 : -1; } - -#endif - }; - - struct EGIsoParticle : public EGParticle, public EGIsoPV { -#ifdef L1Trigger_Phase2L1ParticleFlow_DiscretePFInputs_MORE - // NOTE: this is needed because of CMSSW requirements - // i.e. we need to put the EG object and the TkEm and TkEle ones at the same time to have a valid ref - int ele_idx; -#endif - }; - - struct InputRegion { - float etaCenter, etaMin, etaMax, phiCenter, phiHalfWidth; - float etaExtra, phiExtra; - std::vector calo; - std::vector emcalo; - std::vector track; - std::vector muon; - - InputRegion() - : etaCenter(), - etaMin(), - etaMax(), - phiCenter(), - phiHalfWidth(), - etaExtra(), - phiExtra(), - calo(), - emcalo(), - track(), - muon() {} - InputRegion( - float etacenter, float etamin, float etamax, float phicenter, float phihalfwidth, float etaextra, float phiextra) - : etaCenter(etacenter), - etaMin(etamin), - etaMax(etamax), - phiCenter(phicenter), - phiHalfWidth(phihalfwidth), - etaExtra(etaextra), - phiExtra(phiextra), - calo(), - emcalo(), - track(), - muon() {} - }; - -} // namespace l1tpf_impl -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/L1TCorrelatorLayer1PatternFileWriter.h b/L1Trigger/Phase2L1ParticleFlow/interface/L1TCorrelatorLayer1PatternFileWriter.h new file mode 100644 index 0000000000000..15004caa2f803 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/L1TCorrelatorLayer1PatternFileWriter.h @@ -0,0 +1,76 @@ +#ifndef L1Trigger_Phase2L1ParticleFlow_L1TCorrelatorLayer1PatternFileWriter_h +#define L1Trigger_Phase2L1ParticleFlow_L1TCorrelatorLayer1PatternFileWriter_h + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "L1Trigger/DemonstratorTools/interface/BoardDataWriter.h" +#include "L1Trigger/DemonstratorTools/interface/utilities.h" + +#include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" + +class L1TCorrelatorLayer1PatternFileWriter { +public: + L1TCorrelatorLayer1PatternFileWriter(const edm::ParameterSet& iConfig, const l1ct::Event& eventTemplate); + ~L1TCorrelatorLayer1PatternFileWriter(); + + void write(const l1ct::Event& event); + void flush(); + +private: + enum class Partition { Barrel, HGCal, HGCalNoTk, HF }; + + Partition partition_; + const unsigned int tmuxFactor_ = 6; // not really configurable in current architecture + bool writeInputs_, writeOutputs_; + std::map> channelIdsInput_, channelIdsOutput_; + std::map channelSpecsInput_, channelSpecsOutput_; + + const unsigned int tfTimeslices_ = 3, tfLinksFactor_ = 1; // not really configurable in current architecture + const unsigned int hgcTimeslices_ = 3, hgcLinksFactor_ = 4; // not really configurable in current architecture + const unsigned int gctTimeslices_ = 1, gctSectors_ = 3; // not really configurable in current architecture + const unsigned int gctLinksEcal_ = 1, gctLinksHad_ = 2; // could be made configurable later + const unsigned int gmtTimeslices_ = 3, gmtLinksFactor_ = 1; // not really configurable in current architecture + const unsigned int gttTimeslices_ = 1, gttLinksFactor_ = 1; // not really configurable in current architecture + uint32_t gmtNumberOfMuons_; + uint32_t gttNumberOfPVs_; + uint32_t gttLatency_; + + std::vector outputRegions_, outputLinksPuppi_; + unsigned int nPuppiFramesPerRegion_; + int32_t outputBoard_, outputLinkEgamma_; + uint32_t nEgammaObjectsOut_; + + // Common stuff related to the format + uint32_t nInputFramesPerBX_, nOutputFramesPerBX_; + const std::string fileFormat_; + + // final helper + const uint32_t eventsPerFile_; + uint32_t eventIndex_; + std::unique_ptr inputFileWriter_, outputFileWriter_; + + static Partition parsePartition(const std::string& partition); + + void configTimeSlices(const edm::ParameterSet& iConfig, + const std::string& prefix, + unsigned int nSectors, + unsigned int nTimeSlices, + unsigned int linksFactor); + void configSectors(const edm::ParameterSet& iConfig, + const std::string& prefix, + unsigned int nSectors, + unsigned int linksFactor); + void configLinks(const edm::ParameterSet& iConfig, + const std::string& prefix, + unsigned int linksFactor, + unsigned int offset); + + void writeTF(const l1ct::Event& event, l1t::demo::EventData& out); + void writeBarrelGCT(const l1ct::Event& event, l1t::demo::EventData& out); + void writeHGC(const l1ct::Event& event, l1t::demo::EventData& out); + void writeGMT(const l1ct::Event& event, l1t::demo::EventData& out); + void writeGTT(const l1ct::Event& event, l1t::demo::EventData& out); + void writePuppi(const l1ct::Event& event, l1t::demo::EventData& out); + void writeEgamma(const l1ct::Event& event, l1t::demo::EventData& out); +}; + +#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo2HGC.h b/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo2HGC.h deleted file mode 100644 index 92c9432344a0f..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo2HGC.h +++ /dev/null @@ -1,76 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_PFAlgo2HGC_h -#define L1Trigger_Phase2L1ParticleFlow_PFAlgo2HGC_h - -#include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgoBase.h" - -namespace l1tpf_impl { - class PFAlgo2HGC : public PFAlgoBase { - public: - PFAlgo2HGC(const edm::ParameterSet &); - void runPF(Region &r) const override; - - protected: - float drMatchMu_; - enum class MuMatchMode { BoxBestByPtRatio, DrBestByPtRatio, DrBestByPtDiff } muMatchMode_; - float drMatch_, ptMatchLow_, ptMatchHigh_, maxInvisiblePt_; - bool useTrackCaloSigma_, rescaleUnmatchedTrack_, caloTrkWeightedAverage_; - enum class TkCaloLinkMetric { BestByDR = 0, BestByDRPt = 1, BestByDR2Pt2 = 2 }; - TkCaloLinkMetric tkCaloLinkMetric_; - bool caloReLinkStep_; - float caloReLinkDr_, caloReLinkThreshold_; - bool rescaleTracks_, sumTkCaloErr2_, ecalPriority_, trackEmUseAlsoTrackSigma_, emCaloUseAlsoCaloSigma_; - float tightTrackMaxInvisiblePt_; - enum GoodTrackStatus { GoodTK_Calo_TkPt = 0, GoodTK_Calo_TkCaloPt = 1, GoodTk_Calo_CaloPt = 2, GoodTK_NoCalo = 3 }; - enum BadTrackStatus { BadTK_NoCalo = 1 }; - bool sortInputs_; - - /// do muon track linking (also sets track.muonLink) - void link_tk2mu(Region &r, std::vector &tk2mu, std::vector &mu2tk) const; - - /// track to calo matching - // tk2calo[itk] = icalo or -1 - void link_tk2calo(Region &r, std::vector &tk2calo) const; - - /// for each calo, compute the sum of the track pt - void sum_tk2calo(Region &r, - const std::vector &tk2calo, - std::vector &calo2ntk, - std::vector &calo2sumtkpt, - std::vector &calo2sumtkpterr) const; - - /// promote unlinked low pt tracks to hadrons - void unlinkedtk_algo(Region &r, const std::vector &tk2calo) const; - - /// try to recover split hadron showers (v1.0): - // take hadrons that are not track matched, close by a hadron which has an excess of track pt vs calo pt - // add this pt to the calo pt of the other cluster - // off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events - void calo_relink(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr) const; - - /// process matched calo clusters, compare energy to sum track pt, compute track rescaling factor if needed - // alpha[icalo] = x < 1 if all tracks linked to icalo must have their pt rescaled by x - void linkedcalo_algo(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr, - std::vector &calo2alpha) const; - - /// process matched tracks, if necessary rescale or average - void linkedtk_algo(Region &r, - const std::vector &tk2calo, - const std::vector &calo2ntk, - const std::vector &calo2alpha) const; - - /// process unmatched calo clusters - void unlinkedcalo_algo(Region &r) const; - - /// save muons in output list - void save_muons(Region &r, const std::vector &tk2mu) const; - }; - -} // namespace l1tpf_impl - -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo3.h b/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo3.h deleted file mode 100644 index 931fe0bc7be83..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo3.h +++ /dev/null @@ -1,109 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_PFAlgo3_h -#define L1Trigger_Phase2L1ParticleFlow_PFAlgo3_h - -#include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgoBase.h" - -namespace l1tpf_impl { - class PFAlgo3 : public PFAlgoBase { - public: - PFAlgo3(const edm::ParameterSet &); - void runPF(Region &r) const override; - - protected: - float drMatchMu_; - enum class MuMatchMode { BoxBestByPtRatio, DrBestByPtRatio, DrBestByPtDiff } muMatchMode_; - float drMatch_, ptMatchLow_, ptMatchHigh_, maxInvisiblePt_; - bool useTrackCaloSigma_, rescaleUnmatchedTrack_, caloTrkWeightedAverage_; - enum class TkCaloLinkMetric { BestByDR = 0, BestByDRPt = 1, BestByDR2Pt2 = 2 }; - float drMatchEm_, ptMinFracMatchEm_, drMatchEmHad_; - float emHadSubtractionPtSlope_; - TkCaloLinkMetric tkCaloLinkMetric_; - bool caloReLinkStep_; - float caloReLinkDr_, caloReLinkThreshold_; - bool rescaleTracks_, sumTkCaloErr2_, ecalPriority_, trackEmUseAlsoTrackSigma_, trackEmMayUseCaloMomenta_, - emCaloUseAlsoCaloSigma_; - float tightTrackMaxInvisiblePt_; - enum GoodTrackStatus { GoodTK_Calo_TkPt = 0, GoodTK_Calo_TkCaloPt = 1, GoodTk_Calo_CaloPt = 2, GoodTK_NoCalo = 3 }; - enum BadTrackStatus { BadTK_NoCalo = 1 }; - bool sortInputs_; - - /// do muon track linking (also sets track.muonLink) - void link_tk2mu(Region &r, std::vector &tk2mu, std::vector &mu2tk) const; - - /// match all tracks to the closest EM cluster - // tk2em[itrack] = iem, or -1 if unmatched - void link_tk2em(Region &r, std::vector &tk2em) const; - - /// match all em to the closest had (can happen in parallel to the above) - // em2calo[iem] = icalo or -1 - void link_em2calo(Region &r, std::vector &em2calo) const; - - /// for each EM cluster, count and add up the pt of all the corresponding tracks (skipping muons) - void sum_tk2em(Region &r, - const std::vector &tk2em, - std::vector &em2ntk, - std::vector &em2sumtkpt, - std::vector &em2sumtkpterr) const; - - /// process ecal clusters after linking - void emcalo_algo(Region &r, - const std::vector &em2ntk, - const std::vector &em2sumtkpt, - const std::vector &em2sumtkpterr) const; - - /// promote all flagged tracks to electrons - void emtk_algo(Region &r, - const std::vector &tk2em, - const std::vector &em2ntk, - const std::vector &em2sumtkpterr) const; - - /// subtract EM component from Calo clusters for all photons and electrons (within tracker coverage) - void sub_em2calo(Region &r, const std::vector &em2calo) const; - - /// track to calo matching - // tk2calo[itk] = icalo or -1 - void link_tk2calo(Region &r, std::vector &tk2calo) const; - - /// for each calo, compute the sum of the track pt - void sum_tk2calo(Region &r, - const std::vector &tk2calo, - std::vector &calo2ntk, - std::vector &calo2sumtkpt, - std::vector &calo2sumtkpterr) const; - - /// promote unlinked low pt tracks to hadrons - void unlinkedtk_algo(Region &r, const std::vector &tk2calo) const; - - /// try to recover split hadron showers (v1.0): - // take hadrons that are not track matched, close by a hadron which has an excess of track pt vs calo pt - // add this pt to the calo pt of the other cluster - // off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events - void calo_relink(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr) const; - - /// process matched calo clusters, compare energy to sum track pt, compute track rescaling factor if needed - // alpha[icalo] = x < 1 if all tracks linked to icalo must have their pt rescaled by x - void linkedcalo_algo(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr, - std::vector &calo2alpha) const; - - /// process matched tracks, if necessary rescale or average - void linkedtk_algo(Region &r, - const std::vector &tk2calo, - const std::vector &calo2ntk, - const std::vector &calo2alpha) const; - - /// process unmatched calo clusters - void unlinkedcalo_algo(Region &r) const; - - /// save muons in output list - void save_muons(Region &r, const std::vector &tk2mu) const; - }; - -} // namespace l1tpf_impl - -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgoBase.h b/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgoBase.h deleted file mode 100644 index 188d1746bed5f..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/PFAlgoBase.h +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_PFAlgoBase_h -#define L1Trigger_Phase2L1ParticleFlow_PFAlgoBase_h - -#include - -#include "L1Trigger/Phase2L1ParticleFlow/interface/Region.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -namespace l1tpf_impl { - - class PFAlgoBase { - public: - PFAlgoBase(const edm::ParameterSet &); - virtual ~PFAlgoBase(); - virtual void runPF(Region &r) const = 0; - - protected: - int debug_; - void initRegion(Region &r, bool doSort = true) const; - PFParticle &addTrackToPF(Region &r, const PropagatedTrack &tk) const { return addTrackToPF(r.pf, tk); } - PFParticle &addCaloToPF(Region &r, const CaloCluster &calo) const { return addCaloToPF(r.pf, calo); } - PFParticle &addTrackToPF(std::vector &pfs, const PropagatedTrack &tk) const; - PFParticle &addCaloToPF(std::vector &pfs, const CaloCluster &calo) const; - }; - -} // namespace l1tpf_impl - -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/PFTkEGAlgo.h b/L1Trigger/Phase2L1ParticleFlow/interface/PFTkEGAlgo.h deleted file mode 100644 index 9b4e613d36405..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/PFTkEGAlgo.h +++ /dev/null @@ -1,147 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_PFTkEGAlgo_h -#define L1Trigger_Phase2L1ParticleFlow_PFTkEGAlgo_h - -#include -#include - -#include "L1Trigger/Phase2L1ParticleFlow/interface/Region.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -namespace l1tpf_impl { - - class PFTkEGAlgo { - public: - PFTkEGAlgo(const edm::ParameterSet &); - virtual ~PFTkEGAlgo(); - void runTkEG(Region &r) const; - void runTkIso(Region &r, const float z0) const; - void runPFIso(Region &r, const float z0) const; - - bool writeEgSta() const { return writeEgSta_; } - - protected: - struct IsoParameters { - IsoParameters(const edm::ParameterSet &); - float tkQualityPtMin; - float dZ; - float dRMin; - float dRMax; - float tkQualityChi2Max; - float dRMin2; - float dRMax2; - ; - }; - - int debug_; - bool doBremRecovery_; - bool doTkIsolation_; - bool filterHwQuality_; - int caloHwQual_; - float dEtaMaxBrem_; - float dPhiMaxBrem_; - std::vector absEtaBoundaries_; - std::vector dEtaValues_; - std::vector dPhiValues_; - float caloEtMin_; - float trkQualityPtMin_; - float trkQualityChi2_; - bool writeEgSta_; - IsoParameters tkIsoParametersTkEm_; - IsoParameters tkIsoParametersTkEle_; - IsoParameters pfIsoParametersTkEm_; - IsoParameters pfIsoParametersTkEle_; - - void initRegion(Region &r) const; - void link_emCalo2emCalo(const Region &r, std::vector &emCalo2emCalo) const; - void link_emCalo2tk(const Region &r, std::vector &emCalo2tk) const; - - template - void compute_isolation_tkEm( - Region &r, const std::vector &objects, const IsoParameters ¶ms, const float z0, bool isPF) const { - for (int ic = 0, nc = r.egphotons.size(); ic < nc; ++ic) { - auto &egphoton = r.egphotons[ic]; - - float sumPt = 0.; - float sumPtPV = 0.; - - for (int itk = 0, ntk = objects.size(); itk < ntk; ++itk) { - const auto &tk = objects[itk]; - - if (tk.floatPt() < params.tkQualityPtMin) - continue; - - // FIXME: we compare Tk at vertex against the calo variable....shall we correct for the PV position ? - float d_phi = deltaPhi(tk.floatVtxPhi(), egphoton.floatPhi()); - float d_eta = tk.floatVtxEta() - egphoton.floatEta(); - float dR2 = d_phi * d_phi + d_eta * d_eta; - - if (dR2 > params.dRMin2 && dR2 < params.dRMax2) { - sumPt += tk.floatPt(); - // PF neutrals are not constrained by PV (since their Z0 is 0 by design) - if (tk.intCharge() == 0 || std::abs(tk.floatDZ() - z0) < params.dZ) - sumPtPV += tk.floatPt(); - } - } - if (isPF) { - egphoton.setPFIso(sumPt / egphoton.floatPt()); - egphoton.setPFIsoPV(sumPtPV / egphoton.floatPt()); - } else { - egphoton.setIso(sumPt / egphoton.floatPt()); - egphoton.setIsoPV(sumPtPV / egphoton.floatPt()); - } - } - } - - template - void compute_isolation_tkEle( - Region &r, const std::vector &objects, const IsoParameters ¶ms, const float z0, bool isPF) const { - for (int ic = 0, nc = r.egeles.size(); ic < nc; ++ic) { - auto &egele = r.egeles[ic]; - - float sumPt = 0.; - - for (int itk = 0, ntk = objects.size(); itk < ntk; ++itk) { - const auto &tk = objects[itk]; - - if (tk.floatPt() < params.tkQualityPtMin) - continue; - - // we check the DZ only for charged PFParticles for which Z0 is assigned to (0,0,0) - if (tk.intCharge() != 0 && std::abs(tk.floatDZ() - egele.floatDZ()) > params.dZ) - continue; - - float d_phi = deltaPhi(tk.floatVtxPhi(), egele.floatVtxPhi()); - float d_eta = tk.floatVtxEta() - egele.floatVtxEta(); - float dR2 = d_phi * d_phi + d_eta * d_eta; - - if (dR2 > params.dRMin2 && dR2 < params.dRMax2) { - sumPt += tk.floatPt(); - } - } - if (isPF) { - egele.setPFIso(sumPt / egele.floatPt()); - } else { - egele.setIso(sumPt / egele.floatPt()); - } - } - } - - void eg_algo(Region &r, const std::vector &emCalo2emCalo, const std::vector &emCalo2tk) const; - - void addEgObjsToPF(Region &r, const int calo_idx, const int hwQual, const float ptCorr, const int tk_idx = -1) const; - - EGIsoParticle &addEGIsoToPF(std::vector &egobjs, - const CaloCluster &calo, - const int hwQual, - const float ptCorr) const; - - EGIsoEleParticle &addEGIsoEleToPF(std::vector &egobjs, - const CaloCluster &calo, - const PropagatedTrack &track, - const int hwQual, - const float ptCorr) const; - }; - -} // namespace l1tpf_impl - -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/Region.h b/L1Trigger/Phase2L1ParticleFlow/interface/Region.h deleted file mode 100644 index 4c01068b1360e..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/Region.h +++ /dev/null @@ -1,121 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_Region_h -#define L1Trigger_Phase2L1ParticleFlow_Region_h - -#include "L1Trigger/Phase2L1ParticleFlow/interface/DiscretePFInputs.h" -#include "DataFormats/Math/interface/deltaPhi.h" - -namespace l1tpf_impl { - - struct Region : public InputRegion { - std::vector pf; - std::vector puppi; - std::vector egeles; - std::vector egphotons; - - unsigned int caloOverflow, emcaloOverflow, trackOverflow, muonOverflow, pfOverflow, puppiOverflow; - - const bool relativeCoordinates; // whether the eta,phi in each region are global or relative to the region center - const unsigned int ncaloMax, nemcaloMax, ntrackMax, nmuonMax, npfMax, npuppiMax; - Region(float etamin, - float etamax, - float phicenter, - float phiwidth, - float etaextra, - float phiextra, - bool useRelativeCoordinates, - unsigned int ncalomax, - unsigned int nemcalomax, - unsigned int ntrackmax, - unsigned int nmuonmax, - unsigned int npfmax, - unsigned int npuppimax) - : InputRegion(0.5 * (etamin + etamax), etamin, etamax, phicenter, 0.5 * phiwidth, etaextra, phiextra), - pf(), - puppi(), - egeles(), - egphotons(), - caloOverflow(), - emcaloOverflow(), - trackOverflow(), - muonOverflow(), - pfOverflow(), - puppiOverflow(), - relativeCoordinates(useRelativeCoordinates), - ncaloMax(ncalomax), - nemcaloMax(nemcalomax), - ntrackMax(ntrackmax), - nmuonMax(nmuonmax), - npfMax(npfmax), - npuppiMax(npuppimax) {} - - enum InputType { calo_type = 0, emcalo_type = 1, track_type = 2, l1mu_type = 3, n_input_types = 4 }; - static const char* inputTypeName(int inputType); - - enum OutputType { - any_type = 0, - charged_type = 1, - neutral_type = 2, - electron_type = 3, - pfmuon_type = 4, - charged_hadron_type = 5, - neutral_hadron_type = 6, - photon_type = 7, - n_output_types = 8 - }; - static const char* outputTypeName(int outputType); - - unsigned int nInput(InputType type) const; - unsigned int nOutput(OutputType type, bool puppi, bool fiducial = true) const; - - // global coordinates - bool contains(float eta, float phi) const { - float dphi = deltaPhi(phiCenter, phi); - return (etaMin - etaExtra < eta && eta <= etaMax + etaExtra && -phiHalfWidth - phiExtra < dphi && - dphi <= phiHalfWidth + phiExtra); - } - // global coordinates - bool fiducial(float eta, float phi) const { - float dphi = deltaPhi(phiCenter, phi); - return (etaMin < eta && eta <= etaMax && -phiHalfWidth < dphi && dphi <= phiHalfWidth); - } - // possibly local coordinates - bool fiducialLocal(float localEta, float localPhi) const { - if (relativeCoordinates) { - float dphi = deltaPhi(0.f, localPhi); - return (etaMin < localEta + etaCenter && localEta + etaCenter <= etaMax && -phiHalfWidth < dphi && - dphi <= phiHalfWidth); - } - float dphi = deltaPhi(phiCenter, localPhi); - return (etaMin < localEta && localEta <= etaMax && -phiHalfWidth < dphi && dphi <= phiHalfWidth); - } - float regionAbsEta() const { return std::abs(etaCenter); } - float globalAbsEta(float localEta) const { return std::abs(relativeCoordinates ? localEta + etaCenter : localEta); } - float globalEta(float localEta) const { return relativeCoordinates ? localEta + etaCenter : localEta; } - float globalPhi(float localPhi) const { return relativeCoordinates ? localPhi + phiCenter : localPhi; } - float localEta(float globalEta) const { return relativeCoordinates ? globalEta - etaCenter : globalEta; } - float localPhi(float globalPhi) const { return relativeCoordinates ? deltaPhi(globalPhi, phiCenter) : globalPhi; } - - void zero() { - calo.clear(); - emcalo.clear(); - track.clear(); - muon.clear(); - pf.clear(); - puppi.clear(); - egeles.clear(); - egphotons.clear(); - caloOverflow = 0; - emcaloOverflow = 0; - trackOverflow = 0; - muonOverflow = 0; - pfOverflow = 0; - puppiOverflow = 0; - } - - void inputCrop(bool doSort); - void outputCrop(bool doSort); - }; - -} // namespace l1tpf_impl - -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/RegionMapper.h b/L1Trigger/Phase2L1ParticleFlow/interface/RegionMapper.h deleted file mode 100644 index 68455b9e13f0a..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/interface/RegionMapper.h +++ /dev/null @@ -1,66 +0,0 @@ -#ifndef L1Trigger_Phase2L1ParticleFlow_RegionMapper_h -#define L1Trigger_Phase2L1ParticleFlow_RegionMapper_h - -#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" -#include "L1Trigger/Phase2L1ParticleFlow/interface/Region.h" -#include "DataFormats/Math/interface/deltaPhi.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "DataFormats/L1TCorrelator/interface/TkMuon.h" -#include "DataFormats/L1TCorrelator/interface/TkMuonFwd.h" - -#include - -namespace edm { - class Event; -} - -namespace l1tpf_impl { - class RegionMapper { - // This does the input and filling of regions. - public: - RegionMapper(const edm::ParameterSet &); - - // add object, without tracking references - void addTrack(const l1t::PFTrack &t); - void addMuon(const l1t::Muon &t); - void addMuon(const l1t::TkMuon &t); - void addCalo(const l1t::PFCluster &t); - void addEmCalo(const l1t::PFCluster &t); - - // add object, tracking references - void addTrack(const l1t::PFTrack &t, l1t::PFTrackRef ref); - void addCalo(const l1t::PFCluster &t, l1t::PFClusterRef ref); - void addEmCalo(const l1t::PFCluster &t, l1t::PFClusterRef ref); - - void clear(); - std::vector ®ions() { return regions_; } - - std::unique_ptr fetch(bool puppi = true, float ptMin = 0.01) const; - std::unique_ptr fetchCalo(float ptMin = 0.01, bool emcalo = false) const; - std::unique_ptr fetchTracks(float ptMin = 0.01, bool fromPV = false) const; - - void putEgObjects(edm::Event &iEvent, - const bool writeEgSta, - const std::string &egLablel, - const std::string &tkEmLabel, - const std::string &tkEleLabel, - const float ptMin = 0.01) const; - - std::pair totAndMaxInput(/*Region::InputType*/ int type) const; - std::pair totAndMaxOutput(/*Region::OutputType*/ int type, bool puppi) const; - std::unique_ptr> vecInput(int type) const; - std::unique_ptr> vecOutput(int type, bool puppi) const; - - protected: - std::vector regions_; - bool useRelativeRegionalCoordinates_; // whether the eta,phi in each region are global or relative to the region center - enum class TrackAssoMode { atVertex, atCalo, any = 999 } trackRegionMode_; - - // these are used to link items back - std::unordered_map clusterRefMap_; - std::unordered_map trackRefMap_; - }; - -} // namespace l1tpf_impl -#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egencoder_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egencoder_ref.h index c42a7c709f15b..010fa8801d692 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egencoder_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egencoder_ref.h @@ -3,15 +3,8 @@ #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" #include "DataFormats/L1TParticleFlow/interface/egamma.h" - -#ifdef CMSSW_GIT_HASH #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" -#else -#include "../../utils/dbgPrintf.h" - -#endif - namespace edm { class ParameterSet; } @@ -19,7 +12,6 @@ namespace edm { namespace l1ct { struct L2EgEncoderEmulator { - public: L2EgEncoderEmulator(unsigned int nTKELE_OUT, unsigned int nTKPHO_OUT) : nTkEleOut_(nTKELE_OUT), nTkPhoOut_(nTKPHO_OUT), nEncodedWords_(nTKELE_OUT * 1.5 + nTKPHO_OUT * 1.5) { assert(nTkEleOut_ % 2 == 0); diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egsorter_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egsorter_ref.h index 9804008b137bb..c5b75b37384e6 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egsorter_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/l2egsorter_ref.h @@ -4,18 +4,10 @@ #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" #include "DataFormats/L1TParticleFlow/interface/egamma.h" #include "DataFormats/L1TParticleFlow/interface/pf.h" - -#ifdef CMSSW_GIT_HASH #include "L1Trigger/Phase2L1ParticleFlow/interface/common/bitonic_hybrid_sort_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" -#else -#include "L1Trigger/Phase2L1ParticleFlow/interface/common/bitonic_hybrid_sort_ref.h" -#include "../../utils/dbgPrintf.h" - -#endif #include -#include namespace edm { class ParameterSet; @@ -93,7 +85,7 @@ namespace l1ct { void print_objects(const std::vector &objs, const std::string &label) const { for (unsigned int i = 0; i < objs.size(); ++i) { dbgCout() << label << " [" << i << "] pt: " << objs[i].hwPt << " eta: " << objs[i].hwEta - << " phi: " << objs[i].hwPhi << " qual: " << std::bitset<4>(objs[i].hwQual) << std::endl; + << " phi: " << objs[i].hwPhi << " qual: " << objs[i].hwQual.to_string(2) << std::endl; } } diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index b4a4caaf009a8..b75be595ee9df 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -21,6 +21,7 @@ namespace l1ct { bool doBremRecovery; bool writeBeforeBremRecovery; int caloHwQual; + bool doEndcapHwQual; float emClusterPtMin; // GeV float dEtaMaxBrem; float dPhiMaxBrem; @@ -63,6 +64,7 @@ namespace l1ct { bool doBremRecovery, bool writeBeforeBremRecovery = false, int caloHwQual = 4, + bool doEndcapHwQual = false, float emClusterPtMin = 2., float dEtaMaxBrem = 0.02, float dPhiMaxBrem = 0.1, @@ -78,7 +80,8 @@ namespace l1ct { bool doTkIso = true, bool doPfIso = false, EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, - EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV) + EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, + int debug = 0) : nTRACK(nTrack), nTRACK_EGIN(nTrack_in), @@ -88,6 +91,7 @@ namespace l1ct { doBremRecovery(doBremRecovery), writeBeforeBremRecovery(writeBeforeBremRecovery), caloHwQual(caloHwQual), + doEndcapHwQual(doEndcapHwQual), emClusterPtMin(emClusterPtMin), dEtaMaxBrem(dEtaMaxBrem), dPhiMaxBrem(dPhiMaxBrem), @@ -103,7 +107,8 @@ namespace l1ct { doTkIso(doTkIso), doPfIso(doPfIso), hwIsoTypeTkEle(hwIsoTypeTkEle), - hwIsoTypeTkEm(hwIsoTypeTkEm) {} + hwIsoTypeTkEm(hwIsoTypeTkEm), + debug(debug) {} }; class PFTkEGAlgoEmulator { @@ -157,26 +162,26 @@ namespace l1ct { const std::vector &emcalo, const std::vector &track, const int calo_idx, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr, const int tk_idx, const std::vector &components = {}) const; EGObjEmu &addEGStaToPF(std::vector &egobjs, const EmCaloObjEmu &calo, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr, const std::vector &components) const; EGIsoObjEmu &addEGIsoToPF(std::vector &egobjs, const EmCaloObjEmu &calo, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr) const; EGIsoEleObjEmu &addEGIsoEleToPF(std::vector &egobjs, const EmCaloObjEmu &calo, const TkObjEmu &track, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr) const; // FIXME: reimplemented from PFAlgoEmulatorBase diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegsorter_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegsorter_ref.h index 827d756f90c7e..e0f59a04b2f31 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegsorter_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegsorter_ref.h @@ -5,12 +5,10 @@ #include #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" +#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" #ifdef CMSSW_GIT_HASH #include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" -#else -#include "../../../utils/dbgPrintf.h" #endif namespace edm { diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFHtEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFHtEmulator.h index 6a17677de5027..68e575e491cfa 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFHtEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1PFHtEmulator.h @@ -3,13 +3,10 @@ #include "DataFormats/L1TParticleFlow/interface/jets.h" #include "DataFormats/L1TParticleFlow/interface/sums.h" - -#ifdef CMSSW_GIT_HASH -#include "./L1SeedConePFJetEmulator.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" -#else -#include "../../seededcone/ref/L1SeedConePFJetEmulator.h" -#include "../../../utils/dbgPrintf.h" +#include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h" + +#ifndef CMSSW_GIT_HASH #include "hls_math.h" #endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h index 437baa940182f..00be9c8ee7d69 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h @@ -9,11 +9,7 @@ #include #include -#ifdef CMSSW_GIT_HASH #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" -#else -#include "../../../utils/dbgPrintf.h" -#endif class L1SCJetEmu { public: diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_common_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_common_ref.h index dd04494e30dd5..eacb4ee7ff178 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_common_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_common_ref.h @@ -2,12 +2,7 @@ #define PFALGO_COMMON_REF_H #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" - -#ifdef CMSSW_GIT_HASH #include "pfalgo_types.h" -#else -#include "../firmware/pfalgo_types.h" -#endif #include #include diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/regionizer/multififo_regionizer_elements_ref.icc b/L1Trigger/Phase2L1ParticleFlow/interface/regionizer/multififo_regionizer_elements_ref.icc index 488f704317687..66b0538ac38fe 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/regionizer/multififo_regionizer_elements_ref.icc +++ b/L1Trigger/Phase2L1ParticleFlow/interface/regionizer/multififo_regionizer_elements_ref.icc @@ -53,10 +53,12 @@ void l1ct::multififo_regionizer::RegionBuffer::initFifos(unsigned int nfifos) for (auto& t : queues_.back().second) t.clear(); } - if (!(nfifos == 1 || nfifos == 2 || nfifos == 4 || nfifos == 6 || nfifos == 8 || nfifos == 12)) { + bool isGood = + (nfifos == 1 || nfifos == 2 || nfifos == 3 || nfifos == 4 || nfifos == 6 || nfifos == 8 || nfifos == 12); + if (!isGood) { dbgCerr() << "Error, created regionizer for nfifos == " << nfifos << ", not supported." << std::endl; } - assert(nfifos == 1 || nfifos == 2 || nfifos == 4 || nfifos == 6 || nfifos == 8 || nfifos == 12); + assert(isGood); } template @@ -79,7 +81,7 @@ void l1ct::multififo_regionizer::RegionBuffer::maybe_push(int fifo, const T& template T l1ct::multififo_regionizer::RegionBuffer::pop() { - if (nfifos_ <= 2) // probably works for 3 as well, but not tested + if (nfifos_ <= 3) return pop_next_trivial_(); assert(!queues_.empty()); for (unsigned int istep = 0, nsteps = queues_.size(); istep < nsteps; ++istep) { diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index 27c160c85b633..3eee113d3b15a 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -31,6 +31,7 @@ #include "L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo_common_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegsorter_ref.h" +#include "L1Trigger/Phase2L1ParticleFlow/interface/L1TCorrelatorLayer1PatternFileWriter.h" #include "DataFormats/L1TCorrelator/interface/TkElectron.h" #include "DataFormats/L1TCorrelator/interface/TkElectronFwd.h" @@ -76,6 +77,8 @@ class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { const std::string regionDumpName_; bool writeRawHgcalCluster_; std::fstream fRegionDump_; + const std::vector patternWriterConfigs_; + std::vector> patternWriters_; // region of interest debugging float debugEta_, debugPhi_, debugR_; @@ -87,6 +90,7 @@ class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { // main methods void beginStream(edm::StreamID) override; + void endStream() override; void produce(edm::Event &, const edm::EventSetup &) override; void addUInt(unsigned int value, std::string iLabel, edm::Event &iEvent); @@ -156,6 +160,8 @@ L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet l1tkegsorter_(nullptr), regionDumpName_(iConfig.getUntrackedParameter("dumpFileName", "")), writeRawHgcalCluster_(iConfig.getUntrackedParameter("writeRawHgcalCluster", false)), + patternWriterConfigs_(iConfig.getUntrackedParameter>( + "patternWriters", std::vector())), debugEta_(iConfig.getUntrackedParameter("debugEta", 0)), debugPhi_(iConfig.getUntrackedParameter("debugPhi", 0)), debugR_(iConfig.getUntrackedParameter("debugR", -1)) { @@ -275,6 +281,22 @@ void L1TCorrelatorLayer1Producer::beginStream(edm::StreamID id) { << "Job running with multiple streams, but dump file will have only events on stream zero."; } } + if (!patternWriterConfigs_.empty()) { + if (id == 0) { + for (const auto &pset : patternWriterConfigs_) { + patternWriters_.emplace_back(std::make_unique(pset, event_)); + } + } else { + edm::LogWarning("L1TCorrelatorLayer1Producer") + << "Job running with multiple streams, but pattern files will be written only on stream zero."; + } + } +} + +void L1TCorrelatorLayer1Producer::endStream() { + for (auto &writer : patternWriters_) { + writer->flush(); + } } // ------------ method called to produce the data ------------ @@ -457,6 +479,9 @@ void L1TCorrelatorLayer1Producer::produce(edm::Event &iEvent, const edm::EventSe if (fRegionDump_.is_open()) { event_.write(fRegionDump_); } + for (auto &writer : patternWriters_) { + writer->write(event_); + } // finally clear the regions event_.clear(); diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc index 4d7248b4712b6..75265d23765a0 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc @@ -49,7 +49,7 @@ class L1TCtL2EgProducer : public edm::global::EDProducer<> { BXVector>> oldRefs; std::map>, edm::Ref>> old2newRefMap; - std::vector &, edm::Ptr>> origRefAndPtr; + std::vector, edm::Ptr>> origRefAndPtr; }; void convertToEmu(const l1t::TkElectron &tkele, RefRemapper &refRemapper, l1ct::OutputBoard &boarOut) const; @@ -332,8 +332,9 @@ void L1TCtL2EgProducer::convertToEmu(const l1t::TkElectron &tkele, } refRemapper.origRefAndPtr.push_back(std::make_pair(refEg, tkele.trkPtr())); emu.sta_idx = refRemapper.origRefAndPtr.size() - 1; - emu.setHwIso(EGIsoEleObjEmu::IsoType::TkIso, l1ct::Scales::makeIso(tkele.trkIsol())); - emu.setHwIso(EGIsoEleObjEmu::IsoType::PfIso, l1ct::Scales::makeIso(tkele.pfIsol())); + // NOTE: The emulator and FW data-format stores absolute iso while the CMSSW object stores relative iso + emu.setHwIso(EGIsoEleObjEmu::IsoType::TkIso, l1ct::Scales::makeIso(tkele.trkIsol() * tkele.pt())); + emu.setHwIso(EGIsoEleObjEmu::IsoType::PfIso, l1ct::Scales::makeIso(tkele.pfIsol() * tkele.pt())); // std::cout << "[convertToEmu] TkEle pt: " << emu.hwPt << " eta: " << emu.hwEta << " phi: " << emu.hwPhi << " staidx: " << emu.sta_idx << std::endl; boarOut.egelectron.push_back(emu); @@ -352,44 +353,48 @@ void L1TCtL2EgProducer::convertToEmu(const l1t::TkEm &tkem, } refRemapper.origRefAndPtr.push_back(std::make_pair(refEg, edm::Ptr(nullptr, 0))); emu.sta_idx = refRemapper.origRefAndPtr.size() - 1; - emu.setHwIso(EGIsoObjEmu::IsoType::TkIso, l1ct::Scales::makeIso(tkem.trkIsol())); - emu.setHwIso(EGIsoObjEmu::IsoType::PfIso, l1ct::Scales::makeIso(tkem.pfIsol())); - emu.setHwIso(EGIsoObjEmu::IsoType::TkIsoPV, l1ct::Scales::makeIso(tkem.trkIsolPV())); - emu.setHwIso(EGIsoObjEmu::IsoType::PfIsoPV, l1ct::Scales::makeIso(tkem.pfIsolPV())); + // NOTE: The emulator and FW data-format stores absolute iso while the CMSSW object stores relative iso + emu.setHwIso(EGIsoObjEmu::IsoType::TkIso, l1ct::Scales::makeIso(tkem.trkIsol() * tkem.pt())); + emu.setHwIso(EGIsoObjEmu::IsoType::PfIso, l1ct::Scales::makeIso(tkem.pfIsol() * tkem.pt())); + emu.setHwIso(EGIsoObjEmu::IsoType::TkIsoPV, l1ct::Scales::makeIso(tkem.trkIsolPV() * tkem.pt())); + emu.setHwIso(EGIsoObjEmu::IsoType::PfIsoPV, l1ct::Scales::makeIso(tkem.pfIsolPV() * tkem.pt())); // std::cout << "[convertToEmu] TkEM pt: " << emu.hwPt << " eta: " << emu.hwEta << " phi: " << emu.hwPhi << " staidx: " << emu.sta_idx << std::endl; boarOut.egphoton.push_back(emu); } l1t::TkEm L1TCtL2EgProducer::convertFromEmu(const l1ct::EGIsoObjEmu &egiso, const RefRemapper &refRemapper) const { // std::cout << "[convertFromEmu] TkEm pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi << " staidx: " << egiso.sta_idx << std::endl; - - reco::Candidate::PolarLorentzVector mom(egiso.floatPt(), egiso.floatEta(), egiso.floatPhi(), 0.); + // NOTE: the TkEM object is created with the accuracy as in GT object (not the Correlator internal one)! + const auto gteg = egiso.toGT(); + reco::Candidate::PolarLorentzVector mom( + l1gt::Scales::floatPt(gteg.v3.pt), l1gt::Scales::floatEta(gteg.v3.eta), l1gt::Scales::floatPhi(gteg.v3.phi), 0.); + // NOTE: The emulator and FW data-format stores absolute iso while the CMSSW object stores relative iso l1t::TkEm tkem(reco::Candidate::LorentzVector(mom), refRemapper.origRefAndPtr[egiso.sta_idx].first, egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::TkIso), egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::TkIsoPV)); - // FIXME: need to define a global quality (barrel+endcap) or add a bit to distibguish them? - tkem.setHwQual(egiso.hwQual); + tkem.setHwQual(gteg.quality); tkem.setPFIsol(egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::PfIso)); tkem.setPFIsolPV(egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::PfIsoPV)); - tkem.setEgBinaryWord(egiso.toGT().pack()); + tkem.setEgBinaryWord(gteg.pack()); return tkem; } l1t::TkElectron L1TCtL2EgProducer::convertFromEmu(const l1ct::EGIsoEleObjEmu &egele, const RefRemapper &refRemapper) const { // std::cout << "[convertFromEmu] TkEle pt: " << egele.hwPt << " eta: " << egele.hwEta << " phi: " << egele.hwPhi << " staidx: " << egele.sta_idx << std::endl; - - reco::Candidate::PolarLorentzVector mom(egele.floatPt(), egele.hwEta, egele.hwPhi, 0.); - + // NOTE: the TkElectron object is created with the accuracy as in GT object (not the Correlator internal one)! + const auto gteg = egele.toGT(); + reco::Candidate::PolarLorentzVector mom( + l1gt::Scales::floatPt(gteg.v3.pt), l1gt::Scales::floatEta(gteg.v3.eta), l1gt::Scales::floatPhi(gteg.v3.phi), 0.); + // NOTE: The emulator and FW data-format stores absolute iso while the CMSSW object stores relative iso l1t::TkElectron tkele(reco::Candidate::LorentzVector(mom), refRemapper.origRefAndPtr[egele.sta_idx].first, refRemapper.origRefAndPtr[egele.sta_idx].second, egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::TkIso)); - // FIXME: need to define a global quality (barrel+endcap)? - tkele.setHwQual(egele.hwQual); + tkele.setHwQual(gteg.quality); tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso)); - tkele.setEgBinaryWord(egele.toGT().pack()); + tkele.setEgBinaryWord(gteg.pack()); return tkele; } diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index 6c9faa22e5531..3d1621b655d56 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -10,6 +10,7 @@ writeBeforeBremRecovery=cms.bool(True), filterHwQuality=cms.bool(False), caloHwQual=cms.int32(4), + doEndcapHwQual=cms.bool(False), dEtaMaxBrem=cms.double(0.02), dPhiMaxBrem=cms.double(0.1), absEtaBoundaries=cms.vdouble(0.0, 0.9, 1.5), diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py index 355f9517fa8fd..9005fdf100820 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py @@ -112,7 +112,7 @@ caloSectors = cms.VPSet( cms.PSet( etaBoundaries = cms.vdouble(-1.5, 1.5), - phiSlices = cms.uint32(6), + phiSlices = cms.uint32(3), phiZero = cms.double(0), ) ), @@ -141,12 +141,12 @@ cms.PSet( etaBoundaries = cms.vdouble(-3.0, -1.5), phiSlices = cms.uint32(3), - phiZero = cms.double(math.pi/6) # L1 TrackFinder phi sector and HGCal sectors shifted by 30deg, + phiZero = cms.double(math.pi/2) # The edge of the 0th HGCal sectors is at 30 deg, the center at 30+120/2=90 = pi/2 ), cms.PSet( etaBoundaries = cms.vdouble(+1.5, +3.0), phiSlices = cms.uint32(3), - phiZero = cms.double(math.pi/6) # L1 TrackFinder phi sector and HGCal sectors shifted by 30deg, + phiZero = cms.double(math.pi/2) # As above ) ) @@ -255,6 +255,7 @@ nEMCALO_EGIN = 10, nEM_EGOUT = 5, doBremRecovery=True, + doEndcapHwQual=True, writeBeforeBremRecovery=False, writeEGSta=True), tkEgSorterParameters=tkEgSorterParameters.clone( @@ -340,6 +341,7 @@ nEMCALO_EGIN = 10, nEM_EGOUT = 5, doBremRecovery=True, + doEndcapHwQual=True, writeBeforeBremRecovery=False, writeEGSta=True), tkEgSorterParameters=tkEgSorterParameters.clone( diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_patternWriters_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_patternWriters_cff.py new file mode 100644 index 0000000000000..6a856e48db502 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_patternWriters_cff.py @@ -0,0 +1,165 @@ +import FWCore.ParameterSet.Config as cms + +eventsPerFile_ = 12 +gttLatency_ = 156+120 +gttNumberOfPVs_ = 10 + +##################################################################################################################### +## Barrel configurations: 54 regions, 6 puppi output links, only write out the layer 1 outputs for now +barrelWriterOutputOnly_ = cms.PSet( + partition = cms.string("Barrel"), + outputLinksPuppi = cms.vuint32(*range(6)), + outputLinkEgamma = cms.int32(6), + nEgammaObjectsOut = cms.uint32(16), + nOutputFramesPerBX = cms.uint32(9), + fileFormat = cms.string("EMP"), + maxLinesPerOutputFile = cms.uint32(1024), + eventsPerFile = cms.uint32(eventsPerFile_), +) +## Barrel (54) split in 3 eta slices +barrelWriterOutputOnlyEtaConfigs = [ + barrelWriterOutputOnly_.clone( + outputRegions = cms.vuint32(*[18*ie+i for i in range(18)]), + outputFileName = cms.string("l1BarrelEta%d-outputs-ideal" % (ie+1)), + outputBoard = cms.int32(-1), ## can't output e/gamma in eta split regions + outputLinkEgamma = cms.int32(-1), ## since the boards are defined in phi regions + ) for ie in range(3) +] +## Barrel (54) split in 3 phi slices +barrelWriterOutputOnlyPhiConfigs = [ + barrelWriterOutputOnly_.clone( + outputRegions = cms.vuint32(*[3*ip+9*ie+i for ie in range(6) for i in range(3) ]), + outputBoard = cms.int32(ip), + outputFileName = cms.string("l1BarrelPhi%d-outputs-ideal" % (ip+1)) + ) for ip in range(3) +] +## Barrel9 (27) split in phi eta slices +barrel9WriterOutputOnlyPhiConfigs = [ + barrelWriterOutputOnly_.clone( + outputRegions = cms.vuint32(*[3*ip+9*ie+i for ie in range(3) for i in range(3) ]), + outputBoard = cms.int32(ip), + outputFileName = cms.string("l1Barrel9Phi%d-outputs-ideal" % (ip+1)) + ) for ip in range(3) +] + +barrelWriterConfigs = barrelWriterOutputOnlyPhiConfigs # + barrelWriterOutputOnlyEtaConfigs +barrel9WriterConfigs = [] #barrel9WriterOutputOnlyPhiConfigs + + +##################################################################################################################### +## HGcal configuration: write out both inputs and outputs +hgcalWriterConfig_ = cms.PSet( + partition = cms.string("HGCal"), + outputRegions = cms.vuint32(*[i+9 for i in range(9)]), + outputBoard = cms.int32(1), + nEgammaObjectsOut = cms.uint32(16), + nInputFramesPerBX = cms.uint32(9), + nOutputFramesPerBX = cms.uint32(9), + fileFormat = cms.string("EMP"), + maxLinesPerInputFile = cms.uint32(1024), + maxLinesPerOutputFile = cms.uint32(1024), + eventsPerFile = cms.uint32(eventsPerFile_), + tfTimeSlices = cms.VPSet(*[cms.PSet(tfSectors = cms.VPSet()) for i in range(3)]), + hgcTimeSlices = cms.VPSet(*[cms.PSet(hgcSectors = cms.VPSet()) for i in range(3)]), + gmtTimeSlices = cms.VPSet(cms.PSet(),cms.PSet(),cms.PSet()), + gmtNumberOfMuons = cms.uint32(12), + gttLink = cms.int32(-1), + gttLatency = cms.uint32(gttLatency_), + gttNumberOfPVs = cms.uint32(gttNumberOfPVs_) +) +## Ideal configuration: 27 input links from tf, 36 from hgc, 3 from gmt, 1 from gtt, in this order; output 3 puppi + 1 e/gamma +hgcalPosIdealWriterConfig = hgcalWriterConfig_.clone() +for t in range(3): + hgcalPosIdealWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(-1)) for i in range(9) ] # neg + hgcalPosIdealWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(3*i+t)) for i in range(9) ] # pos + hgcalPosIdealWriterConfig.hgcTimeSlices[t].hgcSectors += [ cms.PSet(hgcLinks = cms.vint32(-1,-1,-1,-1)) for i in range(3) ] # neg + hgcalPosIdealWriterConfig.hgcTimeSlices[t].hgcSectors += [ cms.PSet(hgcLinks = cms.vint32(*[27+12*i+4*t+j for j in range(4)])) for i in range(3) ] # pos + hgcalPosIdealWriterConfig.gmtTimeSlices[t].gmtLink = cms.int32(27+36+t) +hgcalPosIdealWriterConfig.gttLink = 27+36+3 +hgcalPosIdealWriterConfig.outputLinksPuppi = cms.vuint32(0,1,2) +hgcalPosIdealWriterConfig.outputLinkEgamma = cms.int32(5) +hgcalPosIdealWriterConfig.inputFileName = cms.string("l1HGCalPos-inputs-ideal") +hgcalPosIdealWriterConfig.outputFileName = cms.string("l1HGCalPos-outputs-ideal") +hgcalNegIdealWriterConfig = hgcalPosIdealWriterConfig.clone( + inputFileName = "", + outputFileName = "l1HGCalNeg-outputs-ideal", + outputRegions = [i for i in range(9)], + outputBoard = 0, +) +## Current configuration for VU9P at B904 for layer1 - layer2 tests with puppi and e/gamma outputs on links 56-59 +hgcalPosVU9PB904egWriterConfig = hgcalWriterConfig_.clone() +for t in range(3): + hgcalPosVU9PB904egWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(-1)) for i in range(9) ] # neg + hgcalPosVU9PB904egWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(3*i+t+4*2)) for i in range(4) ] # pos, left quads + hgcalPosVU9PB904egWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(3*i+t+4*25)) for i in range(5) ] # pos, right quads + hgcalPosVU9PB904egWriterConfig.hgcTimeSlices[t].hgcSectors += [ cms.PSet(hgcLinks = cms.vint32(-1,-1,-1,-1)) for i in range(3) ] # neg + hgcalPosVU9PB904egWriterConfig.hgcTimeSlices[t].hgcSectors += [ cms.PSet(hgcLinks = cms.vint32(*[4*11+12*i+4*t+j for j in range(4)])) for i in range(3) ] # pos + hgcalPosVU9PB904egWriterConfig.gmtTimeSlices[t].gmtLink = cms.int32(4+t) +hgcalPosVU9PB904egWriterConfig.gttLink = 4+3 +hgcalPosVU9PB904egWriterConfig.outputLinksPuppi = cms.vuint32(56,57,58) +hgcalPosVU9PB904egWriterConfig.outputLinkEgamma = cms.int32(59) +hgcalPosVU9PB904egWriterConfig.inputFileName = cms.string("l1HGCalPos-inputs-vu9p_B904eg") +hgcalPosVU9PB904egWriterConfig.outputFileName = cms.string("l1HGCalPos-outputs-vu9p_B904eg") +## Current configuration for VU13P +hgcalPosVU13PWriterConfig = hgcalWriterConfig_.clone() +for t in range(3): + hgcalPosVU13PWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(-1)) for i in range(9) ] # neg + hgcalPosVU13PWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(3*i+t+4*0)) for i in range(5) ] # pos, left quads + hgcalPosVU13PWriterConfig.tfTimeSlices[t].tfSectors += [ cms.PSet(tfLink = cms.int32(3*i+t+4*28)) for i in range(4) ] # pos, right quads + hgcalPosVU13PWriterConfig.hgcTimeSlices[t].hgcSectors += [ cms.PSet(hgcLinks = cms.vint32(-1,-1,-1,-1)) for i in range(3) ] # neg + for isec,q0 in (0,12),(1,17),(2,20): + hgcalPosVU13PWriterConfig.hgcTimeSlices[t].hgcSectors += [ cms.PSet(hgcLinks = cms.vint32(*[4*q0+4*t+j for j in range(4)])) ] # pos + hgcalPosVU13PWriterConfig.gmtTimeSlices[t].gmtLink = cms.int32(4*27+t) +hgcalPosVU13PWriterConfig.gttLink = 4*27+3 +hgcalPosVU13PWriterConfig.outputLinksPuppi = cms.vuint32(0,1,2) +hgcalPosVU13PWriterConfig.outputLinkEgamma = cms.int32(3) +hgcalPosVU13PWriterConfig.inputFileName = cms.string("l1HGCalPos-inputs-vu13p") +hgcalPosVU13PWriterConfig.outputFileName = cms.string("l1HGCalPos-outputs-vu13p") +## Enable both + +## Enable both + +hgcalWriterConfigs = [ + hgcalPosIdealWriterConfig, + hgcalNegIdealWriterConfig, + hgcalPosVU9PB904egWriterConfig, + hgcalPosVU13PWriterConfig +] + +##################################################################################################################### +## Forward HGCal configuration: only outputs for now, 18 regions, 12 candidates x region, 4 output fibers +hgcalNoTKWriterOutputOnlyConfig = cms.PSet( + partition = cms.string("HGCalNoTk"), + outputRegions = cms.vuint32(*range(18)), + nOutputFramesPerBX = cms.uint32(9), + fileFormat = cms.string("EMP"), + maxLinesPerOutputFile = cms.uint32(1024), + eventsPerFile = cms.uint32(eventsPerFile_), + outputLinksPuppi = cms.vuint32(0,1,2,4), + outputFileName = cms.string("l1HGCalNoTk-outputs-ideal") +) + +hgcalNoTKWriterConfigs = [ + hgcalNoTKWriterOutputOnlyConfig +] + +##################################################################################################################### +## HF configuration: not enabled for the moment +##################################################################################################################### +## HF configuration not realistic, 3 links per endcap, write out the layer 1 outputs for now +hfWriterOutputOnly_ = cms.PSet( + partition = cms.string("HF"), + outputLinksPuppi = cms.vuint32(*range(3)), + nOutputFramesPerBX = cms.uint32(9), + fileFormat = cms.string("EMP"), + maxLinesPerOutputFile = cms.uint32(1024), + eventsPerFile = cms.uint32(eventsPerFile_), +) +hfWriterConfigs = [ + hfWriterOutputOnly_.clone( + outputRegions = cms.vuint32(*[9*ie+i for i in range(9)]), + outputFileName = cms.string("l1HF%s-outputs-ideal" % ("Pos" if ie else "Neg")), + ) for ie in range(2) +] + + diff --git a/L1Trigger/Phase2L1ParticleFlow/src/L1TCorrelatorLayer1PatternFileWriter.cc b/L1Trigger/Phase2L1ParticleFlow/src/L1TCorrelatorLayer1PatternFileWriter.cc new file mode 100644 index 0000000000000..018f156978ca0 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/src/L1TCorrelatorLayer1PatternFileWriter.cc @@ -0,0 +1,355 @@ +#include "L1Trigger/Phase2L1ParticleFlow/interface/L1TCorrelatorLayer1PatternFileWriter.h" +#include "FWCore/Utilities/interface/Exception.h" +#include + +L1TCorrelatorLayer1PatternFileWriter::L1TCorrelatorLayer1PatternFileWriter(const edm::ParameterSet& iConfig, + const l1ct::Event& eventTemplate) + : partition_(parsePartition(iConfig.getParameter("partition"))), + writeInputs_(iConfig.existsAs("inputFileName") && + !iConfig.getParameter("inputFileName").empty()), + writeOutputs_(iConfig.existsAs("outputFileName") && + !iConfig.getParameter("outputFileName").empty()), + outputBoard_(-1), + outputLinkEgamma_(-1), + fileFormat_(iConfig.getParameter("fileFormat")), + eventsPerFile_(iConfig.getParameter("eventsPerFile")), + eventIndex_(0) { + if (writeInputs_) { + nInputFramesPerBX_ = iConfig.getParameter("nInputFramesPerBX"); + + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal) { + configTimeSlices(iConfig, "tf", eventTemplate.raw.track.size(), tfTimeslices_, tfLinksFactor_); + channelSpecsInput_["tf"] = {tmuxFactor_ * tfTimeslices_, tfTimeslices_}; + } + if (partition_ == Partition::Barrel) { + auto sectorConfig = iConfig.getParameter>("gctSectors"); + if (sectorConfig.size() != gctSectors_) + throw cms::Exception("Configuration", "Bad number of GCT sectors"); + for (unsigned int iS = 0; iS < gctSectors_; ++iS) { + auto linksEcal = sectorConfig[iS].getParameter>("gctLinksEcal"); + auto linksHad = sectorConfig[iS].getParameter>("gctLinksHad"); + if (linksEcal.size() != gctLinksEcal_ || linksHad.size() != gctLinksHad_) + throw cms::Exception("Configuration", "Bad number of GCT links"); + unsigned int iLink = 0; + for (unsigned int i = 0; i < gctLinksHad_; ++i, ++iLink) { + if (linksHad[i] != -1) + channelIdsInput_[l1t::demo::LinkId{"gct", iLink + 10 * iS}].push_back(linksHad[i]); + } + for (unsigned int i = 0; i < gctLinksEcal_; ++i) { + if (linksEcal[i] != -1) + channelIdsInput_[l1t::demo::LinkId{"gct", iLink + 10 * iS}].push_back(linksEcal[i]); + } + channelSpecsInput_["gct"] = {tmuxFactor_ * gctTimeslices_, 0}; + } + } + if (partition_ == Partition::HGCal || partition_ == Partition::HGCalNoTk) { + configTimeSlices(iConfig, "hgc", eventTemplate.raw.hgcalcluster.size(), hgcTimeslices_, hgcLinksFactor_); + channelSpecsInput_["hgc"] = {tmuxFactor_ * hgcTimeslices_, hgcTimeslices_}; + } + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal || partition_ == Partition::HGCalNoTk) { + configTimeSlices(iConfig, "gmt", 1, gmtTimeslices_, gmtLinksFactor_); + gmtNumberOfMuons_ = iConfig.getParameter("gmtNumberOfMuons"); + channelSpecsInput_["gmt"] = {tmuxFactor_ * gmtTimeslices_, + gmtTimeslices_ * nInputFramesPerBX_ * tmuxFactor_ - gmtNumberOfMuons_}; + } + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal) { + configTimeSlices(iConfig, "gtt", 1, gttTimeslices_, gttLinksFactor_); + gttLatency_ = iConfig.getParameter("gttLatency"); + gttNumberOfPVs_ = iConfig.getParameter("gttNumberOfPVs"); + channelSpecsInput_["gtt"] = l1t::demo::ChannelSpec{tmuxFactor_, gttTimeslices_, gttLatency_}; + } + inputFileWriter_ = + std::make_unique(l1t::demo::parseFileFormat(fileFormat_), + iConfig.getParameter("inputFileName"), + nInputFramesPerBX_, + tmuxFactor_, + iConfig.getParameter("maxLinesPerInputFile"), + channelIdsInput_, + channelSpecsInput_); + } + + if (writeOutputs_) { + nOutputFramesPerBX_ = iConfig.getParameter("nOutputFramesPerBX"); + + outputRegions_ = iConfig.getParameter>("outputRegions"); + outputLinksPuppi_ = iConfig.getParameter>("outputLinksPuppi"); + for (unsigned int i = 0; i < outputLinksPuppi_.size(); ++i) { + channelIdsOutput_[l1t::demo::LinkId{"puppi", i}].push_back(outputLinksPuppi_[i]); + } + channelSpecsOutput_["puppi"] = {tmuxFactor_, 0}; + nPuppiFramesPerRegion_ = (nOutputFramesPerBX_ * tmuxFactor_) / outputRegions_.size(); + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal) { + outputBoard_ = iConfig.getParameter("outputBoard"); + outputLinkEgamma_ = iConfig.getParameter("outputLinkEgamma"); + nEgammaObjectsOut_ = iConfig.getParameter("nEgammaObjectsOut"); + if (outputLinkEgamma_ != -1) { + channelIdsOutput_[l1t::demo::LinkId{"egamma", 0}].push_back(outputLinkEgamma_); + channelSpecsOutput_["egamma"] = {tmuxFactor_, nOutputFramesPerBX_ * tmuxFactor_ - 3 * nEgammaObjectsOut_}; + } + } + if ((outputBoard_ == -1) != (outputLinkEgamma_ == -1)) { + throw cms::Exception("Configuration", "Inconsistent configuration of outputLinkEgamma, outputBoard"); + } + outputFileWriter_ = + std::make_unique(l1t::demo::parseFileFormat(fileFormat_), + iConfig.getParameter("outputFileName"), + nOutputFramesPerBX_, + tmuxFactor_, + iConfig.getParameter("maxLinesPerOutputFile"), + channelIdsOutput_, + channelSpecsOutput_); + } +} + +L1TCorrelatorLayer1PatternFileWriter::~L1TCorrelatorLayer1PatternFileWriter() {} + +void L1TCorrelatorLayer1PatternFileWriter::write(const l1ct::Event& event) { + if (writeInputs_) { + l1t::demo::EventData inputs; + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal) { + writeTF(event, inputs); + } + if (partition_ == Partition::Barrel) { + writeBarrelGCT(event, inputs); + } + if (partition_ == Partition::HGCal || partition_ == Partition::HGCalNoTk) { + writeHGC(event, inputs); + } + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal || partition_ == Partition::HGCalNoTk) { + writeGMT(event, inputs); + } + if (partition_ == Partition::Barrel || partition_ == Partition::HGCal) { + writeGTT(event, inputs); + } + inputFileWriter_->addEvent(inputs); + } + + if (writeOutputs_) { + l1t::demo::EventData outputs; + writePuppi(event, outputs); + if (outputLinkEgamma_ != -1) + writeEgamma(event, outputs); + outputFileWriter_->addEvent(outputs); + } + + eventIndex_++; + if (eventIndex_ % eventsPerFile_ == 0) { + if (writeInputs_) + inputFileWriter_->flush(); + if (writeOutputs_) + outputFileWriter_->flush(); + } +} + +L1TCorrelatorLayer1PatternFileWriter::Partition L1TCorrelatorLayer1PatternFileWriter::parsePartition( + const std::string& partition) { + if (partition == "Barrel") + return Partition::Barrel; + if (partition == "HGCal") + return Partition::HGCal; + if (partition == "HGCalNoTk") + return Partition::HGCalNoTk; + if (partition == "HF") + return Partition::HF; + throw cms::Exception("Configuration", "Unsupported partition_ '" + partition + "'\n"); +} + +void L1TCorrelatorLayer1PatternFileWriter::configTimeSlices(const edm::ParameterSet& iConfig, + const std::string& prefix, + unsigned int nSectors, + unsigned int nTimeSlices, + unsigned int linksFactor) { + if (nTimeSlices > 1) { + auto timeSliceConfig = iConfig.getParameter>(prefix + "TimeSlices"); + if (timeSliceConfig.size() != nTimeSlices) + throw cms::Exception("Configuration") + << "Mismatched number of " << prefix << "TimeSlices, expected " << nTimeSlices << std::endl; + for (unsigned int iT = 0; iT < nTimeSlices; ++iT) { + configSectors(timeSliceConfig[iT], prefix, nSectors, linksFactor); + } + } else { + configSectors(iConfig, prefix, nSectors, linksFactor); + } +} + +void L1TCorrelatorLayer1PatternFileWriter::configSectors(const edm::ParameterSet& iConfig, + const std::string& prefix, + unsigned int nSectors, + unsigned int linksFactor) { + if (nSectors > 1) { + auto sectorConfig = iConfig.getParameter>(prefix + "Sectors"); + if (sectorConfig.size() != nSectors) + throw cms::Exception("Configuration") + << "Mismatched number of " << prefix << "Sectors, expected " << nSectors << std::endl; + for (unsigned int iS = 0; iS < nSectors; ++iS) { + configLinks(sectorConfig[iS], prefix, linksFactor, linksFactor > 1 ? iS * 10 : iS); + } + } else { + configLinks(iConfig, prefix, linksFactor, 0); + } +} +void L1TCorrelatorLayer1PatternFileWriter::configLinks(const edm::ParameterSet& iConfig, + const std::string& prefix, + unsigned int linksFactor, + unsigned int offset) { + if (linksFactor > 1) { + auto links = iConfig.getParameter>(prefix + "Links"); + if (links.size() != linksFactor) + throw cms::Exception("Configuration") + << "Mismatched number of " << prefix << "Links, expected " << linksFactor << std::endl; + for (unsigned int i = 0; i < linksFactor; ++i) { + if (links[i] != -1) { + channelIdsInput_[l1t::demo::LinkId{prefix, i + offset}].push_back(links[i]); + } + } + } else { + auto link = iConfig.getParameter(prefix + "Link"); + if (link != -1) { + channelIdsInput_[l1t::demo::LinkId{prefix, offset}].push_back(link); + } + } +} + +void L1TCorrelatorLayer1PatternFileWriter::writeTF(const l1ct::Event& event, l1t::demo::EventData& out) { + for (unsigned int iS = 0, nS = event.raw.track.size(); iS < nS; ++iS) { + l1t::demo::LinkId key{"tf", iS}; + if (channelIdsInput_.count(key) == 0) + continue; + std::vector> ret; + std::vector> tracks = event.raw.track[iS].obj; + if (tracks.empty()) + tracks.emplace_back(0); + for (unsigned int i = 0, n = tracks.size(); i < n; ++i) { + const ap_uint<96>& packedtk = tracks[i]; + if (i % 2 == 0) { + ret.emplace_back(packedtk(63, 0)); + ret.emplace_back(0); + ret.back()(31, 0) = packedtk(95, 64); + } else { + ret.back()(63, 32) = packedtk(31, 0); + ret.emplace_back(packedtk(95, 32)); + } + } + out.add(key, ret); + } +} + +void L1TCorrelatorLayer1PatternFileWriter::writeHGC(const l1ct::Event& event, l1t::demo::EventData& out) { + assert(hgcLinksFactor_ == 4); // this piece of code won't really work otherwise + std::vector> ret[hgcLinksFactor_]; + for (unsigned int iS = 0, nS = event.raw.hgcalcluster.size(); iS < nS; ++iS) { + l1t::demo::LinkId key0{"hgc", iS * 10}; + if (channelIdsInput_.count(key0) == 0) + continue; + for (unsigned int il = 0; il < hgcLinksFactor_; ++il) { + // put header word and (dummy) towers + ret[il].resize(31); + ap_uint<64>& head64 = ret[il][0]; + head64(63, 48) = 0xABC0; // Magic + head64(47, 38) = 0; // Opaque + head64(39, 32) = (eventIndex_ % 3) * 6; // TM slice + head64(31, 24) = iS; // Sector + head64(23, 16) = il; // link + head64(15, 0) = eventIndex_ % 3564; // BX + for (unsigned int j = 0; j < 30; ++j) { + ret[il][j + 1] = 4 * j + il; + } + } + for (auto clust : event.raw.hgcalcluster[iS].obj) { + for (unsigned int il = 0; il < hgcLinksFactor_; ++il) { + ret[il].push_back(clust(64 * il + 63, 64 * il)); + } + } + for (unsigned int il = 0; il < hgcLinksFactor_; ++il) { + out.add(l1t::demo::LinkId{"hgc", iS * 10 + il}, ret[il]); + } + } +} + +void L1TCorrelatorLayer1PatternFileWriter::writeBarrelGCT(const l1ct::Event& event, l1t::demo::EventData& out) { + std::vector> ret; + for (unsigned int iS = 0; iS < gctSectors_; ++iS) { + l1t::demo::LinkId key0{"gct", iS * 10}; + if (channelIdsInput_.count(key0) == 0) + continue; + const auto& had = event.decoded.hadcalo[iS]; + const auto& ecal = event.decoded.emcalo[iS]; + unsigned int iLink = 0, nHad = had.size(), nEcal = ecal.size(); + for (unsigned int i = 0; i < gctLinksHad_; ++i, ++iLink) { + ret.clear(); + for (unsigned int iHad = i; iHad < nHad; iHad += gctLinksHad_) { + ret.emplace_back(had[iHad].pack()); + } + if (ret.empty()) + ret.emplace_back(0); + out.add(l1t::demo::LinkId{"gct", iS * 10 + iLink}, ret); + } + for (unsigned int i = 0; i < gctLinksEcal_; ++i, ++iLink) { + ret.clear(); + for (unsigned int iEcal = i; iEcal < nEcal; iEcal += gctLinksEcal_) { + ret.emplace_back(ecal[iEcal].pack()); + } + if (ret.empty()) + ret.emplace_back(0); + out.add(l1t::demo::LinkId{"gct", iS * 10 + iLink}, ret); + } + } +} + +void L1TCorrelatorLayer1PatternFileWriter::writeGMT(const l1ct::Event& event, l1t::demo::EventData& out) { + l1t::demo::LinkId key{"gmt", 0}; + if (channelIdsInput_.count(key) == 0) + return; + std::vector> muons = event.raw.muon.obj; + muons.resize(gmtNumberOfMuons_, ap_uint<64>(0)); + out.add(key, muons); +} + +void L1TCorrelatorLayer1PatternFileWriter::writeGTT(const l1ct::Event& event, l1t::demo::EventData& out) { + l1t::demo::LinkId key{"gtt", 0}; + if (channelIdsInput_.count(key) == 0) + return; + std::vector> pvs = event.pvs_emu; + pvs.resize(gttNumberOfPVs_, ap_uint<64>(0)); + out.add(key, pvs); +} + +void L1TCorrelatorLayer1PatternFileWriter::writePuppi(const l1ct::Event& event, l1t::demo::EventData& out) { + unsigned int n = outputLinksPuppi_.size(); + std::vector>> links(n); + for (auto ir : outputRegions_) { + auto puppi = event.out[ir].puppi; + unsigned int npuppi = puppi.size(); + for (unsigned int i = 0; i < n * nPuppiFramesPerRegion_; ++i) { + links[i / nPuppiFramesPerRegion_].push_back(i < npuppi ? puppi[i].pack() : ap_uint(0)); + } + } + for (unsigned int i = 0; i < n; ++i) { + out.add(l1t::demo::LinkId{"puppi", i}, links[i]); + } +} + +void L1TCorrelatorLayer1PatternFileWriter::writeEgamma(const l1ct::Event& event, l1t::demo::EventData& out) { + std::vector> ret; + const auto& pho = event.board_out[outputBoard_].egphoton; + const auto& ele = event.board_out[outputBoard_].egelectron; + ret.reserve(3 * nEgammaObjectsOut_); + for (const auto& p : pho) { + ret.emplace_back(p.pack()); + } + ret.resize(nEgammaObjectsOut_, ap_uint<64>(0)); + for (const auto& p : ele) { + ap_uint<128> dword = p.pack(); + ret.push_back(dword(63, 0)); + ret.push_back(dword(127, 64)); + } + ret.resize(3 * nEgammaObjectsOut_, ap_uint<64>(0)); + out.add(l1t::demo::LinkId{"egamma", 0}, ret); +} + +void L1TCorrelatorLayer1PatternFileWriter::flush() { + if (inputFileWriter_) + inputFileWriter_->flush(); + if (outputFileWriter_) + outputFileWriter_->flush(); +} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/PFAlgo2HGC.cc b/L1Trigger/Phase2L1ParticleFlow/src/PFAlgo2HGC.cc deleted file mode 100644 index e5304af104e3c..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/PFAlgo2HGC.cc +++ /dev/null @@ -1,652 +0,0 @@ -#include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo2HGC.h" -#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" - -#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" - -#include "FWCore/Utilities/interface/Exception.h" - -#include "DataFormats/Math/interface/deltaR.h" - -namespace { - template - float floatDR(const T1 &t1, const T2 &t2) { - return deltaR(t1.floatEta(), t1.floatPhi(), t2.floatEta(), t2.floatPhi()); - } -} // namespace - -using namespace l1tpf_impl; - -PFAlgo2HGC::PFAlgo2HGC(const edm::ParameterSet &iConfig) : PFAlgoBase(iConfig) { - debug_ = iConfig.getUntrackedParameter("debugPFAlgo2HGC", iConfig.getUntrackedParameter("debug", 0)); - edm::ParameterSet linkcfg = iConfig.getParameter("linking"); - drMatchMu_ = linkcfg.getParameter("trackMuDR"); - - std::string muMatchMode = linkcfg.getParameter("trackMuMatch"); - if (muMatchMode == "boxBestByPtRatio") - muMatchMode_ = MuMatchMode::BoxBestByPtRatio; - else if (muMatchMode == "drBestByPtRatio") - muMatchMode_ = MuMatchMode::DrBestByPtRatio; - else if (muMatchMode == "drBestByPtDiff") - muMatchMode_ = MuMatchMode::DrBestByPtDiff; - else - throw cms::Exception("Configuration", "bad value for trackMuMatch configurable"); - - std::string tkCaloLinkMetric = linkcfg.getParameter("trackCaloLinkMetric"); - if (tkCaloLinkMetric == "bestByDR") - tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDR; - else if (tkCaloLinkMetric == "bestByDRPt") - tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDRPt; - else if (tkCaloLinkMetric == "bestByDR2Pt2") - tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDR2Pt2; - else - throw cms::Exception("Configuration", "bad value for tkCaloLinkMetric configurable"); - - drMatch_ = linkcfg.getParameter("trackCaloDR"); - ptMatchLow_ = linkcfg.getParameter("trackCaloNSigmaLow"); - ptMatchHigh_ = linkcfg.getParameter("trackCaloNSigmaHigh"); - useTrackCaloSigma_ = linkcfg.getParameter("useTrackCaloSigma"); - maxInvisiblePt_ = linkcfg.getParameter("maxInvisiblePt"); - - caloReLinkStep_ = linkcfg.getParameter("caloReLink"); - caloReLinkDr_ = linkcfg.getParameter("caloReLinkDR"); - caloReLinkThreshold_ = linkcfg.getParameter("caloReLinkThreshold"); - rescaleTracks_ = linkcfg.getParameter("rescaleTracks"); - caloTrkWeightedAverage_ = linkcfg.getParameter("useCaloTrkWeightedAverage"); - sumTkCaloErr2_ = linkcfg.getParameter("sumTkCaloErr2"); - ecalPriority_ = linkcfg.getParameter("ecalPriority"); - tightTrackMaxInvisiblePt_ = linkcfg.getParameter("tightTrackMaxInvisiblePt"); - sortInputs_ = iConfig.getParameter("sortInputs"); -} - -void PFAlgo2HGC::runPF(Region &r) const { - initRegion(r, sortInputs_); - - /// ------------- first step (can all go in parallel) ---------------- - - if (debug_) { - dbgPrintf( - "PFAlgo2HGC\nPFAlgo2HGC region eta [ %+5.2f , %+5.2f ], phi [ %+5.2f , %+5.2f ], fiducial eta [ %+5.2f , " - "%+5.2f ], phi [ %+5.2f , %+5.2f ]\n", - r.etaMin - r.etaExtra, - r.etaMax + r.etaExtra, - r.phiCenter - r.phiHalfWidth - r.phiExtra, - r.phiCenter + r.phiHalfWidth + r.phiExtra, - r.etaMin, - r.etaMax, - r.phiCenter - r.phiHalfWidth, - r.phiCenter + r.phiHalfWidth); - dbgPrintf( - "PFAlgo2HGC \t N(track) %3lu N(calo) %3lu N(mu) %3lu\n", r.track.size(), r.calo.size(), r.muon.size()); - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - dbgPrintf( - "PFAlgo2HGC \t track %3d: pt %7.2f +- %5.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi " - "%+5.2f fid %1d calo ptErr %7.2f stubs %2d chi2 %7.1f quality %d\n", - itk, - tk.floatPt(), - tk.floatPtErr(), - tk.floatVtxEta(), - tk.floatVtxPhi(), - tk.floatEta(), - tk.floatPhi(), - int(r.fiducialLocal(tk.floatEta(), tk.floatPhi())), - tk.floatCaloPtErr(), - int(tk.hwStubs), - tk.hwChi2 * 0.1f, - int(tk.hwFlags)); - } - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - dbgPrintf( - "PFAlgo2HGC \t calo %3d: pt %7.2f +- %5.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi " - "%+5.2f fid %1d calo ptErr %7.2f em pt %7.2f isEM %1d \n", - ic, - calo.floatPt(), - calo.floatPtErr(), - calo.floatEta(), - calo.floatPhi(), - calo.floatEta(), - calo.floatPhi(), - int(r.fiducialLocal(calo.floatEta(), calo.floatPhi())), - calo.floatPtErr(), - calo.floatEmPt(), - calo.isEM); - } - for (int im = 0, nm = r.muon.size(); im < nm; ++im) { - auto &mu = r.muon[im]; - dbgPrintf( - "PFAlgo2HGC \t muon %3d: pt %7.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi " - "%+5.2f fid %1d\n", - im, - mu.floatPt(), - mu.floatEta(), - mu.floatPhi(), - mu.floatEta(), - mu.floatPhi(), - int(r.fiducialLocal(mu.floatEta(), mu.floatPhi()))); - } - } - - std::vector tk2mu(r.track.size(), -1), mu2tk(r.muon.size(), -1); - link_tk2mu(r, tk2mu, mu2tk); - - // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later) - std::vector tk2calo(r.track.size(), -1); - link_tk2calo(r, tk2calo); - - /// ------------- next step (needs the previous) ---------------- - // for each calo, compute the sum of the track pt - std::vector calo2ntk(r.calo.size(), 0); - std::vector calo2sumtkpt(r.calo.size(), 0); - std::vector calo2sumtkpterr(r.calo.size(), 0); - sum_tk2calo(r, tk2calo, calo2ntk, calo2sumtkpt, calo2sumtkpterr); - - // in the meantime, promote unlinked low pt tracks to hadrons - unlinkedtk_algo(r, tk2calo); - - /// ------------- next step (needs the previous) ---------------- - /// OPTIONAL STEP: try to recover split hadron showers (v1.0): - // off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events - if (caloReLinkStep_) - calo_relink(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr); - - /// ------------- next step (needs the previous) ---------------- - // process matched calo clusters, compare energy to sum track pt - std::vector calo2alpha(r.calo.size(), 1); - linkedcalo_algo(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr, calo2alpha); - - /// ------------- next step (needs the previous) ---------------- - /// process matched tracks, if necessary rescale or average - linkedtk_algo(r, tk2calo, calo2ntk, calo2alpha); - // process unmatched calo clusters - unlinkedcalo_algo(r); - // finally do muons - save_muons(r, tk2mu); -} - -void PFAlgo2HGC::link_tk2mu(Region &r, std::vector &tk2mu, std::vector &mu2tk) const { - // do a rectangular match for the moment; make a box of the same are as a 0.2 cone - int intDrMuonMatchBox = std::ceil(drMatchMu_ * CaloCluster::ETAPHI_SCALE * std::sqrt(M_PI / 4)); - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - tk2mu[itk] = false; - } - for (int imu = 0, nmu = r.muon.size(); imu < nmu; ++imu) { - const auto &mu = r.muon[imu]; - if (debug_) - dbgPrintf("PFAlgo2HGC \t muon %3d (pt %7.2f, eta %+5.2f, phi %+5.2f) \n", - imu, - mu.floatPt(), - mu.floatEta(), - mu.floatPhi()); - float minDistance = 9e9; - switch (muMatchMode_) { - case MuMatchMode::BoxBestByPtRatio: - minDistance = 4.; - break; - case MuMatchMode::DrBestByPtRatio: - minDistance = 4.; - break; - case MuMatchMode::DrBestByPtDiff: - minDistance = 0.5 * mu.floatPt(); - break; - } - int imatch = -1; - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - int deta = std::abs(mu.hwEta - tk.hwEta); - int dphi = std::abs((mu.hwPhi - tk.hwPhi) % CaloCluster::PHI_WRAP); - float dr = floatDR(mu, tk); - float dpt = std::abs(mu.floatPt() - tk.floatPt()); - float dptr = (mu.hwPt > tk.hwPt ? mu.floatPt() / tk.floatPt() : tk.floatPt() / mu.floatPt()); - bool ok = false; - float distance = 9e9; - switch (muMatchMode_) { - case MuMatchMode::BoxBestByPtRatio: - ok = (deta < intDrMuonMatchBox) && (dphi < intDrMuonMatchBox); - distance = dptr; - break; - case MuMatchMode::DrBestByPtRatio: - ok = (dr < drMatchMu_); - distance = dptr; - break; - case MuMatchMode::DrBestByPtDiff: - ok = (dr < drMatchMu_); - distance = dpt; - break; - } - if (debug_ && dr < 0.4) { - dbgPrintf( - "PFAlgo2HGC \t\t possible match with track %3d (pt %7.2f, caloeta %+5.2f, calophi %+5.2f, dr %.2f, eta " - "%+5.2f, phi %+5.2f, dr %.2f): angular %1d, distance %.3f (vs %.3f)\n", - itk, - tk.floatPt(), - tk.floatEta(), - tk.floatPhi(), - dr, - tk.floatVtxEta(), - tk.floatVtxPhi(), - deltaR(mu.floatEta(), mu.floatPhi(), tk.floatVtxEta(), tk.floatVtxPhi()), - (ok ? 1 : 0), - distance, - minDistance); - } - if (!ok) - continue; - // FIXME for the moment, we do the floating point matching in pt - if (distance < minDistance) { - minDistance = distance; - imatch = itk; - } - } - if (debug_ && imatch > -1) - dbgPrintf("PFAlgo2HGC \t muon %3d (pt %7.2f) linked to track %3d (pt %7.2f)\n", - imu, - mu.floatPt(), - imatch, - r.track[imatch].floatPt()); - if (debug_ && imatch == -1) - dbgPrintf("PFAlgo2HGC \t muon %3d (pt %7.2f) not linked to any track\n", imu, mu.floatPt()); - mu2tk[imu] = imatch; - if (imatch > -1) { - tk2mu[imatch] = imu; - r.track[imatch].muonLink = true; - } - } -} - -void PFAlgo2HGC::link_tk2calo(Region &r, std::vector &tk2calo) const { - // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later) - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - if (tk.muonLink || tk.used) - continue; // not necessary but just a waste of CPU otherwise - float drbest = drMatch_, dptscale = 0; - switch (tkCaloLinkMetric_) { - case TkCaloLinkMetric::BestByDR: - drbest = drMatch_; - break; - case TkCaloLinkMetric::BestByDRPt: - drbest = 1.0; - dptscale = drMatch_ / tk.floatCaloPtErr(); - break; - case TkCaloLinkMetric::BestByDR2Pt2: - drbest = 1.0; - dptscale = drMatch_ / tk.floatCaloPtErr(); - break; - } - float minCaloPt = tk.floatPt() - ptMatchLow_ * tk.floatCaloPtErr(); - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t track %3d (pt %7.2f) to be matched to calo, min pT %7.2f\n", itk, tk.floatPt(), minCaloPt); - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo.used || calo.floatPt() <= minCaloPt) - continue; - float dr = floatDR(tk, calo), dq; - switch (tkCaloLinkMetric_) { - case TkCaloLinkMetric::BestByDR: - if (dr < drbest) { - tk2calo[itk] = ic; - drbest = dr; - } - break; - case TkCaloLinkMetric::BestByDRPt: - dq = dr + std::max(tk.floatPt() - calo.floatPt(), 0.) * dptscale; - if (debug_ > 2 && dr < 0.3) - dbgPrintf("PFAlgo2HGC \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n", - itk, - tk.floatPt(), - ic, - calo.floatPt(), - dr, - dq); - if (dr < drMatch_ && dq < drbest) { - tk2calo[itk] = ic; - drbest = dq; - } - break; - case TkCaloLinkMetric::BestByDR2Pt2: - dq = hypot(dr, std::max(tk.floatPt() - calo.floatPt(), 0.) * dptscale); - if (debug_ > 2 && dr < 0.3) - dbgPrintf("PFAlgo2HGC \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n", - itk, - tk.floatPt(), - ic, - calo.floatPt(), - dr, - dq); - if (dr < drMatch_ && dq < drbest) { - tk2calo[itk] = ic; - drbest = dq; - } - break; - } - } - if (debug_ && tk2calo[itk] != -1) - dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) matches to calo %3d (pt %7.2f) with dist %.3f (dr %.3f)\n", - itk, - tk.floatPt(), - tk2calo[itk], - r.calo[tk2calo[itk]].floatPt(), - drbest, - floatDR(tk, r.calo[tk2calo[itk]])); - // now we re-do this for debugging sake, it may be done for real later - if (debug_ && tk2calo[itk] == -1) { - int ibest = -1; - drbest = 0.3; - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo.used) - continue; - float dr = floatDR(tk, calo); - if (dr < drbest) { - ibest = ic; - drbest = dr; - } - } - if (ibest != -1) - dbgPrintf( - "PFAlgo2HGC \t track %3d (pt %7.2f) would match to calo %3d (pt %7.2f) with dr %.3f if the pt min and dr " - "requirement had been relaxed\n", - itk, - tk.floatPt(), - ibest, - r.calo[ibest].floatPt(), - drbest); - } - } -} - -void PFAlgo2HGC::sum_tk2calo(Region &r, - const std::vector &tk2calo, - std::vector &calo2ntk, - std::vector &calo2sumtkpt, - std::vector &calo2sumtkpterr) const { - // for each calo, compute the sum of the track pt - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - const auto &calo = r.calo[ic]; - if (r.globalAbsEta(calo.floatEta()) > 2.5) - continue; - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - if (tk2calo[itk] == ic) { - const auto &tk = r.track[itk]; - if (tk.muonLink || tk.used) - continue; - calo2ntk[ic]++; - calo2sumtkpt[ic] += tk.floatPt(); - calo2sumtkpterr[ic] += std::pow(tk.floatCaloPtErr(), sumTkCaloErr2_ ? 2 : 1); - } - } - if (sumTkCaloErr2_ && calo2sumtkpterr[ic] > 0) - calo2sumtkpterr[ic] = std::sqrt(calo2sumtkpterr[ic]); - } -} - -void PFAlgo2HGC::unlinkedtk_algo(Region &r, const std::vector &tk2calo) const { - // in the meantime, promote unlinked low pt tracks to hadrons - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - auto &tk = r.track[itk]; - if (tk2calo[itk] != -1 || tk.muonLink || tk.used || !tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - float maxPt = tk.quality(l1tpf_impl::InputTrack::PFTIGHT) ? tightTrackMaxInvisiblePt_ : maxInvisiblePt_; - if (tk.floatPt() < maxPt) { - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t track %3d (pt %7.2f) not matched to calo, kept as charged hadron\n", itk, tk.floatPt()); - auto &p = addTrackToPF(r, tk); - p.hwStatus = GoodTK_NoCalo; - tk.used = true; - } else { - if (debug_) - dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) not matched to calo, dropped\n", itk, tk.floatPt()); - } - } -} - -void PFAlgo2HGC::calo_relink(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr) const { - /// OPTIONAL STEP: try to recover split hadron showers (v1.0): - // take hadrons that are not track matched, close by a hadron which has an excess of track pt vs calo pt - // add this pt to the calo pt of the other cluster - // off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events - std::vector addtopt(r.calo.size(), 0); - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo2ntk[ic] != 0 || calo.used || r.globalAbsEta(calo.floatEta()) > 2.5) - continue; - int i2best = -1; - float drbest = caloReLinkDr_; - for (int ic2 = 0; ic2 < nc; ++ic2) { - const auto &calo2 = r.calo[ic2]; - if (calo2ntk[ic2] == 0 || calo2.used || r.globalAbsEta(calo2.floatEta()) > 2.5) - continue; - float dr = floatDR(calo, calo2); - //// uncomment below for more verbose debugging - //if (debug_ && dr < 0.5) dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f) with no tracks is at dr %.3f from calo %3d with pt %7.2f (sum tk pt %7.2f), track excess %7.2f +- %7.2f\n", ic, calo.floatPt(), dr, ic2, calo2.floatPt(), calo2sumtkpt[ic2], calo2sumtkpt[ic2] - calo2.floatPt(), useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr()); - if (dr < drbest) { - float ptdiff = - calo2sumtkpt[ic2] - calo2.floatPt() + (useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr()); - if (ptdiff >= caloReLinkThreshold_ * calo.floatPt()) { - i2best = ic2; - drbest = dr; - } - } - } - if (i2best != -1) { - const auto &calo2 = r.calo[i2best]; - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t calo %3d (pt %7.2f) with no tracks matched within dr %.3f with calo %3d with pt %7.2f (sum " - "tk pt %7.2f), track excess %7.2f +- %7.2f\n", - ic, - calo.floatPt(), - drbest, - i2best, - calo2.floatPt(), - calo2sumtkpt[i2best], - calo2sumtkpt[i2best] - calo2.floatPt(), - useTrackCaloSigma_ ? calo2sumtkpterr[i2best] : calo2.floatPtErr()); - calo.used = true; - addtopt[i2best] += calo.floatPt(); - } - } - // we do this at the end, so that the above loop is parallelizable - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - if (addtopt[ic]) { - auto &calo = r.calo[ic]; - if (debug_) - dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f, sum tk pt %7.2f) is increased to pt %7.2f after merging\n", - ic, - calo.floatPt(), - calo2sumtkpt[ic], - calo.floatPt() + addtopt[ic]); - calo.setFloatPt(calo.floatPt() + addtopt[ic]); - } - } -} - -void PFAlgo2HGC::linkedcalo_algo(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr, - std::vector &calo2alpha) const { - /// ------------- next step (needs the previous) ---------------- - // process matched calo clusters, compare energy to sum track pt - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo2ntk[ic] == 0 || calo.used) - continue; - float ptdiff = calo.floatPt() - calo2sumtkpt[ic]; - float pterr = useTrackCaloSigma_ ? calo2sumtkpterr[ic] : calo.floatPtErr(); - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t calo %3d (pt %7.2f +- %7.2f, empt %7.2f) has %2d tracks (sumpt %7.2f, sumpterr %7.2f), ptdif " - "%7.2f +- %7.2f\n", - ic, - calo.floatPt(), - calo.floatPtErr(), - calo.floatEmPt(), - calo2ntk[ic], - calo2sumtkpt[ic], - calo2sumtkpterr[ic], - ptdiff, - pterr); - if (ptdiff > +ptMatchHigh_ * pterr) { - if (ecalPriority_) { - if (calo.floatEmPt() > 1) { - float emptdiff = std::min(ptdiff, calo.floatEmPt()); - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t calo %3d (pt %7.2f, empt %7.2f) ---> make photon with pt %7.2f, reduce ptdiff to " - "%7.2f +- %7.2f\n", - ic, - calo.floatPt(), - calo.floatEmPt(), - emptdiff, - ptdiff - emptdiff, - pterr); - auto &p = addCaloToPF(r, calo); - p.setFloatPt(emptdiff); - p.hwId = l1t::PFCandidate::Photon; - ptdiff -= emptdiff; - } - if (ptdiff > 2) { - if (debug_) - dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f, empt %7.2f) ---> make also neutral hadron with pt %7.2f\n", - ic, - calo.floatPt(), - calo.floatEmPt(), - ptdiff); - auto &p = addCaloToPF(r, calo); - p.setFloatPt(ptdiff); - p.hwId = l1t::PFCandidate::NeutralHadron; - } - } else { - if (debug_) - dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f) ---> promoted to neutral with pt %7.2f\n", - ic, - calo.floatPt(), - ptdiff); - auto &p = addCaloToPF(r, calo); - p.setFloatPt(ptdiff); - calo.hwFlags = 0; - } - } else if (ptdiff > -ptMatchLow_ * pterr) { - // nothing to do (weighted average happens when we process the tracks) - calo.hwFlags = 1; - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t calo %3d (pt %7.2f) ---> to be deleted, will use tracks instead\n", ic, calo.floatPt()); - } else { - // tracks overshoot, rescale to tracks to calo - calo2alpha[ic] = rescaleTracks_ ? calo.floatPt() / calo2sumtkpt[ic] : 1.0; - calo.hwFlags = 2; - if (debug_ && rescaleTracks_) - dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f) ---> tracks overshoot and will be scaled down by %.4f\n", - ic, - calo.floatPt(), - calo2alpha[ic]); - if (debug_ && !rescaleTracks_) - dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f) ---> tracks overshoot by %.4f\n", - ic, - calo.floatPt(), - calo2sumtkpt[ic] / calo.floatPt()); - } - calo.used = true; - } -} - -void PFAlgo2HGC::linkedtk_algo(Region &r, - const std::vector &tk2calo, - const std::vector &calo2ntk, - const std::vector &calo2alpha) const { - // process matched tracks, if necessary rescale or average - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - auto &tk = r.track[itk]; - if (tk2calo[itk] == -1 || tk.muonLink || tk.used) - continue; - auto &p = addTrackToPF(r, tk); - tk.used = true; - const auto &calo = r.calo[tk2calo[itk]]; - if (calo.isEM) - p.hwId = l1t::PFCandidate::Electron; - p.cluster.src = calo.src; - if (calo.hwFlags == 1) { - // can do weighted average if there's just one track - if (calo2ntk[tk2calo[itk]] == 1 && caloTrkWeightedAverage_) { - p.hwStatus = GoodTK_Calo_TkPt; - float ptavg = tk.floatPt(); - if (tk.floatPtErr() > 0) { - float wcalo = 1.0 / std::pow(tk.floatCaloPtErr(), 2); - float wtk = 1.0 / std::pow(tk.floatPtErr(), 2); - ptavg = (calo.floatPt() * wcalo + tk.floatPt() * wtk) / (wcalo + wtk); - p.hwStatus = GoodTK_Calo_TkCaloPt; - } - p.setFloatPt(ptavg); - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t track %3d (pt %7.2f +- %7.2f) combined with calo %3d (pt %7.2f +- %7.2f (from tk) " - "yielding candidate of pt %7.2f\n", - itk, - tk.floatPt(), - tk.floatPtErr(), - tk2calo[itk], - calo.floatPt(), - tk.floatCaloPtErr(), - ptavg); - } else { - p.hwStatus = GoodTK_Calo_TkPt; - if (debug_) - dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) linked to calo %3d promoted to %s\n", - itk, - tk.floatPt(), - tk2calo[itk], - (p.hwId == l1t::PFCandidate::Electron ? "electron" : "charged hadron")); - } - } else if (calo.hwFlags == 2) { - // must rescale - p.setFloatPt(tk.floatPt() * calo2alpha[tk2calo[itk]]); - p.hwStatus = GoodTk_Calo_CaloPt; - if (debug_) - dbgPrintf( - "PFAlgo2HGC \t track %3d (pt %7.2f, stubs %2d chi2 %7.1f) linked to calo %3d promoted to %s with pt %7.2f " - "after maybe rescaling\n", - itk, - tk.floatPt(), - int(tk.hwStubs), - tk.hwChi2 * 0.1f, - tk2calo[itk], - (p.hwId == l1t::PFCandidate::Electron ? "electron" : "charged hadron"), - p.floatPt()); - } - } -} - -void PFAlgo2HGC::unlinkedcalo_algo(Region &r) const { - // process unmatched calo clusters - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - if (!r.calo[ic].used) { - addCaloToPF(r, r.calo[ic]); - if (debug_) - dbgPrintf("PFAlgo2HGC \t calo %3d (pt %7.2f) not linked, promoted to neutral\n", ic, r.calo[ic].floatPt()); - } - } -} - -void PFAlgo2HGC::save_muons(Region &r, const std::vector &tk2mu) const { - // finally do muons - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - if (r.track[itk].muonLink) { - auto &p = addTrackToPF(r, r.track[itk]); - p.muonsrc = r.muon[tk2mu[itk]].src; - if (debug_) - dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) promoted to muon.\n", itk, r.track[itk].floatPt()); - } - } -} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/PFAlgo3.cc b/L1Trigger/Phase2L1ParticleFlow/src/PFAlgo3.cc deleted file mode 100644 index abd94c6c08f1b..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/PFAlgo3.cc +++ /dev/null @@ -1,916 +0,0 @@ -#include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo3.h" -#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" - -#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" - -#include "FWCore/Utilities/interface/Exception.h" - -#include "DataFormats/Math/interface/deltaR.h" - -namespace { - template - float floatDR(const T1 &t1, const T2 &t2) { - return deltaR(t1.floatEta(), t1.floatPhi(), t2.floatEta(), t2.floatPhi()); - } -} // namespace - -using namespace l1tpf_impl; - -PFAlgo3::PFAlgo3(const edm::ParameterSet &iConfig) : PFAlgoBase(iConfig) { - debug_ = iConfig.getUntrackedParameter("debugPFAlgo3", iConfig.getUntrackedParameter("debug", 0)); - edm::ParameterSet linkcfg = iConfig.getParameter("linking"); - drMatchMu_ = linkcfg.getParameter("trackMuDR"); - - std::string muMatchMode = linkcfg.getParameter("trackMuMatch"); - if (muMatchMode == "boxBestByPtRatio") - muMatchMode_ = MuMatchMode::BoxBestByPtRatio; - else if (muMatchMode == "drBestByPtRatio") - muMatchMode_ = MuMatchMode::DrBestByPtRatio; - else if (muMatchMode == "drBestByPtDiff") - muMatchMode_ = MuMatchMode::DrBestByPtDiff; - else - throw cms::Exception("Configuration", "bad value for trackMuMatch configurable"); - - std::string tkCaloLinkMetric = linkcfg.getParameter("trackCaloLinkMetric"); - if (tkCaloLinkMetric == "bestByDR") - tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDR; - else if (tkCaloLinkMetric == "bestByDRPt") - tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDRPt; - else if (tkCaloLinkMetric == "bestByDR2Pt2") - tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDR2Pt2; - else - throw cms::Exception("Configuration", "bad value for tkCaloLinkMetric configurable"); - - drMatch_ = linkcfg.getParameter("trackCaloDR"); - ptMatchLow_ = linkcfg.getParameter("trackCaloNSigmaLow"); - ptMatchHigh_ = linkcfg.getParameter("trackCaloNSigmaHigh"); - useTrackCaloSigma_ = linkcfg.getParameter("useTrackCaloSigma"); - maxInvisiblePt_ = linkcfg.getParameter("maxInvisiblePt"); - - drMatchEm_ = linkcfg.getParameter("trackEmDR"); - trackEmUseAlsoTrackSigma_ = linkcfg.getParameter("trackEmUseAlsoTrackSigma"); - trackEmMayUseCaloMomenta_ = linkcfg.getParameter("trackEmMayUseCaloMomenta"); - emCaloUseAlsoCaloSigma_ = linkcfg.getParameter("emCaloUseAlsoCaloSigma"); - ptMinFracMatchEm_ = linkcfg.getParameter("caloEmPtMinFrac"); - drMatchEmHad_ = linkcfg.getParameter("emCaloDR"); - emHadSubtractionPtSlope_ = linkcfg.getParameter("emCaloSubtractionPtSlope"); - caloReLinkStep_ = linkcfg.getParameter("caloReLink"); - caloReLinkDr_ = linkcfg.getParameter("caloReLinkDR"); - caloReLinkThreshold_ = linkcfg.getParameter("caloReLinkThreshold"); - rescaleTracks_ = linkcfg.getParameter("rescaleTracks"); - caloTrkWeightedAverage_ = linkcfg.getParameter("useCaloTrkWeightedAverage"); - sumTkCaloErr2_ = linkcfg.getParameter("sumTkCaloErr2"); - ecalPriority_ = linkcfg.getParameter("ecalPriority"); - tightTrackMaxInvisiblePt_ = linkcfg.getParameter("tightTrackMaxInvisiblePt"); - sortInputs_ = iConfig.getParameter("sortInputs"); -} - -void PFAlgo3::runPF(Region &r) const { - initRegion(r, sortInputs_); - - /// ------------- first step (can all go in parallel) ---------------- - - if (debug_) { - dbgPrintf( - "PFAlgo3\nPFAlgo3 region eta [ %+5.2f , %+5.2f ], phi [ %+5.2f , %+5.2f ], fiducial eta [ %+5.2f , %+5.2f ], " - "phi [ %+5.2f , %+5.2f ]\n", - r.etaMin - r.etaExtra, - r.etaMax + r.etaExtra, - r.phiCenter - r.phiHalfWidth - r.phiExtra, - r.phiCenter + r.phiHalfWidth + r.phiExtra, - r.etaMin, - r.etaMax, - r.phiCenter - r.phiHalfWidth, - r.phiCenter + r.phiHalfWidth); - dbgPrintf("PFAlgo3 \t N(track) %3lu N(em) %3lu N(calo) %3lu N(mu) %3lu\n", - r.track.size(), - r.emcalo.size(), - r.calo.size(), - r.muon.size()); - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - dbgPrintf( - "PFAlgo3 \t track %3d: pt %7.2f +- %5.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi %+5.2f " - "fid %1d calo ptErr %7.2f stubs %2d chi2 %7.1f quality %d\n", - itk, - tk.floatPt(), - tk.floatPtErr(), - tk.floatVtxEta(), - tk.floatVtxPhi(), - tk.floatEta(), - tk.floatPhi(), - int(r.fiducialLocal(tk.floatEta(), tk.floatPhi())), - tk.floatCaloPtErr(), - int(tk.hwStubs), - tk.hwChi2 * 0.1f, - int(tk.hwFlags)); - } - for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) { - const auto &em = r.emcalo[iem]; - dbgPrintf( - "PFAlgo3 \t EM %3d: pt %7.2f +- %5.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi %+5.2f " - "fid %1d calo ptErr %7.2f\n", - iem, - em.floatPt(), - em.floatPtErr(), - em.floatEta(), - em.floatPhi(), - em.floatEta(), - em.floatPhi(), - int(r.fiducialLocal(em.floatEta(), em.floatPhi())), - em.floatPtErr()); - } - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - dbgPrintf( - "PFAlgo3 \t calo %3d: pt %7.2f +- %5.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi %+5.2f " - "fid %1d calo ptErr %7.2f em pt %7.2f \n", - ic, - calo.floatPt(), - calo.floatPtErr(), - calo.floatEta(), - calo.floatPhi(), - calo.floatEta(), - calo.floatPhi(), - int(r.fiducialLocal(calo.floatEta(), calo.floatPhi())), - calo.floatPtErr(), - calo.floatEmPt()); - } - for (int im = 0, nm = r.muon.size(); im < nm; ++im) { - auto &mu = r.muon[im]; - dbgPrintf( - "PFAlgo3 \t muon %3d: pt %7.2f vtx eta %+5.2f vtx phi %+5.2f calo eta %+5.2f calo phi %+5.2f " - "fid %1d \n", - im, - mu.floatPt(), - mu.floatEta(), - mu.floatPhi(), - mu.floatEta(), - mu.floatPhi(), - int(r.fiducialLocal(mu.floatEta(), mu.floatPhi()))); - } - } - - std::vector tk2mu(r.track.size(), -1), mu2tk(r.muon.size(), -1); - link_tk2mu(r, tk2mu, mu2tk); - - // match all tracks to the closest EM cluster - std::vector tk2em(r.track.size(), -1); - link_tk2em(r, tk2em); - - // match all em to the closest had (can happen in parallel to the above) - std::vector em2calo(r.emcalo.size(), -1); - link_em2calo(r, em2calo); - - /// ------------- next step (needs the previous) ---------------- - // for each EM cluster, count and add up the pt of all the corresponding tracks (skipping muons) - std::vector em2ntk(r.emcalo.size(), 0); - std::vector em2sumtkpt(r.emcalo.size(), 0); - std::vector em2sumtkpterr(r.emcalo.size(), 0); - sum_tk2em(r, tk2em, em2ntk, em2sumtkpt, em2sumtkpterr); - - /// ------------- next step (needs the previous) ---------------- - // process ecal clusters after linking - emcalo_algo(r, em2ntk, em2sumtkpt, em2sumtkpterr); - - /// ------------- next step (needs the previous) ---------------- - // promote all flagged tracks to electrons - emtk_algo(r, tk2em, em2ntk, em2sumtkpterr); - sub_em2calo(r, em2calo); - - /// ------------- next step (needs the previous) ---------------- - // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later) - std::vector tk2calo(r.track.size(), -1); - link_tk2calo(r, tk2calo); - - /// ------------- next step (needs the previous) ---------------- - // for each calo, compute the sum of the track pt - std::vector calo2ntk(r.calo.size(), 0); - std::vector calo2sumtkpt(r.calo.size(), 0); - std::vector calo2sumtkpterr(r.calo.size(), 0); - sum_tk2calo(r, tk2calo, calo2ntk, calo2sumtkpt, calo2sumtkpterr); - - // in the meantime, promote unlinked low pt tracks to hadrons - unlinkedtk_algo(r, tk2calo); - - /// ------------- next step (needs the previous) ---------------- - /// OPTIONAL STEP: try to recover split hadron showers (v1.0): - // off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events - if (caloReLinkStep_) - calo_relink(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr); - - /// ------------- next step (needs the previous) ---------------- - // process matched calo clusters, compare energy to sum track pt - std::vector calo2alpha(r.calo.size(), 1); - linkedcalo_algo(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr, calo2alpha); - - /// ------------- next step (needs the previous) ---------------- - /// process matched tracks, if necessary rescale or average - linkedtk_algo(r, tk2calo, calo2ntk, calo2alpha); - // process unmatched calo clusters - unlinkedcalo_algo(r); - // finally do muons - save_muons(r, tk2mu); -} - -void PFAlgo3::link_tk2mu(Region &r, std::vector &tk2mu, std::vector &mu2tk) const { - // do a rectangular match for the moment; make a box of the same are as a 0.2 cone - int intDrMuonMatchBox = std::ceil(drMatchMu_ * CaloCluster::ETAPHI_SCALE * std::sqrt(M_PI / 4)); - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - tk2mu[itk] = false; - } - for (int imu = 0, nmu = r.muon.size(); imu < nmu; ++imu) { - const auto &mu = r.muon[imu]; - if (debug_) - dbgPrintf( - "PFAlgo3 \t muon %3d (pt %7.2f, eta %+5.2f, phi %+5.2f) \n", imu, mu.floatPt(), mu.floatEta(), mu.floatPhi()); - float minDistance = 9e9; - switch (muMatchMode_) { - case MuMatchMode::BoxBestByPtRatio: - minDistance = 4.; - break; - case MuMatchMode::DrBestByPtRatio: - minDistance = 4.; - break; - case MuMatchMode::DrBestByPtDiff: - minDistance = 0.5 * mu.floatPt(); - break; - } - int imatch = -1; - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - int deta = std::abs(mu.hwEta - tk.hwEta); - int dphi = std::abs((mu.hwPhi - tk.hwPhi) % CaloCluster::PHI_WRAP); - float dr = floatDR(mu, tk); - float dpt = std::abs(mu.floatPt() - tk.floatPt()); - float dptr = (mu.hwPt > tk.hwPt ? mu.floatPt() / tk.floatPt() : tk.floatPt() / mu.floatPt()); - bool ok = false; - float distance = 9e9; - switch (muMatchMode_) { - case MuMatchMode::BoxBestByPtRatio: - ok = (deta < intDrMuonMatchBox) && (dphi < intDrMuonMatchBox); - distance = dptr; - break; - case MuMatchMode::DrBestByPtRatio: - ok = (dr < drMatchMu_); - distance = dptr; - break; - case MuMatchMode::DrBestByPtDiff: - ok = (dr < drMatchMu_); - distance = dpt; - break; - } - if (debug_ && dr < 0.4) { - dbgPrintf( - "PFAlgo3 \t\t possible match with track %3d (pt %7.2f, caloeta %+5.2f, calophi %+5.2f, dr %.2f, eta " - "%+5.2f, phi %+5.2f, dr %.2f): angular %1d, distance %.3f (vs %.3f)\n", - itk, - tk.floatPt(), - tk.floatEta(), - tk.floatPhi(), - dr, - tk.floatVtxEta(), - tk.floatVtxPhi(), - deltaR(mu.floatEta(), mu.floatPhi(), tk.floatVtxEta(), tk.floatVtxPhi()), - (ok ? 1 : 0), - distance, - minDistance); - } - if (!ok) - continue; - // FIXME for the moment, we do the floating point matching in pt - if (distance < minDistance) { - minDistance = distance; - imatch = itk; - } - } - if (debug_ && imatch > -1) - dbgPrintf("PFAlgo3 \t muon %3d (pt %7.2f) linked to track %3d (pt %7.2f)\n", - imu, - mu.floatPt(), - imatch, - r.track[imatch].floatPt()); - if (debug_ && imatch == -1) - dbgPrintf("PFAlgo3 \t muon %3d (pt %7.2f) not linked to any track\n", imu, mu.floatPt()); - mu2tk[imu] = imatch; - if (imatch > -1) { - tk2mu[imatch] = imu; - r.track[imatch].muonLink = true; - } - } -} - -void PFAlgo3::link_tk2em(Region &r, std::vector &tk2em) const { - // match all tracks to the closest EM cluster - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - //if (tk.muonLink) continue; // not necessary I think - float drbest = drMatchEm_; - for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) { - const auto &em = r.emcalo[iem]; - float dr = floatDR(tk, em); - if (dr < drbest) { - tk2em[itk] = iem; - drbest = dr; - } - } - if (debug_ && tk2em[itk] != -1) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) matches to EM %3d (pt %7.2f) with dr %.3f\n", - itk, - tk.floatPt(), - tk2em[itk], - tk2em[itk] == -1 ? 0.0 : r.emcalo[tk2em[itk]].floatPt(), - drbest); - } -} - -void PFAlgo3::link_em2calo(Region &r, std::vector &em2calo) const { - // match all em to the closest had (can happen in parallel to the above) - for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) { - const auto &em = r.emcalo[iem]; - float drbest = drMatchEmHad_; - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - const auto &calo = r.calo[ic]; - if (calo.floatEmPt() < ptMinFracMatchEm_ * em.floatPt()) - continue; - float dr = floatDR(calo, em); - if (dr < drbest) { - em2calo[iem] = ic; - drbest = dr; - } - } - if (debug_ && em2calo[iem] != -1) - dbgPrintf("PFAlgo3 \t EM %3d (pt %7.2f) matches to calo %3d (pt %7.2f, empt %7.2f) with dr %.3f\n", - iem, - em.floatPt(), - em2calo[iem], - em2calo[iem] == -1 ? 0.0 : r.calo[em2calo[iem]].floatPt(), - em2calo[iem] == -1 ? 0.0 : r.calo[em2calo[iem]].floatEmPt(), - drbest); - } -} - -void PFAlgo3::sum_tk2em(Region &r, - const std::vector &tk2em, - std::vector &em2ntk, - std::vector &em2sumtkpt, - std::vector &em2sumtkpterr) const { - // for each EM cluster, count and add up the pt of all the corresponding tracks (skipping muons) - for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) { - const auto &em = r.emcalo[iem]; - if (r.globalAbsEta(em.floatEta()) > 2.5) - continue; - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - if (tk2em[itk] == iem) { - const auto &tk = r.track[itk]; - if (tk.muonLink) - continue; - em2ntk[iem]++; - em2sumtkpt[iem] += tk.floatPt(); - em2sumtkpterr[iem] += tk.floatPtErr(); - } - } - } -} - -void PFAlgo3::emcalo_algo(Region &r, - const std::vector &em2ntk, - const std::vector &em2sumtkpt, - const std::vector &em2sumtkpterr) const { - // process ecal clusters after linking - for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) { - auto &em = r.emcalo[iem]; - em.isEM = false; - em.used = false; - em.hwFlags = 0; - if (r.globalAbsEta(em.floatEta()) > 2.5) - continue; - if (debug_) - dbgPrintf("PFAlgo3 \t EM %3d (pt %7.2f) has %2d tracks (sumpt %7.2f, sumpterr %7.2f), ptdif %7.2f +- %7.2f\n", - iem, - em.floatPt(), - em2ntk[iem], - em2sumtkpt[iem], - em2sumtkpterr[iem], - em.floatPt() - em2sumtkpt[iem], - std::max(em2sumtkpterr[iem], em.floatPtErr())); - if (em2ntk[iem] == 0) { // Photon - em.isEM = true; - addCaloToPF(r, em); - em.used = true; - if (debug_) - dbgPrintf("PFAlgo3 \t EM %3d (pt %7.2f) ---> promoted to photon\n", iem, em.floatPt()); - continue; - } - float ptdiff = em.floatPt() - em2sumtkpt[iem]; - float pterr = trackEmUseAlsoTrackSigma_ ? std::max(em2sumtkpterr[iem], em.floatPtErr()) : em.floatPtErr(); - // avoid "pt = inf +- inf" track to become an electron. - if (pterr > 2 * em.floatPt()) { - pterr = 2 * em.floatPt(); - if (debug_) - dbgPrintf("PFAlgo3 \t EM %3d (pt %7.2f) ---> clamp pterr ---> new ptdiff %7.2f +- %7.2f\n", - iem, - em.floatPt(), - ptdiff, - pterr); - } - - if (ptdiff > -ptMatchLow_ * pterr) { - em.isEM = true; - em.used = true; - // convert leftover to a photon if significant - if (ptdiff > +ptMatchHigh_ * pterr) { - auto &p = addCaloToPF(r, em); - p.setFloatPt(ptdiff); - if (debug_) - dbgPrintf("PFAlgo3 \t EM %3d (pt %7.2f) ---> promoted to electron(s) + photon (pt %7.2f)\n", - iem, - em.floatPt(), - ptdiff); - } else { - em.hwFlags = 1; // may use calo momentum - if (debug_) - dbgPrintf("PFAlgo3 \t EM %3d (pt %7.2f) ---> promoted to electron(s)\n", iem, em.floatPt()); - } - } else { - em.isEM = false; - em.used = false; - em.hwFlags = 0; - //discardCalo(r, em, 2); - } - } -} - -void PFAlgo3::emtk_algo(Region &r, - const std::vector &tk2em, - const std::vector &em2ntk, - const std::vector &em2sumtkpterr) const { - // promote all flagged tracks to electrons - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - auto &tk = r.track[itk]; - if (tk2em[itk] == -1 || tk.muonLink) - continue; - const auto &em = r.emcalo[tk2em[itk]]; - if (em.isEM) { - auto &p = addTrackToPF(r, tk); - p.cluster.src = em.src; - // FIXME to check if this is useful - if (trackEmMayUseCaloMomenta_ && em2ntk[tk2em[itk]] == 1 && em.hwFlags == 1) { - if (em.floatPtErr() < em2sumtkpterr[tk2em[itk]]) { - p.setFloatPt(em.floatPt()); - } - } - if (debug_) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) matched to EM %3d (pt %7.2f) promoted to electron with pt %7.2f\n", - itk, - tk.floatPt(), - tk2em[itk], - em.floatPt(), - p.floatPt()); - p.hwId = l1t::PFCandidate::Electron; - tk.used = true; - } - } -} - -void PFAlgo3::sub_em2calo(Region &r, const std::vector &em2calo) const { - // subtract EM component from Calo clusters for all photons and electrons (within tracker coverage) - // kill clusters that end up below their own uncertainty, or that loose 90% of the energy, - // unless they still have live EM clusters pointing to them - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - float pt0 = calo.floatPt(), ept0 = calo.floatEmPt(), pt = pt0, ept = ept0; - bool keepme = false; - for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) { - if (em2calo[iem] == ic) { - const auto &em = r.emcalo[iem]; - if (em.isEM) { - if (debug_) - dbgPrintf( - "PFAlgo3 \t EM %3d (pt %7.2f) is subtracted from calo %3d (pt %7.2f) scaled by %.3f (deltaPt = " - "%7.2f)\n", - iem, - em.floatPt(), - ic, - calo.floatPt(), - emHadSubtractionPtSlope_, - emHadSubtractionPtSlope_ * em.floatPt()); - pt -= emHadSubtractionPtSlope_ * em.floatPt(); - ept -= em.floatPt(); - } else { - keepme = true; - if (debug_) - dbgPrintf( - "PFAlgo3 \t EM %3d (pt %7.2f) not subtracted from calo %3d (pt %7.2f), and calo marked to be kept " - "after EM subtraction\n", - iem, - em.floatPt(), - ic, - calo.floatPt()); - } - } - } - if (pt < pt0) { - if (debug_) - dbgPrintf( - "PFAlgo3 \t calo %3d (pt %7.2f +- %7.2f) has a subtracted pt of %7.2f, empt %7.2f -> %7.2f, isem %d\n", - ic, - calo.floatPt(), - calo.floatPtErr(), - pt, - ept0, - ept, - calo.isEM); - calo.setFloatPt(pt); - calo.setFloatEmPt(ept); - if (!keepme && - ((emCaloUseAlsoCaloSigma_ ? pt < calo.floatPtErr() : false) || pt <= 0.125 * pt0 || - (calo.isEM && ept <= 0.125 * ept0))) { // the <= is important since in firmware the pt0/8 can be zero - if (debug_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f) ----> discarded\n", ic, calo.floatPt()); - calo.used = true; - calo.setFloatPt(pt0); //discardCalo(r, calo, 1); // log this as discarded, for debugging - } - } - } -} - -void PFAlgo3::link_tk2calo(Region &r, std::vector &tk2calo) const { - // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later) - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - if (tk.muonLink || tk.used) - continue; // not necessary but just a waste of CPU otherwise - float drbest = drMatch_, dptscale = 0; - switch (tkCaloLinkMetric_) { - case TkCaloLinkMetric::BestByDR: - drbest = drMatch_; - break; - case TkCaloLinkMetric::BestByDRPt: - drbest = 1.0; - dptscale = drMatch_ / tk.floatCaloPtErr(); - break; - case TkCaloLinkMetric::BestByDR2Pt2: - drbest = 1.0; - dptscale = drMatch_ / tk.floatCaloPtErr(); - break; - } - float minCaloPt = tk.floatPt() - ptMatchLow_ * tk.floatCaloPtErr(); - if (debug_) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) to be matched to calo, min pT %7.2f\n", itk, tk.floatPt(), minCaloPt); - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo.used || calo.floatPt() <= minCaloPt) - continue; - float dr = floatDR(tk, calo), dq; - switch (tkCaloLinkMetric_) { - case TkCaloLinkMetric::BestByDR: - if (dr < drbest) { - tk2calo[itk] = ic; - drbest = dr; - } - break; - case TkCaloLinkMetric::BestByDRPt: - dq = dr + std::max(tk.floatPt() - calo.floatPt(), 0.) * dptscale; - //if (debug_ && dr < 0.2) dbgPrintf("PFAlgo3 \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n", itk, tk.floatPt(), ic, calo.floatPt(), dr, dq); - if (dr < drMatch_ && dq < drbest) { - tk2calo[itk] = ic; - drbest = dq; - } - break; - case TkCaloLinkMetric::BestByDR2Pt2: - dq = hypot(dr, std::max(tk.floatPt() - calo.floatPt(), 0.) * dptscale); - //if (debug_ && dr < 0.2) dbgPrintf("PFAlgo3 \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n", itk, tk.floatPt(), ic, calo.floatPt(), dr, dq); - if (dr < drMatch_ && dq < drbest) { - tk2calo[itk] = ic; - drbest = dq; - } - break; - } - } - if (debug_ && tk2calo[itk] != -1) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) matches to calo %3d (pt %7.2f) with dist %.3f\n", - itk, - tk.floatPt(), - tk2calo[itk], - tk2calo[itk] == -1 ? 0.0 : r.calo[tk2calo[itk]].floatPt(), - drbest); - // now we re-do this for debugging sake, it may be done for real later - if (debug_ && tk2calo[itk] == -1) { - int ibest = -1; - drbest = 0.3; - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo.used) - continue; - float dr = floatDR(tk, calo); - if (dr < drbest) { - ibest = ic; - drbest = dr; - } - } - if (ibest != -1) - dbgPrintf( - "PFAlgo3 \t track %3d (pt %7.2f) would match to calo %3d (pt %7.2f) with dr %.3f if the pt min and dr " - "requirement had been relaxed\n", - itk, - tk.floatPt(), - ibest, - r.calo[ibest].floatPt(), - drbest); - } - } -} - -void PFAlgo3::sum_tk2calo(Region &r, - const std::vector &tk2calo, - std::vector &calo2ntk, - std::vector &calo2sumtkpt, - std::vector &calo2sumtkpterr) const { - // for each calo, compute the sum of the track pt - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - const auto &calo = r.calo[ic]; - if (r.globalAbsEta(calo.floatEta()) > 2.5) - continue; - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - if (tk2calo[itk] == ic) { - const auto &tk = r.track[itk]; - if (tk.muonLink || tk.used) - continue; - calo2ntk[ic]++; - calo2sumtkpt[ic] += tk.floatPt(); - calo2sumtkpterr[ic] += std::pow(tk.floatCaloPtErr(), sumTkCaloErr2_ ? 2 : 1); - } - } - if (sumTkCaloErr2_ && calo2sumtkpterr[ic] > 0) - calo2sumtkpterr[ic] = std::sqrt(calo2sumtkpterr[ic]); - } -} - -void PFAlgo3::unlinkedtk_algo(Region &r, const std::vector &tk2calo) const { - // in the meantime, promote unlinked low pt tracks to hadrons - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - auto &tk = r.track[itk]; - if (tk2calo[itk] != -1 || tk.muonLink || tk.used || !tk.quality(l1tpf_impl::InputTrack::PFLOOSE)) - continue; - float maxPt = tk.quality(l1tpf_impl::InputTrack::PFTIGHT) ? tightTrackMaxInvisiblePt_ : maxInvisiblePt_; - if (tk.floatPt() < maxPt) { - if (debug_) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) not matched to calo, kept as charged hadron\n", itk, tk.floatPt()); - auto &p = addTrackToPF(r, tk); - p.hwStatus = GoodTK_NoCalo; - tk.used = true; - } else { - if (debug_) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) not matched to calo, dropped\n", itk, tk.floatPt()); - //discardTrack(r, tk, BadTK_NoCalo); // log this as discarded, for debugging - } - } -} - -void PFAlgo3::calo_relink(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr) const { - /// OPTIONAL STEP: try to recover split hadron showers (v1.0): - // take hadrons that are not track matched, close by a hadron which has an excess of track pt vs calo pt - // add this pt to the calo pt of the other cluster - // off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events - std::vector addtopt(r.calo.size(), 0); - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo2ntk[ic] != 0 || calo.used || r.globalAbsEta(calo.floatEta()) > 2.5) - continue; - int i2best = -1; - float drbest = caloReLinkDr_; - for (int ic2 = 0; ic2 < nc; ++ic2) { - const auto &calo2 = r.calo[ic2]; - if (calo2ntk[ic2] == 0 || calo2.used || r.globalAbsEta(calo2.floatEta()) > 2.5) - continue; - float dr = floatDR(calo, calo2); - //// uncomment below for more verbose debugging - //if (debug_ && dr < 0.5) dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f) with no tracks is at dr %.3f from calo %3d with pt %7.2f (sum tk pt %7.2f), track excess %7.2f +- %7.2f\n", ic, calo.floatPt(), dr, ic2, calo2.floatPt(), calo2sumtkpt[ic2], calo2sumtkpt[ic2] - calo2.floatPt(), useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr()); - if (dr < drbest) { - float ptdiff = - calo2sumtkpt[ic2] - calo2.floatPt() + (useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr()); - if (ptdiff >= caloReLinkThreshold_ * calo.floatPt()) { - i2best = ic2; - drbest = dr; - } - } - } - if (i2best != -1) { - const auto &calo2 = r.calo[i2best]; - if (debug_) - dbgPrintf( - "PFAlgo3 \t calo %3d (pt %7.2f) with no tracks matched within dr %.3f with calo %3d with pt %7.2f (sum tk " - "pt %7.2f), track excess %7.2f +- %7.2f\n", - ic, - calo.floatPt(), - drbest, - i2best, - calo2.floatPt(), - calo2sumtkpt[i2best], - calo2sumtkpt[i2best] - calo2.floatPt(), - useTrackCaloSigma_ ? calo2sumtkpterr[i2best] : calo2.floatPtErr()); - calo.used = true; - addtopt[i2best] += calo.floatPt(); - } - } - // we do this at the end, so that the above loop is parallelizable - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - if (addtopt[ic]) { - auto &calo = r.calo[ic]; - if (debug_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f, sum tk pt %7.2f) is increased to pt %7.2f after merging\n", - ic, - calo.floatPt(), - calo2sumtkpt[ic], - calo.floatPt() + addtopt[ic]); - calo.setFloatPt(calo.floatPt() + addtopt[ic]); - } - } -} - -void PFAlgo3::linkedcalo_algo(Region &r, - const std::vector &calo2ntk, - const std::vector &calo2sumtkpt, - const std::vector &calo2sumtkpterr, - std::vector &calo2alpha) const { - /// ------------- next step (needs the previous) ---------------- - // process matched calo clusters, compare energy to sum track pt - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - auto &calo = r.calo[ic]; - if (calo2ntk[ic] == 0 || calo.used) - continue; - float ptdiff = calo.floatPt() - calo2sumtkpt[ic]; - float pterr = useTrackCaloSigma_ ? calo2sumtkpterr[ic] : calo.floatPtErr(); - if (debug_) - dbgPrintf( - "PFAlgo3 \t calo %3d (pt %7.2f +- %7.2f, empt %7.2f) has %2d tracks (sumpt %7.2f, sumpterr %7.2f), ptdif " - "%7.2f +- %7.2f\n", - ic, - calo.floatPt(), - calo.floatPtErr(), - calo.floatEmPt(), - calo2ntk[ic], - calo2sumtkpt[ic], - calo2sumtkpterr[ic], - ptdiff, - pterr); - if (ptdiff > +ptMatchHigh_ * pterr) { - if (ecalPriority_) { - if (calo.floatEmPt() > 1) { - float emptdiff = std::min(ptdiff, calo.floatEmPt()); - if (debug_) - dbgPrintf( - "PFAlgo3 \t calo %3d (pt %7.2f, empt %7.2f) ---> make photon with pt %7.2f, reduce ptdiff to %7.2f " - "+- %7.2f\n", - ic, - calo.floatPt(), - calo.floatEmPt(), - emptdiff, - ptdiff - emptdiff, - pterr); - auto &p = addCaloToPF(r, calo); - p.setFloatPt(emptdiff); - p.hwId = l1t::PFCandidate::Photon; - ptdiff -= emptdiff; - } - if (ptdiff > 2) { - if (debug_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f, empt %7.2f) ---> make also neutral hadron with pt %7.2f\n", - ic, - calo.floatPt(), - calo.floatEmPt(), - ptdiff); - auto &p = addCaloToPF(r, calo); - p.setFloatPt(ptdiff); - p.hwId = l1t::PFCandidate::NeutralHadron; - } - } else { - if (debug_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f) ---> promoted to neutral with pt %7.2f\n", - ic, - calo.floatPt(), - ptdiff); - auto &p = addCaloToPF(r, calo); - p.setFloatPt(ptdiff); - calo.hwFlags = 0; - } - } else if (ptdiff > -ptMatchLow_ * pterr) { - // nothing to do (weighted average happens when we process the tracks) - calo.hwFlags = 1; - if (debug_) - dbgPrintf( - "PFAlgo3 \t calo %3d (pt %7.2f) ---> to be deleted, will use tracks instead\n", ic, calo.floatPt()); - //discardCalo(r, calo, 0); // log this as discarded, for debugging - } else { - // tracks overshoot, rescale to tracks to calo - calo2alpha[ic] = rescaleTracks_ ? calo.floatPt() / calo2sumtkpt[ic] : 1.0; - calo.hwFlags = 2; - if (debug_ && rescaleTracks_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f) ---> tracks overshoot and will be scaled down by %.4f\n", - ic, - calo.floatPt(), - calo2alpha[ic]); - if (debug_ && !rescaleTracks_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f) ---> tracks overshoot by %.4f\n", - ic, - calo.floatPt(), - calo2sumtkpt[ic] / calo.floatPt()); - } - calo.used = true; - } -} - -void PFAlgo3::linkedtk_algo(Region &r, - const std::vector &tk2calo, - const std::vector &calo2ntk, - const std::vector &calo2alpha) const { - // process matched tracks, if necessary rescale or average - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - auto &tk = r.track[itk]; - if (tk2calo[itk] == -1 || tk.muonLink || tk.used) - continue; - auto &p = addTrackToPF(r, tk); - tk.used = true; - const auto &calo = r.calo[tk2calo[itk]]; - p.cluster.src = calo.src; - if (calo.hwFlags == 1) { - // can do weighted average if there's just one track - if (calo2ntk[tk2calo[itk]] == 1 && caloTrkWeightedAverage_) { - p.hwStatus = GoodTK_Calo_TkPt; - float ptavg = tk.floatPt(); - if (tk.floatPtErr() > 0) { - float wcalo = 1.0 / std::pow(tk.floatCaloPtErr(), 2); - float wtk = 1.0 / std::pow(tk.floatPtErr(), 2); - ptavg = (calo.floatPt() * wcalo + tk.floatPt() * wtk) / (wcalo + wtk); - p.hwStatus = GoodTK_Calo_TkCaloPt; - } - p.setFloatPt(ptavg); - if (debug_) - dbgPrintf( - "PFAlgo3 \t track %3d (pt %7.2f +- %7.2f) combined with calo %3d (pt %7.2f +- %7.2f (from tk) yielding " - "candidate of pt %7.2f\n", - itk, - tk.floatPt(), - tk.floatPtErr(), - tk2calo[itk], - calo.floatPt(), - tk.floatCaloPtErr(), - ptavg); - } else { - p.hwStatus = GoodTK_Calo_TkPt; - if (debug_) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) linked to calo %3d promoted to charged hadron\n", - itk, - tk.floatPt(), - tk2calo[itk]); - } - } else if (calo.hwFlags == 2) { - // must rescale - p.setFloatPt(tk.floatPt() * calo2alpha[tk2calo[itk]]); - p.hwStatus = GoodTk_Calo_CaloPt; - if (debug_) - dbgPrintf( - "PFAlgo3 \t track %3d (pt %7.2f, stubs %2d chi2 %7.1f) linked to calo %3d promoted to charged hadron with " - "pt %7.2f after maybe rescaling\n", - itk, - tk.floatPt(), - int(tk.hwStubs), - tk.hwChi2 * 0.1f, - tk2calo[itk], - p.floatPt()); - } - } -} - -void PFAlgo3::unlinkedcalo_algo(Region &r) const { - // process unmatched calo clusters - for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) { - if (!r.calo[ic].used) { - addCaloToPF(r, r.calo[ic]); - if (debug_) - dbgPrintf("PFAlgo3 \t calo %3d (pt %7.2f) not linked, promoted to neutral\n", ic, r.calo[ic].floatPt()); - } - } -} - -void PFAlgo3::save_muons(Region &r, const std::vector &tk2mu) const { - // finally do muons - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - if (r.track[itk].muonLink) { - auto &p = addTrackToPF(r, r.track[itk]); - p.muonsrc = r.muon[tk2mu[itk]].src; - if (debug_) - dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) promoted to muon.\n", itk, r.track[itk].floatPt()); - } - } -} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/PFAlgoBase.cc b/L1Trigger/Phase2L1ParticleFlow/src/PFAlgoBase.cc deleted file mode 100644 index 35b8bd6d20a3c..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/PFAlgoBase.cc +++ /dev/null @@ -1,57 +0,0 @@ -#include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgoBase.h" - -#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" - -using namespace l1tpf_impl; - -PFAlgoBase::PFAlgoBase(const edm::ParameterSet &iConfig) : debug_(iConfig.getUntrackedParameter("debug", 0)) {} - -PFAlgoBase::~PFAlgoBase() {} - -void PFAlgoBase::initRegion(Region &r, bool doSort) const { - r.inputCrop(doSort); - r.pf.clear(); - r.puppi.clear(); - for (auto &c : r.calo) - c.used = false; - for (auto &c : r.emcalo) - c.used = false; - for (auto &t : r.track) { - t.used = false; - t.muonLink = false; - } -} - -PFParticle &PFAlgoBase::addTrackToPF(std::vector &pfs, const PropagatedTrack &tk) const { - PFParticle pf; - pf.hwPt = tk.hwPt; - pf.hwEta = tk.hwEta; - pf.hwPhi = tk.hwPhi; - pf.hwVtxEta = tk.hwEta; // FIXME: get from the track - pf.hwVtxPhi = tk.hwPhi; // before propagation - pf.track = tk; - pf.cluster.hwPt = 0; - pf.cluster.src = nullptr; - pf.muonsrc = nullptr; - pf.hwId = (tk.muonLink ? l1t::PFCandidate::Muon : l1t::PFCandidate::ChargedHadron); - pf.hwStatus = 0; - pfs.push_back(pf); - return pfs.back(); -} - -PFParticle &PFAlgoBase::addCaloToPF(std::vector &pfs, const CaloCluster &calo) const { - PFParticle pf; - pf.hwPt = calo.hwPt; - pf.hwEta = calo.hwEta; - pf.hwPhi = calo.hwPhi; - pf.hwVtxEta = calo.hwEta; - pf.hwVtxPhi = calo.hwPhi; - pf.track.hwPt = 0; - pf.track.src = nullptr; - pf.cluster = calo; - pf.muonsrc = nullptr; - pf.hwId = (calo.isEM ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron); - pf.hwStatus = 0; - pfs.push_back(pf); - return pfs.back(); -} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/PFTkEgAlgo.cc b/L1Trigger/Phase2L1ParticleFlow/src/PFTkEgAlgo.cc deleted file mode 100644 index e71c8aa762e96..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/PFTkEgAlgo.cc +++ /dev/null @@ -1,261 +0,0 @@ -#include "L1Trigger/Phase2L1ParticleFlow/interface/PFTkEGAlgo.h" - -using namespace l1tpf_impl; - -#include "DataFormats/Math/interface/deltaPhi.h" -#include "DataFormats/L1TParticleFlow/interface/PFTrack.h" -#include "DataFormats/Common/interface/RefToPtr.h" -#include -#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" - -namespace { - template - float floatDR(const T1 &t1, const T2 &t2) { - return deltaR(t1.floatEta(), t1.floatPhi(), t2.floatEta(), t2.floatPhi()); - } -} // namespace - -PFTkEGAlgo::PFTkEGAlgo(const edm::ParameterSet &pset) - : debug_(pset.getUntrackedParameter("debug")), - doBremRecovery_(pset.getParameter("doBremRecovery")), - doTkIsolation_(pset.getParameter("doTkIsolation")), - filterHwQuality_(pset.getParameter("filterHwQuality")), - caloHwQual_(pset.getParameter("caloHwQual")), - dEtaMaxBrem_(pset.getParameter("dEtaMaxBrem")), - dPhiMaxBrem_(pset.getParameter("dPhiMaxBrem")), - absEtaBoundaries_(pset.getParameter>("absEtaBoundaries")), - dEtaValues_(pset.getParameter>("dEtaValues")), - dPhiValues_(pset.getParameter>("dPhiValues")), - caloEtMin_(pset.getParameter("caloEtMin")), - trkQualityPtMin_(pset.getParameter("trkQualityPtMin")), - trkQualityChi2_(pset.getParameter("trkQualityChi2")), - writeEgSta_(pset.getParameter("writeEgSta")), - tkIsoParametersTkEm_(pset.getParameter("tkIsoParametersTkEm")), - tkIsoParametersTkEle_(pset.getParameter("tkIsoParametersTkEle")), - pfIsoParametersTkEm_(pset.getParameter("pfIsoParametersTkEm")), - pfIsoParametersTkEle_(pset.getParameter("pfIsoParametersTkEle")) {} - -PFTkEGAlgo::~PFTkEGAlgo() {} - -PFTkEGAlgo::IsoParameters::IsoParameters(const edm::ParameterSet &pset) - : tkQualityPtMin(pset.getParameter("tkQualityPtMin")), - dZ(pset.getParameter("dZ")), - dRMin(pset.getParameter("dRMin")), - dRMax(pset.getParameter("dRMax")), - tkQualityChi2Max(pset.getParameter("tkQualityChi2Max")), - dRMin2(dRMin * dRMin), - dRMax2(dRMax * dRMax) {} - -void PFTkEGAlgo::initRegion(Region &r) const { - // NOTE: assume imput is sorted already - r.egphotons.clear(); - r.egeles.clear(); -} - -void PFTkEGAlgo::runTkEG(Region &r) const { - initRegion(r); - - if (debug_ > 0) { - edm::LogInfo("PFTkEGAlgo") << " # emCalo: " << r.emcalo.size() << " # tk: " << r.track.size(); - } - if (debug_ > 0) { - for (int ic = 0, nc = r.emcalo.size(); ic < nc; ++ic) { - const auto &calo = r.emcalo[ic]; - edm::LogInfo("PFTkEGAlgo") << "[OLD] IN calo[" << ic << "] pt: " << calo.floatPt() - << " eta: " << r.globalEta(calo.floatEta()) << " phi: " << r.globalPhi(calo.floatPhi()) - << " hwEta: " << calo.hwEta << " hwPhi: " << calo.hwPhi - << " src pt: " << calo.src->pt() << " src eta: " << calo.src->eta() - << " src phi: " << calo.src->phi(); - } - } - - // NOTE: we run this step for all clusters (before matching) as it is done in the pre-PF EG algorithm - std::vector emCalo2emCalo(r.emcalo.size(), -1); - if (doBremRecovery_) - link_emCalo2emCalo(r, emCalo2emCalo); - - // track to EM calo matching - std::vector emCalo2tk(r.emcalo.size(), -1); - link_emCalo2tk(r, emCalo2tk); - - // add EG objects to region; - eg_algo(r, emCalo2emCalo, emCalo2tk); -} - -void PFTkEGAlgo::eg_algo(Region &r, const std::vector &emCalo2emCalo, const std::vector &emCalo2tk) const { - for (int ic = 0, nc = r.emcalo.size(); ic < nc; ++ic) { - auto &calo = r.emcalo[ic]; - if (filterHwQuality_ && calo.hwFlags != caloHwQual_) - continue; - - int itk = emCalo2tk[ic]; - - // 1. create EG objects before brem recovery - addEgObjsToPF(r, ic, calo.hwFlags, calo.floatPt(), itk); - - // check if brem recovery is on - if (!doBremRecovery_) - continue; - - // check if the cluster has already been used in a brem reclustering - if (emCalo2emCalo[ic] != -1) - continue; - - float ptBremReco = calo.floatPt(); - - for (int jc = ic; jc < nc; ++jc) { - if (emCalo2emCalo[jc] == ic) { - auto &otherCalo = r.emcalo[jc]; - ptBremReco += otherCalo.floatPt(); - } - } - - // 2. create EG objects with brem recovery - // FIXME: duplicating the object is suboptimal but this is done for keeping things as in TDR code... - addEgObjsToPF(r, ic, calo.hwFlags + 1, ptBremReco, itk); - } -} - -EGIsoParticle &PFTkEGAlgo::addEGIsoToPF(std::vector &egobjs, - const CaloCluster &calo, - const int hwQual, - const float ptCorr) const { - EGIsoParticle egiso; - egiso.setFloatPt(ptCorr); - egiso.hwEta = calo.hwEta; - egiso.hwPhi = calo.hwPhi; - egiso.cluster = calo; - - egiso.hwQual = hwQual; - - egiso.hwIso = 0; - egiso.hwIsoPV = 0; - egiso.hwPFIso = 0; - egiso.hwPFIsoPV = 0; - egiso.ele_idx = -1; - egobjs.push_back(egiso); - return egobjs.back(); -} - -EGIsoEleParticle &PFTkEGAlgo::addEGIsoEleToPF(std::vector &egobjs, - const CaloCluster &calo, - const PropagatedTrack &track, - const int hwQual, - const float ptCorr) const { - EGIsoEleParticle egiso; - egiso.setFloatPt(ptCorr); - egiso.hwEta = calo.hwEta; - egiso.hwPhi = calo.hwPhi; - egiso.cluster = calo; - - egiso.hwVtxEta = track.hwVtxEta; - egiso.hwVtxPhi = track.hwVtxPhi; - egiso.hwZ0 = track.hwZ0; - egiso.hwCharge = track.hwCharge; - egiso.track = track; - - egiso.hwQual = hwQual; - egiso.hwIso = 0; - egiso.hwPFIso = 0; - egobjs.push_back(egiso); - return egobjs.back(); -} - -void PFTkEGAlgo::addEgObjsToPF( - Region &r, const int calo_idx, const int hwQual, const float ptCorr, const int tk_idx) const { - EGIsoParticle &egobj = addEGIsoToPF(r.egphotons, r.emcalo[calo_idx], hwQual, ptCorr); - if (tk_idx != -1) { - egobj.ele_idx = r.egeles.size(); - addEGIsoEleToPF(r.egeles, r.emcalo[calo_idx], r.track[tk_idx], hwQual, ptCorr); - } -} - -void PFTkEGAlgo::link_emCalo2emCalo(const Region &r, std::vector &emCalo2emCalo) const { - // NOTE: we assume the input to be sorted!!! - for (int ic = 0, nc = r.emcalo.size(); ic < nc; ++ic) { - auto &calo = r.emcalo[ic]; - if (filterHwQuality_ && calo.hwFlags != caloHwQual_) - continue; - - if (emCalo2emCalo[ic] != -1) - continue; - - for (int jc = ic + 1; jc < nc; ++jc) { - if (emCalo2emCalo[jc] != -1) - continue; - - auto &otherCalo = r.emcalo[jc]; - if (filterHwQuality_ && otherCalo.hwFlags != caloHwQual_) - continue; - - if (fabs(otherCalo.floatEta() - calo.floatEta()) < dEtaMaxBrem_ && - fabs(deltaPhi(otherCalo.floatPhi(), calo.floatPhi())) < dPhiMaxBrem_) { - emCalo2emCalo[jc] = ic; - } - } - } -} - -void PFTkEGAlgo::link_emCalo2tk(const Region &r, std::vector &emCalo2tk) const { - for (int ic = 0, nc = r.emcalo.size(); ic < nc; ++ic) { - auto &calo = r.emcalo[ic]; - - if (filterHwQuality_ && calo.hwFlags != caloHwQual_) - continue; - // compute elliptic matching - auto eta_index = - std::distance( - absEtaBoundaries_.begin(), - std::lower_bound(absEtaBoundaries_.begin(), absEtaBoundaries_.end(), r.globalAbsEta(calo.floatEta()))) - - 1; - float dEtaMax = dEtaValues_[eta_index]; - float dPhiMax = dPhiValues_[eta_index]; - - if (debug_ > 4) - edm::LogInfo("PFTkEGAlgo") << "idx: " << eta_index << " deta: " << dEtaMax << " dphi: " << dPhiMax; - - float dPtMin = 999; - if (debug_ > 3) - edm::LogInfo("PFTkEGAlgo") << "--- calo: pt: " << calo.floatPt() << " eta: " << calo.floatEta() - << " phi: " << calo.floatPhi(); - for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) { - const auto &tk = r.track[itk]; - if (debug_ > 3) - edm::LogInfo("PFTkEGAlgo") << " - tk: pt: " << tk.floatPt() << " eta: " << tk.floatEta() - << " phi: " << tk.floatPhi(); - - if (tk.floatPt() < trkQualityPtMin_) - continue; - - float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi()); - float d_eta = tk.floatEta() - calo.floatEta(); // We only use it squared - - if (debug_ > 3) - edm::LogInfo("PFTkEGAlgo") << " deta: " << fabs(d_eta) << " dphi: " << d_phi << " ell: " - << ((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax)); - - if ((((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) < 1.) { - if (debug_ > 3) - edm::LogInfo("PFTkEGAlgo") << " pass elliptic "; - // NOTE: for now we implement only best pt match. This is NOT what is done in the L1TkElectronTrackProducer - if (fabs(tk.floatPt() - calo.floatPt()) < dPtMin) { - if (debug_ > 3) - edm::LogInfo("PFTkEGAlgo") << " best pt match: " << fabs(tk.floatPt() - calo.floatPt()); - emCalo2tk[ic] = itk; - dPtMin = fabs(tk.floatPt() - calo.floatPt()); - } - } - } - } -} - -void PFTkEGAlgo::runTkIso(Region &r, const float z0) const { - compute_isolation_tkEle(r, r.track, tkIsoParametersTkEle_, z0, false); - compute_isolation_tkEm(r, r.track, tkIsoParametersTkEm_, z0, false); -} - -void PFTkEGAlgo::runPFIso(Region &r, const float z0) const { - compute_isolation_tkEle(r, r.pf, pfIsoParametersTkEle_, z0, true); - compute_isolation_tkEm(r, r.pf, pfIsoParametersTkEm_, z0, true); -} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/Region.cc b/L1Trigger/Phase2L1ParticleFlow/src/Region.cc deleted file mode 100644 index d3e1f26a04b4e..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/Region.cc +++ /dev/null @@ -1,151 +0,0 @@ -#include "L1Trigger/Phase2L1ParticleFlow/interface/Region.h" -#include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" -#include -#include - -const char *l1tpf_impl::Region::inputTypeName(int type) { - switch (InputType(type)) { - case calo_type: - return "Calo"; - case emcalo_type: - return "EmCalo"; - case track_type: - return "TK"; - case l1mu_type: - return "Mu"; - case n_input_types: - throw cms::Exception( - "LogicError", "n_input_types is not a type to be used, but only a compile-time const for iterating on types"); - } - return "NO_SUCH_INPUT_TYPE"; -} -const char *l1tpf_impl::Region::outputTypeName(int type) { - switch (OutputType(type)) { - case any_type: - return ""; - case charged_type: - return "Charged"; - case neutral_type: - return "Neutral"; - case electron_type: - return "Electron"; - case pfmuon_type: - return "Muon"; - case charged_hadron_type: - return "ChargedHadron"; - case neutral_hadron_type: - return "NeutralHadron"; - case photon_type: - return "Photon"; - case n_output_types: - throw cms::Exception( - "LogicError", - "n_output_types is not a type to be used, but only a compile-time const for iterating on types"); - } - return "NO_SUCH_OUTPUT_TYPE"; -} - -unsigned int l1tpf_impl::Region::nInput(InputType type) const { - switch (type) { - case calo_type: - return calo.size(); - case emcalo_type: - return emcalo.size(); - case track_type: - return track.size(); - case l1mu_type: - return muon.size(); - case n_input_types: - throw cms::Exception( - "LogicError", "n_input_types is not a type to be used, but only a compile-time const for iterating on types"); - } - return 9999; -} - -unsigned int l1tpf_impl::Region::nOutput(OutputType type, bool usePuppi, bool fiducial) const { - unsigned int ret = 0; - for (const auto &p : (usePuppi ? puppi : pf)) { - if (p.hwPt <= 0) - continue; - if (fiducial && !fiducialLocal(p.floatEta(), p.floatPhi())) - continue; - switch (type) { - case any_type: - ret++; - break; - case charged_type: - if (p.intCharge() != 0) - ret++; - break; - case neutral_type: - if (p.intCharge() == 0) - ret++; - break; - case electron_type: - if (p.hwId == l1t::PFCandidate::Electron) - ret++; - break; - case pfmuon_type: - if (p.hwId == l1t::PFCandidate::Muon) - ret++; - break; - case charged_hadron_type: - if (p.hwId == l1t::PFCandidate::ChargedHadron) - ret++; - break; - case neutral_hadron_type: - if (p.hwId == l1t::PFCandidate::NeutralHadron) - ret++; - break; - case photon_type: - if (p.hwId == l1t::PFCandidate::Photon) - ret++; - break; - case n_output_types: - throw cms::Exception( - "LogicError", - "n_output_types is not a type to be used, but only a compile-time const for iterating on types"); - } - } - return ret; -} - -void l1tpf_impl::Region::inputCrop(bool doSort) { - if (doSort) { - std::sort(calo.begin(), calo.end()); - std::sort(emcalo.begin(), emcalo.end()); - std::sort(track.begin(), track.end()); - std::sort(muon.begin(), muon.end()); - } - if (ncaloMax > 0 && calo.size() > ncaloMax) { - caloOverflow = calo.size() - ncaloMax; - calo.resize(ncaloMax); - } - if (nemcaloMax > 0 && emcalo.size() > nemcaloMax) { - emcaloOverflow = emcalo.size() - nemcaloMax; - emcalo.resize(nemcaloMax); - } - if (ntrackMax > 0 && track.size() > ntrackMax) { - trackOverflow = track.size() - ntrackMax; - track.resize(ntrackMax); - } - if (nmuonMax > 0 && muon.size() > nmuonMax) { - muonOverflow = muon.size() - nmuonMax; - muon.resize(nmuonMax); - } -} - -void l1tpf_impl::Region::outputCrop(bool doSort) { - if (doSort) { - std::sort(puppi.begin(), puppi.end()); - std::sort(pf.begin(), pf.end()); - } - if (npuppiMax > 0 && puppi.size() > npuppiMax) { - puppiOverflow = puppi.size() - npuppiMax; - puppi.resize(npuppiMax); - } - if (npfMax > 0 && pf.size() > npfMax) { - pfOverflow = pf.size() - npfMax; - pf.resize(npfMax); - } -} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/RegionMapper.cc b/L1Trigger/Phase2L1ParticleFlow/src/RegionMapper.cc deleted file mode 100644 index 0134067cc20b4..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/RegionMapper.cc +++ /dev/null @@ -1,405 +0,0 @@ -#include "L1Trigger/Phase2L1ParticleFlow/interface/RegionMapper.h" -#include "FWCore/Framework/interface/Event.h" -#include "DataFormats/Common/interface/RefToPtr.h" - -#include "DataFormats/L1TCorrelator/interface/TkElectron.h" -#include "DataFormats/L1TCorrelator/interface/TkElectronFwd.h" -#include "DataFormats/L1Trigger/interface/EGamma.h" -#include "DataFormats/L1TCorrelator/interface/TkEm.h" -#include "DataFormats/L1TCorrelator/interface/TkEmFwd.h" - -using namespace l1tpf_impl; - -RegionMapper::RegionMapper(const edm::ParameterSet &iConfig) : useRelativeRegionalCoordinates_(false) { - if (iConfig.existsAs>("regions")) { - useRelativeRegionalCoordinates_ = iConfig.getParameter("useRelativeRegionalCoordinates"); - for (const edm::ParameterSet &preg : iConfig.getParameter>("regions")) { - std::vector etaBoundaries = preg.getParameter>("etaBoundaries"); - unsigned int phiSlices = preg.getParameter("phiSlices"); - float etaExtra = preg.getParameter("etaExtra"); - float phiExtra = preg.getParameter("phiExtra"); - float phiWidth = 2 * M_PI / phiSlices; - unsigned int ncalomax = 0, nemcalomax = 0, ntrackmax = 0, nmuonmax = 0, npfmax = 0, npuppimax = 0; - if (preg.existsAs("caloNMax")) - ncalomax = preg.getParameter("caloNMax"); - if (preg.existsAs("emcaloNMax")) - nemcalomax = preg.getParameter("emcaloNMax"); - if (preg.existsAs("trackNMax")) - ntrackmax = preg.getParameter("trackNMax"); - if (preg.existsAs("muonNMax")) - nmuonmax = preg.getParameter("muonNMax"); - if (preg.existsAs("pfNMax")) - npfmax = preg.getParameter("pfNMax"); - if (preg.existsAs("puppiNMax")) - npuppimax = preg.getParameter("puppiNMax"); - for (unsigned int ieta = 0, neta = etaBoundaries.size() - 1; ieta < neta; ++ieta) { - for (unsigned int iphi = 0; iphi < phiSlices; ++iphi) { - //float phiCenter = (iphi + 0.5) * phiWidth - M_PI; - float phiCenter = reco::reduceRange(iphi * phiWidth); //align with L1 TrackFinder phi sector indexing for now - regions_.push_back(Region(etaBoundaries[ieta], - etaBoundaries[ieta + 1], - phiCenter, - phiWidth, - etaExtra, - phiExtra, - useRelativeRegionalCoordinates_, - ncalomax, - nemcalomax, - ntrackmax, - nmuonmax, - npfmax, - npuppimax)); - } - } - } - std::string trackRegionMode = "TrackAssoMode::any"; - if (iConfig.existsAs("trackRegionMode")) - trackRegionMode = iConfig.getParameter("trackRegionMode"); - if (trackRegionMode == "atVertex") - trackRegionMode_ = TrackAssoMode::atVertex; - else if (trackRegionMode == "atCalo") - trackRegionMode_ = TrackAssoMode::atCalo; - else if (trackRegionMode == "any") - trackRegionMode_ = TrackAssoMode::any; - else - throw cms::Exception( - "Configuration", - "Unsupported value for trackRegionMode: " + trackRegionMode + " (allowed are 'atVertex', 'atCalo', 'any')"); - } else { - // start off with a dummy region - unsigned int ncalomax = 0, nemcalomax = 0, ntrackmax = 0, nmuonmax = 0, npfmax = 0, npuppimax = 0; - regions_.emplace_back(-5.5, - 5.5, - 0, - 2 * M_PI, - 0.5, - 0.5, - useRelativeRegionalCoordinates_, - ncalomax, - nemcalomax, - ntrackmax, - nmuonmax, - npfmax, - npuppimax); - } -} - -void RegionMapper::clear() { - for (Region &r : regions_) - r.zero(); - clusterRefMap_.clear(); - trackRefMap_.clear(); -} - -void RegionMapper::addTrack(const l1t::PFTrack &t) { - // now let's be optimistic and make things very simple - // we propagate in floating point the track to the calo - // we add the track to the region corresponding to its vertex (eta,phi) coordinates AND its (eta,phi) calo coordinates - for (Region &r : regions_) { - bool inside = true; - switch (trackRegionMode_) { - case TrackAssoMode::atVertex: - inside = r.contains(t.eta(), t.phi()); - break; - case TrackAssoMode::atCalo: - inside = r.contains(t.caloEta(), t.caloPhi()); - break; - case TrackAssoMode::any: - inside = r.contains(t.eta(), t.phi()) || r.contains(t.caloEta(), t.caloPhi()); - break; - } - if (inside) { - PropagatedTrack prop; - prop.fillInput(t.pt(), r.localEta(t.eta()), r.localPhi(t.phi()), t.charge(), t.vertex().Z(), t.quality(), &t); - prop.fillPropagated(t.pt(), - t.trkPtError(), - t.caloPtError(), - r.localEta(t.caloEta()), - r.localPhi(t.caloPhi()), - t.quality(), - t.isMuon()); - prop.hwStubs = t.nStubs(); - prop.hwChi2 = round(t.chi2() * 10); - r.track.push_back(prop); - } - } -} -void RegionMapper::addTrack(const l1t::PFTrack &t, l1t::PFTrackRef ref) { - addTrack(t); - trackRefMap_[&t] = ref; -} - -void RegionMapper::addMuon(const l1t::Muon &mu) { - // now let's be optimistic and make things very simple - // we don't propagate anything - for (Region &r : regions_) { - if (r.contains(mu.eta(), mu.phi())) { - Muon prop; - prop.fill(mu.pt(), r.localEta(mu.eta()), r.localPhi(mu.phi()), mu.charge(), mu.hwQual(), &mu); - r.muon.push_back(prop); - } - } -} - -void RegionMapper::addMuon(const l1t::TkMuon &mu) { - // now let's be optimistic and make things very simple - // we don't propagate anything - for (Region &r : regions_) { - if (r.contains(mu.eta(), mu.phi())) { - Muon prop; - prop.fill(mu.pt(), r.localEta(mu.eta()), r.localPhi(mu.phi()), mu.charge(), mu.hwQual()); - r.muon.push_back(prop); - } - } -} - -void RegionMapper::addCalo(const l1t::PFCluster &p) { - if (p.pt() == 0) - return; - for (Region &r : regions_) { - if (r.contains(p.eta(), p.phi())) { - CaloCluster calo; - calo.fill(p.pt(), p.emEt(), p.ptError(), r.localEta(p.eta()), r.localPhi(p.phi()), p.isEM(), 0, &p); - r.calo.push_back(calo); - } - } -} -void RegionMapper::addCalo(const l1t::PFCluster &p, l1t::PFClusterRef ref) { - addCalo(p); - clusterRefMap_[&p] = ref; -} - -void RegionMapper::addEmCalo(const l1t::PFCluster &p) { - if (p.pt() == 0) - return; - for (Region &r : regions_) { - if (r.contains(p.eta(), p.phi())) { - CaloCluster calo; - calo.fill(p.pt(), p.emEt(), p.ptError(), r.localEta(p.eta()), r.localPhi(p.phi()), p.isEM(), p.hwQual(), &p); - r.emcalo.push_back(calo); - } - } -} -void RegionMapper::addEmCalo(const l1t::PFCluster &p, l1t::PFClusterRef ref) { - addEmCalo(p); - clusterRefMap_[&p] = ref; -} - -std::unique_ptr RegionMapper::fetch(bool puppi, float ptMin) const { - auto ret = std::make_unique(); - for (const Region &r : regions_) { - for (const PFParticle &p : (puppi ? r.puppi : r.pf)) { - bool inside = true; - switch (trackRegionMode_) { - case TrackAssoMode::atVertex: - inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi()); - break; - case TrackAssoMode::atCalo: - inside = r.fiducialLocal(p.floatEta(), p.floatPhi()); - break; - case TrackAssoMode::any: - inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi()); - break; // WARNING: this may not be the best choice - } - if (!inside) - continue; - if (p.floatPt() > ptMin) { - reco::Particle::PolarLorentzVector p4( - p.floatPt(), r.globalEta(p.floatVtxEta()), r.globalPhi(p.floatVtxPhi()), 0.13f); - ret->emplace_back( - l1t::PFCandidate::ParticleType(p.hwId), p.intCharge(), p4, p.floatPuppiW(), p.hwPt, p.hwEta, p.hwPhi); - ret->back().setZ0(p.floatDZ()); - ret->back().setHwZ0(p.track.hwZ0); - ret->back().setStatus(p.hwStatus); - if (p.cluster.src) { - auto match = clusterRefMap_.find(p.cluster.src); - if (match == clusterRefMap_.end()) { - throw cms::Exception("CorruptData") << "Invalid cluster pointer in PF candidate id " << p.hwId << " pt " - << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi(); - } - ret->back().setPFCluster(match->second); - } - if (p.track.src) { - auto match = trackRefMap_.find(p.track.src); - if (match == trackRefMap_.end()) { - throw cms::Exception("CorruptData") << "Invalid track pointer in PF candidate id " << p.hwId << " pt " - << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi(); - } - ret->back().setPFTrack(match->second); - } - } - } - } - return ret; -} - -std::unique_ptr RegionMapper::fetchCalo(float ptMin, bool emcalo) const { - auto ret = std::make_unique(); - for (const Region &r : regions_) { - for (const CaloCluster &p : (emcalo ? r.emcalo : r.calo)) { - if (!r.fiducialLocal(p.floatEta(), p.floatPhi())) - continue; - if (p.floatPt() > ptMin) { - reco::Particle::PolarLorentzVector p4(p.floatPt(), r.globalEta(p.floatEta()), r.globalPhi(p.floatPhi()), 0.13f); - l1t::PFCandidate::ParticleType kind = - (p.isEM || emcalo) ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron; - ret->emplace_back(kind, /*charge=*/0, p4, /*puppiW=*/1, p.hwPt, p.hwEta, p.hwPhi); - if (p.src) { - auto match = clusterRefMap_.find(p.src); - if (match == clusterRefMap_.end()) { - throw cms::Exception("CorruptData") - << "Invalid cluster pointer in cluster pt " << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi(); - } - ret->back().setPFCluster(match->second); - } - } - } - } - return ret; -} - -std::unique_ptr RegionMapper::fetchTracks(float ptMin, bool fromPV) const { - auto ret = std::make_unique(); - for (const Region &r : regions_) { - for (const PropagatedTrack &p : r.track) { - if (fromPV && !p.fromPV) - continue; - bool inside = true; - switch (trackRegionMode_) { - case TrackAssoMode::atVertex: - inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi()); - break; - case TrackAssoMode::atCalo: - inside = r.fiducialLocal(p.floatEta(), p.floatPhi()); - break; - case TrackAssoMode::any: - inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi()); - break; // WARNING: this may not be the best choice - } - if (!inside) - continue; - if (p.floatPt() > ptMin) { - reco::Particle::PolarLorentzVector p4( - p.floatVtxPt(), r.globalEta(p.floatVtxEta()), r.globalPhi(p.floatVtxPhi()), 0.13f); - l1t::PFCandidate::ParticleType kind = p.muonLink ? l1t::PFCandidate::Muon : l1t::PFCandidate::ChargedHadron; - ret->emplace_back(kind, p.intCharge(), p4, /*puppiW=*/float(p.fromPV), p.hwPt, p.hwEta, p.hwPhi); - ret->back().setZ0(p.floatDZ()); - ret->back().setHwZ0(p.hwZ0); - if (p.src) { - auto match = trackRefMap_.find(p.src); - if (match == trackRefMap_.end()) { - throw cms::Exception("CorruptData") - << "Invalid track pointer in PF track pt " << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi(); - } - ret->back().setPFTrack(match->second); - } - } - } - } - return ret; -} - -void RegionMapper::putEgObjects(edm::Event &iEvent, - const bool writeEgSta, - const std::string &egLablel, - const std::string &tkEmLabel, - const std::string &tkEleLabel, - const float ptMin) const { - auto egs = std::make_unique>(); - auto tkems = std::make_unique(); - auto tkeles = std::make_unique(); - - edm::RefProd> ref_egs; - if (writeEgSta) - ref_egs = iEvent.getRefBeforePut>(egLablel); - - edm::Ref>::key_type idx = 0; - - for (const Region &r : regions_) { - for (const auto &egphoton : r.egphotons) { - if (egphoton.floatPt() < ptMin) - continue; - - if (!r.fiducialLocal(egphoton.floatEta(), egphoton.floatPhi())) - continue; - - edm::Ref> reg; - auto mom = reco::Candidate::PolarLorentzVector( - egphoton.floatPt(), r.globalEta(egphoton.floatEta()), r.globalPhi(egphoton.floatPhi()), 0.); - if (writeEgSta) { - l1t::EGamma eg(mom); - eg.setHwQual(egphoton.hwQual); - egs->push_back(0, eg); - reg = edm::Ref>(ref_egs, idx++); - } else { - auto egptr = egphoton.cluster.src->constituentsAndFractions()[0].first; - reg = edm::Ref>(egptr.id(), dynamic_cast(egptr.get()), egptr.key()); - } - - l1t::TkEm tkem(reco::Candidate::LorentzVector(mom), reg, egphoton.floatIso(), egphoton.floatIsoPV()); - tkem.setHwQual(egphoton.hwQual); - tkem.setPFIsol(egphoton.floatPFIso()); - tkem.setPFIsolPV(egphoton.floatPFIsoPV()); - tkems->push_back(tkem); - - if (egphoton.ele_idx == -1) - continue; - - const auto &egele = r.egeles[egphoton.ele_idx]; - - if (!r.fiducialLocal(egele.floatEta(), egele.floatPhi())) - continue; - - auto mom_ele = reco::Candidate::PolarLorentzVector( - egele.floatPt(), r.globalEta(egele.floatEta()), r.globalPhi(egele.floatPhi()), 0.); - - l1t::TkElectron tkele( - reco::Candidate::LorentzVector(mom_ele), reg, edm::refToPtr(egele.track.src->track()), egele.floatIso()); - tkele.setHwQual(egele.hwQual); - tkele.setPFIsol(egele.floatPFIso()); - tkeles->push_back(tkele); - } - } - if (writeEgSta) - iEvent.put(std::move(egs), egLablel); - iEvent.put(std::move(tkems), tkEmLabel); - iEvent.put(std::move(tkeles), tkEleLabel); -} - -std::pair RegionMapper::totAndMaxInput(int type) const { - unsigned ntot = 0, nmax = 0; - for (const auto &r : regions_) { - unsigned int ni = r.nInput(Region::InputType(type)); - ntot += ni; - nmax = std::max(nmax, ni); - } - return std::make_pair(ntot, nmax); -} - -std::unique_ptr> RegionMapper::vecInput(int type) const { - auto v = std::make_unique>(); - for (const auto &r : regions_) { - unsigned ni = r.nInput(Region::InputType(type)); - v->push_back(ni); - } - return v; -} - -std::pair RegionMapper::totAndMaxOutput(int type, bool puppi) const { - unsigned ntot = 0, nmax = 0; - for (const auto &r : regions_) { - unsigned int ni = r.nOutput(Region::OutputType(type), puppi); - ntot += ni; - nmax = std::max(nmax, ni); - } - return std::make_pair(ntot, nmax); -} - -std::unique_ptr> RegionMapper::vecOutput(int type, bool puppi) const { - auto v = std::make_unique>(); - for (const auto &r : regions_) { - unsigned ni = r.nOutput(Region::OutputType(type), puppi); - v->push_back(ni); - } - return v; -} diff --git a/L1Trigger/Phase2L1ParticleFlow/src/deregionizer/deregionizer_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/deregionizer/deregionizer_ref.cpp index 6ddcb13f2e97a..7b900c4b0b2ea 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/deregionizer/deregionizer_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/deregionizer/deregionizer_ref.cpp @@ -3,9 +3,10 @@ #include #include +#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" + #ifdef CMSSW_GIT_HASH #include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" l1ct::DeregionizerEmulator::DeregionizerEmulator(const edm::ParameterSet &iConfig) : DeregionizerEmulator(iConfig.getParameter("nPuppiFinalBuffer"), @@ -15,8 +16,6 @@ l1ct::DeregionizerEmulator::DeregionizerEmulator(const edm::ParameterSet &iConfi iConfig.getParameter("nPuppiThirdBuffers")) { debug_ = iConfig.getUntrackedParameter("debug", false); } -#else -#include "../../utils/dbgPrintf.h" #endif l1ct::DeregionizerEmulator::DeregionizerEmulator(const unsigned int nPuppiFinalBuffer /*=128*/, diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index a100230f080a9..2fc881879455d 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -6,6 +6,7 @@ #include #include #include +#include using namespace l1ct; @@ -21,6 +22,7 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) doBremRecovery(pset.getParameter("doBremRecovery")), writeBeforeBremRecovery(pset.getParameter("writeBeforeBremRecovery")), caloHwQual(pset.getParameter("caloHwQual")), + doEndcapHwQual(pset.getParameter("doEndcapHwQual")), emClusterPtMin(pset.getParameter("caloEtMin")), dEtaMaxBrem(pset.getParameter("dEtaMaxBrem")), dPhiMaxBrem(pset.getParameter("dPhiMaxBrem")), @@ -165,7 +167,8 @@ void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const { if (calo.hwPt > 0) dbgCout() << "[REF] IN calo[" << ic << "] pt: " << calo.hwPt << " eta: " << calo.hwEta << " (glb eta: " << in.region.floatGlbEta(calo.hwEta) << ") phi: " << calo.hwPhi - << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << calo.hwEmID << std::endl; + << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << std::bitset<4>(calo.hwEmID) + << std::endl; } } @@ -221,7 +224,15 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion, // check if brem recovery is on if (!cfg.doBremRecovery || cfg.writeBeforeBremRecovery) { // 1. create EG objects before brem recovery - addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, calo.hwPt, itk); + unsigned int egQual = calo.hwEmID; + // If we write both objects with and without brem-recovery + // bit 3 is used for the brem-recovery bit: if set = no recovery + // (for consistency with the barrel hwQual where by default the brem recovery is done upstream) + if (cfg.writeBeforeBremRecovery && cfg.doBremRecovery) { + egQual = calo.hwEmID | 0x8; + } + + addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk); } if (!cfg.doBremRecovery) @@ -244,13 +255,13 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion, // 2. create EG objects with brem recovery // NOTE: duplicating the object is suboptimal but this is done for keeping things as in TDR code... - addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID + 2, ptBremReco, itk, components); + addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, components); } } EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector &egobjs, const EmCaloObjEmu &calo, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr, const std::vector &components) const { EGObjEmu egsta; @@ -260,25 +271,36 @@ EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector &egobjs, egsta.hwPhi = calo.hwPhi; egsta.hwQual = hwQual; egobjs.push_back(egsta); + + if (debug_ > 2) + dbgCout() << "[REF] EGSta pt: " << egsta.hwPt << " eta: " << egsta.hwEta << " phi: " << egsta.hwPhi + << " qual: " << std::bitset<4>(egsta.hwQual) << " packed: " << egsta.pack().to_string(16) << std::endl; + return egobjs.back(); } EGIsoObjEmu &PFTkEGAlgoEmulator::addEGIsoToPF(std::vector &egobjs, const EmCaloObjEmu &calo, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr) const { EGIsoObjEmu egiso; egiso.clear(); egiso.hwPt = ptCorr; egiso.hwEta = calo.hwEta; egiso.hwPhi = calo.hwPhi; - egiso.hwQual = hwQual; + unsigned int egHwQual = hwQual; + if (cfg.doEndcapHwQual) { + // 1. zero-suppress the loose EG-ID (bit 1) + // 2. for now use the standalone tight definition (bit 0) to set the tight point for photons (bit 2) + egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 2); + } + egiso.hwQual = egHwQual; egiso.srcCluster = calo.src; egobjs.push_back(egiso); if (debug_ > 2) dbgCout() << "[REF] EGIsoObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi - << " qual: " << egiso.hwQual << " packed: " << egiso.pack().to_string(16) << std::endl; + << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl; return egobjs.back(); } @@ -286,14 +308,20 @@ EGIsoObjEmu &PFTkEGAlgoEmulator::addEGIsoToPF(std::vector &egobjs, EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector &egobjs, const EmCaloObjEmu &calo, const TkObjEmu &track, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr) const { EGIsoEleObjEmu egiso; egiso.clear(); egiso.hwPt = ptCorr; egiso.hwEta = calo.hwEta; egiso.hwPhi = calo.hwPhi; - egiso.hwQual = hwQual; + unsigned int egHwQual = hwQual; + if (cfg.doEndcapHwQual) { + // 1. zero-suppress the loose EG-ID (bit 1) + // 2. for now use the standalone tight definition (bit 0) to set the tight point for eles (bit 1) + egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1); + } + egiso.hwQual = egHwQual; egiso.hwDEta = track.hwVtxEta() - egiso.hwEta; egiso.hwDPhi = abs(track.hwVtxPhi() - egiso.hwPhi); egiso.hwZ0 = track.hwZ0; @@ -304,7 +332,7 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector if (debug_ > 2) dbgCout() << "[REF] EGIsoEleObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi - << " qual: " << egiso.hwQual << " packed: " << egiso.pack().to_string(16) << std::endl; + << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl; return egobjs.back(); } @@ -315,7 +343,7 @@ void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector &egstas, const std::vector &emcalo, const std::vector &track, const int calo_idx, - const int hwQual, + const unsigned int hwQual, const pt_t ptCorr, const int tk_idx, const std::vector &components) const { diff --git a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/muonGmtToL1ct_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/muonGmtToL1ct_ref.cpp index b6ad2863ffa46..0c3ddda6a1c43 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/muonGmtToL1ct_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/muonGmtToL1ct_ref.cpp @@ -23,20 +23,20 @@ l1ct::MuObjEmu l1ct::GMTMuonDecoderEmulator::decode(const ap_uint<64> &in) const const int z0_scale = std::round(z0Scale_ / l1ct::Scales::Z0_LSB); const int dxy_scale = std::round(dxyScale_ / l1ct::Scales::DXY_LSB); - bool gmt_chg = in[55]; - ap_uint<13> gmt_ipt = in(15, 0); - ap_int<13> gmt_phi = in(28, 16); - ap_int<14> gmt_eta = in(42, 29); - ap_int<5> gmt_z0 = in(47, 43); - ap_int<7> gmt_d0 = in(54, 48); - ap_uint<4> gmt_qual = in(59, 56); + bool gmt_valid = in[0], gmt_chg = in[56]; + ap_uint<13> gmt_ipt = in(16, 1); + ap_int<13> gmt_phi = in(29, 17); + ap_int<14> gmt_eta = in(43, 30); + ap_int<5> gmt_z0 = in(48, 44); + ap_int<7> gmt_d0 = in(55, 49); + ap_uint<4> gmt_qual = in(60, 57); gmt_pt_t gmt_pt; gmt_pt(gmt_pt_t::width - 1, 0) = gmt_ipt(gmt_pt_t::width - 1, 0); // copy the bits l1ct::MuObjEmu out; out.clear(); - if (gmt_pt != 0) { + if (gmt_valid && gmt_pt != 0) { // add a shift in order to get the proper rounding out.hwPt = gmt_pt + gmt_pt_t(l1ct::Scales::INTPT_LSB / 2); diff --git a/L1Trigger/Phase2L1ParticleFlow/src/pf/pfalgo_common_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/pf/pfalgo_common_ref.cpp index 869c5736e946a..1e7fa983f559a 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/pf/pfalgo_common_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/pf/pfalgo_common_ref.cpp @@ -81,7 +81,7 @@ void l1ct::PFAlgoEmulatorBase::pfalgo_mu_ref(const PFInputRegion &in, OutputRegi for (unsigned int im = 0; im < nMU; ++im) { if (in.muon[im].hwPt > 0) { int ibest = -1; - pt_t dptmin = in.muon[im].hwPt >> 1; + pt_t dptmin = (in.muon[im].hwPt << 1) + in.muon[im].hwPt; for (unsigned int it = 0; it < nTRACK; ++it) { if (!in.track[it].isPFLoose()) continue; diff --git a/L1Trigger/Phase2L1ParticleFlow/src/puppi/linpuppi_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/puppi/linpuppi_ref.cpp index 950acb338adb6..8ef2679747d0d 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/puppi/linpuppi_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/puppi/linpuppi_ref.cpp @@ -5,12 +5,7 @@ #include "L1Trigger/Phase2L1ParticleFlow/interface/common/bitonic_hybrid_sort_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/common/bitonic_sort_ref.h" - -#ifdef CMSSW_GIT_HASH #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h" -#else -#include "../utils/dbgPrintf.h" -#endif #ifdef CMSSW_GIT_HASH #include "FWCore/ParameterSet/interface/ParameterSet.h" diff --git a/L1Trigger/Phase2L1ParticleFlow/src/regionizer/multififo_regionizer_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/regionizer/multififo_regionizer_ref.cpp index c91abe9d7db8f..feb87ce4b9dfe 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/regionizer/multififo_regionizer_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/regionizer/multififo_regionizer_ref.cpp @@ -39,7 +39,7 @@ l1ct::MultififoRegionizerEmulator::MultififoRegionizerEmulator(unsigned int nend NTK_SECTORS(9), NCALO_SECTORS(3), NTK_LINKS(2), - NCALO_LINKS(2), + NCALO_LINKS(3), HCAL_LINKS(0), ECAL_LINKS(0), NMU_LINKS(1), @@ -72,15 +72,13 @@ l1ct::MultififoRegionizerEmulator::MultififoRegionizerEmulator(unsigned int nend // hgcal assert(NCALO_SECTORS == 3 && NTK_SECTORS == 9); // otherwise math below is broken, but it's hard to make it generic for (unsigned int ie = 0; ie < nendcaps; ++ie) { - for (unsigned int is = 0; is < NCALO_SECTORS; ++is) { // NCALO_SECTORS sectors - for (unsigned int il = 0; il < NCALO_LINKS; ++il) { // max clusters per sector per clock - for (unsigned int j = 0; j < 3; ++j) { // PF REGION - caloRoutes_.emplace_back(is + 3 * ie, il, 3 * is + j + 9 * ie, il); - if (j == 0 || j == 2) { - int other = (j == 0) ? 2 : 1; // pf region 0, takes from prev. pf region 2 takes from next - // from sector , from link, to region, to fifo - caloRoutes_.emplace_back( - (is + other) % 3 + 3 * ie, il, 3 * is + j + 9 * ie, il + 2); //last 2 = NCALOFIBERS + for (unsigned int is = 0; is < NCALO_SECTORS; ++is) { // NCALO_SECTORS sectors + for (unsigned int il = 0; il < NCALO_LINKS; ++il) { // max clusters per sector per clock + for (unsigned int j = 0; j < 3; ++j) { // PF REGION + caloRoutes_.emplace_back(is + 3 * ie, il, 3 * is + j + 9 * ie, il); // 4 args are: sector, link, region, fifo + if (j != 2) { // pf regions 0 and 1 take also from previous sector + int isprev = (is > 0 ? is - 1 : NCALO_SECTORS - 1); + caloRoutes_.emplace_back(isprev + 3 * ie, il, 3 * is + j + 9 * ie, il + NCALO_LINKS); } } } diff --git a/L1Trigger/Phase2L1ParticleFlow/test/.gitignore b/L1Trigger/Phase2L1ParticleFlow/test/.gitignore new file mode 100644 index 0000000000000..26f0a8dafc043 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/test/.gitignore @@ -0,0 +1,2 @@ +*.txt +*.dump diff --git a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py index 924c8d01a6bd0..cf10f843d9629 100644 --- a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py +++ b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py @@ -1,8 +1,6 @@ import FWCore.ParameterSet.Config as cms from Configuration.StandardSequences.Eras import eras -produceEGStage2Pattern = False - process = cms.Process("RESP", eras.Phase2C9) process.load('Configuration.StandardSequences.Services_cff') @@ -17,7 +15,8 @@ inputCommands = cms.untracked.vstring("keep *", "drop l1tPFClusters_*_*_*", "drop l1tPFTracks_*_*_*", - "drop l1tPFCandidates_*_*_*") + "drop l1tPFCandidates_*_*_*", + "drop l1tTkPrimaryVertexs_*_*_*") ) process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') @@ -29,7 +28,6 @@ from Configuration.AlCa.GlobalTag import GlobalTag process.GlobalTag = GlobalTag(process.GlobalTag, '123X_mcRun4_realistic_v3', '') -process.load("L1Trigger.Phase2L1ParticleFlow.l1ParticleFlow_cff") process.load('L1Trigger.Phase2L1ParticleFlow.l1ctLayer1_cff') process.load('L1Trigger.Phase2L1ParticleFlow.l1ctLayer2EG_cff') process.load('L1Trigger.L1TTrackMatch.l1tGTTInputProducer_cfi') @@ -41,39 +39,33 @@ process.l1tSAMuonsGmt = l1tStandaloneMuons.clone() process.l1tLayer1Barrel9 = process.l1tLayer1Barrel.clone() +process.l1tLayer1Barrel9.puAlgo.nFinalSort = 32 process.l1tLayer1Barrel9.regions[0].etaBoundaries = [ -1.5, -0.5, 0.5, 1.5 ] process.l1tLayer1Barrel9.boards=cms.VPSet( cms.PSet( - regions=cms.vuint32(list(range(0, 3)) + [x+9 for x in range(0, 3)] + [x+18 for x in range(0, 3)])), + regions=cms.vuint32(*[0+9*ie+i for ie in range(3) for i in range(3)])), cms.PSet( - regions=cms.vuint32(list(range(3, 6)) + [x+9 for x in range(3, 6)] + [x+18 for x in range(3, 6)])), + regions=cms.vuint32(*[3+9*ie+i for ie in range(3) for i in range(3)])), cms.PSet( - regions=cms.vuint32(list(range(6, 9)) + [x+9 for x in range(6, 9)] + [x+18 for x in range(6, 9)])), + regions=cms.vuint32(*[6+9*ie+i for ie in range(3) for i in range(3)])), ) process.runPF = cms.Path( process.l1tSAMuonsGmt + process.l1tGTTInputProducer + process.l1tVertexFinderEmulator + - process.l1tPFTracksFromL1Tracks + - process.l1tParticleFlow_calo + process.l1tLayer1Barrel + process.l1tLayer1Barrel9 + process.l1tLayer1HGCal + process.l1tLayer1HGCalNoTK + - process.l1tLayer1HF + - process.l1tLayer1 + - process.l1tLayer2EG - ) + process.l1tLayer1HF +) +process.runPF.associate(process.l1tLayer1TaskInputsTask) + -if produceEGStage2Pattern: - process.l1tLayer2EG.writeInPattern = True - process.l1tLayer2EG.writeOutPattern = True +for det in "Barrel", "Barrel9", "HGCal", "HGCalNoTK", "HF": + l1pf = getattr(process, 'l1ctLayer1'+det) + l1pf.dumpFileName = cms.untracked.string("TTbar_PU200_"+det+".dump") process.source.fileNames = [ '/store/cmst3/group/l1tr/gpetrucc/11_1_0/NewInputs110X/110121.done/TTbar_PU200/inputs110X_%d.root' % i for i in (1,3,7,8,9) ] process.l1tPFClustersFromCombinedCaloHCal.phase2barrelCaloTowers = [cms.InputTag("l1tEGammaClusterEmuProducer",)] - - -for det in "Barrel", "Barrel9", "HGCal", "HGCalNoTK", "HF": - l1pf = getattr(process, 'l1tLayer1'+det) - l1pf.dumpFileName = cms.untracked.string("TTbar_PU200_110X_"+det+".dump") diff --git a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_patternFiles_fromRAW_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_fromRAW_cfg.py similarity index 100% rename from L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_patternFiles_fromRAW_cfg.py rename to L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_fromRAW_cfg.py diff --git a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_patternFiles_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_patternFiles_cfg.py new file mode 100644 index 0000000000000..e55a90022aa50 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_patternFiles_cfg.py @@ -0,0 +1,98 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras + +process = cms.Process("RESP", eras.Phase2C9) + +process.load('Configuration.StandardSequences.Services_cff') +process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") +process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True), allowUnscheduled = cms.untracked.bool(False) ) +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000)) +process.MessageLogger.cerr.FwkReport.reportEvery = 1 + +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('file:inputs110X.root'), + inputCommands = cms.untracked.vstring("keep *", + "drop l1tPFClusters_*_*_*", + "drop l1tPFTracks_*_*_*", + "drop l1tPFCandidates_*_*_*", + "drop l1tTkPrimaryVertexs_*_*_*") +) + +process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D49_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('SimCalorimetry.HcalTrigPrimProducers.hcaltpdigi_cff') # needed to read HCal TPs +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, '123X_mcRun4_realistic_v3', '') + +process.load('L1Trigger.Phase2L1ParticleFlow.l1ctLayer1_cff') +process.load('L1Trigger.Phase2L1ParticleFlow.l1ctLayer2EG_cff') +process.load('L1Trigger.L1TTrackMatch.l1tGTTInputProducer_cfi') +process.load('L1Trigger.VertexFinder.l1tVertexProducer_cfi') +process.l1tVertexFinderEmulator = process.VertexProducer.clone() +process.l1tVertexFinderEmulator.VertexReconstruction.Algorithm = "fastHistoEmulation" +process.l1tVertexFinderEmulator.l1TracksInputTag = cms.InputTag("l1tGTTInputProducer", "Level1TTTracksConverted") +from L1Trigger.Phase2L1GMT.gmt_cfi import l1tStandaloneMuons +process.l1tSAMuonsGmt = l1tStandaloneMuons.clone() + +from L1Trigger.Phase2L1ParticleFlow.l1tSeedConePFJetProducer_cfi import l1tSeedConePFJetEmulatorProducer +from L1Trigger.Phase2L1ParticleFlow.l1tDeregionizerProducer_cfi import l1tDeregionizerProducer +from L1Trigger.Phase2L1ParticleFlow.l1tJetFileWriter_cfi import l1tSeededConeJetFileWriter +process.l1tLayer2Deregionizer = l1tDeregionizerProducer.clone() +process.l1tLayer2SeedConeJets = l1tSeedConePFJetEmulatorProducer.clone(L1PFObject = cms.InputTag('l1tLayer2Deregionizer', 'Puppi')) +process.l1tLayer2SeedConeJetWriter = l1tSeededConeJetFileWriter.clone(jets = "l1tLayer2SeedConeJets") + +process.l1tLayer1Barrel9 = process.l1tLayer1Barrel.clone() +process.l1tLayer1Barrel9.puAlgo.nFinalSort = 32 +process.l1tLayer1Barrel9.regions[0].etaBoundaries = [ -1.5, -0.5, 0.5, 1.5 ] +process.l1tLayer1Barrel9.boards=cms.VPSet( + cms.PSet( + regions=cms.vuint32(*[0+9*ie+i for ie in range(3) for i in range(3)])), + cms.PSet( + regions=cms.vuint32(*[3+9*ie+i for ie in range(3) for i in range(3)])), + cms.PSet( + regions=cms.vuint32(*[6+9*ie+i for ie in range(3) for i in range(3)])), + ) + +from L1Trigger.Phase2L1ParticleFlow.l1ctLayer1_patternWriters_cff import * +process.l1tLayer1Barrel.patternWriters = cms.untracked.VPSet(*barrelWriterConfigs) +#process.l1tLayer1Barrel9.patternWriters = cms.untracked.VPSet(*barrel9WriterConfigs) # not enabled for now +process.l1tLayer1HGCal.patternWriters = cms.untracked.VPSet(*hgcalWriterConfigs) +process.l1tLayer1HGCalNoTK.patternWriters = cms.untracked.VPSet(*hgcalNoTKWriterConfigs) +process.l1tLayer1HF.patternWriters = cms.untracked.VPSet(*hfWriterConfigs) + +process.runPF = cms.Path( + process.l1tSAMuonsGmt + + process.l1tGTTInputProducer + + process.l1tVertexFinderEmulator + + process.l1tLayer1Barrel + + #process.l1tLayer1Barrel9 + + process.l1tLayer1HGCal + + process.l1tLayer1HGCalNoTK + + process.l1tLayer1HF + + process.l1tLayer1 + + process.l1tLayer2EG + + process.l1tLayer2Deregionizer + + process.l1tLayer2SeedConeJets + + process.l1tLayer2SeedConeJetWriter + ) +process.runPF.associate(process.l1tLayer1TaskInputsTask) + + +##################################################################################################################### +## Layer 2 e/gamma + +process.l1tLayer2EG.writeInPattern = True +process.l1tLayer2EG.writeOutPattern = True +process.l1tLayer2EG.inPatternFile.maxLinesPerFile = eventsPerFile_*54 +process.l1tLayer2EG.outPatternFile.maxLinesPerFile = eventsPerFile_*54 + +##################################################################################################################### +## Layer 2 seeded-cone jets +process.l1tLayer2SeedConeJetWriter.maxLinesPerFile = eventsPerFile_*54 + +process.source.fileNames = [ '/store/cmst3/group/l1tr/gpetrucc/11_1_0/NewInputs110X/110121.done/TTbar_PU200/inputs110X_%d.root' % i for i in (1,3,7,8,9) ] +process.pfClustersFromCombinedCaloHCal.phase2barrelCaloTowers = [cms.InputTag("l1tEGammaClusterEmuProducer",)] diff --git a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_patternFiles_fromRAW_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_patternFiles_fromRAW_cfg.py new file mode 100644 index 0000000000000..50e3db4e23812 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_patternFiles_fromRAW_cfg.py @@ -0,0 +1,100 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras + +process = cms.Process("RESP", eras.Phase2C9) + +process.load('Configuration.StandardSequences.Services_cff') +process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") +process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.options = cms.untracked.PSet( + wantSummary = cms.untracked.bool(True), + numberOfThreads = cms.untracked.uint32(2), + numberOfStreams = cms.untracked.uint32(1), +) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(100)) +process.MessageLogger.cerr.FwkReport.reportEvery = 1 + +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/mc/Phase2HLTTDRWinter20DIGI/TT_TuneCP5_14TeV-powheg-pythia8/GEN-SIM-DIGI-RAW/PU200_110X_mcRun4_realistic_v3-v2/110000/005E74D6-B50E-674E-89E6-EAA9A617B476.root',) +) + +process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D49_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, '123X_mcRun4_realistic_v3', '') + +process.load('SimCalorimetry.HcalTrigPrimProducers.hcaltpdigi_cff') # needed to read HCal TPs +process.load('CalibCalorimetry.CaloTPG.CaloTPGTranscoder_cfi') +process.load('Configuration.StandardSequences.SimL1Emulator_cff') +process.load('L1Trigger.TrackTrigger.TrackTrigger_cff') +process.load("L1Trigger.TrackFindingTracklet.L1HybridEmulationTracks_cff") +process.load("L1Trigger.TrackerDTC.ProducerES_cff") +process.load("L1Trigger.TrackerDTC.ProducerED_cff") +process.load("RecoVertex.BeamSpotProducer.BeamSpot_cfi") + +from L1Trigger.Phase2L1ParticleFlow.l1tSeedConePFJetProducer_cfi import l1tSeedConePFJetEmulatorProducer +from L1Trigger.Phase2L1ParticleFlow.l1tDeregionizerProducer_cfi import l1tDeregionizerProducer +from L1Trigger.Phase2L1ParticleFlow.l1tJetFileWriter_cfi import l1tSeededConeJetFileWriter +process.l1tLayer2Deregionizer = l1tDeregionizerProducer.clone() +process.l1tLayer2SeedConeJets = l1tSeedConePFJetEmulatorProducer.clone(L1PFObject = cms.InputTag('l1tLayer2Deregionizer', 'Puppi')) +process.l1tLayer2SeedConeJetWriter = l1tSeededConeJetFileWriter.clone(jets = "l1tLayer2SeedConeJets") + +process.l1tLayer1Barrel9 = process.l1tLayer1Barrel.clone() +process.l1tLayer1Barrel9.puAlgo.nFinalSort = 32 +process.l1tLayer1Barrel9.regions[0].etaBoundaries = [ -1.5, -0.5, 0.5, 1.5 ] +process.l1tLayer1Barrel9.boards=cms.VPSet( + cms.PSet( + regions=cms.vuint32(*[0+9*ie+i for ie in range(3) for i in range(3)])), + cms.PSet( + regions=cms.vuint32(*[3+9*ie+i for ie in range(3) for i in range(3)])), + cms.PSet( + regions=cms.vuint32(*[6+9*ie+i for ie in range(3) for i in range(3)])), + ) + +from L1Trigger.Phase2L1ParticleFlow.l1tLayer1_patternWriters_cff import * +process.l1tLayer1Barrel.patternWriters = cms.untracked.VPSet(*barrelWriterConfigs) +#process.l1tLayer1Barrel9.patternWriters = cms.untracked.VPSet(*barrel9WriterConfigs) # not enabled for now +process.l1tLayer1HGCal.patternWriters = cms.untracked.VPSet(*hgcalWriterConfigs) +process.l1tLayer1HGCalNoTK.patternWriters = cms.untracked.VPSet(*hgcalNoTKWriterConfigs) +process.l1tLayer1HF.patternWriters = cms.untracked.VPSet(*hfWriterConfigs) + +process.L1TPFInputsTask = cms.Task( + process.TTClustersFromPhase2TrackerDigis, + process.TTStubsFromPhase2TrackerDigis, + process.TrackerDTCProducer, + process.offlineBeamSpot, + process.l1tTTTracksFromTrackletEmulation, + process.SimL1EmulatorTask +) +process.runPF = cms.Path( + process.l1tLayer1Barrel + + #process.l1tLayer1Barrel9 + + process.l1tLayer1HGCal + + process.l1tLayer1HGCalNoTK + + process.l1tLayer1HF + + process.l1tLayer1 + + process.l1tLayer2EG + + process.l1tLayer2Deregionizer + + process.l1tLayer2SeedConeJets + + process.l1tLayer2SeedConeJetWriter + ) +process.runPF.associate(process.L1TPFInputsTask) +process.schedule = cms.Schedule(process.runPF) + + +##################################################################################################################### +## Layer 2 e/gamma + +process.l1tLayer2EG.writeInPattern = True +process.l1tLayer2EG.writeOutPattern = True +process.l1tLayer2EG.inPatternFile.maxLinesPerFile = eventsPerFile_*54 +process.l1tLayer2EG.outPatternFile.maxLinesPerFile = eventsPerFile_*54 + +##################################################################################################################### +## Layer 2 seeded-cone jets +process.l1tLayer2SeedConeJetWriter.maxLinesPerFile = eventsPerFile_*54