L1TStage2CaloLayer2Comp.cc

Go to the documentation of this file.

#ifndef L1Trigger_L1TCalorimeter_L1TStage2CaloLayer2Comp
#define L1Trigger_L1TCalorimeter_L1TStage2CaloLayer2Comp

#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/L1Trigger/interface/EGamma.h"
#include "DataFormats/L1Trigger/interface/Jet.h"
#include "DataFormats/L1Trigger/interface/EtSum.h"
#include "DataFormats/L1Trigger/interface/Tau.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/EDProducer.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "CondFormats/L1TObjects/interface/CaloParams.h"
#include "CondFormats/DataRecord/interface/L1TCaloParamsRcd.h"

#include "DataFormats/L1TCalorimeter/interface/CaloTower.h"
#include "DataFormats/L1TCalorimeter/interface/CaloCluster.h"
#include "DataFormats/L1Trigger/interface/EGamma.h"
#include "DataFormats/L1Trigger/interface/Tau.h"
#include "DataFormats/L1Trigger/interface/Jet.h"
#include "DataFormats/L1Trigger/interface/EtSum.h"

// classes required for storage of objects into a separate edm file
#include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"

// classes required for sorting of jets
#include "L1Trigger/L1TCalorimeter/interface/Stage2Layer2JetAlgorithmFirmware.h"
#include "L1Trigger/L1TCalorimeter/interface/AccumulatingSort.h"
#include "L1Trigger/L1TCalorimeter/interface/BitonicSort.h"
#include "DataFormats/Math/interface/LorentzVector.h"


namespace l1t {
  bool operator > ( const l1t::Jet& a, l1t::Jet& b ) {
    return  a.hwPt() > b.hwPt();
  }
  bool operator > ( const l1t::EGamma& a, l1t::EGamma& b ) {
    return  a.hwPt() > b.hwPt();
  }
  bool operator > ( const l1t::Tau& a, l1t::Tau& b ) {
    return  a.hwPt() > b.hwPt();
  }
}


#include "TH1F.h"
#include "TH2F.h"


class  L1TStage2CaloLayer2Comp : public edm::EDProducer {

public:
  L1TStage2CaloLayer2Comp (const edm::ParameterSet & ps);


protected:

  void produce (edm::Event&, const edm::EventSetup&) override;

private:

  bool compareJets(const edm::Handle<l1t::JetBxCollection> & dataCol,
                   const edm::Handle<l1t::JetBxCollection> & emulCol);

  bool compareEGs(const edm::Handle<l1t::EGammaBxCollection> & dataCol,
                  const edm::Handle<l1t::EGammaBxCollection> & emulCol);

  bool compareTaus(const edm::Handle<l1t::TauBxCollection> & dataCol,
                   const edm::Handle<l1t::TauBxCollection> & emulCol);

  bool compareSums(const edm::Handle<l1t::EtSumBxCollection> & dataCol,
                   const edm::Handle<l1t::EtSumBxCollection> & emulCol);

  void accuSort(std::vector<l1t::Jet> *jets);
  void accuSort(std::vector<l1t::Tau> *taus);
  void accuSort(std::vector<l1t::EGamma> *egs);
  void accuSort(std::vector<l1t::L1Candidate> & jets);

  void dumpEventToFile();
  void dumpEventToEDM(edm::Event &e);

  // collections to hold entities reconstructed from data and emulation
  edm::EDGetToken calol2JetCollectionData;
  edm::EDGetToken calol2JetCollectionEmul;
  edm::EDGetToken calol2EGammaCollectionData;
  edm::EDGetToken calol2EGammaCollectionEmul;
  edm::EDGetToken calol2TauCollectionData;
  edm::EDGetToken calol2TauCollectionEmul;
  edm::EDGetToken calol2EtSumCollectionData;
  edm::EDGetToken calol2EtSumCollectionEmul;
  edm::EDGetToken calol2CaloTowerCollectionData;
  edm::EDGetToken calol2CaloTowerCollectionEmul;

  // define collections to hold lists of objects in event
  edm::Handle<l1t::JetBxCollection> jetDataCol;
  edm::Handle<l1t::JetBxCollection> jetEmulCol;
  edm::Handle<l1t::EGammaBxCollection> egDataCol;
  edm::Handle<l1t::EGammaBxCollection> egEmulCol;
  edm::Handle<l1t::TauBxCollection> tauDataCol;
  edm::Handle<l1t::TauBxCollection> tauEmulCol;
  edm::Handle<l1t::EtSumBxCollection> sumDataCol;
  edm::Handle<l1t::EtSumBxCollection> sumEmulCol;
  edm::Handle<l1t::CaloTowerBxCollection> caloTowerDataCol;
  edm::Handle<l1t::CaloTowerBxCollection> caloTowerEmulCol;

  const unsigned int currBx = 0;

  bool dumpTowers = false;
  bool dumpWholeEvent = false;
};

