File indexing completed on 2025-08-06 08:17:47
0001 #include "TowerJetInput.h"
0002
0003 #include "Jet.h"
0004 #include "Jetv2.h"
0005
0006 #include <calobase/RawTower.h>
0007 #include <calobase/RawTowerContainer.h>
0008 #include <calobase/RawTowerDefs.h> // for encode_towerid
0009 #include <calobase/RawTowerGeom.h>
0010 #include <calobase/RawTowerGeomContainer.h>
0011 #include <calobase/TowerInfo.h>
0012 #include <calobase/TowerInfoContainer.h>
0013 #include <globalvertex/GlobalVertex.h>
0014 #include <globalvertex/GlobalVertexMap.h>
0015
0016 #include <phool/getClass.h>
0017
0018 #include <cassert>
0019 #include <cmath> // for asinh, atan2, cos, cosh
0020 #include <iostream>
0021 #include <map> // for _Rb_tree_const_iterator
0022 #include <utility> // for pair
0023 #include <vector>
0024
0025 TowerJetInput::TowerJetInput(Jet::SRC input, const std::string &prefix)
0026 : m_input(input)
0027 , m_towerNodePrefix(prefix)
0028 {
0029 }
0030
0031 void TowerJetInput::identify(std::ostream &os)
0032 {
0033 os << " TowerJetInput: ";
0034 if (m_input == Jet::CEMC_TOWER)
0035 {
0036 os << "TOWER_CEMC to Jet::CEMC_TOWER";
0037 }
0038 else if (m_input == Jet::EEMC_TOWER)
0039 {
0040 os << "TOWER_EEMC to Jet::EEMC_TOWER";
0041 }
0042 else if (m_input == Jet::HCALIN_TOWER)
0043 {
0044 os << "TOWER_HCALIN to Jet::HCALIN_TOWER";
0045 }
0046 else if (m_input == Jet::HCALOUT_TOWER)
0047 {
0048 os << "TOWER_HCALOUT to Jet::HCALOUT_TOWER";
0049 }
0050 else if (m_input == Jet::FEMC_TOWER)
0051 {
0052 os << "TOWER_FEMC to Jet::FEMC_TOWER";
0053 }
0054 else if (m_input == Jet::FHCAL_TOWER)
0055 {
0056 os << "TOWER_FHCAL to Jet::FHCAL_TOWER";
0057 }
0058 os << std::endl;
0059 }
0060
0061 std::vector<Jet *> TowerJetInput::get_input(PHCompositeNode *topNode)
0062 {
0063 if (Verbosity() > 0)
0064 {
0065 std::cout << "TowerJetInput::process_event -- entered" << std::endl;
0066 }
0067 float vtxz = 0;
0068 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0069 if (!vertexmap)
0070 {
0071 std::cout << "TowerJetInput::get_input - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0072 assert(vertexmap);
0073
0074 return std::vector<Jet *>();
0075 }
0076 if (vertexmap->empty())
0077 {
0078 if (Verbosity() > 0)
0079 {
0080 std::cout << "TowerJetInput::get_input - empty vertex map, continuing as if zvtx = 0" << std::endl;
0081 }
0082 }
0083 else
0084 {
0085 GlobalVertex *vtx = vertexmap->begin()->second;
0086 if (vtx)
0087 {
0088 if (m_use_vertextype)
0089 {
0090 auto typeStartIter = vtx->find_vertexes(m_vertex_type);
0091 auto typeEndIter = vtx->end_vertexes();
0092 for (auto iter = typeStartIter; iter != typeEndIter; ++iter)
0093 {
0094 const auto &[type, vertexVec] = *iter;
0095 if (type != m_vertex_type)
0096 {
0097 continue;
0098 }
0099 for (const auto *vertex : vertexVec)
0100 {
0101 if (!vertex)
0102 {
0103 continue;
0104 }
0105 vtxz = vertex->get_z();
0106 }
0107 }
0108 }
0109 else
0110 {
0111 vtxz = vtx->get_z();
0112 }
0113 }
0114 }
0115 if (std::isnan(vtxz))
0116 {
0117 static bool once = true;
0118 if (once)
0119 {
0120 once = false;
0121 std::cout << "TowerJetInput::get_input - WARNING - vertex is NAN. Continue with zvtx = 0 (further vertex warning will be suppressed)." << std::endl;
0122 }
0123 vtxz = 0;
0124 }
0125
0126 if (std::abs(vtxz) > 1e3)
0127 {
0128 static bool once = true;
0129 if (once)
0130 {
0131 once = false;
0132
0133 std::cout << "TowerJetInput::get_input - WARNING - vertex is " << vtxz << ". Set vtxz = 0 (further vertex warning will be suppressed)." << std::endl;
0134 }
0135 vtxz = 0;
0136 }
0137
0138 m_use_towerinfo = false;
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163 RawTowerContainer *towers = nullptr;
0164 TowerInfoContainer *towerinfos = nullptr;
0165 RawTowerGeomContainer *geom = nullptr;
0166 RawTowerGeomContainer *EMCal_geom = nullptr;
0167
0168 if (m_input == Jet::CEMC_TOWER)
0169 {
0170 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0171 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0172 if ((!towers) || !geom)
0173 {
0174 return std::vector<Jet *>();
0175 }
0176 }
0177 else if (m_input == Jet::CEMC_TOWERINFO)
0178 {
0179 m_use_towerinfo = true;
0180 towerName = m_towerNodePrefix + "_CEMC";
0181 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0182 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0183 geocaloid = RawTowerDefs::CalorimeterId::CEMC;
0184 if ((!towerinfos) || !geom)
0185 {
0186 return std::vector<Jet *>();
0187 }
0188 }
0189 else if (m_input == Jet::CEMC_TOWERINFO_EMBED)
0190 {
0191 m_use_towerinfo = true;
0192 towerName = m_towerNodePrefix + "_EMBED_CEMC";
0193 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0194 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0195 geocaloid = RawTowerDefs::CalorimeterId::CEMC;
0196 if ((!towerinfos) || !geom)
0197 {
0198 return std::vector<Jet *>();
0199 }
0200 }
0201 else if (m_input == Jet::CEMC_TOWERINFO_SIM)
0202 {
0203 m_use_towerinfo = true;
0204 towerName = m_towerNodePrefix + "_SIM_CEMC";
0205 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0206 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0207 geocaloid = RawTowerDefs::CalorimeterId::CEMC;
0208 if ((!towerinfos) || !geom)
0209 {
0210 return std::vector<Jet *>();
0211 }
0212 }
0213 else if (m_input == Jet::EEMC_TOWER)
0214 {
0215 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_EEMC");
0216 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_EEMC");
0217 if ((!towers && !towerinfos) || !geom)
0218 {
0219 return std::vector<Jet *>();
0220 }
0221 }
0222 else if (m_input == Jet::HCALIN_TOWER)
0223 {
0224 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0225 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0226 if ((!towers) || !geom)
0227 {
0228 return std::vector<Jet *>();
0229 }
0230 }
0231 else if (m_input == Jet::HCALIN_TOWERINFO)
0232 {
0233 m_use_towerinfo = true;
0234 towerName = m_towerNodePrefix + "_HCALIN";
0235 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0236 geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0237 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0238 if ((!towerinfos) || !geom)
0239 {
0240 return std::vector<Jet *>();
0241 }
0242 }
0243 else if (m_input == Jet::HCALIN_TOWERINFO_EMBED)
0244 {
0245 m_use_towerinfo = true;
0246 towerName = m_towerNodePrefix + "_EMBED_HCALIN";
0247 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0248 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0249 geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0250 if ((!towerinfos) || !geom)
0251 {
0252 return std::vector<Jet *>();
0253 }
0254 }
0255 else if (m_input == Jet::HCALIN_TOWERINFO_SIM)
0256 {
0257 m_use_towerinfo = true;
0258 towerName = m_towerNodePrefix + "_SIM_HCALIN";
0259 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0260 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0261 geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0262 if ((!towerinfos) || !geom)
0263 {
0264 return std::vector<Jet *>();
0265 }
0266 }
0267 else if (m_input == Jet::HCALOUT_TOWER)
0268 {
0269 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0270 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0271 if ((!towers) || !geom)
0272 {
0273 return std::vector<Jet *>();
0274 }
0275 }
0276 else if (m_input == Jet::HCALOUT_TOWERINFO)
0277 {
0278 m_use_towerinfo = true;
0279 towerName = m_towerNodePrefix + "_HCALOUT";
0280 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0281 geocaloid = RawTowerDefs::CalorimeterId::HCALOUT;
0282 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0283 if ((!towerinfos) || !geom)
0284 {
0285 return std::vector<Jet *>();
0286 }
0287 }
0288 else if (m_input == Jet::HCALOUT_TOWERINFO_EMBED)
0289 {
0290 m_use_towerinfo = true;
0291 towerName = m_towerNodePrefix + "_EMBED_HCALOUT";
0292 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0293 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0294 geocaloid = RawTowerDefs::CalorimeterId::HCALOUT;
0295 if ((!towerinfos) || !geom)
0296 {
0297 return std::vector<Jet *>();
0298 }
0299 }
0300 else if (m_input == Jet::HCALOUT_TOWERINFO_SIM)
0301 {
0302 m_use_towerinfo = true;
0303 towerName = m_towerNodePrefix + "_SIM_HCALOUT";
0304 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0305 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0306 geocaloid = RawTowerDefs::CalorimeterId::HCALOUT;
0307 if ((!towerinfos) || !geom)
0308 {
0309 return std::vector<Jet *>();
0310 }
0311 }
0312
0313 else if (m_input == Jet::FEMC_TOWER)
0314 {
0315 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_FEMC");
0316 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_FEMC");
0317 if ((!towers) || !geom)
0318 {
0319 return std::vector<Jet *>();
0320 }
0321 }
0322 else if (m_input == Jet::FHCAL_TOWER)
0323 {
0324 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_FHCAL");
0325 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_FHCAL");
0326 if ((!towers) || !geom)
0327 {
0328 return std::vector<Jet *>();
0329 }
0330 }
0331 else if (m_input == Jet::CEMC_TOWER_RETOWER)
0332 {
0333 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0334 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0335 if ((!towers) || !geom)
0336 {
0337 return std::vector<Jet *>();
0338 }
0339 }
0340 else if (m_input == Jet::CEMC_TOWERINFO_RETOWER)
0341 {
0342 m_use_towerinfo = true;
0343 towerName = m_towerNodePrefix + "_CEMC_RETOWER";
0344 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0345 geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0346 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0347 if ((!towerinfos) || !geom)
0348 {
0349 return std::vector<Jet *>();
0350 }
0351 }
0352 else if (m_input == Jet::CEMC_TOWER_SUB1)
0353 {
0354 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
0355 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0356 if ((!towers) || !geom)
0357 {
0358 return std::vector<Jet *>();
0359 }
0360 }
0361 else if (m_input == Jet::CEMC_TOWERINFO_SUB1)
0362 {
0363 m_use_towerinfo = true;
0364 towerName = m_towerNodePrefix + "_CEMC_RETOWER_SUB1";
0365 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0366 geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0367 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0368 if ((!towerinfos) || !geom)
0369 {
0370 return std::vector<Jet *>();
0371 }
0372 }
0373 else if (m_input == Jet::HCALIN_TOWER_SUB1)
0374 {
0375 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
0376 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0377 if ((!towers) || !geom)
0378 {
0379 return std::vector<Jet *>();
0380 }
0381 }
0382 else if (m_input == Jet::HCALIN_TOWERINFO_SUB1)
0383 {
0384 m_use_towerinfo = true;
0385 towerName = m_towerNodePrefix + "_HCALIN_SUB1";
0386 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0387 geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0388 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0389 if ((!towerinfos) || !geom)
0390 {
0391 return std::vector<Jet *>();
0392 }
0393 }
0394 else if (m_input == Jet::HCALOUT_TOWER_SUB1)
0395 {
0396 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
0397 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0398 if ((!towers) || !geom)
0399 {
0400 return std::vector<Jet *>();
0401 }
0402 }
0403 else if (m_input == Jet::HCALOUT_TOWERINFO_SUB1)
0404 {
0405 m_use_towerinfo = true;
0406 towerName = m_towerNodePrefix + "_HCALOUT_SUB1";
0407 towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerName);
0408 geocaloid = RawTowerDefs::CalorimeterId::HCALOUT;
0409 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0410 if ((!towerinfos) || !geom)
0411 {
0412 return std::vector<Jet *>();
0413 }
0414 }
0415 else if (m_input == Jet::CEMC_TOWER_SUB1CS)
0416 {
0417 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1CS");
0418 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0419 if ((!towers) || !geom)
0420 {
0421 return std::vector<Jet *>();
0422 }
0423 }
0424 else if (m_input == Jet::HCALIN_TOWER_SUB1CS)
0425 {
0426 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1CS");
0427 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0428 if ((!towers) || !geom)
0429 {
0430 return std::vector<Jet *>();
0431 }
0432 }
0433 else if (m_input == Jet::HCALOUT_TOWER_SUB1CS)
0434 {
0435 towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1CS");
0436 geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0437 if ((!towers) || !geom)
0438 {
0439 return std::vector<Jet *>();
0440 }
0441 }
0442 else
0443 {
0444 return std::vector<Jet *>();
0445 }
0446
0447
0448 if (m_input == Jet::CEMC_TOWER_RETOWER || m_input == Jet::CEMC_TOWERINFO_RETOWER || m_input == Jet::CEMC_TOWER_SUB1 || m_input == Jet::CEMC_TOWERINFO_SUB1 || m_input == Jet::CEMC_TOWER_SUB1CS)
0449 {
0450 EMCal_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0451 if (!EMCal_geom)
0452 {
0453 return std::vector<Jet *>();
0454 }
0455 }
0456
0457
0458
0459 std::vector<Jet *> pseudojets;
0460 if (m_use_towerinfo)
0461 {
0462 if (!towerinfos)
0463 {
0464 return std::vector<Jet *>();
0465 }
0466
0467 unsigned int nchannels = towerinfos->size();
0468 for (unsigned int channel = 0; channel < nchannels; channel++)
0469 {
0470 TowerInfo *tower = towerinfos->get_tower_at_channel(channel);
0471 assert(tower);
0472
0473 unsigned int calokey = towerinfos->encode_key(channel);
0474 int ieta = towerinfos->getTowerEtaBin(calokey);
0475 int iphi = towerinfos->getTowerPhiBin(calokey);
0476 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(geocaloid, ieta, iphi);
0477
0478 if (tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2())
0479 {
0480 continue;
0481 }
0482 if (std::isnan(tower->get_energy()))
0483 {
0484 continue;
0485 }
0486 RawTowerGeom *tower_geom = geom->get_tower_geometry(key);
0487 assert(tower_geom);
0488
0489 double r = tower_geom->get_center_radius();
0490 if (m_input == Jet::CEMC_TOWER_RETOWER || m_input == Jet::CEMC_TOWERINFO_RETOWER || m_input == Jet::CEMC_TOWER_SUB1 || m_input == Jet::CEMC_TOWERINFO_SUB1 || m_input == Jet::CEMC_TOWER_SUB1CS)
0491 {
0492 const RawTowerDefs::keytype EMCal_key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, 0, 0);
0493 RawTowerGeom *EMCal_tower_geom = EMCal_geom->get_tower_geometry(EMCal_key);
0494 assert(EMCal_tower_geom);
0495 r = EMCal_tower_geom->get_center_radius();
0496 }
0497 double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0498 double towereta = tower_geom->get_eta();
0499 double z0 = sinh(towereta) * r;
0500 double z = z0 - vtxz;
0501 double eta = asinh(z / r);
0502 double pt = tower->get_energy() / cosh(eta);
0503 double e = tower->get_energy();
0504 double px = pt * cos(phi);
0505 double py = pt * sin(phi);
0506 double pz = pt * sinh(eta);
0507
0508 Jet *jet = new Jetv2();
0509 jet->set_px(px);
0510 jet->set_py(py);
0511 jet->set_pz(pz);
0512 jet->set_e(e);
0513 jet->insert_comp(m_input, channel);
0514 pseudojets.push_back(jet);
0515 }
0516 }
0517 else
0518 {
0519 RawTowerContainer::ConstRange begin_end = towers->getTowers();
0520 RawTowerContainer::ConstIterator rtiter;
0521 for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0522 {
0523 RawTower *tower = rtiter->second;
0524
0525 RawTowerGeom *tower_geom = geom->get_tower_geometry(tower->get_key());
0526 assert(tower_geom);
0527
0528 double r = tower_geom->get_center_radius();
0529 if (m_input == Jet::CEMC_TOWER_RETOWER || m_input == Jet::CEMC_TOWERINFO_RETOWER || m_input == Jet::CEMC_TOWER_SUB1 || m_input == Jet::CEMC_TOWERINFO_SUB1 || m_input == Jet::CEMC_TOWER_SUB1CS)
0530 {
0531 const RawTowerDefs::keytype EMCal_key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, 0, 0);
0532 RawTowerGeom *EMCal_tower_geom = EMCal_geom->get_tower_geometry(EMCal_key);
0533 assert(EMCal_tower_geom);
0534 r = EMCal_tower_geom->get_center_radius();
0535 }
0536 double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0537 double towereta = tower_geom->get_eta();
0538 double z0 = sinh(towereta) * r;
0539 double z = z0 - vtxz;
0540 double eta = asinh(z / r);
0541 double pt = tower->get_energy() / cosh(eta);
0542 double px = pt * cos(phi);
0543 double py = pt * sin(phi);
0544 double pz = pt * sinh(eta);
0545
0546 Jet *jet = new Jetv2();
0547 jet->set_px(px);
0548 jet->set_py(py);
0549 jet->set_pz(pz);
0550 jet->set_e(tower->get_energy());
0551 jet->insert_comp(m_input, tower->get_id());
0552 pseudojets.push_back(jet);
0553 }
0554 }
0555 if (Verbosity() > 0)
0556 {
0557 std::cout << "TowerJetInput::process_event -- exited" << std::endl;
0558 }
0559 return pseudojets;
0560 }
