#include "l1l2_evt/L1L2Chunk.hpp"
namespace {
// One L1 (2x2) tower
struct L1L2 {
// tower eta (-4.:4.)
float l1cal_tot_eta;
// tower phi (0.:2pi)
float l1cal_tot_phi;
// tower energy (EM+FH)
float l1cal_tot;
};
}
void analyzeEvent(const edm::Event& event) {
bool isMC = fwk::HistorySelector::is_monte_carlo(event);
// extract L1L2 chunk with L1 calorimeter towers
TKey l1l2Key;
THandle l1l2Chunk = l1l2Key.find(event);
if ( !l1l2Chunk.isValid() ) {
// if chunk is not available we can not apply new JetID
return;
}
// reserve an array to be filled with L1 cal towers
// Do it only once per event
std::vector l1infos;
// maximum number of towers is 1280
l1infos.reserve(1280);
l1cal_reco l1cal = (l1l2Chunk.ptr())->Ret_l1cal();
int tot_size = l1cal.l1cal_emtot_size();
// towers in l1cal_reco are sorted by energy
// We loop over only 100 the most energetic towers
// We fill l1infos with L1 towers.
for(int i=0 ; i < tot_size && i < 100 ; ++i ) {
L1L2 l1_tower;
l1_tower.l1cal_tot_eta = l1cal.l1cal_tot_eta(i);
l1_tower.l1cal_tot_phi = l1cal.l1cal_tot_phi(i);
l1_tower.l1cal_tot = l1cal.l1cal_tot(i);
l1infos.push_back(l1_tower);
}
// Here I don't show how to construct JetChunkSelector _jetCSel
// it depends on what kind of jets you want
const edm::TKey jetKey(_jetCSel);
edm::THandle jetHandle = jetKey.find(event);
if ( !jetHandle.isValid() ) {
return;
}
// Get list of your jets
const std::vector& jets = *(jetChunk->getParticles());
// Now we loop over list of jets.
for(int ijet = 0 ; ijet < jets.size() ; ++ijet ) {
const Jet& jet = jets[ijet];
//calculate l1set == sum of L1 tower energies
// for towers with 0.5 of the Jet direction
float l1set=0.;
for(int il1t = 0; il1t < l1infos.size() ; il1t++)
{
const L1L2& l1info = l1infos[il1t];
float de = fabs( jet.detEta()*0.1 - l1info.l1cal_tot_eta);
float DetPhi = (jet.detPhi() + 0.05) /64.0*2.0*kinem::PI;
float dp = kinem::delta_phi(DetPhi, l1info.l1cal_tot_phi);
float dr = sqrt (de*de+dp*dp);
if (dr < 0.5)
l1set += l1info.l1cal_tot;
}
// l1set is calculated
// determine if Jet is CC,ICD or EC jet
bool CC=false,EC=false,ICD=false;
if(abs(jet.detEta()*0.1) <= 0.8 ) CC=true;
if(abs(jet.detEta()*0.1) > 0.8 ) ICD=true;
if(abs(jet.detEta()*0.1) > 1.5 ) {EC=true;ICD=false;}
// if JES are applied, unapply
float C=1.;
if( fabs(jet.JES_C())>1e-7 ) C = jet.JES_C();
// jtpt is not corrected ( no JES) Pt of the Jet without
// CH energy
float jtpt = jet.pT()*(1.-jet.chETfraction())/C;
// float ratio -- is ratio between sum of L1 towers
// within 0.5 cone of the jet to the JetPT ( uncorrected, no CH)
float ratio = 0.;
if( jtpt > 0. )
ratio = l1set/jtpt;
// jet's standard quality criteria
if (jet.emETfraction() >= 0.95 ) return continue;
if (jet.emETfraction() <= 0.05 ) return continue;
if (jet.chETfraction() >= 0.4 ) return continue;
if (jet.hotcellratio() >= 10. ) return continue;
if (jet.n90() <= 0 ) return continue;
if (jet.E() <= 0.0) return continue;
// L1 confirmation. We don't apply L1 confirmation for MC
if( CC && ratio < 0.4 && !MC ) return continue;
if( ICD && ratio < 0.2 && !MC ) return continue;
if( EC && ratio < 0.4 && !MC ) return continue;
// O.K. Jets they made it here are "good" jets
}
}