L1TStage2CaloLayer2Comp::L1TStage2CaloLayer2Comp (const edm::ParameterSet& ps)
  : calol2JetCollectionData(consumes <l1t::JetBxCollection>(
                  ps.getParameter<edm::InputTag>(
                "calol2JetCollectionData"))),
    calol2JetCollectionEmul(consumes <l1t::JetBxCollection>(
                  ps.getParameter<edm::InputTag>(
                "calol2JetCollectionEmul"))),
    calol2EGammaCollectionData(consumes <l1t::EGammaBxCollection>(
                 ps.getParameter<edm::InputTag>(
                   "calol2EGammaCollectionData"))),
    calol2EGammaCollectionEmul(consumes <l1t::EGammaBxCollection>(
                 ps.getParameter<edm::InputTag>(
                   "calol2EGammaCollectionEmul"))),
    calol2TauCollectionData(consumes <l1t::TauBxCollection>(
                  ps.getParameter<edm::InputTag>(
                "calol2TauCollectionData"))),
    calol2TauCollectionEmul(consumes <l1t::TauBxCollection>(
                  ps.getParameter<edm::InputTag>(
                "calol2TauCollectionEmul"))),
    calol2EtSumCollectionData(consumes <l1t::EtSumBxCollection>(
                ps.getParameter<edm::InputTag>(
                  "calol2EtSumCollectionData"))),
    calol2EtSumCollectionEmul(consumes <l1t::EtSumBxCollection>(
                ps.getParameter<edm::InputTag>(
                  "calol2EtSumCollectionEmul"))),
    calol2CaloTowerCollectionData(consumes <l1t::CaloTowerBxCollection>(
                    ps.getParameter<edm::InputTag>(
                      "calol2CaloTowerCollectionData"))),
    calol2CaloTowerCollectionEmul(consumes <l1t::CaloTowerBxCollection>(
                    ps.getParameter<edm::InputTag>(
                      "calol2CaloTowerCollectionEmul")))
{
  produces<l1t::JetBxCollection> ("dataJet").setBranchAlias("dataJets");
  produces<l1t::JetBxCollection> ("emulJet").setBranchAlias("emulJets");
  produces<l1t::EGammaBxCollection> ("dataEg").setBranchAlias("dataEgs");
  produces<l1t::EGammaBxCollection> ("emulEg").setBranchAlias("emulEgs");
  produces<l1t::TauBxCollection> ("dataTau").setBranchAlias("dataTaus");
  produces<l1t::TauBxCollection> ("emulTau").setBranchAlias("emulTaus");
  produces<l1t::EtSumBxCollection> ("dataEtSum").setBranchAlias("dataEtSums");
  produces<l1t::EtSumBxCollection> ("emulEtSum").setBranchAlias("emulEtSums");
  produces<l1t::CaloTowerBxCollection> ("dataCaloTower").setBranchAlias("dataCaloTowers");
  produces<l1t::CaloTowerBxCollection> ("emulCaloTower").setBranchAlias("emulCaloTowers");

  dumpTowers = ps.getParameter<bool>("dumpTowers");
  dumpWholeEvent = ps.getParameter<bool>("dumpWholeEvent");
}

void L1TStage2CaloLayer2Comp::produce (
  edm::Event& e,
  const edm::EventSetup & c) {

  bool eventGood = true;

  // map event contents to above collections
  e.getByToken(calol2JetCollectionData, jetDataCol);
  e.getByToken(calol2JetCollectionEmul, jetEmulCol);
  e.getByToken(calol2EGammaCollectionData, egDataCol);
  e.getByToken(calol2EGammaCollectionEmul, egEmulCol);
  e.getByToken(calol2TauCollectionData, tauDataCol);
  e.getByToken(calol2TauCollectionEmul, tauEmulCol);
  e.getByToken(calol2EtSumCollectionData, sumDataCol);
  e.getByToken(calol2EtSumCollectionEmul, sumEmulCol);
  e.getByToken(calol2CaloTowerCollectionData, caloTowerDataCol);
  e.getByToken(calol2CaloTowerCollectionEmul, caloTowerEmulCol);

  edm::LogProblem("l1tcalol2ebec")
    << "Processing event " << e.id() << std::endl;

  if(caloTowerDataCol->size()==0){
    edm::LogProblem("l1tcalol2ebec")
      << "Empty caloTowers. Skipping event " << std::endl;
    return;
  }

  if (!compareJets(jetDataCol, jetEmulCol)) {
    eventGood = false;
  }

  if (!compareEGs(egDataCol, egEmulCol)) {
    eventGood = false;
  }

  if (!compareTaus(tauDataCol, tauEmulCol)) {
    eventGood = false;
  }

  if (!compareSums(sumDataCol, sumEmulCol)) {
    eventGood = false;
  }

  if (!eventGood) {

    if (dumpWholeEvent) {
      // store all contents of event to log file
      dumpEventToFile();
    }

    // store all contents of event to edm file:
    dumpEventToEDM(e);

    edm::LogProblem("l1tcalol2ebec")
    << "Bad event! " << std::endl;

  }
}

