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..5c2f5de0d5dd1
--- /dev/null
+++ b/NtupleProducer/interface/CaloClusterer.h
@@ -0,0 +1,185 @@
+#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_;
+ bool useCorrectedEcal_;
+ };
+
+} // namespace
+
+#endif
diff --git a/NtupleProducer/plugins/BuildFile.xml b/NtupleProducer/plugins/BuildFile.xml
index 759c6413898df..7110d784e45d8 100644
--- a/NtupleProducer/plugins/BuildFile.xml
+++ b/NtupleProducer/plugins/BuildFile.xml
@@ -21,4 +21,6 @@
+
+
diff --git a/NtupleProducer/plugins/NtupleProducer.cc b/NtupleProducer/plugins/NtupleProducer.cc
index 31b36d1a01486..eb93a4a11d4d6 100644
--- a/NtupleProducer/plugins/NtupleProducer.cc
+++ b/NtupleProducer/plugins/NtupleProducer.cc
@@ -41,6 +41,7 @@
#include "FastPUPPI/NtupleProducer/interface/isoanalyzer.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"
@@ -102,14 +103,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);
@@ -138,6 +138,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;
@@ -158,7 +161,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;
@@ -202,6 +205,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),
@@ -270,6 +276,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);
@@ -330,69 +338,69 @@ 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 = 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) {
- 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));
- }
- // 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 );
+ for (const l1tpf::Particle & it : *classicecals) ecalClusterer_.add(it);
+ for (const l1tpf::Particle & it : *hgecals ) ecalClusterer_.add(it);
+
+ ecalClusterer_.run();
+ //// ---- 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()) {
+ 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();
+ }
}
- 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;
@@ -401,103 +409,89 @@ 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);
- 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;
- }
- 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;
+ 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();
+ //
+ //// ---- 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 & 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 : *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;
+
+ // Get particles from the clusterer
+ std::vector RawCaloCands = caloLinker_.fetch(false);
+ std::vector CaloCands = caloLinker_.fetch(true);
+ // 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);
+ }
+ for (const l1tpf::Particle & calo : RawCaloCands) {
+ rawconnector_->addCalo(calo);
}
+
std::vector lRawCalo = rawconnector_->candidates();
std::vector lCorrCalo = connector_->candidates();
connector_->link(metRate_);
@@ -643,64 +637,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()
@@ -735,7 +672,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");
@@ -753,21 +689,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/plugins/ResponseNTuplizer.cc b/NtupleProducer/plugins/ResponseNTuplizer.cc
new file mode 100644
index 0000000000000..b43adee07d72b
--- /dev/null
+++ b/NtupleProducer/plugins/ResponseNTuplizer.cc
@@ -0,0 +1,305 @@
+// -*- 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
+#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_;
+ bool isParticleGun_;
+ TTree *tree_;
+ uint32_t run_, lumi_; uint64_t event_;
+ 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"))),
+ 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) {
+ 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)
+{
+ run_ = iEvent.id().run();
+ lumi_ = iEvent.id().luminosityBlock();
+ event_ = iEvent.id().event();
+
+ 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 (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) {
+ 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/ntupleProducer_cfi.py b/NtupleProducer/python/ntupleProducer_cfi.py
new file mode 100644
index 0000000000000..b3897b213504d
--- /dev/null
+++ b/NtupleProducer/python/ntupleProducer_cfi.py
@@ -0,0 +1,45 @@
+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),
+ useCorrectedEcal = cms.bool(True), # use corrected ecal enery in linking
+ ),
+ ),
+ outputName = cms.untracked.string("ntuple.root"),
+ debug = cms.untracked.int32(0),
+)
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 bfd50bc6851fc..39619454af973 100644
--- a/NtupleProducer/python/runNtupleProducer_cfg.py
+++ b/NtupleProducer/python/runNtupleProducer_cfg.py
@@ -30,36 +30,11 @@
'/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')
-
-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)
@@ -67,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/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
new file mode 100644
index 0000000000000..1458d266ea75f
--- /dev/null
+++ b/NtupleProducer/python/runRespNTupler.py
@@ -0,0 +1,64 @@
+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.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'),
+ 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",),
+ 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)
+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 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)
diff --git a/NtupleProducer/src/CaloClusterer.cc b/NtupleProducer/src/CaloClusterer.cc
new file mode 100644
index 0000000000000..c8c5e2b80124c
--- /dev/null
+++ b/NtupleProducer/src/CaloClusterer.cc
@@ -0,0 +1,304 @@
+#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);
+ // }
+ //}
+}
+
+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;
+ 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);
+ 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")),
+ useCorrectedEcal_(pset.getParameter("useCorrectedEcal"))
+{
+ 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 ? (useCorrectedEcal_ ? ecals[i].et_corr : 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 = (useCorrectedEcal_ ? ecals[i].et_corr : ecals[i].et);
+ cluster_[i].hcal_et = hcals[i].et;
+ 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
+ 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 = (useCorrectedEcal_ ? ecals[i].et_corr : ecals[i].et);
+ cluster_[i].hcal_et = hraw[i];
+ 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
+ }
+ }
+
+}
+
+
+
+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);
+ ret.back().setCaloEtaPhi(cluster_[i].eta, cluster_[i].phi);
+ }
+ }
+ }
+ return ret;
+}