File indexing completed on 2025-07-06 08:34:51
0001
0002
0003
0004
0005
0006
0007 #include <vector>
0008 #include <TVector3.h>
0009
0010 struct towersStrct {
0011 towersStrct()
0012 : energy(0)
0013 , time(0)
0014 , posx(0)
0015 , posy(0)
0016 , posz(0)
0017 , cellID(0)
0018 , cellIDx(-1)
0019 , cellIDy(-1)
0020 , cellIDz(-1)
0021 , tower_trueID(-10000)
0022 , tower_clusterIDA(-1)
0023 , tower_clusterIDB(-1) {}
0024 float energy;
0025 float time;
0026 float posx;
0027 float posy;
0028 float posz;
0029 int cellID;
0030 int cellIDx;
0031 int cellIDy;
0032 int cellIDz;
0033 int tower_trueID;
0034 int tower_clusterIDA;
0035 int tower_clusterIDB;
0036 };
0037
0038 bool acompare(towersStrct lhs, towersStrct rhs) { return lhs.energy > rhs.energy; }
0039
0040 struct clustersStrct {
0041 clustersStrct()
0042 : cluster_E(0.)
0043 , cluster_seed(0.)
0044 , cluster_Eta(-10.)
0045 , cluster_Phi(-10.)
0046 , cluster_X(0.)
0047 , cluster_Y(0.)
0048 , cluster_Z(0.)
0049 , cluster_M02(0.)
0050 , cluster_M20(0.)
0051 , cluster_NTowers(0)
0052 , cluster_trueID(-10000)
0053 , cluster_NtrueID(0) {}
0054 float cluster_E;
0055 float cluster_seed;
0056 float cluster_Eta;
0057 float cluster_Phi;
0058 float cluster_X;
0059 float cluster_Y;
0060 float cluster_Z;
0061 float cluster_M02;
0062 float cluster_M20;
0063 int cluster_NTowers;
0064 int cluster_trueID;
0065 int cluster_NtrueID;
0066 std::vector<towersStrct> cluster_towers;
0067 };
0068
0069 bool acompareCl(clustersStrct lhs, clustersStrct rhs) { return lhs.cluster_E > rhs.cluster_E; }
0070
0071
0072
0073
0074
0075
0076 clustersStrct findMACluster(
0077 float seed,
0078 float ,
0079 std::vector<towersStrct>& input_towers_temp,
0080 std::vector<towersStrct>& cluster_towers_temp,
0081
0082 float aggMargin = 1.0
0083 ) {
0084 clustersStrct tempstructC;
0085 if (input_towers_temp.at(0).energy > seed) {
0086
0087
0088 tempstructC.cluster_E = input_towers_temp.at(0).energy;
0089 tempstructC.cluster_seed = input_towers_temp.at(0).energy;
0090 tempstructC.cluster_NTowers = 1;
0091 tempstructC.cluster_NtrueID = 1;
0092 tempstructC.cluster_trueID = input_towers_temp.at(0).tower_trueID;
0093 cluster_towers_temp.push_back(input_towers_temp.at(0));
0094
0095
0096
0097
0098
0099 input_towers_temp.erase(input_towers_temp.begin());
0100 for (int tit = 0; tit < (int)cluster_towers_temp.size(); tit++) {
0101
0102 int iEtaTwr = cluster_towers_temp.at(tit).cellIDx;
0103 int iPhiTwr = cluster_towers_temp.at(tit).cellIDy;
0104 int iLTwr = cluster_towers_temp.at(tit).cellIDz;
0105 for (int ait = 0; ait < (int)input_towers_temp.size(); ait++) {
0106 int iEtaTwrAgg = input_towers_temp.at(ait).cellIDx;
0107 int iPhiTwrAgg = input_towers_temp.at(ait).cellIDy;
0108 int iLTwrAgg = input_towers_temp.at(ait).cellIDz;
0109
0110 int deltaL = TMath::Abs(iLTwrAgg - iLTwr);
0111 int deltaPhi = TMath::Abs(iPhiTwrAgg - iPhiTwr);
0112 int deltaEta = TMath::Abs(iEtaTwrAgg - iEtaTwr);
0113 bool neighbor = (deltaL + deltaPhi + deltaEta == 1);
0114 bool corner2D = (deltaL == 0 && deltaPhi == 1 && deltaEta == 1) ||
0115 (deltaL == 1 && deltaPhi == 0 && deltaEta == 1) ||
0116 (deltaL == 1 && deltaPhi == 1 && deltaEta == 0);
0117
0118 if (neighbor || corner2D) {
0119
0120
0121
0122 if (input_towers_temp.at(ait).energy >= (cluster_towers_temp.at(tit).energy + aggMargin))
0123 continue;
0124 tempstructC.cluster_E += input_towers_temp.at(ait).energy;
0125 tempstructC.cluster_NTowers++;
0126 cluster_towers_temp.push_back(input_towers_temp.at(ait));
0127
0128
0129
0130
0131
0132
0133 input_towers_temp.erase(input_towers_temp.begin() + ait);
0134 ait--;
0135 }
0136 }
0137 }
0138 }
0139 return tempstructC;
0140 }
0141
0142
0143 float* CalculateM02andWeightedPosition(std::vector<towersStrct> cluster_towers,
0144 float cluster_E_calc, float weight0) {
0145 static float returnVariables[8];
0146 float w_tot = 0;
0147 std::vector<float> w_i;
0148 TVector3 vecTwr;
0149 TVector3 vecTwrTmp;
0150 float w_0 = weight0;
0151
0152 vecTwr = {0., 0., 0.};
0153
0154 for (int cellI = 0; cellI < (int)cluster_towers.size(); cellI++) {
0155 w_i.push_back(TMath::Max(
0156 (float)0, (float)(w_0 + TMath::Log(cluster_towers.at(cellI).energy / cluster_E_calc))));
0157 w_tot += w_i.at(cellI);
0158 if (w_i.at(cellI) > 0) {
0159 vecTwrTmp = TVector3(cluster_towers.at(cellI).posx, cluster_towers.at(cellI).posy,
0160 cluster_towers.at(cellI).posz);
0161 vecTwr += w_i.at(cellI) * vecTwrTmp;
0162 }
0163 }
0164
0165 returnVariables[2] = vecTwr.Eta();
0166 returnVariables[3] =
0167 vecTwr.Phi();
0168 vecTwr *= 1. / w_tot;
0169
0170 returnVariables[4] = vecTwr.X();
0171 returnVariables[5] = vecTwr.Y();
0172 returnVariables[6] = vecTwr.Z();
0173
0174
0175 float delta_phi_phi[4] = {0};
0176 float delta_eta_eta[4] = {0};
0177 float delta_eta_phi[4] = {0};
0178 float dispersion = 0;
0179
0180 for (int cellI = 0; cellI < (int)cluster_towers.size(); cellI++) {
0181 int iphi = cluster_towers.at(cellI).cellIDy;
0182 int ieta = cluster_towers.at(cellI).cellIDx;
0183 delta_phi_phi[1] += (w_i.at(cellI) * iphi * iphi) / w_tot;
0184 delta_phi_phi[2] += (w_i.at(cellI) * iphi) / w_tot;
0185 delta_phi_phi[3] += (w_i.at(cellI) * iphi) / w_tot;
0186
0187 delta_eta_eta[1] += (w_i.at(cellI) * ieta * ieta) / w_tot;
0188 delta_eta_eta[2] += (w_i.at(cellI) * ieta) / w_tot;
0189 delta_eta_eta[3] += (w_i.at(cellI) * ieta) / w_tot;
0190
0191 delta_eta_phi[1] += (w_i.at(cellI) * ieta * iphi) / w_tot;
0192 delta_eta_phi[2] += (w_i.at(cellI) * iphi) / w_tot;
0193 delta_eta_phi[3] += (w_i.at(cellI) * ieta) / w_tot;
0194
0195 vecTwrTmp = TVector3(cluster_towers.at(cellI).posx, cluster_towers.at(cellI).posy,
0196 cluster_towers.at(cellI).posz);
0197
0198 vecTwr *= std::abs(vecTwrTmp.Z() / vecTwr.Z());
0199 float dx2 = pow(vecTwrTmp.X() - vecTwr.X(), 2);
0200 float dy2 = pow(vecTwrTmp.Y() - vecTwr.Y(), 2);
0201 float dz2 = pow(vecTwrTmp.Z() - vecTwr.Z(), 2);
0202 dispersion += (w_i.at(cellI) * (dx2 + dy2 + dz2)) / w_tot;
0203 }
0204 returnVariables[7] = dispersion;
0205 delta_phi_phi[0] = delta_phi_phi[1] - (delta_phi_phi[2] * delta_phi_phi[3]);
0206 delta_eta_eta[0] = delta_eta_eta[1] - (delta_eta_eta[2] * delta_eta_eta[3]);
0207 delta_eta_phi[0] = delta_eta_phi[1] - (delta_eta_phi[2] * delta_eta_phi[3]);
0208
0209 float calcM02 = 0.5 * (delta_phi_phi[0] + delta_eta_eta[0]) +
0210 TMath::Sqrt(0.25 * TMath::Power((delta_phi_phi[0] - delta_eta_eta[0]), 2) +
0211 TMath::Power(delta_eta_phi[0], 2));
0212 float calcM20 = 0.5 * (delta_phi_phi[0] + delta_eta_eta[0]) -
0213 TMath::Sqrt(0.25 * TMath::Power((delta_phi_phi[0] - delta_eta_eta[0]), 2) +
0214 TMath::Power(delta_eta_phi[0], 2));
0215
0216 returnVariables[0] = calcM02;
0217 returnVariables[1] = calcM20;
0218 return returnVariables;
0219 }