// comparison method for jets
bool L1TStage2CaloLayer2Comp::compareJets(
  const edm::Handle<l1t::JetBxCollection> & dataCol,
  const edm::Handle<l1t::JetBxCollection> & emulCol)
{
  bool eventGood = true;

  std::vector<l1t::Jet> *jets = new std::vector<l1t::Jet>;
  std::vector<l1t::Jet>::iterator dataIt;
  l1t::JetBxCollection::const_iterator dataBxIt = dataCol->begin(currBx);
  l1t::JetBxCollection::const_iterator emulBxIt = emulCol->begin(currBx);

  if (dataCol->size(currBx) > 0) {
    // sort data jets
    while (true) {

      jets->emplace_back(*dataBxIt);

      dataBxIt++;
      if (dataBxIt == dataCol->end(currBx))
    break;
    }

    accuSort(jets);
    dataIt = jets->begin();
  }

  // process jets
  if (jets->size() != emulCol->size(currBx)) {
    edm::LogProblem("l1tcalol2ebec")
      << "Jet collection size difference: "
      << "data = " << jets->size()
      << " emul = " << emulCol->size(currBx)
      << std::endl;
    return false;
  }

  if (dataIt != jets->end() ||
      emulBxIt != emulCol->end(currBx)) {

    int nPos=0;
    int nNeg=0;

    while(true) {

      bool posGood = true;
      bool etGood = true;

      // object pt mismatch
      if (dataIt->hwPt() != emulBxIt->hwPt()) {
        etGood = false;
    eventGood = false;
      }

      // object position mismatch (phi)
      if (dataIt->hwPhi() != emulBxIt->hwPhi()){
    posGood = false;
      }

      // object position mismatch (eta)
      if (dataIt->hwEta() != emulBxIt->hwEta()) {
    posGood = false;
      }

      //bypass sorting bug
      if(etGood && !posGood){
        l1t::JetBxCollection::const_iterator emulItCheckSort;
        l1t::JetBxCollection::const_iterator dataItCheckSort;
        for(emulItCheckSort = emulCol->begin(currBx) ; emulItCheckSort != emulCol->end(currBx) ; ++emulItCheckSort) {
          for (dataItCheckSort = jets->begin(); dataItCheckSort != jets->end(); ++dataItCheckSort){
            
            if(dataItCheckSort->hwEta() > 0) ++nPos;
            if(dataItCheckSort->hwEta() < 0) ++nNeg;
            
            if(dataIt->hwPt() == emulItCheckSort->hwPt()
               && dataIt->hwPhi() == emulItCheckSort->hwPhi()
               && dataIt->hwEta() == emulItCheckSort->hwEta())
              posGood = true;
          }
        }
      }

      if(etGood && dataIt->hwEta() > 0 && ((distance(dataIt, jets->end())-nNeg) < 5))
    posGood = true;
      if(etGood && dataIt->hwEta() < 0 && (distance(dataIt, jets->end()) < 5))
    posGood = true;

      if(!posGood) eventGood = false;

      // if both position and energy agree, object is good
      if (!etGood || !posGood) {
    edm::LogProblem("l1tcalol2ebec")
      << "Jet Problem (data emul): "
      << "\tEt = " << dataIt->hwPt() << " " << emulBxIt->hwPt()
      << "\teta = " << dataIt->hwEta() << " " << emulBxIt->hwEta()
      << "\tphi = " << dataIt->hwPhi() << " " << emulBxIt->hwPhi()
      << std::endl;
      }

      // increment position of pointers
      ++dataIt;
      ++emulBxIt;

      if (dataIt == jets->end() ||
      emulBxIt == emulCol->end(currBx))
    break;
    }
  } else {
    if (!jets->empty() || emulCol->size(currBx) != 0)
      return false;
  }

  // return a boolean that states whether the jet data in the event are in
  // agreement
  return eventGood;
}

