GlobalCaloTrigger/src/L1GctTdrJetFinder.cc

0001 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctTdrJetFinder.h"
0002
0003 using namespace std;
0004
0005 //DEFINE STATICS
0006 // *** Note the following definition in terms of COL_OFFSET appears not to work ***
0007 // ***          for some deep C++ reason that I don't understand - GPH          ***
0008 // const unsigned int L1GctTdrJetFinder::MAX_REGIONS_IN = L1GctJetFinderBase::COL_OFFSET*L1GctTdrJetFinder::N_COLS;
0009 // ***                        So - use the following instead                    ***
0010 const unsigned int L1GctTdrJetFinder::MAX_REGIONS_IN =
0011     (((L1CaloRegionDetId::N_ETA) / 2) + 1) * L1GctTdrJetFinder::N_COLS;
0012
0013 const unsigned int L1GctTdrJetFinder::N_COLS = 4;
0014 const unsigned int L1GctTdrJetFinder::CENTRAL_COL0 = 1;
0015
0016 L1GctTdrJetFinder::L1GctTdrJetFinder(int id) : L1GctJetFinderBase(id) {
0017   this->reset();
0018   // Initialise parameters for Region input calculations in the
0019   // derived class so we get the right values of constants.
0020   static const unsigned NPHI = L1CaloRegionDetId::N_PHI;
0021   m_minColThisJf = (NPHI + m_id * 2 - CENTRAL_COL0) % NPHI;
0022 }
0023
0024 L1GctTdrJetFinder::~L1GctTdrJetFinder() {}
0025
0026 ostream& operator<<(ostream& os, const L1GctTdrJetFinder& algo) {
0027   os << "===L1GctTdrJetFinder===" << endl;
0028   const L1GctJetFinderBase* temp = &algo;
0029   os << *temp;
0030   return os;
0031 }
0032
0033 void L1GctTdrJetFinder::fetchInput() {}
0034
0035 void L1GctTdrJetFinder::process() {
0036   if (setupOk()) {
0037     findJets();
0038     sortJets();
0039     doEnergySums();
0040   }
0041 }
0042
0043 /// HERE IS THE JETFINDER CODE
0044
0045 void L1GctTdrJetFinder::findJets() {
0046   UShort jetNum = 0;  //holds the number of jets currently found
0047   UShort centreIndex = COL_OFFSET * this->centralCol0();
0048   for (UShort column = 0; column < 2; ++column)  //Find jets in the central search region
0049   {
0050     //don't include row zero as it is not in the search region
0051     ++centreIndex;
0052     for (UShort row = 1; row < COL_OFFSET; ++row) {
0053       //Determine if we are at end of the HF or not (so need 3*2 window)
0054       bool hfBoundary = (row == COL_OFFSET - 1);
0055       //Determine if we are at the end of the endcap HCAL regions, so need boundary condition tauveto
0056       bool heBoundary = (row == COL_OFFSET - 5);
0057
0058       //debug checks for improper input indices
0059       if ((centreIndex % COL_OFFSET != 0)  //Don't want the 4 regions from other half of detector
0060           && (centreIndex >= COL_OFFSET)   //Don't want the shared column to left of jet finding area
0061           && (centreIndex < (MAX_REGIONS_IN - COL_OFFSET))) {  //Don't want column to the right either
0062
0063         if (detectJet(centreIndex, hfBoundary)) {
0064           if (jetNum < MAX_JETS_OUT) {
0065             m_outputJets.at(jetNum).setRawsum(calcJetEnergy(centreIndex, hfBoundary));
0066             m_outputJets.at(jetNum).setDetId(calcJetPosition(centreIndex));
0067             m_outputJets.at(jetNum).setBx(m_inputRegions.at(centreIndex).bx());
0068             if (row < COL_OFFSET - 4)  //if we are not in the HF, perform tauVeto analysis
0069             {
0070               m_outputJets.at(jetNum).setForward(false);
0071               m_outputJets.at(jetNum).setTauVeto(calcJetTauVeto(centreIndex, heBoundary));
0072             } else  //can't be a tau jet because we are in the HF
0073             {
0074               m_outputJets.at(jetNum).setForward(true);
0075               m_outputJets.at(jetNum).setTauVeto(true);
0076             }
0077             ++jetNum;
0078           }
0079         }
0080         ++centreIndex;
0081       }
0082     }
0083   }
0084 }
0085
0086 // Returns true if region index is the centre of a jet. Set boundary = true if at edge of HCAL.
0087 bool L1GctTdrJetFinder::detectJet(const UShort centreIndex, const bool boundary) const {
0088   if (!boundary)  //Not at boundary, so use 3*3 window of regions to determine if a jet
0089   {
0090     // Get the energy of the central region
0091     ULong testEt = m_inputRegions.at(centreIndex).et();
0092
0093     //Test if our region qualifies as a jet by comparing its energy with the energies of the
0094     //surrounding eight regions.  In the event of neighbouring regions with identical energy,
0095     //this will locate the jet in the lower-most (furthest away from eta=0), left-most (least phi) region.
0096     if (testEt > m_inputRegions.at(centreIndex - 1 - COL_OFFSET).et() &&
0097         testEt > m_inputRegions.at(centreIndex - COL_OFFSET).et() &&
0098         testEt > m_inputRegions.at(centreIndex + 1 - COL_OFFSET).et() &&
0099
0100         testEt >= m_inputRegions.at(centreIndex - 1).et() && testEt > m_inputRegions.at(centreIndex + 1).et() &&
0101
0102         testEt >= m_inputRegions.at(centreIndex - 1 + COL_OFFSET).et() &&
0103         testEt >= m_inputRegions.at(centreIndex + COL_OFFSET).et() &&
0104         testEt >= m_inputRegions.at(centreIndex + 1 + COL_OFFSET).et()) {
0105       return true;
0106     }
0107     //USE THIS BLOCK INSTEAD IF YOU WANT OVERFLOW BIT FUNCTIONALITY
0108     //*** BUT IT WILL NEED MODIFICATION SINCE L1GctRegion IS OBSOLETE ***
0109     /*    // Get the energy of the central region & OR the overflow bit to become the MSB
0110     ULong testEt = (m_inputRegions.at(centreIndex).et() | (m_inputRegions.at(centreIndex).getOverFlow() << L1GctRegion::ET_BITWIDTH));
0111
0112     //Test if our region qualifies as a jet by comparing its energy with the energies of the
0113     //surrounding eight regions.  In the event of neighbouring regions with identical energy,
0114     //this will locate the jet in the lower-most (furthest away from eta=0), left-most (least phi) region.
0115     if(testEt >  (m_inputRegions.at(centreIndex-1-COL_OFFSET).et() | (m_inputRegions.at(centreIndex-1-COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0116        testEt >  (m_inputRegions.at(centreIndex - COL_OFFSET).et() | (m_inputRegions.at(centreIndex - COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0117        testEt >  (m_inputRegions.at(centreIndex+1-COL_OFFSET).et() | (m_inputRegions.at(centreIndex+1-COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0118
0119        testEt >= (m_inputRegions.at(centreIndex - 1).et() | (m_inputRegions.at(centreIndex - 1).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0120        testEt >  (m_inputRegions.at(centreIndex + 1).et() | (m_inputRegions.at(centreIndex + 1).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0121
0122        testEt >= (m_inputRegions.at(centreIndex-1+COL_OFFSET).et() | (m_inputRegions.at(centreIndex-1+COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0123        testEt >= (m_inputRegions.at(centreIndex + COL_OFFSET).et() | (m_inputRegions.at(centreIndex + COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
0124        testEt >= (m_inputRegions.at(centreIndex+1+COL_OFFSET).et() | (m_inputRegions.at(centreIndex+1+COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)))
0125     {
0126       return true;
0127     }
0128 */  //END OVERFLOW FUNCTIONALITY
0129   } else  //...so only test surround 5 regions in our jet testing.
0130   {
0131     // Get the energy of the central region
0132     // Don't need all the overflow bit adjustments as above, since we are in the HF here
0133     ULong testEt = m_inputRegions.at(centreIndex).et();
0134
0135     if (testEt > m_inputRegions.at(centreIndex - 1 - COL_OFFSET).et() &&
0136         testEt > m_inputRegions.at(centreIndex - COL_OFFSET).et() &&
0137
0138         testEt >= m_inputRegions.at(centreIndex - 1).et() &&
0139
0140         testEt >= m_inputRegions.at(centreIndex - 1 + COL_OFFSET).et() &&
0141         testEt >= m_inputRegions.at(centreIndex + COL_OFFSET).et()) {
0142       return true;
0143     }
0144   }
0145   return false;
0146 }
0147
0148 //returns the energy sum of the nine regions centred (physically) about centreIndex
0149 L1GctJetFinderBase::ULong L1GctTdrJetFinder::calcJetEnergy(const UShort centreIndex, const bool boundary) const {
0150   ULong energy = 0;
0151
0152   if (!boundary) {
0153     for (int column = -1; column <= +1; ++column) {
0154       energy += m_inputRegions.at(centreIndex - 1 + (column * COL_OFFSET)).et() +
0155                 m_inputRegions.at(centreIndex + (column * COL_OFFSET)).et() +
0156                 m_inputRegions.at(centreIndex + 1 + (column * COL_OFFSET)).et();
0157     }
0158   } else {
0159     for (int column = -1; column <= +1; ++column) {
0160       energy += m_inputRegions.at(centreIndex - 1 + (column * COL_OFFSET)).et() +
0161                 m_inputRegions.at(centreIndex + (column * COL_OFFSET)).et();
0162     }
0163   }
0164
0165   return energy;
0166 }
0167
0168 // returns the encoded (eta, phi) position of the centre region
0169 L1CaloRegionDetId L1GctTdrJetFinder::calcJetPosition(const UShort centreIndex) const {
0170   return m_inputRegions.at(centreIndex).id();
0171 }
0172
0173 // returns the combined tauveto of the nine regions centred (physically) about centreIndex. Set boundary = true if at edge of Endcap.
0174 bool L1GctTdrJetFinder::calcJetTauVeto(const UShort centreIndex, const bool boundary) const {
0175   bool partial[3] = {false, false, false};
0176
0177   if (!boundary) {
0178     for (int column = -1; column <= +1; ++column) {
0179       partial[column + 1] = m_inputRegions.at(centreIndex - 1 + (column * COL_OFFSET)).tauVeto() ||
0180                             m_inputRegions.at(centreIndex + (column * COL_OFFSET)).tauVeto() ||
0181                             m_inputRegions.at(centreIndex + 1 + (column * COL_OFFSET)).tauVeto();
0182     }
0183   } else {
0184     for (int column = -1; column <= +1; ++column) {
0185       partial[column + 1] = m_inputRegions.at(centreIndex - 1 + (column * COL_OFFSET)).tauVeto() ||
0186                             m_inputRegions.at(centreIndex + (column * COL_OFFSET)).tauVeto();
0187     }
0188   }
0189   return partial[0] || partial[1] || partial[2];
0190 }