Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:45

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   // Initial values for emean_, emean2, esigma_, ntowers
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 }