// comparison method for e/gammas
bool L1TStage2CaloLayer2Comp::compareEGs(
  const edm::Handle<l1t::EGammaBxCollection> & dataCol,
  const edm::Handle<l1t::EGammaBxCollection> & emulCol)
{
  bool eventGood = true;

  std::vector<l1t::EGamma> *egs = new std::vector<l1t::EGamma>;
  std::vector<l1t::EGamma>::iterator dataIt;
  l1t::EGammaBxCollection::const_iterator dataBxIt = dataCol->begin(currBx);
  l1t::EGammaBxCollection::const_iterator emulBxIt = emulCol->begin(currBx);

  if (dataCol->size(currBx) > 0) {
    // sort data egs
    while (true) {

      egs->emplace_back(*dataBxIt);

      dataBxIt++;
      if (dataBxIt == dataCol->end(currBx))
    break;
    }
    accuSort(egs);
    dataIt = egs->begin();
  }

  // check length of collections
  if (egs->size() != emulCol->size(currBx)) {
    edm::LogProblem("l1tcalol2ebec")
      << "EG collection size difference: "
      << "data = " << egs->size()
      << " emul = " << emulCol->size(currBx)
      << std::endl;
    return false;
  }

  // processing continues only of length of data collections is the same
  if (dataIt != egs->end() ||
      emulBxIt != emulCol->end(currBx)) {

    int nPos=0;
    int nNeg=0;

    while(true) {

      bool posGood = true;
      bool etGood = true;

      // object pt mismatch
      if (dataIt->hwPt() != emulBxIt->hwPt()) {
    etGood = false;
    eventGood = false;
      }

      // object position mismatch (phi)
      if (dataIt->hwPhi() != emulBxIt->hwPhi()) {
    posGood = false;
      }

      // object position mismatch (eta)
      if (dataIt->hwEta() != emulBxIt->hwEta()) {
    posGood = false;
      }

      //bypass sorting bug
      if(etGood && !posGood){
        l1t::EGammaBxCollection::const_iterator emulItCheckSort;
        l1t::EGammaBxCollection::const_iterator dataItCheckSort;
        for(emulItCheckSort = emulCol->begin(currBx); emulItCheckSort != emulCol->end(currBx); ++emulItCheckSort) {
          for (dataItCheckSort = egs->begin(); dataItCheckSort != egs->end(); ++dataItCheckSort){
            
            if(dataItCheckSort->hwEta() > 0) ++nPos;
            if(dataItCheckSort->hwEta() < 0) ++nNeg;
            
            if(dataIt->hwPt() == emulItCheckSort->hwPt()
               && dataIt->hwPhi() == emulItCheckSort->hwPhi()
               && dataIt->hwEta() == emulItCheckSort->hwEta())
              posGood = true;
          }
        }
      }

      if(etGood && dataIt->hwEta() > 0 && ((distance(dataIt, egs->end())-nNeg) < 5))
    posGood = true;
      if(etGood && dataIt->hwEta() < 0 && (distance(dataIt, egs->end()) < 5))
    posGood = true;

      if(!posGood) eventGood = false;

      // if both position and energy agree, object is good
      if (!posGood || !etGood) {
    edm::LogProblem("l1tcalol2ebec")
      << "EG Problem (data emul): "
      << "\tEt = " << dataIt->hwPt() << " " << emulBxIt->hwPt()
      << "\teta = " << dataIt->hwEta() << " " << emulBxIt->hwEta()
      << "\tphi = " << dataIt->hwPhi() << " " << emulBxIt->hwPhi()
      << std::endl;
      }

      // increment position of pointers
      ++dataIt;
      ++emulBxIt;

      if (dataIt == egs->end() ||
      emulBxIt == emulCol->end(currBx))
    break;
    }
  } else {
    if (!egs->empty() || emulCol->size(currBx) != 0)
      return false;
  }

  // return a boolean that states whether the eg data in the event are in
  // agreement
  return eventGood;
}

// comparison method for taus
bool L1TStage2CaloLayer2Comp::compareTaus(
  const edm::Handle<l1t::TauBxCollection> & dataCol,
  const edm::Handle<l1t::TauBxCollection> & emulCol)
{
  bool eventGood = true;

  std::vector<l1t::Tau> *taus = new std::vector<l1t::Tau>;
  std::vector<l1t::Tau>::iterator dataIt;
  l1t::TauBxCollection::const_iterator dataBxIt = dataCol->begin(currBx);
  l1t::TauBxCollection::const_iterator emulBxIt = emulCol->begin(currBx);

  if (dataCol->size(currBx) > 0) {
    // sort data taus
    while (true) {

      taus->emplace_back(*dataBxIt);

      dataBxIt++;
      if (dataBxIt == dataCol->end(currBx))
    break;
    }
    accuSort(taus);
    dataIt = taus->begin();
  }

  // check length of collections
  if (taus->size() != emulCol->size(currBx)) {
    edm::LogProblem("l1tcalol2ebec")
      << "Tau collection size difference: "
      << "data = " << taus->size()
      << " emul = " << emulCol->size(currBx)
      << std::endl;
    return false;
  }

  // processing continues only of length of data collections is the same
  if (dataIt != taus->end() ||
      emulBxIt != emulCol->end(currBx)) {

    int nPos=0;
    int nNeg=0;

    while(true) {

      bool posGood = true;
      bool etGood = true;

      // object Et mismatch
      if (dataIt->hwPt() != emulBxIt->hwPt()) {
    etGood = false;
    eventGood = false;
      }

      // object position mismatch (phi)
      if (dataIt->hwPhi() != emulBxIt->hwPhi()) {
    posGood = false;
      }

      // object position mismatch (eta)
      if (dataIt->hwEta() != emulBxIt->hwEta()) {
    posGood = false;
      }

      //bypass sorting bug
      if(etGood && !posGood){
        l1t::TauBxCollection::const_iterator emulItCheckSort;
        l1t::TauBxCollection::const_iterator dataItCheckSort;
        for(emulItCheckSort = emulCol->begin(currBx); emulItCheckSort != emulCol->end(currBx); ++emulItCheckSort) {
          for (dataItCheckSort = taus->begin(); dataItCheckSort != taus->end(); ++dataItCheckSort){
          
            if(dataItCheckSort->hwEta() > 0) ++nPos;
            if(dataItCheckSort->hwEta() < 0) ++nNeg;
            
            if(dataIt->hwPt() == emulItCheckSort->hwPt()
               && dataIt->hwPhi() == emulItCheckSort->hwPhi()
               && dataIt->hwEta() == emulItCheckSort->hwEta())
              posGood = true;
          }
        }
      }

      if(etGood && dataIt->hwEta() > 0 && ((distance(dataIt, taus->end())-nNeg) < 5))
    posGood = true;
      if(etGood && dataIt->hwEta() < 0 && (distance(dataIt, taus->end()) < 5))
    posGood = true;

      if(!posGood) eventGood = false;

      // if both position and energy agree, object is good
      if (!posGood || !etGood) {
    edm::LogProblem("l1tcalol2ebec")
      << "Tau Problem (data emul): "
      << "\tEt = " << dataIt->hwPt() << " " << emulBxIt->hwPt()
      << "\teta = " << dataIt->hwEta() << " " << emulBxIt->hwEta()
      << "\tphi = " << dataIt->hwPhi() << " " << emulBxIt->hwPhi()
      << std::endl;
      }

      // increment position of pointers
      ++dataIt;
      ++emulBxIt;

      if (dataIt == taus->end() ||
      emulBxIt == emulCol->end(currBx))
    break;
    }
  } else {
    if (!taus->empty() || emulCol->size(currBx) != 0)
      return false;
  }

  // return a boolean that states whether the tau data in the event are in
  // agreement
  return eventGood;
}

