L1TCaloLayer1/test/testUCTLayer1.cpp

0001 #include <iostream>
0002 #include <iomanip>
0003 #include <vector>
0004 #include <stdlib.h>
0005 #include <stdint.h>
0006 #include <math.h>
0007
0008 using namespace std;
0009
0010 #include "L1Trigger/L1TCaloLayer1/src/UCTLayer1.hh"
0011 #include "L1Trigger/L1TCaloLayer1/src/UCTCrate.hh"
0012 #include "L1Trigger/L1TCaloLayer1/src/UCTCard.hh"
0013 #include "L1Trigger/L1TCaloLayer1/src/UCTRegion.hh"
0014 #include "L1Trigger/L1TCaloLayer1/src/UCTTower.hh"
0015
0016 #include "L1Trigger/L1TCaloLayer1/src/UCTGeometry.hh"
0017
0018 double flatRandom(double min, double max) {
0019   static double rMax = (double)0x7FFFFFFF;
0020   uint32_t r = random();
0021   double d = (double)r;
0022   double f = min + ((max - min) * d / rMax);
0023   if (f < min)
0024     f = min;
0025   if (f > max)
0026     f = max;
0027   return f;
0028 }
0029
0030 double poissonRandom(double mean) {
0031   static double oldMean = -1;
0032   static double g;
0033   if (mean != oldMean) {
0034     oldMean = mean;
0035     if (mean == 0) {
0036       g = 0;
0037     } else {
0038       g = exp(-mean);
0039     }
0040   }
0041   double em = -1;
0042   double t = 1;
0043   do {
0044     em++;
0045     t *= flatRandom(0., 1.);
0046   } while (t > g);
0047   return em;
0048 }
0049
0050 void print(UCTLayer1& uct) {
0051   vector<UCTCrate*> crates = uct.getCrates();
0052   for (uint32_t crt = 0; crt < crates.size(); crt++) {
0053     vector<UCTCard*> cards = crates[crt]->getCards();
0054     for (uint32_t crd = 0; crd < cards.size(); crd++) {
0055       vector<UCTRegion*> regions = cards[crd]->getRegions();
0056       for (uint32_t rgn = 0; rgn < regions.size(); rgn++) {
0057         if (regions[rgn]->et() > 0) {
0058           int hitEta = regions[rgn]->hitCaloEta();
0059           int hitPhi = regions[rgn]->hitCaloPhi();
0060           auto const& towers = regions[rgn]->getTowers();
0061           for (uint32_t twr = 0; twr < towers.size(); twr++) {
0062             if (towers[twr]->caloPhi() == hitPhi && towers[twr]->caloEta() == hitEta) {
0063               std::cout << "*";
0064             }
0065             std::cout << towers[twr];
0066           }
0067           std::cout << regions[rgn];
0068         }
0069       }
0070       std::cout << cards[crd];
0071     }
0072     std::cout << crates[crt];
0073   }
0074   std::cout << uct;
0075 }
0076
0077 int main(int argc, char** argv) {
0078   int nEvents = 10000;
0079   if (argc == 1)
0080     std::cout << "Running on " << nEvents << std::endl;
0081   else if (argc == 2)
0082     nEvents = atoi(argv[1]);
0083   else {
0084     std::cout << "Command syntax: testUCTLayer1 [nEvents]" << std::endl;
0085     return 1;
0086   }
0087
0088   UCTLayer1 uctLayer1;
0089
0090   // Event loop for test
0091   for (int event = 0; event < nEvents; event++) {
0092     if (!uctLayer1.clearEvent()) {
0093       std::cerr << "UCT: Failed to clear event" << std::endl;
0094       exit(1);
0095     }
0096
0097     // Put a random number of towers in the UCT
0098
0099     uint32_t expectedTotalET = 0;
0100
0101     // ECAL TPGs - set a mean of 100 random ECAL towers!
0102     uint32_t nHitTowers = poissonRandom(100.);
0103     for (uint32_t i = 0; i < nHitTowers; i++) {
0104       uint32_t et = (random() & 0xFF);      // Random energy up to the maximum allowed
0105       bool fg = ((random() % 100) < 95);    // 5% of the time eleFG Veto should "kill" electron
0106       int caloEta = ((random() + 1) % 28);  // Distribute uniformly in +/- eta within acceptance
0107       while (caloEta < 1 || caloEta > 28)
0108         caloEta = ((random() + 1) % 28);
0109       if ((random() & 0x1) != 0)
0110         caloEta = -caloEta;
0111       int caloPhi = ((random() + 1) % 72);  // Distribute uniformly in all phi
0112       while (caloPhi < 1 || caloPhi > 72)
0113         caloPhi = ((random() + 1) % 72);
0114       UCTTowerIndex t = UCTTowerIndex(caloEta, caloPhi);
0115       if (!uctLayer1.setECALData(t, fg, et)) {
0116         std::cerr << "UCT: Failed loading an ECAL tower" << std::endl;
0117         exit(1);
0118       }
0119       expectedTotalET += et;
0120     }
0121
0122     // HCAL TPGs - set a mean of 100 random HCAL towers!
0123     nHitTowers = poissonRandom(100.);
0124     for (uint32_t i = 0; i < nHitTowers; i++) {
0125       uint32_t et = (random() & 0xFF);      // Random energy up to the maximum allowed
0126       uint32_t fb = (random() & 0x1F);      // Set random five bits - this is true emulation!
0127       int caloEta = ((random() + 1) % 28);  // Distribute uniformly in +/- eta within acceptance
0128       while (caloEta < 1 || caloEta > 28)
0129         caloEta = ((random() + 1) % 28);
0130       if ((random() & 0x1) != 0)
0131         caloEta = -caloEta;
0132       int caloPhi = ((random() + 1) % 72);  // Distribute uniformly in all phi
0133       while (caloPhi < 1 || caloPhi > 72)
0134         caloPhi = ((random() + 1) % 72);
0135       UCTTowerIndex t = UCTTowerIndex(caloEta, caloPhi);
0136       if (!uctLayer1.setHCALData(t, fb, et)) {
0137         std::cerr << "UCT: Failed loading an HCAL tower" << std::endl;
0138         exit(1);
0139       }
0140       expectedTotalET += et;
0141     }
0142
0143     // Process
0144     if (!uctLayer1.process()) {
0145       std::cerr << "UCT: Failed to process layer 1" << std::endl;
0146       exit(1);
0147     }
0148
0149     // Crude check if total ET is approximately OK!
0150     // We can't expect exact match as there is region level saturation to 10-bits
0151     // 10% is good enough
0152     int diff = uctLayer1.et() - expectedTotalET;
0153     if (diff < 0)
0154       diff = diff * -1;
0155     if (diff > (0.10 * expectedTotalET)) {
0156       print(uctLayer1);
0157       std::cout << "Expected " << std::showbase << std::internal << std::setfill('0') << std::setw(10) << std::hex
0158                 << expectedTotalET << std::endl;
0159     }
0160   }
0161
0162   return 0;
0163 }