35 #define TRG_SHORT_NAMES
36 #define TRGECLCLUSTER_SHORT_NAMES
38 #include <framework/datastore/StoreArray.h>
41 #include <trg/ecl/TrgEclCluster.h>
42 #include "trg/ecl/dataobjects/TRGECLCluster.h"
50 TrgEclCluster::TrgEclCluster(): _BRICN(0), _FWDICN(0), _BWDICN(0), _BRNofCluster(0), _FWDNofCluster(0), _BWDNofCluster(0),
51 _EventId(0), _Method(1), _LimitNCluster(6), _Position(1)
76 _Quadrant.resize(3, std::vector<int>(4, 0.0));
160 _Quadrant.resize(3, std::vector<int>(4, 0.0));
198 for (
int iposition = 0; iposition < 3 ; iposition ++) {
200 for (
int icluster = 0; icluster < Ncluster; icluster++) {
212 ClusterArray[m_hitNum]->setEventId(m_nEvent);
213 ClusterArray[m_hitNum]->setClusterId(clusterId);
214 ClusterArray[m_hitNum]->setMaxTCId(
MaxTCId[iposition][icluster]);
217 ClusterArray[m_hitNum]->setNofTCinCluster(
NofTCinCluster[iposition][icluster]);
218 ClusterArray[m_hitNum]->setEnergyDep(
ClusterEnergy[iposition][icluster]);
219 ClusterArray[m_hitNum]->setTimeAve(
ClusterTiming[iposition][icluster]);
221 ClusterArray[m_hitNum]->setPositionX(
ClusterPositionX[iposition][icluster]);
222 ClusterArray[m_hitNum]->setPositionY(
ClusterPositionY[iposition][icluster]);
223 ClusterArray[m_hitNum]->setPositionZ(
ClusterPositionZ[iposition][icluster]);
233 std::vector<int> TCFire;
234 std::vector<double> TCFireEnergy;
235 std::vector<double> TCFireTiming;
236 std::vector<std::vector<double>> TCFirePosition;
239 TCFireEnergy.clear();
240 TCFireTiming.clear();
241 TCFirePosition.clear();
243 TCFire.resize(432, 0);
244 TCFireEnergy.resize(432, 0.);
245 TCFireTiming.resize(432, 0.);
246 TCFirePosition.resize(432, std::vector<double>(3, 0.));
251 const int hit_size =
TCId.size();
252 for (
int ihit = 0 ; ihit < hit_size ; ihit++) {
253 if (
TCId[ihit] >= 81 &&
TCId[ihit] <= 512) {
254 TCFire[
TCId[ihit] - 81] =
TCId[ihit];
268 int tc_upper_right = 0;
270 int tc_lower_right = 0;
272 int tc_lower_left = 0;
274 int tc_upper_left = 0;
276 for (
int iii = 0 ; iii < 432 ; iii++) {
277 if (TCFire[iii] == 0) {
continue; }
280 tc_upper = TCFire[iii + 420] ;
281 tc_upper_right = TCFire[iii + 419] ;
282 tc_right = TCFire[iii - 1] ;
283 tc_lower_right = TCFire[iii + 11] ;
284 tc_lower = TCFire[iii + 12] ;
285 tc_lower_left = TCFire[iii + 13] ;
286 tc_left = TCFire[iii + 1] ;
287 tc_upper_left = TCFire[iii + 421] ;
294 if (iii % 12 == 11) {
301 if (iii > 11 && iii < 420) {
302 tc_upper = TCFire[iii - 12] ;
303 tc_upper_right = TCFire[iii - 13] ;
304 tc_right = TCFire[iii - 1] ;
305 tc_lower_right = TCFire[iii + 11] ;
306 tc_lower = TCFire[iii + 12] ;
307 tc_lower_left = TCFire[iii + 13] ;
308 tc_left = TCFire[iii + 1] ;
309 tc_upper_left = TCFire[iii - 11] ;
316 if (iii % 12 == 11) {
325 tc_upper = TCFire[iii - 12] ;
326 tc_upper_right = TCFire[iii - 13] ;
327 tc_right = TCFire[iii - 1] ;
328 tc_lower_right = TCFire[iii - 421] ;
329 tc_lower = TCFire[iii - 420] ;
330 tc_lower_left = TCFire[iii - 419];
331 tc_left = TCFire[iii + 1] ;
332 tc_upper_left = TCFire[iii - 11] ;
338 if (iii % 12 == 11) {
358 if (!(tc_upper != 0 || tc_left != 0)) {
359 if (!(tc_lower != 0 && tc_lower_left != 0)) {
364 double maxTCEnergy = 0;
365 for (
int iTC = 0; iTC < 9; iTC++) {
367 if (maxTCEnergy < TCFireEnergy[
TempCluster[iTC] - 81]) {
386 if ((maxTCid - 81) % 12 == 0) {
393 if ((maxTCid - 81) % 12 == 11) {
401 if (maxTCid > 92 && maxTCid < 501) {
412 if ((maxTCid - 81) % 12 == 0) {
418 if ((maxTCid - 81) % 12 == 11) {
439 if ((maxTCid - 81) % 12 == 0) {
445 if ((maxTCid - 81) % 12 == 11) {
458 for (
int iNearTC = 1; iNearTC < 9; iNearTC ++) {
459 for (
int jNearTC = 1; jNearTC < 9; jNearTC ++) {
461 if (iNearTC == jNearTC) {
continue;}
469 double clusterenergy = 0;
470 double clustertiming = 0;
471 double clusterpositionX = 0;
472 double clusterpositionY = 0;
473 double clusterpositionZ = 0;
474 int noftcincluster = 0;
475 for (
int iNearTC = 0; iNearTC < 9; iNearTC ++) {
477 else {noftcincluster++;}
478 clusterenergy += TCFireEnergy[
TempCluster[iNearTC] - 81];
480 clusterpositionX += TCFireEnergy[
TempCluster[iNearTC] - 81] * TCFirePosition[
TempCluster[iNearTC] - 81][0];
481 clusterpositionY += TCFireEnergy[
TempCluster[iNearTC] - 81] * TCFirePosition[
TempCluster[iNearTC] - 81][1];
482 clusterpositionZ += TCFireEnergy[
TempCluster[iNearTC] - 81] * TCFirePosition[
TempCluster[iNearTC] - 81][2];
490 clustertiming /= clusterenergy;
492 clusterpositionX /= clusterenergy;
493 clusterpositionY /= clusterenergy;
494 clusterpositionZ /= clusterenergy;
496 clustertiming = TCFireTiming[maxTCId - 81];
497 clusterpositionX = TCFirePosition[maxTCId - 81][0];
498 clusterpositionY = TCFirePosition[maxTCId - 81][1];
499 clusterpositionZ = TCFirePosition[maxTCId - 81][2];
501 if (clustertiming == 0 && clusterenergy == 0) {
continue;}
520 std::vector<int> TCFire;
521 std::vector<double> TCFireEnergy;
522 std::vector<double> TCFireTiming;
523 std::vector<std::vector<double>> TCFirePosition;
525 std::vector<double> TempClusterEnergy;
526 std::vector<double> TempClusterTiming;
527 std::vector<double> TempClusterPositionX;
528 std::vector<double> TempClusterPositionY;
529 std::vector<double> TempClusterPositionZ;
530 std::vector<double> Sort1D;
531 std::vector<std::vector<double>> Sort2D;
533 std::vector<int> TempNofTCinCluster;
534 std::vector<int> TempMaxTCId;
536 int TempICNTCId = 0;;
538 TempClusterEnergy.clear();
539 TempClusterTiming.clear();
540 TempClusterPositionX.clear();
541 TempClusterPositionY.clear();
542 TempClusterPositionZ.clear();
543 TempNofTCinCluster.clear();
554 TCFireEnergy.clear();
555 TCFireTiming.clear();
556 TCFirePosition.clear();
558 TCFire.resize(96, 0);
559 TCFireEnergy.resize(96, 0.);
560 TCFireTiming.resize(96, 0.);
562 TCFirePosition.resize(96, std::vector<double>(3, 0.));
565 const int hit_size =
TCId.size();
566 for (
int ihit = 0 ; ihit < hit_size ; ihit++) {
567 if (
TCId[ihit] > 80) {
continue;}
593 int iTCId0 =
TCId[ihit] - 1;
595 if (iTCId0 % 5 == 0) {
596 kkk = (iTCId0 / 5) * 2;
597 TCFire[kkk] =
TCId[ihit];
598 TCFire[kkk + 1] =
TCId[ihit];
601 switch (iTCId0 % 5) {
603 TCFire[32 + 2 * kkk] =
TCId[ihit];
break;
605 TCFire[64 + 2 * kkk] =
TCId[ihit];
break;
607 TCFire[64 + 2 * kkk + 1] =
TCId[ihit];
break;
609 TCFire[32 + 2 * kkk + 1] =
TCId[ihit];
break;
616 for (
int iii = 0 ; iii < 96 ; iii++) {
623 for (
int iinit = 0; iinit < 9; iinit ++) {
TempCluster[iinit] = 0;}
624 if (TCFire[iii] == 0) {
continue; }
637 }
else if (iii == 30) {
668 }
else if (iii > 31 && iii < 64) {
680 }
else if (iii == 33) {
694 else if (iii == 62) {
706 }
else if (iii == 63) {
757 }
else if (iii == 95) {
790 TempICNTCId = TCFire[iii];
794 double maxTCEnergy = 0;
795 for (
int iTC = 0; iTC < 9; iTC++) {
797 if (maxTCEnergy < TCFireEnergy[
TempCluster[iTC] - 1]) {
812 kkk = 64 + 1 + 2 * kkk;
814 kkk = 32 + 1 + 2 * kkk;
829 }
else if (kkk == 30) {
853 }
else if (kkk > 31 && kkk < 64) {
865 }
else if (kkk == 33) {
879 else if (kkk == 62) {
891 }
else if (kkk == 63) {
935 }
else if (kkk == 95) {
960 for (
int iNearTC = 1; iNearTC < 9; iNearTC ++) {
961 for (
int jNearTC = 1; jNearTC < 9; jNearTC ++) {
963 if (iNearTC == jNearTC)
continue;
971 double clusterenergy = 0;
972 double clustertiming = 0;
973 double clusterpositionX = 0;
974 double clusterpositionY = 0;
975 double clusterpositionZ = 0;
976 int noftcincluster = 0;
977 for (
int iNearTC = 0; iNearTC < 9; iNearTC ++) {
979 else {noftcincluster++;}
980 clusterenergy += TCFireEnergy[
TempCluster[iNearTC] - 1];
991 clustertiming /= clusterenergy;
993 clusterpositionX /= clusterenergy;
994 clusterpositionY /= clusterenergy;
995 clusterpositionZ /= clusterenergy;
997 clustertiming = TCFireTiming[maxTCId - 1];
998 clusterpositionX = TCFirePosition[maxTCId - 1][0];
999 clusterpositionY = TCFirePosition[maxTCId - 1][1];
1000 clusterpositionZ = TCFirePosition[maxTCId - 1][2];
1002 if (clustertiming == 0 && clusterenergy == 0) {
continue;}
1004 TempClusterEnergy.push_back(clusterenergy);
1005 TempClusterTiming.push_back(clustertiming);
1006 TempClusterPositionX.push_back(clusterpositionX);
1007 TempClusterPositionY.push_back(clusterpositionY);
1008 TempClusterPositionZ.push_back(clusterpositionZ);
1009 TempNofTCinCluster.push_back(noftcincluster);
1010 TempMaxTCId.push_back(maxTCId);
1013 Sort1D.push_back(TempICNTCId);
1014 Sort1D.push_back(clusterenergy);
1015 Sort1D.push_back(clustertiming);
1016 Sort1D.push_back(clusterpositionX);
1017 Sort1D.push_back(clusterpositionY);
1018 Sort1D.push_back(clusterpositionZ);
1019 Sort1D.push_back(noftcincluster);
1020 Sort1D.push_back(maxTCId);
1022 Sort2D.push_back(Sort1D);
1031 sort(Sort2D.begin(), Sort2D.end(),
1032 [](
const vector<double>& aa1,
const vector<double>& aa2) {return aa1[0] < aa2[0];});
1034 const int clustersize = Sort2D.size();
1035 for (
int jtc = 0; jtc < clustersize; jtc++) {
1043 MaxTCId[1].push_back(Sort2D[jtc][7]);
1059 std::vector<int> TCFire;
1060 std::vector<double> TCFireEnergy;
1061 std::vector<double> TCFireTiming;
1062 std::vector<std::vector<double>> TCFirePosition;
1065 std::vector<double> TempClusterEnergy;
1066 std::vector<double> TempClusterTiming;
1067 std::vector<double> TempClusterPositionX;
1068 std::vector<double> TempClusterPositionY;
1069 std::vector<double> TempClusterPositionZ;
1070 std::vector<int> TempNofTCinCluster;
1071 std::vector<int> TempMaxTCId;
1073 std::vector<double> Sort1D;
1074 std::vector<std::vector<double>> Sort2D;
1076 int TempICNTCId = 0;
1078 TempClusterEnergy.clear();
1079 TempClusterTiming.clear();
1080 TempClusterPositionX.clear();
1081 TempClusterPositionY.clear();
1082 TempClusterPositionZ.clear();
1083 TempNofTCinCluster.clear();
1084 TempMaxTCId.clear();
1091 TCFireEnergy.clear();
1092 TCFireTiming.clear();
1093 TCFirePosition.clear();
1095 TCFire.resize(64, 0);
1096 TCFireEnergy.resize(64, 0.);
1097 TCFireTiming.resize(64, 0.);
1098 TCFirePosition.resize(64, std::vector<double>(3, 0.));
1102 const int hit_size =
TCId.size();
1103 for (
int ihit = 0 ; ihit < hit_size ; ihit++) {
1104 if (
TCId[ihit] < 513) {
continue;}
1105 TCFireEnergy[
TCId[ihit] - 513] =
Energy[ihit];
1106 TCFireTiming[
TCId[ihit] - 513] =
Timing[ihit];
1127 int iTCId0 =
TCId[ihit] - 1;
1129 if ((iTCId0 - 512) % 4 == 2) {
1130 kkk = (iTCId0 - 512) / 2 - 1;
1132 if ((iTCId0 - 512) % 4 == 1) {
1133 kkk = ((iTCId0 - 512) + 1) / 2;
1135 if ((iTCId0 - 512) % 4 == 3) {
1136 kkk = 32 + ((iTCId0 - 512) - 3) / 2;
1138 if ((iTCId0 - 512) % 4 == 0) {
1139 kkk = 33 + ((iTCId0 - 512)) / 2;
1142 TCFire[kkk] = iTCId0 + 1;
1147 for (
int iii = 0 ; iii < 64 ; iii ++) {
1149 if (TCFire[iii] == 0) {
continue; }
1151 for (
int iinit = 0; iinit < 9; iinit ++) {
TempCluster[iinit] = 0;}
1163 }
else if (iii == 31) {
1184 }
else if (iii > 31) {
1195 }
else if (iii == 63) {
1221 TempICNTCId = TCFire[iii];
1225 double maxTCEnergy = 0;
1226 for (
int iTC = 0; iTC < 9; iTC++) {
1228 if (maxTCEnergy < TCFireEnergy[
TempCluster[iTC] - 513]) {
1229 maxTCEnergy = TCFireEnergy[
TempCluster[iTC] - 513];
1259 }
else if (kkk == 31) {
1280 }
else if (kkk > 31) {
1291 }
else if (kkk == 63) {
1316 for (
int iNearTC = 1; iNearTC < 9; iNearTC ++) {
1317 for (
int jNearTC = 1; jNearTC < 9; jNearTC ++) {
1319 if (iNearTC == jNearTC) {
continue;}
1327 double clusterenergy = 0;
1328 double clustertiming = 0;
1329 double clusterpositionX = 0;
1330 double clusterpositionY = 0;
1331 double clusterpositionZ = 0;
1332 int noftcincluster = 0;
1333 for (
int iNearTC = 0; iNearTC < 9; iNearTC ++) {
1335 else {noftcincluster++;}
1337 clusterenergy += TCFireEnergy[
TempCluster[iNearTC] - 513];
1339 clusterpositionX += TCFireEnergy[
TempCluster[iNearTC] - 513] * TCFirePosition[
TempCluster[iNearTC] - 513][0];
1340 clusterpositionY += TCFireEnergy[
TempCluster[iNearTC] - 513] * TCFirePosition[
TempCluster[iNearTC] - 513][1];
1341 clusterpositionZ += TCFireEnergy[
TempCluster[iNearTC] - 513] * TCFirePosition[
TempCluster[iNearTC] - 513][2];
1349 clustertiming /= clusterenergy;
1351 clusterpositionX /= clusterenergy;
1352 clusterpositionY /= clusterenergy;
1353 clusterpositionZ /= clusterenergy;
1355 clustertiming = TCFireTiming[maxTCId - 513];
1356 clusterpositionX = TCFirePosition[maxTCId - 513][0];
1357 clusterpositionY = TCFirePosition[maxTCId - 513][1];
1358 clusterpositionZ = TCFirePosition[maxTCId - 513][2];
1361 if (clustertiming == 0 && clusterenergy == 0) {
continue;}
1362 TempClusterEnergy.push_back(clusterenergy);
1363 TempClusterTiming.push_back(clustertiming);
1364 TempClusterPositionX.push_back(clusterpositionX);
1365 TempClusterPositionY.push_back(clusterpositionY);
1366 TempClusterPositionZ.push_back(clusterpositionZ);
1367 TempNofTCinCluster.push_back(noftcincluster);
1368 TempMaxTCId.push_back(maxTCId);
1371 Sort1D.push_back(TempICNTCId);
1372 Sort1D.push_back(clusterenergy);
1373 Sort1D.push_back(clustertiming);
1374 Sort1D.push_back(clusterpositionX);
1375 Sort1D.push_back(clusterpositionY);
1376 Sort1D.push_back(clusterpositionZ);
1377 Sort1D.push_back(noftcincluster);
1378 Sort1D.push_back(maxTCId);
1380 Sort2D.push_back(Sort1D);
1388 sort(Sort2D.begin(), Sort2D.end(),
1389 [](
const vector<double>& aa1,
const vector<double>& aa2) {return aa1[0] < aa2[0];});
1391 const int clustersize = Sort2D.size();
1392 for (
int jtc = 0; jtc < clustersize; jtc++) {
1400 MaxTCId[2].push_back(Sort2D[jtc][7]);
1418 std::vector<int> TCFire;
1422 TCFire.resize(432, 0);
1424 const int hit_size =
TCId.size();
1425 for (
int ihit = 0 ; ihit < hit_size ; ihit++) {
1426 if (
TCId[ihit] >= 81 &&
TCId[ihit] <= 512) {
1427 TCFire[
TCId[ihit] - 81] =
TCId[ihit];
1435 int tc_upper_right = 0;
1437 int tc_lower_right = 0;
1439 int tc_lower_left = 0;
1441 int tc_upper_left = 0;
1445 for (
int iii = 0 ; iii < 432 ; iii++) {
1446 if (TCFire[iii] == 0) {
continue; }
1449 tc_upper = TCFire[iii + 420] ;
1450 tc_upper_right = TCFire[iii + 419] ;
1451 tc_right = TCFire[iii - 1] ;
1452 tc_lower_right = TCFire[iii + 11] ;
1453 tc_lower = TCFire[iii + 12] ;
1454 tc_lower_left = TCFire[iii + 13] ;
1455 tc_left = TCFire[iii + 1] ;
1456 tc_upper_left = TCFire[iii + 421] ;
1457 if (iii % 12 == 0) {
1463 if (iii % 12 == 11) {
1470 if (iii > 11 && iii < 420) {
1471 tc_upper = TCFire[iii - 12] ;
1472 tc_upper_right = TCFire[iii - 13] ;
1473 tc_right = TCFire[iii - 1] ;
1474 tc_lower_right = TCFire[iii + 11] ;
1475 tc_lower = TCFire[iii + 12] ;
1476 tc_lower_left = TCFire[iii + 13] ;
1477 tc_left = TCFire[iii + 1] ;
1478 tc_upper_left = TCFire[iii - 11] ;
1479 if (iii % 12 == 0) {
1485 if (iii % 12 == 11) {
1494 tc_upper = TCFire[iii - 12] ;
1495 tc_upper_right = TCFire[iii - 13] ;
1496 tc_right = TCFire[iii - 1] ;
1497 tc_lower_right = TCFire[iii - 421] ;
1498 tc_lower = TCFire[iii - 420] ;
1499 tc_lower_left = TCFire[iii - 419];
1500 tc_left = TCFire[iii + 1] ;
1501 tc_upper_left = TCFire[iii - 11] ;
1502 if (iii % 12 == 0) {
1507 if (iii % 12 == 11) {
1528 if (!(tc_upper != 0 || tc_left != 0)) {
1529 if (!(tc_lower != 0 && tc_lower_left != 0)) {
1531 int phiid =
_TCMap -> getTCPhiIdFromTCId(iii + 80 + 1);
1533 if (phiid == 36 || (phiid > 0 && phiid < 11)) {
1536 if (phiid > 8 && phiid < 20) {
1539 if (phiid > 17 && phiid < 29) {
1542 if ((phiid > 26 && phiid < 37) || phiid == 1) {
1559 std::vector<int> TCFire;
1564 TCFire.resize(96, 0);
1567 const int hit_size =
TCId.size();
1568 for (
int ihit = 0 ; ihit < hit_size ; ihit++) {
1569 if (
TCId[ihit] > 80) {
continue;}
1588 int iTCId0 =
TCId[ihit] - 1;
1590 if (iTCId0 % 5 == 0) {
1591 kkk = (iTCId0 / 5) * 2;
1592 TCFire[kkk] =
TCId[ihit];
1593 TCFire[kkk + 1] =
TCId[ihit];
1596 switch (iTCId0 % 5) {
1598 TCFire[32 + 2 * kkk] =
TCId[ihit];
break;
1600 TCFire[64 + 2 * kkk] =
TCId[ihit];
break;
1602 TCFire[64 + 2 * kkk + 1] =
TCId[ihit];
break;
1604 TCFire[32 + 2 * kkk + 1] =
TCId[ihit];
break;
1612 for (
int iii = 32 ; iii < 96 ; iii++) {
1613 for (
int iinit = 0; iinit < 9; iinit ++) {
TempCluster[iinit] = 0;}
1614 if (TCFire[iii] == 0) {
continue; }
1626 }
else if (iii == 63) {
1658 }
else if (iii == 95) {
1685 int phiid =
_TCMap -> getTCPhiIdFromTCId(TCFire[iii]);
1687 if (phiid == 32 || (phiid > 0 && phiid < 10)) {
1690 if (phiid > 7 && phiid < 18) {
1693 if (phiid > 15 && phiid < 26) {
1696 if ((phiid > 22 && phiid < 33) || phiid == 1) {
1711 std::vector<int> TCFire;
1716 TCFire.resize(64, 0);
1720 const int hit_size =
TCId.size();
1721 for (
int ihit = 0 ; ihit < hit_size ; ihit++) {
1722 if (
TCId[ihit] < 513) {
continue;}
1741 int iTCId0 =
TCId[ihit] - 1;
1743 if ((iTCId0 - 512) % 4 == 2) {
1744 kkk = (iTCId0 - 512) / 2 - 1;
1746 if ((iTCId0 - 512) % 4 == 1) {
1747 kkk = ((iTCId0 - 512) + 1) / 2;
1749 if ((iTCId0 - 512) % 4 == 3) {
1750 kkk = 32 + ((iTCId0 - 512) - 3) / 2;
1752 if ((iTCId0 - 512) % 4 == 0) {
1753 kkk = 33 + ((iTCId0 - 512)) / 2;
1756 TCFire[kkk] = iTCId0 + 1;
1761 for (
int iii = 0 ; iii < 64 ; iii ++) {
1763 if (TCFire[iii] == 0) {
continue; }
1765 for (
int iinit = 0; iinit < 9; iinit ++) {
TempCluster[iinit] = 0;}
1777 }
else if (iii == 31) {
1798 }
else if (iii > 31) {
1809 }
else if (iii == 63) {
1836 int phiid =
_TCMap -> getTCPhiIdFromTCId(TCFire[iii]);
1838 if (phiid == 32 || (phiid > 0 && phiid < 10)) {
1841 if (phiid > 7 && phiid < 18) {
1844 if (phiid > 15 && phiid < 26) {
1847 if ((phiid > 22 && phiid < 33) || phiid == 1) {