JetAlgorithms/src/FixedGridEnergyDensity.cc

0001 #include "RecoJets/JetAlgorithms/interface/FixedGridEnergyDensity.h"
0002
0003 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include "TMath.h"
0006
0007 using namespace reco;
0008 using namespace edm;
0009 using namespace std;
0010
0011 float FixedGridEnergyDensity::fixedGridRho(EtaRegion etaRegion) {
0012   //define the phi bins
0013   vector<float> phibins;
0014   phibins.reserve(10);
0015   for (int i = 0; i < 10; i++)
0016     phibins.push_back(-TMath::Pi() + (2 * i + 1) * TMath::TwoPi() / 20.);
0017   //define the eta bins
0018   vector<float> etabins;
0019   if (etaRegion == Central) {
0020     for (int i = 0; i < 8; ++i)
0021       etabins.push_back(-2.1 + 0.6 * i);
0022   } else if (etaRegion == Forward) {
0023     for (int i = 0; i < 10; ++i) {
0024       if (i < 5)
0025         etabins.push_back(-5.1 + 0.6 * i);
0026       else
0027         etabins.push_back(2.7 + 0.6 * (i - 5));
0028     }
0029   } else if (etaRegion == All) {
0030     for (int i = 0; i < 18; ++i)
0031       etabins.push_back(-5.1 + 0.6 * i);
0032   }
0033   return fixedGridRho(etabins, phibins);
0034 }
0035
0036 float FixedGridEnergyDensity::fixedGridRho(std::vector<float>& etabins, std::vector<float>& phibins) {
0037   float etadist = etabins[1] - etabins[0];
0038   float phidist = phibins[1] - phibins[0];
0039   float etahalfdist = (etabins[1] - etabins[0]) / 2.;
0040   float phihalfdist = (phibins[1] - phibins[0]) / 2.;
0041   vector<float> sumPFNallSMDQ;
0042   sumPFNallSMDQ.reserve(etabins.size() * phibins.size());
0043   for (unsigned int ieta = 0; ieta < etabins.size(); ++ieta) {
0044     for (unsigned int iphi = 0; iphi < phibins.size(); ++iphi) {
0045       float pfniso_ieta_iphi = 0;
0046       for (PFCandidateCollection::const_iterator pf_it = pfCandidates->begin(); pf_it != pfCandidates->end(); pf_it++) {
0047         if (fabs(etabins[ieta] - pf_it->eta()) > etahalfdist)
0048           continue;
0049         if (fabs(reco::deltaPhi(phibins[iphi], pf_it->phi())) > phihalfdist)
0050           continue;
0051         pfniso_ieta_iphi += pf_it->pt();
0052       }
0053       sumPFNallSMDQ.push_back(pfniso_ieta_iphi);
0054     }
0055   }
0056   float evt_smdq = 0;
0057   sort(sumPFNallSMDQ.begin(), sumPFNallSMDQ.end());
0058   if (sumPFNallSMDQ.size() % 2)
0059     evt_smdq = sumPFNallSMDQ[(sumPFNallSMDQ.size() - 1) / 2];
0060   else
0061     evt_smdq = (sumPFNallSMDQ[sumPFNallSMDQ.size() / 2] + sumPFNallSMDQ[(sumPFNallSMDQ.size() - 2) / 2]) / 2.;
0062   return evt_smdq / (etadist * phidist);
0063 }