diff --git a/NtupleProducer/plugins/L1PFDecodedCaloTableProducer.cc b/NtupleProducer/plugins/L1PFDecodedCaloTableProducer.cc index 807ab49b..a9d7e850 100644 --- a/NtupleProducer/plugins/L1PFDecodedCaloTableProducer.cc +++ b/NtupleProducer/plugins/L1PFDecodedCaloTableProducer.cc @@ -17,7 +17,7 @@ #include "CommonTools/Utils/interface/StringCutObjectSelector.h" #include "CommonTools/Utils/interface/StringObjectFunction.h" - +#include "L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h" #include @@ -67,14 +67,19 @@ L1PFDecodedCaloTableProducer::produce(edm::StreamID id, edm::Event& iEvent, cons std::vector vals_empt, vals_srrTot, vals_hwSrrTot, vals_meanz, vals_hwMeanZ, vals_hoe, vals_piIdProb, vals_PuIdProb, vals_EmIdProb, - vals_caloIso, vals_showerShape; - std::vector vals_showerlength, + vals_caloIso, vals_caloShowerShape, vals_hwEmID, vals_showerShape, + vals_hwShowerShape, vals_hwRelIso, vals_relIso; + + std::vector vals_showerlength, vals_coreshowerlength, vals_emf, vals_hw_emf, vals_abseta, vals_hw_abseta, vals_hw_meanz, vals_sigmaetaeta, vals_hw_sigmaetaeta, vals_sigmaphiphi, vals_hw_sigmaphiphi, vals_sigmazz, vals_hw_sigmazz; + std::vector vals_caloPt, vals_caloEta, vals_caloPhi; + std::vector vals_relIsoHack, vals_relIsoOld; + vals_empt.resize(ncands); vals_srrTot.resize(ncands); vals_hwSrrTot.resize(ncands); @@ -85,7 +90,12 @@ L1PFDecodedCaloTableProducer::produce(edm::StreamID id, edm::Event& iEvent, cons vals_PuIdProb.resize(ncands); vals_EmIdProb.resize(ncands); vals_caloIso.resize(ncands); + vals_caloShowerShape.resize(ncands); + vals_hwEmID.resize(ncands); vals_showerShape.resize(ncands); + vals_hwShowerShape.resize(ncands); + vals_hwRelIso.resize(ncands); + vals_relIso.resize(ncands); vals_showerlength.resize(ncands); vals_coreshowerlength.resize(ncands); vals_emf.resize(ncands); @@ -99,7 +109,11 @@ L1PFDecodedCaloTableProducer::produce(edm::StreamID id, edm::Event& iEvent, cons vals_hw_sigmaphiphi.resize(ncands); vals_sigmazz.resize(ncands); vals_hw_sigmazz.resize(ncands); - + vals_caloPt.resize(ncands); + vals_caloEta.resize(ncands); + vals_caloPhi.resize(ncands); + vals_relIsoHack.resize(ncands); + vals_relIsoOld.resize(ncands); for (unsigned int i = 0; i < ncands; ++i) { const auto cand = selected[i]; @@ -115,10 +129,42 @@ L1PFDecodedCaloTableProducer::produce(edm::StreamID id, edm::Event& iEvent, cons vals_PuIdProb[i] = digi->floatPuProb(); vals_EmIdProb[i] = digi->floatEmProb(); + vals_hwEmID[i] = digi->hwEmID; + vals_showerShape[i] = digi->floatShowerShape(); + vals_hwShowerShape[i] = digi->hwShowerShape; + vals_hwRelIso[i] = digi->hwRelIso; + vals_relIso[i] = digi->floatRelIso(); + + const l1tp2::CaloCrystalCluster *crycl = dynamic_cast(cand->constituentsAndFractions().front().first.get()); if(crycl) { vals_caloIso[i] = crycl->isolation(); - vals_showerShape[i] = crycl->e2x5() / crycl->e5x5(); + vals_caloShowerShape[i] = crycl->e2x5() / crycl->e5x5(); + vals_caloPt[i] = crycl->pt(); + vals_caloEta[i] = crycl->eta(); + vals_caloPhi[i] = crycl->phi(); + + ap_ufixed<16, 0> calo_invPt = l1ct::invert_with_shift, 1024>(digi->hwPt); + + float reliso_hack = digi->hwPt == 1 ? digi->floatRelIso() : digi->floatRelIso() * calo_invPt.to_float(); + l1ct::rel_iso_t hw_reliso_hack = l1ct::Scales::makeRelIso(reliso_hack); + + vals_relIsoHack[i] = l1ct::Scales::floatRelIso(hw_reliso_hack); + + + vals_relIsoOld[i] = l1ct::Scales::floatRelIso(l1ct::Scales::makeRelIso(crycl->isolation() / digi->floatPt())); + // std::cout << " reliso (float): " << crycl->isolation() + // << " reliso (emu): " << digi->floatRelIso() + // << " reliso buggy (float): " << crycl->isolation() / digi->floatPt() + // << " reliso buggy (emu): " << vals_relIsoOld[i] + // << " pt (float): " << digi->floatPt() + // << " pt (emu): " << digi->hwPt.to_float() + // << " invPt (float): " << 1./digi->floatPt() + // << " invPt (emu): " << calo_invPt.to_float() + // << " reliso hack (float): " << vals_relIsoHack[i] + // << std::endl; + + } } else if(auto digi = std::get_if(&obj)){ vals_empt[i] = digi->floatEmPt(); @@ -162,6 +208,13 @@ L1PFDecodedCaloTableProducer::produce(edm::StreamID id, edm::Event& iEvent, cons vals_hw_sigmaphiphi[i] = w_sigmaphiphi.to_float(); vals_sigmazz[i] = w_sigmazz * SIGMAZZ_LSB; vals_hw_sigmazz[i] = w_sigmazz.to_float(); + vals_caloPt[i] = hgcalcl->pt(); + vals_caloEta[i] = hgcalcl->eta(); + vals_caloPhi[i] = hgcalcl->phi(); + } else { + vals_caloPt[i] =cand->constituentsAndFractions().front().first.get()->pt(); + vals_caloEta[i] =cand->constituentsAndFractions().front().first.get()->eta(); + vals_caloPhi[i] =cand->constituentsAndFractions().front().first.get()->phi(); } } } @@ -176,23 +229,37 @@ L1PFDecodedCaloTableProducer::produce(edm::StreamID id, edm::Event& iEvent, cons out->addColumn("piIdProb", vals_piIdProb, ""); out->addColumn("PuIdProb", vals_PuIdProb, ""); out->addColumn("EmIdProb", vals_EmIdProb, ""); - out->addColumn("caloIso", vals_caloIso, ""); - out->addColumn("showerShape", vals_showerShape, ""); + out->addColumn("clRelIso", vals_caloIso, ""); + out->addColumn("clRelIsoHack", vals_relIsoHack, ""); + out->addColumn("clRelIsoOld", vals_relIsoOld, ""); - out->addColumn("showerlength", vals_showerlength, ""); - out->addColumn("coreshowerlength", vals_coreshowerlength, ""); - out->addColumn("emf", vals_emf, ""); - out->addColumn("hwEmf", vals_hw_emf, ""); - out->addColumn("abseta", vals_abseta, ""); - out->addColumn("hwAbseta", vals_hw_abseta, ""); - out->addColumn("hwFPMeanz", vals_hw_meanz, ""); - out->addColumn("sigmaetaeta", vals_sigmaetaeta, ""); - out->addColumn("hwSigmaetaeta", vals_hw_sigmaetaeta, ""); - out->addColumn("sigmaphiphi", vals_sigmaphiphi, ""); - out->addColumn("hwSigmaphiphi", vals_hw_sigmaphiphi, ""); - out->addColumn("sigmazz", vals_sigmazz, ""); - out->addColumn("hwSigmazz", vals_hw_sigmazz, ""); + + out->addColumn("clShowerShape", vals_caloShowerShape, ""); + out->addColumn("hwEmID", vals_hwEmID, ""); + out->addColumn("showerShape", vals_showerShape, ""); + out->addColumn("hwShowerShape", vals_hwShowerShape, ""); + out->addColumn("hwRelIso", vals_hwRelIso, ""); + out->addColumn("relIso", vals_relIso, ""); + + + + out->addColumn("clShowerlength", vals_showerlength, ""); + out->addColumn("clCoreshowerlength", vals_coreshowerlength, ""); + out->addColumn("clEmf", vals_emf, ""); + out->addColumn("clHwEmf", vals_hw_emf, ""); + out->addColumn("clAbseta", vals_abseta, ""); + out->addColumn("clHwAbseta", vals_hw_abseta, ""); + out->addColumn("clHwMeanz", vals_hw_meanz, ""); + out->addColumn("clSigmaetaeta", vals_sigmaetaeta, ""); + out->addColumn("clHwSigmaetaeta", vals_hw_sigmaetaeta, ""); + out->addColumn("clSigmaphiphi", vals_sigmaphiphi, ""); + out->addColumn("clHwSigmaphiphi", vals_hw_sigmaphiphi, ""); + out->addColumn("clSigmazz", vals_sigmazz, ""); + out->addColumn("clHwSigmazz", vals_hw_sigmazz, ""); + out->addColumn("clPt", vals_caloPt, ""); + out->addColumn("clEta", vals_caloEta, ""); + out->addColumn("clPhi", vals_caloPhi, ""); // save to the event branches iEvent.put(std::move(out)); diff --git a/NtupleProducer/plugins/L1PFTkEleTableProducer.cc b/NtupleProducer/plugins/L1PFTkEleTableProducer.cc new file mode 100644 index 00000000..8de9842c --- /dev/null +++ b/NtupleProducer/plugins/L1PFTkEleTableProducer.cc @@ -0,0 +1,107 @@ +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/Common/interface/Handle.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "DataFormats/NanoAOD/interface/FlatTable.h" + +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" + +#include "DataFormats/L1TCorrelator/interface/TkElectron.h" +#include "DataFormats/L1TCorrelator/interface/TkElectronFwd.h" +#include "DataFormats/L1TParticleFlow/interface/egamma.h" + +#include + +class L1PFTkEleTableProducer : public edm::global::EDProducer<> { +public: + explicit L1PFTkEleTableProducer(const edm::ParameterSet&); + ~L1PFTkEleTableProducer() override = default; + +private: + void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; + + std::string name_; + edm::EDGetTokenT> src_; + StringCutObjectSelector sel_; +}; + +L1PFTkEleTableProducer::L1PFTkEleTableProducer(const edm::ParameterSet& iConfig) + : name_(iConfig.getParameter("name")), + src_(consumes>(iConfig.getParameter("src"))), + sel_(iConfig.getParameter("cut"), true) { + produces(); +} + +void L1PFTkEleTableProducer::produce(edm::StreamID id, + edm::Event& iEvent, + const edm::EventSetup& iSetup) const { + edm::Handle> src; + iEvent.getByToken(src_, src); + + std::vector selected; + for (const auto& ele : *src) + if (sel_(ele)) + selected.push_back(&ele); + + unsigned int ncands = selected.size(); + auto out = std::make_unique(ncands, name_, false, true); + + std::vector hwPt(ncands, 0), hwEta(ncands, 0), hwPhi(ncands, 0), + hwQual(ncands, 0), hwIso(ncands, 0), + hwDEta(ncands, 0), hwDPhi(ncands, 0), hwZ0(ncands, 0), + hwCharge(ncands, 0), hwIDScore(ncands, 0), + hwTkRedChi2RPhi(ncands, 0), hwTkCaloDphi(ncands, 0), + hwCaloShowerShape(ncands, 0), hwCaloTkPtRatio(ncands, 0); + + for (unsigned int i = 0; i < ncands; ++i) { + const auto& ele = *selected[i]; + if (ele.encoding() != l1t::TkEm::HWEncoding::CT) + continue; + + auto word = ele.egBinaryWord(); + l1ct::EGIsoEleObj ct; + ct.initFromBits(word); + + hwPt[i] = ct.hwPt.to_float(); + hwEta[i] = ct.hwEta.to_float(); + hwPhi[i] = ct.hwPhi.to_float(); + hwQual[i] = float(ct.hwQual); + hwIso[i] = ct.hwIso.to_float(); + hwDEta[i] = ct.hwDEta.to_float(); + hwDPhi[i] = ct.hwDPhi.to_float(); + hwZ0[i] = ct.hwZ0.to_float(); + hwCharge[i] = float(ct.hwCharge); + hwIDScore[i] = ct.hwIDScore.to_float(); + hwTkRedChi2RPhi[i] = ct.hwTkRedChi2RPhi.to_float(); + hwTkCaloDphi[i] = ct.hwTkCaloDphi.to_float(); + hwCaloShowerShape[i] = ct.hwCaloShowerShape.to_float(); + hwCaloTkPtRatio[i] = ct.hwCaloTkPtRatio.to_float(); + } + + out->addColumn("hwPt", hwPt, "CT hardware pT [LSB = 0.25 GeV]"); + out->addColumn("hwEta", hwEta, "CT hardware eta at calo face [LSB = pi/720]"); + out->addColumn("hwPhi", hwPhi, "CT hardware phi at calo face [LSB = pi/720]"); + out->addColumn("hwQual", hwQual, "CT hardware quality flags"); + out->addColumn("hwIso", hwIso, "CT hardware isolation pT [LSB = 0.25 GeV]"); + out->addColumn("hwDEta", hwDEta, "CT hardware track-calo delta eta [LSB = pi/720]"); + out->addColumn("hwDPhi", hwDPhi, "CT hardware track-calo delta phi [LSB = pi/720]"); + out->addColumn("hwZ0", hwZ0, "CT hardware track z0 [LSB = 0.05 cm]"); + out->addColumn("hwCharge", hwCharge, "CT hardware charge (1=positive, 0=negative)"); + out->addColumn("hwIDScore", hwIDScore, "CT hardware ID score [-1, 1]"); + out->addColumn("hwTkRedChi2RPhi", hwTkRedChi2RPhi, "CT hardware track reduced chi2 in R-phi (4-bit bin)"); + out->addColumn("hwTkCaloDphi", hwTkCaloDphi, "CT hardware track-calo delta phi (7 bits)"); + out->addColumn("hwCaloShowerShape", hwCaloShowerShape, "CT hardware calo shower shape (6 bits)"); + out->addColumn("hwCaloTkPtRatio", hwCaloTkPtRatio, "CT hardware calo/track pT ratio (10 bits)"); + + iEvent.put(std::move(out)); +} + +//define this as a plug-in +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(L1PFTkEleTableProducer); diff --git a/NtupleProducer/python/runPerformanceNTuple.py b/NtupleProducer/python/runPerformanceNTuple.py index 17347dd6..5ad9019e 100644 --- a/NtupleProducer/python/runPerformanceNTuple.py +++ b/NtupleProducer/python/runPerformanceNTuple.py @@ -16,7 +16,7 @@ def LazyVar(expr, valtype, doc=None, precision=-1): process.MessageLogger.cerr.FwkReport.reportEvery = 1 process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring('file:inputs125X.root'), + fileNames = cms.untracked.vstring('file:inputs140X.root'), inputCommands = cms.untracked.vstring("keep *", "drop l1tPFClusters_*_*_*", "drop l1tPFTracks_*_*_*", @@ -52,6 +52,10 @@ def LazyVar(expr, valtype, doc=None, precision=-1): from L1Trigger.L1CaloTrigger.l1tPhase2GCTBarrelToCorrelatorLayer1Emulator_cfi import l1tPhase2GCTBarrelToCorrelatorLayer1Emulator process.l1tPhase2GCTBarrelToCorrelatorLayer1Emulator = l1tPhase2GCTBarrelToCorrelatorLayer1Emulator.clone() +from L1Trigger.L1CaloTrigger.l1tPhase2CaloToCorrelatorTM18_cfi import l1tPhase2CaloToCorrelatorTM18 +process.l1tPhase2CaloToCorrelatorTM18 = l1tPhase2CaloToCorrelatorTM18.clone() + + from L1Trigger.Phase2L1ParticleFlow.L1NNTauProducer_cff import l1tNNTauProducerPuppi process.l1tNNTauProducerPuppi = l1tNNTauProducerPuppi.clone() @@ -63,6 +67,7 @@ def LazyVar(expr, valtype, doc=None, precision=-1): process.l1tPhase2L1CaloEGammaEmulator, process.l1tPhase2CaloPFClusterEmulator, process.l1tPhase2GCTBarrelToCorrelatorLayer1Emulator, + process.l1tPhase2CaloToCorrelatorTM18, process.l1tSAMuonsGmt, process.l1tGTTInputProducer, process.l1tTrackSelectionProducer, @@ -253,7 +258,10 @@ def addCalib(): process.l1tLayer1BarrelRaw = process.l1tLayer1Barrel.clone( gctEmInputConversionParameters = process.l1tLayer1Barrel.gctEmInputConversionParameters.clone( gctEmCorrector = cms.string("") - ) + ), + gctHadInputConversionParameters = process.l1tLayer1Barrel.gctHadInputConversionParameters.clone( + gctHadCorrector = cms.string("") + ), ) @@ -276,7 +284,7 @@ def addCalib(): #uncalibrated process.ntuple.objects.L1RawBarrelEcal = cms.VInputTag('l1tLayer1BarrelRaw:DecodedEmClusters') - process.ntuple.objects.L1RawBarrelCalo = cms.VInputTag('l1tPFClustersFromCombinedCaloHCal:uncalibrated') + process.ntuple.objects.L1RawBarrelCalo = cms.VInputTag('ll1tLayer1BarrelRaw:DecodedHadClusters') process.ntuple.objects.L1RawHGCal = cms.VInputTag('l1tLayer1HGCalRaw:DecodedHadClusters', 'l1tLayer1HGCalNoTKRaw:DecodedHadClusters')#use only this, try to understand if you have to use emf or emf_tot process.ntuple.objects.L1RawHGCalEM = cms.VInputTag('l1tLayer1HGCalRaw:DecodedEmClusters', 'l1tLayer1HGCalNoTKRaw:DecodedEmClusters') process.ntuple.objects.L1RawHFCalo = cms.VInputTag('l1tPFClustersFromCombinedCaloHF:uncalibrated') @@ -378,6 +386,7 @@ def addGen(pdgs): vz = Var("vz", float,precision=8), charge = Var("charge", int, doc="charge id"), prompt = Var("2*statusFlags().isPrompt() + statusFlags().isDirectPromptTauDecayProduct()", int, doc="Particle status."), + pdgId = Var("pdgId", int, doc="PDG id") ) ) genLepTableExt = cms.EDProducer("L1PFGenTableProducer", @@ -420,6 +429,17 @@ def addGen(pdgs): name = process.genPiTable.name ) process.extraPFStuff.add(process.genPiTable, process.genPiExtTable) + elif pdgId == 130: + process.genK0LTable = genLepTable.clone( + cut = cms.string("abs(pdgId) == %d && status == 1 && pt > 2" % pdgId), + name = cms.string("GenK0L")) + process.genK0LExtTable = genLepTableExt.clone( + cut = process.genK0LTable.cut, + name = process.genK0LTable.name + ) + process.extraPFStuff.add(process.genK0LTable, process.genK0LExtTable) + else: + raise ValueError("pdgId %d not supported in addGen" % pdgId) def addGenPi(pdgs=[211]): addGen(pdgs) @@ -627,19 +647,27 @@ def getTkEgTables(slice, postfix, tkem_inputtag, tkele_inputtag): tkEleTable.variables.caloEta = LazyVar("egCaloPtr.eta", float,precision=8) tkEleTable.variables.caloPhi = LazyVar("egCaloPtr.phi", float,precision=8) - return tkEmTable, tkEleTable + tkEleTableExt = cms.EDProducer("L1PFTkEleTableProducer", + src = cms.InputTag(tkele_inputtag), + name = cms.string("TkEle"+slice+postfix), + cut = cms.string(""),) + + return tkEmTable, tkEleTable, tkEleTableExt if doL1: for w in "EB","EE": - tkEmTable, tkEleTable = getTkEgTables(w, postfix, f"l1tLayer1EG{postfix}:L1TkEm{w}", f'l1tLayer1EG{postfix}:L1TkEle{w}') + tkEmTable, tkEleTable, tkEleTableExt = getTkEgTables(w, postfix, f"l1tLayer1EG{postfix}:L1TkEm{w}", f'l1tLayer1EG{postfix}:L1TkEle{w}') setattr(process, "TkEm%s%sTable" % (w,postfix), tkEmTable) setattr(process, "TkEle%s%sTable" % (w,postfix), tkEleTable) - process.extraPFStuff.add(tkEmTable,tkEleTable) + setattr(process, "TkEle%s%sExtTable" % (w,postfix), tkEleTableExt) + + process.extraPFStuff.add(tkEmTable,tkEleTable,tkEleTableExt) if doL2: - tkEmTable, tkEleTable = getTkEgTables('L2', postfix, f"l1tLayer2EG:L1CtTkEm", f'l1tLayer2EG:L1CtTkElectron') + tkEmTable, tkEleTable, tkEleTableExt = getTkEgTables('L2', postfix, f"l1tLayer2EG:L1CtTkEm", f'l1tLayer2EG:L1CtTkElectron') setattr(process, "TkEmL2%sTable" % (postfix), tkEmTable) setattr(process, "TkEleL2%sTable" % (postfix), tkEleTable) + # setattr(process, "TkEleL2%sExtTable" % (postfix), tkEleTableExt) process.extraPFStuff.add(tkEmTable,tkEleTable) @@ -732,7 +760,7 @@ def addDecodedCalo(types=['Had', 'Em'], regs=['HGCal','Barrel','HGCalNoTK']): hwPhi = LazyVar("hwPhi", int, doc="hwPhi"), ) ) - + decCaloTableExt = cms.EDProducer("L1PFDecodedCaloTableProducer", src = cms.InputTag("l1tLayer1"+reg, f'Decoded{tp}Clusters'), name = cms.string(""), @@ -842,3 +870,6 @@ def saveGenCands(): ), ) process.p += process.gencandTable + + +# addDecodedCalo() \ No newline at end of file