From f3f32cac95d829f790e805705dde297273abefa7 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 24 Apr 2017 12:26:38 +0200 Subject: [PATCH 1/7] Fix memory leak --- NtupleProducer/plugins/NtupleProducer.cc | 29 ++++++++++++------------ 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/NtupleProducer/plugins/NtupleProducer.cc b/NtupleProducer/plugins/NtupleProducer.cc index 71b101a9e430e..3426f0bf1a977 100644 --- a/NtupleProducer/plugins/NtupleProducer.cc +++ b/NtupleProducer/plugins/NtupleProducer.cc @@ -334,16 +334,16 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { iEvent.getByToken(TokEcalTPTag_, classicecals); edm::Handle hgecals; iEvent.getByToken(TokHGEcalTPTag_, hgecals); - L1PFCollection * ecals = new L1PFCollection; - for (const l1tpf::Particle & it : *classicecals) ecals->push_back(it); - for (const l1tpf::Particle & it : *hgecals ) ecals->push_back(it); - for (const l1tpf::Particle & it : *ecals) { + L1PFCollection ecals; + for (const l1tpf::Particle & it : *classicecals) ecals.push_back(it); + for (const l1tpf::Particle & it : *hgecals ) ecals.push_back(it); + for (const l1tpf::Particle & it : ecals) { lEEt[it.aEta()][it.aPhi()] = it.pt(); } std::vector pClust; int ne = 0; - for (const l1tpf::Particle & it : *ecals) { + for (const l1tpf::Particle & it : ecals) { double et = it.pt(); ecal_subdet = 0;// it->id().subDet(); // FIXME: missing in Particle class ecal_ieta = it.iEta(); @@ -397,13 +397,13 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { iEvent.getByToken(TokBHHcalTPTag_, bhhcals); edm::Handle hfhcals; iEvent.getByToken(TokHFTPTag_, hfhcals); - L1PFCollection * hcals = new L1PFCollection; - for (const l1tpf::Particle & it : *classichcals) hcals->push_back(it); - for (const l1tpf::Particle & it : *hfhcals ) hcals->push_back(it); + L1PFCollection hcals; + for (const l1tpf::Particle & it : *classichcals) hcals.push_back(it); + for (const l1tpf::Particle & it : *hfhcals ) hcals.push_back(it); for (const l1tpf::Particle & it : *hghcals ) { //combine duplicates for towers we use bool lFound = false; - for (l1tpf::Particle & it2 : *hcals ) { + for (l1tpf::Particle & it2 : hcals ) { if(it2.eta() == it.eta() && it2.phi() == it.phi() && it2.pt() == it.pt()) {std::cout << "copied!!!!!!!!!!!!!!!!!" << std::endl; break;} if(it2.iEta() != it.iEta() || it2.iPhi() != it.iPhi()) continue; TLorentzVector lVec = it.tp4(); @@ -412,12 +412,12 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { lFound = true; break; } - if(!lFound) hcals->push_back(it); + if(!lFound) hcals.push_back(it); } for (const l1tpf::Particle & it : *bhhcals ) { //Flatten the BH to the HG bool lFound = false; - for (l1tpf::Particle & it2 : *hcals) { + for (l1tpf::Particle & it2 : hcals) { if(it2.iPhi() == it.iPhi() && it2.iEta() == it.iEta()) { TLorentzVector lVec = it.tp4(); lVec += it2.tp4(); @@ -426,11 +426,12 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { break; } } - if(!lFound) hcals->push_back(it); + if(!lFound) hcals.push_back(it); } - for (const l1tpf::Particle & it : *hcals) lHEt[it.aEta()][it.aPhi()] += it.pt(); + for (const l1tpf::Particle & it : hcals) lHEt[it.aEta()][it.aPhi()] += it.pt(); + unsigned nh = 0; - for (const l1tpf::Particle & it : *hcals) { + for (const l1tpf::Particle & it : hcals) { hcal_subdet = 0; // it->id().subdet(); // FIXME: missing in Particle class hcal_ieta = it.iEta(); hcal_iphi = it.iPhi(); From f0f2c12bcf1e74a1c07c6187c446e047448b27ae Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 24 Apr 2017 12:29:16 +0200 Subject: [PATCH 2/7] New calo clustering code (without energy corrections for the moment) --- NtupleProducer/BuildFile.xml | 2 +- NtupleProducer/interface/CaloClusterer.h | 184 ++++++++++++ NtupleProducer/plugins/BuildFile.xml | 1 + NtupleProducer/src/CaloClusterer.cc | 340 +++++++++++++++++++++++ 4 files changed, 526 insertions(+), 1 deletion(-) create mode 100644 NtupleProducer/interface/CaloClusterer.h create mode 100644 NtupleProducer/src/CaloClusterer.cc diff --git a/NtupleProducer/BuildFile.xml b/NtupleProducer/BuildFile.xml index 5602bb9931a25..97f033001a362 100644 --- a/NtupleProducer/BuildFile.xml +++ b/NtupleProducer/BuildFile.xml @@ -1,4 +1,4 @@ - + diff --git a/NtupleProducer/interface/CaloClusterer.h b/NtupleProducer/interface/CaloClusterer.h new file mode 100644 index 0000000000000..dd474958d34f2 --- /dev/null +++ b/NtupleProducer/interface/CaloClusterer.h @@ -0,0 +1,184 @@ +#ifndef FASTPUPPI_NTUPLERPRODUCER_CALOCLUSTERER_H +#define FASTPUPPI_NTUPLERPRODUCER_CALOCLUSTERER_H +/** + * Classes for new calorimetric clustering + * + * */ + +// fwd declarations +namespace edm { class ParameterSet; } + +// real includes +#include +#include +#include +#include +#include +#include +#include "L1TPFParticle.h" + +namespace l1pf_calo { + class Grid { + public: + virtual ~Grid() {} + unsigned int size() const { return ncells_; } + virtual int find_cell(float eta, float phi) const = 0; + int neighbour(int icell, unsigned int idx) const { return neighbours_[icell][idx]; } + float eta(int icell) const { return eta_[icell]; } + float phi(int icell) const { return phi_[icell]; } + float etaWidth(int icell) const { return etaWidth_[icell]; } + float phiWidth(int icell) const { return phiWidth_[icell]; } + virtual int ieta(int icell) const { return 0; } + virtual int iphi(int icell) const { return 0; } + protected: + Grid(unsigned int size) : ncells_(size), eta_(size), etaWidth_(size), phi_(size), phiWidth_(size), neighbours_(size) {} + unsigned int ncells_; + std::vector eta_, etaWidth_, phi_, phiWidth_; + std::vector> neighbours_; // indices of the neigbours, -1 = none + }; + + class Stage1Grid : public Grid { + public: + Stage1Grid() ; + virtual int ieta(int icell) const override { return ieta_[icell]; } + virtual int iphi(int icell) const override { return iphi_[icell]; } + virtual int find_cell(float eta, float phi) const override ; + int ifind_cell(int ieta, int iphi) const { return cell_map_[(ieta+nEta_) + 2*nEta_*(iphi-1)]; } + protected: + static const int nEta_ = 41, nPhi_ = 72, ietaCoarse_ = 29, ietaVeryCoarse_ = 40; + static const float towerEtas_[nEta_]; + std::vector ieta_, iphi_; + std::vector cell_map_; + // valid ieta, iphi (does not check for outside bounds, only for non-existence of ieta=0, iphi=0, and coarser towers at high eta + bool valid_ieta_iphi(int ieta, int iphi) const { + if (ieta == 0 || iphi == 0) return false; + if (std::abs(ieta) >= ietaVeryCoarse_ && (iphi % 4 != 1)) return false; + if (std::abs(ieta) >= ietaCoarse_ && (iphi % 2 != 1)) return false; + return true; + } + // move by +/-1 around a cell; return icell or -1 if not available + int imove(int ieta, int iphi, int deta, int dphi) ; + + }; + + template + class GridData { + public: + GridData() : grid_(nullptr), data_(), empty_() {} + GridData(const Grid &grid) : grid_(&grid), data_(grid.size()), empty_() {} + + T & operator()(float eta, float phi) { return data_[grid_->find_cell(eta,phi)]; } + const T & operator()(float eta, float phi) const { return data_[grid_->find_cell(eta,phi)]; } + + const Grid & grid() const { return *grid_; } + + unsigned int size() const { return data_.size(); } + + float eta(int icell) const { return grid().eta(icell); } + float phi(int icell) const { return grid().phi(icell); } + int ieta(int icell) const { return grid().ieta(icell); } + int iphi(int icell) const { return grid().iphi(icell); } + + T & operator[](int icell) { return data_[icell]; } + const T & operator[](int icell) const { return data_[icell]; } + + const T & neigh(int icell, unsigned int idx) const { + int ineigh = grid_->neighbour(icell, idx); + return (ineigh < 0 ? empty_ : data_[ineigh]); + } + + // always defined + void fill(const T &val) { std::fill(data_.begin(), data_.end(), val); } + void zero() { fill(T()); } + + // defined only if T has a 'clear' method + void clear() { for (T & t : data_) t.clear(); } + + private: + const Grid * grid_; + std::vector data_; + const T empty_; + }; + typedef GridData EtGrid; + + struct PreCluster { + PreCluster() : ptLocalMax(0), ptOverNeighLocalMaxSum(0) {} + float ptLocalMax; // pt if it's a local max, zero otherwise + float ptOverNeighLocalMaxSum; // pt / (sum of ptLocalMax of neighbours); zero if no neighbours + void clear() { ptLocalMax = ptOverNeighLocalMaxSum = 0; } + }; + typedef GridData PreClusterGrid; + + struct Cluster { + Cluster() : et(0), et_corr(0), eta(0), phi(0) {} + float et, et_corr, eta, phi; + void clear() { et = et_corr = eta = phi = 0; } + }; + typedef GridData ClusterGrid; + + struct CombinedCluster : public Cluster { + float ecal_et, hcal_et; + void clear() { et = et_corr = ecal_et = hcal_et = eta = phi = 0; } + }; + typedef GridData CombinedClusterGrid; + + class SingleCaloClusterer { + public: + SingleCaloClusterer(const edm::ParameterSet &pset) ; + ~SingleCaloClusterer() ; + void clear() { rawet_.zero(); } + void add(const l1tpf::Particle &particle) { + if (particle.pt() > 0) { + rawet_(particle.eta(), particle.phi()) += particle.pt(); + } + } + void run() ; + const EtGrid & raw() const { return rawet_; } + const ClusterGrid & clusters() const { return cluster_; } + + // for the moment, generic interface that takes a cluster and corrects it + template + void correct(const Corrector & corrector) { + for (unsigned int i = 0, ncells = grid_->size(); i < ncells; ++i) { + if (cluster_[i].et > 0) { + cluster_[i].et_corr = corrector(cluster_[i], grid_->ieta(i), grid_->iphi(i)); + } + } + } + private: + std::unique_ptr grid_; + EtGrid rawet_; + PreClusterGrid precluster_; + ClusterGrid cluster_; + float zsEt_, seedEt_, minClusterEt_; + }; + + class SimpleCaloLinker { + public: + SimpleCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer & ecal, const SingleCaloClusterer & hcal) ; + ~SimpleCaloLinker() ; + void run() ; + const CombinedClusterGrid & clusters() const { return cluster_; } + + // for the moment, generic interface that takes a cluster and corrects it + template + void correct(const Corrector & corrector) { + for (unsigned int i = 0, ncells = grid_->size(); i < ncells; ++i) { + if (cluster_[i].et > 0) { + cluster_[i].et_corr = corrector(cluster_[i], grid_->ieta(i), grid_->iphi(i)); + } + } + } + + std::vector fetch(bool corrected=true) const ; + private: + std::unique_ptr grid_; + const SingleCaloClusterer & ecal_, & hcal_; + PreClusterGrid ecalToHCal_; + CombinedClusterGrid cluster_; + float hoeCut_, minPhotonEt_, minHadronEt_; + }; + +} // namespace + +#endif diff --git a/NtupleProducer/plugins/BuildFile.xml b/NtupleProducer/plugins/BuildFile.xml index 759c6413898df..0acc75836cf80 100644 --- a/NtupleProducer/plugins/BuildFile.xml +++ b/NtupleProducer/plugins/BuildFile.xml @@ -21,4 +21,5 @@ + diff --git a/NtupleProducer/src/CaloClusterer.cc b/NtupleProducer/src/CaloClusterer.cc new file mode 100644 index 0000000000000..99f443571a5c4 --- /dev/null +++ b/NtupleProducer/src/CaloClusterer.cc @@ -0,0 +1,340 @@ +#include "../interface/CaloClusterer.h" +#include +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +const float l1pf_calo::Stage1Grid::towerEtas_[l1pf_calo::Stage1Grid::nEta_] = {0,0.087,0.174,0.261,0.348,0.435,0.522,0.609,0.696,0.783,0.870,0.957,1.044,1.131,1.218,1.305,1.392,1.479,1.566,1.653,1.740,1.830,1.930,2.043,2.172,2.322,2.5,2.650,2.853,3.139,3.314,3.489,3.664,3.839,4.013,4.191,4.363,4.538,4.716,4.889,5.191}; + +l1pf_calo::Stage1Grid::Stage1Grid() : + Grid(2*((ietaCoarse_-1)*nPhi_ + (ietaVeryCoarse_-ietaCoarse_)*(nPhi_/2) + (nEta_-ietaVeryCoarse_+1) * (nPhi_/4))), + ieta_(ncells_), iphi_(ncells_), + cell_map_(2*nEta_*nPhi_, -1) +{ + int icell = 0; + for (int ie = -nEta_; ie <= nEta_; ++ie) { + int absie = std::abs(ie); + for (int iph = 1; iph <= nPhi_; ++iph) { + if (!valid_ieta_iphi(ie,iph)) continue; + ieta_[icell] = ie; + iphi_[icell] = iph; + eta_[icell] = (ie > 0 ? 0.5 : -0.5)*(towerEtas_[absie-1] + towerEtas_[absie]); + etaWidth_[icell] = (towerEtas_[absie] - towerEtas_[absie-1]); + phiWidth_[icell] = 2*M_PI/nPhi_; + if (absie >= ietaVeryCoarse_) phiWidth_[icell] *= 4; + else if (absie >= ietaCoarse_) phiWidth_[icell] *= 2; + phi_[icell] = (iph-1)*2*M_PI/nPhi_ + 0.5*phiWidth_[icell]; + if (phi_[icell] > M_PI) phi_[icell] -= 2*M_PI; + std::fill(neighbours_[icell].begin(), neighbours_[icell].end(), -1); + cell_map_[(ie+nEta_) + 2 * nEta_*(iph-1)] = icell; + icell++; + } + } + assert(unsigned(icell) == ncells_); + // now link the cells + for (icell = 0; icell < int(ncells_); ++icell) { + int ie = ieta_[icell], iph = iphi_[icell]; + int ineigh = 0; + for (int deta = -1; deta <= +1; ++deta) { + for (int dphi = -1; dphi <= +1; ++dphi) { + if (deta == 0 && dphi == 0) continue; + neighbours_[icell][ineigh++] = imove(ie, iph, deta, dphi); + } + } + } + // consistency check 1: check that find_cell works + for (float teta = 0; teta <= 5.0; teta += 0.02) { + for (float tphi = -M_PI; tphi <= M_PI; tphi += 0.02) { + find_cell(+teta, tphi); + find_cell(-teta, tphi); + } + } + // consistency check: neighbours for cells at small abs(ieta) + for (int aie = 1; aie < ietaCoarse_-1; ++aie) { + for (int is = -1; is <= +1; is += 2) { + int ie = aie*is; + for (int iph = 1; iph <= nPhi_; ++iph) { + int icell = ifind_cell(ie,iph); + int nvalid = 0; + for (int in = 0; in < 8; ++in) { + int ineigh = neighbour(icell, in); + if (ineigh != -1) { + if ((eta(icell) == eta(ineigh)) || + (std::abs(std::abs(eta(icell)-eta(ineigh)) - 0.5*(etaWidth(icell)+etaWidth(ineigh))) < 0.01)) { + if ((phi(icell) == phi(ineigh)) || + (std::abs(std::abs(deltaPhi(phi(icell),phi(ineigh))) - 0.5*(phiWidth(icell)+phiWidth(ineigh))) < 0.01)) { + nvalid++; + } else { + printf("Error in neighbour cell for ieta %+3d iphi %2d eta %+7.4f +- %.4f phi %+7.4f +- %.4f, got neigh ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f\n", + ie, iph, eta(icell), etaWidth(icell), phi(icell), phiWidth(icell), + ieta(ineigh), iphi(ineigh), eta(ineigh), etaWidth(ineigh), phi(ineigh), phiWidth(ineigh)); + } + } else { + printf("Error in neighbour cell for ieta %+3d iphi %2d eta %+7.4f +- %.4f phi %+7.4f +- %.4f, got neigh ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f\n", + ie, iph, eta(icell), etaWidth(icell), phi(icell), phiWidth(icell), + ieta(ineigh), iphi(ineigh), eta(ineigh), etaWidth(ineigh), phi(ineigh), phiWidth(ineigh)); + } + } else { + printf("Mismatch cell ieta %+3d iphi %2d neigh %d is null\n", ie,iph,in); + } + } + } + } + } +} + +int l1pf_calo::Stage1Grid::find_cell(float eta, float phi) const { + int ieta = (eta != 0) ? std::distance(towerEtas_, std::lower_bound(towerEtas_, towerEtas_+nEta_, std::abs(eta))) : 1; + if (!(ieta > 0 && ieta < nEta_)) { + printf("Error in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d which is not valid\n", eta, phi, ieta); + } + assert(ieta > 0 && ieta < nEta_); + if (ieta > nEta_) ieta = nEta_; + if (eta < 0) ieta = -ieta; + if (phi > 2*M_PI) phi -= 2*M_PI; + if (phi < 0) phi += 2*M_PI; + int iphi = std::floor(phi * nPhi_/(2*M_PI)); + if (phi >= 2*M_PI) iphi = nPhi_-1; // fix corner case due to roundings etc + assert(iphi < nPhi_); + if (std::abs(ieta) >= ietaVeryCoarse_) iphi -= (iphi%4); + else if (std::abs(ieta) >= ietaCoarse_) iphi -= (iphi%2); + iphi += 1; + if (!valid_ieta_iphi(ieta,iphi)) { + printf("Error in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which is not valid\n", + eta, phi, ieta, iphi); + } + assert(valid_ieta_iphi(ieta,iphi)); + int icell = ifind_cell(ieta,iphi); + assert(icell != -1); + + if (std::abs(eta - eta_[icell]) > 0.501*etaWidth_[icell] || std::abs(deltaPhi(phi, phi_[icell])) > 0.501*phiWidth_[icell]) { + printf("Mismatch in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f ; deta = %+7.4f dphi = %+7.4f\n", + eta, phi, ieta, iphi, eta_[icell], etaWidth_[icell], phi_[icell], phiWidth_[icell], eta - eta_[icell], deltaPhi(phi, phi_[icell])); + } + //assert(std::abs(eta - eta_[icell]) <= 0.5*etaWidth_[icell]); + //assert(std::abs(deltaPhi(phi, phi_[icell])) <= 0.5*phiWidth_[icell]); + return icell; +} + +int l1pf_calo::Stage1Grid::imove(int ieta, int iphi, int deta, int dphi) { + int ie = ieta, iph = iphi; + switch (deta) { + case -1: ie = (ie == -nEta_ ? 0 : (ie == +1 ? -1 : ie-1)); break; + case +1: ie = (ie == +nEta_ ? 0 : (ie == -1 ? +1 : ie+1)); break; + case 0: break; + default: assert(false); + }; + if (ie == 0) return -1; + switch (dphi) { + case -1: iph = (iph == 1 ? nPhi_ : iph-1); break; + case +1: iph = (iph == nPhi_ ? 1 : iph+1); break; + case 0: break; + default: assert(false); + }; + if (!valid_ieta_iphi(ie,iph)) return -1; + int icell = ifind_cell(ie,iph); + if ((ie == ieta && iph == iphi)) { + printf("Error bogus move %+3d %2d + (%+1d %+1d) = %+3d %2d \n", ieta, iphi, deta, dphi, ie, iph); + } + assert(!(ie == ieta && iph == iphi)); + assert(icell != -1); + assert(icell != ifind_cell(ieta,iphi)); + return icell; +} + +l1pf_calo::SingleCaloClusterer::SingleCaloClusterer(const edm::ParameterSet &pset) : + grid_(new Stage1Grid()), + rawet_(*grid_), + precluster_(*grid_), + cluster_(*grid_), + zsEt_(pset.getParameter("zsEt")), + seedEt_(pset.getParameter("seedEt")), + minClusterEt_(pset.getParameter("minClusterEt")) +{ +} + +l1pf_calo::SingleCaloClusterer::~SingleCaloClusterer() +{ +} + +void l1pf_calo::SingleCaloClusterer::run() +{ + unsigned int i, ncells = grid_->size(); + + // kill zeros + for (i = 0; i < ncells; ++i) { + if (rawet_[i] < zsEt_) rawet_[i] = 0; + } + + precluster_.clear(); + // pre-cluster step 1: at each cell, set the value equal to itself if it's a local maxima, zero otherwise + // can be done in parallel on all cells + for (i = 0; i < ncells; ++i) { + if (rawet_[i] > seedEt_) { + precluster_[i].ptLocalMax = rawet_[i]; + //printf(" candidate precluster pt %7.2f at %4d (ieta %+3d iphi %2d)\n", rawet_[i], i, grid_->ieta(i), grid_->iphi(i)); + for (int ineigh = 0; ineigh <= 3; ++ineigh) { + if (rawet_.neigh(i,ineigh) > rawet_[i]) precluster_[i].ptLocalMax = 0; + //int ncell = grid_->neighbour(i,ineigh); + //if (ncell == -1) printf(" \t neigh %d is null\n", ineigh); + //else printf(" \t neigh %d at %4d (ieta %+3d iphi %2d) has pt %7.2f: comparison %1d \n", ineigh, ncell, grid_->ieta(ncell), grid_->iphi(ncell), rawet_[ncell], precluster_[i].ptLocalMax > 0); + } + for (int ineigh = 4; ineigh < 8; ++ineigh) { + if (rawet_.neigh(i,ineigh) >= rawet_[i]) precluster_[i].ptLocalMax = 0; + //int ncell = grid_->neighbour(i,ineigh); + //if (ncell == -1) printf(" \t neigh %d is null\n", ineigh); + //else printf(" \t neigh %d at %4d (ieta %+3d iphi %2d) has pt %7.2f: comparison %1d \n", ineigh, ncell, grid_->ieta(ncell), grid_->iphi(ncell), rawet_[ncell], precluster_[i].ptLocalMax > 0); + } + } + } + // pre-cluster step 2: at each non-max cell, add up ptLocalMax of neighbours + // could be done for all cells, but it's not needed for max cells since they can't have ptLocalMax around them + for (i = 0; i < ncells; ++i) { + if (precluster_[i].ptLocalMax == 0) { + float tot = 0; + for (int ineigh = 0; ineigh < 8; ++ineigh) { + tot += precluster_.neigh(i,ineigh).ptLocalMax; + } + precluster_[i].ptOverNeighLocalMaxSum = tot ? rawet_[i]/tot : 0; + } + } + + cluster_.clear(); + // cluster: at each localMax cell, take itself plus the weighted contributions of the neighbours + for (i = 0; i < ncells; ++i) { + if (precluster_[i].ptLocalMax > 0) { + float myet = rawet_[i]; + float tot = myet; + float avg_eta = 0; + float avg_phi = 0; + for (int ineigh = 0; ineigh < 8; ++ineigh) { + int ineighcell = grid_->neighbour(i, ineigh); + if (ineighcell == -1) continue; // skip dummy cells + float fracet = myet * precluster_.neigh(i,ineigh).ptOverNeighLocalMaxSum; + tot += fracet; + avg_eta += fracet * (grid_->eta(ineighcell) - grid_->eta(i)); + avg_phi += fracet * deltaPhi(grid_->phi(ineighcell), grid_->phi(i)); + } + if (tot > minClusterEt_) { + cluster_[i].et = tot; + cluster_[i].eta = grid_->eta(i) + avg_eta / tot; + cluster_[i].phi = grid_->phi(i) + avg_phi / tot; + // wrap around phi + if (cluster_[i].phi > M_PI) cluster_[i].phi -= 2*M_PI; + if (cluster_[i].phi < -M_PI) cluster_[i].phi += 2*M_PI; + } + } + } + + +} + +l1pf_calo::SimpleCaloLinker::SimpleCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer & ecal, const SingleCaloClusterer & hcal) : + grid_(new Stage1Grid()), + ecal_(ecal), hcal_(hcal), + ecalToHCal_(*grid_), + cluster_(*grid_), + hoeCut_(pset.getParameter("hoeCut")), + minPhotonEt_(pset.getParameter("minPhotonEt")), + minHadronEt_(pset.getParameter("minHadronEt")) +{ + assert(grid_->size() == ecal.raw().grid().size()); + assert(grid_->size() == hcal.raw().grid().size()); +} + +l1pf_calo::SimpleCaloLinker::~SimpleCaloLinker() +{ +} + +void l1pf_calo::SimpleCaloLinker::run() +{ + unsigned int i, ncells = grid_->size(); + + const EtGrid & hraw = hcal_.raw(); + const ClusterGrid & ecals = ecal_.clusters(); + const ClusterGrid & hcals = hcal_.clusters(); + + // for each ECal cluster, get the corresponding HCal cluster and the sum of the neighbour HCal clusters + ecalToHCal_.clear(); + for (i = 0; i < ncells; ++i) { + if (ecals[i].et > 0) { + if (hcals[i].et > 0) { + ecalToHCal_[i].ptLocalMax = hcals[i].et; + } else { + float tot = 0; + for (int ineigh = 0; ineigh < 8; ++ineigh) { + tot += hcals.neigh(i,ineigh).et; + } + ecalToHCal_[i].ptOverNeighLocalMaxSum = tot ? ecals[i].et/tot : 0; + } + } + } + + cluster_.clear(); + // promote HCal clusters to final clusters + for (i = 0; i < ncells; ++i) { + if (hcals[i].et > 0) { + if (ecalToHCal_[i].ptLocalMax > 0) { + // direct linking is easy + cluster_[i].ecal_et = ecals[i].et; + cluster_[i].hcal_et = hcals[i].et; + cluster_[i].et = hcals[i].et + ecals[i].et; + float wecal = ecals[i].et/cluster_[i].et, whcal = 1.0 - wecal; + cluster_[i].eta = ecals[i].eta * wecal + hcals[i].eta * whcal; + cluster_[i].phi = ecals[i].phi * wecal + hcals[i].phi * whcal; + // wrap around phi + if (cluster_[i].phi > M_PI) cluster_[i].phi -= 2*M_PI; + if (cluster_[i].phi < -M_PI) cluster_[i].phi += 2*M_PI; + } else { + // sidewas linking is more annonying + float myet = hcals[i].et; + float etot = 0; + float avg_eta = 0; + float avg_phi = 0; + for (int ineigh = 0; ineigh < 8; ++ineigh) { + int ineighcell = grid_->neighbour(i, ineigh); + if (ineighcell == -1) continue; // skip dummy cells + float fracet = myet * ecalToHCal_.neigh(i,ineigh).ptOverNeighLocalMaxSum; + etot += fracet; + avg_eta += fracet * (grid_->eta(ineighcell) - grid_->eta(i)); + avg_phi += fracet * deltaPhi(grid_->phi(ineighcell), grid_->phi(i)); + } + cluster_[i].hcal_et = hcals[i].et; + cluster_[i].ecal_et = etot; + cluster_[i].et = myet + etot; + cluster_[i].eta = hcals[i].eta + avg_eta / cluster_[i].et; + cluster_[i].phi = hcals[i].phi + avg_phi / cluster_[i].et; + // wrap around phi + if (cluster_[i].phi > M_PI) cluster_[i].phi -= 2*M_PI; + if (cluster_[i].phi < -M_PI) cluster_[i].phi += 2*M_PI; + } + } + } + + // promote Unlinked ECal clusters to final clusters + for (i = 0; i < ncells; ++i) { + if (ecals[i].et > 0 && ecalToHCal_[i].ptLocalMax == 0 && ecalToHCal_[i].ptOverNeighLocalMaxSum == 0) { + // direct linking is easy + cluster_[i].ecal_et = ecals[i].et; + cluster_[i].hcal_et = hraw[i]; + cluster_[i].et = hraw[i] + ecals[i].et; + cluster_[i].eta = ecals[i].eta; + cluster_[i].phi = ecals[i].phi; + // no need to wrap around phi + } + } + +} + + + +std::vector l1pf_calo::SimpleCaloLinker::fetch(bool corrected) const { + std::vector ret; + for (unsigned int i = 0, ncells = grid_->size(); i < ncells; ++i) { + if (cluster_[i].et > 0) { + bool photon = (cluster_[i].hcal_et < hoeCut_* cluster_[i].ecal_et); + if (cluster_[i].et > (photon ? minPhotonEt_ : minHadronEt_)) { + ret.emplace_back(corrected ? cluster_[i].et_corr : cluster_[i].et, cluster_[i].eta, cluster_[i].phi, 0.0, photon ? 3 : 2); + } + } + } + return ret; +} From 1778cf9aa0bd4b311d3179f20dc49d4017f93079 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 24 Apr 2017 16:14:47 +0200 Subject: [PATCH 3/7] Cmssw module for more response studies --- NtupleProducer/plugins/BuildFile.xml | 1 + NtupleProducer/plugins/ResponseNTuplizer.cc | 289 ++++++++++++++++++++ NtupleProducer/python/runRespNTupler.py | 43 +++ 3 files changed, 333 insertions(+) create mode 100644 NtupleProducer/plugins/ResponseNTuplizer.cc create mode 100644 NtupleProducer/python/runRespNTupler.py diff --git a/NtupleProducer/plugins/BuildFile.xml b/NtupleProducer/plugins/BuildFile.xml index 0acc75836cf80..7110d784e45d8 100644 --- a/NtupleProducer/plugins/BuildFile.xml +++ b/NtupleProducer/plugins/BuildFile.xml @@ -21,5 +21,6 @@ + diff --git a/NtupleProducer/plugins/ResponseNTuplizer.cc b/NtupleProducer/plugins/ResponseNTuplizer.cc new file mode 100644 index 0000000000000..254591f950395 --- /dev/null +++ b/NtupleProducer/plugins/ResponseNTuplizer.cc @@ -0,0 +1,289 @@ +// -*- C++ -*- +// +// Package: Giovanni/NTuplizer +// Class: NTuplizer +// +/**\class NTuplizer NTuplizer.cc Giovanni/NTuplizer/plugins/NTuplizer.cc + + Description: [one line class summary] + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Giovanni Petrucciani +// Created: Thu, 01 Sep 2016 11:30:38 GMT +// +// + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Common/interface/View.h" + +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/JetReco/interface/GenJet.h" + +#include "DataFormats/Math/interface/deltaR.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +//#include "CommonTools/Utils/interface/StringCutObjectSelector.h" +//#include "CommonTools/Utils/interface/StringObjectFunction.h" + +#include + + +namespace { + struct SimpleObject { + float pt, eta, phi; + SimpleObject(float apt, float aneta, float aphi) : pt(apt), eta(aneta), phi(aphi) {} + bool operator<(const SimpleObject &other) const { return eta < other.eta; } + bool operator<(const float &other) const { return eta < other; } + }; + class MultiCollection { + public: + MultiCollection(const edm::ParameterSet &iConfig, const std::string &name, edm::ConsumesCollector && coll) : + name_(name) + { + const std::vector & tags = iConfig.getParameter< std::vector>(name); + for (const auto & tag : tags) tokens_.push_back(coll.consumes(tag)); + } + const std::string & name() const { return name_; } + void get(const edm::Event &iEvent) { + edm::Handle handle; + for (const auto & token : tokens_) { + iEvent.getByToken(token, handle); + for (const reco::Candidate & c : *handle) { + objects_.emplace_back(c.pt(), c.eta(), c.phi()); + } + } + std::sort(objects_.begin(), objects_.end()); + } + const std::vector & objects() const { return objects_; } + void clear() { objects_.clear(); } + private: + std::string name_; + std::vector> tokens_; + std::vector objects_; + }; + class InCone { + public: + InCone(const std::vector & objects, float eta, float phi, float dr) { + auto first = std::lower_bound(objects.begin(), objects.end(), eta-dr-0.01f); // small offset to avoid dealing with == + auto end = std::lower_bound(objects.begin(), objects.end(), eta+dr+0.01f); + float dr2 = dr*dr; + for (auto it = first; it < end; ++it) { + float mydr2 = ::deltaR2(eta,phi, it->eta,it->phi); + if (mydr2 < dr2) ptdr2.emplace_back(it->pt, mydr2); + } + } + float sum() const { + float mysum = 0; + for (const auto & p : ptdr2) mysum += p.first; + return mysum; + } + float sum(float dr) const { + float dr2 = dr*dr; + float mysum = 0; + for (const auto & p : ptdr2) { + if (p.second < dr2) mysum += p.first; + } + return mysum; + } + float nearest() const { + std::pair best(0,9999); + for (const auto & p : ptdr2) { + if (p.second < best.second) best = p; + } + return best.first; + } + float max() const { + float best = 0; + for (const auto & p : ptdr2) { + if (p.first > best) best = p.first; + } + return best; + } + + private: + std::vector> ptdr2; + }; + +} +class ResponseNTuplizer : public edm::one::EDAnalyzer { + public: + explicit ResponseNTuplizer(const edm::ParameterSet&); + ~ResponseNTuplizer(); + + private: + virtual void beginJob() override; + virtual void analyze(const edm::Event&, const edm::EventSetup&) override; + + edm::EDGetTokenT> genjets_; + edm::EDGetTokenT> genparticles_; + TTree *tree_; + struct McVars { + float pt, pt02, eta, phi, iso02, iso04; + int id; + void makeBranches(TTree *tree) { + tree->Branch("mc_pt", &pt, "mc_pt/F"); + tree->Branch("mc_pt02", &pt02, "mc_pt02/F"); + tree->Branch("mc_eta", &eta, "mc_eta/F"); + tree->Branch("mc_phi", &phi, "mc_phi/F"); + tree->Branch("mc_iso02", &iso02, "mc_iso02/F"); + tree->Branch("mc_iso04", &iso04, "mc_iso04/F"); + tree->Branch("mc_id", &id, "mc_id/I"); + } + void fillP4(const reco::Candidate &c) { + pt = c.pt(); eta = c.eta(); phi = c.phi(); + } + } mc_; + struct RecoVars { + float pt, pt02, ptbest, pthighest; + void makeBranches(const std::string &prefix, TTree *tree) { + tree->Branch((prefix+"_pt").c_str(), &pt, (prefix+"_pt/F").c_str()); + tree->Branch((prefix+"_pt02").c_str(), &pt02, (prefix+"_pt02/F").c_str()); + tree->Branch((prefix+"_ptbest").c_str(), &ptbest, (prefix+"_ptbest/F").c_str()); + tree->Branch((prefix+"_pthighest").c_str(), &pthighest, (prefix+"_pthighest/F").c_str()); + } + void fill(const std::vector<::SimpleObject> & objects, float eta, float phi) { + ::InCone incone(objects, eta, phi, 0.4); + pt = incone.sum(); + pt02 = incone.sum(0.2); + ptbest = incone.nearest(); + pthighest = incone.max(); + } + }; + std::vector> reco_; +}; + +ResponseNTuplizer::ResponseNTuplizer(const edm::ParameterSet& iConfig) : + genjets_(consumes>(iConfig.getParameter("genJets"))), + genparticles_(consumes>(iConfig.getParameter("genParticles"))) +{ + usesResource("TFileService"); + edm::Service fs; + tree_ = fs->make("tree","tree"); + + auto reconames = iConfig.getParameterNamesForType>(); + for (const std::string & name : reconames) { + reco_.emplace_back(::MultiCollection(iConfig,name,consumesCollector()),RecoVars()); + } +} + +ResponseNTuplizer::~ResponseNTuplizer() { } + +// ------------ method called once each job just before starting event loop ------------ +void +ResponseNTuplizer::beginJob() +{ + mc_.makeBranches(tree_); + for (auto & pair : reco_) { + pair.second.makeBranches(pair.first.name(), tree_); + } +} + + +// ------------ method called for each event ------------ +void +ResponseNTuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::Handle> genjets; + edm::Handle> genparticles; + iEvent.getByToken(genjets_, genjets); + iEvent.getByToken(genparticles_, genparticles); + + std::vector prompts, taus; + for (const reco::GenParticle &gen : *genparticles) { + if ((gen.isPromptFinalState() || gen.isDirectPromptTauDecayProductFinalState()) && (std::abs(gen.pdgId()) == 11 || std::abs(gen.pdgId()) == 13) && gen.pt() > 5) { + prompts.push_back(&gen); + } else if (gen.isPromptFinalState() && std::abs(gen.pdgId()) == 22 && gen.pt() > 10) { + prompts.push_back(&gen); + } else if (abs(gen.pdgId()) == 15 && gen.isPromptDecayed()) { + taus.push_back(&gen); + } + } + + for (auto & recopair : reco_) { + recopair.first.get(iEvent); + } + for (const reco::GenJet & j : *genjets) { + bool ok = true; + const reco::Candidate * match = nullptr; + for (const reco::GenParticle * p : prompts) { + if (::deltaR2(*p, j) < 0.16f) { + if (match != nullptr) { ok = false; break; } + else { match = p; } + } + } + if (!ok) continue; + if (!match) { + // look for a tau + for (const reco::GenParticle * p : taus) { + if (::deltaR2(*p, j) < 0.16f) { + if (match != nullptr) { ok = false; break; } + else { match = p; } + } + } + if (!ok) continue; + if (match != nullptr && match->numberOfDaughters() == 2 && std::abs(match->daughter(0)->pdgId()) + std::abs(match->daughter(1)->pdgId()) == 211+16) { + // one-prong tau, consider it a pion + match = (std::abs(match->daughter(0)->pdgId()) == 211 ? match->daughter(0) : match->daughter(1)); + } + } + if (match != nullptr) { + if (std::abs(match->pdgId()) == 15) { + reco::Particle::LorentzVector pvis; + for (unsigned int i = 0, n = match->numberOfDaughters(); i < n; ++i) { + const reco::Candidate *dau = match->daughter(i); + if (std::abs(dau->pdgId()) == 12 || std::abs(dau->pdgId()) == 14 || std::abs(dau->pdgId()) == 16) { + continue; + } + pvis += dau->p4(); + } + mc_.pt = pvis.Pt(); + mc_.eta = pvis.Eta(); + mc_.phi = pvis.Phi(); + } else { + mc_.fillP4(*match); + } + mc_.id = std::abs(match->pdgId()); + mc_.iso04 = j.pt()/mc_.pt - 1; + mc_.iso02 = 0; + for (const auto &dptr : j.daughterPtrVector()) { + if (::deltaR2(*dptr, *match) < 0.04f) { + mc_.iso02 += dptr->pt(); + } + } + mc_.iso02 = mc_.iso02/mc_.pt - 1; + } else { + if (j.pt() < 20) continue; + mc_.fillP4(j); + mc_.id = 0; + mc_.iso02 = 0; + mc_.iso04 = 0; + } + for (auto & recopair : reco_) { + recopair.second.fill(recopair.first.objects(), mc_.eta, mc_.phi); + } + tree_->Fill(); + } + for (auto & recopair : reco_) { + recopair.first.clear(); + } + +} + +//define this as a plug-in +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(ResponseNTuplizer); diff --git a/NtupleProducer/python/runRespNTupler.py b/NtupleProducer/python/runRespNTupler.py new file mode 100644 index 0000000000000..953f77a2fe981 --- /dev/null +++ b/NtupleProducer/python/runRespNTupler.py @@ -0,0 +1,43 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("RESP") + +process.load('Configuration.StandardSequences.Services_cff') +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(-1)) +process.MessageLogger.cerr.FwkReport.reportEvery = 1000 + +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('file:l1pf_out.root') +) + +process.ntuple = cms.EDAnalyzer("ResponseNTuplizer", + genJets = cms.InputTag("ak4GenJetsNoNu"), + genParticles = cms.InputTag("genParticles"), + # -- inputs -- + Ecal = cms.VInputTag('l1tPFEcalProducerFromOfflineRechits:towers','l1tPFHGCalEEProducerFromOfflineRechits:towers', 'l1tPFHFProducerFromOfflineRechits:towers'), + Hcal = cms.VInputTag('l1tPFHcalProducerFromOfflineRechits:towers','l1tPFHGCalFHProducerFromOfflineRechits:towers', 'l1tPFHGCalBHProducerFromOfflineRechits:towers', 'l1tPFHFProducerFromOfflineRechits:towers'), + Calo = cms.VInputTag('l1tPFEcalProducerFromOfflineRechits:towers','l1tPFHGCalEEProducerFromOfflineRechits:towers', 'l1tPFHcalProducerFromOfflineRechits:towers', 'l1tPFHGCalFHProducerFromOfflineRechits:towers', 'l1tPFHGCalBHProducerFromOfflineRechits:towers', 'l1tPFHFProducerFromOfflineRechits:towers'), + TK = cms.VInputTag('l1tPFTkProducersFromOfflineTracksStrips'), + # -- processed -- + L1RawCalo = cms.VInputTag("InfoOut:RawCalo",), + L1NewCalo = cms.VInputTag("InfoOut:NewCalo",), + L1Calo = cms.VInputTag("InfoOut:Calo",), + L1TK = cms.VInputTag("InfoOut:TK",), + L1PF = cms.VInputTag("InfoOut:PF",), + L1Puppi = cms.VInputTag("InfoOut:Puppi",), + ## -- clustered -- + #L1ak4RawCalo = cms.VInputTag("ak4L1RawCalo",), + #L1ak4Calo = cms.VInputTag("ak4L1Calo",), + #L1ak4TK = cms.VInputTag("ak4L1TK",), + #L1ak4PF = cms.VInputTag("ak4L1PF",), + #L1ak4Puppi = cms.VInputTag("ak4L1Puppi",), +) +process.p = cms.Path(process.ntuple) +process.TFileService = cms.Service("TFileService", fileName = cms.string("respTupleNew.root")) +if True: + process.load('FastPUPPI.NtupleProducer.ntupleProducer_cfi') + process.InfoOut.outputName = ""; # turn off Ntuples + process.p = cms.Path(process.InfoOut + process.ntuple) + From 9b3aa15fb85878e920a01a896eeee8feec8d3560 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 24 Apr 2017 16:16:49 +0200 Subject: [PATCH 4/7] New clusterer --- NtupleProducer/plugins/NtupleProducer.cc | 339 +++++++----------- NtupleProducer/python/ntupleProducer_cfi.py | 44 +++ .../python/runNtupleProducer_cfg.py | 28 +- 3 files changed, 168 insertions(+), 243 deletions(-) create mode 100644 NtupleProducer/python/ntupleProducer_cfi.py diff --git a/NtupleProducer/plugins/NtupleProducer.cc b/NtupleProducer/plugins/NtupleProducer.cc index 3426f0bf1a977..a5a6b00016f7b 100644 --- a/NtupleProducer/plugins/NtupleProducer.cc +++ b/NtupleProducer/plugins/NtupleProducer.cc @@ -40,6 +40,7 @@ #include "FastPUPPI/NtupleProducer/interface/jetanalyzer.hh" #include "FastPUPPI/NtupleProducer/interface/L1TPFUtils.h" #include "FastPUPPI/NtupleProducer/interface/DiscretePF.h" +#include "FastPUPPI/NtupleProducer/interface/CaloClusterer.h" #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h" #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h" @@ -101,14 +102,13 @@ class NtupleProducer : public edm::EDProducer { int ieta, iphi; float et, corr_et, eta, phi; MyEcalCluster(int iIeta, int iIphi, float iEt, float iCorr_et, float iEta, float iPhi) : ieta(iIeta), iphi(iIphi), et(iEt), corr_et(iCorr_et), eta(iEta), phi(iPhi) {} + bool operator<(const MyEcalCluster &other) const { return eta < other.eta; } }; virtual void beginJob() override; virtual void produce(edm::Event&, const edm::EventSetup&) override; virtual void endJob() override; void genMatch(std::vector &iGenVars,int iType,double iEta,double iPhi,double iPt,const reco::GenParticleCollection &iGenParticles); - TLorentzVector getVector(double iPt[][73],int iEta,int iPhi,int iEta0,int iPhi0,double iNSigma=2,double iENoise=0.2); - void simpleCluster(std::vector &iClusters,double iEta,double iPhi,double iPt[][73],double iNSigma=2,double iENoise=0.2); virtual void beginRun(edm::Run const&, edm::EventSetup const&) override; void addPF(std::vector &iCandidates,std::string iLabel,edm::Event& iEvent); @@ -136,6 +136,9 @@ class NtupleProducer : public edm::EDProducer { // discretized version l1tpf_int::RegionMapper l1regions_; l1tpf_int::PFAlgo l1pfalgo_; + // new calo clusterer (float) + l1pf_calo::SingleCaloClusterer ecalClusterer_, hcalClusterer_; + l1pf_calo::SimpleCaloLinker caloLinker_; // debug flag int fDebug; @@ -156,7 +159,7 @@ class NtupleProducer : public edm::EDProducer { float ecal_subdet, ecal_ieta, ecal_iphi, ecal_curTwrEta, ecal_curTwrPhi, ecal_et, ecal_num,ecal_dr; float hcal_subdet, hcal_ieta, hcal_iphi, hcal_TwrR, hcal_et, hcal_num, hcal_ecal_et,hcal_ecal_etcorr,hcal_ecal_eta,hcal_ecal_phi,hcal_dr; float ecal_clust_et,ecal_clust_eta,ecal_clust_phi,ecal_corr_et; - float hcal_clust_et,hcal_clust_eta,hcal_clust_phi,hcal_corr_et,hcal_corr_emf; + float hcal_clust_et,hcal_clust_eta,hcal_clust_phi,hcal_clust_emf,hcal_corr_et,hcal_corr_emf; float ecal_genPt, ecal_genEta, ecal_genPhi, ecal_genId; float hcal_genPt, hcal_genEta, hcal_genPhi, hcal_genId; //virtual void endRun(edm::Run const&, edm::EventSetup const&) override; @@ -200,6 +203,9 @@ NtupleProducer::NtupleProducer(const edm::ParameterSet& iConfig): vtxRes_ (iConfig.getParameter ("vtxRes")), l1regions_ (iConfig), l1pfalgo_ (iConfig), + ecalClusterer_ (iConfig.getParameter("caloClusterer").getParameter("ecal")), + hcalClusterer_ (iConfig.getParameter("caloClusterer").getParameter("hcal")), + caloLinker_ (iConfig.getParameter("caloClusterer").getParameter("linker"), ecalClusterer_, hcalClusterer_), fDebug (iConfig.getUntrackedParameter("debug",0)), fOutputName (iConfig.getUntrackedParameter("outputName", "ntuple.root")), fOutputFile (0), @@ -226,6 +232,7 @@ NtupleProducer::NtupleProducer(const edm::ParameterSet& iConfig): produces("Calo"); produces("PF"); produces("Puppi"); + produces("NewCalo"); produces("L1TK"); produces("L1Calo"); produces("L1PF"); @@ -266,6 +273,8 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { connector_->clear(); rawconnector_->clear(); l1regions_.clear(); + ecalClusterer_.clear(); + hcalClusterer_.clear(); /// ----------------TRACK INFO------------------- edm::Handle> l1tks; iEvent.getByLabel(L1TrackTag_, l1tks); @@ -326,69 +335,58 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { } /// ----------------ECAL INFO------------------- - std::vector lEcals; - double lEEt[l1tpf::towerNEta()][73]; - double lHEt[l1tpf::towerNEta()][73]; - for(int i0 = 0; i0 < l1tpf::towerNEta(); i0++) { for(int i1 = 0; i1 classicecals; iEvent.getByToken(TokEcalTPTag_, classicecals); edm::Handle hgecals; iEvent.getByToken(TokHGEcalTPTag_, hgecals); - L1PFCollection ecals; - for (const l1tpf::Particle & it : *classicecals) ecals.push_back(it); - for (const l1tpf::Particle & it : *hgecals ) ecals.push_back(it); - for (const l1tpf::Particle & it : ecals) { - lEEt[it.aEta()][it.aPhi()] = it.pt(); - } - std::vector pClust; - int ne = 0; - - for (const l1tpf::Particle & it : ecals) { - double et = it.pt(); - ecal_subdet = 0;// it->id().subDet(); // FIXME: missing in Particle class - ecal_ieta = it.iEta(); - ecal_iphi = it.iPhi(); - ecal_et = et; - ecal_curTwrEta = it.caloEta(); - ecal_curTwrPhi = it.caloPhi(); - ecal_num = ne; - simpleCluster(pClust,it.aEta(),it.aPhi(),lEEt); - // FIXME Gio: shouldn't we continue to the next TP here if we don't find a cluster? - // Phil : Yes, except this ntuple also contains debug info so we can look at energy lost from the cluster (You are right we don't use it) - ecal_clust_et=0,ecal_clust_eta=0,ecal_clust_phi =0; - if(pClust.size() > 0) ecal_clust_et =pClust[0].Pt(); - if(pClust.size() > 0) ecal_clust_eta=pClust[0].Eta(); - if(pClust.size() > 0) ecal_clust_phi=pClust[0].Phi(); - pClust.clear(); - ecal_corr_et = ecorrector_->correct(0.,double(ecal_clust_et),ecal_ieta,ecal_iphi); - //if(ecal_clust_et > 0) std::cout << "=ECalo=> " << ecal_clust_et << "---> " << ecal_corr_et << " ---> " << ecal_et << " -- " << ecal_clust_eta << " -- " << ecal_clust_phi << std::endl; - if(ecal_clust_et > 0) { - lEcals.push_back(MyEcalCluster(it.iEta(), it.iPhi(), ecal_clust_et, ecal_corr_et, ecal_clust_eta, ecal_clust_phi)); + for (const l1tpf::Particle & it : *classicecals) ecalClusterer_.add(it); + for (const l1tpf::Particle & it : *hgecals ) ecalClusterer_.add(it); + + ecalClusterer_.run(); + ///=== FIXME calibration goes here ==== + //ecal_corr_et = ecorrector_->correct(0.,double(ecal_clust_et),ecal_ieta,ecal_iphi); + // very naive calibration here + ecalClusterer_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { return c.et; } ); + //ecalClusteter_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { + // return c.et/(1+0.20*(std::abs(c.eta)<3)-0.15*(std::abs(c.eta)<1.5)); } ); + + // write debug output tree + if (!fOutputName.empty()) { + unsigned int ne = 0; + const auto & ecraw = ecalClusterer_.raw(); + const auto & ecals = ecalClusterer_.clusters(); + for (unsigned int i = 0, ncells = ecals.size(); i < ncells; ++i) { + if (ecals[i].et == 0) continue; + ecal_num = ne++; + ecal_et = ecraw[i]; + ecal_ieta = ecals.ieta(i); + ecal_iphi = ecals.iphi(i); + ecal_curTwrEta = ecals.eta(i); + ecal_curTwrPhi = ecals.phi(i); + ecal_clust_et = ecals[i].et; + ecal_clust_eta = ecals[i].eta; + ecal_clust_phi = ecals[i].phi; + ecal_corr_et = ecals[i].et_corr; + + std::vector lGenVars; + genMatch(lGenVars,1,ecal_clust_eta,ecal_clust_phi,ecal_clust_et,genParticles); + ecal_genPt=0; ecal_genEta=0; ecal_genPhi=0; ecal_genId=0; ecal_dr = 0; + if(lGenVars.size() > 3) { + ecal_genPt = float(lGenVars[0]); + ecal_genEta = float(lGenVars[1]); + ecal_genPhi = float(lGenVars[2]); + ecal_genId = float(lGenVars[3]); + ecal_dr = reco::deltaR( ecal_genEta, ecal_genPhi, ecal_clust_eta,ecal_clust_phi ); + } + if (zeroSuppress_) { + if(ecal_genPt > 1.) fEcalInfoTree->Fill(); + } else { + if(ecal_et > 1.) fEcalInfoTree->Fill(); + } } - // fill tree - if (fOutputName.empty()) continue; - std::vector lGenVars; - genMatch(lGenVars,1,double(it.caloEta()),double(it.caloPhi()),double(et),genParticles); - ecal_genPt=0; ecal_genEta=0; ecal_genPhi=0; ecal_genId=0; ecal_dr = 0; - if(lGenVars.size() > 3) { - ecal_genPt = float(lGenVars[0]); - ecal_genEta = float(lGenVars[1]); - ecal_genPhi = float(lGenVars[2]); - ecal_genId = float(lGenVars[3]); - ecal_dr = reco::deltaR( ecal_genEta, ecal_genPhi, ecal_clust_eta,ecal_clust_phi ); - } - if(ecal_genPt > 1. && zeroSuppress_) fEcalInfoTree->Fill(); - //if(!zeroSuppress_) fEcalInfoTree->Fill(); - if(!zeroSuppress_ && ecal_et > 1.) fEcalInfoTree->Fill(); - ne++; - // if (ne > 20) break; } + // / ----------------HCAL INFO------------------- - // / Stealing some other code! - // / https://github.com/cms-sw/cmssw/blob/0397259dd747cee94b68928f17976224c037057a/L1Trigger/L1TNtuples/src/L1AnalysisCaloTP.cc#L40 - // / 72 in phi and 56 in eta = 4032 - // / 144 TP is HF: 4 (in eta) * 18 (in phi) * 2 (sides) - // / Hcal TPs edm::Handle classichcals; iEvent.getByToken(TokHcalTPTag_, classichcals); edm::Handle hghcals; @@ -397,104 +395,70 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { iEvent.getByToken(TokBHHcalTPTag_, bhhcals); edm::Handle hfhcals; iEvent.getByToken(TokHFTPTag_, hfhcals); - L1PFCollection hcals; - for (const l1tpf::Particle & it : *classichcals) hcals.push_back(it); - for (const l1tpf::Particle & it : *hfhcals ) hcals.push_back(it); - for (const l1tpf::Particle & it : *hghcals ) { - //combine duplicates for towers we use - bool lFound = false; - for (l1tpf::Particle & it2 : hcals ) { - if(it2.eta() == it.eta() && it2.phi() == it.phi() && it2.pt() == it.pt()) {std::cout << "copied!!!!!!!!!!!!!!!!!" << std::endl; break;} - if(it2.iEta() != it.iEta() || it2.iPhi() != it.iPhi()) continue; - TLorentzVector lVec = it.tp4(); - lVec += it2.tp4(); - it2.setPtEtaPhiM(lVec.Pt(),lVec.Eta(),lVec.Phi(),lVec.M()); - lFound = true; - break; + for (const l1tpf::Particle & it : *classichcals) hcalClusterer_.add(it); + for (const l1tpf::Particle & it : *hfhcals ) hcalClusterer_.add(it); + for (const l1tpf::Particle & it : *hghcals ) hcalClusterer_.add(it); + for (const l1tpf::Particle & it : *bhhcals ) hcalClusterer_.add(it); + hcalClusterer_.run(); + + // Calorimeter linking + caloLinker_.run(); + ///=== FIXME calibration goes here ==== + /// if(hcal_clust_et > 0) hcal_corr_et = corrector_ ->correct(double(hcal_clust_et),double(hcal_ecal_et),hcal_ieta,hcal_iphi); + /// if(hcal_corr_et > 0) hcal_corr_et = corrector2_->correct(double(hcal_corr_et) ,double(hcal_clust_et),double(hcal_ecal_et),hcal_ieta,hcal_iphi); + /// if(hcal_corr_et < 0) hcal_corr_et = 0; + /// if(hcal_corr_et) hcal_corr_emf = corrector_->ecalFrac(); + caloLinker_.correct( [](const l1pf_calo::CombinedCluster &c, int ieta, int iphi) -> double { return c.et; } ); + + // write debug output tree + if (!fOutputName.empty()) { + const auto & hcraw = hcalClusterer_.raw(); + const auto & clusters = caloLinker_.clusters(); + unsigned int nh = 0; + for (unsigned int i = 0, ncells = clusters.size(); i < ncells; ++i) { + if (clusters[i].et == 0) continue; + hcal_num = nh++; + hcal_et = clusters[i].hcal_et; + hcal_ecal_et = clusters[i].ecal_et; + hcal_ecal_eta = -999; // FIXME missing + hcal_ecal_phi = -999; // FIXME missing + hcal_ieta = clusters.ieta(i); + hcal_iphi = clusters.iphi(i); + hcal_clust_et = clusters[i].et; // all these are of the combined + hcal_clust_eta = clusters[i].eta; // cluster (ecal+hcal) + hcal_clust_phi = clusters[i].phi; + hcal_clust_emf = clusters[i].ecal_et / clusters[i].et; // note: this is raw EMF + hcal_corr_et = clusters[i].et_corr; + hcal_corr_emf = -999; // FIXME + std::vector lGenVars; + genMatch(lGenVars,1,hcal_clust_eta,hcal_clust_phi,hcal_clust_et,genParticles); + hcal_genPt=0; hcal_genEta=0; hcal_genPhi=0; hcal_genId=0; hcal_dr = 0; + if(lGenVars.size() > 3) { + hcal_genPt = float(lGenVars[0]); + hcal_genEta = float(lGenVars[1]); + hcal_genPhi = float(lGenVars[2]); + hcal_genId = float(lGenVars[3]); + hcal_dr = reco::deltaR( hcal_genEta, hcal_genPhi, hcal_clust_eta,hcal_clust_phi ); + } + if (zeroSuppress_) { + if(hcal_genPt > 1.) fHcalInfoTree->Fill(); + } else { + if(hcal_et > 1.) fHcalInfoTree->Fill(); + } } - if(!lFound) hcals.push_back(it); - } - for (const l1tpf::Particle & it : *bhhcals ) { - //Flatten the BH to the HG - bool lFound = false; - for (l1tpf::Particle & it2 : hcals) { - if(it2.iPhi() == it.iPhi() && it2.iEta() == it.iEta()) { - TLorentzVector lVec = it.tp4(); - lVec += it2.tp4(); - it2.setPtEtaPhiM(lVec.Pt(),lVec.Eta(),lVec.Phi(),lVec.M()); - lFound = true; - break; - } - } - if(!lFound) hcals.push_back(it); } - for (const l1tpf::Particle & it : hcals) lHEt[it.aEta()][it.aPhi()] += it.pt(); - - unsigned nh = 0; - for (const l1tpf::Particle & it : hcals) { - hcal_subdet = 0; // it->id().subdet(); // FIXME: missing in Particle class - hcal_ieta = it.iEta(); - hcal_iphi = it.iPhi(); - hcal_TwrR = 0; // towerR; // FIXME: missing in Particle class - hcal_num = nh; - hcal_et = it.pt(); - simpleCluster(pClust,it.aEta(),it.aPhi(),lHEt); - // FIXME GP: shouldn't there be a iEta < 31 as for ECAL, since lHEt is not fixed outside it? - //PH : If I recall by construction it should be limited, but let me follow up. - // FIXME Gio: shouldn't we continue to the next TP here if we don't find a cluster? - // PH : Unfortunately the code is written such that ecal is dealt with here, so a continue will ignore H/E = 0 clusters. I agree its confusing. - hcal_clust_et=0,hcal_clust_eta=0,hcal_clust_phi=0; - if(pClust.size() > 0) hcal_clust_et =pClust[0].Pt(); - if(pClust.size() > 0) hcal_clust_eta=pClust[0].Eta(); - if(pClust.size() > 0) hcal_clust_phi=pClust[0].Phi(); - pClust.clear(); - std::vector lGenVars; - genMatch(lGenVars,1,double(hcal_clust_eta),double(hcal_clust_phi),double(it.pt()),genParticles); - hcal_genPt=0; hcal_genEta=0; hcal_genPhi=0; hcal_genId=0; hcal_dr=0; - if(lGenVars.size() > 3) { - hcal_genPt = float(lGenVars[0]); - hcal_genEta = float(lGenVars[1]); - hcal_genPhi = float(lGenVars[2]); - hcal_genId = float(lGenVars[3]); - hcal_dr = reco::deltaR( hcal_genEta, hcal_genPhi, hcal_clust_eta,hcal_clust_phi ); - } - //Linking Ecal and Hcal by Delta 1 in iEta iPhi - bool isZero = true; if(it.pt() > 2. * 0.5) isZero = false; - bool link = false; if(!isZero && hcal_clust_et > 0) link = true; - hcal_ecal_et = 0; hcal_ecal_etcorr = 0; hcal_ecal_eta = 0; hcal_ecal_phi = 0; - for(unsigned int i1 = 0; i1 < lEcals.size(); i1++) { - if(!isZero && !link) continue; //Avoid double counting ecal - if( isZero && (it.iEta() != lEcals[i1].eta || it.iPhi() != lEcals[i1].phi)) continue; - if(abs(it.iEta()-lEcals[i1].ieta) > 1) continue; - int pDPhi = abs(it.iPhi()-lEcals[i1].iphi); if(pDPhi > l1tpf::towerNPhi(it.iEta())-pDPhi) pDPhi = l1tpf::towerNPhi(it.iEta())-pDPhi; - if(abs(pDPhi) > 1) continue; - if(lEcals[i1].et < 0.2) continue; - hcal_ecal_et = lEcals[i1].et; - hcal_ecal_etcorr = lEcals[i1].corr_et; - hcal_ecal_eta = lEcals[i1].eta; - hcal_ecal_phi = lEcals[i1].phi; - break; - } - hcal_corr_et = 0; - if(hcal_clust_et > 0) hcal_corr_et = corrector_ ->correct(double(hcal_clust_et),double(hcal_ecal_et),hcal_ieta,hcal_iphi); - if(hcal_corr_et > 0) hcal_corr_et = corrector2_->correct(double(hcal_corr_et) ,double(hcal_clust_et),double(hcal_ecal_et),hcal_ieta,hcal_iphi); - if(hcal_corr_et < 0) hcal_corr_et = 0; - if(hcal_corr_et) hcal_corr_emf = corrector_->ecalFrac(); - if (hcal_corr_et >= 1 || hcal_ecal_etcorr >= 1) { - l1tpf::Particle calo = connector_->makeCalo(double(hcal_corr_et) ,double(hcal_ecal_etcorr),double(hcal_clust_eta),double(hcal_clust_phi),double(hcal_ecal_eta),double(hcal_ecal_phi)); - connector_->addCalo(calo); - l1regions_.addCalo(calo); - } - if (hcal_clust_et >= 1 || hcal_ecal_et >= 1) { - rawconnector_->addCalo(connector_->makeCalo(double(hcal_clust_et),double(hcal_ecal_et), double(hcal_clust_eta),double(hcal_clust_phi),double(hcal_ecal_eta),double(hcal_ecal_phi))); - } - // fill tree - if (fOutputName.empty()) continue; - if(hcal_genPt > 1. && zeroSuppress_) fHcalInfoTree->Fill(); - if(!zeroSuppress_ && hcal_et > 1.) fHcalInfoTree->Fill(); - nh++; - if (nh > 99999) break; + + // Now feed calo into PF + std::vector RawCaloCands = caloLinker_.fetch(false); + std::vector CaloCands = caloLinker_.fetch(true); + for (const l1tpf::Particle & calo : caloLinker_.fetch()) { + connector_->addCalo(calo); + l1regions_.addCalo(calo); + } + for (const l1tpf::Particle & calo : RawCaloCands) { + rawconnector_->addCalo(calo); } + std::vector lRawCalo = rawconnector_->candidates(); std::vector lCorrCalo = connector_->candidates(); connector_->link(metRate_); @@ -628,64 +592,7 @@ void NtupleProducer::genMatch(std::vector &iGenVars,int iType,double iEt if(lVec.Pt() > 0) iGenVars.push_back(lVec.Phi()); if(lVec.Pt() > 0) iGenVars.push_back(lId); } -//3x3 clusterizer with potential to grow clusters commented out -TLorentzVector NtupleProducer::getVector(double iPt[][73],int iEta,int iPhi,int iEta0,int iPhi0,double iNSigma,double iENoise) { //iEtaC,iPhiC - double lPtTot = 0; - std::vector > lGrow; - for(int i0=-1; i0 < 2; i0++) { - for(int i1=-1; i1 < 2; i1++) { - if(i0 == 0 && i1 == 0) continue; - if(iEta+i0 < 0 || iEta+i0 > l1tpf::towerNEta()) continue; - int pPhi = iPhi+i1; - if(pPhi < 0) pPhi=l1tpf::towerNPhi(iEta+i0)+pPhi; - if(pPhi > l1tpf::towerNPhi(iEta+i0)) pPhi=pPhi-l1tpf::towerNPhi(iEta+i0); - if(( i1+i0 > -1 && i1 > -1) && iPt[iEta][iPhi] < iPt[iEta+i0][pPhi]) lPtTot += iPt[iEta+i0][pPhi]; - if(!(i1+i0 > -1 && i1 > -1) && iPt[iEta][iPhi] <= iPt[iEta+i0][pPhi]) lPtTot += iPt[iEta+i0][pPhi]; - /* - if(( i1+i0 > -1 && i1 > -1) && iPt[iEta][iPhi] < iPt[iEta+i0][iPhi+i1]) continue; - if(!(i1+i0 > -1 && i1 > -1) && iPt[iEta][iPhi] <= iPt[iEta+i0][iPhi+i1]) continue; - double lDEta0 = fabs(iEtaC-iEta); - double lDPhi0 = fabs(iPhiC-iPhi); - double lDEta1 = fabs(iEtaC-iEta-i0); - double lDPhi1 = fabs(iPhiC-iPhi-i1); - if(sqrt(lDEta0*lEta0+lDPhi0*lDPhi0) > sqrt(lDEta1*lEta1+lDPhi1*lDPhi1)) continue; - if(sqrt(lDEta1*lEta1+lDPhi1*lDPhi1) > fDistance) continue; - lGrow.push_back(std::pair(iEta+i0,iPhi+i1)); - */ - } - } - double lPt = iPt[iEta][iPhi]; - if(lPtTot > iNSigma * iENoise && !(iEta0 == iEta && iPhi0 == iPhi)) lPt = iPt[iEta0][iPhi0]/lPtTot * lPt; - if(lPtTot > iNSigma * iENoise && iEta0 == iEta && iPhi0 == iPhi ) lPt = 0; - TLorentzVector lVec; - lVec.SetPtEtaPhiM(lPt,l1tpf::towerEta(l1tpf::translateAEta(iEta,true)),l1tpf::towerPhi(l1tpf::translateAEta(iEta,true),l1tpf::translateAPhi(iPhi,true)),0); - //if(lPt > 0) for(unsigned int i0 = 0; i0 < lGrow.size(); i0++) { lVec += getVector(iPt,lGrow[i0].first,lGrow[i0].second,iEta,iPhi,iEta,iPhiC);} - return lVec; -} -//--- Simple Clustering Start with Local maxima (in 3x3) keep adding neighbors 2sigma above threshold require bottom left to be equal or greater (to avoid points with the same value) -void NtupleProducer::simpleCluster(std::vector &iClusters,double iEta,double iPhi,double iPt[][73],double iNSigma,double iENoise) { - for(int i0 = 0; i0 < l1tpf::towerNEta()-1; i0++) { - for(int i1 = 0; i1 < l1tpf::towerNPhi(iEta); i1++) { - if(i0 != iEta || i1 != iPhi) continue; - if (iPt[i0][i1] < iNSigma * iENoise) continue; - //Max requirement - TLorentzVector pVec = getVector(iPt,i0,i1,i0,i1,iNSigma,iENoise); - if(pVec.Pt() == 0) continue; - for(int i2=-1; i2 < 2; i2++) { - for(int i3=-1; i3 < 2; i3++) { - if(i2 == 0 && i3 == 0) continue; - if(i0+i2 < 0 || i0+i2 > l1tpf::towerNEta()-1) continue; - int pPhi = i1+i3; - if(pPhi < 0) pPhi=l1tpf::towerNPhi(iEta)+pPhi; - if(pPhi > l1tpf::towerNPhi(i0+i2)) pPhi=pPhi-l1tpf::towerNPhi(i0+i2); - pVec += getVector(iPt,i0+i2,pPhi,i0,i1,iNSigma,iENoise); - } - } - iClusters.push_back(pVec); - } - } -} // ------------ method called once each job just before starting event loop ------------ void NtupleProducer::beginJob() @@ -720,7 +627,6 @@ NtupleProducer::beginJob() fTrkInfoTree->Branch("genPhi", &genPhi, "genPhi/F"); fTrkInfoTree->Branch("genid", &genId, "genid/F"); - fEcalInfoTree->Branch("ecal_subdet", &ecal_subdet, "ecal_subdet/F"); fEcalInfoTree->Branch("ecal_ieta", &ecal_ieta, "ecal_ieta/F"); fEcalInfoTree->Branch("ecal_iphi", &ecal_iphi, "ecal_iphi/F"); fEcalInfoTree->Branch("ecal_curTwrEta", &ecal_curTwrEta, "ecal_curTwrEta/F"); @@ -738,21 +644,22 @@ NtupleProducer::beginJob() fEcalInfoTree->Branch("gendr", &ecal_dr, "ecal_dr/F"); - fHcalInfoTree->Branch("hcal_subdet", &hcal_subdet, "hcal_subdet/F"); + //fHcalInfoTree->Branch("hcal_subdet", &hcal_subdet, "hcal_subdet/F"); fHcalInfoTree->Branch("hcal_ieta", &hcal_ieta, "hcal_ieta/F"); fHcalInfoTree->Branch("hcal_iphi", &hcal_iphi, "hcal_iphi/F"); - fHcalInfoTree->Branch("hcal_TwrR", &hcal_TwrR, "hcal_TwrR/F"); + //fHcalInfoTree->Branch("hcal_TwrR", &hcal_TwrR, "hcal_TwrR/F"); fHcalInfoTree->Branch("hcal_num", &hcal_num, "hcal_num/F"); fHcalInfoTree->Branch("hcal_et", &hcal_et, "hcal_et/F"); fHcalInfoTree->Branch("hcal_clust_et", &hcal_clust_et , "hcal_clust_et/F"); fHcalInfoTree->Branch("hcal_clust_eta", &hcal_clust_eta, "hcal_clust_eta/F"); fHcalInfoTree->Branch("hcal_clust_phi", &hcal_clust_phi, "hcal_clust_phi/F"); + fHcalInfoTree->Branch("hcal_clust_emf" , &hcal_clust_emf , "hcal_clust_emf/F"); fHcalInfoTree->Branch("hcal_corr_et" , &hcal_corr_et , "hcal_corr_et/F"); fHcalInfoTree->Branch("hcal_corr_emf" , &hcal_corr_emf , "hcal_corr_emf/F"); fHcalInfoTree->Branch("ecal_et", &hcal_ecal_et, "hcal_ecal_et/F"); fHcalInfoTree->Branch("ecal_etcorr", &hcal_ecal_etcorr, "hcal_ecal_etcorr/F"); - fHcalInfoTree->Branch("ecal_eta", &hcal_ecal_eta, "hcal_ecal_eta/F"); - fHcalInfoTree->Branch("ecal_phi", &hcal_ecal_phi, "hcal_ecal_phi/F"); + //fHcalInfoTree->Branch("ecal_eta", &hcal_ecal_eta, "hcal_ecal_eta/F"); // FIXME + //fHcalInfoTree->Branch("ecal_phi", &hcal_ecal_phi, "hcal_ecal_phi/F"); // FIXME fHcalInfoTree->Branch("genPt", &hcal_genPt ,"hcal_genPt/F"); fHcalInfoTree->Branch("genEta", &hcal_genEta,"hcal_genEta/F"); fHcalInfoTree->Branch("genPhi", &hcal_genPhi,"hcal_genPhi/F"); diff --git a/NtupleProducer/python/ntupleProducer_cfi.py b/NtupleProducer/python/ntupleProducer_cfi.py new file mode 100644 index 0000000000000..917ce9bc9a0ef --- /dev/null +++ b/NtupleProducer/python/ntupleProducer_cfi.py @@ -0,0 +1,44 @@ +import FWCore.ParameterSet.Config as cms + +InfoOut = cms.EDProducer('NtupleProducer', + L1TrackTag = cms.InputTag('l1tPFTkProducersFromOfflineTracksStrips'), + EcalTPTag = cms.InputTag('l1tPFEcalProducerFromOfflineRechits','towers'), + HGEcalTPTag = cms.InputTag('l1tPFHGCalEEProducerFromOfflineRechits','towers'), + HcalTPTag = cms.InputTag('l1tPFHcalProducerFromOfflineRechits','towers'), + HGHcalTPTag = cms.InputTag('l1tPFHGCalFHProducerFromOfflineRechits','towers'), + BHHcalTPTag = cms.InputTag('l1tPFHGCalBHProducerFromOfflineRechits','towers'), + HFTPTag = cms.InputTag('l1tPFHFProducerFromOfflineRechits','towers'), + MuonTPTag = cms.InputTag('gmtStage2Digis','Muon'), + genParTag = cms.InputTag('genParticles'), + zeroSuppress = cms.bool(False), + corrector = cms.string("FastPUPPI/NtupleProducer/data/pion_eta_phi.root"), + corrector2 = cms.string("FastPUPPI/NtupleProducer/data/pion_eta_phi_res_old.root"), + ecorrector = cms.string("FastPUPPI/NtupleProducer/data/ecorr.root"), + trackres = cms.string("FastPUPPI/NtupleProducer/data/tkres.root"), + eleres = cms.string("FastPUPPI/NtupleProducer/data/eres.root"), + pionres = cms.string("FastPUPPI/NtupleProducer/data/pionres.root"), + trkPtCut = cms.double(2.0), + metRate = cms.bool(False), + etaCharged = cms.double(2.5), + puppiPtCut = cms.double(4.0), + vtxRes = cms.double(0.333), + caloClusterer = cms.PSet( + ecal = cms.PSet( + zsEt = cms.double(0.4), + seedEt = cms.double(0.5), + minClusterEt = cms.double(0.5), + ), + hcal = cms.PSet( + zsEt = cms.double(0.4), + seedEt = cms.double(0.5), + minClusterEt = cms.double(0.8), + ), + linker = cms.PSet( + hoeCut = cms.double(0.1), + minPhotonEt = cms.double(1.0), + minHadronEt = cms.double(1.0), + ), + ), + outputName = cms.untracked.string("ntuple.root"), + debug = cms.untracked.int32(0), +) diff --git a/NtupleProducer/python/runNtupleProducer_cfg.py b/NtupleProducer/python/runNtupleProducer_cfg.py index c82f0b4d5fa8d..3c8d26d8f180a 100644 --- a/NtupleProducer/python/runNtupleProducer_cfg.py +++ b/NtupleProducer/python/runNtupleProducer_cfg.py @@ -33,33 +33,7 @@ process.load('FastPUPPI.NtupleProducer.l1tPFCaloProducersFromOfflineRechits_cff') process.load('FastPUPPI.NtupleProducer.l1tPFTkProducersFromOfflineTracks_cfi') - -process.InfoOut = cms.EDProducer('NtupleProducer', - L1TrackTag = cms.InputTag('l1tPFTkProducersFromOfflineTracksStrips'), - EcalTPTag = cms.InputTag('l1tPFEcalProducerFromOfflineRechits','towers'), - HGEcalTPTag = cms.InputTag('l1tPFHGCalEEProducerFromOfflineRechits','towers'), - HcalTPTag = cms.InputTag('l1tPFHcalProducerFromOfflineRechits','towers'), - HGHcalTPTag = cms.InputTag('l1tPFHGCalFHProducerFromOfflineRechits','towers'), - BHHcalTPTag = cms.InputTag('l1tPFHGCalBHProducerFromOfflineRechits','towers'), - HFTPTag = cms.InputTag('l1tPFHFProducerFromOfflineRechits','towers'), - MuonTPTag = cms.InputTag('gmtStage2Digis','Muon'), - genParTag = cms.InputTag('genParticles'), - zeroSuppress = cms.bool(False), - corrector = cms.string("FastPUPPI/NtupleProducer/data/pion_eta_phi.root"), - corrector2 = cms.string("FastPUPPI/NtupleProducer/data/pion_eta_phi_res_old.root"), - ecorrector = cms.string("FastPUPPI/NtupleProducer/data/ecorr.root"), - trackres = cms.string("FastPUPPI/NtupleProducer/data/tkres.root"), - eleres = cms.string("FastPUPPI/NtupleProducer/data/eres.root"), - pionres = cms.string("FastPUPPI/NtupleProducer/data/pionres.root"), - trkPtCut = cms.double(2.0), - metRate = cms.bool(False), - etaCharged = cms.double(2.5), - puppiPtCut = cms.double(4.0), - vtxRes = cms.double(0.333), - debug = cms.untracked.int32(1), - outputName = cms.untracked.string("ntuple.root"), - ) - +process.load('FastPUPPI.NtupleProducer.ntupleProducer_cfi') process.l1Puppi = cms.Sequence(process.l1tPFCaloProducersFromOfflineRechits+process.l1tPFTkProducersFromOfflineTracksStrips) From e5c7fbb1d118099aabe212cf3b6cc2016d94e583 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Tue, 25 Apr 2017 00:45:55 +0200 Subject: [PATCH 5/7] Working version --- NtupleProducer/interface/CaloClusterer.h | 1 + NtupleProducer/plugins/NtupleProducer.cc | 51 ++++++++++++++---- NtupleProducer/plugins/ResponseNTuplizer.cc | 18 ++++++- NtupleProducer/python/ntupleProducer_cfi.py | 1 + NtupleProducer/python/pfdiff.py | 54 +++++++++++++------ .../python/runNtupleProducer_cfg.py | 2 + NtupleProducer/python/runRespNTupler.py | 25 ++++++++- NtupleProducer/src/CaloClusterer.cc | 16 +++--- 8 files changed, 130 insertions(+), 38 deletions(-) diff --git a/NtupleProducer/interface/CaloClusterer.h b/NtupleProducer/interface/CaloClusterer.h index dd474958d34f2..5c2f5de0d5dd1 100644 --- a/NtupleProducer/interface/CaloClusterer.h +++ b/NtupleProducer/interface/CaloClusterer.h @@ -177,6 +177,7 @@ namespace l1pf_calo { PreClusterGrid ecalToHCal_; CombinedClusterGrid cluster_; float hoeCut_, minPhotonEt_, minHadronEt_; + bool useCorrectedEcal_; }; } // namespace diff --git a/NtupleProducer/plugins/NtupleProducer.cc b/NtupleProducer/plugins/NtupleProducer.cc index a5a6b00016f7b..eebc8ea7969f6 100644 --- a/NtupleProducer/plugins/NtupleProducer.cc +++ b/NtupleProducer/plugins/NtupleProducer.cc @@ -232,7 +232,6 @@ NtupleProducer::NtupleProducer(const edm::ParameterSet& iConfig): produces("Calo"); produces("PF"); produces("Puppi"); - produces("NewCalo"); produces("L1TK"); produces("L1Calo"); produces("L1PF"); @@ -343,12 +342,23 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { for (const l1tpf::Particle & it : *hgecals ) ecalClusterer_.add(it); ecalClusterer_.run(); - ///=== FIXME calibration goes here ==== - //ecal_corr_et = ecorrector_->correct(0.,double(ecal_clust_et),ecal_ieta,ecal_iphi); - // very naive calibration here - ecalClusterer_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { return c.et; } ); - //ecalClusteter_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { - // return c.et/(1+0.20*(std::abs(c.eta)<3)-0.15*(std::abs(c.eta)<1.5)); } ); + //// ---- Dummy calibration == no calibration + // ecalClusterer_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { return c.et; } ); + // + //// ---- Trivial calibration by hand + ecalClusterer_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { + if (std::abs(c.eta)<1.5) { + return c.et - (3.0 - std::abs(c.eta)); // it looks like otherwise there's an offset + } else if (std::abs(c.eta)<3) { + return c.et/1.2; // HGCal scale is off by ~1.2% + } else { + return c.et; + } + } ); + //// ---- Calibration from Phil's workflow + //ecalClusterer_.correct( [](const l1pf_calo::Cluster &c, int ieta, int iphi) -> double { + // return ecorrector_->correct(0., c.et, ieta, iphi); // FIXME is ieta or abs(ieta) ? + //} ); // write debug output tree if (!fOutputName.empty()) { @@ -408,11 +418,22 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { /// if(hcal_corr_et > 0) hcal_corr_et = corrector2_->correct(double(hcal_corr_et) ,double(hcal_clust_et),double(hcal_ecal_et),hcal_ieta,hcal_iphi); /// if(hcal_corr_et < 0) hcal_corr_et = 0; /// if(hcal_corr_et) hcal_corr_emf = corrector_->ecalFrac(); - caloLinker_.correct( [](const l1pf_calo::CombinedCluster &c, int ieta, int iphi) -> double { return c.et; } ); + // + //// ---- Trivial calibration by hand + // + caloLinker_.correct( [](const l1pf_calo::CombinedCluster &c, int ieta, int iphi) -> double { + if (std::abs(c.eta)<3.0) { + return c.ecal_et + c.hcal_et * 1.25; + } else { + return c.et; + } + } ); + // + //// ---- Dummy calibration (no calibration at all) + // caloLinker_.correct( [](const l1pf_calo::CombinedCluster &c, int ieta, int iphi) -> double { return c.et; } ); // write debug output tree if (!fOutputName.empty()) { - const auto & hcraw = hcalClusterer_.raw(); const auto & clusters = caloLinker_.clusters(); unsigned int nh = 0; for (unsigned int i = 0, ncells = clusters.size(); i < ncells; ++i) { @@ -448,10 +469,18 @@ NtupleProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { } } - // Now feed calo into PF + // Get particles from the clusterer std::vector RawCaloCands = caloLinker_.fetch(false); std::vector CaloCands = caloLinker_.fetch(true); - for (const l1tpf::Particle & calo : caloLinker_.fetch()) { + // FIXME the sigma is known to the combiner, not the calo clusterer, at the moment + for (l1tpf::Particle & calo : CaloCands) { + calo.setSigma(calo.pdgId() == combiner::GAMMA ? + connector_->getEleRes(calo.pt(), calo.eta(), calo.phi()) : + connector_->getPionRes(calo.pt(), calo.eta(), calo.phi())); + } + + // pass to the PF algo + for (const l1tpf::Particle & calo : CaloCands) { connector_->addCalo(calo); l1regions_.addCalo(calo); } diff --git a/NtupleProducer/plugins/ResponseNTuplizer.cc b/NtupleProducer/plugins/ResponseNTuplizer.cc index 254591f950395..b43adee07d72b 100644 --- a/NtupleProducer/plugins/ResponseNTuplizer.cc +++ b/NtupleProducer/plugins/ResponseNTuplizer.cc @@ -40,6 +40,7 @@ //#include "CommonTools/Utils/interface/StringCutObjectSelector.h" //#include "CommonTools/Utils/interface/StringObjectFunction.h" +#include #include @@ -131,7 +132,9 @@ class ResponseNTuplizer : public edm::one::EDAnalyzer edm::EDGetTokenT> genjets_; edm::EDGetTokenT> genparticles_; + bool isParticleGun_; TTree *tree_; + uint32_t run_, lumi_; uint64_t event_; struct McVars { float pt, pt02, eta, phi, iso02, iso04; int id; @@ -169,11 +172,16 @@ class ResponseNTuplizer : public edm::one::EDAnalyzer ResponseNTuplizer::ResponseNTuplizer(const edm::ParameterSet& iConfig) : genjets_(consumes>(iConfig.getParameter("genJets"))), - genparticles_(consumes>(iConfig.getParameter("genParticles"))) + genparticles_(consumes>(iConfig.getParameter("genParticles"))), + isParticleGun_(iConfig.getParameter("isParticleGun")), + random_(new TRandom3()) { usesResource("TFileService"); edm::Service fs; tree_ = fs->make("tree","tree"); + tree_->Branch("run", &run_, "run/i"); + tree_->Branch("lumi", &lumi_, "lumi/i"); + tree_->Branch("event", &event_, "event/l"); auto reconames = iConfig.getParameterNamesForType>(); for (const std::string & name : reconames) { @@ -198,6 +206,10 @@ ResponseNTuplizer::beginJob() void ResponseNTuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + run_ = iEvent.id().run(); + lumi_ = iEvent.id().luminosityBlock(); + event_ = iEvent.id().event(); + edm::Handle> genjets; edm::Handle> genparticles; iEvent.getByToken(genjets_, genjets); @@ -205,6 +217,10 @@ ResponseNTuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSet std::vector prompts, taus; for (const reco::GenParticle &gen : *genparticles) { + if (isParticleGun_) { + if (gen.statusFlags().isPrompt() == 1) prompts.push_back(&gen); + continue; + } if ((gen.isPromptFinalState() || gen.isDirectPromptTauDecayProductFinalState()) && (std::abs(gen.pdgId()) == 11 || std::abs(gen.pdgId()) == 13) && gen.pt() > 5) { prompts.push_back(&gen); } else if (gen.isPromptFinalState() && std::abs(gen.pdgId()) == 22 && gen.pt() > 10) { diff --git a/NtupleProducer/python/ntupleProducer_cfi.py b/NtupleProducer/python/ntupleProducer_cfi.py index 917ce9bc9a0ef..b3897b213504d 100644 --- a/NtupleProducer/python/ntupleProducer_cfi.py +++ b/NtupleProducer/python/ntupleProducer_cfi.py @@ -37,6 +37,7 @@ hoeCut = cms.double(0.1), minPhotonEt = cms.double(1.0), minHadronEt = cms.double(1.0), + useCorrectedEcal = cms.bool(True), # use corrected ecal enery in linking ), ), outputName = cms.untracked.string("ntuple.root"), diff --git a/NtupleProducer/python/pfdiff.py b/NtupleProducer/python/pfdiff.py index 4cb7fc9979960..acb1561c7c0bf 100644 --- a/NtupleProducer/python/pfdiff.py +++ b/NtupleProducer/python/pfdiff.py @@ -8,14 +8,20 @@ import os from sys import argv + from math import * from DataFormats.FWLite import Handle, Events -from PhysicsTools.HeppyCore.statistics.tree import Tree from PhysicsTools.HeppyCore.utils.deltar import * -tree = Tree("t","t",defaultFloatType="F") -events = Events(argv[1]) +from optparse import OptionParser +parser = OptionParser("%(prog) infile [ src [ dst ] ]") +parser.add_option("-n", "--events", type=int, nargs=1, default=3, help="Number of events to consider") +parser.add_option("-E", "--select-events", type=str, nargs=1, action="append", default=[], help="Specific event to select (run:lumi:event)") +options, args = parser.parse_args() +while len(args) < 3: args.append("") + +events = Events(args[0]) pf1 = Handle("std::vector") pf2 = Handle("std::vector") @@ -23,18 +29,19 @@ genp = Handle("std::vector") for iev,event in enumerate(events): - print "\nEvent %1d %5d %12d" % ( - event.eventAuxiliary().run(), - event.eventAuxiliary().luminosityBlock(), - event.eventAuxiliary().event()) + if iev >= options.events: break + idev = "%d:%d:%d" % ( event.eventAuxiliary().run(), event.eventAuxiliary().luminosityBlock(), event.eventAuxiliary().event()) + if options.select_events: + if idev not in options.select_events: continue + print "Event %s" % idev - if len(argv) == 2: + if not args[1] or ("gen:" in args[2]): event.getByLabel("ak4GenJetsNoNu", genj) allGenJ = [ g for g in genj.product() ] event.getByLabel("genParticles", genp) - genLeptons = [ p for p in genp.product() if p.status() == 1 and abs(p.pdgId()) in (11,13) and (p.isPromptFinalState() or p.isDirectPromptTauDecayProductFinalState()) and p.pt() > 5 ] + genLeptons = [ p for p in genp.product() if p.status() == 1 and abs(p.pdgId()) in (11,13,211) and (p.isPromptFinalState() or p.isDirectPromptTauDecayProductFinalState()) and p.pt() > 5 ] genTau = [ p for p in genp.product() if abs(p.pdgId()) == 15 and p.isLastCopy() and p.isPromptDecayed() and p.pt() > 5 ] - genGamma = [ p for p in genp.product() if p.status() == 1 and abs(p.pdgId()) == 22 and p.isPromptFinalState() and p.pt() > 10 ] + genGamma = [ p for p in genp.product() if p.status() == 1 and abs(p.pdgId()) == 22 and p.isPromptFinalState() and p.pt() > 5 ] genItems = genLeptons + genTau + genGamma for j in allGenJ: incone = [ i for i in genItems if deltaR(i,j) < 0.4 ] @@ -45,18 +52,33 @@ else: j.mcId = 99 allGenJ.sort(key = lambda j : - j.pt()) + if not args[1]: print "Gen jets in the event:" for j in allGenJ: if abs(j.eta()) > 5: continue if j.mcId in (0, 1, 99) and j.pt() < 5: continue print " pt %7.2f eta %+5.2f phi %+5.2f id %d" % (j.pt(), j.eta(), j.phi(), j.mcId) print "" + elif not args[2]: + event.getByLabel(args[1], pf1) + o1 = [ o for o in sorted(pf1.product(), key = lambda p : -p.pt()) ] + for p1 in o1: + print "pdgId %+4d pt %7.2f eta %+5.2f phi %+5.2f \n" % (p1.pdgId(), p1.pt(), p1.eta(), p1.phi()), + print "" else: - print "Comparing %s to %s" % (argv[2], argv[3]) - event.getByLabel(argv[2], pf1) - event.getByLabel(argv[3], pf2) + print "Comparing %s to %s" % (args[1], args[2]) + event.getByLabel(args[1], pf1) o1 = [ o for o in sorted(pf1.product(), key = lambda p : -p.pt()) ] - o2 = [ o for o in sorted(pf2.product(), key = lambda p : -p.pt()) ] + if "gen:" in args[2]: + (g,gwhat) = args[2].split(":") + if gwhat == "e" : o2 = [ g for g in genLeptons if abs(g.pdgId()) == 11 ] + elif gwhat == "mu": o2 = [ g for g in genLeptons if abs(g.pdgId()) == 13 ] + elif gwhat == "pi": o2 = [ g for g in genLeptons if abs(g.pdgId()) == 211 ] + elif gwhat == "gamma": o2 = genGamma + else: raise RuntimeError, "Not yet supported" + else: + event.getByLabel(args[2], pf2) + o2 = [ o for o in sorted(pf2.product(), key = lambda p : -p.pt()) ] match = matchObjectCollection3(o1, o2, deltaRMax=0.3) for p2 in o2: p2.seen = False for p1 in o1: @@ -76,6 +98,4 @@ for p2 in o2: if p2.seen: continue print " ... -> pdgId %+4d pt %7.2f eta %+5.2f phi %+5.2f " % (p2.pdgId(), p2.pt(), p2.eta(), p2.phi()) - if len(argv) > 4 and iev < int(argv[4]): - print "\n\n"; continue - exit() + print "" diff --git a/NtupleProducer/python/runNtupleProducer_cfg.py b/NtupleProducer/python/runNtupleProducer_cfg.py index 3c8d26d8f180a..505c408484f24 100644 --- a/NtupleProducer/python/runNtupleProducer_cfg.py +++ b/NtupleProducer/python/runNtupleProducer_cfg.py @@ -30,6 +30,7 @@ '/store/relval/CMSSW_9_1_0_pre1/RelValZMM_14/GEN-SIM-RECO/90X_upgrade2023_realistic_v9_D12-v1/00000/0A8B6505-A215-E711-B25B-0CC47A4C8F26.root' ) ) +process.source.duplicateCheckMode = cms.untracked.string("noDuplicateCheck") process.load('FastPUPPI.NtupleProducer.l1tPFCaloProducersFromOfflineRechits_cff') process.load('FastPUPPI.NtupleProducer.l1tPFTkProducersFromOfflineTracks_cfi') @@ -41,6 +42,7 @@ if False: # turn on CMSSW downstream processing and output process.InfoOut.outputName = ""; # turn off Ntuples + process.InfoOut.debug = 0; from RecoJets.JetProducers.ak4PFJets_cfi import ak4PFJets process.ak4L1RawCalo = ak4PFJets.clone(src = 'InfoOut:RawCalo') diff --git a/NtupleProducer/python/runRespNTupler.py b/NtupleProducer/python/runRespNTupler.py index 953f77a2fe981..10515105de404 100644 --- a/NtupleProducer/python/runRespNTupler.py +++ b/NtupleProducer/python/runRespNTupler.py @@ -11,10 +11,12 @@ process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring('file:l1pf_out.root') ) +process.source.duplicateCheckMode = cms.untracked.string("noDuplicateCheck") process.ntuple = cms.EDAnalyzer("ResponseNTuplizer", genJets = cms.InputTag("ak4GenJetsNoNu"), genParticles = cms.InputTag("genParticles"), + isParticleGun = cms.bool(False), # -- inputs -- Ecal = cms.VInputTag('l1tPFEcalProducerFromOfflineRechits:towers','l1tPFHGCalEEProducerFromOfflineRechits:towers', 'l1tPFHFProducerFromOfflineRechits:towers'), Hcal = cms.VInputTag('l1tPFHcalProducerFromOfflineRechits:towers','l1tPFHGCalFHProducerFromOfflineRechits:towers', 'l1tPFHGCalBHProducerFromOfflineRechits:towers', 'l1tPFHFProducerFromOfflineRechits:towers'), @@ -22,7 +24,6 @@ TK = cms.VInputTag('l1tPFTkProducersFromOfflineTracksStrips'), # -- processed -- L1RawCalo = cms.VInputTag("InfoOut:RawCalo",), - L1NewCalo = cms.VInputTag("InfoOut:NewCalo",), L1Calo = cms.VInputTag("InfoOut:Calo",), L1TK = cms.VInputTag("InfoOut:TK",), L1PF = cms.VInputTag("InfoOut:PF",), @@ -40,4 +41,24 @@ process.load('FastPUPPI.NtupleProducer.ntupleProducer_cfi') process.InfoOut.outputName = ""; # turn off Ntuples process.p = cms.Path(process.InfoOut + process.ntuple) - +if False: + process.out = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string("l1pf_remade.root"), + outputCommands = cms.untracked.vstring("drop *", + "keep *_genParticles_*_*", + "keep *_ak4GenJetsNoNu_*_*", + "keep *_InfoOut_*_*", + ) + ) + process.e = cms.EndPath(process.out) +if False: + process.MessageLogger.cerr.FwkReport.reportEvery = 1 + process.maxEvents.input = 3 + process.InfoOut.debug = 2 + if True: + process.filter = cms.EDFilter("CandViewSelector", + src = cms.InputTag("genParticles"), + cut = cms.string("pt > 10 && abs(eta) < 1.5"), + filter = cms.bool(True), + ) + process.p = cms.Path(process.filter + process.InfoOut) diff --git a/NtupleProducer/src/CaloClusterer.cc b/NtupleProducer/src/CaloClusterer.cc index 99f443571a5c4..935adfc561051 100644 --- a/NtupleProducer/src/CaloClusterer.cc +++ b/NtupleProducer/src/CaloClusterer.cc @@ -234,7 +234,8 @@ l1pf_calo::SimpleCaloLinker::SimpleCaloLinker(const edm::ParameterSet &pset, con cluster_(*grid_), hoeCut_(pset.getParameter("hoeCut")), minPhotonEt_(pset.getParameter("minPhotonEt")), - minHadronEt_(pset.getParameter("minHadronEt")) + minHadronEt_(pset.getParameter("minHadronEt")), + useCorrectedEcal_(pset.getParameter("useCorrectedEcal")) { assert(grid_->size() == ecal.raw().grid().size()); assert(grid_->size() == hcal.raw().grid().size()); @@ -263,7 +264,7 @@ void l1pf_calo::SimpleCaloLinker::run() for (int ineigh = 0; ineigh < 8; ++ineigh) { tot += hcals.neigh(i,ineigh).et; } - ecalToHCal_[i].ptOverNeighLocalMaxSum = tot ? ecals[i].et/tot : 0; + ecalToHCal_[i].ptOverNeighLocalMaxSum = tot ? (useCorrectedEcal_ ? ecals[i].et_corr : ecals[i].et)/tot : 0; } } } @@ -274,10 +275,10 @@ void l1pf_calo::SimpleCaloLinker::run() if (hcals[i].et > 0) { if (ecalToHCal_[i].ptLocalMax > 0) { // direct linking is easy - cluster_[i].ecal_et = ecals[i].et; + cluster_[i].ecal_et = (useCorrectedEcal_ ? ecals[i].et_corr : ecals[i].et); cluster_[i].hcal_et = hcals[i].et; - cluster_[i].et = hcals[i].et + ecals[i].et; - float wecal = ecals[i].et/cluster_[i].et, whcal = 1.0 - wecal; + cluster_[i].et = cluster_[i].ecal_et + cluster_[i].hcal_et; + float wecal = cluster_[i].ecal_et/cluster_[i].et, whcal = 1.0 - wecal; cluster_[i].eta = ecals[i].eta * wecal + hcals[i].eta * whcal; cluster_[i].phi = ecals[i].phi * wecal + hcals[i].phi * whcal; // wrap around phi @@ -313,9 +314,9 @@ void l1pf_calo::SimpleCaloLinker::run() for (i = 0; i < ncells; ++i) { if (ecals[i].et > 0 && ecalToHCal_[i].ptLocalMax == 0 && ecalToHCal_[i].ptOverNeighLocalMaxSum == 0) { // direct linking is easy - cluster_[i].ecal_et = ecals[i].et; + cluster_[i].ecal_et = (useCorrectedEcal_ ? ecals[i].et_corr : ecals[i].et); cluster_[i].hcal_et = hraw[i]; - cluster_[i].et = hraw[i] + ecals[i].et; + cluster_[i].et = cluster_[i].ecal_et + cluster_[i].hcal_et; cluster_[i].eta = ecals[i].eta; cluster_[i].phi = ecals[i].phi; // no need to wrap around phi @@ -333,6 +334,7 @@ std::vector l1pf_calo::SimpleCaloLinker::fetch(bool corrected) bool photon = (cluster_[i].hcal_et < hoeCut_* cluster_[i].ecal_et); if (cluster_[i].et > (photon ? minPhotonEt_ : minHadronEt_)) { ret.emplace_back(corrected ? cluster_[i].et_corr : cluster_[i].et, cluster_[i].eta, cluster_[i].phi, 0.0, photon ? 3 : 2); + ret.back().setCaloEtaPhi(cluster_[i].eta, cluster_[i].phi); } } } From 35d5878503ccb00868d4f5f2140d0ef78ffeeeb6 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Tue, 25 Apr 2017 01:03:17 +0200 Subject: [PATCH 6/7] use cfi --- .../python/runNtupleProducer_cfg_tmp.py | 30 ++----------------- NtupleProducer/python/runRespNTupler.py | 4 +-- 2 files changed, 5 insertions(+), 29 deletions(-) diff --git a/NtupleProducer/python/runNtupleProducer_cfg_tmp.py b/NtupleProducer/python/runNtupleProducer_cfg_tmp.py index 8f9782a78747b..8d3d4237e123c 100644 --- a/NtupleProducer/python/runNtupleProducer_cfg_tmp.py +++ b/NtupleProducer/python/runNtupleProducer_cfg_tmp.py @@ -33,33 +33,9 @@ process.load('FastPUPPI.NtupleProducer.l1tPFCaloProducersFromOfflineRechits_cff') process.load('FastPUPPI.NtupleProducer.l1tPFTkProducersFromOfflineTracks_cfi') - -process.InfoOut = cms.EDProducer('NtupleProducer', - L1TrackTag = cms.InputTag('l1tPFTkProducersFromOfflineTracksStrips'), - EcalTPTag = cms.InputTag('l1tPFEcalProducerFromOfflineRechits','towers'), - HGEcalTPTag = cms.InputTag('l1tPFHGCalEEProducerFromOfflineRechits','towers'), - HcalTPTag = cms.InputTag('l1tPFHcalProducerFromOfflineRechits','towers'), - HGHcalTPTag = cms.InputTag('l1tPFHGCalFHProducerFromOfflineRechits','towers'), - BHHcalTPTag = cms.InputTag('l1tPFHGCalBHProducerFromOfflineRechits','towers'), - HFTPTag = cms.InputTag('l1tPFHFProducerFromOfflineRechits','towers'), - MuonTPTag = cms.InputTag('gmtStage2Digis','Muon'), - genParTag = cms.InputTag('genParticles'), - zeroSuppress = cms.bool(False), - corrector = cms.string("/afs/cern.ch/work/n/ntran/private/Correlator/analysis/go15-gpetruc/CMSSW_9_1_0_pre2/src/FastPUPPI/NtupleProducer/data/pion_eta_phi.root"), - corrector2 = cms.string("/afs/cern.ch/work/n/ntran/private/Correlator/analysis/go15-gpetruc/CMSSW_9_1_0_pre2/src/FastPUPPI/NtupleProducer/data/pion_eta_phi_res_old.root"), - ecorrector = cms.string("/afs/cern.ch/work/n/ntran/private/Correlator/analysis/go15-gpetruc/CMSSW_9_1_0_pre2/src/FastPUPPI/NtupleProducer/data/ecorr.root"), - trackres = cms.string("/afs/cern.ch/work/n/ntran/private/Correlator/analysis/go15-gpetruc/CMSSW_9_1_0_pre2/src/FastPUPPI/NtupleProducer/data/tkres.root"), - eleres = cms.string("/afs/cern.ch/work/n/ntran/private/Correlator/analysis/go15-gpetruc/CMSSW_9_1_0_pre2/src/FastPUPPI/NtupleProducer/data/eres.root"), - pionres = cms.string("/afs/cern.ch/work/n/ntran/private/Correlator/analysis/go15-gpetruc/CMSSW_9_1_0_pre2/src/FastPUPPI/NtupleProducer/data/pionres.root"), - trkPtCut = cms.double(YYYY), - metRate = cms.bool(ZZZZ), - etaCharged = cms.double(2.5), - puppiPtCut = cms.double(4.0), - vtxRes = cms.double(0.333), - debug = cms.untracked.int32(1), - outputName = cms.untracked.string("ntuple.root"), - ) - +process.load('FastPUPPI.NtupleProducer.ntupleProducer_cfi') +process.InfoOut.trkPtCut = cms.double(YYYY) +process.InfoOut.metRate = cms.bool(ZZZZ) process.l1Puppi = cms.Sequence(process.l1tPFCaloProducersFromOfflineRechits+process.l1tPFTkProducersFromOfflineTracksStrips) diff --git a/NtupleProducer/python/runRespNTupler.py b/NtupleProducer/python/runRespNTupler.py index 10515105de404..1458d266ea75f 100644 --- a/NtupleProducer/python/runRespNTupler.py +++ b/NtupleProducer/python/runRespNTupler.py @@ -55,10 +55,10 @@ process.MessageLogger.cerr.FwkReport.reportEvery = 1 process.maxEvents.input = 3 process.InfoOut.debug = 2 - if True: + if False: process.filter = cms.EDFilter("CandViewSelector", src = cms.InputTag("genParticles"), cut = cms.string("pt > 10 && abs(eta) < 1.5"), filter = cms.bool(True), ) - process.p = cms.Path(process.filter + process.InfoOut) + process.p = cms.Path(process.filter + process.InfoOut) From d29511de00c00ef21a653b37b26cdf5a302c2b22 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Tue, 25 Apr 2017 01:10:46 +0200 Subject: [PATCH 7/7] Cleanup some debug code --- NtupleProducer/src/CaloClusterer.cc | 68 +++++++---------------------- 1 file changed, 15 insertions(+), 53 deletions(-) diff --git a/NtupleProducer/src/CaloClusterer.cc b/NtupleProducer/src/CaloClusterer.cc index 935adfc561051..c8c5e2b80124c 100644 --- a/NtupleProducer/src/CaloClusterer.cc +++ b/NtupleProducer/src/CaloClusterer.cc @@ -40,52 +40,17 @@ l1pf_calo::Stage1Grid::Stage1Grid() : } } } - // consistency check 1: check that find_cell works - for (float teta = 0; teta <= 5.0; teta += 0.02) { - for (float tphi = -M_PI; tphi <= M_PI; tphi += 0.02) { - find_cell(+teta, tphi); - find_cell(-teta, tphi); - } - } - // consistency check: neighbours for cells at small abs(ieta) - for (int aie = 1; aie < ietaCoarse_-1; ++aie) { - for (int is = -1; is <= +1; is += 2) { - int ie = aie*is; - for (int iph = 1; iph <= nPhi_; ++iph) { - int icell = ifind_cell(ie,iph); - int nvalid = 0; - for (int in = 0; in < 8; ++in) { - int ineigh = neighbour(icell, in); - if (ineigh != -1) { - if ((eta(icell) == eta(ineigh)) || - (std::abs(std::abs(eta(icell)-eta(ineigh)) - 0.5*(etaWidth(icell)+etaWidth(ineigh))) < 0.01)) { - if ((phi(icell) == phi(ineigh)) || - (std::abs(std::abs(deltaPhi(phi(icell),phi(ineigh))) - 0.5*(phiWidth(icell)+phiWidth(ineigh))) < 0.01)) { - nvalid++; - } else { - printf("Error in neighbour cell for ieta %+3d iphi %2d eta %+7.4f +- %.4f phi %+7.4f +- %.4f, got neigh ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f\n", - ie, iph, eta(icell), etaWidth(icell), phi(icell), phiWidth(icell), - ieta(ineigh), iphi(ineigh), eta(ineigh), etaWidth(ineigh), phi(ineigh), phiWidth(ineigh)); - } - } else { - printf("Error in neighbour cell for ieta %+3d iphi %2d eta %+7.4f +- %.4f phi %+7.4f +- %.4f, got neigh ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f\n", - ie, iph, eta(icell), etaWidth(icell), phi(icell), phiWidth(icell), - ieta(ineigh), iphi(ineigh), eta(ineigh), etaWidth(ineigh), phi(ineigh), phiWidth(ineigh)); - } - } else { - printf("Mismatch cell ieta %+3d iphi %2d neigh %d is null\n", ie,iph,in); - } - } - } - } - } + //// consistency check 1: check that find_cell works + //for (float teta = 0; teta <= 5.0; teta += 0.02) { + // for (float tphi = -M_PI; tphi <= M_PI; tphi += 0.02) { + // find_cell(+teta, tphi); + // find_cell(-teta, tphi); + // } + //} } int l1pf_calo::Stage1Grid::find_cell(float eta, float phi) const { int ieta = (eta != 0) ? std::distance(towerEtas_, std::lower_bound(towerEtas_, towerEtas_+nEta_, std::abs(eta))) : 1; - if (!(ieta > 0 && ieta < nEta_)) { - printf("Error in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d which is not valid\n", eta, phi, ieta); - } assert(ieta > 0 && ieta < nEta_); if (ieta > nEta_) ieta = nEta_; if (eta < 0) ieta = -ieta; @@ -97,18 +62,18 @@ int l1pf_calo::Stage1Grid::find_cell(float eta, float phi) const { if (std::abs(ieta) >= ietaVeryCoarse_) iphi -= (iphi%4); else if (std::abs(ieta) >= ietaCoarse_) iphi -= (iphi%2); iphi += 1; - if (!valid_ieta_iphi(ieta,iphi)) { - printf("Error in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which is not valid\n", - eta, phi, ieta, iphi); - } + //if (!valid_ieta_iphi(ieta,iphi)) { + // printf("Error in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which is not valid\n", + // eta, phi, ieta, iphi); + //} assert(valid_ieta_iphi(ieta,iphi)); int icell = ifind_cell(ieta,iphi); assert(icell != -1); - if (std::abs(eta - eta_[icell]) > 0.501*etaWidth_[icell] || std::abs(deltaPhi(phi, phi_[icell])) > 0.501*phiWidth_[icell]) { - printf("Mismatch in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f ; deta = %+7.4f dphi = %+7.4f\n", - eta, phi, ieta, iphi, eta_[icell], etaWidth_[icell], phi_[icell], phiWidth_[icell], eta - eta_[icell], deltaPhi(phi, phi_[icell])); - } + //if (std::abs(eta - eta_[icell]) > 0.501*etaWidth_[icell] || std::abs(deltaPhi(phi, phi_[icell])) > 0.501*phiWidth_[icell]) { + // printf("Mismatch in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f ; deta = %+7.4f dphi = %+7.4f\n", + // eta, phi, ieta, iphi, eta_[icell], etaWidth_[icell], phi_[icell], phiWidth_[icell], eta - eta_[icell], deltaPhi(phi, phi_[icell])); + //} //assert(std::abs(eta - eta_[icell]) <= 0.5*etaWidth_[icell]); //assert(std::abs(deltaPhi(phi, phi_[icell])) <= 0.5*phiWidth_[icell]); return icell; @@ -131,9 +96,6 @@ int l1pf_calo::Stage1Grid::imove(int ieta, int iphi, int deta, int dphi) { }; if (!valid_ieta_iphi(ie,iph)) return -1; int icell = ifind_cell(ie,iph); - if ((ie == ieta && iph == iphi)) { - printf("Error bogus move %+3d %2d + (%+1d %+1d) = %+3d %2d \n", ieta, iphi, deta, dphi, ie, iph); - } assert(!(ie == ieta && iph == iphi)); assert(icell != -1); assert(icell != ifind_cell(ieta,iphi));