RegionalCaloTrigger/src/L1RCT.cc

0001 #include "L1Trigger/RegionalCaloTrigger/interface/L1RCT.h"
0002
0003 #include <vector>
0004 using std::vector;
0005
0006 #include <fstream>
0007 #include <string>
0008
0009 #include <iostream>
0010 using std::cerr;
0011 using std::cout;
0012 using std::endl;
0013 using std::ostream;
0014
0015 #include <iomanip>
0016
0017 //#include "DataFormats/L1CaloTrigger/interface/L1CaloEmCand.h"
0018
0019 #include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
0020 #include "CondFormats/L1TObjects/interface/L1RCTParameters.h"
0021 #include "L1Trigger/RegionalCaloTrigger/interface/L1RCTLookupTables.h"
0022
0023 // Main method to process a single event, hence the name.
0024 // First it sets up all the neighbors, sharing the pointers to the proper
0025 // regions.  This is done via the neighborMap auxiliary class, which
0026 // is very dry and contains the proper mappings from crate,card,and
0027 // region numbers to the crate,card, and region numbers of the neighbors.
0028 // The next step is to pass along the pointers for the regions
0029 // to their corresponding Electron Isolation Card
0030 // This is done in the crate method fillElectronIsolationCards
0031 // Then the actual processing of the data begins with the
0032 // processReceiverCards and processElectronIsolationCards methods.
0033 // Next the region sums, tau bits, mip bits, and electron
0034 // candidates are passed onto the Jet Summary Card, and that's where
0035 // the data flow ends for the Regional CaloTrigger.
0036 void L1RCT::processEvent() {
0037   for (int i = 0; i < 18; i++)
0038     crates.at(i).processReceiverCards();
0039   shareNeighbors();
0040   for (int i = 0; i < 18; i++) {
0041     crates.at(i).fillElectronIsolationCards();
0042     crates.at(i).processElectronIsolationCards();
0043     crates.at(i).fillJetSummaryCard();
0044     crates.at(i).processJetSummaryCard();
0045   }
0046 }
0047
0048 void L1RCT::makeCrates() {
0049   for (int i = 0; i < 18; i++) {
0050     L1RCTCrate c(i, rctLookupTables_);
0051     crates.push_back(c);
0052   }
0053 }
0054
0055 L1RCT::L1RCT(const L1RCTLookupTables *rctLookupTables)
0056     : rctLookupTables_(rctLookupTables),
0057       empty(),
0058       neighborMap(),
0059       barrel(18, std::vector<std::vector<unsigned short>>(7, std::vector<unsigned short>(64))),
0060       hf(18, std::vector<unsigned short>(8)) {
0061   makeCrates();
0062 }
0063
0064 void L1RCT::input() {
0065   for (int i = 0; i < 18; i++) {
0066     crates.at(i).input(barrel.at(i), hf.at(i));
0067   }
0068 }
0069
0070 void L1RCT::input(const std::vector<std::vector<std::vector<unsigned short>>> &barrelIn,
0071                   const std::vector<std::vector<unsigned short>> &hfIn) {
0072   for (int i = 0; i < 18; i++) {
0073     crates.at(i).input(barrelIn.at(i), hfIn.at(i));
0074   }
0075 }
0076
0077 // This is a method for taking input from a file.  Any entries in excess
0078 // of 18*7*64 will simply be ignored.  This *only* fills input for a single
0079 // event.  At the moment you cannot put a ton of data and have it be
0080 // read as separate events.
0081 void L1RCT::fileInput(const char *filename) {  // added "const" also in .h
0082   unsigned short x;
0083   std::ifstream instream(filename);
0084   if (instream) {
0085     for (int i = 0; i < 18; i++) {
0086       for (int j = 0; j < 7; j++) {
0087         for (int k = 0; k < 64; k++) {
0088           if (instream >> x) {
0089             unsigned short bit = x / 256;             // added J.Leonard Aug. 16 06
0090             unsigned short energy = x & 255;          //
0091             unsigned short input = energy * 2 + bit;  //
0092             barrel.at(i).at(j).at(k) = input;
0093           } else
0094             break;
0095         }
0096       }
0097       for (int j = 0; j < 8; j++) {
0098         if (instream >> x) {
0099           hf.at(i).at(j) = x;
0100         } else
0101           break;
0102       }
0103     }
0104   }
0105   input();
0106 }
0107
0108 // takes hcal and ecal digi input, including HF
0109 void L1RCT::digiInput(const EcalTrigPrimDigiCollection &ecalCollection,
0110                       const HcalTrigPrimDigiCollection &hcalCollection) {
0111   // fills input vectors with 0's in case ecal or hcal collection not used
0112   for (int i = 0; i < 18; i++) {
0113     for (int j = 0; j < 7; j++) {
0114       for (int k = 0; k < 64; k++) {
0115         barrel.at(i).at(j).at(k) = 0;
0116       }
0117     }
0118     for (int j = 0; j < 8; j++) {
0119       hf.at(i).at(j) = 0;
0120     }
0121   }
0122
0123   int nEcalDigi = ecalCollection.size();
0124   if (nEcalDigi > 4032) {
0125     nEcalDigi = 4032;
0126   }
0127   for (int i = 0; i < nEcalDigi; i++) {
0128     short ieta = (short)ecalCollection[i].id().ieta();
0129     // Note absIeta counts from 1-28 (not 0-27)
0130     unsigned short absIeta = (unsigned short)abs(ieta);
0131     unsigned short cal_iphi = (unsigned short)ecalCollection[i].id().iphi();
0132     unsigned short iphi = (72 + 18 - cal_iphi) % 72;  // transform TOWERS (not regions) into local
0133                                                       // rct (intuitive) phi bins
0134
0135     // map digis to crates, cards, and towers
0136     unsigned short crate = 999, card = 999, tower = 999;
0137     crate = rctLookupTables_->rctParameters()->calcCrate(iphi, ieta);
0138     card = rctLookupTables_->rctParameters()->calcCard(iphi, absIeta);
0139     tower = rctLookupTables_->rctParameters()->calcTower(iphi, absIeta);
0140
0141     unsigned short energy = ecalCollection[i].compressedEt();
0142     unsigned short fineGrain = (unsigned short)ecalCollection[i].fineGrain();  // 0 or 1
0143     unsigned short ecalInput = energy * 2 + fineGrain;
0144
0145     // put input into correct crate/card/tower of barrel
0146     if ((crate < 18) && (card < 7) && (tower < 32)) {   // changed 64 to 32 Sept. 19 J. Leonard, rm -1 7Nov07
0147       barrel.at(crate).at(card).at(tower) = ecalInput;  // rm -1
0148     } else {
0149       std::cerr << "L1RCT: ecal out of range! tower = " << tower << " iphi is " << iphi << " absieta is " << absIeta
0150                 << std::endl;
0151     }
0152   }
0153
0154   // same for hcal, once we get the hcal digis, just need to add 32 to towers:
0155   // just copied and pasted and changed names where necessary
0156   int nHcalDigi = hcalCollection.size();
0157   // if (nHcalDigi != 4176){ std::cout << "L1RCT: Warning: There are " <<
0158   // nHcalDigi << "hcal digis instead of 4176!" << std::endl;}
0159   // incl HF 4032 + 144 = 4176
0160   for (int i = 0; i < nHcalDigi; i++) {
0161     if (hcalCollection[i].id().version() != 0) {
0162       continue;
0163     }
0164     short ieta = (short)hcalCollection[i].id().ieta();
0165     unsigned short absIeta = (unsigned short)abs(ieta);
0166     unsigned short cal_iphi = (unsigned short)hcalCollection[i].id().iphi();
0167     // All Hcal primitives (including HF) are reported
0168     // with phi bin numbering in the range 0-72.
0169     unsigned short iphi = (72 + 18 - cal_iphi) % 72;
0170     // transform Hcal TOWERS (1-72)into local rct (intuitive) phi bins (72 bins)
0171     // 0-71 Use local iphi to work out the region and crate (for HB/HE and HF)
0172     // HF regions need to have local iphi 0-17
0173     if (absIeta >= 29) {
0174       iphi = iphi / 4;
0175     }
0176
0177     // map digis to crates, cards, and towers
0178     unsigned short crate = 999, card = 999, tower = 999;
0179     crate = rctLookupTables_->rctParameters()->calcCrate(iphi, ieta);
0180     if (absIeta < 29) {
0181       card = rctLookupTables_->rctParameters()->calcCard(iphi, absIeta);
0182     }
0183     tower = rctLookupTables_->rctParameters()->calcTower(iphi, absIeta);
0184
0185     unsigned short energy = hcalCollection[i].SOI_compressedEt();  // access only sample of interest
0186     unsigned short fineGrain = (unsigned short)hcalCollection[i].SOI_fineGrain();
0187     unsigned short hcalInput = energy * 2 + fineGrain;
0188     if (absIeta <= 28) {
0189       // put input into correct crate/card/tower of barrel
0190       if ((crate < 18) && (card < 7) && (tower < 32)) {        // changed 64 to 32 Sept. 19 J. Leonard, rm -1 7Nov07
0191         barrel.at(crate).at(card).at(tower + 32) = hcalInput;  // hcal towers are ecal + 32 see RC.cc; rm -1 7Nov07
0192       } else {
0193         std::cout << "L1RCT: hcal out of range!  tower = " << tower << std::endl;
0194       }
0195     } else if ((absIeta >= 29) && (absIeta <= 32)) {
0196       // put input into correct crate/region of HF
0197       if ((crate < 18) && (tower < 8)) {
0198         hf.at(crate).at(tower) = hcalInput;
0199       } else {
0200         std::cout << "L1RCT: hf out of range!  region = " << tower << std::endl;
0201       }
0202     }
0203   }
0204
0205   input();
0206
0207   return;
0208 }
0209
0210 // As the name implies, it will randomly generate input for the
0211 // regional calotrigger.
0212 void L1RCT::randomInput() {
0213   for (int i = 0; i < 18; i++) {
0214     for (int j = 0; j < 7; j++) {
0215       for (int k = 0; k < 64; k++) {
0216         barrel.at(i).at(j).at(k) = rand() % 511;
0217       }
0218     }
0219     for (int j = 0; j < 8; j++) {
0220       hf.at(i).at(j) = rand() % 255;  // changed from 1023 (10 bits)
0221     }
0222   }
0223   input();
0224   return;
0225 }
0226
0227 // This method handles the bulk of the pointer passing, giving
0228 // to each region pointers to its neighbors.  If it does *not*
0229 // have a neighbor in that direction then it passes it a pointer
0230 // to an empty region that contains no data and is disconnected
0231 // from anything else.  This makes the electron finding algorithm simpler
0232 // as then all regions can be treated equally.
0233 void L1RCT::shareNeighbors() {
0234   L1RCTRegion *north;
0235   L1RCTRegion *south;
0236   L1RCTRegion *west;
0237   L1RCTRegion *east;
0238   L1RCTRegion *se;
0239   L1RCTRegion *sw;
0240   L1RCTRegion *nw;
0241   L1RCTRegion *ne;
0242   L1RCTRegion *primary;
0243
0244   for (int i = 0; i < 18; i++) {
0245     for (int j = 0; j < 7; j++) {
0246       for (int k = 0; k < 2; k++) {
0247         primary = crates.at(i).getReceiverCard(j)->getRegion(k);
0248
0249         vector<int> northIndices = neighborMap.north(i, j, k);
0250         if (northIndices.at(0) != -1)
0251           north = crates.at(northIndices.at(0)).getReceiverCard(northIndices.at(1))->getRegion(northIndices.at(2));
0252         else
0253           north = &empty;
0254
0255         vector<int> southIndices = neighborMap.south(i, j, k);
0256         if (southIndices.at(0) != -1)
0257           south = crates.at(southIndices.at(0)).getReceiverCard(southIndices.at(1))->getRegion(southIndices.at(2));
0258         else
0259           south = &empty;
0260
0261         vector<int> westIndices = neighborMap.west(i, j, k);
0262         if (westIndices.at(0) != -1)
0263           west = crates.at(westIndices.at(0)).getReceiverCard(westIndices.at(1))->getRegion(westIndices.at(2));
0264         else
0265           west = &empty;
0266
0267         vector<int> eastIndices = neighborMap.east(i, j, k);
0268         if (eastIndices.at(0) != -1)
0269           east = crates.at(eastIndices.at(0)).getReceiverCard(eastIndices.at(1))->getRegion(eastIndices.at(2));
0270         else
0271           east = &empty;
0272
0273         vector<int> seIndices = neighborMap.se(i, j, k);
0274         if (seIndices.at(0) != -1)
0275           se = crates.at(seIndices.at(0)).getReceiverCard(seIndices.at(1))->getRegion(seIndices.at(2));
0276         else
0277           se = &empty;
0278
0279         vector<int> swIndices = neighborMap.sw(i, j, k);
0280         if (swIndices.at(0) != -1)
0281           sw = crates.at(swIndices.at(0)).getReceiverCard(swIndices.at(1))->getRegion(swIndices.at(2));
0282         else
0283           sw = &empty;
0284
0285         vector<int> neIndices = neighborMap.ne(i, j, k);
0286         if (neIndices.at(0) != -1)
0287           ne = crates.at(neIndices.at(0)).getReceiverCard(neIndices.at(1))->getRegion(neIndices.at(2));
0288         else
0289           ne = &empty;
0290
0291         vector<int> nwIndices = neighborMap.nw(i, j, k);
0292         if (nwIndices.at(0) != -1)
0293           nw = crates.at(nwIndices.at(0)).getReceiverCard(nwIndices.at(1))->getRegion(nwIndices.at(2));
0294         else
0295           nw = &empty;
0296
0297         primary->setNorthEt(north->giveNorthEt());
0298         primary->setNorthHE_FG(north->giveNorthHE_FG());
0299         primary->setSouthEt(south->giveSouthEt());
0300         primary->setSouthHE_FG(south->giveSouthHE_FG());
0301         primary->setEastEt(east->giveEastEt());
0302         primary->setEastHE_FG(east->giveEastHE_FG());
0303         primary->setWestEt(west->giveWestEt());
0304         primary->setWestHE_FG(west->giveWestHE_FG());
0305         primary->setSEEt(se->giveSEEt());
0306         primary->setSEHE_FG(se->giveSEHE_FG());
0307         primary->setSWEt(sw->giveSWEt());
0308         primary->setSWHE_FG(sw->giveSWHE_FG());
0309         primary->setNWEt(nw->giveNWEt());
0310         primary->setNWHE_FG(nw->giveNWHE_FG());
0311         primary->setNEEt(ne->giveNEEt());
0312         primary->setNEHE_FG(ne->giveNEHE_FG());
0313       }
0314     }
0315   }
0316 }
0317
0318 void L1RCT::print() {
0319   for (int i = 0; i < 18; i++) {
0320     std::cout << "Crate " << i << std::endl;
0321     crates.at(i).print();
0322   }
0323 }
0324
0325 // Returns the top four isolated electrons from given crate
0326 // in a vector of L1CaloEmCands
0327 L1CaloEmCollection L1RCT::getIsolatedEGObjects(unsigned crate) {
0328   std::vector<unsigned short> isoEmObjects = crates.at(crate).getIsolatedEGObjects();
0329   L1CaloEmCollection isoEmCands;
0330   for (uint16_t i = 0; i < 4; i++) {
0331     unsigned rgn = ((isoEmObjects.at(i)) & 1);
0332     unsigned crd = (((isoEmObjects.at(i)) / 2) & 7);
0333     unsigned energy = ((isoEmObjects.at(i)) / 16);
0334     unsigned rank = rctLookupTables_->emRank(energy);
0335     L1CaloEmCand isoCand(rank, rgn, crd, crate, true, i,
0336                          0);  // includes emcand index
0337     isoEmCands.push_back(isoCand);
0338   }
0339   return isoEmCands;
0340 }
0341
0342 // Returns the top four nonisolated electrons from the given crate
0343 // in a vector of L1CaloEmCands
0344 L1CaloEmCollection L1RCT::getNonisolatedEGObjects(unsigned crate) {
0345   std::vector<unsigned short> nonIsoEmObjects = crates.at(crate).getNonisolatedEGObjects();
0346   L1CaloEmCollection nonIsoEmCands;
0347   for (uint16_t i = 0; i < 4; i++) {
0348     unsigned rgn = ((nonIsoEmObjects.at(i)) & 1);
0349     unsigned crd = (((nonIsoEmObjects.at(i)) / 2) & 7);
0350     unsigned energy = ((nonIsoEmObjects.at(i)) / 16);
0351     unsigned rank = rctLookupTables_->emRank(energy);
0352     L1CaloEmCand nonIsoCand(rank, rgn, crd, crate, false, i,
0353                             0);  // includes emcand index
0354     nonIsoEmCands.push_back(nonIsoCand);
0355   }
0356   return nonIsoEmCands;
0357 }
0358
0359 vector<L1CaloRegion> L1RCT::getRegions(unsigned crate) {
0360   // barrel regions
0361   std::bitset<14> taus((long)crates.at(crate).getTauBits());
0362   std::bitset<14> mips((long)crates.at(crate).getMIPBits());
0363   std::bitset<14> quiets((long)crates.at(crate).getQuietBits());
0364   std::bitset<14> overflows((long)crates.at(crate).getOverFlowBits());
0365   std::vector<unsigned short> barrelEnergies = crates.at(crate).getBarrelRegions();
0366   std::vector<L1CaloRegion> regionCollection;
0367   for (unsigned card = 0; card < 7; card++) {
0368     for (unsigned rgn = 0; rgn < 2; rgn++) {
0369       bool tau = taus[card * 2 + rgn];
0370       bool mip = mips[card * 2 + rgn];
0371       bool quiet = quiets[card * 2 + rgn];
0372       bool overflow = overflows[card * 2 + rgn];
0373       unsigned barrelEnergy = barrelEnergies.at(card * 2 + rgn);
0374       L1CaloRegion region(barrelEnergy, overflow, tau, mip, quiet, crate, card, rgn);
0375       regionCollection.push_back(region);
0376     }
0377   }
0378
0379   // hf regions
0380   std::vector<unsigned short> hfEnergies = crates.at(crate).getHFRegions();
0381   // fine grain bits -- still have to work out digi input
0382   std::vector<unsigned short> hfFineGrainBits = crates.at(crate).getHFFineGrainBits();
0383   for (unsigned hfRgn = 0; hfRgn < 8; hfRgn++) {  // region number, see diagram on paper.  make sure know how hf
0384                                                   // regions come in.
0385     unsigned energy = hfEnergies.at(hfRgn);
0386     bool fineGrain = hfFineGrainBits.at(hfRgn);
0387     L1CaloRegion hfRegion(energy, fineGrain, crate, hfRgn);
0388     regionCollection.push_back(hfRegion);
0389   }
0390   return regionCollection;
0391 }