GlobalCaloTrigger/test/gctTestEnergyAlgos.cc

0001 #include "L1Trigger/GlobalCaloTrigger/test/gctTestEnergyAlgos.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 #include "CondFormats/L1TObjects/interface/L1GctChannelMask.h"
0006
0007 #include "L1Trigger/GlobalCaloTrigger/interface/L1GlobalCaloTrigger.h"
0008 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctGlobalEnergyAlgos.h"
0009 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctJetFinderBase.h"
0010
0011 #include <math.h>
0012 #include <algorithm>
0013 #include <cassert>
0014 #include <iostream>
0015
0016 using namespace std;
0017
0018 //=================================================================================================================
0019 //
0020 /// Constructor and destructor
0021
0022 gctTestEnergyAlgos::gctTestEnergyAlgos()
0023     : m_bxStart(0), m_numOfBx(1), etStripSums(36), inMinusOvrFlow(1), inPlusOverFlow(1) {}
0024 gctTestEnergyAlgos::~gctTestEnergyAlgos() {}
0025
0026 //=================================================================================================================
0027 //
0028 /// Load another event into the gct. Overloaded for the various ways of doing this.
0029 //  Here's a random event generator. Loads isolated input regions to check the energy sums.
0030 std::vector<L1CaloRegion> gctTestEnergyAlgos::loadEvent(L1GlobalCaloTrigger*& gct,
0031                                                         const bool simpleEvent,
0032                                                         const int16_t bx) {
0033   static const unsigned MAX_ET_CENTRAL = 0x3ff;
0034   static const unsigned MAX_ET_FORWARD = 0x0ff;
0035   std::vector<L1CaloRegion> inputRegions;
0036
0037   // For initial tests just try things out with one region input
0038   // Then test with summing multiple regions. Choose one value
0039   // of energy and phi for each eta to avoid trying to set the
0040   // same region several times.
0041   for (unsigned i = 0; i < (simpleEvent ? 1 : L1CaloRegionDetId::N_ETA); i++) {
0042     etmiss_vec etVector = randomMissingEtVector();
0043     //    cout << "Region et " << etVector.mag << " phi " << etVector.phi << endl;
0044     // Set a single region input
0045     unsigned etaRegion = i;
0046     unsigned phiRegion = etVector.phi / 4;
0047
0048     bool regionOf = false;
0049     bool regionFg = false;
0050
0051     // Central or forward region?
0052     if (etaRegion < 4 || etaRegion >= 18) {
0053       // forward
0054       etVector.mag = etVector.mag >> 2;
0055       if (etVector.mag >= MAX_ET_FORWARD) {
0056         // deal with et values in overflow
0057         etVector.mag = MAX_ET_FORWARD;
0058       }
0059
0060     } else {
0061       // central
0062       regionOf = etVector.mag > MAX_ET_CENTRAL;
0063       regionFg = true;
0064     }
0065     //    cout << "Region et " << etVector.mag << " eta " << etaRegion << " phi " << etVector.phi << " bx " << bx << endl;
0066     // Arguments to named ctor are (et, overflow, finegrain, mip, quiet, eta, phi)
0067     L1CaloRegion temp =
0068         L1CaloRegion::makeRegionFromGctIndices(etVector.mag, regionOf, regionFg, false, false, etaRegion, phiRegion);
0069     temp.setBx(bx);
0070     inputRegions.push_back(temp);
0071   }
0072
0073   loadInputRegions(gct, inputRegions, bx);
0074   return inputRegions;
0075 }
0076
0077 // This method reads the gct input data for jetfinding from a file
0078 // as an array of region energies, one for each eta-phi bin
0079 std::vector<L1CaloRegion> gctTestEnergyAlgos::loadEvent(L1GlobalCaloTrigger*& gct,
0080                                                         const std::string& fileName,
0081                                                         bool& endOfFile,
0082                                                         const int16_t bx) {
0083   std::vector<L1CaloRegion> inputRegions;
0084
0085   //Open the file
0086   if (!regionEnergyMapInputFile.is_open()) {
0087     regionEnergyMapInputFile.open(fileName.c_str(), ios::in);
0088   }
0089
0090   //Error message and abandon ship if we can't read the file
0091   if (!regionEnergyMapInputFile.good()) {
0092     throw cms::Exception("fileReadError")
0093         << " in gctTestEnergyAlgos::loadEvent(L1GlobalCaloTrigger *&, const std::string)\n"
0094         << "Couldn't read data from file " << fileName << "!";
0095   }
0096
0097   // Here we read the energy map from the file.
0098   // Set each region at the input to the gct, and fill the expected strip sums etc.
0099   for (unsigned jphi = 0; jphi < L1CaloRegionDetId::N_PHI; ++jphi) {
0100     unsigned iphi = (L1CaloRegionDetId::N_PHI + 4 - jphi) % L1CaloRegionDetId::N_PHI;
0101     for (unsigned ieta = 0; ieta < L1CaloRegionDetId::N_ETA; ++ieta) {
0102       L1CaloRegion temp = nextRegionFromFile(ieta, iphi, bx);
0103       inputRegions.push_back(temp);
0104     }
0105   }
0106   endOfFile = regionEnergyMapInputFile.eof();
0107
0108   loadInputRegions(gct, inputRegions, bx);
0109   return inputRegions;
0110 }
0111
0112 // Load a vector of regions found elsewhere
0113 std::vector<L1CaloRegion> gctTestEnergyAlgos::loadEvent(L1GlobalCaloTrigger*& gct,
0114                                                         const std::vector<L1CaloRegion>& inputRegions,
0115                                                         const int16_t bx) {
0116   loadInputRegions(gct, inputRegions, bx);
0117   return inputRegions;
0118 }
0119
0120 //=================================================================================================================
0121 //
0122 /// Sends input regions to the gct and remembers strip sums for checking
0123 /// This routine is called from all of the above input methods
0124 void gctTestEnergyAlgos::loadInputRegions(L1GlobalCaloTrigger*& gct,
0125                                           const std::vector<L1CaloRegion>& inputRegions,
0126                                           const int16_t bx) {
0127   int bxRel = bx - m_bxStart;
0128   int base = 36 * bxRel;
0129   assert(((base >= 0) && (base + 36) <= (int)etStripSums.size()));
0130
0131   //Initialise our local sums etc for this event
0132   for (int i = 0; i < 36; i++) {
0133     etStripSums.at(i + base) = 0;
0134   }
0135   inMinusOvrFlow.at(bxRel) = false;
0136   inPlusOverFlow.at(bxRel) = false;
0137
0138   gct->fillRegions(inputRegions);
0139
0140   // Now add up the input regions
0141   for (std::vector<L1CaloRegion>::const_iterator reg = inputRegions.begin(); reg != inputRegions.end(); ++reg) {
0142     assert(reg->bx() == bx);
0143
0144     static const unsigned MAX_ET_CENTRAL = 0x3ff;
0145     static const unsigned MAX_ET_FORWARD = 0x0ff;
0146     // Check the channel masking for Et sums
0147     if (!m_chanMask->totalEtMask(reg->gctEta())) {
0148       if (reg->id().ieta() < (L1CaloRegionDetId::N_ETA / 2)) {
0149         unsigned strip = reg->id().iphi();
0150         if (reg->overFlow() || ((reg->rctEta() >= 7) && (reg->et() == MAX_ET_FORWARD))) {
0151           etStripSums.at(strip + base) += MAX_ET_CENTRAL;
0152           inMinusOvrFlow.at(bxRel) = true;
0153         } else {
0154           etStripSums.at(strip + base) += reg->et();
0155         }
0156       } else {
0157         unsigned strip = reg->id().iphi() + L1CaloRegionDetId::N_PHI;
0158         if (reg->overFlow() || ((reg->rctEta() >= 7) && (reg->et() == MAX_ET_FORWARD))) {
0159           etStripSums.at(strip + base) += MAX_ET_CENTRAL;
0160           inPlusOverFlow.at(bxRel) = true;
0161         } else {
0162           etStripSums.at(strip + base) += reg->et();
0163         }
0164       }
0165     }
0166   }
0167 }
0168
0169 //=================================================================================================================
0170 //
0171 /// Set array sizes for the number of bunch crossings
0172 void gctTestEnergyAlgos::setBxRange(const int bxStart, const int numOfBx) {
0173   // Allow the start of the bunch crossing range to be altered
0174   // without losing previously stored etStripSums
0175   for (int bx = bxStart; bx < m_bxStart; bx++) {
0176     etStripSums.insert(etStripSums.begin(), 36, (unsigned)0);
0177     inMinusOvrFlow.insert(inMinusOvrFlow.begin(), false);
0178     inPlusOverFlow.insert(inPlusOverFlow.begin(), false);
0179   }
0180
0181   m_bxStart = bxStart;
0182
0183   // Resize the vectors without clearing previously stored values
0184   etStripSums.resize(36 * numOfBx);
0185   inMinusOvrFlow.resize(numOfBx);
0186   inPlusOverFlow.resize(numOfBx);
0187   m_numOfBx = numOfBx;
0188 }
0189
0190 //=================================================================================================================
0191 //
0192 /// Check the energy sums algorithms
0193 bool gctTestEnergyAlgos::checkEnergySums(const L1GlobalCaloTrigger* gct) const {
0194   bool testPass = true;
0195   L1GctGlobalEnergyAlgos* myGlobalEnergy = gct->getEnergyFinalStage();
0196
0197   for (int bx = 0; bx < m_numOfBx; bx++) {
0198     int exMinusVl = 0;
0199     int eyMinusVl = 0;
0200     unsigned etMinusVl = 0;
0201     bool exMinusOvrFlow = inMinusOvrFlow.at(bx);
0202     bool eyMinusOvrFlow = inMinusOvrFlow.at(bx);
0203     bool etMinusOvrFlow = inMinusOvrFlow.at(bx);
0204
0205     int exPlusVal = 0;
0206     int eyPlusVal = 0;
0207     unsigned etPlusVal = 0;
0208     bool exPlusOverFlow = inPlusOverFlow.at(bx);
0209     bool eyPlusOverFlow = inPlusOverFlow.at(bx);
0210     bool etPlusOverFlow = inPlusOverFlow.at(bx);
0211
0212     unsigned rctStrip = 0;
0213     for (unsigned leaf = 0; leaf < 3; leaf++) {
0214       int exMinusSm = 0;
0215       int eyMinusSm = 0;
0216       unsigned etMinusSm = 0;
0217       int exPlusSum = 0;
0218       int eyPlusSum = 0;
0219       unsigned etPlusSum = 0;
0220
0221       unsigned etm0 = 0, etm1 = 0, etp0 = 0, etp1 = 0;
0222
0223       for (unsigned col = 0; col < 6; col++) {
0224         unsigned strip = (22 - rctStrip) % 18 + 36 * bx;
0225         unsigned etm = etStripSums.at(strip) % 4096;
0226         unsigned etp = etStripSums.at(strip + 18) % 4096;
0227         etMinusSm += etm;
0228         etPlusSum += etp;
0229
0230         if (col % 2 == 0) {
0231           etm0 = etm;
0232           etp0 = etp;
0233         } else {
0234           etm1 = etm;
0235           etp1 = etp;
0236
0237           int exm = etComponent(etm0, ((2 * strip + 11) % 36), etm1, ((2 * strip + 9) % 36));
0238           int eym = etComponent(etm0, ((2 * strip + 2) % 36), etm1, ((2 * strip) % 36));
0239
0240           int exp = etComponent(etp0, ((2 * strip + 11) % 36), etp1, ((2 * strip + 9) % 36));
0241           int eyp = etComponent(etp0, ((2 * strip + 2) % 36), etp1, ((2 * strip) % 36));
0242
0243           exMinusSm += exm;
0244           eyMinusSm += eym;
0245
0246           exPlusSum += exp;
0247           eyPlusSum += eyp;
0248         }
0249         rctStrip++;
0250       }
0251       // Check overflow for each leaf card
0252       if (exMinusSm < -65535) {
0253         exMinusSm += 131072;
0254         exMinusOvrFlow = true;
0255       }
0256       if (exMinusSm >= 65535) {
0257         exMinusSm -= 131072;
0258         exMinusOvrFlow = true;
0259       }
0260       if (eyMinusSm < -65535) {
0261         eyMinusSm += 131072;
0262         eyMinusOvrFlow = true;
0263       }
0264       if (eyMinusSm >= 65535) {
0265         eyMinusSm -= 131072;
0266         eyMinusOvrFlow = true;
0267       }
0268       if (etMinusSm >= 4096) {
0269         etMinusOvrFlow = true;
0270       }
0271       exMinusVl += exMinusSm;
0272       eyMinusVl += eyMinusSm;
0273       etMinusVl += etMinusSm;
0274       if (exPlusSum < -65535) {
0275         exPlusSum += 131072;
0276         exPlusOverFlow = true;
0277       }
0278       if (exPlusSum >= 65535) {
0279         exPlusSum -= 131072;
0280         exPlusOverFlow = true;
0281       }
0282       if (eyPlusSum < -65535) {
0283         eyPlusSum += 131072;
0284         eyPlusOverFlow = true;
0285       }
0286       if (eyPlusSum >= 65535) {
0287         eyPlusSum -= 131072;
0288         eyPlusOverFlow = true;
0289       }
0290       if (etPlusSum >= 4096) {
0291         etPlusOverFlow = true;
0292       }
0293       exPlusVal += exPlusSum;
0294       eyPlusVal += eyPlusSum;
0295       etPlusVal += etPlusSum;
0296     }
0297     // Check overflow for the overall sums
0298     if (exMinusVl < -65535) {
0299       exMinusVl += 131072;
0300       exMinusOvrFlow = true;
0301     }
0302     if (exMinusVl >= 65535) {
0303       exMinusVl -= 131072;
0304       exMinusOvrFlow = true;
0305     }
0306     if (eyMinusVl < -65535) {
0307       eyMinusVl += 131072;
0308       eyMinusOvrFlow = true;
0309     }
0310     if (eyMinusVl >= 65535) {
0311       eyMinusVl -= 131072;
0312       eyMinusOvrFlow = true;
0313     }
0314     if (etMinusVl >= 4096 || etMinusOvrFlow) {
0315       etMinusVl = 4095;
0316       etMinusOvrFlow = true;
0317     }
0318
0319     if (exPlusVal < -65535) {
0320       exPlusVal += 131072;
0321       exPlusOverFlow = true;
0322     }
0323     if (exPlusVal >= 65535) {
0324       exPlusVal -= 131072;
0325       exPlusOverFlow = true;
0326     }
0327     if (eyPlusVal < -65535) {
0328       eyPlusVal += 131072;
0329       eyPlusOverFlow = true;
0330     }
0331     if (eyPlusVal >= 65535) {
0332       eyPlusVal -= 131072;
0333       eyPlusOverFlow = true;
0334     }
0335     if (etPlusVal >= 4096 || etPlusOverFlow) {
0336       etPlusVal = 4095;
0337       etPlusOverFlow = true;
0338     }
0339
0340     int exTotal = exMinusVl + exPlusVal;
0341     int eyTotal = eyMinusVl + eyPlusVal;
0342     unsigned etTotal = etMinusVl + etPlusVal;
0343
0344     bool exTotalOvrFlow = exMinusOvrFlow || exPlusOverFlow;
0345     bool eyTotalOvrFlow = eyMinusOvrFlow || eyPlusOverFlow;
0346     bool etTotalOvrFlow = etMinusOvrFlow || etPlusOverFlow;
0347
0348     if (exTotal < -65535) {
0349       exTotal += 131072;
0350       exTotalOvrFlow = true;
0351     }
0352     if (exTotal >= 65535) {
0353       exTotal -= 131072;
0354       exTotalOvrFlow = true;
0355     }
0356     if (eyTotal < -65535) {
0357       eyTotal += 131072;
0358       eyTotalOvrFlow = true;
0359     }
0360     if (eyTotal >= 65535) {
0361       eyTotal -= 131072;
0362       eyTotalOvrFlow = true;
0363     }
0364     if (etTotal >= 4096 || etTotalOvrFlow) {
0365       etTotal = 4095;
0366       etTotalOvrFlow = true;
0367     }
0368
0369     etmiss_vec etResult = trueMissingEt(-exTotal / 2, -eyTotal / 2);
0370
0371     bool etMissOverFlow = exTotalOvrFlow || eyTotalOvrFlow;
0372     if (etMissOverFlow) {
0373       etResult.mag = 4095;
0374       etResult.phi = 45;
0375     }
0376     if (etResult.mag >= 4095) {
0377       etResult.mag = 4095;
0378       etMissOverFlow = true;
0379     }
0380
0381     //
0382     // Check the input to the final GlobalEnergyAlgos is as expected
0383     //--------------------------------------------------------------------------------------
0384     //
0385     if ((myGlobalEnergy->getInputExVlMinusWheel().at(bx).overFlow() != exMinusOvrFlow) ||
0386         (myGlobalEnergy->getInputExVlMinusWheel().at(bx).value() != exMinusVl)) {
0387       cout << "ex Minus at GlobalEnergy input " << exMinusVl << (exMinusOvrFlow ? " overflow " : "  ") << " from Gct "
0388            << myGlobalEnergy->getInputExVlMinusWheel().at(bx) << endl;
0389       testPass = false;
0390     }
0391     if ((myGlobalEnergy->getInputExValPlusWheel().at(bx).overFlow() != exPlusOverFlow) ||
0392         (myGlobalEnergy->getInputExValPlusWheel().at(bx).value() != exPlusVal)) {
0393       cout << "ex Plus at GlobalEnergy input " << exPlusVal << (exPlusOverFlow ? " overflow " : "  ") << " from Gct "
0394            << myGlobalEnergy->getInputExValPlusWheel().at(bx) << endl;
0395       testPass = false;
0396     }
0397     if ((myGlobalEnergy->getInputEyVlMinusWheel().at(bx).overFlow() != eyMinusOvrFlow) ||
0398         (myGlobalEnergy->getInputEyVlMinusWheel().at(bx).value() != eyMinusVl)) {
0399       cout << "ey Minus at GlobalEnergy input " << eyMinusVl << (eyMinusOvrFlow ? " overflow " : "  ") << " from Gct "
0400            << myGlobalEnergy->getInputEyVlMinusWheel().at(bx) << endl;
0401       testPass = false;
0402     }
0403     if ((myGlobalEnergy->getInputEyValPlusWheel().at(bx).overFlow() != eyPlusOverFlow) ||
0404         (myGlobalEnergy->getInputEyValPlusWheel().at(bx).value() != eyPlusVal)) {
0405       cout << "ey Plus at GlobalEnergy input " << eyPlusVal << (eyPlusOverFlow ? " overflow " : "  ") << " from Gct "
0406            << myGlobalEnergy->getInputEyValPlusWheel().at(bx) << endl;
0407       testPass = false;
0408     }
0409     if ((myGlobalEnergy->getInputEtVlMinusWheel().at(bx).overFlow() != etMinusOvrFlow) ||
0410         (myGlobalEnergy->getInputEtVlMinusWheel().at(bx).value() != etMinusVl)) {
0411       cout << "et Minus at GlobalEnergy input " << etMinusVl << (etMinusOvrFlow ? " overflow " : "  ") << " from Gct "
0412            << myGlobalEnergy->getInputEtVlMinusWheel().at(bx) << endl;
0413       testPass = false;
0414     }
0415     if ((myGlobalEnergy->getInputEtValPlusWheel().at(bx).overFlow() != etPlusOverFlow) ||
0416         (myGlobalEnergy->getInputEtValPlusWheel().at(bx).value() != etPlusVal)) {
0417       cout << "et Plus at GlobalEnergy input " << etPlusVal << (etPlusOverFlow ? " overflow " : "  ") << " from Gct "
0418            << myGlobalEnergy->getInputEtValPlusWheel().at(bx) << endl;
0419       testPass = false;
0420     }
0421
0422     //
0423     // Now check the processing in the final stage GlobalEnergyAlgos
0424     //--------------------------------------------------------------------------------------
0425     //
0426     // Check the missing Et calculation. Allow some margin for the
0427     // integer calculation of missing Et.
0428     unsigned etDiff, phDiff;
0429     unsigned etMargin, phMargin;
0430
0431     etDiff = (unsigned)abs((long int)etResult.mag - (long int)myGlobalEnergy->getEtMissColl().at(bx).value());
0432     phDiff = (unsigned)abs((long int)etResult.phi - (long int)myGlobalEnergy->getEtMissPhiColl().at(bx).value());
0433     if (etDiff > 2000) {
0434       etDiff = 4096 - etDiff;
0435       etMissOverFlow = true;
0436     }
0437     if (phDiff > 60) {
0438       phDiff = 72 - phDiff;
0439     }
0440     //
0441     etMargin = (etMissOverFlow ? 40 : max((etResult.mag / 100), (unsigned)1) + 2);
0442     if (etResult.mag == 0) {
0443       phMargin = 72;
0444     } else {
0445       phMargin = (30 / etResult.mag) + 1;
0446     }
0447     if ((etDiff > etMargin) || (phDiff > phMargin)) {
0448       cout << "Algo etMiss diff " << etDiff << " phi diff " << phDiff << endl;
0449       testPass = false;
0450       cout << " exTotal " << exTotal << " eyTotal " << eyTotal << endl;
0451       cout << "etMiss mag " << etResult.mag << " phi " << etResult.phi << "; from Gct mag "
0452            << myGlobalEnergy->getEtMissColl().at(bx).value() << " phi "
0453            << myGlobalEnergy->getEtMissPhiColl().at(bx).value() << endl;
0454       cout << "ex Minus " << exMinusVl << (exMinusOvrFlow ? " overflow " : "  ") << " from Gct "
0455            << myGlobalEnergy->getInputExVlMinusWheel().at(bx) << endl;
0456       cout << "ex Plus " << exPlusVal << (exPlusOverFlow ? " overflow " : "  ") << " from Gct "
0457            << myGlobalEnergy->getInputExValPlusWheel().at(bx) << endl;
0458       cout << "ey Minus " << eyMinusVl << (eyMinusOvrFlow ? " overflow " : "  ") << " from Gct "
0459            << myGlobalEnergy->getInputEyVlMinusWheel().at(bx) << endl;
0460       cout << "ey Plus " << eyPlusVal << (eyPlusOverFlow ? " overflow " : "  ") << " from Gct "
0461            << myGlobalEnergy->getInputEyValPlusWheel().at(bx) << endl;
0462       rctStrip = 0;
0463       for (unsigned i = 0; i < gct->getJetLeafCards().size(); i++) {
0464         cout << "Leaf card " << i << " ex " << gct->getJetLeafCards().at(i)->getAllOutputEx().at(bx) << " ey "
0465              << gct->getJetLeafCards().at(i)->getAllOutputEy().at(bx) << endl;
0466         cout << "strip sums ";
0467         for (unsigned col = 0; col < 6; col++) {
0468           unsigned strip = (40 - rctStrip) % 18 + 36 * bx;
0469           cout << " s " << strip << " e " << etStripSums.at(strip);
0470           rctStrip++;
0471         }
0472         cout << endl;
0473       }
0474     }
0475
0476     if (etMissOverFlow != myGlobalEnergy->getEtMissColl().at(bx).overFlow()) {
0477       cout << "etMiss overFlow " << (etMissOverFlow ? "expected but not found in Gct" : "found in Gct but not expected")
0478            << std::endl;
0479       unsigned etm0 = myGlobalEnergy->getEtMissColl().at(bx).value();
0480       unsigned etm1 = etResult.mag;
0481       cout << "etMiss value from Gct " << etm0 << "; expected " << etm1 << endl;
0482       if ((etm0 > 4090) && (etm1 > 4090) && (etm0 < 4096) && (etm1 < 4096)) {
0483         cout << "Known effect - continue testing" << endl;
0484       } else {
0485         testPass = false;
0486       }
0487     }
0488     // Check the total Et calculation
0489     if (!myGlobalEnergy->getEtSumColl().at(bx).overFlow() && !etTotalOvrFlow &&
0490         (myGlobalEnergy->getEtSumColl().at(bx).value() != etTotal)) {
0491       cout << "Algo etSum" << endl;
0492       testPass = false;
0493     }
0494     // end of loop over bunch crossings
0495   }
0496   return testPass;
0497 }
0498
0499 //=================================================================================================================
0500 //
0501 // PRIVATE MEMBER FUNCTIONS
0502 //
0503 // Function definitions for event generation
0504 //=========================================================================
0505 // Generates 2-d missing Et vector
0506 gctTestEnergyAlgos::etmiss_vec gctTestEnergyAlgos::randomMissingEtVector() const {
0507   // This produces random variables distributed as a 2-d Gaussian
0508   // with a standard deviation of sigma for each component,
0509   // and magnitude ranging up to 5*sigma.
0510   //
0511   // With sigma set to 400 we will always be in the range
0512   // of 12-bit signed integers (-2048 to 2047).
0513   // With sigma set to 200 we are always in the range
0514   // of 10-bit input region Et values.
0515   const float sigma = 400.;
0516
0517   // rmax controls the magnitude range
0518   // Chosen as a power of two conveniently close to
0519   // exp(5*5/2) to give a 5*sigma range.
0520   const unsigned rmax = 262144;
0521
0522   // Generate a pair of uniform pseudo-random integers
0523   vector<unsigned> components = randomTestData((int)2, rmax);
0524
0525   // Exclude the value zero for the first random integer
0526   // (Alternatively, return an overflow bit)
0527   while (components[0] == 0) {
0528     components = randomTestData((int)2, rmax);
0529   }
0530
0531   // Convert to the 2-d Gaussian
0532   float p, r, s;
0533   unsigned Emag, Ephi;
0534
0535   const float nbins = 18.;
0536
0537   r = float(rmax);
0538   s = r / float(components[0]);
0539   p = float(components[1]) / r;
0540   // Force phi value into the centre of a bin
0541   Emag = int(sigma * sqrt(2. * log(s)));
0542   Ephi = int(nbins * p);
0543   // Copy to the output
0544   etmiss_vec Et;
0545   Et.mag = Emag;
0546   Et.phi = (4 * Ephi);
0547   return Et;
0548 }
0549
0550 /// Generates test data consisting of a vector of energies
0551 /// uniformly distributed between zero and max
0552 vector<unsigned> gctTestEnergyAlgos::randomTestData(const int size, const unsigned max) const {
0553   vector<unsigned> energies;
0554   int r, e;
0555   float p, q, s, t;
0556
0557   p = float(max);
0558   q = float(RAND_MAX);
0559   for (int i = 0; i < size; i++) {
0560     r = rand();
0561     s = float(r);
0562     t = s * p / q;
0563     e = int(t);
0564
0565     energies.push_back(e);
0566   }
0567   return energies;
0568 }
0569
0570 // Loads test input regions from a text file.
0571 L1CaloRegion gctTestEnergyAlgos::nextRegionFromFile(const unsigned ieta, const unsigned iphi, const int16_t bx) {
0572   // The file just contains lists of region energies
0573   unsigned et;
0574   regionEnergyMapInputFile >> et;
0575   L1CaloRegion temp = L1CaloRegion::makeRegionFromGctIndices(et, false, true, false, false, ieta, iphi);
0576   temp.setBx(bx);
0577   return temp;
0578 }
0579
0580 //
0581 // Function definitions for energy sum checking
0582 //=========================================================================
0583 int gctTestEnergyAlgos::etComponent(const unsigned Emag0,
0584                                     const unsigned fact0,
0585                                     const unsigned Emag1,
0586                                     const unsigned fact1) const {
0587   // Copy the Ex, Ey conversion from the hardware emulation
0588   const unsigned sinFact[10] = {0, 2845, 5603, 8192, 10531, 12550, 14188, 15395, 16134, 16383};
0589   unsigned myFact;
0590   bool neg0 = false, neg1 = false, negativeResult;
0591   int res0 = 0, res1 = 0, result;
0592   unsigned Emag, fact;
0593
0594   for (int i = 0; i < 2; i++) {
0595     if (i == 0) {
0596       Emag = Emag0;
0597       fact = fact0;
0598     } else {
0599       Emag = Emag1;
0600       fact = fact1;
0601     }
0602
0603     switch (fact / 9) {
0604       case 0:
0605         myFact = sinFact[fact];
0606         negativeResult = false;
0607         break;
0608       case 1:
0609         myFact = sinFact[(18 - fact)];
0610         negativeResult = false;
0611         break;
0612       case 2:
0613         myFact = sinFact[(fact - 18)];
0614         negativeResult = true;
0615         break;
0616       case 3:
0617         myFact = sinFact[(36 - fact)];
0618         negativeResult = true;
0619         break;
0620       default:
0621         cout << "Invalid factor " << fact << endl;
0622         return 0;
0623     }
0624     result = static_cast<int>(Emag * myFact);
0625     if (i == 0) {
0626       res0 = result;
0627       neg0 = negativeResult;
0628     } else {
0629       res1 = result;
0630       neg1 = negativeResult;
0631     }
0632   }
0633   if (neg0 == neg1) {
0634     result = res0 + res1;
0635     negativeResult = neg0;
0636   } else {
0637     if (res0 >= res1) {
0638       result = res0 - res1;
0639       negativeResult = neg0;
0640     } else {
0641       result = res1 - res0;
0642       negativeResult = neg1;
0643     }
0644   }
0645   // Divide by 8192 using bit-shift; but emulate
0646   // twos-complement arithmetic for negative numbers
0647   if (negativeResult) {
0648     result = (1 << 28) - result;
0649     result = (result + 0x1000) >> 13;
0650     result = result - (1 << 15);
0651   } else {
0652     result = (result + 0x1000) >> 13;
0653   }
0654   return result;
0655 }
0656
0657 /// Calculate the expected missing Et vector for a given
0658 /// ex and ey sum, for comparison with the hardware
0659 gctTestEnergyAlgos::etmiss_vec gctTestEnergyAlgos::trueMissingEt(const int ex, const int ey) const {
0660   etmiss_vec result;
0661
0662   double fx = static_cast<double>(ex);
0663   double fy = static_cast<double>(ey);
0664   double fmag = sqrt(fx * fx + fy * fy);
0665   double fphi = 36. * atan2(fy, fx) / 3.1415927;
0666
0667   result.mag = static_cast<int>(fmag);
0668   if (fphi >= 0) {
0669     result.phi = static_cast<int>(fphi);
0670   } else {
0671     result.phi = static_cast<int>(fphi + 72.);
0672   }
0673
0674   return result;
0675 }