File indexing completed on 2024-04-06 12:19:54
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
0021
0022 gctTestEnergyAlgos::gctTestEnergyAlgos()
0023 : m_bxStart(0), m_numOfBx(1), etStripSums(36), inMinusOvrFlow(1), inPlusOverFlow(1) {}
0024 gctTestEnergyAlgos::~gctTestEnergyAlgos() {}
0025
0026
0027
0028
0029
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
0038
0039
0040
0041 for (unsigned i = 0; i < (simpleEvent ? 1 : L1CaloRegionDetId::N_ETA); i++) {
0042 etmiss_vec etVector = randomMissingEtVector();
0043
0044
0045 unsigned etaRegion = i;
0046 unsigned phiRegion = etVector.phi / 4;
0047
0048 bool regionOf = false;
0049 bool regionFg = false;
0050
0051
0052 if (etaRegion < 4 || etaRegion >= 18) {
0053
0054 etVector.mag = etVector.mag >> 2;
0055 if (etVector.mag >= MAX_ET_FORWARD) {
0056
0057 etVector.mag = MAX_ET_FORWARD;
0058 }
0059
0060 } else {
0061
0062 regionOf = etVector.mag > MAX_ET_CENTRAL;
0063 regionFg = true;
0064 }
0065
0066
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
0078
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
0086 if (!regionEnergyMapInputFile.is_open()) {
0087 regionEnergyMapInputFile.open(fileName.c_str(), ios::in);
0088 }
0089
0090
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
0098
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
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
0123
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
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
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
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
0172 void gctTestEnergyAlgos::setBxRange(const int bxStart, const int numOfBx) {
0173
0174
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
0184 etStripSums.resize(36 * numOfBx);
0185 inMinusOvrFlow.resize(numOfBx);
0186 inPlusOverFlow.resize(numOfBx);
0187 m_numOfBx = numOfBx;
0188 }
0189
0190
0191
0192
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
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
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
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
0424
0425
0426
0427
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
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
0495 }
0496 return testPass;
0497 }
0498
0499
0500
0501
0502
0503
0504
0505
0506 gctTestEnergyAlgos::etmiss_vec gctTestEnergyAlgos::randomMissingEtVector() const {
0507
0508
0509
0510
0511
0512
0513
0514
0515 const float sigma = 400.;
0516
0517
0518
0519
0520 const unsigned rmax = 262144;
0521
0522
0523 vector<unsigned> components = randomTestData((int)2, rmax);
0524
0525
0526
0527 while (components[0] == 0) {
0528 components = randomTestData((int)2, rmax);
0529 }
0530
0531
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
0541 Emag = int(sigma * sqrt(2. * log(s)));
0542 Ephi = int(nbins * p);
0543
0544 etmiss_vec Et;
0545 Et.mag = Emag;
0546 Et.phi = (4 * Ephi);
0547 return Et;
0548 }
0549
0550
0551
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
0571 L1CaloRegion gctTestEnergyAlgos::nextRegionFromFile(const unsigned ieta, const unsigned iphi, const int16_t bx) {
0572
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
0582
0583 int gctTestEnergyAlgos::etComponent(const unsigned Emag0,
0584 const unsigned fact0,
0585 const unsigned Emag1,
0586 const unsigned fact1) const {
0587
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
0646
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
0658
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 }