L1TCaloLayer1/test/testUCTLayer1HF.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 gaussianRandom(double mean, double standardDeviation) { return 0; }
0031
0032 double poissonRandom(double mean) {
0033   static double oldMean = -1;
0034   static double g;
0035   if (mean != oldMean) {
0036     oldMean = mean;
0037     if (mean == 0) {
0038       g = 0;
0039     } else {
0040       g = exp(-mean);
0041     }
0042   }
0043   double em = -1;
0044   double t = 1;
0045   do {
0046     em++;
0047     t *= flatRandom(0., 1.);
0048   } while (t > g);
0049   return em;
0050 }
0051
0052 void print(UCTLayer1& uct) {
0053   vector<UCTCrate*> crates = uct.getCrates();
0054   for (uint32_t crt = 0; crt < crates.size(); crt++) {
0055     vector<UCTCard*> cards = crates[crt]->getCards();
0056     for (uint32_t crd = 0; crd < cards.size(); crd++) {
0057       vector<UCTRegion*> regions = cards[crd]->getRegions();
0058       for (uint32_t rgn = 0; rgn < regions.size(); rgn++) {
0059         vector<UCTTower*> towers = regions[rgn]->getTowers();
0060         for (uint32_t twr = 0; twr < towers.size(); twr++) {
0061           if (towers[twr]->et() != 0)
0062             std::cout << *towers[twr] << std::endl;
0063         }
0064         if (regions[rgn]->et() != 0)
0065           std::cout << *regions[rgn] << std::endl;
0066       }
0067       std::cout << *cards[crd] << std::endl;
0068     }
0069     std::cout << *crates[crt] << std::endl;
0070   }
0071   std::cout << uct << std::endl;
0072 }
0073
0074 int main(int argc, char** argv) {
0075   int nEvents = 10000;
0076   if (argc == 1)
0077     std::cout << "Running on " << nEvents << std::endl;
0078   else if (argc == 2)
0079     nEvents = atoi(argv[1]);
0080   else {
0081     std::cout << "Command syntax: testUCTLayer1 [nEvents]" << std::endl;
0082     return 1;
0083   }
0084
0085   UCTLayer1 uctLayer1;
0086
0087   // Event loop for test
0088   for (int event = 0; event < nEvents; event++) {
0089     if (!uctLayer1.clearEvent()) {
0090       std::cerr << "UCT: Failed to clear event" << std::endl;
0091       exit(1);
0092     }
0093
0094     // Put a random number of towers in the UCT HF region
0095
0096     uint32_t expectedTotalET = 0;
0097     uint32_t nHitTowers = poissonRandom(3.);
0098     for (uint32_t i = 0; i < nHitTowers; i++) {
0099       uint32_t et = (random() & 0xFF);     // Random energy up to the maximum allowed
0100       bool fg = ((random() % 100) < 95);   // 5% of the time eleFG Veto should "kill" electron
0101       int caloEta = (30 + random() % 12);  // Distribute uniformly in +/- eta within acceptance
0102       while (caloEta < 30 || caloEta > 41)
0103         caloEta = (30 + random() % 12);
0104       int caloPhi;
0105       caloPhi = (1 + random() % 72);  // Distribute uniformly in all phi
0106       while (caloPhi < 1 || caloPhi > 72)
0107         caloPhi = (1 + random() % 72);
0108       caloPhi = caloPhi - ((caloPhi - 1) % 4);
0109       if ((random() & 0x1) != 0)
0110         caloEta = -caloEta;
0111       UCTTowerIndex t = UCTTowerIndex(caloEta, caloPhi);
0112       if (!uctLayer1.setHCALData(t, fg, et)) {
0113         std::cerr << "UCT: Failed loading an HF tower" << std::endl;
0114         exit(1);
0115       }
0116       expectedTotalET += et;
0117     }
0118
0119     // Process
0120     if (!uctLayer1.process()) {
0121       std::cerr << "UCT: Failed to process layer 1" << std::endl;
0122       exit(1);
0123     }
0124
0125     // Crude check if total ET is approximately OK!
0126     // We can't expect exact match as there is region level saturation to 10-bits
0127     // 10% is good enough
0128     int diff = uctLayer1.et() - expectedTotalET;
0129     if (diff < 0)
0130       diff = diff * -1;
0131     double fraction = 0.10;
0132     if (expectedTotalET < 32)
0133       fraction = 0.20;
0134     if (diff > (fraction * expectedTotalET)) {
0135       std::cout << std::dec << std::endl;
0136       print(uctLayer1);
0137       std::cout << std::endl
0138                 << "Expected " << std::showbase << std::internal << std::setfill('0') << std::setw(10) << std::hex
0139                 << expectedTotalET << std::dec << std::endl;
0140     }
0141   }
0142
0143   return 0;
0144 }