Go to the documentation of this file.00001 #include "RecoJets/JetAlgorithms/interface/FixedGridEnergyDensity.h"
00002
00003 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00004 #include "DataFormats/Math/interface/deltaPhi.h"
00005 #include "TMath.h"
00006
00007 using namespace reco;
00008 using namespace edm;
00009 using namespace std;
00010
00011 float FixedGridEnergyDensity::fixedGridRho(EtaRegion etaRegion){
00012
00013 vector<float> phibins;
00014 for (int i=0;i<10;i++) phibins.push_back(-TMath::Pi()+(2*i+1)*TMath::TwoPi()/20.);
00015
00016 vector<float> etabins;
00017 if (etaRegion==Central) {
00018 for (int i=0;i<8;++i) etabins.push_back(-2.1+0.6*i);
00019 } else if (etaRegion==Forward) {
00020 for (int i=0;i<10;++i) {
00021 if (i<5) etabins.push_back(-5.1+0.6*i);
00022 else etabins.push_back(2.7+0.6*(i-5));
00023 }
00024 } else if (etaRegion==All) {
00025 for (int i=0;i<18;++i) etabins.push_back(-5.1+0.6*i);
00026 }
00027 return fixedGridRho(etabins,phibins);
00028 }
00029
00030
00031 float FixedGridEnergyDensity::fixedGridRho(std::vector<float>& etabins,std::vector<float>& phibins) {
00032 float etadist = etabins[1]-etabins[0];
00033 float phidist = phibins[1]-phibins[0];
00034 float etahalfdist = (etabins[1]-etabins[0])/2.;
00035 float phihalfdist = (phibins[1]-phibins[0])/2.;
00036 vector<float> sumPFNallSMDQ;
00037 sumPFNallSMDQ.reserve(etabins.size()*phibins.size());
00038 for (unsigned int ieta=0;ieta<etabins.size();++ieta) {
00039 for (unsigned int iphi=0;iphi<phibins.size();++iphi) {
00040 float pfniso_ieta_iphi = 0;
00041 for(PFCandidateCollection::const_iterator pf_it = pfCandidates->begin(); pf_it != pfCandidates->end(); pf_it++) {
00042 if (fabs(etabins[ieta]-pf_it->eta())>etahalfdist) continue;
00043 if (fabs(reco::deltaPhi(phibins[iphi],pf_it->phi()))>phihalfdist) continue;
00044 pfniso_ieta_iphi+=pf_it->pt();
00045 }
00046 sumPFNallSMDQ.push_back(pfniso_ieta_iphi);
00047 }
00048 }
00049 float evt_smdq = 0;
00050 sort(sumPFNallSMDQ.begin(),sumPFNallSMDQ.end());
00051 if (sumPFNallSMDQ.size()%2) evt_smdq = sumPFNallSMDQ[(sumPFNallSMDQ.size()-1)/2];
00052 else evt_smdq = (sumPFNallSMDQ[sumPFNallSMDQ.size()/2]+sumPFNallSMDQ[(sumPFNallSMDQ.size()-2)/2])/2.;
00053 return evt_smdq/(etadist*phidist);
00054 }
1.7.3