// comparison method for sums
bool L1TStage2CaloLayer2Comp::compareSums(
  const edm::Handle<l1t::EtSumBxCollection> & dataCol,
  const edm::Handle<l1t::EtSumBxCollection> & emulCol)
{
  bool eventGood = true;

  double dataEt = 0;
  double emulEt = 0;

  // if either data or emulator collections are empty, or they have different
  // size, mark the event as bad (this should never occur in normal running)
  if ((dataCol->isEmpty(currBx)) || (emulCol->isEmpty(currBx)) ||
      (dataCol->size(currBx) != emulCol->size(currBx))) {
    edm::LogProblem("l1tcalol2ebec")
      << "EtSum collection size difference: "
      << "data = " << dataCol->size(currBx)
      << " emul = " << emulCol->size(currBx)
      << std::endl;

    return false;
  }

  l1t::EtSum dataSum;
  l1t::EtSum emulSum;
  for (unsigned int i = 0; i < emulCol->size(currBx); i++) {

    emulSum = emulCol->at(currBx, i);
    dataSum = dataCol->at(currBx, l1t::CaloTools::emul_to_data_sum_index_map[i]);

    if (emulSum.getType() != dataSum.getType()) {
      edm::LogProblem("l1tcalol2ebec")
    << "EtSum type problem (data emul): "
    << dataSum.getType() << " " << emulSum.getType()
    << std::endl;
    }

    dataEt = dataSum.hwPt();
    emulEt = emulSum.hwPt();

    if (dataEt != emulEt) {
      eventGood = false;
      edm::LogProblem("l1tcalol2ebec")
    << "EtSum problem (data emul):\tType = " << emulSum.getType()
    << "\tEt = " << dataEt << " " << emulEt
    << std::endl;
    }
  }

  // return a boolean that states whether the sum data in the event are in
  // agreement
  return eventGood;
}

// sort for jets
void L1TStage2CaloLayer2Comp::accuSort(std::vector<l1t::Jet> *jets){

  math::PtEtaPhiMLorentzVector emptyP4;
  l1t::Jet tempJet (emptyP4, 0, 0, 0, 0);
  std::vector< std::vector<l1t::Jet> > jetEtaPos( 41 , std::vector<l1t::Jet>(18, tempJet));
  std::vector< std::vector<l1t::Jet> > jetEtaNeg( 41 , std::vector<l1t::Jet>(18, tempJet));

  for (unsigned int iJet = 0; iJet < jets->size(); iJet++) {
    if (jets->at(iJet).hwEta() > 0) jetEtaPos.at(jets->at(iJet).hwEta()-1).at((jets->at(iJet).hwPhi()-1)/4) = jets->at(iJet);
    else  jetEtaNeg.at(-(jets->at(iJet).hwEta()+1)).at((jets->at(iJet).hwPhi()-1)/4) = jets->at(iJet);
  }

  AccumulatingSort <l1t::Jet> etaPosSorter(7);
  AccumulatingSort <l1t::Jet> etaNegSorter(7);
  std::vector<l1t::Jet> accumEtaPos;
  std::vector<l1t::Jet> accumEtaNeg;

  for( int ieta = 0 ; ieta < 41 ; ++ieta) {
    // eta +
    std::vector<l1t::Jet>::iterator start_, end_;
    start_ = jetEtaPos.at(ieta).begin();
    end_   = jetEtaPos.at(ieta).end();
    BitonicSort<l1t::Jet>(down, start_, end_);
    etaPosSorter.Merge( jetEtaPos.at(ieta) , accumEtaPos );

    // eta -
    start_ = jetEtaNeg.at(ieta).begin();
    end_   = jetEtaNeg.at(ieta).end();
    BitonicSort<l1t::Jet>(down, start_, end_);
    etaNegSorter.Merge( jetEtaNeg.at(ieta) , accumEtaNeg );

  }

  //check for 6 & 7th jets with same et and eta. Keep jet with larger phi
  if(accumEtaPos.at(6).hwPt()==accumEtaPos.at(5).hwPt() && accumEtaPos.at(6).hwEta()==accumEtaPos.at(5).hwEta()
     && accumEtaPos.at(6).hwPhi() > accumEtaPos.at(5).hwPhi()){
    accumEtaPos.at(5)=accumEtaPos.at(6);
  }
  if(accumEtaNeg.at(6).hwPt()==accumEtaNeg.at(5).hwPt() && accumEtaNeg.at(6).hwEta()==accumEtaNeg.at(5).hwEta()
     && accumEtaNeg.at(6).hwPhi() > accumEtaNeg.at(5).hwPhi()){
    accumEtaNeg.at(5)=accumEtaNeg.at(6);
  }

  //truncate
  accumEtaPos.resize(6);
  accumEtaNeg.resize(6);

  // put all 12 candidates in the original jet vector, removing zero energy ones
  jets->clear();
  for (l1t::Jet accjet : accumEtaPos) {
    if (accjet.hwPt() > 0) jets->push_back(accjet);
  }
  for (l1t::Jet accjet : accumEtaNeg) {
    if (accjet.hwPt() > 0) jets->push_back(accjet);
  }
}

