L1RCT.cc

Go to the documentation of this file.

#include "L1Trigger/RegionalCaloTrigger/interface/L1RCT.h"
 
#include <vector>
using std::vector;
 
#include <fstream>
#include <string>
 
#include <iostream>
using std::cerr;
using std::cout;
using std::endl;
using std::ostream;
 
#include <iomanip>
 
//#include "DataFormats/L1CaloTrigger/interface/L1CaloEmCand.h"
 
#include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
#include "CondFormats/L1TObjects/interface/L1RCTParameters.h"
#include "L1Trigger/RegionalCaloTrigger/interface/L1RCTLookupTables.h"
 
// Main method to process a single event, hence the name.
// First it sets up all the neighbors, sharing the pointers to the proper
// regions.  This is done via the neighborMap auxiliary class, which
// is very dry and contains the proper mappings from crate,card,and
// region numbers to the crate,card, and region numbers of the neighbors.
// The next step is to pass along the pointers for the regions
// to their corresponding Electron Isolation Card
// This is done in the crate method fillElectronIsolationCards
// Then the actual processing of the data begins with the
// processReceiverCards and processElectronIsolationCards methods.
// Next the region sums, tau bits, mip bits, and electron
// candidates are passed onto the Jet Summary Card, and that's where
// the data flow ends for the Regional CaloTrigger.
void L1RCT::processEvent() {
  for (int i = 0; i < 18; i++)
    crates.at(i).processReceiverCards();
  shareNeighbors();
  for (int i = 0; i < 18; i++) {
    crates.at(i).fillElectronIsolationCards();
    crates.at(i).processElectronIsolationCards();
    crates.at(i).fillJetSummaryCard();
    crates.at(i).processJetSummaryCard();
  }
}
 
void L1RCT::makeCrates() {
  for (int i = 0; i < 18; i++) {
    L1RCTCrate c(i, rctLookupTables_);
    crates.push_back(c);
  }
}
 
L1RCT::L1RCT(const L1RCTLookupTables *rctLookupTables)
    : rctLookupTables_(rctLookupTables),
      empty(),
      neighborMap(),
      barrel(18, std::vector<std::vector<unsigned short>>(7, std::vector<unsigned short>(64))),
      hf(18, std::vector<unsigned short>(8)) {
  makeCrates();
}
 
void L1RCT::input() {
  for (int i = 0; i < 18; i++) {
    crates.at(i).input(barrel.at(i), hf.at(i));
  }
}
 
void L1RCT::input(const std::vector<std::vector<std::vector<unsigned short>>> &barrelIn,
                  const std::vector<std::vector<unsigned short>> &hfIn) {
  for (int i = 0; i < 18; i++) {
    crates.at(i).input(barrelIn.at(i), hfIn.at(i));
  }
}
 
// This is a method for taking input from a file.  Any entries in excess
// of 18*7*64 will simply be ignored.  This *only* fills input for a single
// event.  At the moment you cannot put a ton of data and have it be
// read as separate events.
void L1RCT::fileInput(const char *filename) {  // added "const" also in .h
  unsigned short x;
  std::ifstream instream(filename);
  if (instream) {
    for (int i = 0; i < 18; i++) {
      for (int j = 0; j < 7; j++) {
        for (int k = 0; k < 64; k++) {
          if (instream >> x) {
            unsigned short bit = x / 256;             // added J.Leonard Aug. 16 06
            unsigned short energy = x & 255;          //
            unsigned short input = energy * 2 + bit;  //
            barrel.at(i).at(j).at(k) = input;
          } else
            break;
        }
      }
      for (int j = 0; j < 8; j++) {
        if (instream >> x) {
          hf.at(i).at(j) = x;
        } else
          break;
      }
    }
  }
  input();
}
 
