HiJetAlgos/src/ParametrizedSubtractor.cc

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   // The function below that is commented out was deleted from
0051   // DataFormats/HeavyIonEvent/src/Centrality.cc
0052   // in June 2015. See comments associated with that commit.
0053   //   if(!cbins_) getCentralityBinsFromDB(iSetup);
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   //double e = hEta[bin_]->GetBinContent(hEta[bin_]->FindBin(ieta));
0220   //double c = fPU->Eval(centrality_);
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   //   return e*c;
0255   return addMean * em * cm + addSigma * nSigmaPU_ * er * cr;
0256 }