// sort for eg
void L1TStage2CaloLayer2Comp::accuSort(std::vector<l1t::EGamma> *egs){

  math::PtEtaPhiMLorentzVector emptyP4;
  l1t::EGamma tempEGamma (emptyP4, 0, 0, 0, 0);
  std::vector< std::vector<l1t::EGamma> > jetEtaPos( 41 , std::vector<l1t::EGamma>(18, tempEGamma));
  std::vector< std::vector<l1t::EGamma> > jetEtaNeg( 41 , std::vector<l1t::EGamma>(18, tempEGamma));

  for (unsigned int iEGamma = 0; iEGamma < egs->size(); iEGamma++) {
    if (egs->at(iEGamma).hwEta() > 0) jetEtaPos.at(egs->at(iEGamma).hwEta()-1).at((egs->at(iEGamma).hwPhi()-1)/4) = egs->at(iEGamma);
    else  jetEtaNeg.at(-(egs->at(iEGamma).hwEta()+1)).at((egs->at(iEGamma).hwPhi()-1)/4) = egs->at(iEGamma);
  }

  AccumulatingSort <l1t::EGamma> etaPosSorter(7);
  AccumulatingSort <l1t::EGamma> etaNegSorter(7);
  std::vector<l1t::EGamma> accumEtaPos;
  std::vector<l1t::EGamma> accumEtaNeg;

  for( int ieta = 0 ; ieta < 41 ; ++ieta) {
    // eta +
    std::vector<l1t::EGamma>::iterator start_, end_;
    start_ = jetEtaPos.at(ieta).begin();
    end_   = jetEtaPos.at(ieta).end();
    BitonicSort<l1t::EGamma>(down, start_, end_);
    etaPosSorter.Merge( jetEtaPos.at(ieta) , accumEtaPos );

    // eta -
    start_ = jetEtaNeg.at(ieta).begin();
    end_   = jetEtaNeg.at(ieta).end();
    BitonicSort<l1t::EGamma>(down, start_, end_);
    etaNegSorter.Merge( jetEtaNeg.at(ieta) , accumEtaNeg );

  }

  //check for 6 & 7th egs with same et and eta. Keep jet with larger phi
  if(accumEtaPos.at(6).hwPt()==accumEtaPos.at(5).hwPt() && accumEtaPos.at(6).hwEta()==accumEtaPos.at(5).hwEta()
     && accumEtaPos.at(6).hwPhi() > accumEtaPos.at(5).hwPhi()){
    accumEtaPos.at(5)=accumEtaPos.at(6);
  }
  if(accumEtaNeg.at(6).hwPt()==accumEtaNeg.at(5).hwPt() && accumEtaNeg.at(6).hwEta()==accumEtaNeg.at(5).hwEta()
     && accumEtaNeg.at(6).hwPhi() > accumEtaNeg.at(5).hwPhi()){
    accumEtaNeg.at(5)=accumEtaNeg.at(6);
  }

  //truncate
  accumEtaPos.resize(6);
  accumEtaNeg.resize(6);

  // put all 12 candidates in the original jet vector, removing zero energy ones
  egs->clear();
  for (l1t::EGamma accjet : accumEtaPos) {
    if (accjet.hwPt() > 0) egs->push_back(accjet);
  }
  for (l1t::EGamma accjet : accumEtaNeg) {
    if (accjet.hwPt() > 0) egs->push_back(accjet);
  }
}