// takes hcal and ecal digi input, including HF
void L1RCT::digiInput(const EcalTrigPrimDigiCollection &ecalCollection,
                      const HcalTrigPrimDigiCollection &hcalCollection) {
  // fills input vectors with 0's in case ecal or hcal collection not used
  for (int i = 0; i < 18; i++) {
    for (int j = 0; j < 7; j++) {
      for (int k = 0; k < 64; k++) {
        barrel.at(i).at(j).at(k) = 0;
      }
    }
    for (int j = 0; j < 8; j++) {
      hf.at(i).at(j) = 0;
    }
  }
 
  int nEcalDigi = ecalCollection.size();
  if (nEcalDigi > 4032) {
    nEcalDigi = 4032;
  }
  for (int i = 0; i < nEcalDigi; i++) {
    short ieta = (short)ecalCollection[i].id().ieta();
    // Note absIeta counts from 1-28 (not 0-27)
    unsigned short absIeta = (unsigned short)abs(ieta);
    unsigned short cal_iphi = (unsigned short)ecalCollection[i].id().iphi();
    unsigned short iphi = (72 + 18 - cal_iphi) % 72;  // transform TOWERS (not regions) into local
                                                      // rct (intuitive) phi bins
 
    // map digis to crates, cards, and towers
    unsigned short crate = 999, card = 999, tower = 999;
    crate = rctLookupTables_->rctParameters()->calcCrate(iphi, ieta);
    card = rctLookupTables_->rctParameters()->calcCard(iphi, absIeta);
    tower = rctLookupTables_->rctParameters()->calcTower(iphi, absIeta);
 
    unsigned short energy = ecalCollection[i].compressedEt();
    unsigned short fineGrain = (unsigned short)ecalCollection[i].fineGrain();  // 0 or 1
    unsigned short ecalInput = energy * 2 + fineGrain;
 
    // put input into correct crate/card/tower of barrel
    if ((crate < 18) && (card < 7) && (tower < 32)) {   // changed 64 to 32 Sept. 19 J. Leonard, rm -1 7Nov07
      barrel.at(crate).at(card).at(tower) = ecalInput;  // rm -1
    } else {
      std::cerr << "L1RCT: ecal out of range! tower = " << tower << " iphi is " << iphi << " absieta is " << absIeta
                << std::endl;
    }
  }
 
  // same for hcal, once we get the hcal digis, just need to add 32 to towers:
  // just copied and pasted and changed names where necessary
  int nHcalDigi = hcalCollection.size();
  // if (nHcalDigi != 4176){ std::cout << "L1RCT: Warning: There are " <<
  // nHcalDigi << "hcal digis instead of 4176!" << std::endl;}
  // incl HF 4032 + 144 = 4176
  for (int i = 0; i < nHcalDigi; i++) {
    if (hcalCollection[i].id().version() != 0) {
      continue;
    }
    short ieta = (short)hcalCollection[i].id().ieta();
    unsigned short absIeta = (unsigned short)abs(ieta);
    unsigned short cal_iphi = (unsigned short)hcalCollection[i].id().iphi();
    // All Hcal primitives (including HF) are reported
    // with phi bin numbering in the range 0-72.
    unsigned short iphi = (72 + 18 - cal_iphi) % 72;
    // transform Hcal TOWERS (1-72)into local rct (intuitive) phi bins (72 bins)
    // 0-71 Use local iphi to work out the region and crate (for HB/HE and HF)
    // HF regions need to have local iphi 0-17
    if (absIeta >= 29) {
      iphi = iphi / 4;
    }
 
    // map digis to crates, cards, and towers
    unsigned short crate = 999, card = 999, tower = 999;
    crate = rctLookupTables_->rctParameters()->calcCrate(iphi, ieta);
    if (absIeta < 29) {
      card = rctLookupTables_->rctParameters()->calcCard(iphi, absIeta);
    }
    tower = rctLookupTables_->rctParameters()->calcTower(iphi, absIeta);
 
    unsigned short energy = hcalCollection[i].SOI_compressedEt();  // access only sample of interest
    unsigned short fineGrain = (unsigned short)hcalCollection[i].SOI_fineGrain();
    unsigned short hcalInput = energy * 2 + fineGrain;
    if (absIeta <= 28) {
      // put input into correct crate/card/tower of barrel
      if ((crate < 18) && (card < 7) && (tower < 32)) {        // changed 64 to 32 Sept. 19 J. Leonard, rm -1 7Nov07
        barrel.at(crate).at(card).at(tower + 32) = hcalInput;  // hcal towers are ecal + 32 see RC.cc; rm -1 7Nov07
      } else {
        std::cout << "L1RCT: hcal out of range!  tower = " << tower << std::endl;
      }
    } else if ((absIeta >= 29) && (absIeta <= 32)) {
      // put input into correct crate/region of HF
      if ((crate < 18) && (tower < 8)) {
        hf.at(crate).at(tower) = hcalInput;
      } else {
        std::cout << "L1RCT: hf out of range!  region = " << tower << std::endl;
      }
    }
  }
 
  input();
 
  return;
}
 
