From 547875f1faa4f4abce58dfb1b1f2be8fececa62a Mon Sep 17 00:00:00 2001
From: Martin Wessels <martin.wessels@cern.ch>
Date: Wed, 30 Oct 2024 14:52:00 +0100
Subject: [PATCH] Improved calculation of LUT offset

---
 .../src/L1TriggerTowerToolRun3.cxx            | 67 +++++++++----------
 .../src/L1TriggerTowerToolRun3.h              |  2 +
 2 files changed, 35 insertions(+), 34 deletions(-)

diff --git a/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.cxx b/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.cxx
index cc84d6d1edb4..93e3b3d87ec0 100644
--- a/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.cxx
+++ b/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.cxx
@@ -221,7 +221,6 @@ void L1TriggerTowerToolRun3::simulateChannel(const xAOD::TriggerTower& tt, std::
 
   //If we have 80 MHz readout, we need to extract the 40 MHz samples. The central 80 MHz sample is always a 40 MHz sample. We use the cool database (runParameters folder) to understand if we are in 80MHz readout
   
-  
   SG::ReadCondHandle<L1CaloRunParametersContainer> runParameters( m_runParametersContainer);
   unsigned int readoutConfigID   = runParameters->runParameters(1)->readoutConfigID(); 
   ATH_MSG_DEBUG("RunParameters:: readoutConfigID " <<  readoutConfigID);
@@ -689,7 +688,6 @@ void L1TriggerTowerToolRun3::cpLut(const std::vector<int> &fir, const L1CaloCool
   int startBit = 0;
   int strategy = 0;
   int offset   = 0;
-  double offsetReal = 0;
   int slope    = 0;
   int cut      = 0;
   unsigned short scale_menu = 0;
@@ -717,17 +715,11 @@ void L1TriggerTowerToolRun3::cpLut(const std::vector<int> &fir, const L1CaloCool
       auto l1Menu = getL1Menu(ctx);
       scale_menu = l1Menu->thrExtraInfo().EM().emScale(); // Retrieve scale param from menu
     
-      
       for( auto &coeffs :  *hwCoeffs) { 
 	hwCoeffSum += coeffs;
       }
-      if (strategy == 0){
-        offsetReal = pedMean * hwCoeffSum / pow(2.,startBit);
-      }
-      else{
-	offsetReal = pedMean * static_cast<double>(hwCoeffSum) * static_cast<double>(slope) / pow(2.,static_cast<double>(startBit)) - static_cast<double>(slope)/2.;
-      }
-      offset = static_cast<unsigned short>( offsetReal < 0. ? 0 : offsetReal + 0.5 );
+      
+      offset = this->getLutOffset(pedMean, startBit, *hwCoeffs, slope, strategy);
 
       ATH_MSG_DEBUG( "::cpLut: Offset: offset/strategy/pedMean/firCoeffSum/startBit/slope: "
 		       << offset << " " << strategy << " " << " " << pedMean << " " << hwCoeffSum << " " << startBit << " " << slope );
@@ -764,7 +756,6 @@ void L1TriggerTowerToolRun3::jepLut(const std::vector<int> &fir, const L1CaloCoo
   int startBit = 0;
   int strategy   = 0;
   int offset     = 0;
-  double offsetReal = 0;
   int slope      = 0;
   int cut        = 0;
   unsigned short scale_db   = 0;
@@ -813,13 +804,8 @@ void L1TriggerTowerToolRun3::jepLut(const std::vector<int> &fir, const L1CaloCoo
       for( auto &coeffs :  *hwCoeffs) { 
 	 hwCoeffSum += coeffs;
       }
-      if (strategy == 0){
-        offsetReal = pedMean * hwCoeffSum / pow(2.,startBit);
-      }
-      else{
-        offsetReal = pedMean * static_cast<double>(hwCoeffSum) * static_cast<double>(slope) / pow(2.,static_cast<double>(startBit)) - static_cast<double>(slope)/2.;
-      }
-      offset = static_cast<unsigned short>( offsetReal < 0. ? 0 : offsetReal + 0.5 );
+
+      offset = this->getLutOffset(pedMean, startBit, *hwCoeffs, slope, strategy);
 
       ATH_MSG_VERBOSE( "::jepLut: Offset: offset/strategy/pedMean/firCoeffSum/startBit/slope: "
 		       << offset << " " << strategy << " " << " " << pedMean << " " << hwCoeffSum << " " << startBit << " " << slope );
@@ -1118,7 +1104,6 @@ void L1TriggerTowerToolRun3::cpLutParams(const L1CaloCoolChannelId& channelId, i
   startBit = 0;
   strategy = 0;
   offset   = 0;
-  double offsetReal = 0;
   slope    = 0;
   cut      = 0;
   pedValue = 0;
@@ -1126,6 +1111,7 @@ void L1TriggerTowerToolRun3::cpLutParams(const L1CaloCoolChannelId& channelId, i
   disabled = true;
   int hwCoeffSum = 0;
   const std::vector<short int>* hwCoeffs;
+
   
   if(!isRun2()) {
     // assert instead ?!
@@ -1149,13 +1135,8 @@ void L1TriggerTowerToolRun3::cpLutParams(const L1CaloCoolChannelId& channelId, i
       for( auto &coeffs :  *hwCoeffs) { 
 	hwCoeffSum += coeffs;
       }
-      if (strategy == 0){
-	offsetReal = pedMean * hwCoeffSum / pow(2.,startBit);
-      }
-      else{
-	offsetReal = pedMean * hwCoeffSum * slope / pow(2.,startBit) - slope/2.;
-      }
-      offset = static_cast<unsigned short>( offsetReal < 0. ? 0 : offsetReal + 0.5 );
+
+      offset  = this->getLutOffset(pedMean, startBit, *hwCoeffs, slope, strategy);
       
       ATH_MSG_VERBOSE( "::jepLutParams: Offset: offset/strategy/pedMean/firCoeffSum/startBit/slope: "
 		     << offset << " " << strategy << " " << " " << pedMean << " " << hwCoeffSum << " " << startBit << " " << slope );
@@ -1175,7 +1156,6 @@ void L1TriggerTowerToolRun3::jepLutParams(const L1CaloCoolChannelId& channelId,
   startBit = 0;
   strategy = 0;
   offset   = 0;
-  double offsetReal = 0;
   slope    = 0;
   cut      = 0;
   pedValue = 0;
@@ -1205,13 +1185,8 @@ void L1TriggerTowerToolRun3::jepLutParams(const L1CaloCoolChannelId& channelId,
       for( auto &coeffs :  *hwCoeffs) { 
 	hwCoeffSum += coeffs;
       }
-      if (strategy == 0){
-	offsetReal = pedMean * hwCoeffSum / pow(2.,startBit);
-      }
-      else{
-	offsetReal = pedMean * hwCoeffSum * slope / pow(2.,startBit) - slope/2.;
-      }
-      offset = static_cast<unsigned short>( offsetReal < 0. ? 0 : offsetReal + 0.5 );
+
+      offset = this->getLutOffset(pedMean, startBit, *hwCoeffs, slope, strategy);
       
       ATH_MSG_VERBOSE( "::jepLutParams: Offset: offset/strategy/pedMean/firCoeffSum/startBit/slope: "
 		       << offset << " " << strategy << " " << " " << pedMean << " " << hwCoeffSum << " " << startBit << " " << slope );
@@ -1510,7 +1485,31 @@ bool L1TriggerTowerToolRun3::isRun2() const
   return false;
 }
 
+unsigned int L1TriggerTowerToolRun3::getLutOffset(const double &pedMean, const unsigned int &firStartBit, const std::vector<short int> &firCoeff, const unsigned int &lutSlope, const unsigned int &lutStrategy) const
+{
+  unsigned int lutOffset = 0;
+  // essential to save in long long to avoid rounding errors
+  long long int lutOffsetLong = 0;
+  long long int lutSlopeLong = lutSlope;
+  long long int firStartBitLong = firStartBit;
+  long long int pedMeanLong = std::lround(pedMean * 10000.);
+  long long int firCoeffSum = 0;
+  
+  for (unsigned int i=0; i<firCoeff.size(); i++) {
+    firCoeffSum += firCoeff.at(i);
+  }
+  
+  if ( lutStrategy == 0 ) {
+    lutOffsetLong = ((pedMeanLong*firCoeffSum) >> firStartBitLong);
+  }
+  else {
+    lutOffsetLong = ((pedMeanLong*firCoeffSum*lutSlopeLong) >> firStartBitLong) - ((lutSlopeLong * 10000) >> 1);    
+  }
 
+  lutOffsetLong = (lutOffsetLong + (10000-1))/10000;
+  lutOffset = static_cast<unsigned int>( lutOffsetLong < 0 ? 0 : lutOffsetLong );
+  return lutOffset;
+}
 
 const TrigConf::L1Menu* L1TriggerTowerToolRun3::getL1Menu(const EventContext& ctx) const {
   const TrigConf::L1Menu* menu = nullptr;
diff --git a/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.h b/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.h
index 86dc381f1c0a..9bdd1c1412a2 100644
--- a/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.h
+++ b/Trigger/TrigT1/TrigT1CaloTools/src/L1TriggerTowerToolRun3.h
@@ -126,6 +126,8 @@ namespace LVL1
       /** Get extra noise cut with disabled channel */
       bool disabledChannel(const L1CaloCoolChannelId& channelId, unsigned int& noiseCut) const;
 
+      // calculate the LUT offset from DB parameters
+      unsigned int getLutOffset(const double &pedMean, const unsigned int &firStartBit, const std::vector<short int> &firCoeff, const unsigned int &lutSlope, const unsigned int &lutStrategy) const;
 
       /// Id managers
       const CaloIdManager* m_caloMgr;
-- 
GitLab