// sort for tau
void L1TStage2CaloLayer2Comp::accuSort(std::vector<l1t::Tau> *taus){

  math::PtEtaPhiMLorentzVector emptyP4;
  l1t::Tau tempTau (emptyP4, 0, 0, 0, 0);
  std::vector< std::vector<l1t::Tau> > jetEtaPos( 41 , std::vector<l1t::Tau>(18, tempTau));
  std::vector< std::vector<l1t::Tau> > jetEtaNeg( 41 , std::vector<l1t::Tau>(18, tempTau));

  for (unsigned int iTau = 0; iTau < taus->size(); iTau++) {
    if (taus->at(iTau).hwEta() > 0) jetEtaPos.at(taus->at(iTau).hwEta()-1).at((taus->at(iTau).hwPhi()-1)/4) = taus->at(iTau);
    else  jetEtaNeg.at(-(taus->at(iTau).hwEta()+1)).at((taus->at(iTau).hwPhi()-1)/4) = taus->at(iTau);
  }

  AccumulatingSort <l1t::Tau> etaPosSorter(7);
  AccumulatingSort <l1t::Tau> etaNegSorter(7);
  std::vector<l1t::Tau> accumEtaPos;
  std::vector<l1t::Tau> accumEtaNeg;

  for( int ieta = 0 ; ieta < 41 ; ++ieta) {
    // eta +
    std::vector<l1t::Tau>::iterator start_, end_;
    start_ = jetEtaPos.at(ieta).begin();
    end_   = jetEtaPos.at(ieta).end();
    BitonicSort<l1t::Tau>(down, start_, end_);
    etaPosSorter.Merge( jetEtaPos.at(ieta) , accumEtaPos );

    // eta -
    start_ = jetEtaNeg.at(ieta).begin();
    end_   = jetEtaNeg.at(ieta).end();
    BitonicSort<l1t::Tau>(down, start_, end_);
    etaNegSorter.Merge( jetEtaNeg.at(ieta) , accumEtaNeg );

  }

  //check for 6 & 7th taus with same et and eta. Keep jet with larger phi
  if(accumEtaPos.at(6).hwPt()==accumEtaPos.at(5).hwPt() && accumEtaPos.at(6).hwEta()==accumEtaPos.at(5).hwEta()
     && accumEtaPos.at(6).hwPhi() > accumEtaPos.at(5).hwPhi()){
    accumEtaPos.at(5)=accumEtaPos.at(6);
  }
  if(accumEtaNeg.at(6).hwPt()==accumEtaNeg.at(5).hwPt() && accumEtaNeg.at(6).hwEta()==accumEtaNeg.at(5).hwEta()
     && accumEtaNeg.at(6).hwPhi() > accumEtaNeg.at(5).hwPhi()){
    accumEtaNeg.at(5)=accumEtaNeg.at(6);
  }

  //truncate
  accumEtaPos.resize(6);
  accumEtaNeg.resize(6);

  // put all 12 candidates in the original jet vector, removing zero energy ones
  taus->clear();
  for (l1t::Tau accjet : accumEtaPos) {
    if (accjet.hwPt() > 0) taus->push_back(accjet);
  }
  for (l1t::Tau accjet : accumEtaNeg) {
    if (accjet.hwPt() > 0) taus->push_back(accjet);
  }
}

void L1TStage2CaloLayer2Comp::dumpEventToFile()
{
    edm::LogProblem("l1tcalol2ebec")
      << "==== Problems found, dumping full event contents ====" << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== Event contents in data: ====" << std::endl;

    if (dumpTowers) {
      edm::LogProblem("l1tcalol2ebec")
        << "==== Towers: ====" << std::endl;

      for (auto tower = caloTowerDataCol->begin(currBx); tower != caloTowerDataCol->end(currBx); ++tower)
    edm::LogProblem("l1tcalol2ebec")
      << "Tower: Et = " << tower->hwPt() << ", "
      << "eta = " << tower->hwEta() << ", "
      << "phi = " << tower->hwPhi() << std::endl;
    }

    edm::LogProblem("l1tcalol2ebec")
      << "==== Jets: ====" << std::endl;
    for (auto jet = jetDataCol->begin(currBx); jet != jetDataCol->end(currBx); ++jet)
      edm::LogProblem("l1tcalol2ebec")
    << "Jet: Et = " << jet->hwPt() << ", "
    << "eta = " << jet->hwEta() << ", "
    << "phi = " << jet->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== EGs: ====" << std::endl;
    for (auto eg = egDataCol->begin(currBx); eg != egDataCol->end(currBx); ++eg)
      edm::LogProblem("l1tcalol2ebec")
    << "EG: Et = " << eg->hwPt() << ", "
    << "eta = " << eg->hwEta() << ", "
    << "phi = " << eg->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== Taus: ====" << std::endl;
    for (auto tau = tauDataCol->begin(currBx); tau != tauDataCol->end(currBx); ++tau)
      edm::LogProblem("l1tcalol2ebec")
    << "Tau: Et = " << tau->hwPt() << ", "
    << "eta = " << tau->hwEta() << ", "
    << "phi = " << tau->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== Sums: ====" << std::endl;
    for (auto sum = sumDataCol->begin(currBx); sum != sumDataCol->end(currBx); ++sum)
      edm::LogProblem("l1tcalol2ebec")
    << "Sum: type = " << sum->getType() << " "
    << "Et = " << sum->hwPt() << ", "
    << "eta = " << sum->hwEta() << ", "
    << "phi = " << sum->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== Event contents in emul: ====" << std::endl;

    if (dumpTowers) {
      edm::LogProblem("l1tcalol2ebec")
    << "==== Towers: ====" << std::endl;

      for (auto tower = caloTowerEmulCol->begin(currBx); tower != caloTowerEmulCol->end(currBx); ++tower)
    edm::LogProblem("l1tcalol2ebec")
      << "Tower: Et = " << tower->hwPt() << ", "
      << "eta = " << tower->hwEta() << ", "
      << "phi = " << tower->hwPhi() << std::endl;
    }

    edm::LogProblem("l1tcalol2ebec")
      << "==== Jets: ====" << std::endl;
    for (auto jet = jetEmulCol->begin(currBx); jet != jetEmulCol->end(currBx); ++jet)
      edm::LogProblem("l1tcalol2ebec")
    << "Jet: Et = " << jet->hwPt() << ", "
    << "eta = " << jet->hwEta() << ", "
    << "phi = " << jet->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== EGs: ====" << std::endl;
    for (auto eg = egEmulCol->begin(currBx); eg != egEmulCol->end(currBx); ++eg)
      edm::LogProblem("l1tcalol2ebec")
    << "EG: Et = " << eg->hwPt() << ", "
    << "eta = " << eg->hwEta() << ", "
    << "phi = " << eg->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== Taus: ====" << std::endl;
    for (auto tau = tauEmulCol->begin(currBx); tau != tauEmulCol->end(currBx); ++tau)
      edm::LogProblem("l1tcalol2ebec")
    << "Tau: Et = " << tau->hwPt() << ", "
    << "eta = " << tau->hwEta() << ", "
    << "phi = " << tau->hwPhi() << std::endl;

    edm::LogProblem("l1tcalol2ebec")
      << "==== Sums: ====" << std::endl;
    for (auto sum = sumEmulCol->begin(currBx); sum != sumEmulCol->end(currBx); ++sum)
      edm::LogProblem("l1tcalol2ebec")
    << "Sum: type = " << sum->getType() << " "
    << "Et = " << sum->hwPt() << ", "
    << "eta = " << sum->hwEta() << ", "
    << "phi = " << sum->hwPhi() << std::endl;

}