// As the name implies, it will randomly generate input for the
// regional calotrigger.
void L1RCT::randomInput() {
  for (int i = 0; i < 18; i++) {
    for (int j = 0; j < 7; j++) {
      for (int k = 0; k < 64; k++) {
        barrel.at(i).at(j).at(k) = rand() % 511;
      }
    }
    for (int j = 0; j < 8; j++) {
      hf.at(i).at(j) = rand() % 255;  // changed from 1023 (10 bits)
    }
  }
  input();
  return;
}
 
// This method handles the bulk of the pointer passing, giving
// to each region pointers to its neighbors.  If it does *not*
// have a neighbor in that direction then it passes it a pointer
// to an empty region that contains no data and is disconnected
// from anything else.  This makes the electron finding algorithm simpler
// as then all regions can be treated equally.
void L1RCT::shareNeighbors() {
  L1RCTRegion *north;
  L1RCTRegion *south;
  L1RCTRegion *west;
  L1RCTRegion *east;
  L1RCTRegion *se;
  L1RCTRegion *sw;
  L1RCTRegion *nw;
  L1RCTRegion *ne;
  L1RCTRegion *primary;
 
  for (int i = 0; i < 18; i++) {
    for (int j = 0; j < 7; j++) {
      for (int k = 0; k < 2; k++) {
        primary = crates.at(i).getReceiverCard(j)->getRegion(k);
 
        vector<int> northIndices = neighborMap.north(i, j, k);
        if (northIndices.at(0) != -1)
          north = crates.at(northIndices.at(0)).getReceiverCard(northIndices.at(1))->getRegion(northIndices.at(2));
        else
          north = &empty;
 
        vector<int> southIndices = neighborMap.south(i, j, k);
        if (southIndices.at(0) != -1)
          south = crates.at(southIndices.at(0)).getReceiverCard(southIndices.at(1))->getRegion(southIndices.at(2));
        else
          south = &empty;
 
        vector<int> westIndices = neighborMap.west(i, j, k);
        if (westIndices.at(0) != -1)
          west = crates.at(westIndices.at(0)).getReceiverCard(westIndices.at(1))->getRegion(westIndices.at(2));
        else
          west = &empty;
 
        vector<int> eastIndices = neighborMap.east(i, j, k);
        if (eastIndices.at(0) != -1)
          east = crates.at(eastIndices.at(0)).getReceiverCard(eastIndices.at(1))->getRegion(eastIndices.at(2));
        else
          east = &empty;
 
        vector<int> seIndices = neighborMap.se(i, j, k);
        if (seIndices.at(0) != -1)
          se = crates.at(seIndices.at(0)).getReceiverCard(seIndices.at(1))->getRegion(seIndices.at(2));
        else
          se = &empty;
 
        vector<int> swIndices = neighborMap.sw(i, j, k);
        if (swIndices.at(0) != -1)
          sw = crates.at(swIndices.at(0)).getReceiverCard(swIndices.at(1))->getRegion(swIndices.at(2));
        else
          sw = &empty;
 
        vector<int> neIndices = neighborMap.ne(i, j, k);
        if (neIndices.at(0) != -1)
          ne = crates.at(neIndices.at(0)).getReceiverCard(neIndices.at(1))->getRegion(neIndices.at(2));
        else
          ne = &empty;
 
        vector<int> nwIndices = neighborMap.nw(i, j, k);
        if (nwIndices.at(0) != -1)
          nw = crates.at(nwIndices.at(0)).getReceiverCard(nwIndices.at(1))->getRegion(nwIndices.at(2));
        else
          nw = &empty;
 
        primary->setNorthEt(north->giveNorthEt());
        primary->setNorthHE_FG(north->giveNorthHE_FG());
        primary->setSouthEt(south->giveSouthEt());
        primary->setSouthHE_FG(south->giveSouthHE_FG());
        primary->setEastEt(east->giveEastEt());
        primary->setEastHE_FG(east->giveEastHE_FG());
        primary->setWestEt(west->giveWestEt());
        primary->setWestHE_FG(west->giveWestHE_FG());
        primary->setSEEt(se->giveSEEt());
        primary->setSEHE_FG(se->giveSEHE_FG());
        primary->setSWEt(sw->giveSWEt());
        primary->setSWHE_FG(sw->giveSWHE_FG());
        primary->setNWEt(nw->giveNWEt());
        primary->setNWHE_FG(nw->giveNWHE_FG());
        primary->setNEEt(ne->giveNEEt());
        primary->setNEHE_FG(ne->giveNEHE_FG());
      }
    }
  }
}
 
