File indexing completed on 2024-11-01 06:11:54
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
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
0098
0099 uint32_t expectedTotalET = 0;
0100
0101
0102 uint32_t nHitTowers = poissonRandom(100.);
0103 for (uint32_t i = 0; i < nHitTowers; i++) {
0104 uint32_t et = (random() & 0xFF);
0105 bool fg = ((random() % 100) < 95);
0106 int caloEta = ((random() + 1) % 28);
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);
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
0123 nHitTowers = poissonRandom(100.);
0124 for (uint32_t i = 0; i < nHitTowers; i++) {
0125 uint32_t et = (random() & 0xFF);
0126 uint32_t fb = (random() & 0x1F);
0127 int caloEta = ((random() + 1) % 28);
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);
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
0144 if (!uctLayer1.process()) {
0145 std::cerr << "UCT: Failed to process layer 1" << std::endl;
0146 exit(1);
0147 }
0148
0149
0150
0151
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 }