void L1TStage2CaloLayer2Comp::dumpEventToEDM(edm::Event &e)
{
    // store all jets to an edm file
  std::unique_ptr<l1t::JetBxCollection> mpjets_data (new l1t::JetBxCollection(0, currBx, currBx));
  std::unique_ptr<l1t::JetBxCollection> mpjets_emul (new l1t::JetBxCollection(0, currBx, currBx));

    for (auto jet = jetDataCol->begin(currBx); jet != jetDataCol->end(currBx); ++jet)
      mpjets_data->push_back(0, (*jet));
    for (auto jet = jetEmulCol->begin(currBx); jet != jetEmulCol->end(currBx); ++jet)
      mpjets_emul->push_back(0, (*jet));

    e.put(std::move(mpjets_data), "dataJet");
    e.put(std::move(mpjets_emul), "emulJet");

    // store all egs to an edm file
    std::unique_ptr<l1t::EGammaBxCollection> mpEGammas_data (new l1t::EGammaBxCollection(0, currBx, currBx));
    std::unique_ptr<l1t::EGammaBxCollection> mpEGammas_emul (new l1t::EGammaBxCollection(0, currBx, currBx));

    for (auto eg = egDataCol->begin(currBx); eg != egDataCol->end(currBx); ++eg)
      mpEGammas_data->push_back(0, (*eg));
    for (auto eg = egEmulCol->begin(currBx); eg != egEmulCol->end(currBx); ++eg)
      mpEGammas_emul->push_back(0, (*eg));

    e.put(std::move(mpEGammas_data), "dataEg");
    e.put(std::move(mpEGammas_emul), "emulEg");

    // store all taus to an edm file
    std::unique_ptr<l1t::TauBxCollection> mptaus_data (new l1t::TauBxCollection(0, currBx, currBx));
    std::unique_ptr<l1t::TauBxCollection> mptaus_emul (new l1t::TauBxCollection(0, currBx, currBx));

    for (auto tau = tauDataCol->begin(currBx); tau != tauDataCol->end(currBx); ++tau)
      mptaus_data->push_back(0, (*tau));
    for (auto tau = tauEmulCol->begin(currBx); tau != tauEmulCol->end(currBx); ++tau)
      mptaus_emul->push_back(0, (*tau));

    e.put(std::move(mptaus_data), "dataTau");
    e.put(std::move(mptaus_emul), "emulTau");

    // store all sums to an edm file
    std::unique_ptr<l1t::EtSumBxCollection> mpsums_data (new l1t::EtSumBxCollection(0, currBx, currBx));
    std::unique_ptr<l1t::EtSumBxCollection> mpsums_emul (new l1t::EtSumBxCollection(0, currBx, currBx));

    for (auto sum = sumDataCol->begin(currBx); sum != sumDataCol->end(currBx); ++sum)
      mpsums_data->push_back(0, (*sum));
    for (auto sum = sumEmulCol->begin(currBx); sum != sumEmulCol->end(currBx); ++sum)
      mpsums_emul->push_back(0, (*sum));

    e.put(std::move(mpsums_data), "dataEtSum");
    e.put(std::move(mpsums_emul), "emulEtSum");

    // store calorimeter towers
    std::unique_ptr<l1t::CaloTowerBxCollection> mptowers_data (new l1t::CaloTowerBxCollection(0, currBx, currBx));
    std::unique_ptr<l1t::CaloTowerBxCollection> mptowers_emul (new l1t::CaloTowerBxCollection(0, currBx, currBx));

    for (auto tower = caloTowerDataCol->begin(currBx); tower != caloTowerDataCol->end(currBx); ++tower)
      mptowers_data->push_back(0, (*tower));
    for (auto tower = caloTowerEmulCol->begin(currBx); tower != caloTowerEmulCol->end(currBx); ++tower)
      mptowers_emul->push_back(0, (*tower));

    e.put(std::move(mptowers_data), "dataCaloTower");
    e.put(std::move(mptowers_emul), "emulCaloTower");
}

DEFINE_FWK_MODULE (L1TStage2CaloLayer2Comp);

#endif