void L1RCT::print() {
  for (int i = 0; i < 18; i++) {
    std::cout << "Crate " << i << std::endl;
    crates.at(i).print();
  }
}
 
// Returns the top four isolated electrons from given crate
// in a vector of L1CaloEmCands
L1CaloEmCollection L1RCT::getIsolatedEGObjects(unsigned crate) {
  std::vector<unsigned short> isoEmObjects = crates.at(crate).getIsolatedEGObjects();
  L1CaloEmCollection isoEmCands;
  for (uint16_t i = 0; i < 4; i++) {
    unsigned rgn = ((isoEmObjects.at(i)) & 1);
    unsigned crd = (((isoEmObjects.at(i)) / 2) & 7);
    unsigned energy = ((isoEmObjects.at(i)) / 16);
    unsigned rank = rctLookupTables_->emRank(energy);
    L1CaloEmCand isoCand(rank, rgn, crd, crate, true, i,
                         0);  // includes emcand index
    isoEmCands.push_back(isoCand);
  }
  return isoEmCands;
}
 
// Returns the top four nonisolated electrons from the given crate
// in a vector of L1CaloEmCands
L1CaloEmCollection L1RCT::getNonisolatedEGObjects(unsigned crate) {
  std::vector<unsigned short> nonIsoEmObjects = crates.at(crate).getNonisolatedEGObjects();
  L1CaloEmCollection nonIsoEmCands;
  for (uint16_t i = 0; i < 4; i++) {
    unsigned rgn = ((nonIsoEmObjects.at(i)) & 1);
    unsigned crd = (((nonIsoEmObjects.at(i)) / 2) & 7);
    unsigned energy = ((nonIsoEmObjects.at(i)) / 16);
    unsigned rank = rctLookupTables_->emRank(energy);
    L1CaloEmCand nonIsoCand(rank, rgn, crd, crate, false, i,
                            0);  // includes emcand index
    nonIsoEmCands.push_back(nonIsoCand);
  }
  return nonIsoEmCands;
}
 
vector<L1CaloRegion> L1RCT::getRegions(unsigned crate) {
  // barrel regions
  std::bitset<14> taus((long)crates.at(crate).getTauBits());
  std::bitset<14> mips((long)crates.at(crate).getMIPBits());
  std::bitset<14> quiets((long)crates.at(crate).getQuietBits());
  std::bitset<14> overflows((long)crates.at(crate).getOverFlowBits());
  std::vector<unsigned short> barrelEnergies = crates.at(crate).getBarrelRegions();
  std::vector<L1CaloRegion> regionCollection;
  for (unsigned card = 0; card < 7; card++) {
    for (unsigned rgn = 0; rgn < 2; rgn++) {
      bool tau = taus[card * 2 + rgn];
      bool mip = mips[card * 2 + rgn];
      bool quiet = quiets[card * 2 + rgn];
      bool overflow = overflows[card * 2 + rgn];
      unsigned barrelEnergy = barrelEnergies.at(card * 2 + rgn);
      L1CaloRegion region(barrelEnergy, overflow, tau, mip, quiet, crate, card, rgn);
      regionCollection.push_back(region);
    }
  }
 
  // hf regions
  std::vector<unsigned short> hfEnergies = crates.at(crate).getHFRegions();
  // fine grain bits -- still have to work out digi input
  std::vector<unsigned short> hfFineGrainBits = crates.at(crate).getHFFineGrainBits();
  for (unsigned hfRgn = 0; hfRgn < 8; hfRgn++) {  // region number, see diagram on paper.  make sure know how hf
                                                  // regions come in.
    unsigned energy = hfEnergies.at(hfRgn);
    bool fineGrain = hfFineGrainBits.at(hfRgn);
    L1CaloRegion hfRegion(energy, fineGrain, crate, hfRgn);
    regionCollection.push_back(hfRegion);
  }
  return regionCollection;
}