File indexing completed on 2023-03-17 11:18:13
0001 #include "RecoHI/HiJetAlgos/interface/ReflectedIterator.h"
0002 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0003 #include "DataFormats/Candidate/interface/Candidate.h"
0004
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 using namespace std;
0007
0008 void ReflectedIterator::rescaleRMS(double s) {
0009 for (std::map<int, double>::iterator iter = esigma_.begin(); iter != esigma_.end(); ++iter) {
0010 iter->second = s * (iter->second);
0011 }
0012 }
0013
0014 void ReflectedIterator::offsetCorrectJets() {
0015 LogDebug("PileUpSubtractor") << "The subtractor correcting jets...\n";
0016 jetOffset_.clear();
0017
0018 using namespace reco;
0019
0020 (*fjInputs_) = fjOriginalInputs_;
0021 rescaleRMS(nSigmaPU_);
0022 subtractPedestal(*fjInputs_);
0023 const fastjet::JetDefinition& def = fjClusterSeq_->jet_def();
0024 if (!doAreaFastjet_ && !doRhoFastjet_) {
0025 fastjet::ClusterSequence newseq(*fjInputs_, def);
0026 (*fjClusterSeq_) = newseq;
0027 } else {
0028 fastjet::ClusterSequenceArea newseq(*fjInputs_, def, *fjActiveArea_);
0029 (*fjClusterSeq_) = newseq;
0030 }
0031
0032 (*fjJets_) = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
0033
0034 jetOffset_.reserve(fjJets_->size());
0035
0036 vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin(), jetsEnd = fjJets_->end();
0037 for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
0038 int ijet = pseudojetTMP - fjJets_->begin();
0039 jetOffset_[ijet] = 0;
0040
0041 std::vector<fastjet::PseudoJet> towers = sorted_by_pt(fjClusterSeq_->constituents(*pseudojetTMP));
0042
0043 double newjetet = 0.;
0044 for (vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(), towEnd = towers.end(); ito != towEnd; ++ito) {
0045 const reco::CandidatePtr& originalTower = (*inputs_)[ito->user_index()];
0046 int it = ieta(originalTower);
0047 double Original_Et = originalTower->et();
0048 double etnew = Original_Et - (*emean_.find(-it)).second - (*esigma_.find(-it)).second;
0049 if (etnew < 0.)
0050 etnew = 0;
0051 newjetet = newjetet + etnew;
0052 jetOffset_[ijet] += Original_Et - etnew;
0053 }
0054 }
0055 }
0056
0057 void ReflectedIterator::subtractPedestal(vector<fastjet::PseudoJet>& coll) {
0058 LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
0059
0060 int it = -100;
0061
0062 vector<fastjet::PseudoJet> newcoll;
0063
0064 for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0065 input_object != fjInputsEnd;
0066 ++input_object) {
0067 reco::CandidatePtr const& itow = (*inputs_)[input_object->user_index()];
0068
0069 it = ieta(itow);
0070 iphi(itow);
0071
0072 double Original_Et = itow->et();
0073 if (sumRecHits_) {
0074 Original_Et = getEt(itow);
0075 }
0076
0077 double etnew = Original_Et - (*(emean_.find(-it))).second - (*(esigma_.find(-it))).second;
0078 float mScale = etnew / input_object->Et();
0079 if (etnew < 0.)
0080 mScale = 0.;
0081
0082 math::XYZTLorentzVectorD towP4(input_object->px() * mScale,
0083 input_object->py() * mScale,
0084 input_object->pz() * mScale,
0085 input_object->e() * mScale);
0086
0087 int index = input_object->user_index();
0088 input_object->reset(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
0089 input_object->set_user_index(index);
0090
0091 if (etnew > 0. && dropZeroTowers_)
0092 newcoll.push_back(*input_object);
0093 }
0094
0095 if (dropZeroTowers_)
0096 coll = newcoll;
0097 }
0098
0099 void ReflectedIterator::calculatePedestal(vector<fastjet::PseudoJet> const& coll) {
0100 LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
0101
0102 map<int, double> emean2;
0103 map<int, int> ntowers;
0104
0105 int ietaold = -10000;
0106 int ieta0 = -100;
0107
0108
0109
0110 for (int i = ietamin_; i < ietamax_ + 1; i++) {
0111 emean_[i] = 0.;
0112 emean2[i] = 0.;
0113 esigma_[i] = 0.;
0114 ntowers[i] = 0;
0115 }
0116
0117 for (vector<fastjet::PseudoJet>::const_iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0118 input_object != fjInputsEnd;
0119 ++input_object) {
0120 const reco::CandidatePtr& originalTower = (*inputs_)[input_object->user_index()];
0121 ieta0 = ieta(originalTower);
0122 double Original_Et = originalTower->et();
0123 if (sumRecHits_) {
0124 Original_Et = getEt(originalTower);
0125 }
0126
0127 if (ieta0 - ietaold != 0) {
0128 emean_[ieta0] = emean_[ieta0] + Original_Et;
0129 emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
0130 ntowers[ieta0] = 1;
0131 ietaold = ieta0;
0132 } else {
0133 emean_[ieta0] = emean_[ieta0] + Original_Et;
0134 emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
0135 ntowers[ieta0]++;
0136 }
0137 }
0138
0139 for (map<int, int>::const_iterator gt = geomtowers_.begin(); gt != geomtowers_.end(); gt++) {
0140 int it = (*gt).first;
0141
0142 double e1 = (*(emean_.find(it))).second;
0143 double e2 = (*emean2.find(it)).second;
0144 int nt = (*gt).second - (*(ntowersWithJets_.find(it))).second;
0145
0146 LogDebug("PileUpSubtractor") << " ieta : " << it << " number of towers : " << nt << " e1 : " << e1 << " e2 : " << e2
0147 << "\n";
0148
0149 if (nt > 0) {
0150 emean_[it] = e1 / nt;
0151 double eee = e2 / nt - e1 * e1 / (nt * nt);
0152 if (eee < 0.)
0153 eee = 0.;
0154 esigma_[it] = nSigmaPU_ * sqrt(eee);
0155 } else {
0156 emean_[it] = 0.;
0157 esigma_[it] = 0.;
0158 }
0159 LogDebug("PileUpSubtractor") << " ieta : " << it << " Pedestals : " << emean_[it] << " " << esigma_[it] << "\n";
0160 }
0161 }
0162
0163 double ReflectedIterator::getEt(const reco::CandidatePtr& in) const {
0164 const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0165 const GlobalPoint& pos = geo_->getPosition(ctc->id());
0166 double energy = ctc->emEnergy() + ctc->hadEnergy();
0167 double et = energy * sin(pos.theta());
0168 return et;
0169 }
0170
0171 double ReflectedIterator::getEta(const reco::CandidatePtr& in) const {
0172 const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0173 const GlobalPoint& pos = geo_->getPosition(ctc->id());
0174 double eta = pos.eta();
0175 return eta;
0176 }
