File indexing completed on 2024-04-06 12:25:20
0001 #include "RecoHI/HiJetAlgos/interface/ParametrizedSubtractor.h"
0002 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0003 #include "DataFormats/Candidate/interface/Candidate.h"
0004 #include "FWCore/ParameterSet/interface/FileInPath.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "TFile.h"
0012
0013 #include <string>
0014 #include <iostream>
0015 using namespace std;
0016
0017 void ParametrizedSubtractor::rescaleRMS(double s) {
0018 for (std::map<int, double>::iterator iter = esigma_.begin(); iter != esigma_.end(); ++iter) {
0019 iter->second = s * (iter->second);
0020 }
0021 }
0022
0023 ParametrizedSubtractor::ParametrizedSubtractor(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
0024 : PileUpSubtractor(iConfig, std::move(iC)),
0025 dropZeroTowers_(iConfig.getUntrackedParameter<bool>("dropZeroTowers", true)),
0026 cbins_(nullptr) {
0027 centTag_ = iC.consumes<reco::Centrality>(
0028 iConfig.getUntrackedParameter<edm::InputTag>("centTag", edm::InputTag("hiCentrality", "", "RECO")));
0029
0030 interpolate_ = iConfig.getParameter<bool>("interpolate");
0031 sumRecHits_ = iConfig.getParameter<bool>("sumRecHits");
0032
0033 std::string ifname = "RecoHI/HiJetAlgos/data/PU_DATA.root";
0034 TFile* inf = new TFile(edm::FileInPath(ifname).fullPath().data());
0035 fPU = (TF1*)inf->Get("fPU");
0036 fMean = (TF1*)inf->Get("fMean");
0037 fRMS = (TF1*)inf->Get("fRMS");
0038 hC = (TH1D*)inf->Get("hC");
0039
0040 for (int i = 0; i < 40; ++i) {
0041 hEta.push_back((TH1D*)inf->Get(Form("hEta_%d", i)));
0042 hEtaMean.push_back((TH1D*)inf->Get(Form("hEtaMean_%d", i)));
0043 hEtaRMS.push_back((TH1D*)inf->Get(Form("hEtaRMS_%d", i)));
0044 }
0045 }
0046
0047 void ParametrizedSubtractor::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0048 LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
0049
0050
0051
0052
0053
0054
0055 edm::Handle<reco::Centrality> cent;
0056 iEvent.getByToken(centTag_, cent);
0057
0058 centrality_ = cent->EtHFhitSum();
0059 bin_ = 40 - hC->FindBin(centrality_);
0060 if (bin_ > 39)
0061 bin_ = 39;
0062 if (bin_ < 0)
0063 bin_ = 0;
0064
0065 PileUpSubtractor::setupGeometryMap(iEvent, iSetup);
0066 for (int i = ietamin_; i < ietamax_ + 1; i++) {
0067 emean_[i] = 0.;
0068 esigma_[i] = 0.;
0069 }
0070 }
0071
0072 void ParametrizedSubtractor::calculatePedestal(vector<fastjet::PseudoJet> const& coll) { return; }
0073
0074 void ParametrizedSubtractor::subtractPedestal(vector<fastjet::PseudoJet>& coll) {
0075 if (false) {
0076 return;
0077 } else {
0078 LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
0079
0080 int it = -100;
0081 vector<fastjet::PseudoJet> newcoll;
0082
0083 for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0084 input_object != fjInputsEnd;
0085 ++input_object) {
0086 reco::CandidatePtr const& itow = (*inputs_)[input_object->user_index()];
0087
0088 it = ieta(itow);
0089 iphi(itow);
0090
0091 double Original_Et = itow->et();
0092 if (sumRecHits_) {
0093 Original_Et = getEt(itow);
0094 }
0095
0096 double etnew = Original_Et - getPU(it, true, true);
0097 float mScale = etnew / input_object->Et();
0098 if (etnew < 0.)
0099 mScale = 0.;
0100
0101 math::XYZTLorentzVectorD towP4(input_object->px() * mScale,
0102 input_object->py() * mScale,
0103 input_object->pz() * mScale,
0104 input_object->e() * mScale);
0105
0106 int index = input_object->user_index();
0107 input_object->reset(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
0108 input_object->set_user_index(index);
0109 if (etnew > 0. && dropZeroTowers_)
0110 newcoll.push_back(*input_object);
0111 }
0112 if (dropZeroTowers_)
0113 coll = newcoll;
0114 }
0115 }
0116
0117 void ParametrizedSubtractor::calculateOrphanInput(vector<fastjet::PseudoJet>& orphanInput) { orphanInput = *fjInputs_; }
0118
0119 void ParametrizedSubtractor::offsetCorrectJets() {
0120 LogDebug("PileUpSubtractor") << "The subtractor correcting jets...\n";
0121 jetOffset_.clear();
0122
0123 using namespace reco;
0124
0125 (*fjInputs_) = fjOriginalInputs_;
0126 rescaleRMS(nSigmaPU_);
0127 subtractPedestal(*fjInputs_);
0128
0129 if (false) {
0130 const fastjet::JetDefinition& def = fjClusterSeq_->jet_def();
0131 if (!doAreaFastjet_ && !doRhoFastjet_) {
0132 fastjet::ClusterSequence newseq(*fjInputs_, def);
0133 (*fjClusterSeq_) = newseq;
0134 } else {
0135 fastjet::ClusterSequenceArea newseq(*fjInputs_, def, *fjActiveArea_);
0136 (*fjClusterSeq_) = newseq;
0137 }
0138
0139 (*fjJets_) = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
0140 }
0141
0142 jetOffset_.reserve(fjJets_->size());
0143
0144 vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin(), jetsEnd = fjJets_->end();
0145 for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
0146 int ijet = pseudojetTMP - fjJets_->begin();
0147 jetOffset_[ijet] = 0;
0148
0149 std::vector<fastjet::PseudoJet> towers = sorted_by_pt(fjClusterSeq_->constituents(*pseudojetTMP));
0150
0151 double newjetet = 0.;
0152 for (vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(), towEnd = towers.end(); ito != towEnd; ++ito) {
0153 const reco::CandidatePtr& originalTower = (*inputs_)[ito->user_index()];
0154 int it = ieta(originalTower);
0155 double Original_Et = originalTower->et();
0156
0157 if (sumRecHits_) {
0158 Original_Et = getEt(originalTower);
0159 }
0160
0161 double etnew = Original_Et - getPU(it, true, true);
0162 if (etnew < 0.)
0163 etnew = 0;
0164 newjetet = newjetet + etnew;
0165 jetOffset_[ijet] += Original_Et - etnew;
0166 }
0167
0168 if (sumRecHits_) {
0169 double mScale = newjetet / pseudojetTMP->Et();
0170 int cshist = pseudojetTMP->cluster_hist_index();
0171 pseudojetTMP->reset(pseudojetTMP->px() * mScale,
0172 pseudojetTMP->py() * mScale,
0173 pseudojetTMP->pz() * mScale,
0174 pseudojetTMP->e() * mScale);
0175 pseudojetTMP->set_cluster_hist_index(cshist);
0176 }
0177 }
0178 }
0179
0180 double ParametrizedSubtractor::getEt(const reco::CandidatePtr& in) const {
0181 const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0182 const GlobalPoint& pos = geo_->getPosition(ctc->id());
0183 double energy = ctc->emEnergy() + ctc->hadEnergy();
0184
0185 if (false) {
0186 energy = 0;
0187 const std::vector<DetId>& hitids = ctc->constituents();
0188 for (unsigned int i = 0; i < hitids.size(); ++i) {
0189 }
0190 }
0191
0192 double et = energy * sin(pos.theta());
0193 return et;
0194 }
0195
0196 double ParametrizedSubtractor::getEta(const reco::CandidatePtr& in) const {
0197 const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0198 const GlobalPoint& pos = geo_->getPosition(ctc->id());
0199 double eta = pos.eta();
0200 return eta;
0201 }
0202
0203 double ParametrizedSubtractor::getMeanAtTower(const reco::CandidatePtr& in) const {
0204 int it = ieta(in);
0205 return getPU(it, true, false);
0206 }
0207
0208 double ParametrizedSubtractor::getSigmaAtTower(const reco::CandidatePtr& in) const {
0209 int it = ieta(in);
0210 return getPU(it, false, true);
0211 }
0212
0213 double ParametrizedSubtractor::getPileUpAtTower(const reco::CandidatePtr& in) const {
0214 int it = ieta(in);
0215 return getPU(it, true, true);
0216 }
0217
0218 double ParametrizedSubtractor::getPU(int ieta, bool addMean, bool addSigma) const {
0219
0220
0221
0222 double em = hEtaMean[bin_]->GetBinContent(hEtaMean[bin_]->FindBin(ieta));
0223 double cm = fMean->Eval(centrality_);
0224
0225 double er = hEtaRMS[bin_]->GetBinContent(hEtaRMS[bin_]->FindBin(ieta));
0226 double cr = fRMS->Eval(centrality_);
0227
0228 if (interpolate_) {
0229 double n = 0;
0230 int hbin = 40 - bin_;
0231 double centerweight = (centrality_ - hC->GetBinCenter(hbin));
0232 double lowerweight = (centrality_ - hC->GetBinLowEdge(hbin));
0233 double upperweight = (centrality_ - hC->GetBinLowEdge(hbin + 1));
0234
0235 em *= lowerweight * upperweight;
0236 er *= lowerweight * upperweight;
0237 n += lowerweight * upperweight;
0238
0239 if (bin_ > 0) {
0240 em += upperweight * centerweight * hEtaMean[bin_]->GetBinContent(hEtaMean[bin_ - 1]->FindBin(ieta));
0241 er += upperweight * centerweight * hEtaRMS[bin_]->GetBinContent(hEtaRMS[bin_ - 1]->FindBin(ieta));
0242 n += upperweight * centerweight;
0243 }
0244
0245 if (bin_ < 39) {
0246 em += lowerweight * centerweight * hEtaMean[bin_]->GetBinContent(hEtaMean[bin_ + 1]->FindBin(ieta));
0247 er += lowerweight * centerweight * hEtaRMS[bin_]->GetBinContent(hEtaRMS[bin_ + 1]->FindBin(ieta));
0248 n += lowerweight * centerweight;
0249 }
0250 em /= n;
0251 er /= n;
0252 }
0253
0254
0255 return addMean * em * cm + addSigma * nSigmaPU_ * er * cr;
0256 }