File indexing completed on 2025-12-17 09:19:52
0001 #include "RawClusterv1.h"
0002
0003 #include "RawTowerDefs.h"
0004
0005 #include <phool/phool.h>
0006
0007 #include <cstdlib>
0008 #include <limits>
0009 #include <string>
0010 #include <vector>
0011
0012 RawClusterv1::RawClusterv1(const RawCluster& cluster)
0013 {
0014 set_id(cluster.get_id());
0015 set_energy(cluster.get_energy());
0016 set_r(cluster.get_r());
0017 set_phi(cluster.get_phi());
0018 set_z(cluster.get_z());
0019
0020 for (const auto& tower : cluster.get_towermap())
0021 {
0022 addTower(tower.first, tower.second);
0023 }
0024
0025 static const PROPERTY kKnownProperties[] = {
0026 prop_ecore,
0027 prop_prob,
0028 prop_chi2,
0029 prop_merged_cluster_prob,
0030 prop_et_iso_calotower_sub_R01,
0031 prop_et_iso_calotower_R01,
0032 prop_et_iso_calotower_sub_R02,
0033 prop_et_iso_calotower_R02,
0034 prop_et_iso_calotower_sub_R03,
0035 prop_et_iso_calotower_R03,
0036 prop_et_iso_calotower_sub_R04,
0037 prop_et_iso_calotower_R04,
0038
0039
0040 prop_tower_x_raw,
0041 prop_tower_y_raw,
0042 prop_tower_x_corr,
0043 prop_tower_y_corr,
0044 prop_tower_t_mean,
0045
0046
0047 prop_incidence_alpha_phi,
0048 prop_incidence_alpha_eta};
0049
0050 for (const auto prop_id : kKnownProperties)
0051 {
0052 copy_property_from_cluster(cluster, prop_id);
0053 }
0054 }
0055
0056 void RawClusterv1::copy_property_from_cluster(const RawCluster& source, const PROPERTY prop_id)
0057 {
0058 if (!source.has_property(prop_id))
0059 {
0060 return;
0061 }
0062
0063 std::pair<const std::string, PROPERTY_TYPE> property_info = RawCluster::get_property_info(prop_id);
0064 switch (property_info.second)
0065 {
0066 case type_int:
0067 set_property(prop_id, source.get_property_int(prop_id));
0068 break;
0069 case type_uint:
0070 set_property(prop_id, source.get_property_uint(prop_id));
0071 break;
0072 case type_float:
0073 set_property(prop_id, source.get_property_float(prop_id));
0074 break;
0075 default:
0076 break;
0077 }
0078 }
0079
0080 void RawClusterv1::Reset()
0081 {
0082 clusterid = 0;
0083 _z = std::numeric_limits<float>::quiet_NaN();
0084 _r = std::numeric_limits<float>::quiet_NaN();
0085 _phi = std::numeric_limits<float>::quiet_NaN();
0086 _energy = std::numeric_limits<float>::quiet_NaN();
0087 prop_map.clear();
0088 towermap.clear();
0089 }
0090
0091 void RawClusterv1::addTower(const RawClusterDefs::keytype twrid, const float etower)
0092 {
0093 if (towermap.contains(twrid))
0094 {
0095 std::cout << "tower 0x" << std::hex << twrid << ", dec: " << std::dec
0096 << twrid << " already exists, that is bad" << std::endl;
0097 exit(1);
0098 }
0099 towermap[twrid] = etower;
0100 }
0101
0102 void RawClusterv1::identify(std::ostream& os) const
0103 {
0104 os << "RawClusterv1 ID " << get_id() << " consist of " << getNTowers() << " towers with total energy of " << get_energy() << " GeV ";
0105 os << "@ (r,phi,z) = (" << get_r() << ", " << get_phi() << ", " << get_z() << "), (x,y,z) = (" << get_x() << ", " << get_y() << ", " << get_z() << ")";
0106
0107 for (auto i : prop_map)
0108 {
0109 PROPERTY prop_id = static_cast<PROPERTY>(i.first);
0110 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0111 os << "\t" << prop_id << ":\t" << property_info.first << " = \t";
0112 switch (property_info.second)
0113 {
0114 case type_int:
0115 os << get_property_int(prop_id);
0116 break;
0117 case type_uint:
0118 os << get_property_uint(prop_id);
0119 break;
0120 case type_float:
0121 os << get_property_float(prop_id);
0122 break;
0123 default:
0124 os << " unknown type ";
0125 break;
0126 }
0127 os << std::endl;
0128 }
0129 }
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145 bool RawClusterv1::has_property(const PROPERTY prop_id) const
0146 {
0147 prop_map_t::const_iterator i = prop_map.find(prop_id);
0148 return i != prop_map.end();
0149 }
0150
0151 float RawClusterv1::get_property_float(const PROPERTY prop_id) const
0152 {
0153 if (!check_property(prop_id, type_float))
0154 {
0155 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0156 std::cout << PHWHERE << " Property " << property_info.first << " with id "
0157 << prop_id << " is of type " << get_property_type(property_info.second)
0158 << " not " << get_property_type(type_float) << std::endl;
0159 exit(1);
0160 }
0161 prop_map_t::const_iterator i = prop_map.find(prop_id);
0162
0163 if (i != prop_map.end())
0164 {
0165 return u_property(i->second).fdata;
0166 }
0167
0168 return std::numeric_limits<float>::quiet_NaN();
0169 }
0170
0171 int RawClusterv1::get_property_int(const PROPERTY prop_id) const
0172 {
0173 if (!check_property(prop_id, type_int))
0174 {
0175 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0176 std::cout << PHWHERE << " Property " << property_info.first << " with id "
0177 << prop_id << " is of type " << get_property_type(property_info.second)
0178 << " not " << get_property_type(type_int) << std::endl;
0179 exit(1);
0180 }
0181 prop_map_t::const_iterator i = prop_map.find(prop_id);
0182
0183 if (i != prop_map.end())
0184 {
0185 return u_property(i->second).idata;
0186 }
0187
0188 return std::numeric_limits<int>::min();
0189 }
0190
0191 unsigned int
0192 RawClusterv1::get_property_uint(const PROPERTY prop_id) const
0193 {
0194 if (!check_property(prop_id, type_uint))
0195 {
0196 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0197 std::cout << PHWHERE << " Property " << property_info.first << " with id "
0198 << prop_id << " is of type " << get_property_type(property_info.second)
0199 << " not " << get_property_type(type_uint) << std::endl;
0200 exit(1);
0201 }
0202 prop_map_t::const_iterator i = prop_map.find(prop_id);
0203
0204 if (i != prop_map.end())
0205 {
0206 return u_property(i->second).uidata;
0207 }
0208
0209 return std::numeric_limits<unsigned int>::max();
0210 }
0211
0212 void RawClusterv1::set_property(const PROPERTY prop_id, const float value)
0213 {
0214 if (!check_property(prop_id, type_float))
0215 {
0216 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0217 std::cout << PHWHERE << " Property " << property_info.first << " with id "
0218 << prop_id << " is of type " << get_property_type(property_info.second)
0219 << " not " << get_property_type(type_float) << std::endl;
0220 exit(1);
0221 }
0222 prop_map[prop_id] = u_property(value).get_storage();
0223 }
0224
0225 void RawClusterv1::set_property(const PROPERTY prop_id, const int value)
0226 {
0227 if (!check_property(prop_id, type_int))
0228 {
0229 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0230 std::cout << PHWHERE << " Property " << property_info.first << " with id "
0231 << prop_id << " is of type " << get_property_type(property_info.second)
0232 << " not " << get_property_type(type_int) << std::endl;
0233 exit(1);
0234 }
0235 prop_map[prop_id] = u_property(value).get_storage();
0236 }
0237
0238 void RawClusterv1::set_property(const PROPERTY prop_id, const unsigned int value)
0239 {
0240 if (!check_property(prop_id, type_uint))
0241 {
0242 std::pair<const std::string, PROPERTY_TYPE> property_info = get_property_info(prop_id);
0243 std::cout << PHWHERE << " Property " << property_info.first << " with id "
0244 << prop_id << " is of type " << get_property_type(property_info.second)
0245 << " not " << get_property_type(type_uint) << std::endl;
0246 exit(1);
0247 }
0248 prop_map[prop_id] = u_property(value).get_storage();
0249 }
0250 void RawClusterv1::set_et_iso(const float et_iso, const int radiusx10, bool subtracted, bool clusterTower = true)
0251 {
0252 if (clusterTower)
0253 {
0254 if (subtracted)
0255 {
0256 switch (radiusx10)
0257 {
0258 case 1:
0259 set_property(prop_et_iso_calotower_sub_R01, et_iso);
0260 break;
0261 case 2:
0262 set_property(prop_et_iso_calotower_sub_R02, et_iso);
0263 break;
0264 case 3:
0265 set_property(prop_et_iso_calotower_sub_R03, et_iso);
0266 break;
0267 case 4:
0268 set_property(prop_et_iso_calotower_sub_R04, et_iso);
0269 break;
0270 default:
0271 std::string warning = "set_et_iso(const int radiusx10, bool subtracted, bool clusterTower) - radius:" + std::to_string(radiusx10) + " has not been defined";
0272 PHOOL_VIRTUAL_WARN(warning.c_str());
0273 break;
0274 }
0275 }
0276 else
0277 {
0278 switch (radiusx10)
0279 {
0280 case 1:
0281 set_property(prop_et_iso_calotower_R01, et_iso);
0282 break;
0283 case 2:
0284 set_property(prop_et_iso_calotower_R02, et_iso);
0285 break;
0286 case 3:
0287 set_property(prop_et_iso_calotower_R03, et_iso);
0288 break;
0289 case 4:
0290 set_property(prop_et_iso_calotower_R04, et_iso);
0291 break;
0292 default:
0293 std::string warning = "set_et_iso(const int radiusx10, bool subtracted, bool clusterTower) - radius:" + std::to_string(radiusx10) + " has not been defined";
0294 PHOOL_VIRTUAL_WARN(warning.c_str());
0295 break;
0296 }
0297 }
0298 }
0299 else
0300 {
0301 PHOOL_VIRTUAL_WARN("set_et_iso(const int radiusx10, bool subtracted, bool clusterTower) - nonclusterTower algorithms have not been defined");
0302 }
0303 }
0304
0305 unsigned int
0306 RawClusterv1::get_property_nocheck(const PROPERTY prop_id) const
0307 {
0308 prop_map_t::const_iterator iter = prop_map.find(prop_id);
0309 if (iter != prop_map.end())
0310 {
0311 return iter->second;
0312 }
0313 return std::numeric_limits<unsigned int>::max();
0314 }
0315
0316 float RawClusterv1::get_et_iso(const int radiusx10 = 3, bool subtracted = false, bool clusterTower = true) const
0317 {
0318 float r;
0319 if (clusterTower)
0320 {
0321 if (subtracted)
0322 {
0323 switch (radiusx10)
0324 {
0325 case 1:
0326 r = get_property_float(prop_et_iso_calotower_sub_R01);
0327 break;
0328 case 2:
0329 r = get_property_float(prop_et_iso_calotower_sub_R02);
0330 break;
0331 case 3:
0332 r = get_property_float(prop_et_iso_calotower_sub_R03);
0333 break;
0334 case 4:
0335 r = get_property_float(prop_et_iso_calotower_sub_R04);
0336 break;
0337 default:
0338 std::string warning = "get_et_iso(const int radiusx10, bool subtracted, bool clusterTower) - radius:" + std::to_string(radiusx10) + " has not been defined";
0339 PHOOL_VIRTUAL_WARN(warning.c_str());
0340 r = std::numeric_limits<float>::quiet_NaN();
0341 break;
0342 }
0343 }
0344 else
0345 {
0346 switch (radiusx10)
0347 {
0348 case 1:
0349 r = get_property_float(prop_et_iso_calotower_R01);
0350 break;
0351 case 2:
0352 r = get_property_float(prop_et_iso_calotower_R02);
0353 break;
0354 case 3:
0355 r = get_property_float(prop_et_iso_calotower_R03);
0356 break;
0357 case 4:
0358 r = get_property_float(prop_et_iso_calotower_R04);
0359 break;
0360 default:
0361 std::string warning = "get_et_iso(const int radiusx10, bool subtracted, bool clusterTower) - radius:" + std::to_string(radiusx10) + " has not been defined";
0362 PHOOL_VIRTUAL_WARN(warning.c_str());
0363 r = std::numeric_limits<float>::quiet_NaN();
0364 break;
0365 }
0366 }
0367 }
0368 else
0369 {
0370 PHOOL_VIRTUAL_WARN("get_et_iso(const int radiusx10, bool subtracted, bool clusterTower) - nonclusterTower algorithms have not been defined");
0371 r = std::numeric_limits<float>::quiet_NaN();
0372 }
0373 return r;
0374 }
0375
0376 std::vector<float> RawClusterv1::get_shower_shapes(float tower_thresh) const
0377 {
0378
0379 RawTowerDefs::keytype maxtowerkey = 0;
0380 float maxtowerE = 0;
0381 for (const auto& iter : towermap)
0382 {
0383 if (iter.second > maxtowerE)
0384 {
0385 maxtowerE = iter.second;
0386 maxtowerkey = iter.first;
0387 }
0388 }
0389
0390
0391 if (towermap.empty() || maxtowerE == 0)
0392 {
0393 std::vector<float> v;
0394 return v;
0395 }
0396 int maxtowerieta = RawTowerDefs::decode_index1(maxtowerkey);
0397 int maxtoweriphi = RawTowerDefs::decode_index2(maxtowerkey);
0398
0399 float deta1 = 0;
0400 float dphi1 = 0;
0401
0402 float deta2 = 0;
0403 float dphi2 = 0;
0404
0405 float totalE = 0;
0406
0407
0408 int totalphibins = 256;
0409 auto dphiwrap = [totalphibins](float towerphi, float maxiphi)
0410 {
0411 float idphi = towerphi - maxiphi;
0412 float idphiwrap = totalphibins - std::abs(idphi);
0413 if (std::abs(idphiwrap) < std::abs(idphi))
0414 {
0415 return (idphi > 0) ? -idphiwrap : idphiwrap;
0416 }
0417 return idphi;
0418 };
0419 for (const auto& iter : towermap)
0420 {
0421 RawTowerDefs::keytype towerkey = iter.first;
0422 float E = iter.second;
0423 float eta = RawTowerDefs::decode_index1(towerkey);
0424 float phi = RawTowerDefs::decode_index2(towerkey);
0425 float deta = eta - maxtowerieta;
0426
0427 float dphi = dphiwrap(phi, maxtoweriphi);
0428 if (E > tower_thresh)
0429 {
0430 totalE += E;
0431 deta1 += E * deta;
0432 dphi1 += E * dphi;
0433 deta2 += E * deta * deta;
0434 dphi2 += E * dphi * dphi;
0435 }
0436 }
0437 deta1 /= totalE;
0438 dphi1 /= totalE;
0439 deta2 = deta2 / totalE - (deta1 * deta1);
0440 dphi2 = dphi2 / totalE - (dphi1 * dphi1);
0441
0442 int centertowerieta = std::floor(deta1 + 0.5);
0443 int centertoweriphi = std::floor(dphi1 + 0.5);
0444 int etashift = (deta1 - centertowerieta) < 0 ? -1 : 1;
0445 int phishift = (dphi1 - centertoweriphi) < 0 ? -1 : 1;
0446
0447 auto wraptowerphi = [totalphibins](int phi) -> int
0448 {
0449 if (phi < 0)
0450 {
0451 return phi + totalphibins;
0452 }
0453 if (phi >= totalphibins)
0454 {
0455 return phi - totalphibins;
0456 }
0457 return phi;
0458 };
0459
0460 int eta1 = centertowerieta + maxtowerieta;
0461 int eta2 = centertowerieta + etashift + maxtowerieta;
0462 int phi1 = wraptowerphi(centertoweriphi + maxtoweriphi);
0463 int phi2 = wraptowerphi(centertoweriphi + phishift + maxtoweriphi);
0464
0465 auto gettowerenergy = [totalphibins, tower_thresh, this](int ieta, int phi) -> float
0466 {
0467 if (phi < 0)
0468 {
0469 return phi + totalphibins;
0470 }
0471 if (phi >= totalphibins)
0472 {
0473 return phi - totalphibins;
0474 }
0475 RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, phi);
0476 if (towermap.contains(key))
0477 {
0478 if (towermap.at(key) > tower_thresh)
0479 {
0480 return towermap.at(key);
0481 }
0482 }
0483 return 0;
0484 };
0485
0486 float e1 = gettowerenergy(eta1, phi1);
0487 float e2 = gettowerenergy(eta1, phi2);
0488 float e3 = gettowerenergy(eta2, phi2);
0489 float e4 = gettowerenergy(eta2, phi1);
0490
0491 float e1t = (e1 + e2 + e3 + e4) / totalE;
0492 float e2t = (e1 + e2 - e3 - e4) / totalE;
0493 float e3t = (e1 - e2 - e3 + e4) / totalE;
0494 float e4t = (e3) / totalE;
0495
0496 std::vector<float> v = {e1t, e2t, e3t, e4t, deta1 + maxtowerieta, dphi1 + maxtoweriphi, deta2, dphi2, e1, e2, e3, e4, totalE};
0497 return v;
0498 }
0499
0500 std::pair<int, int> RawClusterv1::get_lead_tower() const
0501 {
0502 if (towermap.empty())
0503 {
0504 return {0, 0};
0505 }
0506
0507 int phi_max = 0;
0508 int eta_max = 0;
0509 float e_max = std::numeric_limits<float>::lowest();
0510
0511 for (const auto& iter : towermap)
0512 {
0513 RawTowerDefs::keytype towerkey = iter.first;
0514 float E = iter.second;
0515 float eta = RawTowerDefs::decode_index1(towerkey);
0516 float phi = RawTowerDefs::decode_index2(towerkey);
0517 if (E > e_max)
0518 {
0519 phi_max = phi;
0520 eta_max = eta;
0521 e_max = E;
0522 }
0523 }
0524 return {eta_max, phi_max};
0525 }
