Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:34 

0001 
0002 #include "RecoJets/JetProducers/interface/PileUpSubtractor.h"
0003 
0004 #include "DataFormats/Common/interface/View.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 
0009 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0010 #include "DataFormats/Candidate/interface/Candidate.h"
0011 #include "DataFormats/Math/interface/deltaR.h"
0012 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0013 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014 
0015 #include <map>
0016 #include <memory>
0017 
0018 using namespace std;
0019 
0020 PileUpSubtractor::PileUpSubtractor(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC) {
0021   geo_ = nullptr;
0022   geoToken_ = iC.esConsumes();
0023   doAreaFastjet_ = iConfig.getParameter<bool>("doAreaFastjet");
0024   doRhoFastjet_ = iConfig.getParameter<bool>("doRhoFastjet");
0025   nSigmaPU_ = iConfig.getParameter<double>("nSigmaPU");
0026   radiusPU_ = iConfig.getParameter<double>("radiusPU");
0027   jetPtMin_ = iConfig.getParameter<double>("jetPtMin");
0028   puPtMin_ = iConfig.getParameter<double>("puPtMin");
0029   ghostEtaMax = iConfig.getParameter<double>("Ghost_EtaMax");
0030   activeAreaRepeats = iConfig.getParameter<int>("Active_Area_Repeats");
0031   ghostArea = iConfig.getParameter<double>("GhostArea");
0032 
0033   if (doAreaFastjet_ || doRhoFastjet_) {
0034     fjActiveArea_ = std::make_shared<fastjet::ActiveAreaSpec>(ghostEtaMax, activeAreaRepeats, ghostArea);
0035     if ((ghostEtaMax < 0) || (activeAreaRepeats < 0) || (ghostArea < 0))
0036       throw cms::Exception("doAreaFastjet or doRhoFastjet")
0037           << "Parameters ghostEtaMax, activeAreaRepeats or ghostArea for doAreaFastjet/doRhoFastjet are not defined."
0038           << std::endl;
0039   }
0040 }
0041 
0042 void PileUpSubtractor::reset(std::vector<edm::Ptr<reco::Candidate> >& input,
0043                              std::vector<fastjet::PseudoJet>& towers,
0044                              std::vector<fastjet::PseudoJet>& output) {
0045   inputs_ = &input;
0046   fjInputs_ = &towers;
0047   fjJets_ = &output;
0048   fjOriginalInputs_ = (*fjInputs_);
0049   for (unsigned int i = 0; i < fjInputs_->size(); ++i) {
0050     fjOriginalInputs_[i].set_user_index((*fjInputs_)[i].user_index());
0051   }
0052 }
0053 
0054 void PileUpSubtractor::setDefinition(JetDefPtr const& jetDef) {
0055   fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(*jetDef);
0056 }
0057 
0058 void PileUpSubtractor::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0059   LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
0060 
0061   if (geo_ == nullptr) {
0062     geo_ = &iSetup.getData(geoToken_);
0063     std::vector<DetId> alldid = geo_->getValidDetIds();
0064 
0065     int ietaold = -10000;
0066     ietamax_ = -10000;
0067     ietamin_ = 10000;
0068     for (std::vector<DetId>::const_iterator did = alldid.begin(); did != alldid.end(); did++) {
0069       if ((*did).det() == DetId::Hcal) {
0070         HcalDetId hid = HcalDetId(*did);
0071         allgeomid_.push_back(*did);
0072 
0073         if (hid.ieta() != ietaold) {
0074           ietaold = hid.ieta();
0075           geomtowers_[hid.ieta()] = 1;
0076           if (hid.ieta() > ietamax_)
0077             ietamax_ = hid.ieta();
0078           if (hid.ieta() < ietamin_)
0079             ietamin_ = hid.ieta();
0080         } else {
0081           geomtowers_[hid.ieta()]++;
0082         }
0083       }
0084     }
0085   }
0086 
0087   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0088     ntowersWithJets_[i] = 0;
0089   }
0090 }
0091 
0092 //
0093 // Calculate mean E and sigma from jet collection "coll".
0094 //
0095 void PileUpSubtractor::calculatePedestal(vector<fastjet::PseudoJet> const& coll) {
0096   LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
0097   map<int, double> emean2;
0098   map<int, int> ntowers;
0099 
0100   int ietaold = -10000;
0101   int ieta0 = -100;
0102 
0103   // Initial values for emean_, emean2, esigma_, ntowers
0104 
0105   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0106     emean_[i] = 0.;
0107     emean2[i] = 0.;
0108     esigma_[i] = 0.;
0109     ntowers[i] = 0;
0110   }
0111 
0112   for (vector<fastjet::PseudoJet>::const_iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0113        input_object != fjInputsEnd;
0114        ++input_object) {
0115     const reco::CandidatePtr& originalTower = (*inputs_)[input_object->user_index()];
0116     ieta0 = ieta(originalTower);
0117     double Original_Et = originalTower->et();
0118     if (ieta0 - ietaold != 0) {
0119       emean_[ieta0] = emean_[ieta0] + Original_Et;
0120       emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
0121       ntowers[ieta0] = 1;
0122       ietaold = ieta0;
0123     } else {
0124       emean_[ieta0] = emean_[ieta0] + Original_Et;
0125       emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
0126       ntowers[ieta0]++;
0127     }
0128   }
0129 
0130   for (map<int, int>::const_iterator gt = geomtowers_.begin(); gt != geomtowers_.end(); gt++) {
0131     int it = (*gt).first;
0132 
0133     double e1 = (*(emean_.find(it))).second;
0134     double e2 = (*emean2.find(it)).second;
0135     int nt = (*gt).second - (*(ntowersWithJets_.find(it))).second;
0136 
0137     LogDebug("PileUpSubtractor") << " ieta : " << it << " number of towers : " << nt << " e1 : " << e1 << " e2 : " << e2
0138                                  << "\n";
0139     if (nt > 0) {
0140       emean_[it] = e1 / nt;
0141       double eee = e2 / nt - e1 * e1 / (nt * nt);
0142       if (eee < 0.)
0143         eee = 0.;
0144       esigma_[it] = nSigmaPU_ * sqrt(eee);
0145     } else {
0146       emean_[it] = 0.;
0147       esigma_[it] = 0.;
0148     }
0149     LogDebug("PileUpSubtractor") << " ieta : " << it << " Pedestals : " << emean_[it] << "  " << esigma_[it] << "\n";
0150   }
0151 }
0152 
0153 //
0154 // Subtract mean and sigma from fjInputs_
0155 //
0156 void PileUpSubtractor::subtractPedestal(vector<fastjet::PseudoJet>& coll) {
0157   LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
0158 
0159   int it = -100;
0160   for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin(), fjInputsEnd = coll.end();
0161        input_object != fjInputsEnd;
0162        ++input_object) {
0163     reco::CandidatePtr const& itow = (*inputs_)[input_object->user_index()];
0164 
0165     it = ieta(itow);
0166 
0167     double etnew = itow->et() - (*(emean_.find(it))).second - (*(esigma_.find(it))).second;
0168     float mScale = etnew / input_object->Et();
0169     if (etnew < 0.)
0170       mScale = 0.;
0171 
0172     math::XYZTLorentzVectorD towP4(input_object->px() * mScale,
0173                                    input_object->py() * mScale,
0174                                    input_object->pz() * mScale,
0175                                    input_object->e() * mScale);
0176 
0177     int index = input_object->user_index();
0178     input_object->reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
0179     input_object->set_user_index(index);
0180   }
0181 }
0182 
0183 void PileUpSubtractor::calculateOrphanInput(vector<fastjet::PseudoJet>& orphanInput) {
0184   LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
0185 
0186   (*fjInputs_) = fjOriginalInputs_;
0187 
0188   vector<int> jettowers;                   // vector of towers indexed by "user_index"
0189   vector<pair<int, int> > excludedTowers;  // vector of excluded ieta, iphi values
0190 
0191   vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin(), fjJetsEnd = fjJets_->end();
0192   for (; pseudojetTMP != fjJetsEnd; ++pseudojetTMP) {
0193     if (pseudojetTMP->perp() < puPtMin_)
0194       continue;
0195 
0196     // find towers within radiusPU_ of this jet
0197     for (vector<HcalDetId>::const_iterator im = allgeomid_.begin(); im != allgeomid_.end(); im++) {
0198       double dr = reco::deltaR(geo_->getPosition((DetId)(*im)), (*pseudojetTMP));
0199       vector<pair<int, int> >::const_iterator exclude =
0200           find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(im->ieta(), im->iphi()));
0201       if (dr < radiusPU_ && exclude == excludedTowers.end()) {
0202         ntowersWithJets_[(*im).ieta()]++;
0203         excludedTowers.push_back(pair<int, int>(im->ieta(), im->iphi()));
0204       }
0205     }
0206     vector<fastjet::PseudoJet>::const_iterator it = fjInputs_->begin(), fjInputsEnd = fjInputs_->end();
0207 
0208     for (; it != fjInputsEnd; ++it) {
0209       int index = it->user_index();
0210       int ie = ieta((*inputs_)[index]);
0211       int ip = iphi((*inputs_)[index]);
0212       vector<pair<int, int> >::const_iterator exclude =
0213           find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(ie, ip));
0214       if (exclude != excludedTowers.end()) {
0215         jettowers.push_back(index);
0216       }  //dr < radiusPU_
0217     }  // initial input collection
0218   }  // pseudojets
0219 
0220   //
0221   // Create a new collections from the towers not included in jets
0222   //
0223   for (vector<fastjet::PseudoJet>::const_iterator it = fjInputs_->begin(), fjInputsEnd = fjInputs_->end();
0224        it != fjInputsEnd;
0225        ++it) {
0226     int index = it->user_index();
0227     vector<int>::const_iterator itjet = find(jettowers.begin(), jettowers.end(), index);
0228     if (itjet == jettowers.end()) {
0229       const reco::CandidatePtr& originalTower = (*inputs_)[index];
0230       fastjet::PseudoJet orphan(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
0231       orphan.set_user_index(index);
0232 
0233       orphanInput.push_back(orphan);
0234     }
0235   }
0236 }
0237 
0238 void PileUpSubtractor::offsetCorrectJets() {
0239   LogDebug("PileUpSubtractor") << "The subtractor correcting jets...\n";
0240   jetOffset_.clear();
0241   using namespace reco;
0242 
0243   //
0244   // Reestimate energy of jet (energy of jet with initial map)
0245   //
0246   jetOffset_.reserve(fjJets_->size());
0247   vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin(), jetsEnd = fjJets_->end();
0248   for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
0249     int ijet = pseudojetTMP - fjJets_->begin();
0250     jetOffset_[ijet] = 0;
0251 
0252     std::vector<fastjet::PseudoJet> towers = fastjet::sorted_by_pt(pseudojetTMP->constituents());
0253     double newjetet = 0.;
0254     for (vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(), towEnd = towers.end(); ito != towEnd; ++ito) {
0255       const reco::CandidatePtr& originalTower = (*inputs_)[ito->user_index()];
0256       int it = ieta(originalTower);
0257       double Original_Et = originalTower->et();
0258       double etnew = Original_Et - (*emean_.find(it)).second - (*esigma_.find(it)).second;
0259       if (etnew < 0.)
0260         etnew = 0;
0261       newjetet = newjetet + etnew;
0262       jetOffset_[ijet] += Original_Et - etnew;
0263     }
0264     double mScale = newjetet / pseudojetTMP->Et();
0265     LogDebug("PileUpSubtractor") << "pseudojetTMP->Et() : " << pseudojetTMP->Et() << '\n';
0266     LogDebug("PileUpSubtractor") << "newjetet : " << newjetet << '\n';
0267     LogDebug("PileUpSubtractor") << "jetOffset_[ijet] : " << jetOffset_[ijet] << '\n';
0268     LogDebug("PileUpSubtractor") << "pseudojetTMP->Et() - jetOffset_[ijet] : " << pseudojetTMP->Et() - jetOffset_[ijet]
0269                                  << '\n';
0270     LogDebug("PileUpSubtractor") << "Scale is : " << mScale << '\n';
0271     int cshist = pseudojetTMP->cluster_hist_index();
0272     pseudojetTMP->reset_momentum(pseudojetTMP->px() * mScale,
0273                                  pseudojetTMP->py() * mScale,
0274                                  pseudojetTMP->pz() * mScale,
0275                                  pseudojetTMP->e() * mScale);
0276     pseudojetTMP->set_cluster_hist_index(cshist);
0277   }
0278 }
0279 
0280 double PileUpSubtractor::getCone(double cone, double eta, double phi, double& et, double& pu) {
0281   pu = 0;
0282 
0283   for (vector<HcalDetId>::const_iterator im = allgeomid_.begin(); im != allgeomid_.end(); im++) {
0284     if (im->depth() != 1)
0285       continue;
0286     const GlobalPoint& point = geo_->getPosition((DetId)(*im));
0287     double dr = reco::deltaR(point.eta(), point.phi(), eta, phi);
0288     if (dr < cone) {
0289       pu += (*emean_.find(im->ieta())).second + (*esigma_.find(im->ieta())).second;
0290     }
0291   }
0292 
0293   return pu;
0294 }
0295 
0296 double PileUpSubtractor::getMeanAtTower(const reco::CandidatePtr& in) const {
0297   int it = ieta(in);
0298   return (*emean_.find(it)).second;
0299 }
0300 
0301 double PileUpSubtractor::getSigmaAtTower(const reco::CandidatePtr& in) const {
0302   int it = ieta(in);
0303   return (*esigma_.find(it)).second;
0304 }
0305 
0306 double PileUpSubtractor::getPileUpAtTower(const reco::CandidatePtr& in) const {
0307   int it = ieta(in);
0308   return (*emean_.find(it)).second + (*esigma_.find(it)).second;
0309 }
0310 
0311 int PileUpSubtractor::getN(const reco::CandidatePtr& in) const {
0312   int it = ieta(in);
0313 
0314   int n = (*(geomtowers_.find(it))).second - (*(ntowersWithJets_.find(it))).second;
0315   return n;
0316 }
0317 
0318 int PileUpSubtractor::getNwithJets(const reco::CandidatePtr& in) const {
0319   int it = ieta(in);
0320   int n = (*(ntowersWithJets_.find(it))).second;
0321   return n;
0322 }
0323 
0324 int PileUpSubtractor::ieta(const reco::CandidatePtr& in) const {
0325   int it = 0;
0326   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0327   if (ctc) {
0328     it = ctc->id().ieta();
0329   } else {
0330     throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type";
0331   }
0332   return it;
0333 }
0334 
0335 int PileUpSubtractor::iphi(const reco::CandidatePtr& in) const {
0336   int it = 0;
0337   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
0338   if (ctc) {
0339     it = ctc->id().iphi();
0340   } else {
0341     throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type";
0342   }
0343   return it;
0344 }
0345 
0346 #include "FWCore/PluginManager/interface/PluginFactory.h"
0347 EDM_REGISTER_PLUGINFACTORY(PileUpSubtractorFactory, "PileUpSubtractorFactory");
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!