packages/CaloBase/RawClusterv1.cc

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