From d2a55048b197db66da065921304c6b21584ec8f6 Mon Sep 17 00:00:00 2001 From: mcyoren Date: Tue, 8 Sep 2026 16:25:39 -0400 Subject: [PATCH 1/6] Move gem currents to InitRun and Keff adjustment by segment --- .../packages/TpcConditions/TpcConditions.h | 33 ++++ .../TpcConditions/TpcConditionsReco.cc | 151 ++++++++++++++++-- .../TpcConditions/TpcConditionsReco.h | 2 + offline/packages/tpctrackreco/Makefile.am | 3 +- .../tpctrackreco/Tpc_PolyClusterizer.cc | 20 +++ .../tpctrackreco/Tpc_PolyClusterizer.h | 2 + 6 files changed, 196 insertions(+), 15 deletions(-) diff --git a/offline/packages/TpcConditions/TpcConditions.h b/offline/packages/TpcConditions/TpcConditions.h index e0ef4386d7..1e1cd67bad 100644 --- a/offline/packages/TpcConditions/TpcConditions.h +++ b/offline/packages/TpcConditions/TpcConditions.h @@ -32,6 +32,17 @@ class TpcConditions void set_LoadNR1(float Value) { m_LoadNR1 = Value; } void set_LoadNR2(float Value) { m_LoadNR2 = Value; } void set_LoadNR3(float Value) { m_LoadNR3 = Value; } + + void set_AverageLoadCurrent(float Value) { m_AverageLoadCurrent = Value; } + void set_AverageLoadNorth(float Value) { m_AverageLoadNorth = Value; } + void set_AverageLoadSouth(float Value) { m_AverageLoadSouth = Value; } + void set_AverageLoadSR1(float Value) { m_AverageLoadSR1 = Value; } + void set_AverageLoadSR2(float Value) { m_AverageLoadSR2 = Value; } + void set_AverageLoadSR3(float Value) { m_AverageLoadSR3 = Value; } + void set_AverageLoadNR1(float Value) { m_AverageLoadNR1 = Value; } + void set_AverageLoadNR2(float Value) { m_AverageLoadNR2 = Value; } + void set_AverageLoadNR3(float Value) { m_AverageLoadNR3 = Value; } + void set_Temperature(float Value) { m_Temperature = Value; } void set_Pressure(float Value) { m_Pressure = Value; } void set_FieldOK(bool Value) { m_FieldOK = Value; } @@ -46,6 +57,17 @@ class TpcConditions float get_LoadNR1() const { return m_LoadNR1; } float get_LoadNR2() const { return m_LoadNR2; } float get_LoadNR3() const { return m_LoadNR3; } + + float get_AverageLoadCurrent() const { return m_AverageLoadCurrent; } + float get_AverageLoadNorth() const { return m_AverageLoadNorth; } + float get_AverageLoadSouth() const { return m_AverageLoadSouth; } + float get_AverageLoadSR1() const { return m_AverageLoadSR1; } + float get_AverageLoadSR2() const { return m_AverageLoadSR2; } + float get_AverageLoadSR3() const { return m_AverageLoadSR3; } + float get_AverageLoadNR1() const { return m_AverageLoadNR1; } + float get_AverageLoadNR2() const { return m_AverageLoadNR2; } + float get_AverageLoadNR3() const { return m_AverageLoadNR3; } + float get_Temperature() const { return m_Temperature; } float get_Pressure() const { return m_Pressure; } bool get_FieldOK() const { return m_FieldOK; } @@ -61,6 +83,17 @@ class TpcConditions float m_LoadNR1{0.0}; float m_LoadNR2{0.0}; float m_LoadNR3{0.0}; + + float m_AverageLoadCurrent{0.0}; + float m_AverageLoadNorth{0.0}; + float m_AverageLoadSouth{0.0}; + float m_AverageLoadSR1{0.0}; + float m_AverageLoadSR2{0.0}; + float m_AverageLoadSR3{0.0}; + float m_AverageLoadNR1{0.0}; + float m_AverageLoadNR2{0.0}; + float m_AverageLoadNR3{0.0}; + float m_Temperature{0.0}; float m_Pressure{0.0}; bool m_FieldOK{false}; diff --git a/offline/packages/TpcConditions/TpcConditionsReco.cc b/offline/packages/TpcConditions/TpcConditionsReco.cc index b5430a54e0..982aeab4e7 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.cc +++ b/offline/packages/TpcConditions/TpcConditionsReco.cc @@ -14,9 +14,12 @@ #include #include +#include + #include #include #include +#include #include #include #include @@ -105,11 +108,76 @@ TpcConditionsReco::~TpcConditionsReco() delete m_tree; } +void TpcConditionsReco::fillConditions(int channel) +{ + m_conditions->set_Temperature(m_tree->GetFloatValue(channel, "gas_temperature")); + m_conditions->set_Pressure(m_tree->GetFloatValue(channel, "gas_pressure")); + m_conditions->set_FieldOK(m_tree->GetFloatValue(channel, "FieldOK") != 0.0); + m_conditions->set_GainOK(m_tree->GetFloatValue(channel, "GainOK") != 0.0); + + m_conditions->set_LoadCurrent(get_MedianCurrent(channel, ALL)); + m_conditions->set_LoadNorth(get_MedianCurrent(channel, NORTH)); + m_conditions->set_LoadSouth(get_MedianCurrent(channel, SOUTH)); + m_conditions->set_LoadSR1(get_MedianCurrent(channel, SR1)); + m_conditions->set_LoadSR2(get_MedianCurrent(channel, SR2)); + m_conditions->set_LoadSR3(get_MedianCurrent(channel, SR3)); + m_conditions->set_LoadNR1(get_MedianCurrent(channel, NR1)); + m_conditions->set_LoadNR2(get_MedianCurrent(channel, NR2)); + m_conditions->set_LoadNR3(get_MedianCurrent(channel, NR3)); +} + +float TpcConditionsReco::get_AverageMedianCurrent( + const std::vector &channels) +{ + if (m_bco_to_channel.empty()) + { + return 0.0F; + } + + double sum = 0.0; + + for (const auto &[bco, channel] : m_bco_to_channel) + { + sum += get_MedianCurrent(channel, channels); + } + + return static_cast(sum / static_cast(m_bco_to_channel.size())); +} + int TpcConditionsReco::InitRun(PHCompositeNode *topNode) { std::cout << "TpcConditionsReco::InitRun(PHCompositeNode *topNode) Initializing" << std::endl; + Gl1Packet *gl1 = findNode::getClass(topNode, "GL1RAWHIT"); + if (!gl1) + { + gl1 = findNode::getClass(topNode, "GL1Packet"); + } + + if (!gl1) + { + std::cout << "TpcConditionsReco::InitRun - could not find GL1 packet" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + const int64_t rawBco = gl1->lValue(0, "BCO"); + + if (rawBco <= 0) + { + std::cout << "TpcConditionsReco::InitRun - invalid GL1 BCO " + << rawBco << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + const uint64_t targetBco = static_cast(rawBco); + + std::cout << "\033[1;31m" + << " targetBCO = " + << targetBco + << "\033[0m" << std::endl; + // Get the TPC conditions payload from CDB std::string calibdir = CDBInterface::instance()->getUrl("TPC_CONDITIONS"); m_tree = new CDBTTree(calibdir); @@ -140,6 +208,75 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) PHDataNode *conditionsNode = new PHDataNode(m_conditions, "TpcConditions", "Data"); parNode->addNode(conditionsNode); + m_conditions->set_AverageLoadCurrent(get_AverageMedianCurrent(ALL)); + m_conditions->set_AverageLoadNorth(get_AverageMedianCurrent(NORTH)); + m_conditions->set_AverageLoadSouth(get_AverageMedianCurrent(SOUTH)); + + m_conditions->set_AverageLoadSR1(get_AverageMedianCurrent(SR1)); + m_conditions->set_AverageLoadSR2(get_AverageMedianCurrent(SR2)); + m_conditions->set_AverageLoadSR3(get_AverageMedianCurrent(SR3)); + + m_conditions->set_AverageLoadNR1(get_AverageMedianCurrent(NR1)); + m_conditions->set_AverageLoadNR2(get_AverageMedianCurrent(NR2)); + m_conditions->set_AverageLoadNR3(get_AverageMedianCurrent(NR3)); + + if (!m_bco_to_channel.empty()) + { + recoConsts *rc = recoConsts::instance(); + const int segment = rc->get_IntFlag("RUNSEGMENT"); + + // Find closest stored conditions BCO. + auto upper = m_bco_to_channel.upper_bound(targetBco); + auto selected = upper; + + if (upper == m_bco_to_channel.end()) + { + selected = std::prev(m_bco_to_channel.end()); + } + else if (upper != m_bco_to_channel.begin()) + { + auto lower = std::prev(upper); + + if ((targetBco - lower->first) <= (upper->first - targetBco)) + { + selected = lower; + } + } + + if (!m_bco_to_channel.empty()) + { + const auto first = m_bco_to_channel.begin(); + const auto last = std::prev(m_bco_to_channel.end()); + const auto middle = std::next( + m_bco_to_channel.begin(), + static_cast(m_bco_to_channel.size() / 2)); + + std::cout << "\033[1;31m" + << "R1 conditions:" + << "\n FIRST BCO=" << first->first + << " NR1=" << get_MedianCurrent(first->second, NR1) + << " SR1=" << get_MedianCurrent(first->second, SR1) + << "\n MIDDLE BCO=" << middle->first + << " NR1=" << get_MedianCurrent(middle->second, NR1) + << " SR1=" << get_MedianCurrent(middle->second, SR1) + << "\n LAST BCO=" << last->first + << " NR1=" << get_MedianCurrent(last->second, NR1) + << " SR1=" << get_MedianCurrent(last->second, SR1) + << "\n AVERAGE NR1=" << m_conditions->get_AverageLoadNR1() + << " SR1=" << m_conditions->get_AverageLoadSR1() + << "\033[0m" << std::endl; + } + + fillConditions(selected->second); + + std::cout << "TpcConditionsReco::InitRun - segment " << segment + << ", target BCO " << targetBco + << ", conditions BCO " << selected->first + << ", LoadNR1 " << m_conditions->get_LoadNR1() + << ", LoadSR1 " << m_conditions->get_LoadSR1() + << std::endl; + } + return Fun4AllReturnCodes::EVENT_OK; } @@ -157,7 +294,6 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) uint64_t bco = gl1->getBCO(); - // Find the most recent conditions sample at or before this BCO auto iter = m_bco_to_channel.upper_bound(bco); if (iter == m_bco_to_channel.begin()) { @@ -170,22 +306,9 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) --iter; int channel = iter->second; - // Fill current-event conditions - m_conditions->set_Temperature(m_tree->GetFloatValue(channel, "gas_temperature")); - m_conditions->set_Pressure(m_tree->GetFloatValue(channel, "gas_pressure")); m_conditions->set_FieldOK(m_tree->GetFloatValue(channel, "FieldOK") != 0.0); m_conditions->set_GainOK(m_tree->GetFloatValue(channel, "GainOK") != 0.0); - // Here are the critical LoadCurrent calculations. - m_conditions->set_LoadCurrent(get_MedianCurrent(channel, ALL)); - m_conditions->set_LoadNorth(get_MedianCurrent(channel, NORTH)); - m_conditions->set_LoadSouth(get_MedianCurrent(channel, SOUTH)); - m_conditions->set_LoadSR1(get_MedianCurrent(channel, SR1)); - m_conditions->set_LoadSR2(get_MedianCurrent(channel, SR2)); - m_conditions->set_LoadSR3(get_MedianCurrent(channel, SR3)); - m_conditions->set_LoadNR1(get_MedianCurrent(channel, NR1)); - m_conditions->set_LoadNR2(get_MedianCurrent(channel, NR2)); - m_conditions->set_LoadNR3(get_MedianCurrent(channel, NR3)); // Do or die if (!m_conditions->get_FieldOK() || diff --git a/offline/packages/TpcConditions/TpcConditionsReco.h b/offline/packages/TpcConditions/TpcConditionsReco.h index ffbd2887cb..e3b38a4eff 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.h +++ b/offline/packages/TpcConditions/TpcConditionsReco.h @@ -26,6 +26,8 @@ class TpcConditionsReco : public SubsysReco CDBTTree *m_tree{nullptr}; TpcConditions *m_conditions{nullptr}; float get_MedianCurrent(int channel, const std::vector &channels); + void fillConditions(int channel); + float get_AverageMedianCurrent(const std::vector &channels); std::map m_bco_to_channel; }; diff --git a/offline/packages/tpctrackreco/Makefile.am b/offline/packages/tpctrackreco/Makefile.am index ea5b489293..d199358951 100644 --- a/offline/packages/tpctrackreco/Makefile.am +++ b/offline/packages/tpctrackreco/Makefile.am @@ -141,7 +141,8 @@ libtpctrackreco_la_LIBADD = \ -lgenfit2exp \ -lPHGenFit \ -lphfield \ - -lPHGarfield + -lPHGarfield \ + -lTpcConditions # Rule for generating table CINT dictionaries. %_Dict.cc: %.h %LinkDef.h diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc index 22503a4fb6..efe626ca15 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc @@ -24,6 +24,8 @@ #include #include +#include + #include // for PHG4CylinderGeom #include #include @@ -226,6 +228,17 @@ bool Tpc_PolyClusterizer::load_cdb_inputs() return " [CDB]"; }; + std::cout << "\033[31m" + << "Load NR1 " << m_conditions->get_LoadNR1() + << " Load SR1 " << m_conditions->get_LoadSR1() + << "\033[0m" << std::endl; + + if(true) + { + m_kEffSide0 = m_kEffSide0 * m_conditions->get_LoadSR1()/m_conditions->get_AverageLoadSR1(); + m_kEffSide1 = m_kEffSide1 * m_conditions->get_LoadNR1()/m_conditions->get_AverageLoadNR1(); + } + std::cout << Name() << "::load_cdb_inputs - final kEff values: side0 = " << m_kEffSide0 << keff_source(m_kEffSide0Override) << ", side1 = " << m_kEffSide1 << keff_source(m_kEffSide1Override) << std::endl; return ok; @@ -362,6 +375,13 @@ int Tpc_PolyClusterizer::getNodes(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTRUN; } + m_conditions = findNode::getClass(topNode, "TpcConditions"); + if (!m_conditions) + { + std::cerr << Name() << "::getNodes - missing TpcConditions" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + return Fun4AllReturnCodes::EVENT_OK; } diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h index 1c09259dd2..6e4ea564be 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h @@ -20,6 +20,7 @@ class TpcCrossingDecisionContainer; class TrkrHitSetContainer; class PHG4CylinderGeomContainer; class PHG4TpcGeomContainer; +class TpcConditions; class Tpc_PolyClusterizer : public SubsysReco { @@ -180,6 +181,7 @@ class Tpc_PolyClusterizer : public SubsysReco std::unique_ptr m_garfield{}; PHG4TpcGeomContainer* m_geomContainerTpc{nullptr}; + TpcConditions* m_conditions{nullptr}; std::array m_driftLookup; unsigned int m_event{0}; double m_t0{8}; From 6a3a1011b32b47d5f7fa0905a26ebe76b50f5355 Mon Sep 17 00:00:00 2001 From: mcyoren Date: Wed, 9 Sep 2026 16:37:05 -0400 Subject: [PATCH 2/6] Switch IFC Voltage by runnumber; Add time-dependent TPC kEff correction; Survey Geometry is on by default --- offline/framework/phool/RunnumberRange.h | 1 + .../TpcConditions/TpcConditionsReco.cc | 193 +++++++++++++----- .../TpcConditions/TpcConditionsReco.h | 1 + .../tpctrackreco/Tpc_PolyClusterizer.cc | 72 +++++-- .../tpctrackreco/Tpc_PolyClusterizer.h | 2 + 5 files changed, 207 insertions(+), 62 deletions(-) diff --git a/offline/framework/phool/RunnumberRange.h b/offline/framework/phool/RunnumberRange.h index 7aadfd7d72..29b45bd404 100644 --- a/offline/framework/phool/RunnumberRange.h +++ b/offline/framework/phool/RunnumberRange.h @@ -26,6 +26,7 @@ namespace RunnumberRange constexpr int RUN2AUAU_LAST = 54974; constexpr int RUN3_TPCFW_CLOCK_CHANGE = 58667; constexpr int RUN3AUAU_FIRST = 66457; + constexpr int RUN3AUAU_IFC_V_CHANGE = 75800; // first run after IFC voltage change constexpr int RUN3AUAU_LAST = 78954; constexpr int RUN3PP_FIRST = 79146; // first beam data constexpr int RUN3PP_LAST = 81668; diff --git a/offline/packages/TpcConditions/TpcConditionsReco.cc b/offline/packages/TpcConditions/TpcConditionsReco.cc index 982aeab4e7..ae5e4d3eee 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.cc +++ b/offline/packages/TpcConditions/TpcConditionsReco.cc @@ -19,8 +19,8 @@ #include #include #include -#include #include +#include #include #include @@ -108,6 +108,114 @@ TpcConditionsReco::~TpcConditionsReco() delete m_tree; } +float TpcConditionsReco::get_AverageMedianCurrent( + const std::vector &channels) +{ + if (m_bco_to_channel.empty()) + { + return 0.0F; + } + + double sum = 0.0; + + for (const auto &entry : m_bco_to_channel) + { + sum += get_MedianCurrent(entry.second, channels); + } + + return static_cast(sum / static_cast(m_bco_to_channel.size())); +} + +float TpcConditionsReco::get_InterpolatedMedianCurrent(const uint64_t bco, const std::vector &channels) +{ + if (m_bco_to_channel.empty()) + { + return 0.0F; + } + + if (m_bco_to_channel.size() == 1) + { + return get_MedianCurrent(m_bco_to_channel.begin()->second, channels); + } + + const auto first = m_bco_to_channel.begin(); + const auto last = std::prev(m_bco_to_channel.end()); + + if (bco <= first->first) + { + return get_MedianCurrent(first->second, channels); + } + + if (bco >= last->first) + { + return get_MedianCurrent(last->second, channels); + } + + const auto right = m_bco_to_channel.upper_bound(bco); + const auto left = std::prev(right); + + const double x1 = static_cast(left->first); + const double x2 = static_cast(right->first); + const double y1 = get_MedianCurrent(left->second, channels); + const double y2 = get_MedianCurrent(right->second, channels); + + const double h1 = x2 - x1; + const double d1 = (y2 - y1) / h1; + + double m1 = d1; + double m2 = d1; + + if (left != m_bco_to_channel.begin()) + { + const auto left2 = std::prev(left); + const double x0 = static_cast(left2->first); + const double y0 = get_MedianCurrent(left2->second, channels); + const double h0 = x1 - x0; + const double d0 = (y1 - y0) / h0; + + if (d0 * d1 > 0.0) + { + const double w1 = 2.0 * h1 + h0; + const double w2 = h1 + 2.0 * h0; + m1 = (w1 + w2) / (w1 / d0 + w2 / d1); + } + else + { + m1 = 0.0; + } + } + + const auto right2 = std::next(right); + if (right2 != m_bco_to_channel.end()) + { + const double x3 = static_cast(right2->first); + const double y3 = get_MedianCurrent(right2->second, channels); + const double h2 = x3 - x2; + const double d2 = (y3 - y2) / h2; + + if (d1 * d2 > 0.0) + { + const double w1 = 2.0 * h2 + h1; + const double w2 = h2 + 2.0 * h1; + m2 = (w1 + w2) / (w1 / d1 + w2 / d2); + } + else + { + m2 = 0.0; + } + } + + const double t = (static_cast(bco) - x1) / h1; + const double t2 = t * t; + const double t3 = t2 * t; + + return static_cast( + (2.0 * t3 - 3.0 * t2 + 1.0) * y1 + + (t3 - 2.0 * t2 + t) * h1 * m1 + + (-2.0 * t3 + 3.0 * t2) * y2 + + (t3 - t2) * h1 * m2); +} + void TpcConditionsReco::fillConditions(int channel) { m_conditions->set_Temperature(m_tree->GetFloatValue(channel, "gas_temperature")); @@ -126,24 +234,6 @@ void TpcConditionsReco::fillConditions(int channel) m_conditions->set_LoadNR3(get_MedianCurrent(channel, NR3)); } -float TpcConditionsReco::get_AverageMedianCurrent( - const std::vector &channels) -{ - if (m_bco_to_channel.empty()) - { - return 0.0F; - } - - double sum = 0.0; - - for (const auto &[bco, channel] : m_bco_to_channel) - { - sum += get_MedianCurrent(channel, channels); - } - - return static_cast(sum / static_cast(m_bco_to_channel.size())); -} - int TpcConditionsReco::InitRun(PHCompositeNode *topNode) { std::cout << "TpcConditionsReco::InitRun(PHCompositeNode *topNode) Initializing" @@ -243,37 +333,47 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) } } - if (!m_bco_to_channel.empty()) - { - const auto first = m_bco_to_channel.begin(); - const auto last = std::prev(m_bco_to_channel.end()); - const auto middle = std::next( - m_bco_to_channel.begin(), - static_cast(m_bco_to_channel.size() / 2)); - - std::cout << "\033[1;31m" - << "R1 conditions:" - << "\n FIRST BCO=" << first->first - << " NR1=" << get_MedianCurrent(first->second, NR1) - << " SR1=" << get_MedianCurrent(first->second, SR1) - << "\n MIDDLE BCO=" << middle->first - << " NR1=" << get_MedianCurrent(middle->second, NR1) - << " SR1=" << get_MedianCurrent(middle->second, SR1) - << "\n LAST BCO=" << last->first - << " NR1=" << get_MedianCurrent(last->second, NR1) - << " SR1=" << get_MedianCurrent(last->second, SR1) - << "\n AVERAGE NR1=" << m_conditions->get_AverageLoadNR1() - << " SR1=" << m_conditions->get_AverageLoadSR1() - << "\033[0m" << std::endl; - } + const auto first = m_bco_to_channel.begin(); + const auto last = std::prev(m_bco_to_channel.end()); + const auto middle = std::next( + m_bco_to_channel.begin(), + static_cast(m_bco_to_channel.size() / 2)); + + std::cout << "\033[1;31m" + << "R1 conditions:" + << "\n FIRST BCO=" << first->first + << " NR1=" << get_MedianCurrent(first->second, NR1) + << " SR1=" << get_MedianCurrent(first->second, SR1) + << "\n MIDDLE BCO=" << middle->first + << " NR1=" << get_MedianCurrent(middle->second, NR1) + << " SR1=" << get_MedianCurrent(middle->second, SR1) + << "\n LAST BCO=" << last->first + << " NR1=" << get_MedianCurrent(last->second, NR1) + << " SR1=" << get_MedianCurrent(last->second, SR1) + << "\n AVERAGE NR1=" << m_conditions->get_AverageLoadNR1() + << " SR1=" << m_conditions->get_AverageLoadSR1() + << "\033[0m" << std::endl; fillConditions(selected->second); + // GEM load currents are evaluated smoothly at the actual segment BCO. + m_conditions->set_LoadCurrent(get_InterpolatedMedianCurrent(targetBco, ALL)); + m_conditions->set_LoadNorth(get_InterpolatedMedianCurrent(targetBco, NORTH)); + m_conditions->set_LoadSouth(get_InterpolatedMedianCurrent(targetBco, SOUTH)); + + m_conditions->set_LoadSR1(get_InterpolatedMedianCurrent(targetBco, SR1)); + m_conditions->set_LoadSR2(get_InterpolatedMedianCurrent(targetBco, SR2)); + m_conditions->set_LoadSR3(get_InterpolatedMedianCurrent(targetBco, SR3)); + + m_conditions->set_LoadNR1(get_InterpolatedMedianCurrent(targetBco, NR1)); + m_conditions->set_LoadNR2(get_InterpolatedMedianCurrent(targetBco, NR2)); + m_conditions->set_LoadNR3(get_InterpolatedMedianCurrent(targetBco, NR3)); + std::cout << "TpcConditionsReco::InitRun - segment " << segment << ", target BCO " << targetBco - << ", conditions BCO " << selected->first - << ", LoadNR1 " << m_conditions->get_LoadNR1() - << ", LoadSR1 " << m_conditions->get_LoadSR1() + << ", nearest conditions BCO " << selected->first + << ", interpolated LoadNR1 " << m_conditions->get_LoadNR1() + << ", interpolated LoadSR1 " << m_conditions->get_LoadSR1() << std::endl; } @@ -306,10 +406,11 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) --iter; int channel = iter->second; + m_conditions->set_Temperature(m_tree->GetFloatValue(channel, "gas_temperature")); + m_conditions->set_Pressure(m_tree->GetFloatValue(channel, "gas_pressure")); m_conditions->set_FieldOK(m_tree->GetFloatValue(channel, "FieldOK") != 0.0); m_conditions->set_GainOK(m_tree->GetFloatValue(channel, "GainOK") != 0.0); - // Do or die if (!m_conditions->get_FieldOK() || !m_conditions->get_GainOK()) diff --git a/offline/packages/TpcConditions/TpcConditionsReco.h b/offline/packages/TpcConditions/TpcConditionsReco.h index e3b38a4eff..d2b962ae22 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.h +++ b/offline/packages/TpcConditions/TpcConditionsReco.h @@ -28,6 +28,7 @@ class TpcConditionsReco : public SubsysReco float get_MedianCurrent(int channel, const std::vector &channels); void fillConditions(int channel); float get_AverageMedianCurrent(const std::vector &channels); + float get_InterpolatedMedianCurrent(uint64_t bco, const std::vector &channels); std::map m_bco_to_channel; }; diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc index efe626ca15..06e4a5ea5e 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc @@ -16,7 +16,9 @@ #include #include #include +#include #include +#include #include #include @@ -228,16 +230,27 @@ bool Tpc_PolyClusterizer::load_cdb_inputs() return " [CDB]"; }; - std::cout << "\033[31m" - << "Load NR1 " << m_conditions->get_LoadNR1() - << " Load SR1 " << m_conditions->get_LoadSR1() - << "\033[0m" << std::endl; - - if(true) + if (m_useBCOkEffs) { - m_kEffSide0 = m_kEffSide0 * m_conditions->get_LoadSR1()/m_conditions->get_AverageLoadSR1(); - m_kEffSide1 = m_kEffSide1 * m_conditions->get_LoadNR1()/m_conditions->get_AverageLoadNR1(); - } + const double averageSR1 = m_conditions->get_AverageLoadSR1(); + const double averageNR1 = m_conditions->get_AverageLoadNR1(); + + std::cout << Name() << "::load_cdb_inputs" + << " - SR1=" << m_conditions->get_LoadSR1() + << " avgSR1=" << averageSR1 + << " NR1=" << m_conditions->get_LoadNR1() + << " avgNR1=" << averageNR1 + << std::endl; + if (averageSR1 != 0.0 && averageNR1 != 0.0) + { + m_kEffSide0 *= m_conditions->get_LoadSR1() / averageSR1; + m_kEffSide1 *= m_conditions->get_LoadNR1() / averageNR1; + } + else + { + std::cout << Name() << "::load_cdb_inputs - warning: average SR1 or NR1 is zero, cannot apply BC correction" << std::endl; + } + } std::cout << Name() << "::load_cdb_inputs - final kEff values: side0 = " << m_kEffSide0 << keff_source(m_kEffSide0Override) << ", side1 = " << m_kEffSide1 << keff_source(m_kEffSide1Override) << std::endl; @@ -265,6 +278,13 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) // get layer geometry for layer 20. auto* layergeom = m_geomContainerTpc->GetLayerCellGeom(20); + if (!layergeom) + { + std::cout << Name() << "::InitRun - missing TPC geometry for layer 20" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + if (use_survey_geometry) { // apply survey geometry @@ -332,10 +352,27 @@ void Tpc_PolyClusterizer::configure_garfield(PHGarfield* garfield) const garfield->SetFrameChargeScale(m_frameChargeScale); - garfield->SetUseIFCVoltageDistortion(true); - garfield->SetUseOFCVoltageDistortion(true); + recoConsts* rc = recoConsts::instance(); + int runnumber = rc->get_IntFlag("RUNNUMBER"); + + if (runnumber > RunnumberRange::RUN3AUAU_IFC_V_CHANGE) + { + garfield->SetUseIFCVoltageDistortion(true); + garfield->SetUseOFCVoltageDistortion(true); - garfield->SetFieldCageVoltageOffsets(m_fieldCageVoltageOffsets[0], m_fieldCageVoltageOffsets[1], m_fieldCageVoltageOffsets[2], m_fieldCageVoltageOffsets[3]); + garfield->SetFieldCageVoltageOffsets(m_fieldCageVoltageOffsets[0], m_fieldCageVoltageOffsets[1], m_fieldCageVoltageOffsets[2], m_fieldCageVoltageOffsets[3]); + + std::cout << Name() << "::configure_garfield - using IFC and OFC voltage distortion with offsets: " + << "IFC South = " << m_fieldCageVoltageOffsets[0] << " V, " + << "IFC North = " << m_fieldCageVoltageOffsets[1] << " V, " + << "OFC South = " << m_fieldCageVoltageOffsets[2] << " V, " + << "OFC North = " << m_fieldCageVoltageOffsets[3] << " V" + << std::endl; + } + else + { + std::cout << Name() << "::configure_garfield - not using IFC and OFC voltage distortion for run number " << runnumber << std::endl; + } garfield->MoveTpc(m_tpcMove[0], m_tpcMove[1], m_tpcMove[2]); for (const auto& rotation : m_tpcRotations) @@ -375,11 +412,14 @@ int Tpc_PolyClusterizer::getNodes(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTRUN; } - m_conditions = findNode::getClass(topNode, "TpcConditions"); - if (!m_conditions) + if (m_useBCOkEffs) { - std::cerr << Name() << "::getNodes - missing TpcConditions" << std::endl; - return Fun4AllReturnCodes::ABORTRUN; + m_conditions = findNode::getClass(topNode, "TpcConditions"); + if (!m_conditions) + { + std::cout << Name() << "::getNodes - missing TpcConditions" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } } return Fun4AllReturnCodes::EVENT_OK; diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h index 6e4ea564be..009de3d48b 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h @@ -96,6 +96,7 @@ class Tpc_PolyClusterizer : public SubsysReco void setFrameChargeScale(double v) { m_frameChargeScale = v; } void setFieldCageVoltageOffsets(double ifcSouth, double ifcNorth, double ofcSouth, double ofcNorth) { m_fieldCageVoltageOffsets = {{ifcSouth, ifcNorth, ofcSouth, ofcNorth}}; } void setUse2DElectricFieldMap(bool v) { m_use2DElectricFieldMap = v; } + void setUseBCOkEffs(bool v) { m_useBCOkEffs = v; } private: struct Point @@ -210,6 +211,7 @@ class Tpc_PolyClusterizer : public SubsysReco bool m_framesSide0Override{false}; bool m_framesSide1Override{false}; bool m_use2DElectricFieldMap{false}; + bool m_useBCOkEffs{true}; std::array m_tpcMove{{0.0, 0.0, 0.0}}; //{{-0.16775, -0.0337, -0.71365}}; std::array, 2> m_tpcRotations{{{{0.0, 0.0, 0.0}}, {{0.0, 0.0, 0.0}}}}; //{{{{0.0, 0.01485 / 10.0, 0.0}}, {{0.0298 / 8.0, 0.0, 0.0}}}}; double m_frameChargeScale{-180.0}; From a472c7f60fe1fe9d3b735aca965598cb6110ce98 Mon Sep 17 00:00:00 2001 From: mcyoren Date: Wed, 16 Sep 2026 17:26:04 -0400 Subject: [PATCH 3/6] Move TPC drift configuration into PHGarfield --- offline/packages/PHGarfield/Makefile.am | 1 + offline/packages/PHGarfield/PHGarfield.cc | 357 +++++++++++++++++- offline/packages/PHGarfield/PHGarfield.h | 133 +++++-- .../TpcConditions/TpcConditionsReco.cc | 23 +- .../tpctrackreco/Tpc_PolyClusterizer.cc | 39 +- .../tpctrackreco/Tpc_PolyClusterizer.h | 101 ++++- 6 files changed, 595 insertions(+), 59 deletions(-) diff --git a/offline/packages/PHGarfield/Makefile.am b/offline/packages/PHGarfield/Makefile.am index cbaae76232..e4749a46f0 100644 --- a/offline/packages/PHGarfield/Makefile.am +++ b/offline/packages/PHGarfield/Makefile.am @@ -25,6 +25,7 @@ libPHGarfield_la_LIBADD = \ -lffamodules \ -lffarawobjects \ -lcdbobjects \ + -lg4detectors \ -lEvent \ -lphool \ -lphfield \ diff --git a/offline/packages/PHGarfield/PHGarfield.cc b/offline/packages/PHGarfield/PHGarfield.cc index 011d201e2b..c099587b28 100644 --- a/offline/packages/PHGarfield/PHGarfield.cc +++ b/offline/packages/PHGarfield/PHGarfield.cc @@ -1,12 +1,17 @@ #include "PHGarfield.h" -#include -#include +#include #include +#include #include #include - -#include +#include +#include +#include +#include +#include +#include +#include #include #include @@ -38,6 +43,16 @@ namespace fs = std::filesystem; +namespace +{ + constexpr const char* CdbEField3DSide0 = "Tpc_PolySeeding_EField3D_Side0"; + constexpr const char* CdbEField3DSide1 = "Tpc_PolySeeding_EField3D_Side1"; + constexpr const char* CdbModuleFrames3DSide0 = "Tpc_PolySeeding_ModuleFrames3D_Side0"; + constexpr const char* CdbModuleFrames3DSide1 = "Tpc_PolySeeding_ModuleFrames3D_Side1"; + constexpr const char* CdbKEff = "Tpc_PolyClusterizer_kEff"; + constexpr const char* CdbEField2D = "Tpc_PolySeeding_EField"; +} + PHGarfield::PHGarfield(const std::string& name, const std::string& electricFieldMap, double spaceChargeScale_side0, @@ -55,9 +70,15 @@ PHGarfield::PHGarfield(const std::string& name, : SubsysReco(name) , m_electricFieldMap(electricFieldMap) , m_electricFieldMap3D{{electricFieldMap3D_side0, electricFieldMap3D_side1}} - , m_spaceChargeScale_side0(spaceChargeScale_side0) - , m_spaceChargeScale_side1(spaceChargeScale_side1) + , m_spaceChargeScale_side0(std::isfinite(spaceChargeScale_side0) ? spaceChargeScale_side0 : 1.0) + , m_spaceChargeScale_side1(std::isfinite(spaceChargeScale_side1) ? spaceChargeScale_side1 : 1.0) { + m_electricFieldMapOverride = !electricFieldMap.empty(); + m_electricFieldMap3DOverride[0] = !electricFieldMap3D_side0.empty(); + m_electricFieldMap3DOverride[1] = !electricFieldMap3D_side1.empty(); + + m_spaceChargeScaleOverride[0] = std::isfinite(spaceChargeScale_side0); + m_spaceChargeScaleOverride[1] = std::isfinite(spaceChargeScale_side1); } PHGarfield::~PHGarfield() @@ -76,6 +97,313 @@ PHGarfield::~PHGarfield() ClearFrameElectricFieldCorrections3D(1); } +bool PHGarfield::LoadCDBInputs(PHCompositeNode* topNode) +{ + auto* cdb = CDBInterface::instance(); + + auto resolve_file = [this, cdb](const std::string& payload, std::string& filename, const bool manual_override) -> bool + { + if (manual_override) + { + std::cout << Name() << "::LoadCDBInputs - manual override for " << payload << ": " << filename << std::endl; + + if (filename.empty()) + { + std::cout << Name() << "::LoadCDBInputs - manual filename is empty for " << payload << std::endl; + return false; + } + + return true; + } + + filename = cdb->getUrl(payload); + + if (filename.empty()) + { + std::cout << Name() << "::LoadCDBInputs - CDB payload not found: " << payload << std::endl; + return false; + } + + std::cout << Name() << "::LoadCDBInputs - loaded " << payload << " from CDB: " << filename << std::endl; + + return true; + }; + + bool ok = true; + + if (use_survey_geometry && !m_tpcGeometryOverride) + { + auto* geomContainer = + findNode::getClass(topNode, "TPCGEOMCONTAINER"); + + if (!geomContainer) + { + std::cout << Name() + << "::LoadCDBInputs - missing TPCGEOMCONTAINER" + << std::endl; + return false; + } + + PHG4TpcGeom* layergeom = geomContainer->GetLayerCellGeom(20); + + if (!layergeom) + { + std::cout << Name() + << "::LoadCDBInputs - missing TPC geometry for layer 20" + << std::endl; + return false; + } + + m_tpcpos.SetXYZ(layergeom->get_place_x(), + layergeom->get_place_y(), + layergeom->get_place_z()); + + m_tpcrot = TRotation(); + m_tpcrot.RotateX(layergeom->get_rot_x()); + m_tpcrot.RotateY(layergeom->get_rot_y()); + m_tpcrot.RotateZ(layergeom->get_rot_z()); + } + + // Space-charge field maps. + // A manually selected 2D map takes precedence over automatic 3D CDB maps. + if (m_electricFieldMapOverride) + { + std::cout << Name() << "::LoadCDBInputs - manual 2D electric-field map: " << m_electricFieldMap << std::endl; + } + else if (m_use2DElectricFieldMap) + { + if (!resolve_file(CdbEField2D, m_electricFieldMap, false)) + { + ok = false; + } + } + else + { + if (!resolve_file(CdbEField3DSide0, m_electricFieldMap3D[0], m_electricFieldMap3DOverride[0])) + { + ok = false; + } + + if (!resolve_file(CdbEField3DSide1, m_electricFieldMap3D[1], m_electricFieldMap3DOverride[1])) + { + ok = false; + } + } + + // Frame-charge maps. + // A manually selected 2D frame map takes precedence over automatic 3D CDB maps. + if (m_frameElectricFieldMapOverride) + { + std::cout << Name() << "::LoadCDBInputs - manual 2D frame electric-field map: " + << m_frameElectricFieldMap << std::endl; + } + else + { + if (!resolve_file(CdbModuleFrames3DSide0, m_frameElectricFieldMap3D[0], m_frameElectricFieldMap3DOverride[0])) + { + ok = false; + } + + if (!resolve_file(CdbModuleFrames3DSide1, m_frameElectricFieldMap3D[1], m_frameElectricFieldMap3DOverride[1])) + { + ok = false; + } + } + + // kEff from CDB for any side that was not explicitly set. + if (!m_spaceChargeScaleOverride[0] || !m_spaceChargeScaleOverride[1]) + { + std::string kefffile; + + if (m_field3DCoefficientFileOverride) + { + kefffile = m_field3DCoefficientFile; + std::cout << Name() + << "::LoadCDBInputs - manual kEff coefficient file: " + << kefffile << std::endl; + } + else + { + kefffile = cdb->getUrl(CdbKEff); + std::cout << Name() + << "::LoadCDBInputs - kEff coefficient file from CDB: " + << kefffile << std::endl; + } + + if (kefffile.empty()) + { + std::cout << Name() << "::LoadCDBInputs - kEff coefficient file is empty" << std::endl; + ok = false; + } + else + { + CDBTTree keffcdbtree(kefffile); + keffcdbtree.LoadCalibrations(); + + if (!m_spaceChargeScaleOverride[0]) + { + m_spaceChargeScale_side0 = keffcdbtree.GetSingleFloatValue("keffside0"); + } + + if (!m_spaceChargeScaleOverride[1]) + { + m_spaceChargeScale_side1 = keffcdbtree.GetSingleFloatValue("keffside1"); + } + } + } + + // Print the base kEff source before the current correction. + std::cout << Name() << "::LoadCDBInputs - base kEff:" + << " side0=" << m_spaceChargeScale_side0 + << (m_spaceChargeScaleOverride[0] ? " [manual override]" : " [CDB]") + << ", side1=" << m_spaceChargeScale_side1 + << (m_spaceChargeScaleOverride[1] ? " [manual override]" : " [CDB]") + << std::endl; + + // Apply the same bunch-crossing/current correction used by Tpc_PolyClusterizer, + // but only if TpcConditions exists. + if (!m_useBCOkEffs) + { + std::cout << Name() + << "::LoadCDBInputs - TpcConditions current correction disabled" + << std::endl; + } + else + { + auto* conditions = findNode::getClass(topNode, "TpcConditions"); + + if (conditions) + { + const double averageSR1 = conditions->get_AverageLoadSR1(); + const double averageNR1 = conditions->get_AverageLoadNR1(); + + std::cout << Name() << "::LoadCDBInputs - TpcConditions:" + << " SR1=" << conditions->get_LoadSR1() + << " avgSR1=" << averageSR1 + << " NR1=" << conditions->get_LoadNR1() + << " avgNR1=" << averageNR1 + << std::endl; + + if (m_useBCOkEffs) + { + if (averageSR1 != 0.0) + { + m_spaceChargeScale_side0 *= conditions->get_LoadSR1() / averageSR1; + } + else + { + std::cout << Name() + << "::LoadCDBInputs - warning: average SR1 is zero, " + << "cannot apply side0 current correction" + << std::endl; + } + } + else + { + std::cout << Name() + << "::LoadCDBInputs - side0 kEff manually overridden; " + << "not applying TpcConditions correction" + << std::endl; + } + + if (m_useBCOkEffs) + { + if (averageNR1 != 0.0) + { + m_spaceChargeScale_side1 *= conditions->get_LoadNR1() / averageNR1; + } + else + { + std::cout << Name() + << "::LoadCDBInputs - warning: average NR1 is zero, " + << "cannot apply side1 current correction" + << std::endl; + } + } + else + { + std::cout << Name() + << "::LoadCDBInputs - side1 kEff manually overridden; " + << "not applying TpcConditions correction" + << std::endl; + } + } + else + { + std::cout << Name() + << "::LoadCDBInputs - TpcConditions not found; " + << "no current correction applied" + << std::endl; + } + } + std::cout << Name() << "::LoadCDBInputs - final kEff:" + << " side0=" << m_spaceChargeScale_side0 + << ", side1=" << m_spaceChargeScale_side1 + << std::endl; + + // Keep these prints even though frame scale currently has no CDB source. + // This makes an explicit SetFrameChargeScale(1.0) visible. + std::cout << Name() << "::LoadCDBInputs - frame charge scale:" + << " side0=" << m_frameChargeScale_side0 + << (m_frameChargeScaleOverride[0] ? " [manual override]" : " [default]") + << ", side1=" << m_frameChargeScale_side1 + << (m_frameChargeScaleOverride[1] ? " [manual override]" : " [default]") + << std::endl; + + return ok; +} + +void PHGarfield::ConfigureRunDependentFieldCage() +{ + const int runnumber = recoConsts::instance()->get_IntFlag("RUNNUMBER"); + + if (m_fieldCageConfigurationOverride) + { + std::cout << Name() + << "::ConfigureRunDependentFieldCage - manual override:" + << " useIFC=" << m_useIFCVoltageDistortion + << " useOFC=" << m_useOFCVoltageDistortion + << " IFC South=" << m_ifcVoltageOffset_side0 << " V" + << " IFC North=" << m_ifcVoltageOffset_side1 << " V" + << " OFC South=" << m_ofcVoltageOffset_side0 << " V" + << " OFC North=" << m_ofcVoltageOffset_side1 << " V" + << std::endl; + return; + } + + m_useIFCVoltageDistortion = false; + m_useOFCVoltageDistortion = false; + + m_ifcVoltageOffset_side0 = 0.0; + m_ifcVoltageOffset_side1 = 0.0; + m_ofcVoltageOffset_side0 = 0.0; + m_ofcVoltageOffset_side1 = 0.0; + + if (runnumber > RunnumberRange::RUN3AUAU_IFC_V_CHANGE) + { + m_useIFCVoltageDistortion = true; + m_useOFCVoltageDistortion = true; + + m_ifcVoltageOffset_side0 = 211.0; + + std::cout << Name() + << "::ConfigureRunDependentFieldCage - automatic configuration" + << " for run " << runnumber + << ": IFC South=" << m_ifcVoltageOffset_side0 << " V" + << " IFC North=" << m_ifcVoltageOffset_side1 << " V" + << " OFC South=" << m_ofcVoltageOffset_side0 << " V" + << " OFC North=" << m_ofcVoltageOffset_side1 << " V" + << std::endl; + } + else + { + std::cout << Name() + << "::ConfigureRunDependentFieldCage - IFC/OFC distortion disabled" + << " for run " << runnumber + << std::endl; + } +} + int PHGarfield::InitRun(PHCompositeNode *topNode) { if (Verbosity() > 1) @@ -100,6 +428,15 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) m_cdbTPCMAPttree = new CDBTTree(text); m_cdbTPCMAPttree->LoadCalibrations(); + if (!LoadCDBInputs(topNode)) + { + std::cout << PHWHERE << " Failed to load PHGarfield CDB inputs" + << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + ConfigureRunDependentFieldCage(); + // Load the optional axisymmetric space-charge correction map. // Failure is non-fatal: Garfield then uses only the nominal 400 V/cm field. if (!m_electricFieldMap.empty()) @@ -178,7 +515,11 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) gasfile = m_defaultGasfile; } InitializeGas(gasfile); - + if (!m_GasFilesLoaded) + { + std::cerr << PHWHERE << " Failed to initialize Garfield gas" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } // Diagnostic during code development... FillRadii(); if (Verbosity() > 1) @@ -331,6 +672,7 @@ void PHGarfield::RotateMagnet(double theta_x, double theta_y, double theta_z) void PHGarfield::MoveTpc(double x_cm, double y_cm, double z_cm) { + m_tpcGeometryOverride = true; m_tpcpos.SetXYZ(x_cm, y_cm, z_cm); if (Verbosity() > 0) { @@ -341,6 +683,7 @@ void PHGarfield::MoveTpc(double x_cm, double y_cm, double z_cm) void PHGarfield::RotateTpc(double theta_x, double theta_y, double theta_z) { + m_tpcGeometryOverride = true; m_tpcrot.RotateX(theta_x); m_tpcrot.RotateY(theta_y); m_tpcrot.RotateZ(theta_z); diff --git a/offline/packages/PHGarfield/PHGarfield.h b/offline/packages/PHGarfield/PHGarfield.h index 046e5d3185..a86191d5e7 100644 --- a/offline/packages/PHGarfield/PHGarfield.h +++ b/offline/packages/PHGarfield/PHGarfield.h @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -30,8 +31,8 @@ class PHGarfield : public SubsysReco // Backward-compatible constructor used by existing sPHENIX/Devon TPC code. PHGarfield(const std::string &name = "PHGarfield", const std::string &electricFieldMap = "", - double spaceChargeScale_side0 = 1.0, - double spaceChargeScale_side1 = 1.0); + double spaceChargeScale_side0 = std::numeric_limits::quiet_NaN(), + double spaceChargeScale_side1 = std::numeric_limits::quiet_NaN()); // Extended constructor for direct side-separated 3D field-map configuration. PHGarfield(const std::string &name, @@ -70,8 +71,9 @@ class PHGarfield : public SubsysReco void MoveTpc(double x_cm, double y_cm, double z_cm); void RotateTpc(double theta_x, double theta_y, double theta_z); void SetCMVoltageDefault(double voltage) { m_CMVoltageDefault = voltage; } + double GetCMVoltageDefault() const { return m_CMVoltageDefault; } + bool GetSpaceChargeScaleOverride() const { return m_spaceChargeScaleOverride[0] && m_spaceChargeScaleOverride[1]; } - // These are left in public namespace for easy plotting macros... // The user is encouraged to add more routine to fit their analysis goals... // Existing macros should call this one. Input and returned polyline are in @@ -87,7 +89,11 @@ class PHGarfield : public SubsysReco // Axisymmetric ROOT map must contain QA/hErDefault and QA/hEzDefault. // The histograms are expected in cm on the axes and V/m in the bins. - void SetElectricFieldMap(const std::string &filename) { m_electricFieldMap = filename; } + void SetElectricFieldMap(const std::string &filename) + { + m_electricFieldMap = filename; + m_electricFieldMapOverride = true; + } // Side-separated 3D ROOT maps must contain Field3D/hEx, Field3D/hEy, // and Field3D/hEz. Axes are (r [cm], phi [rad], |z| [cm]); bin contents @@ -96,38 +102,78 @@ class PHGarfield : public SubsysReco { m_electricFieldMap3D[0] = side0_filename; m_electricFieldMap3D[1] = side1_filename; + m_electricFieldMap3DOverride = {{true, true}}; + } + void SetElectricFieldMap3DSide0(const std::string &filename) + { + m_electricFieldMap3D[0] = filename; + m_electricFieldMap3DOverride[0] = true; } - void SetElectricFieldMap3DSide0(const std::string &filename) { m_electricFieldMap3D[0] = filename; } - void SetElectricFieldMap3DSide1(const std::string &filename) { m_electricFieldMap3D[1] = filename; } + void SetElectricFieldMap3DSide1(const std::string &filename) + { + m_electricFieldMap3D[1] = filename; + m_electricFieldMap3DOverride[1] = true; + } // Optional frame-charge correction maps are added on top of the existing // space-charge correction. The 2D format is QA/hErDefault + QA/hEzDefault; // the side-separated 3D format may be either Field3D/hEx + hEy + hEz // or the Rossegger cylindrical format hEr + hEphi + hEz at file root. - void SetFrameElectricFieldMap(const std::string &filename) { m_frameElectricFieldMap = filename; } + void SetFrameElectricFieldMap(const std::string &filename) + { + m_frameElectricFieldMap = filename; + m_frameElectricFieldMapOverride = true; + } + void SetFrameElectricFieldMap3D(const std::string &side0_filename, const std::string &side1_filename) { m_frameElectricFieldMap3D[0] = side0_filename; m_frameElectricFieldMap3D[1] = side1_filename; + m_frameElectricFieldMap3DOverride = {{true, true}}; } - void SetFrameElectricFieldMap3DSide0(const std::string &filename) { m_frameElectricFieldMap3D[0] = filename; } - void SetFrameElectricFieldMap3DSide1(const std::string &filename) { m_frameElectricFieldMap3D[1] = filename; } - void SetFrameChargeScale(double value) + + void SetFrameElectricFieldMap3DSide0(const std::string &filename) { - m_frameChargeScale_side0 = value; - m_frameChargeScale_side1 = value; + m_frameElectricFieldMap3D[0] = filename; + m_frameElectricFieldMap3DOverride[0] = true; + } + + void SetFrameElectricFieldMap3DSide1(const std::string &filename) + { + m_frameElectricFieldMap3D[1] = filename; + m_frameElectricFieldMap3DOverride[1] = true; } - void SetFrameChargeScaleSide0(double value) { m_frameChargeScale_side0 = value; } - void SetFrameChargeScaleSide1(double value) { m_frameChargeScale_side1 = value; } void SetSpaceChargeScale(double value) { m_spaceChargeScale_side0 = value; m_spaceChargeScale_side1 = value; + m_spaceChargeScaleOverride = {{true, true}}; } - void SetSpaceChargeScaleSide0(double value) { m_spaceChargeScale_side0 = value; } - void SetSpaceChargeScaleSide1(double value) { m_spaceChargeScale_side1 = value; } + void SetSpaceChargeScaleSide0(double value) { + m_spaceChargeScale_side0 = value; + m_spaceChargeScaleOverride[0] = true; + } + void SetSpaceChargeScaleSide1(double value) { + m_spaceChargeScale_side1 = value; + m_spaceChargeScaleOverride[1] = true; + } + + void SetFrameChargeScale(double value) + { + m_frameChargeScale_side0 = value; + m_frameChargeScale_side1 = value; + m_frameChargeScaleOverride = {{true, true}}; + } + void SetFrameChargeScaleSide0(double value) { + m_frameChargeScale_side0 = value; + m_frameChargeScaleOverride[0] = true; + } + void SetFrameChargeScaleSide1(double value) { + m_frameChargeScale_side1 = value; + m_frameChargeScaleOverride[1] = true; + } double GetSpaceChargeScaleSide0() const { return m_spaceChargeScale_side0; } double GetSpaceChargeScaleSide1() const { return m_spaceChargeScale_side1; } @@ -139,25 +185,30 @@ class PHGarfield : public SubsysReco // Offsets are endpoint perturbations relative to the nominal resistor-chain // boundary voltage, in volts. IFC and OFC are independent, giving four // tunable parameters: IFC South/North and OFC South/North. - void SetUseIFCVoltageDistortion(bool value) { m_useIFCVoltageDistortion = value; } + void SetUseIFCVoltageDistortion(bool value) { + m_useIFCVoltageDistortion = value; + m_fieldCageConfigurationOverride = true; + } void SetIFCVoltageOffset(double side0_south_v, double side1_north_v) { m_ifcVoltageOffset_side0 = side0_south_v; m_ifcVoltageOffset_side1 = side1_north_v; + m_fieldCageConfigurationOverride = true; } - void SetIFCVoltageOffsetSide0(double value_v) { m_ifcVoltageOffset_side0 = value_v; } - void SetIFCVoltageOffsetSide1(double value_v) { m_ifcVoltageOffset_side1 = value_v; } + void SetIFCVoltageOffsetSide0(double value_v) { m_ifcVoltageOffset_side0 = value_v; m_fieldCageConfigurationOverride = true; } + void SetIFCVoltageOffsetSide1(double value_v) { m_ifcVoltageOffset_side1 = value_v; m_fieldCageConfigurationOverride = true; } double GetIFCVoltageOffsetSide0() const { return m_ifcVoltageOffset_side0; } double GetIFCVoltageOffsetSide1() const { return m_ifcVoltageOffset_side1; } - void SetUseOFCVoltageDistortion(bool value) { m_useOFCVoltageDistortion = value; } + void SetUseOFCVoltageDistortion(bool value) { m_useOFCVoltageDistortion = value; m_fieldCageConfigurationOverride = true; } void SetOFCVoltageOffset(double side0_south_v, double side1_north_v) { m_ofcVoltageOffset_side0 = side0_south_v; m_ofcVoltageOffset_side1 = side1_north_v; + m_fieldCageConfigurationOverride = true; } - void SetOFCVoltageOffsetSide0(double value_v) { m_ofcVoltageOffset_side0 = value_v; } - void SetOFCVoltageOffsetSide1(double value_v) { m_ofcVoltageOffset_side1 = value_v; } + void SetOFCVoltageOffsetSide0(double value_v) { m_ofcVoltageOffset_side0 = value_v; m_fieldCageConfigurationOverride = true; } + void SetOFCVoltageOffsetSide1(double value_v) { m_ofcVoltageOffset_side1 = value_v; m_fieldCageConfigurationOverride = true; } double GetOFCVoltageOffsetSide0() const { return m_ofcVoltageOffset_side0; } double GetOFCVoltageOffsetSide1() const { return m_ofcVoltageOffset_side1; } @@ -183,6 +234,21 @@ class PHGarfield : public SubsysReco m_ifcGridNZ = nz; } + void SetUseSurveyGeometry(bool value) { use_survey_geometry = value; } + bool GetUseSurveyGeometry() const { return use_survey_geometry; } + + void SetUseBCOkEffs(bool value) { m_useBCOkEffs = value; } + bool GetUseBCOkEffs() const { return m_useBCOkEffs; } + + void SetUse2DElectricFieldMap(bool value) { m_use2DElectricFieldMap = value; } + bool GetUse2DElectricFieldMap() const { return m_use2DElectricFieldMap; } + + void SetField3DCoefficientFile(const std::string &filename) + { + m_field3DCoefficientFile = filename; + m_field3DCoefficientFileOverride = true; + } + private: void GetMagneticFieldTesla(double x_cm, double y_cm, double z_cm, double &bx_t, double &by_t, double &bz_t) const; // Feeds magnetic field to Garfield void GetElectricFieldVcm(double x_cm, double y_cm, double z_cm, double &ex_vcm, double &ey_vcm, double &ez_vcm) const; // Feeds electric field to Garfield @@ -212,6 +278,9 @@ class PHGarfield : public SubsysReco void FillRadii(); static double bounder(double phi, double phi_min); + bool LoadCDBInputs(PHCompositeNode *topNode); + void ConfigureRunDependentFieldCage(); + CDBTTree *m_cdbTPCMAPttree{nullptr}; // Locations of the pads from CDB... //PHField3DCartesian *m_field{nullptr}; // The standard sPHENIX field holding container. PHField *m_field{nullptr}; @@ -236,7 +305,7 @@ class PHGarfield : public SubsysReco std::array m_electricFieldMap3D{}; double m_spaceChargeScale_side0{1.0}; // south, z < 0 double m_spaceChargeScale_side1{1.0}; // north, z > 0 - double m_CMVoltageDefault{380.0}; // V/cm, nominal TPC field + double m_CMVoltageDefault{375.0}; // V/cm, nominal TPC field bool m_zerofield{false}; TH2 *m_erCorrection{nullptr}; // radial correction, input bins in V/m TH2 *m_ezCorrection{nullptr}; // local longitudinal correction, input bins in V/m @@ -248,8 +317,8 @@ class PHGarfield : public SubsysReco // over the optional axisymmetric frame map on that side. std::string m_frameElectricFieldMap; std::array m_frameElectricFieldMap3D{}; - double m_frameChargeScale_side0{1.0}; - double m_frameChargeScale_side1{1.0}; + double m_frameChargeScale_side0{-180.0}; + double m_frameChargeScale_side1{-180.0}; TH2 *m_frameErCorrection{nullptr}; TH2 *m_frameEzCorrection{nullptr}; std::array, 2> m_frameField3DCorrection{}; @@ -275,6 +344,20 @@ class PHGarfield : public SubsysReco double m_ifcGridDr_cm{0.0}; double m_ifcGridDz_cm{0.0}; bool m_fieldCageGridReady{false}; + bool m_tpcGeometryOverride{false}; + bool m_electricFieldMapOverride{false}; + std::array m_electricFieldMap3DOverride{{false, false}}; + bool m_frameElectricFieldMapOverride{false}; + std::array m_frameElectricFieldMap3DOverride{{false, false}}; + std::array m_spaceChargeScaleOverride{{false, false}}; + std::array m_frameChargeScaleOverride{{false, false}}; + bool m_fieldCageConfigurationOverride{false}; + bool use_survey_geometry{true}; + bool m_useBCOkEffs{true}; + bool m_use2DElectricFieldMap{false}; + std::string m_field3DCoefficientFile; + bool m_field3DCoefficientFileOverride{false}; + std::vector m_ifcUnitErGrid; // V/cm for +1 V IFC endpoint perturbation std::vector m_ifcUnitEsGrid; // V/cm along +|z| for +1 V IFC endpoint perturbation std::vector m_ofcUnitErGrid; // V/cm for +1 V OFC endpoint perturbation diff --git a/offline/packages/TpcConditions/TpcConditionsReco.cc b/offline/packages/TpcConditions/TpcConditionsReco.cc index ae5e4d3eee..d176b77140 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.cc +++ b/offline/packages/TpcConditions/TpcConditionsReco.cc @@ -263,10 +263,7 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) const uint64_t targetBco = static_cast(rawBco); - std::cout << "\033[1;31m" - << " targetBCO = " - << targetBco - << "\033[0m" << std::endl; + std::cout << " targetBCO = " << targetBco << std::endl; // Get the TPC conditions payload from CDB std::string calibdir = CDBInterface::instance()->getUrl("TPC_CONDITIONS"); @@ -339,8 +336,7 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) m_bco_to_channel.begin(), static_cast(m_bco_to_channel.size() / 2)); - std::cout << "\033[1;31m" - << "R1 conditions:" + std::cout << "R1 conditions:" << "\n FIRST BCO=" << first->first << " NR1=" << get_MedianCurrent(first->second, NR1) << " SR1=" << get_MedianCurrent(first->second, SR1) @@ -352,7 +348,7 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) << " SR1=" << get_MedianCurrent(last->second, SR1) << "\n AVERAGE NR1=" << m_conditions->get_AverageLoadNR1() << " SR1=" << m_conditions->get_AverageLoadSR1() - << "\033[0m" << std::endl; + << std::endl; fillConditions(selected->second); @@ -411,6 +407,19 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) m_conditions->set_FieldOK(m_tree->GetFloatValue(channel, "FieldOK") != 0.0); m_conditions->set_GainOK(m_tree->GetFloatValue(channel, "GainOK") != 0.0); + m_conditions->set_LoadCurrent(get_InterpolatedMedianCurrent(channel, ALL)); + m_conditions->set_LoadNorth(get_InterpolatedMedianCurrent(channel, NORTH)); + m_conditions->set_LoadSouth(get_InterpolatedMedianCurrent(channel, SOUTH)); + + m_conditions->set_LoadSR1(get_InterpolatedMedianCurrent(channel, SR1)); + m_conditions->set_LoadSR2(get_InterpolatedMedianCurrent(channel, SR2)); + m_conditions->set_LoadSR3(get_InterpolatedMedianCurrent(channel, SR3)); + + m_conditions->set_LoadNR1(get_InterpolatedMedianCurrent(channel, NR1)); + m_conditions->set_LoadNR2(get_InterpolatedMedianCurrent(channel, NR2)); + m_conditions->set_LoadNR3(get_InterpolatedMedianCurrent(channel, NR3)); + + // Do or die if (!m_conditions->get_FieldOK() || !m_conditions->get_GainOK()) diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc index 06e4a5ea5e..05b0b54d70 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc @@ -276,7 +276,7 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTRUN; } - // get layer geometry for layer 20. + // get layer geometry for layer 20. auto* layergeom = m_geomContainerTpc->GetLayerCellGeom(20); if (!layergeom) { @@ -285,7 +285,7 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTRUN; } - if (use_survey_geometry) + if (!m_usePHGarfieldDefaults && use_survey_geometry) { // apply survey geometry const double rot_x = layergeom->get_rot_x(); @@ -308,15 +308,40 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) std::cout << Name() << "::InitRun - m_startZSouth: " << m_startZSouth << " cm" << std::endl; std::cout << Name() << "::InitRun - m_startZNorth: " << m_startZNorth << " cm" << std::endl; - if (!load_cdb_inputs()) + if (!m_usePHGarfieldDefaults) { - std::cout << Name() << "::InitRun - failed to load CDB inputs" << std::endl; - return Fun4AllReturnCodes::ABORTRUN; + if (!load_cdb_inputs()) + { + std::cout << Name() << "::InitRun - failed to load CDB inputs" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } } - m_garfield.reset( new PHGarfield(Name() + "_PHGarfield", "", m_kEffSide0, m_kEffSide1) ); + delete m_garfield; + m_garfield = nullptr; + + if (m_usePHGarfieldDefaults) + { + std::cout << Name() + << "::InitRun - using PHGarfield default configuration" + << std::endl; + + m_garfield = new PHGarfield(Name() + "_PHGarfield"); - configure_garfield(m_garfield.get()); + reconfigure_garfield(m_garfield); + } + else + { + std::cout << Name() + << "::InitRun - using PolyClusterizer/manual PHGarfield configuration" + << std::endl; + + m_garfield = new PHGarfield(Name() + "_PHGarfield", "", m_kEffSide0, m_kEffSide1); + + m_garfield->SetUseSurveyGeometry(false); + configure_garfield(m_garfield); + } + if (m_garfield->InitRun(topNode) != Fun4AllReturnCodes::EVENT_OK) { std::cerr << Name() << "::InitRun - PHGarfield InitRun failed" << std::endl; diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h index 009de3d48b..9313358b8f 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h @@ -37,10 +37,26 @@ class Tpc_PolyClusterizer : public SubsysReco void setOutputNodeName(const std::string& n) { m_outputNodeName = n; } void setCrossingDecisionNodeName(const std::string& n) { m_crossingDecisionNodeName = n; } void setMaxAcceptedTier(unsigned char v) { m_maxAcceptedTier = v; } - void setT0(double v) { m_t0 = v; } - void setTpcAdcClock(double v) { m_tpcAdcClock = v; } - void setCrossingPeriodNs(double v) { m_crossingPeriodNs = v; } - void setReverseDriftStepNs(double v) { m_reverseDriftStepNs = v; } + void setT0(double v) + { + m_t0 = v; + } + + void setTpcAdcClock(double v) + { + m_tpcAdcClock = v; + } + + void setCrossingPeriodNs(double v) + { + m_crossingPeriodNs = v; + } + + void setReverseDriftStepNs(double v) + { + m_reverseDriftStepNs = v; + } + void setKEffSide0(double v) { m_kEffSide0 = v; @@ -52,51 +68,101 @@ class Tpc_PolyClusterizer : public SubsysReco m_kEffSide1 = v; m_kEffSide1Override = true; } - void setField3DCoefficientFile(const std::string& n) { + + void setField3DCoefficientFile(const std::string& n) + { m_field3DCoefficientFile = n; m_field3DCoefficientFileOverride = true; } + void setElectricFieldMap(const std::string& n) { m_electricFieldMap = n; m_electricFieldMapOverride = true; } + void setElectricFieldMap3DSide0(const std::string& n) { m_field3DSide0 = n; m_field3DSide0Override = true; } + void setElectricFieldMap3DSide1(const std::string& n) { m_field3DSide1 = n; m_field3DSide1Override = true; } + void setFrameElectricFieldMap3DSide0(const std::string& n) { m_framesSide0 = n; m_framesSide0Override = true; } + void setFrameElectricFieldMap3DSide1(const std::string& n) { m_framesSide1 = n; m_framesSide1Override = true; } - void setCMVoltageDefault(double v) { m_cmVoltageDefault = v; } - void setUseSurveyGeometry(bool v) { use_survey_geometry = v; } - void setMoveTpc(double x, double y, double z) { m_tpcMove = {{x, y, z}}; } + + void setCMVoltageDefault(double v) + { + m_cmVoltageDefault = v; + m_cmVoltageDefaultOverride = true; + } + + void setUseSurveyGeometry(bool v) + { + use_survey_geometry = v; + m_useSurveyGeometryOverride = true; + } + + void setMoveTpc(double x, double y, double z) + { + m_tpcMove = {{x, y, z}}; + m_tpcGeometryOverride = true; + } + void setRotateTpc(unsigned int index, double x, double y, double z) { - if (index < m_tpcRotations.size()) m_tpcRotations[index] = {{x, y, z}}; + if (index < m_tpcRotations.size()) + { + m_tpcRotations[index] = {{x, y, z}}; + m_tpcGeometryOverride = true; + } } + void setStartZ(double south_z, double north_z) { m_startZSouth = south_z; m_startZNorth = north_z; } - void setFrameChargeScale(double v) { m_frameChargeScale = v; } - void setFieldCageVoltageOffsets(double ifcSouth, double ifcNorth, double ofcSouth, double ofcNorth) { m_fieldCageVoltageOffsets = {{ifcSouth, ifcNorth, ofcSouth, ofcNorth}}; } - void setUse2DElectricFieldMap(bool v) { m_use2DElectricFieldMap = v; } - void setUseBCOkEffs(bool v) { m_useBCOkEffs = v; } + + void setFrameChargeScale(double v) + { + m_frameChargeScale = v; + m_frameChargeScaleOverride = true; + } + + void setFieldCageVoltageOffsets(double ifcSouth, double ifcNorth, + double ofcSouth, double ofcNorth) + { + m_fieldCageVoltageOffsets = {{ifcSouth, ifcNorth, ofcSouth, ofcNorth}}; + m_fieldCageVoltageOverride = true; + } + + void setUse2DElectricFieldMap(bool v) + { + m_use2DElectricFieldMap = v; + } + + void setUseBCOkEffs(bool v) + { + m_useBCOkEffs = v; + m_useBCOkEffsOverride = true; + } + + void setUsePHGarfieldDefaults(bool v) { m_usePHGarfieldDefaults = v; } private: struct Point @@ -162,6 +228,7 @@ class Tpc_PolyClusterizer : public SubsysReco ClusterParameters make_cluster_parameters(const std::vector& points, const Centroid& centroid, int side) const; static Centroid make_centroid(const std::vector& points); void configure_garfield(PHGarfield* garfield) const; + void reconfigure_garfield(PHGarfield* garfield) const; static unsigned int drift_lookup_index(unsigned int layer_index, unsigned int side, unsigned int sector, unsigned int sample); std::string m_inputNodeName; std::string m_outputNodeName; @@ -212,6 +279,14 @@ class Tpc_PolyClusterizer : public SubsysReco bool m_framesSide1Override{false}; bool m_use2DElectricFieldMap{false}; bool m_useBCOkEffs{true}; + bool m_usePHGarfieldDefaults{true}; + bool m_useSurveyGeometryOverride{false}; + bool m_cmVoltageDefaultOverride{false}; + bool m_tpcGeometryOverride{false}; + bool m_frameChargeScaleOverride{false}; + bool m_fieldCageVoltageOverride{false}; + bool m_useBCOkEffsOverride{false}; + std::array m_tpcMove{{0.0, 0.0, 0.0}}; //{{-0.16775, -0.0337, -0.71365}}; std::array, 2> m_tpcRotations{{{{0.0, 0.0, 0.0}}, {{0.0, 0.0, 0.0}}}}; //{{{{0.0, 0.01485 / 10.0, 0.0}}, {{0.0298 / 8.0, 0.0, 0.0}}}}; double m_frameChargeScale{-180.0}; From 8ef8f21b24adce891481495a6f8e71058769c3f2 Mon Sep 17 00:00:00 2001 From: mcyoren Date: Thu, 17 Sep 2026 15:28:33 -0400 Subject: [PATCH 4/6] Handle unavailable TPC conditions for kEff correction --- offline/packages/PHGarfield/PHGarfield.cc | 61 +++++----- .../packages/TpcConditions/TpcConditions.h | 4 + .../TpcConditions/TpcConditionsReco.cc | 105 ++++++++++++------ .../tpctrackreco/Tpc_PolyClusterizer.cc | 56 +++++++--- 4 files changed, 151 insertions(+), 75 deletions(-) diff --git a/offline/packages/PHGarfield/PHGarfield.cc b/offline/packages/PHGarfield/PHGarfield.cc index c099587b28..51563a6b07 100644 --- a/offline/packages/PHGarfield/PHGarfield.cc +++ b/offline/packages/PHGarfield/PHGarfield.cc @@ -5,6 +5,8 @@ #include #include #include +#include +#include #include #include #include @@ -272,7 +274,21 @@ bool PHGarfield::LoadCDBInputs(PHCompositeNode* topNode) { auto* conditions = findNode::getClass(topNode, "TpcConditions"); - if (conditions) + if (!conditions) + { + std::cout << Name() + << "::LoadCDBInputs - WARNING: TpcConditions node not found; " + << "using unscaled kEff" + << std::endl; + } + else if (!conditions->get_ConditionsAvailable()) + { + std::cout << Name() + << "::LoadCDBInputs - WARNING: TpcConditions are not available; " + << "using unscaled kEff" + << std::endl; + } + else { const double averageSR1 = conditions->get_AverageLoadSR1(); const double averageNR1 = conditions->get_AverageLoadNR1(); @@ -283,58 +299,41 @@ bool PHGarfield::LoadCDBInputs(PHCompositeNode* topNode) << " NR1=" << conditions->get_LoadNR1() << " avgNR1=" << averageNR1 << std::endl; - - if (m_useBCOkEffs) + + if (averageSR1 != 0.0 && averageNR1 != 0.0) { - if (averageSR1 != 0.0) + if (!m_spaceChargeScaleOverride[0]) { m_spaceChargeScale_side0 *= conditions->get_LoadSR1() / averageSR1; } else { std::cout << Name() - << "::LoadCDBInputs - warning: average SR1 is zero, " - << "cannot apply side0 current correction" + << "::LoadCDBInputs - side0 kEff manually overridden; " + << "not applying TpcConditions correction" << std::endl; } - } - else - { - std::cout << Name() - << "::LoadCDBInputs - side0 kEff manually overridden; " - << "not applying TpcConditions correction" - << std::endl; - } - if (m_useBCOkEffs) - { - if (averageNR1 != 0.0) + if (!m_spaceChargeScaleOverride[1]) { m_spaceChargeScale_side1 *= conditions->get_LoadNR1() / averageNR1; } else { std::cout << Name() - << "::LoadCDBInputs - warning: average NR1 is zero, " - << "cannot apply side1 current correction" + << "::LoadCDBInputs - side1 kEff manually overridden; " + << "not applying TpcConditions correction" << std::endl; } } else { std::cout << Name() - << "::LoadCDBInputs - side1 kEff manually overridden; " - << "not applying TpcConditions correction" + << "::LoadCDBInputs - WARNING: average SR1 or NR1 is zero, " + << "cannot apply current correction; using unscaled kEff" << std::endl; } } - else - { - std::cout << Name() - << "::LoadCDBInputs - TpcConditions not found; " - << "no current correction applied" - << std::endl; - } } std::cout << Name() << "::LoadCDBInputs - final kEff:" << " side0=" << m_spaceChargeScale_side0 @@ -414,6 +413,12 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) // Here we use the CDBInterface to set up the magnetic field map: m_field = PHFieldUtility::GetFieldMapNode(nullptr, topNode); + auto *fieldConfig = PHFieldUtility::GetFieldConfigNode(nullptr, topNode); + + std::cout << "Magnetic field map: " << fieldConfig->get_filename() + << "\nField type: " << fieldConfig->get_field_config_description() + << "\nField rescale: " << fieldConfig->get_magfield_rescale() + << std::endl; if (!m_field) { diff --git a/offline/packages/TpcConditions/TpcConditions.h b/offline/packages/TpcConditions/TpcConditions.h index 1e1cd67bad..f481b62b0c 100644 --- a/offline/packages/TpcConditions/TpcConditions.h +++ b/offline/packages/TpcConditions/TpcConditions.h @@ -73,6 +73,9 @@ class TpcConditions bool get_FieldOK() const { return m_FieldOK; } bool get_GainOK() const { return m_GainOK; } + void set_ConditionsAvailable(bool value) { m_ConditionsAvailable = value; } + bool get_ConditionsAvailable() const { return m_ConditionsAvailable; } + protected: float m_LoadCurrent{0.0}; float m_LoadNorth{0.0}; @@ -98,6 +101,7 @@ class TpcConditions float m_Pressure{0.0}; bool m_FieldOK{false}; bool m_GainOK{false}; + bool m_ConditionsAvailable{false}; }; #endif diff --git a/offline/packages/TpcConditions/TpcConditionsReco.cc b/offline/packages/TpcConditions/TpcConditionsReco.cc index d176b77140..28cc77f06c 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.cc +++ b/offline/packages/TpcConditions/TpcConditionsReco.cc @@ -239,6 +239,29 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) std::cout << "TpcConditionsReco::InitRun(PHCompositeNode *topNode) Initializing" << std::endl; + PHNodeIterator iter(topNode); + PHCompositeNode *parNode = dynamic_cast(iter.findFirst("PHCompositeNode", "PAR")); + + if (!parNode) + { + parNode = new PHCompositeNode("PAR"); + topNode->addNode(parNode); + } + + m_conditions = findNode::getClass(topNode, "TpcConditions"); + if (!m_conditions) + { + m_conditions = new TpcConditions(); + auto *conditionsNode = new PHDataNode(m_conditions, "TpcConditions", "Data"); + parNode->addNode(conditionsNode); + } + + m_conditions->set_ConditionsAvailable(false); + + delete m_tree; + m_tree = nullptr; + m_bco_to_channel.clear(); + Gl1Packet *gl1 = findNode::getClass(topNode, "GL1RAWHIT"); if (!gl1) { @@ -249,7 +272,8 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) { std::cout << "TpcConditionsReco::InitRun - could not find GL1 packet" << std::endl; - return Fun4AllReturnCodes::ABORTRUN; + m_conditions->set_ConditionsAvailable(false); + return Fun4AllReturnCodes::EVENT_OK; } const int64_t rawBco = gl1->lValue(0, "BCO"); @@ -258,7 +282,8 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) { std::cout << "TpcConditionsReco::InitRun - invalid GL1 BCO " << rawBco << std::endl; - return Fun4AllReturnCodes::ABORTRUN; + m_conditions->set_ConditionsAvailable(false); + return Fun4AllReturnCodes::EVENT_OK; } const uint64_t targetBco = static_cast(rawBco); @@ -267,6 +292,16 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) // Get the TPC conditions payload from CDB std::string calibdir = CDBInterface::instance()->getUrl("TPC_CONDITIONS"); + + if (calibdir.empty()) + { + std::cout << "TpcConditionsReco::InitRun - WARNING: " + << "TPC_CONDITIONS not available for this run; " + << "continuing without TPC conditions" + << std::endl; + return Fun4AllReturnCodes::EVENT_OK; + } + m_tree = new CDBTTree(calibdir); m_tree->LoadCalibrations(); @@ -280,21 +315,15 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) m_bco_to_channel[bco] = channel; } - // Find/create RUN node - PHNodeIterator iter(topNode); - PHCompositeNode *parNode = dynamic_cast(iter.findFirst("PHCompositeNode", "PAR")); - - if (!parNode) + if (m_bco_to_channel.empty()) { - parNode = new PHCompositeNode("PAR"); - topNode->addNode(parNode); + std::cout << "TpcConditionsReco::InitRun - WARNING: " + << "TPC_CONDITIONS contains no BCO entries; " + << "continuing without TPC conditions" + << std::endl; + return Fun4AllReturnCodes::EVENT_OK; } - // Create the transient TPC conditions object - m_conditions = new TpcConditions(); - PHDataNode *conditionsNode = new PHDataNode(m_conditions, "TpcConditions", "Data"); - parNode->addNode(conditionsNode); - m_conditions->set_AverageLoadCurrent(get_AverageMedianCurrent(ALL)); m_conditions->set_AverageLoadNorth(get_AverageMedianCurrent(NORTH)); m_conditions->set_AverageLoadSouth(get_AverageMedianCurrent(SOUTH)); @@ -373,11 +402,18 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) << std::endl; } + m_conditions->set_ConditionsAvailable(true); return Fun4AllReturnCodes::EVENT_OK; } int TpcConditionsReco::process_event(PHCompositeNode *topNode) { + + if (!m_tree || m_bco_to_channel.empty()) + { + return Fun4AllReturnCodes::EVENT_OK; + } + // Get the current event BCO Gl1Packet *gl1 = findNode::getClass(topNode, "GL1RAWHIT"); if (!gl1) @@ -385,21 +421,23 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) std::cout << "TpcConditionsReco::process_event - " << "could not find GL1RAWHIT node" << std::endl; - return Fun4AllReturnCodes::ABORTEVENT; + m_conditions->set_ConditionsAvailable(false); + return Fun4AllReturnCodes::EVENT_OK; } uint64_t bco = gl1->getBCO(); auto iter = m_bco_to_channel.upper_bound(bco); + if (iter == m_bco_to_channel.begin()) { - std::cout << "TpcConditionsReco::process_event - " - << "no TPC conditions available for BCO " - << bco << std::endl; - return Fun4AllReturnCodes::ABORTEVENT; + iter = m_bco_to_channel.begin(); } - - --iter; + else + { + --iter; + } + int channel = iter->second; m_conditions->set_Temperature(m_tree->GetFloatValue(channel, "gas_temperature")); @@ -407,22 +445,25 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) m_conditions->set_FieldOK(m_tree->GetFloatValue(channel, "FieldOK") != 0.0); m_conditions->set_GainOK(m_tree->GetFloatValue(channel, "GainOK") != 0.0); - m_conditions->set_LoadCurrent(get_InterpolatedMedianCurrent(channel, ALL)); - m_conditions->set_LoadNorth(get_InterpolatedMedianCurrent(channel, NORTH)); - m_conditions->set_LoadSouth(get_InterpolatedMedianCurrent(channel, SOUTH)); - m_conditions->set_LoadSR1(get_InterpolatedMedianCurrent(channel, SR1)); - m_conditions->set_LoadSR2(get_InterpolatedMedianCurrent(channel, SR2)); - m_conditions->set_LoadSR3(get_InterpolatedMedianCurrent(channel, SR3)); + m_conditions->set_LoadCurrent(get_InterpolatedMedianCurrent(bco, ALL)); + m_conditions->set_LoadNorth(get_InterpolatedMedianCurrent(bco, NORTH)); + m_conditions->set_LoadSouth(get_InterpolatedMedianCurrent(bco, SOUTH)); + + m_conditions->set_LoadSR1(get_InterpolatedMedianCurrent(bco, SR1)); + m_conditions->set_LoadSR2(get_InterpolatedMedianCurrent(bco, SR2)); + m_conditions->set_LoadSR3(get_InterpolatedMedianCurrent(bco, SR3)); - m_conditions->set_LoadNR1(get_InterpolatedMedianCurrent(channel, NR1)); - m_conditions->set_LoadNR2(get_InterpolatedMedianCurrent(channel, NR2)); - m_conditions->set_LoadNR3(get_InterpolatedMedianCurrent(channel, NR3)); + m_conditions->set_LoadNR1(get_InterpolatedMedianCurrent(bco, NR1)); + m_conditions->set_LoadNR2(get_InterpolatedMedianCurrent(bco, NR2)); + m_conditions->set_LoadNR3(get_InterpolatedMedianCurrent(bco, NR3)); + m_conditions->set_ConditionsAvailable(true); // Do or die - if (!m_conditions->get_FieldOK() || - !m_conditions->get_GainOK()) + if( m_conditions->get_ConditionsAvailable() && + (!m_conditions->get_FieldOK() || + !m_conditions->get_GainOK())) { return Fun4AllReturnCodes::ABORTEVENT; } diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc index 05b0b54d70..0f9f62cd9b 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc @@ -232,23 +232,40 @@ bool Tpc_PolyClusterizer::load_cdb_inputs() if (m_useBCOkEffs) { - const double averageSR1 = m_conditions->get_AverageLoadSR1(); - const double averageNR1 = m_conditions->get_AverageLoadNR1(); - - std::cout << Name() << "::load_cdb_inputs" - << " - SR1=" << m_conditions->get_LoadSR1() - << " avgSR1=" << averageSR1 - << " NR1=" << m_conditions->get_LoadNR1() - << " avgNR1=" << averageNR1 - << std::endl; - if (averageSR1 != 0.0 && averageNR1 != 0.0) + if (!m_conditions) + { + std::cout << Name() + << "::load_cdb_inputs - WARNING: TpcConditions node not found; " + << "using unscaled kEff" + << std::endl; + } + else if (!m_conditions->get_ConditionsAvailable()) { - m_kEffSide0 *= m_conditions->get_LoadSR1() / averageSR1; - m_kEffSide1 *= m_conditions->get_LoadNR1() / averageNR1; + std::cout << Name() + << "::load_cdb_inputs - WARNING: TpcConditions are not available; " + << "using unscaled kEff" + << std::endl; } else { - std::cout << Name() << "::load_cdb_inputs - warning: average SR1 or NR1 is zero, cannot apply BC correction" << std::endl; + const double averageSR1 = m_conditions->get_AverageLoadSR1(); + const double averageNR1 = m_conditions->get_AverageLoadNR1(); + + std::cout << Name() << "::load_cdb_inputs" + << " - SR1=" << m_conditions->get_LoadSR1() + << " avgSR1=" << averageSR1 + << " NR1=" << m_conditions->get_LoadNR1() + << " avgNR1=" << averageNR1 + << std::endl; + if (averageSR1 != 0.0 && averageNR1 != 0.0) + { + m_kEffSide0 *= m_conditions->get_LoadSR1() / averageSR1; + m_kEffSide1 *= m_conditions->get_LoadNR1() / averageNR1; + } + else + { + std::cout << Name() << "::load_cdb_inputs - warning: average SR1 or NR1 is zero, cannot apply BC correction" << std::endl; + } } } @@ -442,8 +459,17 @@ int Tpc_PolyClusterizer::getNodes(PHCompositeNode* topNode) m_conditions = findNode::getClass(topNode, "TpcConditions"); if (!m_conditions) { - std::cout << Name() << "::getNodes - missing TpcConditions" << std::endl; - return Fun4AllReturnCodes::ABORTRUN; + std::cout << Name() + << "::getNodes - WARNING: TpcConditions node not found; " + << "continuing with unscaled kEff" + << std::endl; + } + else if (!m_conditions->get_ConditionsAvailable()) + { + std::cout << Name() + << "::getNodes - WARNING: TpcConditions are not available; " + << "continuing with unscaled kEff" + << std::endl; } } From 487b276e3c7f4c7fbac47aed70203582a26a2709 Mon Sep 17 00:00:00 2001 From: mcyoren Date: Fri, 18 Sep 2026 14:32:05 -0400 Subject: [PATCH 5/6] Fix Tpc_PolyClusterizer after rebase --- .../tpctrackreco/Tpc_PolyClusterizer.cc | 162 +++++++++++++++++- 1 file changed, 153 insertions(+), 9 deletions(-) diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc index 0f9f62cd9b..bbba2972c9 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc @@ -293,7 +293,7 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTRUN; } - // get layer geometry for layer 20. + // get layer geometry for layer 20. auto* layergeom = m_geomContainerTpc->GetLayerCellGeom(20); if (!layergeom) { @@ -334,8 +334,7 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) } } - delete m_garfield; - m_garfield = nullptr; + m_garfield.reset(); if (m_usePHGarfieldDefaults) { @@ -343,9 +342,9 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) << "::InitRun - using PHGarfield default configuration" << std::endl; - m_garfield = new PHGarfield(Name() + "_PHGarfield"); + m_garfield = std::make_unique(Name() + "_PHGarfield"); - reconfigure_garfield(m_garfield); + reconfigure_garfield(m_garfield.get()); } else { @@ -353,12 +352,13 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) << "::InitRun - using PolyClusterizer/manual PHGarfield configuration" << std::endl; - m_garfield = new PHGarfield(Name() + "_PHGarfield", "", m_kEffSide0, m_kEffSide1); + m_garfield = std::make_unique( + Name() + "_PHGarfield", "", m_kEffSide0, m_kEffSide1); m_garfield->SetUseSurveyGeometry(false); - configure_garfield(m_garfield); + configure_garfield(m_garfield.get()); } - + if (m_garfield->InitRun(topNode) != Fun4AllReturnCodes::EVENT_OK) { std::cerr << Name() << "::InitRun - PHGarfield InitRun failed" << std::endl; @@ -374,6 +374,150 @@ int Tpc_PolyClusterizer::InitRun(PHCompositeNode* topNode) return Fun4AllReturnCodes::EVENT_OK; } +void Tpc_PolyClusterizer::reconfigure_garfield(PHGarfield* garfield) const +{ + if (!garfield) + { + return; + } + + if (garfield->GetCMVoltageDefault() != m_cmVoltageDefault) + { + std::cout << Name() << "::reconfigure_garfield - CM voltage " + << garfield->GetCMVoltageDefault() << " -> " + << m_cmVoltageDefault << " V/cm" << std::endl; + garfield->SetCMVoltageDefault(m_cmVoltageDefault); + } + + if (garfield->GetUseSurveyGeometry() != use_survey_geometry) + { + std::cout << Name() << "::reconfigure_garfield - survey geometry " + << garfield->GetUseSurveyGeometry() << " -> " + << use_survey_geometry << std::endl; + garfield->SetUseSurveyGeometry(use_survey_geometry); + } + + if (garfield->GetUseBCOkEffs() != m_useBCOkEffs) + { + std::cout << Name() << "::reconfigure_garfield - BCO/current correction " + << garfield->GetUseBCOkEffs() << " -> " + << m_useBCOkEffs << std::endl; + garfield->SetUseBCOkEffs(m_useBCOkEffs); + } + + if (garfield->GetUse2DElectricFieldMap() != m_use2DElectricFieldMap) + { + std::cout << Name() << "::reconfigure_garfield - use 2D electric-field map " + << garfield->GetUse2DElectricFieldMap() << " -> " + << m_use2DElectricFieldMap << std::endl; + garfield->SetUse2DElectricFieldMap(m_use2DElectricFieldMap); + } + + if (m_field3DCoefficientFileOverride) + { + std::cout << Name() << "::reconfigure_garfield - manual kEff coefficient file: " + << m_field3DCoefficientFile << std::endl; + garfield->SetField3DCoefficientFile(m_field3DCoefficientFile); + } + + if (m_electricFieldMapOverride) + { + std::cout << Name() << "::reconfigure_garfield - manual 2D electric-field map: " + << m_electricFieldMap << std::endl; + garfield->SetElectricFieldMap(m_electricFieldMap); + } + + if (m_kEffSide0Override) + { + std::cout << Name() << "::reconfigure_garfield - manual kEff side0: " + << m_kEffSide0 << std::endl; + garfield->SetSpaceChargeScaleSide0(m_kEffSide0); + } + + if (m_kEffSide1Override) + { + std::cout << Name() << "::reconfigure_garfield - manual kEff side1: " + << m_kEffSide1 << std::endl; + garfield->SetSpaceChargeScaleSide1(m_kEffSide1); + } + + if (m_field3DSide0Override) + { + std::cout << Name() << "::reconfigure_garfield - manual 3D field map side0: " + << m_field3DSide0 << std::endl; + garfield->SetElectricFieldMap3DSide0(m_field3DSide0); + } + + if (m_field3DSide1Override) + { + std::cout << Name() << "::reconfigure_garfield - manual 3D field map side1: " + << m_field3DSide1 << std::endl; + garfield->SetElectricFieldMap3DSide1(m_field3DSide1); + } + + if (m_framesSide0Override) + { + std::cout << Name() << "::reconfigure_garfield - manual frame map side0: " + << m_framesSide0 << std::endl; + garfield->SetFrameElectricFieldMap3DSide0(m_framesSide0); + } + + if (m_framesSide1Override) + { + std::cout << Name() << "::reconfigure_garfield - manual frame map side1: " + << m_framesSide1 << std::endl; + garfield->SetFrameElectricFieldMap3DSide1(m_framesSide1); + } + + if (m_frameChargeScaleOverride) + { + std::cout << Name() << "::reconfigure_garfield - manual frame charge scale: " + << m_frameChargeScale << std::endl; + garfield->SetFrameChargeScale(m_frameChargeScale); + } + + if (m_tpcGeometryOverride) + { + std::cout << Name() << "::reconfigure_garfield - manual TPC translation: " + << "(" << m_tpcMove[0] << ", " + << m_tpcMove[1] << ", " + << m_tpcMove[2] << ") cm" << std::endl; + + garfield->MoveTpc(m_tpcMove[0], m_tpcMove[1], m_tpcMove[2]); + + for (std::size_t i = 0; i < m_tpcRotations.size(); ++i) + { + const auto& rotation = m_tpcRotations[i]; + + std::cout << Name() << "::reconfigure_garfield - manual TPC rotation " + << i << ": (" + << rotation[0] << ", " + << rotation[1] << ", " + << rotation[2] << ") rad" << std::endl; + + garfield->RotateTpc(rotation[0], rotation[1], rotation[2]); + } + } + + if (m_fieldCageVoltageOverride) + { + std::cout << Name() << "::reconfigure_garfield - manual field-cage offsets:" + << " IFC South=" << m_fieldCageVoltageOffsets[0] + << " V, IFC North=" << m_fieldCageVoltageOffsets[1] + << " V, OFC South=" << m_fieldCageVoltageOffsets[2] + << " V, OFC North=" << m_fieldCageVoltageOffsets[3] + << " V" << std::endl; + + garfield->SetUseIFCVoltageDistortion(true); + garfield->SetUseOFCVoltageDistortion(true); + garfield->SetFieldCageVoltageOffsets( + m_fieldCageVoltageOffsets[0], + m_fieldCageVoltageOffsets[1], + m_fieldCageVoltageOffsets[2], + m_fieldCageVoltageOffsets[3]); + } +} + void Tpc_PolyClusterizer::configure_garfield(PHGarfield* garfield) const { if (!garfield) @@ -454,7 +598,7 @@ int Tpc_PolyClusterizer::getNodes(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTRUN; } - if (m_useBCOkEffs) + if (!m_usePHGarfieldDefaults && m_useBCOkEffs) { m_conditions = findNode::getClass(topNode, "TpcConditions"); if (!m_conditions) From f0d8c9cefaa22ccdd6203c1956d0889cccbd1b46 Mon Sep 17 00:00:00 2001 From: mcyoren Date: Fri, 18 Sep 2026 15:20:21 -0400 Subject: [PATCH 6/6] Clean up formatting and clang-tidy warnings --- offline/packages/PHGarfield/PHGarfield.cc | 112 +++++++++--------- offline/packages/PHGarfield/PHGarfield.h | 71 +++++++---- .../packages/TpcConditions/TpcConditions.h | 4 +- .../TpcConditions/TpcConditionsReco.cc | 10 +- .../tpctrackreco/Tpc_PolyClusterizer.cc | 2 +- .../tpctrackreco/Tpc_PolyClusterizer.h | 5 +- 6 files changed, 111 insertions(+), 93 deletions(-) diff --git a/offline/packages/PHGarfield/PHGarfield.cc b/offline/packages/PHGarfield/PHGarfield.cc index 51563a6b07..288824197a 100644 --- a/offline/packages/PHGarfield/PHGarfield.cc +++ b/offline/packages/PHGarfield/PHGarfield.cc @@ -3,8 +3,9 @@ #include #include #include +#include +#include #include -#include #include #include #include @@ -12,8 +13,6 @@ #include #include #include -#include -#include #include #include @@ -53,7 +52,7 @@ namespace constexpr const char* CdbModuleFrames3DSide1 = "Tpc_PolySeeding_ModuleFrames3D_Side1"; constexpr const char* CdbKEff = "Tpc_PolyClusterizer_kEff"; constexpr const char* CdbEField2D = "Tpc_PolySeeding_EField"; -} +} // namespace PHGarfield::PHGarfield(const std::string& name, const std::string& electricFieldMap, @@ -74,13 +73,10 @@ PHGarfield::PHGarfield(const std::string& name, , m_electricFieldMap3D{{electricFieldMap3D_side0, electricFieldMap3D_side1}} , m_spaceChargeScale_side0(std::isfinite(spaceChargeScale_side0) ? spaceChargeScale_side0 : 1.0) , m_spaceChargeScale_side1(std::isfinite(spaceChargeScale_side1) ? spaceChargeScale_side1 : 1.0) + , m_electricFieldMapOverride(!electricFieldMap.empty()) + , m_electricFieldMap3DOverride{{!electricFieldMap3D_side0.empty(), !electricFieldMap3D_side1.empty()}} + , m_spaceChargeScaleOverride{{std::isfinite(spaceChargeScale_side0), std::isfinite(spaceChargeScale_side1)}} { - m_electricFieldMapOverride = !electricFieldMap.empty(); - m_electricFieldMap3DOverride[0] = !electricFieldMap3D_side0.empty(); - m_electricFieldMap3DOverride[1] = !electricFieldMap3D_side1.empty(); - - m_spaceChargeScaleOverride[0] = std::isfinite(spaceChargeScale_side0); - m_spaceChargeScaleOverride[1] = std::isfinite(spaceChargeScale_side1); } PHGarfield::~PHGarfield() @@ -299,7 +295,7 @@ bool PHGarfield::LoadCDBInputs(PHCompositeNode* topNode) << " NR1=" << conditions->get_LoadNR1() << " avgNR1=" << averageNR1 << std::endl; - + if (averageSR1 != 0.0 && averageNR1 != 0.0) { if (!m_spaceChargeScaleOverride[0]) @@ -403,7 +399,7 @@ void PHGarfield::ConfigureRunDependentFieldCage() } } -int PHGarfield::InitRun(PHCompositeNode *topNode) +int PHGarfield::InitRun(PHCompositeNode* topNode) { if (Verbosity() > 1) { @@ -413,7 +409,7 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) // Here we use the CDBInterface to set up the magnetic field map: m_field = PHFieldUtility::GetFieldMapNode(nullptr, topNode); - auto *fieldConfig = PHFieldUtility::GetFieldConfigNode(nullptr, topNode); + auto* fieldConfig = PHFieldUtility::GetFieldConfigNode(nullptr, topNode); std::cout << "Magnetic field map: " << fieldConfig->get_filename() << "\nField type: " << fieldConfig->get_field_config_description() @@ -423,8 +419,8 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) if (!m_field) { std::cerr << PHWHERE - << " ERROR: magnetic field map not available" - << std::endl; + << " ERROR: magnetic field map not available" + << std::endl; return Fun4AllReturnCodes::ABORTRUN; } @@ -433,7 +429,7 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) m_cdbTPCMAPttree = new CDBTTree(text); m_cdbTPCMAPttree->LoadCalibrations(); - if (!LoadCDBInputs(topNode)) + if (!LoadCDBInputs(topNode)) { std::cout << PHWHERE << " Failed to load PHGarfield CDB inputs" << std::endl; @@ -441,7 +437,7 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) } ConfigureRunDependentFieldCage(); - + // Load the optional axisymmetric space-charge correction map. // Failure is non-fatal: Garfield then uses only the nominal 400 V/cm field. if (!m_electricFieldMap.empty()) @@ -449,7 +445,7 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) if (!LoadElectricFieldCorrections(m_electricFieldMap)) { std::cout << PHWHERE << " Failed to load electric-field correction map: " - << m_electricFieldMap << std::endl; + << m_electricFieldMap << std::endl; } } @@ -506,9 +502,9 @@ int PHGarfield::InitRun(PHCompositeNode *topNode) // Make the Garfield Component and register the methods that will interface to our fields... m_component = new Garfield::ComponentUser(); m_component->SetMagneticField([this](double x, double y, double z, double& bx, double& by, double& bz) - { GetMagneticFieldTesla(x, y, z, bx, by, bz); }); + { GetMagneticFieldTesla(x, y, z, bx, by, bz); }); m_component->SetElectricField([this](double x, double y, double z, double& ex, double& ey, double& ez) - { GetElectricFieldVcm(x, y, z, ex, ey, ez); }); + { GetElectricFieldVcm(x, y, z, ex, ey, ez); }); // Here we fetch the gas from the CDB std::string gasfile = m_cdb->getUrl("PHGARFIELD_GAS"); @@ -573,18 +569,18 @@ void PHGarfield::PrintGarfield(double x, double y, double z) const GetMagneticFieldTesla(x, y, z, bx, by, bz); m_gas->ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz); std::cout << " x:" << x - << " y:" << y - << " z:" << z - << " ex:" << ex - << " ey:" << ey - << " ez:" << ez - << " bx:" << bx - << " by:" << by - << " bz:" << bz - << " vx:" << vx - << " vy:" << vy - << " vz:" << vz - << std::endl; + << " y:" << y + << " z:" << z + << " ex:" << ex + << " ey:" << ey + << " ez:" << ez + << " bx:" << bx + << " by:" << by + << " bz:" << bz + << " vx:" << vx + << " vy:" << vy + << " vz:" << vz + << std::endl; } void PHGarfield::PrintGasSummary() const @@ -659,7 +655,7 @@ void PHGarfield::MoveMagnet(double x_cm, double y_cm, double z_cm) if (Verbosity() > 0) { std::cout << "PHGarfield: magnetic-field map translation = (" - << x_cm << ", " << y_cm << ", " << z_cm << ") cm" << std::endl; + << x_cm << ", " << y_cm << ", " << z_cm << ") cm" << std::endl; } } @@ -671,7 +667,7 @@ void PHGarfield::RotateMagnet(double theta_x, double theta_y, double theta_z) if (Verbosity() > 0) { std::cout << "PHGarfield: magnetic-field map rotation increment = (" - << theta_x << ", " << theta_y << ", " << theta_z << ") rad" << std::endl; + << theta_x << ", " << theta_y << ", " << theta_z << ") rad" << std::endl; } } @@ -682,7 +678,7 @@ void PHGarfield::MoveTpc(double x_cm, double y_cm, double z_cm) if (Verbosity() > 0) { std::cout << "PHGarfield: TPC translation = (" - << x_cm << ", " << y_cm << ", " << z_cm << ") cm" << std::endl; + << x_cm << ", " << y_cm << ", " << z_cm << ") cm" << std::endl; } } @@ -695,7 +691,7 @@ void PHGarfield::RotateTpc(double theta_x, double theta_y, double theta_z) if (Verbosity() > 0) { std::cout << "PHGarfield: TPC rotation increment = (" - << theta_x << ", " << theta_y << ", " << theta_z << ") rad" << std::endl; + << theta_x << ", " << theta_y << ", " << theta_z << ") rad" << std::endl; } } @@ -743,11 +739,11 @@ void PHGarfield::GetMagneticFieldTesla(double x_cm, double y_cm, double z_cm, do const TVector3 p_map_cm = TpcPointToMagnetFieldMapPoint(x_cm, y_cm, z_cm); double point[4] = - { - p_map_cm.X() * CLHEP::cm, - p_map_cm.Y() * CLHEP::cm, - p_map_cm.Z() * CLHEP::cm, - 0.0}; + { + p_map_cm.X() * CLHEP::cm, + p_map_cm.Y() * CLHEP::cm, + p_map_cm.Z() * CLHEP::cm, + 0.0}; double bfield_map[3] = {0.0, 0.0, 0.0}; @@ -823,9 +819,9 @@ void PHGarfield::GetElectricFieldVcm(double x_cm, double y_cm, double z_cm, doub } const double delta_er_vcm = spaceChargeScale * - InterpolateCorrectionVcm(m_erCorrection, r_cm, abs_z_cm); + InterpolateCorrectionVcm(m_erCorrection, r_cm, abs_z_cm); const double delta_ez_local_vcm = spaceChargeScale * - InterpolateCorrectionVcm(m_ezCorrection, r_cm, abs_z_cm); + InterpolateCorrectionVcm(m_ezCorrection, r_cm, abs_z_cm); // Convert the cylindrical radial correction to Cartesian components. if (r_cm > 0.0) @@ -1095,7 +1091,7 @@ bool PHGarfield::LoadElectricFieldCorrections(const std::string& filename) if (!input || input->IsZombie()) { std::cout << PHWHERE << " Could not open electric-field map: " - << filename << std::endl; + << filename << std::endl; return false; } @@ -1115,8 +1111,8 @@ bool PHGarfield::LoadElectricFieldCorrections(const std::string& filename) if (!er || !ez) { std::cout << PHWHERE - << " Missing QA/hErDefault or QA/hEzDefault in " - << filename << std::endl; + << " Missing QA/hErDefault or QA/hEzDefault in " + << filename << std::endl; return false; } @@ -1138,17 +1134,17 @@ bool PHGarfield::LoadElectricFieldCorrections(const std::string& filename) m_ezCorrection->SetDirectory(nullptr); std::cout << "Loaded axisymmetric electric-field corrections from " - << filename << std::endl; + << filename << std::endl; std::cout << " scale k_eff side0/south/z<0 = " - << m_spaceChargeScale_side0 << std::endl; + << m_spaceChargeScale_side0 << std::endl; std::cout << " scale k_eff side1/north/z>0 = " - << m_spaceChargeScale_side1 << std::endl; + << m_spaceChargeScale_side1 << std::endl; std::cout << " r range [cm] = [" - << m_erCorrection->GetXaxis()->GetXmin() << ", " - << m_erCorrection->GetXaxis()->GetXmax() << "]" << std::endl; + << m_erCorrection->GetXaxis()->GetXmin() << ", " + << m_erCorrection->GetXaxis()->GetXmax() << "]" << std::endl; std::cout << " |z| range [cm] = [" - << m_erCorrection->GetYaxis()->GetXmin() << ", " - << m_erCorrection->GetYaxis()->GetXmax() << "]" << std::endl; + << m_erCorrection->GetYaxis()->GetXmin() << ", " + << m_erCorrection->GetYaxis()->GetXmax() << "]" << std::endl; return true; } @@ -1181,12 +1177,12 @@ double PHGarfield::InterpolateCorrectionVcm( constexpr double epsilon = 1.0e-6; const double r_eval = - std::clamp(r_cm, r_min + epsilon, r_max - epsilon); + std::clamp(r_cm, r_min + epsilon, r_max - epsilon); const double z_eval = - std::clamp(std::abs(abs_z_cm), - z_min + epsilon, - z_max - epsilon); + std::clamp(std::abs(abs_z_cm), + z_min + epsilon, + z_max - epsilon); // Input histogram is in V/m. Garfield expects V/cm. return 0.01 * hist->Interpolate(r_eval, z_eval); @@ -1930,7 +1926,7 @@ TPolyLine3D* PHGarfield::ReverseDriftGlobalCoords(double x_cm, double y_cm, doub } PHGarfield::ReverseDriftStatus PHGarfield::StopHere(const double x, const double y, const double z, - const double zPrevious) + const double zPrevious) { const double r = std::hypot(x, y); diff --git a/offline/packages/PHGarfield/PHGarfield.h b/offline/packages/PHGarfield/PHGarfield.h index a86191d5e7..79db09a4ad 100644 --- a/offline/packages/PHGarfield/PHGarfield.h +++ b/offline/packages/PHGarfield/PHGarfield.h @@ -78,12 +78,12 @@ class PHGarfield : public SubsysReco // The user is encouraged to add more routine to fit their analysis goals... // Existing macros should call this one. Input and returned polyline are in // local TPC/Garfield coordinates. - TPolyLine3D *ReverseDrift(double x_cm, double y_cm, double z_cm, double step_ns = 50.0, ReverseDriftStatus* status = nullptr); + TPolyLine3D *ReverseDrift(double x_cm, double y_cm, double z_cm, double step_ns = 50.0, ReverseDriftStatus *status = nullptr); // Debug/visualization helper. Input and returned polyline are in global // detector coordinates. Internally the drift is still computed in local TPC // coordinates to keep the Garfield gas tables valid. - TPolyLine3D *ReverseDriftGlobalCoords(double x_cm, double y_cm, double z_cm, double step_ns = 50.0, ReverseDriftStatus* status = nullptr); + TPolyLine3D *ReverseDriftGlobalCoords(double x_cm, double y_cm, double z_cm, double step_ns = 50.0, ReverseDriftStatus *status = nullptr); double GetRadius(size_t index) const { return radii.at(index); } @@ -151,28 +151,32 @@ class PHGarfield : public SubsysReco m_spaceChargeScaleOverride = {{true, true}}; } - void SetSpaceChargeScaleSide0(double value) { + void SetSpaceChargeScaleSide0(double value) + { m_spaceChargeScale_side0 = value; - m_spaceChargeScaleOverride[0] = true; + m_spaceChargeScaleOverride[0] = true; } - void SetSpaceChargeScaleSide1(double value) { - m_spaceChargeScale_side1 = value; - m_spaceChargeScaleOverride[1] = true; + void SetSpaceChargeScaleSide1(double value) + { + m_spaceChargeScale_side1 = value; + m_spaceChargeScaleOverride[1] = true; } - + void SetFrameChargeScale(double value) { m_frameChargeScale_side0 = value; m_frameChargeScale_side1 = value; m_frameChargeScaleOverride = {{true, true}}; } - void SetFrameChargeScaleSide0(double value) { + void SetFrameChargeScaleSide0(double value) + { m_frameChargeScale_side0 = value; m_frameChargeScaleOverride[0] = true; } - void SetFrameChargeScaleSide1(double value) { - m_frameChargeScale_side1 = value; - m_frameChargeScaleOverride[1] = true; + void SetFrameChargeScaleSide1(double value) + { + m_frameChargeScale_side1 = value; + m_frameChargeScaleOverride[1] = true; } double GetSpaceChargeScaleSide0() const { return m_spaceChargeScale_side0; } @@ -185,8 +189,9 @@ class PHGarfield : public SubsysReco // Offsets are endpoint perturbations relative to the nominal resistor-chain // boundary voltage, in volts. IFC and OFC are independent, giving four // tunable parameters: IFC South/North and OFC South/North. - void SetUseIFCVoltageDistortion(bool value) { - m_useIFCVoltageDistortion = value; + void SetUseIFCVoltageDistortion(bool value) + { + m_useIFCVoltageDistortion = value; m_fieldCageConfigurationOverride = true; } void SetIFCVoltageOffset(double side0_south_v, double side1_north_v) @@ -195,20 +200,40 @@ class PHGarfield : public SubsysReco m_ifcVoltageOffset_side1 = side1_north_v; m_fieldCageConfigurationOverride = true; } - void SetIFCVoltageOffsetSide0(double value_v) { m_ifcVoltageOffset_side0 = value_v; m_fieldCageConfigurationOverride = true; } - void SetIFCVoltageOffsetSide1(double value_v) { m_ifcVoltageOffset_side1 = value_v; m_fieldCageConfigurationOverride = true; } + void SetIFCVoltageOffsetSide0(double value_v) + { + m_ifcVoltageOffset_side0 = value_v; + m_fieldCageConfigurationOverride = true; + } + void SetIFCVoltageOffsetSide1(double value_v) + { + m_ifcVoltageOffset_side1 = value_v; + m_fieldCageConfigurationOverride = true; + } double GetIFCVoltageOffsetSide0() const { return m_ifcVoltageOffset_side0; } double GetIFCVoltageOffsetSide1() const { return m_ifcVoltageOffset_side1; } - void SetUseOFCVoltageDistortion(bool value) { m_useOFCVoltageDistortion = value; m_fieldCageConfigurationOverride = true; } + void SetUseOFCVoltageDistortion(bool value) + { + m_useOFCVoltageDistortion = value; + m_fieldCageConfigurationOverride = true; + } void SetOFCVoltageOffset(double side0_south_v, double side1_north_v) { m_ofcVoltageOffset_side0 = side0_south_v; m_ofcVoltageOffset_side1 = side1_north_v; m_fieldCageConfigurationOverride = true; } - void SetOFCVoltageOffsetSide0(double value_v) { m_ofcVoltageOffset_side0 = value_v; m_fieldCageConfigurationOverride = true; } - void SetOFCVoltageOffsetSide1(double value_v) { m_ofcVoltageOffset_side1 = value_v; m_fieldCageConfigurationOverride = true; } + void SetOFCVoltageOffsetSide0(double value_v) + { + m_ofcVoltageOffset_side0 = value_v; + m_fieldCageConfigurationOverride = true; + } + void SetOFCVoltageOffsetSide1(double value_v) + { + m_ofcVoltageOffset_side1 = value_v; + m_fieldCageConfigurationOverride = true; + } double GetOFCVoltageOffsetSide0() const { return m_ofcVoltageOffset_side0; } double GetOFCVoltageOffsetSide1() const { return m_ofcVoltageOffset_side1; } @@ -281,8 +306,8 @@ class PHGarfield : public SubsysReco bool LoadCDBInputs(PHCompositeNode *topNode); void ConfigureRunDependentFieldCage(); - CDBTTree *m_cdbTPCMAPttree{nullptr}; // Locations of the pads from CDB... - //PHField3DCartesian *m_field{nullptr}; // The standard sPHENIX field holding container. + CDBTTree *m_cdbTPCMAPttree{nullptr}; // Locations of the pads from CDB... + // PHField3DCartesian *m_field{nullptr}; // The standard sPHENIX field holding container. PHField *m_field{nullptr}; Garfield::ComponentUser *m_component{nullptr}; // This handles the interface of the electric and magnetic fields as handed to Garfield Garfield::MediumMagboltz *m_gas{nullptr}; // This is the pre-tabulated gas properties required by Garfield... @@ -307,8 +332,8 @@ class PHGarfield : public SubsysReco double m_spaceChargeScale_side1{1.0}; // north, z > 0 double m_CMVoltageDefault{375.0}; // V/cm, nominal TPC field bool m_zerofield{false}; - TH2 *m_erCorrection{nullptr}; // radial correction, input bins in V/m - TH2 *m_ezCorrection{nullptr}; // local longitudinal correction, input bins in V/m + TH2 *m_erCorrection{nullptr}; // radial correction, input bins in V/m + TH2 *m_ezCorrection{nullptr}; // local longitudinal correction, input bins in V/m // Component order is Ex, Ey, Ez. Ez is along +|z| in the map. std::array, 2> m_field3DCorrection{}; diff --git a/offline/packages/TpcConditions/TpcConditions.h b/offline/packages/TpcConditions/TpcConditions.h index f481b62b0c..7e8fedf505 100644 --- a/offline/packages/TpcConditions/TpcConditions.h +++ b/offline/packages/TpcConditions/TpcConditions.h @@ -67,7 +67,7 @@ class TpcConditions float get_AverageLoadNR1() const { return m_AverageLoadNR1; } float get_AverageLoadNR2() const { return m_AverageLoadNR2; } float get_AverageLoadNR3() const { return m_AverageLoadNR3; } - + float get_Temperature() const { return m_Temperature; } float get_Pressure() const { return m_Pressure; } bool get_FieldOK() const { return m_FieldOK; } @@ -86,7 +86,7 @@ class TpcConditions float m_LoadNR1{0.0}; float m_LoadNR2{0.0}; float m_LoadNR3{0.0}; - + float m_AverageLoadCurrent{0.0}; float m_AverageLoadNorth{0.0}; float m_AverageLoadSouth{0.0}; diff --git a/offline/packages/TpcConditions/TpcConditionsReco.cc b/offline/packages/TpcConditions/TpcConditionsReco.cc index 28cc77f06c..bdc30d7137 100644 --- a/offline/packages/TpcConditions/TpcConditionsReco.cc +++ b/offline/packages/TpcConditions/TpcConditionsReco.cc @@ -408,7 +408,6 @@ int TpcConditionsReco::InitRun(PHCompositeNode *topNode) int TpcConditionsReco::process_event(PHCompositeNode *topNode) { - if (!m_tree || m_bco_to_channel.empty()) { return Fun4AllReturnCodes::EVENT_OK; @@ -437,7 +436,7 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) { --iter; } - + int channel = iter->second; m_conditions->set_Temperature(m_tree->GetFloatValue(channel, "gas_temperature")); @@ -445,7 +444,6 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) m_conditions->set_FieldOK(m_tree->GetFloatValue(channel, "FieldOK") != 0.0); m_conditions->set_GainOK(m_tree->GetFloatValue(channel, "GainOK") != 0.0); - m_conditions->set_LoadCurrent(get_InterpolatedMedianCurrent(bco, ALL)); m_conditions->set_LoadNorth(get_InterpolatedMedianCurrent(bco, NORTH)); m_conditions->set_LoadSouth(get_InterpolatedMedianCurrent(bco, SOUTH)); @@ -461,9 +459,9 @@ int TpcConditionsReco::process_event(PHCompositeNode *topNode) m_conditions->set_ConditionsAvailable(true); // Do or die - if( m_conditions->get_ConditionsAvailable() && - (!m_conditions->get_FieldOK() || - !m_conditions->get_GainOK())) + if (m_conditions->get_ConditionsAvailable() && + (!m_conditions->get_FieldOK() || + !m_conditions->get_GainOK())) { return Fun4AllReturnCodes::ABORTEVENT; } diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc index bbba2972c9..24564fd84a 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.cc @@ -1295,7 +1295,7 @@ int Tpc_PolyClusterizer::process_event(PHCompositeNode* topNode) const int iphi = static_cast(p.pad); const int it = static_cast(p.tbin); - if(it >= layergeom->get_zbins()) + if (it >= layergeom->get_zbins()) { continue; } diff --git a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h index 9313358b8f..31ea146ae2 100644 --- a/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h +++ b/offline/packages/tpctrackreco/Tpc_PolyClusterizer.h @@ -242,7 +242,7 @@ class Tpc_PolyClusterizer : public SubsysReco unsigned char m_maxAcceptedTier{1}; Tpc_AssembledTrackContainer* m_assembledTracks{nullptr}; Tpc_PolyClusterContainer* m_clusters{nullptr}; - TpcCrossingDecisionContainer* m_crossingDecisions {nullptr}; + TpcCrossingDecisionContainer* m_crossingDecisions{nullptr}; TrkrHitSetContainer* m_hits{nullptr}; IdealPadMap* m_idealPadMap{nullptr}; @@ -254,10 +254,9 @@ class Tpc_PolyClusterizer : public SubsysReco unsigned int m_event{0}; double m_t0{8}; double m_tpcAdcClock{56.881262}; - double m_crossingPeriodNs {106.56}; + double m_crossingPeriodNs{106.56}; double m_reverseDriftStepNs{56.881262}; - //! starting z position for primary electron backward drift /** * quoted values must be kept consistent with _max_driftlength + _CM_halfwidth