Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:19:52 

0001 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctHardwareJetFinder.h"
0002 
0003 //DEFINE STATICS
0004 const unsigned int L1GctHardwareJetFinder::N_COLS = 2;
0005 const unsigned int L1GctHardwareJetFinder::CENTRAL_COL0 = 0;
0006 const unsigned int L1GctHardwareJetFinder::MAX_REGIONS_IN =
0007     (((L1CaloRegionDetId::N_ETA) / 2) + N_EXTRA_REGIONS_ETA00) * L1GctHardwareJetFinder::N_COLS;
0008 
0009 L1GctHardwareJetFinder::L1GctHardwareJetFinder(int id)
0010     : L1GctJetFinderBase(id),
0011       m_localMaxima(MAX_JETS_OUT),
0012       m_clusters(MAX_JETS_OUT),
0013       m_numberOfClusters(0),
0014       m_localMax00(2),
0015       m_cluster00(2) {
0016   this->reset();
0017   // Initialise parameters for Region input calculations in the
0018   // derived class so we get the right values of constants.
0019   static const unsigned NPHI = L1CaloRegionDetId::N_PHI;
0020   m_minColThisJf = (NPHI + m_id * 2 - CENTRAL_COL0) % NPHI;
0021 }
0022 
0023 L1GctHardwareJetFinder::~L1GctHardwareJetFinder() {}
0024 
0025 std::ostream& operator<<(std::ostream& os, const L1GctHardwareJetFinder& algo) {
0026   os << "===L1GctHardwareJetFinder===" << std::endl;
0027   const L1GctJetFinderBase* temp = &algo;
0028   os << *temp;
0029   return os;
0030 }
0031 
0032 void L1GctHardwareJetFinder::reset() { L1GctJetFinderBase::reset(); }
0033 
0034 void L1GctHardwareJetFinder::fetchInput() {
0035   if (setupOk()) {
0036     findProtoJets();
0037   }
0038 }
0039 
0040 void L1GctHardwareJetFinder::process() {
0041   if (setupOk()) {
0042     fetchProtoJetsFromNeighbour(TOPBOT);
0043     findJets();
0044     sortJets();
0045     doEnergySums();
0046   }
0047 }
0048 
0049 /// HERE IS THE JETFINDER CODE
0050 
0051 /// The first stage of clustering, called by fetchInput()
0052 void L1GctHardwareJetFinder::findProtoJets() {
0053   findLocalMaxima();
0054   findProtoClusters();
0055   convertClustersToProtoJets();
0056 }
0057 
0058 /// The second stage of clustering, called by process()
0059 void L1GctHardwareJetFinder::findJets() {
0060   findFinalClusters();
0061   convertClustersToOutputJets();
0062 }
0063 
0064 /// Both clustering stages need to find local maxima in the search array
0065 //  Find the local et maxima in the 2x11 array of regions
0066 void L1GctHardwareJetFinder::findLocalMaxima() {
0067   m_localMaxima.clear();
0068   m_localMaxima.resize(MAX_JETS_OUT);
0069   m_localMax00.clear();
0070   m_localMax00.resize(2);
0071 
0072   UShort jetNum = 0;  //holds the number of jets currently found
0073   UShort centreIndex = COL_OFFSET * this->centralCol0();
0074   for (UShort column = 0; column < 2; ++column)  //Find jets in the central search region
0075   {
0076     // The input regions include two extra bins on the other side of eta=0. This allows "seamless"
0077     // jetfinding across the eta=0 boundary. We skip the first input region in each row. We perform
0078     // the full pre-clustering on the next region but store the resulting clusters separately
0079     // from the main list of output pre-clusters - they will be used in the final cluster stage to
0080     // make sure we do not produce jets in adjacent regions on opposite sides of eta=0.
0081     ++centreIndex;
0082     for (UShort row = 1; row < COL_OFFSET; ++row) {
0083       // Here's the array of greater-than and greater-or-equal tests
0084       // to ensure each localMaximum appears once and only once in the list
0085       // It is different for forward and backward eta.
0086       unsigned JET_THRESHOLD = ((row > m_EtaBoundry) ? m_FwdJetSeed : m_CenJetSeed);
0087       bool localMax = !m_inputRegions.at(centreIndex).empty() && (m_inputRegions.at(centreIndex).et() >= JET_THRESHOLD);
0088       if (m_positiveEtaWheel) {  // Forward eta
0089         localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - 1).et());
0090         if (row < (COL_OFFSET - 1)) {
0091           localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + 1).et());
0092         }
0093         if (column == 0) {
0094           localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + COL_OFFSET).et());
0095           localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + COL_OFFSET - 1).et());
0096           if (row < (COL_OFFSET - 1)) {
0097             localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + COL_OFFSET + 1).et());
0098           }
0099         } else {
0100           localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - COL_OFFSET).et());
0101           localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - COL_OFFSET - 1).et());
0102           if (row < (COL_OFFSET - 1)) {
0103             localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - COL_OFFSET + 1).et());
0104           }
0105         }
0106       } else {  // Backward eta
0107         localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - 1).et());
0108         if (row < (COL_OFFSET - 1)) {
0109           localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + 1).et());
0110         }
0111         if (column == 0) {
0112           localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + COL_OFFSET).et());
0113           localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + COL_OFFSET - 1).et());
0114           if (row < (COL_OFFSET - 1)) {
0115             localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + COL_OFFSET + 1).et());
0116           }
0117         } else {
0118           localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - COL_OFFSET).et());
0119           localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - COL_OFFSET - 1).et());
0120           if (row < (COL_OFFSET - 1)) {
0121             localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - COL_OFFSET + 1).et());
0122           }
0123         }
0124       }
0125       if (localMax) {
0126         if (row > 1) {
0127           if (jetNum < MAX_JETS_OUT) {
0128             m_localMaxima.at(jetNum) = m_inputRegions.at(centreIndex);
0129             ++jetNum;
0130           }
0131         }
0132         // Treat row 1 as a separate case. It's not required for jetfinding but
0133         // is used for vetoing of jets double counted across the eta=0 boundary
0134         else {
0135           unsigned phi = m_inputRegions.at(centreIndex).rctPhi();
0136           m_localMax00.at(phi) = m_inputRegions.at(centreIndex);
0137         }
0138       }
0139       ++centreIndex;
0140     }
0141   }
0142 
0143   m_numberOfClusters = jetNum;
0144 }
0145 
0146 //  For each local maximum, find the cluster et in a 2x3 region.
0147 //  The logic ensures that a given region et cannot be used in more than one cluster.
0148 //  The sorting of the local maxima ensures the highest et maximum has priority.
0149 void L1GctHardwareJetFinder::findProtoClusters() {
0150   m_clusters.clear();
0151   m_clusters.resize(MAX_JETS_OUT);
0152   m_cluster00.clear();
0153   m_cluster00.resize(2);
0154 
0155   RegionsVector topJets(MAX_JETS_OUT), botJets(MAX_JETS_OUT);
0156   std::vector<unsigned> topJetsPosition(MAX_JETS_OUT), botJetsPosition(MAX_JETS_OUT);
0157   unsigned numberOfTopJets = 0, numberOfBotJets = 0;
0158 
0159   // Loop over local maxima
0160   for (unsigned j = 0; j < m_numberOfClusters; ++j) {
0161     // Make a proto-jet cluster
0162     L1GctRegion temp = makeProtoJet(m_localMaxima.at(j));
0163 
0164     if (m_localMaxima.at(j).rctPhi() == 0) {
0165       // Store "top edge" jets
0166       topJets.at(numberOfTopJets) = temp;
0167       topJetsPosition.at(numberOfTopJets) = 0;
0168       for (unsigned k = 0; k < numberOfTopJets; ++k) {
0169         if (topJets.at(numberOfTopJets).et() >= topJets.at(k).et()) {
0170           ++topJetsPosition.at(k);
0171         }
0172         if (topJets.at(numberOfTopJets).et() <= topJets.at(k).et()) {
0173           ++topJetsPosition.at(numberOfTopJets);
0174         }
0175       }
0176       ++numberOfTopJets;
0177     } else {
0178       // Store "bottom edge" jets
0179       botJets.at(numberOfBotJets) = temp;
0180       botJetsPosition.at(numberOfBotJets) = 0;
0181       for (unsigned k = 0; k < numberOfBotJets; ++k) {
0182         if (botJets.at(numberOfBotJets).et() >= botJets.at(k).et()) {
0183           ++botJetsPosition.at(k);
0184         }
0185         if (botJets.at(numberOfBotJets).et() <= botJets.at(k).et()) {
0186           ++botJetsPosition.at(numberOfBotJets);
0187         }
0188       }
0189       ++numberOfBotJets;
0190     }
0191   }
0192   // Now we've found all the proto-jets, copy the best ones to the output array
0193   //
0194   // We fill the first half of the array with "bottom jets"
0195   // and the remainder with "top jets". For cases where
0196   // we have found too many jets in one phi column,
0197   // we keep those with the highest Et.
0198   static const unsigned int MAX_TOPBOT_JETS = MAX_JETS_OUT / 2;
0199   unsigned pos = 0;
0200   for (unsigned j = 0; j < numberOfBotJets; ++j) {
0201     if (botJetsPosition.at(j) < MAX_TOPBOT_JETS) {
0202       m_clusters.at(pos++) = botJets.at(j);
0203     }
0204   }
0205   pos = MAX_TOPBOT_JETS;
0206   for (unsigned j = 0; j < numberOfTopJets; ++j) {
0207     if (topJetsPosition.at(j) < MAX_TOPBOT_JETS) {
0208       m_clusters.at(pos++) = topJets.at(j);
0209     }
0210   }
0211   // Finally, deal with eta00 maxima
0212   if (!m_localMax00.at(0).empty())
0213     m_cluster00.at(0) = makeProtoJet(m_localMax00.at(0));
0214   if (!m_localMax00.at(1).empty())
0215     m_cluster00.at(1) = makeProtoJet(m_localMax00.at(1));
0216 }
0217 
0218 /// Method to make a single proto-jet
0219 L1GctRegion L1GctHardwareJetFinder::makeProtoJet(L1GctRegion localMax) {
0220   unsigned eta = localMax.gctEta();
0221   unsigned phi = localMax.gctPhi();
0222   int16_t bx = localMax.bx();
0223 
0224   unsigned localEta = localMax.rctEta();
0225   unsigned localPhi = localMax.rctPhi();
0226 
0227   unsigned etCluster = 0;
0228   bool ovrFlowOr = false;
0229   bool tauVetoOr = false;
0230   unsigned rgnsAboveIsoThreshold = 0;
0231 
0232   // check for row00
0233   const unsigned midEta = (L1CaloRegionDetId::N_ETA) / 2;
0234   bool wrongEtaWheel = ((!m_positiveEtaWheel) && (eta >= midEta)) || ((m_positiveEtaWheel) && (eta < midEta));
0235 
0236   // Which rows are we looking over?
0237   unsigned rowStart, rowEnd, rowMid;
0238   static const unsigned row0 = N_EXTRA_REGIONS_ETA00 - 1;
0239   if (wrongEtaWheel) {
0240     if (localEta > row0 - 1) {
0241       rowStart = 0;
0242       rowMid = 0;
0243     } else {
0244       rowStart = row0 - 1 - localEta;
0245       rowMid = rowStart + 1;
0246     }
0247     if (localEta > row0 + 2) {  // Shouldn't happen, but big problems if it does
0248       rowEnd = 0;
0249     } else {
0250       rowEnd = row0 + 2 - localEta;
0251     }
0252   } else {
0253     rowStart = row0 + localEta;
0254     rowMid = rowStart + 1;
0255     if (localEta < COL_OFFSET - row0 - 2) {
0256       rowEnd = rowStart + 3;
0257     } else {
0258       rowEnd = COL_OFFSET;
0259     }
0260   }
0261 
0262   for (unsigned row = rowStart; row < rowEnd; ++row) {
0263     for (unsigned column = 0; column < 2; ++column) {
0264       unsigned index = column * COL_OFFSET + row;
0265       etCluster += m_inputRegions.at(index).et();
0266       ovrFlowOr |= m_inputRegions.at(index).overFlow();
0267       // Distinguish between central and tau-flagged jets. Two versions of the algorithm.
0268       if (m_useImprovedTauAlgo) {
0269         //===========================================================================================
0270         // "Old" version of improved tau algorithm tests the tau veto for the central region always
0271         //  if ((row==(localEta+N_EXTRA_REGIONS_ETA00)) && (column==localPhi)) {
0272         //    // central region - check the tau veto
0273         //    tauVetoOr |= m_inputRegions.at(index).tauVeto();
0274         //  } else {
0275         //    // other regions - check the tau veto if required
0276         //    if (!m_ignoreTauVetoBitsForIsolation) {
0277         //      tauVetoOr |= m_inputRegions.at(index).tauVeto();
0278         //    }
0279         //    // check the region energy against the isolation threshold
0280         //    if (m_inputRegions.at(index).et() >= m_tauIsolationThreshold) {
0281         //      rgnsAboveIsoThreshold++;
0282         //    }
0283         //  }
0284         //===========================================================================================
0285 
0286         // In the hardware, the ignoreTauVetoBitsForIsolation switch ignores all the veto bits,
0287         // including the one for the central region.
0288         if (!((row == rowMid) && (column == localPhi))) {
0289           // non-central region - check the region energy against the isolation threshold
0290           if (m_inputRegions.at(index).et() >= m_tauIsolationThreshold) {
0291             rgnsAboveIsoThreshold++;
0292           }
0293         }
0294         // all regions - check the tau veto if required
0295         if (!m_ignoreTauVetoBitsForIsolation) {
0296           tauVetoOr |= m_inputRegions.at(index).tauVeto();
0297         }
0298         // End of improved tau algorithm
0299       } else {
0300         // Original tau algorithm
0301         tauVetoOr |= m_inputRegions.at(index).tauVeto();
0302       }
0303     }
0304   }
0305   // Encode the number of towers over threshold for the isolated tau algorithm
0306   bool tauFeatureBit = false;
0307   if (m_useImprovedTauAlgo) {
0308     tauVetoOr |= (rgnsAboveIsoThreshold > 1);
0309     tauFeatureBit |= (rgnsAboveIsoThreshold == 1);
0310   }
0311 
0312   L1GctRegion temp(L1GctRegion::makeProtoJetRegion(etCluster, ovrFlowOr, tauVetoOr, tauFeatureBit, eta, phi, bx));
0313   return temp;
0314 }
0315 
0316 /// Convert protojets to final jets
0317 void L1GctHardwareJetFinder::findFinalClusters() {
0318   m_clusters.clear();
0319   m_clusters.resize(MAX_JETS_OUT);
0320 
0321   // Loop over proto-jets received from neighbours.
0322   // Form a jet to send to the output if there is no proto-jet nearby in the
0323   // list of jets found locally. If local jets are found nearby, form a jet
0324   // if the received jet has higher Et than any one of the local ones.
0325   for (unsigned j = 0; j < MAX_JETS_OUT; ++j) {
0326     unsigned et0 = m_rcvdProtoJets.at(j).et();
0327     unsigned localEta0 = m_rcvdProtoJets.at(j).rctEta();
0328     unsigned localPhi0 = m_rcvdProtoJets.at(j).rctPhi();
0329     unsigned JET_THRESHOLD = ((localEta0 >= m_EtaBoundry) ? m_FwdJetSeed : m_CenJetSeed);
0330     if (et0 >= JET_THRESHOLD) {
0331       bool storeJet = false;
0332       bool isolated = true;
0333       // eta00 boundary check/veto
0334       if (localEta0 == 0) {
0335         unsigned neighbourEt = m_cluster00.at(1 - localPhi0).et();
0336         isolated &= et0 >= neighbourEt;
0337       }
0338       // If the jet is NOT vetoed, look at the jets found locally (m_keptProtoJets).
0339       // We accept the jet if there are no local jets nearby, or if the local jet
0340       // (there should be no more than one) has lower Et.
0341       if (isolated) {
0342         for (unsigned k = 0; k < MAX_JETS_OUT; ++k) {
0343           unsigned et1 = m_keptProtoJets.at(k).et();
0344           unsigned localEta1 = m_keptProtoJets.at(k).rctEta();
0345           unsigned localPhi1 = m_keptProtoJets.at(k).rctPhi();
0346           if (et1 > 0) {
0347             bool distantJet = ((localPhi0 == localPhi1) || (localEta1 > localEta0 + 1) || (localEta0 > localEta1 + 1));
0348 
0349             isolated &= distantJet;
0350             storeJet |= !distantJet && ((et0 > et1) || ((et0 == et1) && localPhi0 == 1));
0351           }
0352         }
0353       }
0354 
0355       storeJet |= isolated;
0356 
0357       if (storeJet) {
0358         // Start with the et sum, tau veto and overflow flags of the protoJet (2x3 regions)
0359         unsigned etCluster = et0;
0360         bool ovrFlowOr = m_rcvdProtoJets.at(j).overFlow();
0361         bool tauVetoOr = m_rcvdProtoJets.at(j).tauVeto();
0362         unsigned rgnsAboveIsoThreshold = (m_rcvdProtoJets.at(j).featureBit0() ? 1 : 0);
0363 
0364         // Combine with the corresponding regions from
0365         // the local array to make a 3x3 jet cluster
0366         unsigned column = 1 - localPhi0;
0367         // Which rows are we looking over?
0368         unsigned rowStart, rowEnd;
0369         static const unsigned row0 = N_EXTRA_REGIONS_ETA00 - 1;
0370         rowStart = row0 + localEta0;
0371         if (localEta0 < COL_OFFSET - row0 - 2) {
0372           rowEnd = rowStart + 3;
0373         } else {
0374           rowEnd = COL_OFFSET;
0375         }
0376         unsigned index = COL_OFFSET * (this->centralCol0() + column) + rowStart;
0377         for (unsigned row = rowStart; row < rowEnd; ++row) {
0378           etCluster += m_inputRegions.at(index).et();
0379           ovrFlowOr |= m_inputRegions.at(index).overFlow();
0380           if (m_useImprovedTauAlgo) {
0381             if (!m_ignoreTauVetoBitsForIsolation) {
0382               tauVetoOr |= m_inputRegions.at(index).tauVeto();
0383             }
0384             // check the region energy against the isolation threshold
0385             if (m_inputRegions.at(index).et() >= m_tauIsolationThreshold) {
0386               rgnsAboveIsoThreshold++;
0387             }
0388           } else {
0389             tauVetoOr |= m_inputRegions.at(index).tauVeto();
0390           }
0391 
0392           ++index;
0393         }
0394 
0395         // Store the new jet
0396         unsigned eta = m_rcvdProtoJets.at(j).gctEta();
0397         unsigned phi = m_rcvdProtoJets.at(j).gctPhi();
0398         int16_t bx = m_rcvdProtoJets.at(j).bx();
0399 
0400         // Use the number of towers over threshold for the isolated tau algorithm
0401         if (m_useImprovedTauAlgo) {
0402           tauVetoOr |= (rgnsAboveIsoThreshold > 1);
0403         }
0404 
0405         L1GctRegion temp(L1GctRegion::makeFinalJetRegion(etCluster, ovrFlowOr, tauVetoOr, eta, phi, bx));
0406         m_clusters.at(j) = temp;
0407       }
0408     }
0409   }
0410 }
0411 
0412 /// Organise the pre-clustered jets into the ones we keep and those we send to the neighbour
0413 void L1GctHardwareJetFinder::convertClustersToProtoJets() {
0414   for (unsigned j = 0; j < MAX_JETS_OUT; ++j) {
0415     bool isForward = (m_clusters.at(j).rctEta() >= m_EtaBoundry);
0416     unsigned JET_THRESHOLD = (isForward ? m_FwdJetSeed : m_CenJetSeed);
0417     if (m_clusters.at(j).et() >= JET_THRESHOLD) {
0418       m_keptProtoJets.at(j) = m_clusters.at(j);
0419       m_sentProtoJets.at(j) = m_clusters.at(j);
0420     }
0421   }
0422 }
0423 
0424 /// Organise the final clustered jets into L1GctJets
0425 void L1GctHardwareJetFinder::convertClustersToOutputJets() {
0426   for (unsigned j = 0; j < MAX_JETS_OUT; ++j) {
0427     bool isForward = (m_clusters.at(j).rctEta() >= m_EtaBoundry);
0428     unsigned JET_THRESHOLD = (isForward ? m_FwdJetSeed : m_CenJetSeed);
0429     if (m_clusters.at(j).et() >= JET_THRESHOLD) {
0430       L1GctJet temp(m_clusters.at(j).et(),
0431                     m_clusters.at(j).gctEta(),
0432                     m_clusters.at(j).gctPhi(),
0433                     m_clusters.at(j).overFlow(),
0434                     isForward,
0435                     m_clusters.at(j).tauVeto(),
0436                     m_clusters.at(j).bx());
0437       m_outputJets.at(j) = temp;
0438     }
0439   }
0440 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!