#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

   } 

}