diff --git a/offline/packages/CaloReco/CaloTowerStatus.cc b/offline/packages/CaloReco/CaloTowerStatus.cc index a98a64a27b..86eb680e3e 100644 --- a/offline/packages/CaloReco/CaloTowerStatus.cc +++ b/offline/packages/CaloReco/CaloTowerStatus.cc @@ -3,10 +3,6 @@ #include // for TowerInfo #include -#include -#include -#include -#include #include // for CDBTTree @@ -154,7 +150,7 @@ int CaloTowerStatus::InitRun(PHCompositeNode *topNode) } m_calibName_hotMap = m_detector + "nome"; - if (m_dettype == CaloTowerDefs::CEMC) + if (m_dettype == CaloTowerDefs::CEMC || m_dettype == CaloTowerDefs::SEPD) { m_calibName_hotMap = m_detector + "_BadTowerMap"; } diff --git a/offline/packages/eventplaneinfo/EventPlaneRecov2.cc b/offline/packages/eventplaneinfo/EventPlaneRecov2.cc new file mode 100644 index 0000000000..9635e0f3d6 --- /dev/null +++ b/offline/packages/eventplaneinfo/EventPlaneRecov2.cc @@ -0,0 +1,637 @@ +#include "EventPlaneRecov2.h" + +#include "EventplaneinfoMapv1.h" +#include "Eventplaneinfov2.h" + +#include +#include +#include + +// -- Centrality +#include + +// -- sEPD +#include + +#include // for CDBTTree + +// -- event +#include + +#include + +#include + +#include +#include +#include + +// -- root includes -- +#include +#include + +// c++ includes -- +#include +#include +#include +#include +#include +#include + +//____________________________________________________________________________.. +EventPlaneRecov2::EventPlaneRecov2(const std::string &name): + SubsysReco(name) +{ +} + +//____________________________________________________________________________.. +int EventPlaneRecov2::Init(PHCompositeNode *topNode) +{ + std::string calibdir = CDBInterface::instance()->getUrl(m_calibName); + + if (!m_directURL_EventPlaneCalib.empty()) + { + m_cdbttree = new CDBTTree(m_directURL_EventPlaneCalib); + std::cout << PHWHERE << " Custom Event Plane Calib Found: " << m_directURL_EventPlaneCalib << std::endl; + } + else if (!calibdir.empty()) + { + m_cdbttree = new CDBTTree(calibdir); + std::cout << PHWHERE << " Event Plane Calib Found: " << calibdir << std::endl; + } + else if (m_doAbortNoEventPlaneCalib) + { + std::cout << PHWHERE << " Error: No Event Plane Calib Found and m_doAbortNoEventPlaneCalib is true. Aborting." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + else + { + std::cout << PHWHERE << " Error: No Event Plane Calib Found. Skipping Event Plane Calibrations." << std::endl; + m_doNotCalib = true; + } + + if (!m_doNotCalib) + { + LoadCalib(); + } + + if (Verbosity() > 0) + { + print_correction_data(); + } + + CreateNodes(topNode); + + return Fun4AllReturnCodes::EVENT_OK; +} + +std::array, 2> EventPlaneRecov2::calculate_flattening_matrix(double xx, double yy, double xy, int n, int cent_bin, const std::string& det_label) +{ + std::array, 2> mat{}; + + double D_arg = (xx * yy) - (xy * xy); + if (D_arg <= 0) + { + std::cout << PHWHERE << "Invalid D-term " << D_arg << " for n=" << n << ", cent bin=" << cent_bin << ", det=" << det_label << std::endl; + // Return Identity Matrix to preserve Recentered vector + mat[0][0] = 1.0; + mat[1][1] = 1.0; + return mat; + } + double D = std::sqrt(D_arg); + + double N_term = D * (xx + yy + (2 * D)); + if (N_term <= 0) + { + std::cout << PHWHERE << "Invalid N-term " << N_term << " for n=" << n << ", cent bin=" << cent_bin << ", det=" << det_label << std::endl; + // Return Identity Matrix to preserve Recentered vector + mat[0][0] = 1.0; + mat[1][1] = 1.0; + return mat; + } + double inv_sqrt_N = 1.0 / std::sqrt(N_term); + + mat[0][0] = inv_sqrt_N * (yy + D); + mat[0][1] = -inv_sqrt_N * xy; + mat[1][0] = mat[0][1]; + mat[1][1] = inv_sqrt_N * (xx + D); + return mat; +} + +//____________________________________________________________________________.. +void EventPlaneRecov2::LoadCalib() +{ + size_t south_idx = static_cast(Subdetector::S); + size_t north_idx = static_cast(Subdetector::N); + size_t ns_idx = static_cast(Subdetector::NS); + + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + int n = m_harmonics[h_idx]; + + std::string S_x_avg_name = std::format("Q_S_x_{}_avg", n); + std::string S_y_avg_name = std::format("Q_S_y_{}_avg", n); + std::string N_x_avg_name = std::format("Q_N_x_{}_avg", n); + std::string N_y_avg_name = std::format("Q_N_y_{}_avg", n); + + std::string S_xx_avg_name = std::format("Q_S_xx_{}_avg", n); + std::string S_yy_avg_name = std::format("Q_S_yy_{}_avg", n); + std::string S_xy_avg_name = std::format("Q_S_xy_{}_avg", n); + std::string N_xx_avg_name = std::format("Q_N_xx_{}_avg", n); + std::string N_yy_avg_name = std::format("Q_N_yy_{}_avg", n); + std::string N_xy_avg_name = std::format("Q_N_xy_{}_avg", n); + + std::string NS_xx_avg_name = std::format("Q_NS_xx_{}_avg", n); + std::string NS_yy_avg_name = std::format("Q_NS_yy_{}_avg", n); + std::string NS_xy_avg_name = std::format("Q_NS_xy_{}_avg", n); + + for (size_t cent_bin = 0; cent_bin < m_bins_cent; ++cent_bin) + { + int key = cent_bin; + + // South + auto& dataS = m_correction_data[h_idx][cent_bin][south_idx]; + dataS.avg_Q.x = m_cdbttree->GetDoubleValue(key, S_x_avg_name); + dataS.avg_Q.y = m_cdbttree->GetDoubleValue(key, S_y_avg_name); + + dataS.avg_Q_xx = m_cdbttree->GetDoubleValue(key, S_xx_avg_name); + dataS.avg_Q_yy = m_cdbttree->GetDoubleValue(key, S_yy_avg_name); + dataS.avg_Q_xy = m_cdbttree->GetDoubleValue(key, S_xy_avg_name); + + dataS.X_matrix = calculate_flattening_matrix(dataS.avg_Q_xx, dataS.avg_Q_yy, dataS.avg_Q_xy, n, cent_bin, "South"); + + // North + auto& dataN = m_correction_data[h_idx][cent_bin][north_idx]; + dataN.avg_Q.x = m_cdbttree->GetDoubleValue(key, N_x_avg_name); + dataN.avg_Q.y = m_cdbttree->GetDoubleValue(key, N_y_avg_name); + + dataN.avg_Q_xx = m_cdbttree->GetDoubleValue(key, N_xx_avg_name); + dataN.avg_Q_yy = m_cdbttree->GetDoubleValue(key, N_yy_avg_name); + dataN.avg_Q_xy = m_cdbttree->GetDoubleValue(key, N_xy_avg_name); + + dataN.X_matrix = calculate_flattening_matrix(dataN.avg_Q_xx, dataN.avg_Q_yy, dataN.avg_Q_xy, n, cent_bin, "North"); + + // North South + // Note: We do NOT load avg_Q (x,y) for NS because NS is recentered by summing the recentered S and N vectors. + auto& dataNS = m_correction_data[h_idx][cent_bin][ns_idx]; + + dataNS.avg_Q_xx = m_cdbttree->GetDoubleValue(key, NS_xx_avg_name); + dataNS.avg_Q_yy = m_cdbttree->GetDoubleValue(key, NS_yy_avg_name); + dataNS.avg_Q_xy = m_cdbttree->GetDoubleValue(key, NS_xy_avg_name); + + dataNS.X_matrix = calculate_flattening_matrix(dataNS.avg_Q_xx, dataNS.avg_Q_yy, dataNS.avg_Q_xy, n, cent_bin, "NorthSouth"); + } + } + delete m_cdbttree; + m_cdbttree = nullptr; +} + +//____________________________________________________________________________.. +void EventPlaneRecov2::print_correction_data() +{ + std::cout << std::format("\n{:=>60}\n", ""); + std::cout << std::format("{:^60}\n", "EVENT PLANE CORRECTION DATA SUMMARY"); + std::cout << std::format("{:=>60}\n", ""); + + // Iterate through harmonics {2, 3, 4} + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + int n = m_harmonics[h_idx]; + std::cout << std::format("\n>>> HARMONIC n = {} <<<\n", n); + + // Iterate through Centrality Bins (0-7) + for (size_t cent = 0; cent < m_bins_cent; ++cent) + { + std::cout << std::format("\n Centrality Bin: {}\n", cent); + std::cout << std::format(" {:->30}\n", ""); + + // Header with fixed column widths + std::cout << std::format(" {:<12} {:>10} {:>10} {:>10} {:>10} {:>10}\n", + "Detector", "Avg Qx", "Avg Qy", "Avg Qxx", "Avg Qyy", "Avg Qxy"); + + // Iterate through Subdetectors {S, N} + for (size_t det_idx = 0; det_idx < 3; ++det_idx) + { + std::string det_name; + if (det_idx == 0) + { + det_name = "South"; + } + else if (det_idx == 1) + { + det_name = "North"; + } + else + { + det_name = "NorthSouth"; + } + + const auto& data = m_correction_data[h_idx][cent][det_idx]; + + // For NS, Avg Qx/Qy will be 0.0 because they are not loaded from CDB. + // This is expected behavior. + std::cout << std::format(" {:<12} {:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f}\n", + det_name, + data.avg_Q.x, data.avg_Q.y, + data.avg_Q_xx, data.avg_Q_yy, data.avg_Q_xy); + + // Print X-Matrix in a bracketed layout + std::cout << std::format(" X-Matrix: [ {:>8.6f}, {:>8.6f} ]\n", + data.X_matrix[0][0], data.X_matrix[0][1]); + std::cout << std::format(" [ {:>8.6f}, {:>8.6f} ]\n", + data.X_matrix[1][0], data.X_matrix[1][1]); + } + } + } + std::cout << std::format("\n{:=>60}\n", ""); +} + +int EventPlaneRecov2::CreateNodes(PHCompositeNode *topNode) { + PHNodeIterator iter(topNode); + + PHCompositeNode *dstNode = dynamic_cast(iter.findFirst("PHCompositeNode", "DST")); + if (!dstNode) + { + std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + PHCompositeNode *globalNode = dynamic_cast(iter.findFirst("PHCompositeNode", "GLOBAL")); + if (!globalNode) + { + globalNode = new PHCompositeNode("GLOBAL"); + dstNode->addNode(globalNode); + } + + EventplaneinfoMap *eps = findNode::getClass(topNode, "EventplaneinfoMap"); + if (!eps) + { + eps = new EventplaneinfoMapv1(); + PHIODataNode *newNode = new PHIODataNode(eps , "EventplaneinfoMap", "PHObject"); + globalNode->addNode(newNode); + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +//____________________________________________________________________________.. +int EventPlaneRecov2::process_centrality(PHCompositeNode *topNode) +{ + CentralityInfo* centInfo = findNode::getClass(topNode, "CentralityInfo"); + if (!centInfo) + { + std::cout << PHWHERE << " CentralityInfo is missing doing nothing" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + m_cent = centInfo->get_centile(CentralityInfo::PROP::mbd_NS) * 100; + + if (!std::isfinite(m_cent) || m_cent < 0) + { + if (Verbosity() > 1) + { + std::cout << PHWHERE << " Warning Centrality is out of range. Cent: " << m_cent << ". Cannot calibrate Q vector for this event." << std::endl; + } + m_doNotCalibEvent = true; + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +//____________________________________________________________________________.. +int EventPlaneRecov2::process_sEPD(PHCompositeNode* topNode) +{ + TowerInfoContainer* towerinfosEPD = findNode::getClass(topNode, m_inputNode); + if (!towerinfosEPD) + { + std::cout << PHWHERE << " TOWERINFO_CALIB_SEPD is missing doing nothing" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + EpdGeom* epdgeom = findNode::getClass(topNode, "TOWERGEOM_EPD"); + if (!epdgeom) + { + std::cout << PHWHERE << " TOWERGEOM_EPD is missing doing nothing" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + // sepd + unsigned int nchannels_epd = towerinfosEPD->size(); + + double sepd_total_charge_south = 0; + double sepd_total_charge_north = 0; + + for (unsigned int channel = 0; channel < nchannels_epd; ++channel) + { + TowerInfo* tower = towerinfosEPD->get_tower_at_channel(channel); + + unsigned int key = TowerInfoDefs::encode_epd(channel); + double charge = tower->get_energy(); + double phi = epdgeom->get_phi(key); + + // skip bad channels + // skip channels with very low charge + if (!tower->get_isGood() || charge < m_sepd_min_channel_charge) + { + continue; + } + + // arm = 0: South + // arm = 1: North + unsigned int arm = TowerInfoDefs::get_epd_arm(key); + + // sepd charge sums + double& sepd_total_charge = (arm == 0) ? sepd_total_charge_south : sepd_total_charge_north; + + // Compute total charge for the respective sEPD arm + sepd_total_charge += charge; + + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + int n = m_harmonics[h_idx]; + QVec q_n = {charge * std::cos(n * phi), charge * std::sin(n * phi)}; + m_Q_raw[h_idx][arm].x += q_n.x; + m_Q_raw[h_idx][arm].y += q_n.y; + } + } + + // ensure both total charges are nonzero + if (sepd_total_charge_south == 0 || sepd_total_charge_north == 0) + { + if (Verbosity() > 1) + { + std::cout << PHWHERE << " Error: Total sEPD Charge is Zero: " + << "South = " << sepd_total_charge_south + << ", North = " << sepd_total_charge_north << std::endl; + } + + // ensure raw Q vec is reset + m_Q_raw = {}; + m_doNotCalibEvent = true; + return Fun4AllReturnCodes::EVENT_OK; + } + + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + m_Q_raw[h_idx][0].x /= sepd_total_charge_south; + m_Q_raw[h_idx][0].y /= sepd_total_charge_south; + + m_Q_raw[h_idx][1].x /= sepd_total_charge_north; + m_Q_raw[h_idx][1].y /= sepd_total_charge_north; + + // NEW: Calculate Raw NS (Sum of Raw S + Raw N) + m_Q_raw[h_idx][2].x = m_Q_raw[h_idx][0].x + m_Q_raw[h_idx][1].x; + m_Q_raw[h_idx][2].y = m_Q_raw[h_idx][0].y + m_Q_raw[h_idx][1].y; + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +void EventPlaneRecov2::correct_QVecs() +{ + size_t cent_bin = static_cast(m_cent / 10.0); + if (cent_bin >= m_bins_cent) + { + cent_bin = m_bins_cent - 1; // Clamp max + } + + size_t south_idx = static_cast(Subdetector::S); + size_t north_idx = static_cast(Subdetector::N); + size_t ns_idx = static_cast(Subdetector::NS); + + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + auto& dataS = m_correction_data[h_idx][cent_bin][south_idx]; + auto& dataN = m_correction_data[h_idx][cent_bin][north_idx]; + auto& dataNS = m_correction_data[h_idx][cent_bin][ns_idx]; + + double Q_S_x_avg = dataS.avg_Q.x; + double Q_S_y_avg = dataS.avg_Q.y; + double Q_N_x_avg = dataN.avg_Q.x; + double Q_N_y_avg = dataN.avg_Q.y; + + QVec q_S = m_Q_raw[h_idx][south_idx]; + QVec q_N = m_Q_raw[h_idx][north_idx]; + + // Apply Recentering + QVec q_S_recenter = {q_S.x - Q_S_x_avg, q_S.y - Q_S_y_avg}; + QVec q_N_recenter = {q_N.x - Q_N_x_avg, q_N.y - Q_N_y_avg}; + QVec q_NS_recenter = {q_S_recenter.x + q_N_recenter.x, q_S_recenter.y + q_N_recenter.y}; + + m_Q_recentered[h_idx][south_idx] = q_S_recenter; + m_Q_recentered[h_idx][north_idx] = q_N_recenter; + m_Q_recentered[h_idx][ns_idx] = q_NS_recenter; + + // Flattening Matrix + const auto &X_S = dataS.X_matrix; + const auto &X_N = dataN.X_matrix; + const auto &X_NS = dataNS.X_matrix; + + // Apply Flattening + double Q_S_x_flat = X_S[0][0] * q_S_recenter.x + X_S[0][1] * q_S_recenter.y; + double Q_S_y_flat = X_S[1][0] * q_S_recenter.x + X_S[1][1] * q_S_recenter.y; + double Q_N_x_flat = X_N[0][0] * q_N_recenter.x + X_N[0][1] * q_N_recenter.y; + double Q_N_y_flat = X_N[1][0] * q_N_recenter.x + X_N[1][1] * q_N_recenter.y; + + double Q_NS_x_flat = X_NS[0][0] * q_NS_recenter.x + X_NS[0][1] * q_NS_recenter.y; + double Q_NS_y_flat = X_NS[1][0] * q_NS_recenter.x + X_NS[1][1] * q_NS_recenter.y; + + QVec q_S_flat = {Q_S_x_flat, Q_S_y_flat}; + QVec q_N_flat = {Q_N_x_flat, Q_N_y_flat}; + QVec q_NS_flat = {Q_NS_x_flat, Q_NS_y_flat}; + + m_Q_flat[h_idx][south_idx] = q_S_flat; + m_Q_flat[h_idx][north_idx] = q_N_flat; + m_Q_flat[h_idx][ns_idx] = q_NS_flat; + } +} + +void EventPlaneRecov2::print_QVectors() +{ + std::string header_text = std::format("EVENT Q-VECTOR SUMMARY (Event: {}, CENTRALITY: {:.0f}%)", m_globalEvent, m_cent); + + std::cout << std::format("\n{:*>100}\n", ""); + std::cout << std::format("{:^100}\n", header_text); + std::cout << std::format("{:*>100}\n", ""); + + // Table Header + std::cout << std::format(" {:<10} {:<10} | {:>21} | {:>21} | {:>21}\n", + "Harmonic", "Detector", "Raw (x, y)", "Recentered (x, y)", "Flattened (x, y)"); + std::cout << std::format(" {:-<100}\n", ""); + + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + int n = m_harmonics[h_idx]; + + for (size_t det_idx = 0; det_idx < 3; ++det_idx) + { + std::string det_name; + if (det_idx == 0) + { + det_name = "South"; + } + else if (det_idx == 1) + { + det_name = "North"; + } + else + { + det_name = "NorthSouth"; + } + + const auto& raw = m_Q_raw[h_idx][det_idx]; + const auto& rec = m_Q_recentered[h_idx][det_idx]; + const auto& flat = m_Q_flat[h_idx][det_idx]; + + std::string h_label = (det_idx == 0) ? std::format("n={}", n) : ""; + + // Groups x and y into (val, val) pairs for better scannability + std::string raw_str = std::format("({:>8.5f}, {:>8.5f})", raw.x, raw.y); + std::string rec_str = std::format("({:>8.5f}, {:>8.5f})", rec.x, rec.y); + std::string flat_str = std::format("({:>8.5f}, {:>8.5f})", flat.x, flat.y); + + std::cout << std::format(" {:<10} {:<10} | {:<21} | {:<21} | {:10}\n", + h_label, det_name, raw_str, rec_str, flat_str); + } + if (h_idx < m_harmonics.size() - 1) + { + std::cout << std::format(" {:.>100}\n", ""); + } + } + std::cout << std::format("{:*>100}\n\n", ""); +} + +int EventPlaneRecov2::FillNode(PHCompositeNode *topNode) +{ + EventplaneinfoMap *epmap = findNode::getClass(topNode, "EventplaneinfoMap"); + if (!epmap) + { + std::cout << PHWHERE << " EventplaneinfoMap is missing doing nothing" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + size_t vec_size = static_cast(*std::ranges::max_element(m_harmonics)); + + std::vector> south_Qvec_raw(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + std::vector> south_Qvec_recentered(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + std::vector> south_Qvec(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + + std::vector> north_Qvec_raw(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + std::vector> north_Qvec_recentered(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + std::vector> north_Qvec(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + + std::vector> northsouth_Qvec_raw(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + std::vector> northsouth_Qvec_recentered(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + std::vector> northsouth_Qvec(vec_size, {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}); + + for (size_t h_idx = 0; h_idx < m_harmonics.size(); ++h_idx) + { + int n = m_harmonics[h_idx]; + int idx = n - 1; + + // Fallback logic: Use raw if calibration failed or centrality is out of range + const auto& Q_S = (m_doNotCalib || m_doNotCalibEvent) ? m_Q_raw[h_idx][0] : m_Q_flat[h_idx][0]; + const auto& Q_N = (m_doNotCalib || m_doNotCalibEvent) ? m_Q_raw[h_idx][1] : m_Q_flat[h_idx][1]; + const auto& Q_NS = (m_doNotCalib || m_doNotCalibEvent) ? m_Q_raw[h_idx][2] : m_Q_flat[h_idx][2]; + + const auto& Q_S_raw = m_Q_raw[h_idx][0]; + const auto& Q_S_recentered = m_Q_recentered[h_idx][0]; + + const auto& Q_N_raw = m_Q_raw[h_idx][1]; + const auto& Q_N_recentered = m_Q_recentered[h_idx][1]; + + const auto& Q_NS_raw = m_Q_raw[h_idx][2]; + const auto& Q_NS_recentered = m_Q_recentered[h_idx][2]; + + // South + south_Qvec_raw[idx] = {Q_S_raw.x, Q_S_raw.y}; + south_Qvec_recentered[idx] = {Q_S_recentered.x, Q_S_recentered.y}; + south_Qvec[idx] = {Q_S.x, Q_S.y}; + + // North + north_Qvec_raw[idx] = {Q_N_raw.x, Q_N_raw.y}; + north_Qvec_recentered[idx] = {Q_N_recentered.x, Q_N_recentered.y}; + north_Qvec[idx] = {Q_N.x, Q_N.y}; + + // Combined (North + South) + northsouth_Qvec_raw[idx] = {Q_NS_raw.x, Q_NS_raw.y}; + northsouth_Qvec_recentered[idx] = {Q_NS_recentered.x, Q_NS_recentered.y}; + northsouth_Qvec[idx] = {Q_NS.x, Q_NS.y}; + } + + // Helper lambda to fill nodes using the class's GetPsi method + auto create_and_fill = [&](const std::vector>& qvecs_raw, const std::vector>& qvecs_recentered, const std::vector>& qvecs) { + auto node = std::make_unique(); + node->set_qvector_raw(qvecs_raw); + node->set_qvector_recentered(qvecs_recentered); + node->set_qvector(qvecs); + + std::vector psi_vec(vec_size, std::numeric_limits::quiet_NaN()); + for (int n : m_harmonics) { + psi_vec[n-1] = node->GetPsi(qvecs[n-1].first, qvecs[n-1].second, n); + } + node->set_shifted_psi(psi_vec); + return node; + }; + + epmap->insert(create_and_fill(south_Qvec_raw, south_Qvec_recentered, south_Qvec).release(), EventplaneinfoMap::sEPDS); + epmap->insert(create_and_fill(north_Qvec_raw, north_Qvec_recentered, north_Qvec).release(), EventplaneinfoMap::sEPDN); + epmap->insert(create_and_fill(northsouth_Qvec_raw, northsouth_Qvec_recentered, northsouth_Qvec).release(), EventplaneinfoMap::sEPDNS); + + return Fun4AllReturnCodes::EVENT_OK; +} + +//____________________________________________________________________________.. +int EventPlaneRecov2::process_event(PHCompositeNode *topNode) +{ + EventHeader *eventInfo = findNode::getClass(topNode, "EventHeader"); + if (!eventInfo) + { + return Fun4AllReturnCodes::ABORTRUN; + } + + m_globalEvent = eventInfo->get_EvtSequence(); + + int ret = process_centrality(topNode); + if (ret) + { + return ret; + } + + ret = process_sEPD(topNode); + if (ret) + { + return ret; + } + + // Calibrate Q Vectors + if (!m_doNotCalib && !m_doNotCalibEvent) + { + correct_QVecs(); + } + + ret = FillNode(topNode); + if (ret) + { + return ret; + } + + if (Verbosity() > 1) + { + print_QVectors(); + } + + return Fun4AllReturnCodes::EVENT_OK; +} + +//____________________________________________________________________________.. +int EventPlaneRecov2::ResetEvent(PHCompositeNode */*topNode*/) +{ + m_doNotCalibEvent = false; + + m_Q_raw = {}; + m_Q_recentered = {}; + m_Q_flat = {}; + + return Fun4AllReturnCodes::EVENT_OK; +} diff --git a/offline/packages/eventplaneinfo/EventPlaneRecov2.h b/offline/packages/eventplaneinfo/EventPlaneRecov2.h new file mode 100644 index 0000000000..217faddf88 --- /dev/null +++ b/offline/packages/eventplaneinfo/EventPlaneRecov2.h @@ -0,0 +1,134 @@ +#ifndef EVENTPLANEINFO_EVENTPLANERECOV2_H +#define EVENTPLANEINFO_EVENTPLANERECOV2_H + +#include + + +#include +#include +#include + +class CDBTTree; +class PHCompositeNode; + +class EventPlaneRecov2 : public SubsysReco +{ + public: + + explicit EventPlaneRecov2(const std::string &name = "EventPlaneRecov2"); + ~EventPlaneRecov2() override = default; + + // Explicitly disable copying and moving + EventPlaneRecov2(const EventPlaneRecov2&) = delete; + EventPlaneRecov2& operator=(const EventPlaneRecov2&) = delete; + EventPlaneRecov2(EventPlaneRecov2&&) = delete; + EventPlaneRecov2& operator=(EventPlaneRecov2&&) = delete; + + /** Called during initialization. + Typically this is where you can book histograms, and e.g. + register them to Fun4AllServer (so they can be output to file + using Fun4AllServer::dumpHistos() method). + */ + int Init(PHCompositeNode *topNode) override; + + /** Called for each event. + This is where you do the real work. + */ + int process_event(PHCompositeNode *topNode) override; + + /// Clean up internals after each event. + int ResetEvent(PHCompositeNode *topNode) override; + + + void set_inputNode(const std::string &inputNode) + { + m_inputNode = inputNode; + } + + void set_directURL_EventPlaneCalib(const std::string &directURL_EventPlaneCalib) + { + m_directURL_EventPlaneCalib = directURL_EventPlaneCalib; + } + + void set_doAbortNoEventPlaneCalib(bool status = true) + { + m_doAbortNoEventPlaneCalib = status; + } + + void set_sepd_min_channel_charge(double sepd_min_channel_charge) + { + m_sepd_min_channel_charge = sepd_min_channel_charge; + } + + private: + + static int CreateNodes(PHCompositeNode *topNode); + + std::array, 2> calculate_flattening_matrix(double xx, double yy, double xy, int n, int cent_bin, const std::string& det_label); + void LoadCalib(); + + void print_correction_data(); + void print_QVectors(); + + int process_centrality(PHCompositeNode *topNode); + int process_sEPD(PHCompositeNode *topNode); + void correct_QVecs(); + + int FillNode(PHCompositeNode *topNode); + + std::string m_directURL_EventPlaneCalib; + bool m_doAbortNoEventPlaneCalib{false}; + bool m_doNotCalib{false}; + bool m_doNotCalibEvent{false}; + + double m_cent{0.0}; + double m_globalEvent{0}; + double m_sepd_min_channel_charge{0.2}; + + std::string m_calibName{"SEPD_EventPlaneCalib"}; + std::string m_inputNode{"TOWERINFO_CALIB_SEPD"}; + + CDBTTree *m_cdbttree {nullptr}; + + enum class Subdetector + { + S, + N, + NS + }; + + struct QVec + { + double x{0.0}; + double y{0.0}; + }; + + struct CorrectionData + { + // Averages of Qx, Qy, Qx^2, Qy^2, Qxy + QVec avg_Q{}; + double avg_Q_xx{0.0}; + double avg_Q_yy{0.0}; + double avg_Q_xy{0.0}; + + // Correction matrix + std::array, 2> X_matrix{}; + }; + + static constexpr size_t m_bins_cent {8}; + static constexpr std::array m_harmonics = {2, 3, 4}; + + // Holds all correction data + // key: [Harmonic][Cent][Subdetector] + // Harmonics {2,3,4} -> 3 elements + // Subdetectors {S,N,NS} -> 3 elements + std::array, m_bins_cent>, m_harmonics.size()> m_correction_data; + + // sEPD Q Vectors + // key: [Harmonic][Subdetector] + // Subdetectors {S,N,NS} -> 3 elements + std::array, m_harmonics.size()> m_Q_raw{}; + std::array, m_harmonics.size()> m_Q_recentered{}; + std::array, m_harmonics.size()> m_Q_flat{}; +}; +#endif diff --git a/offline/packages/eventplaneinfo/Eventplaneinfo.h b/offline/packages/eventplaneinfo/Eventplaneinfo.h index 1ad2665881..af903b6d06 100644 --- a/offline/packages/eventplaneinfo/Eventplaneinfo.h +++ b/offline/packages/eventplaneinfo/Eventplaneinfo.h @@ -1,12 +1,12 @@ // Tell emacs that this is a C++ source // -*- C++ -*-. -#ifndef EVENTPLANEINFO_H -#define EVENTPLANEINFO_H +#ifndef EVENTPLANEINFO_EVENTPLANEINFO_H +#define EVENTPLANEINFO_EVENTPLANEINFO_H #include -#include #include +#include #include #include @@ -23,17 +23,21 @@ class Eventplaneinfo : public PHObject PHObject* CloneMe() const override { return nullptr; } virtual void set_qvector(std::vector> /*Qvec*/) { return; } + virtual void set_qvector_raw(const std::vector>& /*Qvec*/) { return; } + virtual void set_qvector_recentered(const std::vector>& /*Qvec*/) { return; } virtual void set_shifted_psi(std::vector /*Psi_Shifted*/) { return; } - virtual std::pair get_qvector(int /*order*/) const { return std::make_pair(NAN, NAN); } - virtual double get_psi(int /*order*/) const { return NAN; } - virtual double get_shifted_psi(int /*order*/) const { return NAN; } - virtual double GetPsi(const double /*Qx*/, const double /*Qy*/, const unsigned int /*order*/) const { return NAN; } + virtual std::pair get_qvector(int /*order*/) const { return std::make_pair(std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); } + virtual std::pair get_qvector_raw(int /*order*/) const { return std::make_pair(std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); } + virtual std::pair get_qvector_recentered(int /*order*/) const { return std::make_pair(std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); } + virtual double get_psi(int /*order*/) const { return std::numeric_limits::quiet_NaN(); } + virtual double get_shifted_psi(int /*order*/) const { return std::numeric_limits::quiet_NaN(); } + virtual double GetPsi(const double /*Qx*/, const double /*Qy*/, const unsigned int /*order*/) const { return std::numeric_limits::quiet_NaN(); } virtual void set_ring_qvector(std::vector>> /*RingQvecs*/) { return; } - virtual std::pair get_ring_qvector(int /*rbin*/, int /*order*/) const { return std::make_pair(NAN, NAN); } - virtual double get_ring_psi(int /*rbin*/, int /*order*/) const { return NAN; } + virtual std::pair get_ring_qvector(int /*rbin*/, int /*order*/) const { return std::make_pair(std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); } + virtual double get_ring_psi(int /*rbin*/, int /*order*/) const { return std::numeric_limits::quiet_NaN(); } protected: - Eventplaneinfo() {} + Eventplaneinfo() = default; private: ClassDefOverride(Eventplaneinfo, 1); diff --git a/offline/packages/eventplaneinfo/Eventplaneinfov2.cc b/offline/packages/eventplaneinfo/Eventplaneinfov2.cc new file mode 100644 index 0000000000..7331117a17 --- /dev/null +++ b/offline/packages/eventplaneinfo/Eventplaneinfov2.cc @@ -0,0 +1,23 @@ +#include "Eventplaneinfov2.h" + +#include +#include + +void Eventplaneinfov2::identify(std::ostream& os) const +{ + os << "---------Eventplaneinfov2------------------" << std::endl; + return; +} + +double Eventplaneinfov2::GetPsi(const double Qx, const double Qy, const unsigned int order) const +{ + if (order == 0) + { + return std::numeric_limits::quiet_NaN(); + } + if ((Qx == 0.0) && (Qy == 0.0)) + { + return std::numeric_limits::quiet_NaN(); + } + return std::atan2(Qy, Qx) / order; +} diff --git a/offline/packages/eventplaneinfo/Eventplaneinfov2.h b/offline/packages/eventplaneinfo/Eventplaneinfov2.h new file mode 100644 index 0000000000..bd6cb553c8 --- /dev/null +++ b/offline/packages/eventplaneinfo/Eventplaneinfov2.h @@ -0,0 +1,86 @@ +// Tell emacs that this is a C++ source +// -*- C++ -*-. +#ifndef EVENTPLANEINFO_EVENTPLANEINFOV2_H +#define EVENTPLANEINFO_EVENTPLANEINFOV2_H + +#include "Eventplaneinfo.h" + +#include // for size_t +#include +#include +#include // for pair, make_pair +#include + +class PHObject; + +class Eventplaneinfov2 : public Eventplaneinfo +{ + public: + Eventplaneinfov2() = default; + ~Eventplaneinfov2() override = default; + + Eventplaneinfov2(const Eventplaneinfov2&) = default; + Eventplaneinfov2& operator=(const Eventplaneinfov2&) = default; + Eventplaneinfov2(Eventplaneinfov2&&) = default; + Eventplaneinfov2& operator=(Eventplaneinfov2&&) = default; + + void identify(std::ostream& os = std::cout) const override; + void Reset() override { *this = Eventplaneinfov2(); } + PHObject* CloneMe() const override { return new Eventplaneinfov2(*this); } + + void set_qvector(std::vector> Qvec) override { mQvec = Qvec; } + void set_qvector_raw(const std::vector>& Qvec) override { mQvec_raw = Qvec; } + void set_qvector_recentered(const std::vector>& Qvec) override { mQvec_recentered = Qvec; } + void set_shifted_psi(std::vector Psi_Shifted) override { mPsi_Shifted = Psi_Shifted; } + std::pair get_qvector(int order) const override { return safe_qvec(mQvec, order); } + std::pair get_qvector_raw(int order) const override { return safe_qvec(mQvec_raw, order); } + std::pair get_qvector_recentered(int order) const override { return safe_qvec(mQvec_recentered, order); } + void set_ring_qvector(std::vector>> Qvec) override { ring_Qvec = Qvec; } + std::pair get_ring_qvector(int ring_index, int order) const override + { + if (ring_index < 0 || static_cast(ring_index) >= ring_Qvec.size()) + { + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; + } + return safe_qvec(ring_Qvec[ring_index], order); + } + double get_ring_psi(int ring_index, int order) const override + { + auto q = get_ring_qvector(ring_index, order); + return GetPsi(q.first, q.second, static_cast(order)); + } + + double GetPsi(double Qx, double Qy, unsigned int order) const override; + double get_psi(int order) const override + { + auto q = get_qvector(order); + return GetPsi(q.first, q.second, static_cast(order)); + } + double get_shifted_psi(int order) const override + { + if (order <= 0 || static_cast(order) > mPsi_Shifted.size()) + { + return std::numeric_limits::quiet_NaN(); + } + return mPsi_Shifted[order - 1]; + } + + private: + static std::pair safe_qvec(const std::vector>& v, int order) + { + if (order <= 0 || static_cast(order) > v.size()) + { + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; + } + return v[order - 1]; + } + + std::vector> mQvec; + std::vector> mQvec_raw; + std::vector> mQvec_recentered; + std::vector mPsi_Shifted; + std::vector>> ring_Qvec; + ClassDefOverride(Eventplaneinfov2, 1); +}; + +#endif diff --git a/offline/packages/eventplaneinfo/Eventplaneinfov2LinkDef.h b/offline/packages/eventplaneinfo/Eventplaneinfov2LinkDef.h new file mode 100644 index 0000000000..961c0446cf --- /dev/null +++ b/offline/packages/eventplaneinfo/Eventplaneinfov2LinkDef.h @@ -0,0 +1,5 @@ +#ifdef __CINT__ + +#pragma link C++ class Eventplaneinfov2 + ; + +#endif /* __CINT__ */ diff --git a/offline/packages/eventplaneinfo/Makefile.am b/offline/packages/eventplaneinfo/Makefile.am index 3a1192095d..786a285bee 100644 --- a/offline/packages/eventplaneinfo/Makefile.am +++ b/offline/packages/eventplaneinfo/Makefile.am @@ -14,8 +14,7 @@ AM_LDFLAGS = \ -L$(OFFLINE_MAIN)/lib libeventplaneinfo_io_la_LIBADD = \ - -lphool \ - -lcentrality_io + -lphool libeventplaneinfo_la_LIBADD = \ libeventplaneinfo_io.la \ @@ -24,6 +23,7 @@ libeventplaneinfo_la_LIBADD = \ -lfun4all \ -lffamodules \ -lcalotrigger_io \ + -lcentrality_io \ -lffarawobjects \ -lcdbobjects \ -lglobalvertex_io @@ -32,13 +32,16 @@ pkginclude_HEADERS = \ EventPlaneCalibration.h \ Eventplaneinfo.h \ Eventplaneinfov1.h \ + Eventplaneinfov2.h \ EventplaneinfoMap.h \ EventplaneinfoMapv1.h \ - EventPlaneReco.h + EventPlaneReco.h \ + EventPlaneRecov2.h ROOTDICTS = \ Eventplaneinfo_Dict.cc \ Eventplaneinfov1_Dict.cc \ + Eventplaneinfov2_Dict.cc \ EventplaneinfoMap_Dict.cc \ EventplaneinfoMapv1_Dict.cc @@ -50,12 +53,14 @@ libeventplaneinfo_io_la_SOURCES = \ $(ROOTDICTS) \ Eventplaneinfo.cc \ Eventplaneinfov1.cc \ + Eventplaneinfov2.cc \ EventplaneinfoMap.cc \ EventplaneinfoMapv1.cc libeventplaneinfo_la_SOURCES = \ EventPlaneCalibration.cc \ - EventPlaneReco.cc + EventPlaneReco.cc \ + EventPlaneRecov2.cc # Rule for generating table CINT dictionaries. %_Dict.cc: %.h %LinkDef.h