HiJetAlgos/src/MultipleAlgoIterator.cc

0001 #include <memory>
0002
0003 #include "RecoHI/HiJetAlgos/interface/MultipleAlgoIterator.h"
0004 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0005 #include "DataFormats/Candidate/interface/Candidate.h"
0006 #include "DataFormats/Math/interface/deltaR.h"
0007
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 using namespace std;
0010
0011 MultipleAlgoIterator::MultipleAlgoIterator(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
0012     : PileUpSubtractor(iConfig, std::move(iC)),
0013       minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
0014       sumRecHits_(iConfig.getParameter<bool>("sumRecHits")),
0015       dropZeroTowers_(iConfig.getUntrackedParameter<bool>("dropZeroTowers", true)) {
0016   LogDebug("PileUpSubtractor") << "LIMITING THE MINIMUM TOWERS FRACTION TO : " << minimumTowersFraction_ << endl;
0017 }
0018
0019 void MultipleAlgoIterator::rescaleRMS(double s) {
0020   for (std::map<int, double>::iterator iter = esigma_.begin(); iter != esigma_.end(); ++iter) {
0021     iter->second = s * (iter->second);
0022   }
0023 }
0024
0025 void MultipleAlgoIterator::offsetCorrectJets() {
0026   LogDebug("PileUpSubtractor") << "The subtractor correcting jets...\n";
0027   jetOffset_.clear();
0028
0029   using namespace reco;
0030
0031   (*fjInputs_) = fjOriginalInputs_;
0032   rescaleRMS(nSigmaPU_);
0033   subtractPedestal(*fjInputs_);
0034   const fastjet::JetDefinition& def = *fjJetDefinition_;
0035   if (!doAreaFastjet_ && !doRhoFastjet_) {
0036     fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(*fjInputs_, def);
0037   } else {
0038     fjClusterSeq_ = ClusterSequencePtr(new fastjet::ClusterSequenceArea(*fjInputs_, def, *fjActiveArea_));
0039   }
0040
0041   (*fjJets_) = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
0042
0043   jetOffset_.reserve(fjJets_->size());
0044
0045   vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin(), jetsEnd = fjJets_->end();
0046   for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
0047     int ijet = pseudojetTMP - fjJets_->begin();
0048     jetOffset_[ijet] = 0;
0049
0050     std::vector<fastjet::PseudoJet> towers = sorted_by_pt(pseudojetTMP->constituents());
0051
0052     double newjetet = 0.;
0053     for (vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(), towEnd = towers.end(); ito != towEnd; ++ito) {
0054       const reco::CandidatePtr& originalTower = (*inputs_)[ito->user_index()];
0055       int it = ieta(originalTower);
0056       double Original_Et = originalTower->et();
0057       double etnew = Original_Et - (*emean_.find(it)).second - (*esigma_.find(it)).second;
0058       if (etnew < 0.)
0059         etnew = 0;
0060       newjetet = newjetet + etnew;
0061       jetOffset_[ijet] += Original_Et - etnew;
0062     }
0063   }
0064 }
0065
0066 void MultipleAlgoIterator::subtractPedestal(vector<fastjet::PseudoJet>& coll) {
0067   LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
0068
0069   int it = -100;
0070
0071   vector<fastjet::PseudoJet> newcoll;
0072
0073   for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0074        input_object != fjInputsEnd;
0075        ++input_object) {
0076     reco::CandidatePtr const& itow = (*inputs_)[input_object->user_index()];
0077
0078     it = ieta(itow);
0079     iphi(itow);
0080
0081     double Original_Et = itow->et();
0082     if (sumRecHits_) {
0083       Original_Et = getEt(itow);
0084     }
0085
0086     double etnew = Original_Et - (*(emean_.find(it))).second - (*(esigma_.find(it))).second;
0087     float mScale = etnew / input_object->Et();
0088     if (etnew < 0.)
0089       mScale = 0.;
0090
0091     math::XYZTLorentzVectorD towP4(input_object->px() * mScale,
0092                                    input_object->py() * mScale,
0093                                    input_object->pz() * mScale,
0094                                    input_object->e() * mScale);
0095
0096     int index = input_object->user_index();
0097     input_object->reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
0098     input_object->set_user_index(index);
0099
0100     if (etnew > 0. && dropZeroTowers_)
0101       newcoll.push_back(*input_object);
0102   }
0103
0104   if (dropZeroTowers_)
0105     coll = newcoll;
0106 }
0107
0108 void MultipleAlgoIterator::calculatePedestal(vector<fastjet::PseudoJet> const& coll) {
0109   LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
0110
0111   map<int, double> emean2;
0112   map<int, int> ntowers;
0113
0114   int ietaold = -10000;
0115   int ieta0 = -100;
0116
0117   // Initial values for emean_, emean2, esigma_, ntowers
0118
0119   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0120     emean_[i] = 0.;
0121     emean2[i] = 0.;
0122     esigma_[i] = 0.;
0123     ntowers[i] = 0;
0124   }
0125
0126   for (vector<fastjet::PseudoJet>::const_iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0127        input_object != fjInputsEnd;
0128        ++input_object) {
0129     const reco::CandidatePtr& originalTower = (*inputs_)[input_object->user_index()];
0130     ieta0 = ieta(originalTower);
0131     double Original_Et = originalTower->et();
0132     if (sumRecHits_) {
0133       Original_Et = getEt(originalTower);
0134     }
0135
0136     if (ieta0 - ietaold != 0) {
0137       emean_[ieta0] = emean_[ieta0] + Original_Et;
0138       emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
0139       ntowers[ieta0] = 1;
0140       ietaold = ieta0;
0141     } else {
0142       emean_[ieta0] = emean_[ieta0] + Original_Et;
0143       emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
0144       ntowers[ieta0]++;
0145     }
0146   }
0147
0148   for (map<int, int>::const_iterator gt = geomtowers_.begin(); gt != geomtowers_.end(); gt++) {
0149     int it = (*gt).first;
0150
0151     double e1 = (*(emean_.find(it))).second;
0152     double e2 = (*emean2.find(it)).second;
0153     int nt = (*gt).second - (*(ntowersWithJets_.find(it))).second;
0154
0155     LogDebug("PileUpSubtractor") << " ieta : " << it << " number of towers : " << nt << " e1 : " << e1 << " e2 : " << e2
0156                                  << "\n";
0157
0158     if (nt > 0) {
0159       emean_[it] = e1 / nt;
0160       double eee = e2 / nt - e1 * e1 / (nt * nt);
0161       if (eee < 0.)
0162         eee = 0.;
0163       esigma_[it] = nSigmaPU_ * sqrt(eee);
0164     } else {
0165       emean_[it] = 0.;
0166       esigma_[it] = 0.;
0167     }
0168     LogDebug("PileUpSubtractor") << " ieta : " << it << " Pedestals : " << emean_[it] << "  " << esigma_[it] << "\n";
0169   }
0170 }
0171
0172 void MultipleAlgoIterator::calculateOrphanInput(vector<fastjet::PseudoJet>& orphanInput) {
0173   LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
0174
0175   (*fjInputs_) = fjOriginalInputs_;
0176
0177   vector<int> jettowers;                   // vector of towers indexed by "user_index"
0178   vector<pair<int, int> > excludedTowers;  // vector of excluded ieta, iphi values
0179
0180   for (auto const& pseudojetTMP : *fjJets_) {
0181     if (pseudojetTMP.perp() < puPtMin_)
0182       continue;
0183
0184     // find towers within radiusPU_ of this jet
0185     for (auto const im : allgeomid_) {
0186       double dr = reco::deltaR(geo_->getPosition(im), pseudojetTMP);
0187       auto exclude = find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(im.ieta(), im.iphi()));
0188       if (dr < radiusPU_ && exclude == excludedTowers.end() &&
0189           (geomtowers_[im.ieta()] - ntowersWithJets_[im.ieta()]) > minimumTowersFraction_ * (geomtowers_[im.ieta()])) {
0190         ntowersWithJets_[im.ieta()]++;
0191
0192         excludedTowers.push_back(pair<int, int>(im.ieta(), im.iphi()));
0193       }
0194     }
0195
0196     for (auto const& it : *fjInputs_) {
0197       int index = it.user_index();
0198       int ie = ieta((*inputs_)[index]);
0199       int ip = iphi((*inputs_)[index]);
0200       auto exclude = find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(ie, ip));
0201       if (exclude != excludedTowers.end()) {
0202         jettowers.push_back(index);
0203       }
0204     }  // initial input collection
0205
0206   }  // pseudojets
0207
0208   //
0209   // Create a new collections from the towers not included in jets
0210   //
0211
0212   for (auto const& it : *fjInputs_) {
0213     int index = it.user_index();
0214     vector<int>::const_iterator itjet = find(jettowers.begin(), jettowers.end(), index);
0215     if (itjet == jettowers.end()) {
0216       const reco::CandidatePtr& originalTower = (*inputs_)[index];
0217       fastjet::PseudoJet orphan(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
0218       orphan.set_user_index(index);
0219
0220       orphanInput.push_back(orphan);
0221     }
0222   }
0223 }
0224
0225 double MultipleAlgoIterator::getEt(const reco::CandidatePtr& in) const {
0226   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0227   const GlobalPoint& pos = geo_->getPosition(ctc->id());
0228   double energy = ctc->emEnergy() + ctc->hadEnergy();
0229   double et = energy * sin(pos.theta());
0230   return et;
0231 }
0232
0233 double MultipleAlgoIterator::getEta(const reco::CandidatePtr& in) const {
0234   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0235   const GlobalPoint& pos = geo_->getPosition(ctc->id());
0236   double eta = pos.eta();
0237   return eta;
0238 }