File indexing completed on 2024-09-07 04:36:48
0001 #include "L1Trigger/GlobalCaloTrigger/test/gctTestHt.h"
0002
0003 #include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
0004 #include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
0005
0006 #include "L1Trigger/GlobalCaloTrigger/interface/L1GlobalCaloTrigger.h"
0007 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctGlobalEnergyAlgos.h"
0008 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctJetFinderBase.h"
0009 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctWheelJetFpga.h"
0010
0011 #include <math.h>
0012 #include <iostream>
0013 #include <cassert>
0014
0015 using namespace std;
0016
0017
0018
0019
0020
0021 gctTestHt::gctTestHt()
0022 : m_bxStart(), m_numOfBx(1), minusWheelJetDta(), plusWheelJetData(), m_jetEtScale(), m_htMissScale(), m_jfPars() {}
0023
0024 gctTestHt::~gctTestHt() {}
0025
0026
0027
0028
0029
0030 void gctTestHt::configure(const L1CaloEtScale* jetScale,
0031 const L1CaloEtScale* mhtScale,
0032 const L1GctJetFinderParams* jfPars) {
0033 m_jetEtScale = jetScale;
0034 m_htMissScale = mhtScale;
0035 m_jfPars = jfPars;
0036 }
0037
0038 bool gctTestHt::setupOk() const { return (m_jetEtScale != 0 && m_htMissScale != 0 && m_jfPars != 0); }
0039
0040
0041
0042
0043 void gctTestHt::fillRawJetData(const L1GlobalCaloTrigger* gct) {
0044 minusWheelJetDta.clear();
0045 plusWheelJetData.clear();
0046 minusWheelJetDta.resize(9 * m_numOfBx);
0047 plusWheelJetData.resize(9 * m_numOfBx);
0048
0049 int bx = m_bxStart;
0050 unsigned mPos = 0;
0051 unsigned pPos = 0;
0052 for (int ibx = 0; ibx < m_numOfBx; ibx++) {
0053 int leaf = 0;
0054
0055
0056 for (; leaf < 3; leaf++) {
0057 minusWheelJetDta.at(mPos) = rawJetFinderOutput(gct->getJetLeafCards().at(leaf)->getJetFinderA(), mPos % 9, bx);
0058 mPos++;
0059 minusWheelJetDta.at(mPos) = rawJetFinderOutput(gct->getJetLeafCards().at(leaf)->getJetFinderB(), mPos % 9, bx);
0060 mPos++;
0061 minusWheelJetDta.at(mPos) = rawJetFinderOutput(gct->getJetLeafCards().at(leaf)->getJetFinderC(), mPos % 9, bx);
0062 mPos++;
0063 }
0064
0065
0066 for (; leaf < 6; leaf++) {
0067 plusWheelJetData.at(pPos) = rawJetFinderOutput(gct->getJetLeafCards().at(leaf)->getJetFinderA(), pPos % 9, bx);
0068 pPos++;
0069 plusWheelJetData.at(pPos) = rawJetFinderOutput(gct->getJetLeafCards().at(leaf)->getJetFinderB(), pPos % 9, bx);
0070 pPos++;
0071 plusWheelJetData.at(pPos) = rawJetFinderOutput(gct->getJetLeafCards().at(leaf)->getJetFinderC(), pPos % 9, bx);
0072 pPos++;
0073 }
0074
0075 bx++;
0076 }
0077 }
0078
0079
0080
0081
0082 void gctTestHt::setBxRange(const int bxStart, const int numOfBx) {
0083
0084
0085 rawJetData temp;
0086 for (int bx = bxStart; bx < m_bxStart; bx++) {
0087 minusWheelJetDta.insert(minusWheelJetDta.begin(), 9, temp);
0088 plusWheelJetData.insert(plusWheelJetData.begin(), 9, temp);
0089 }
0090
0091 m_bxStart = bxStart;
0092
0093
0094 minusWheelJetDta.resize(9 * numOfBx);
0095 plusWheelJetData.resize(9 * numOfBx);
0096 m_numOfBx = numOfBx;
0097 }
0098
0099
0100
0101
0102 bool gctTestHt::checkHtSums(const L1GlobalCaloTrigger* gct) const {
0103 if (!setupOk()) {
0104 cout << "checkHtSums setup check failed" << endl;
0105 return false;
0106 }
0107
0108 bool testPass = true;
0109 L1GctGlobalEnergyAlgos* myGlobalEnergy = gct->getEnergyFinalStage();
0110
0111 L1GctInternJetDataCollection internalJets = gct->getInternalJets();
0112
0113 std::vector<rawJetData>::const_iterator mJet = minusWheelJetDta.begin();
0114 std::vector<rawJetData>::const_iterator pJet = plusWheelJetData.begin();
0115
0116 for (int bx = 0; bx < m_numOfBx; bx++) {
0117 unsigned htMinusVl = 0;
0118 unsigned htPlusVal = 0;
0119 int hxMinusVl = 0;
0120 int hxPlusVal = 0;
0121 int hyMinusVl = 0;
0122 int hyPlusVal = 0;
0123 bool httMinusInputOf = false;
0124 bool httPlusInputOvf = false;
0125 bool htmMinusInputOf = false;
0126 bool htmPlusInputOvf = false;
0127
0128
0129
0130
0131
0132
0133 int leaf = 0;
0134 for (; leaf < 3; leaf++) {
0135 unsigned leafHttSum = 0;
0136 int leafHtxSum = 0;
0137 int leafHtySum = 0;
0138 bool leafHttOvf = false;
0139 bool leafHtmOvf = false;
0140
0141 for (int jf = 0; jf < 3; jf++) {
0142 assert(mJet != minusWheelJetDta.end());
0143 leafHttSum += (mJet->httSum);
0144 leafHttOvf |= (mJet->httOverFlow);
0145 leafHtxSum += (mJet->htxSum);
0146 leafHtySum += (mJet->htySum);
0147 leafHtmOvf |= (mJet->htmOverFlow);
0148 mJet++;
0149 }
0150 if (leafHttSum >= 4096 || leafHttOvf) {
0151 leafHttSum = 4095;
0152 leafHttOvf = true;
0153 }
0154 if (leafHttSum == gct->getJetLeafCards().at(leaf)->getAllOutputHt().at(bx).value()) {
0155 htMinusVl += leafHttSum;
0156 } else {
0157 cout << "Ht sum check leaf " << leaf << " bx " << bx << endl;
0158 testPass = false;
0159 }
0160 if (leafHtxSum == gct->getJetLeafCards().at(leaf)->getAllOutputHx().at(bx).value()) {
0161 hxMinusVl += leafHtxSum;
0162 } else {
0163 cout << "Hx sum check leaf " << leaf << " bx " << bx << endl;
0164 testPass = false;
0165 }
0166 if (leafHtySum == gct->getJetLeafCards().at(leaf)->getAllOutputHy().at(bx).value()) {
0167 hyMinusVl += leafHtySum;
0168 } else {
0169 cout << "Hy sum check leaf " << leaf << " bx " << bx << endl;
0170 testPass = false;
0171 }
0172 if ((gct->getJetLeafCards().at(leaf)->getAllOutputHt().at(bx).overFlow() == leafHttOvf) &&
0173 (gct->getJetLeafCards().at(leaf)->getAllOutputHx().at(bx).overFlow() == leafHtmOvf) &&
0174 (gct->getJetLeafCards().at(leaf)->getAllOutputHy().at(bx).overFlow() == leafHtmOvf)) {
0175 httMinusInputOf |= leafHttOvf;
0176 htmMinusInputOf |= leafHtmOvf;
0177 } else {
0178 cout << "Ht minus overflow check leaf " << leaf << " bx " << bx << endl;
0179 testPass = false;
0180 }
0181 }
0182
0183
0184 for (; leaf < 6; leaf++) {
0185 unsigned leafHttSum = 0;
0186 int leafHtxSum = 0;
0187 int leafHtySum = 0;
0188 bool leafHttOvf = false;
0189 bool leafHtmOvf = false;
0190 for (int jf = 0; jf < 3; jf++) {
0191 assert(pJet != plusWheelJetData.end());
0192 leafHttSum += (pJet->httSum);
0193 leafHttOvf |= (pJet->httOverFlow);
0194 leafHtxSum += (pJet->htxSum);
0195 leafHtySum += (pJet->htySum);
0196 leafHtmOvf |= (pJet->htmOverFlow);
0197 pJet++;
0198 }
0199 if (leafHttSum >= 4096 || leafHttOvf) {
0200 leafHttSum = 4095;
0201 leafHttOvf = true;
0202 }
0203 if (leafHttSum == gct->getJetLeafCards().at(leaf)->getAllOutputHt().at(bx).value()) {
0204 htPlusVal += leafHttSum;
0205 } else {
0206 cout << "Ht sum check leaf " << leaf << " bx " << bx << endl;
0207 testPass = false;
0208 }
0209 if (leafHtxSum == gct->getJetLeafCards().at(leaf)->getAllOutputHx().at(bx).value()) {
0210 hxPlusVal += leafHtxSum;
0211 } else {
0212 cout << "Hx sum check leaf " << leaf << " bx " << bx << endl;
0213 testPass = false;
0214 }
0215 if (leafHtySum == gct->getJetLeafCards().at(leaf)->getAllOutputHy().at(bx).value()) {
0216 hyPlusVal += leafHtySum;
0217 } else {
0218 cout << "Hy sum check leaf " << leaf << endl;
0219 testPass = false;
0220 }
0221 if ((gct->getJetLeafCards().at(leaf)->getAllOutputHt().at(bx).overFlow() == leafHttOvf) &&
0222 (gct->getJetLeafCards().at(leaf)->getAllOutputHx().at(bx).overFlow() == leafHtmOvf) &&
0223 (gct->getJetLeafCards().at(leaf)->getAllOutputHy().at(bx).overFlow() == leafHtmOvf)) {
0224 httPlusInputOvf |= leafHttOvf;
0225 htmPlusInputOvf |= leafHtmOvf;
0226 } else {
0227 cout << "Ht plus overflow check leaf " << leaf << " bx " << bx << endl;
0228 testPass = false;
0229 }
0230 }
0231
0232 unsigned htTotal = htMinusVl + htPlusVal;
0233
0234 bool httMinusOvrFlow = (htMinusVl >= 4096) || httMinusInputOf;
0235 bool httPlusOverFlow = (htPlusVal >= 4096) || httPlusInputOvf;
0236
0237 if (httMinusOvrFlow)
0238 htMinusVl = 4095;
0239 if (httPlusOverFlow)
0240 htPlusVal = 4095;
0241
0242 bool httTotalOvrFlow = (htTotal >= 4096) || httMinusOvrFlow || httPlusOverFlow;
0243
0244 if (httTotalOvrFlow)
0245 htTotal = 4095;
0246
0247 int hxTotal = hxMinusVl + hxPlusVal;
0248 int hyTotal = hyMinusVl + hyPlusVal;
0249
0250 bool htmMinusOvrFlow = htmMinusInputOf;
0251 bool htmPlusOverFlow = htmPlusInputOvf;
0252
0253
0254
0255 if (!myGlobalEnergy->getInputHtVlMinusWheel().at(bx).overFlow() && !httMinusOvrFlow &&
0256 (myGlobalEnergy->getInputHtVlMinusWheel().at(bx).value() != htMinusVl)) {
0257 cout << "ht Minus " << htMinusVl << " bx " << bx << endl;
0258 testPass = false;
0259 }
0260
0261 if (myGlobalEnergy->getInputHtVlMinusWheel().at(bx).value() != htMinusVl) {
0262 cout << "ht Minus ovF " << htMinusVl << " found " << myGlobalEnergy->getInputHtVlMinusWheel().at(bx) << " bx "
0263 << bx << endl;
0264 testPass = false;
0265 }
0266
0267 if (!myGlobalEnergy->getInputHtValPlusWheel().at(bx).overFlow() && !httPlusOverFlow &&
0268 (myGlobalEnergy->getInputHtValPlusWheel().at(bx).value() != htPlusVal)) {
0269 cout << "ht Plus " << htPlusVal << " bx " << bx << endl;
0270 testPass = false;
0271 }
0272
0273 if (myGlobalEnergy->getInputHtValPlusWheel().at(bx).value() != htPlusVal) {
0274 cout << "ht Plus ovF " << htPlusVal << " found " << myGlobalEnergy->getInputHtValPlusWheel().at(bx) << " bx "
0275 << bx << endl;
0276 testPass = false;
0277 }
0278
0279 if (myGlobalEnergy->getInputHxVlMinusWheel().at(bx).value() != hxMinusVl) {
0280 cout << "hx Minus " << hxMinusVl << " bx " << bx << endl;
0281 testPass = false;
0282 }
0283 if (myGlobalEnergy->getInputHxValPlusWheel().at(bx).value() != hxPlusVal) {
0284 cout << "hx Plus " << hxPlusVal << " bx " << bx << endl;
0285 testPass = false;
0286 }
0287
0288 if (myGlobalEnergy->getInputHyVlMinusWheel().at(bx).value() != hyMinusVl) {
0289 cout << "hy Minus " << hyMinusVl << " bx " << bx << endl;
0290 testPass = false;
0291 }
0292 if (myGlobalEnergy->getInputHyValPlusWheel().at(bx).value() != hyPlusVal) {
0293 cout << "hy Plus " << hyPlusVal << " bx " << bx << endl;
0294 testPass = false;
0295 }
0296
0297 if ((myGlobalEnergy->getInputHtVlMinusWheel().at(bx).overFlow() == httMinusOvrFlow) &&
0298 (myGlobalEnergy->getInputHxVlMinusWheel().at(bx).overFlow() == htmMinusOvrFlow) &&
0299 (myGlobalEnergy->getInputHyVlMinusWheel().at(bx).overFlow() == htmMinusOvrFlow)) {
0300 } else {
0301 cout << "Ht minus overflow check wheel"
0302 << " bx " << bx << endl;
0303 testPass = false;
0304 }
0305
0306 if ((myGlobalEnergy->getInputHtValPlusWheel().at(bx).overFlow() == httPlusOverFlow) &&
0307 (myGlobalEnergy->getInputHxValPlusWheel().at(bx).overFlow() == htmPlusOverFlow) &&
0308 (myGlobalEnergy->getInputHyValPlusWheel().at(bx).overFlow() == htmPlusOverFlow)) {
0309 } else {
0310 cout << "Ht plus overflow check wheel"
0311 << " bx " << bx << endl;
0312 testPass = false;
0313 }
0314
0315
0316 if (!myGlobalEnergy->getEtHadColl().at(bx).overFlow() && !httTotalOvrFlow &&
0317 (myGlobalEnergy->getEtHadColl().at(bx).value() != htTotal)) {
0318 cout << "Algo etHad"
0319 << " bx " << bx << endl;
0320 testPass = false;
0321 }
0322
0323 if (myGlobalEnergy->getEtHadColl().at(bx).value() != htTotal) {
0324 cout << "Algo etHad ovf"
0325 << " found " << myGlobalEnergy->getEtHadColl().at(bx) << " bx " << bx << endl;
0326 testPass = false;
0327 }
0328
0329
0330 unsigned htMiss = 0;
0331 unsigned htMPhi = 9;
0332
0333 if ((htmMinusOvrFlow || htmPlusOverFlow) || ((abs(hxTotal) > 2047) || (abs(hyTotal) > 2047))) {
0334 hxTotal = 2047;
0335 hyTotal = 2047;
0336 }
0337
0338 if ((((hxTotal) & 0xff0) != 0) || (((hyTotal) & 0xff0) != 0)) {
0339 double dhx = htComponentGeVForHtMiss(hxTotal);
0340 double dhy = htComponentGeVForHtMiss(hyTotal);
0341 double dhm = sqrt(dhx * dhx + dhy * dhy);
0342 double phi = atan2(dhy, dhx) + M_PI;
0343
0344 htMiss = m_htMissScale->rank(dhm);
0345 htMPhi = static_cast<unsigned>(phi / M_PI * 9.);
0346 }
0347
0348 if ((htMiss != myGlobalEnergy->getHtMissColl().at(bx).value()) ||
0349 (htMPhi != myGlobalEnergy->getHtMissPhiColl().at(bx).value())) {
0350 cout << "Missing Ht: expected " << htMiss << " phi " << htMPhi << " from x input " << hxTotal << " y input "
0351 << hyTotal << ", found " << myGlobalEnergy->getHtMissColl().at(bx).value() << " phi "
0352 << myGlobalEnergy->getHtMissPhiColl().at(bx).value() << " bx " << bx << endl;
0353 testPass = false;
0354 }
0355
0356
0357 unsigned htFromInternalJets = 0;
0358 bool htOvfFromInternalJets = false;
0359 for (L1GctInternJetDataCollection::const_iterator jet = internalJets.begin(); jet != internalJets.end(); jet++) {
0360 if ((jet->bx() == bx + m_bxStart) && (jet->et() >= m_jfPars->getHtJetEtThresholdGct())) {
0361 htFromInternalJets += jet->et();
0362 htOvfFromInternalJets |= jet->oflow();
0363 }
0364 }
0365 if ((htFromInternalJets >= 4096) || htOvfFromInternalJets)
0366 htFromInternalJets = 4095;
0367 if (htFromInternalJets != htTotal) {
0368 cout << "Found ht from jets " << htFromInternalJets << " expected " << htTotal << " bx " << bx << endl;
0369 testPass = false;
0370 }
0371
0372
0373 }
0374 return testPass;
0375 }
0376
0377
0378
0379
0380
0381 gctTestHt::rawJetData gctTestHt::rawJetFinderOutput(const L1GctJetFinderBase* jf,
0382 const unsigned phiPos,
0383 const int bx) const {
0384 assert(phiPos < 9);
0385 RawJetsVector jetsFromJf = jf->getRawJets();
0386 RawJetsVector jetList;
0387 unsigned sumHtt = 0;
0388 unsigned sumHtStrip0 = 0;
0389 unsigned sumHtStrip1 = 0;
0390 bool sumHttOvrFlo = false;
0391 bool sumHtmOvrFlo = false;
0392 for (RawJetsVector::const_iterator jet = jetsFromJf.begin(); jet != jetsFromJf.end(); jet++) {
0393 if (jet->bx() == bx && !jet->isNullJet()) {
0394
0395
0396
0397
0398
0399
0400
0401
0402 jetList.push_back(*jet);
0403
0404
0405
0406
0407 unsigned htJet = (jet->rawsum() == 0x3ff ? 0x3ff : jet->rawsum());
0408
0409 if (htJet >= m_jfPars->getHtJetEtThresholdGct()) {
0410 sumHtt += htJet;
0411 sumHttOvrFlo |= (jet->overFlow());
0412 }
0413
0414 if (htJet >= m_jfPars->getMHtJetEtThresholdGct()) {
0415 if (jet->rctPhi() == 0) {
0416 sumHtStrip0 += htJet;
0417 }
0418 if (jet->rctPhi() == 1) {
0419 sumHtStrip1 += htJet;
0420 }
0421 sumHtmOvrFlo |= (jet->overFlow());
0422 }
0423 }
0424 }
0425
0426 unsigned xFact0 = (53 - 4 * phiPos) % 36;
0427 unsigned xFact1 = (59 - 4 * phiPos) % 36;
0428 unsigned yFact0 = (44 - 4 * phiPos) % 36;
0429 unsigned yFact1 = (50 - 4 * phiPos) % 36;
0430
0431 int sumHtx = htComponent(sumHtStrip0, xFact0, sumHtStrip1, xFact1);
0432 int sumHty = htComponent(sumHtStrip0, yFact0, sumHtStrip1, yFact1);
0433
0434
0435 const int maxOutput = 0x800;
0436 while (sumHtx >= maxOutput) {
0437 sumHtx -= maxOutput * 2;
0438 sumHtmOvrFlo = true;
0439 }
0440 while (sumHtx < -maxOutput) {
0441 sumHtx += maxOutput * 2;
0442 sumHtmOvrFlo = true;
0443 }
0444 while (sumHty >= maxOutput) {
0445 sumHty -= maxOutput * 2;
0446 sumHtmOvrFlo = true;
0447 }
0448 while (sumHty < -maxOutput) {
0449 sumHty += maxOutput * 2;
0450 sumHtmOvrFlo = true;
0451 }
0452
0453 rawJetData result(jetList, sumHtt, sumHtx, sumHty, sumHttOvrFlo, sumHtmOvrFlo);
0454 return result;
0455 }
0456
0457 int gctTestHt::htComponent(const unsigned Emag0,
0458 const unsigned fact0,
0459 const unsigned Emag1,
0460 const unsigned fact1) const {
0461
0462 const unsigned sinFact[10] = {0, 2845, 5603, 8192, 10531, 12550, 14188, 15395, 16134, 16383};
0463 unsigned myFact;
0464 bool neg0 = false, neg1 = false, negativeResult;
0465 int res0 = 0, res1 = 0, result;
0466 unsigned Emag, fact;
0467
0468 for (int i = 0; i < 2; i++) {
0469 if (i == 0) {
0470 Emag = Emag0;
0471 fact = fact0;
0472 } else {
0473 Emag = Emag1;
0474 fact = fact1;
0475 }
0476
0477 switch (fact / 9) {
0478 case 0:
0479 myFact = sinFact[fact];
0480 negativeResult = false;
0481 break;
0482 case 1:
0483 myFact = sinFact[(18 - fact)];
0484 negativeResult = false;
0485 break;
0486 case 2:
0487 myFact = sinFact[(fact - 18)];
0488 negativeResult = true;
0489 break;
0490 case 3:
0491 myFact = sinFact[(36 - fact)];
0492 negativeResult = true;
0493 break;
0494 default:
0495 cout << "Invalid factor " << fact << endl;
0496 return 0;
0497 }
0498 result = static_cast<int>(Emag * myFact);
0499 if (i == 0) {
0500 res0 = result;
0501 neg0 = negativeResult;
0502 } else {
0503 res1 = result;
0504 neg1 = negativeResult;
0505 }
0506 }
0507 if (neg0 == neg1) {
0508 result = res0 + res1;
0509 negativeResult = neg0;
0510 } else {
0511 if (res0 >= res1) {
0512 result = res0 - res1;
0513 negativeResult = neg0;
0514 } else {
0515 result = res1 - res0;
0516 negativeResult = neg1;
0517 }
0518 }
0519
0520
0521 if (negativeResult) {
0522 result = (1 << 28) - result;
0523 result = (result + 0x1000) >> 13;
0524 result = result - (1 << 15);
0525 } else {
0526 result = (result + 0x1000) >> 13;
0527 }
0528
0529 return result;
0530 }
0531
0532 double gctTestHt::htComponentGeVForHtMiss(int inputComponent) const {
0533
0534
0535 int truncatedComponent = (((inputComponent + 0x800) >> 4) & 0xff) - 0x80;
0536 return ((static_cast<double>(truncatedComponent) + 0.5) * 8.0 * m_jfPars->getHtLsbGeV());
0537 }