From 8d57d57fc98b59fdeed2e75007acb1b022817e51 Mon Sep 17 00:00:00 2001
From: Luca Martinelli <luca.martinelli@cern.ch>
Date: Tue, 6 Dec 2022 17:02:52 +0100
Subject: [PATCH] New MM thresholds

New MM thresholds
---
 .../MM_Digitization/MM_DigitizationTool.h     |  13 +-
 .../src/MM_DigitizationTool.cxx               | 121 ++++++++++++++----
 Tools/PROCTools/data/q445_AOD_digest.ref      |   2 +-
 Tools/WorkflowTestRunner/python/References.py |   4 +-
 4 files changed, 105 insertions(+), 35 deletions(-)

diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitizationTool.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitizationTool.h
index 63fbd9d387da..76307e37e490 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitizationTool.h
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_DigitizationTool.h
@@ -64,6 +64,7 @@
 #include "StoreGate/ReadCondHandleKey.h"
 #include "xAODEventInfo/EventAuxInfo.h"  // SubEventIterator
 #include "xAODEventInfo/EventInfo.h"     // SubEventIterator
+#include "MuonCondData/NswCalibDbTimeChargeData.h"
 
 /*******************************************************************************/
 
@@ -143,7 +144,7 @@ private:
     Gaudi::Property<bool> m_needsMcEventCollHelper{this, "UseMcEventCollectionHelper", false};
     Gaudi::Property<bool> m_checkMMSimHits{this, "CheckSimHits", true, "Control on the hit validity"};
     Gaudi::Property<bool> m_useTimeWindow{this, "UseTimeWindow", true};
-    Gaudi::Property<bool> m_vmmNeighborLogic{this, "VMMNeighborLogic", false};
+    Gaudi::Property<bool> m_vmmNeighborLogic{this, "VMMNeighborLogic", true};
     Gaudi::Property<bool> m_doSmearing{this, "doSmearing", false,
                                        "set the usage or not of the smearing tool for realistic detector performance"};
     // Constants vars for the MM_StripsResponseSimulation class
@@ -157,7 +158,7 @@ private:
     Gaudi::Property<float> m_crossTalk1{this, "crossTalk1", 0.3, "Strip Cross Talk with Nearest Neighbor"};
     Gaudi::Property<float> m_crossTalk2{this, "crossTalk2", 0.09, "Strip Cross Talk with 2nd Nearest Neighbor"};
 
-    Gaudi::Property<float> m_avalancheGain{this, "AvalancheGain", 8.0e3, "avalanche Gain for rach gas mixture"};
+    Gaudi::Property<float> m_avalancheGain{this, "AvalancheGain", 6.0e3, "avalanche Gain for rach gas mixture"};
 
     // Constants vars for the MM_ElectronicsResponseSimulation
     Gaudi::Property<float> m_electronicsThreshold{this, "electronicsThreshold", 15000,
@@ -176,7 +177,7 @@ private:
                                               "Use conditions data to get thresholds, overrules useThresholdScaling"};
     Gaudi::Property<bool> m_useThresholdScaling{this, "useThresholdScaling", true,
                                                 "Use a strip length dependent threshold in MM digitiation"};
-    Gaudi::Property<float> m_thresholdScaleFactor{this, "thresholdScaleFactor", 9.0,
+    Gaudi::Property<float> m_thresholdScaleFactor{this, "thresholdScaleFactor", 10.0,
                                                   "Use x times the strip length dependent noise as MM threshold"};
     Gaudi::Property<float> m_vmmDeadtime{
         this, "vmmDeadtime", 200, "Specifies how much before the lower time limit the VMM simulation should start evaluating the signal"};
@@ -195,8 +196,10 @@ private:
     std::unique_ptr<MM_StripsResponseSimulation> m_StripsResponseSimulation{};
     std::unique_ptr<MM_ElectronicsResponseSimulation> m_ElectronicsResponseSimulation{};
 
-    double m_noiseSlope{0.};
-    double m_noiseIntercept{0.};
+    using NoiseCalibConstants = NswCalibDbTimeChargeData::CalibConstants;
+    /// Define a map to cache the noise parameters individually
+    /// Key: stationName * std::abs(stationEta)
+    std::map<int, NoiseCalibConstants> m_noiseParams{};
 };
 
 #endif  // MM_DigitizationTool
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx
index 72560ea9305c..4e93a97f4642 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_DigitizationTool.cxx
@@ -61,11 +61,17 @@
 
 namespace {
     // thresholds for the shortest and longest strips
-    // values from https://indico.cern.ch/event/1002207/contributions/4240818
-    // slide 10
+    //values from https://indico.cern.ch/event/1131762/contributions/4749097/attachments/2431773/4164431/MMGcoord2022.04.26.pdf
 
-    constexpr float maxNoise = 2600;
-    constexpr float minNoise = 1000;
+    constexpr float maxNoiseSmall_eta1 = 2100;
+    constexpr float minNoiseSmall_eta1 = 1000;
+    constexpr float maxNoiseSmall_eta2 = 2600;
+    constexpr float minNoiseSmall_eta2 = 2100;
+
+    constexpr float maxNoiseLarge_eta1 = 2500;
+    constexpr float minNoiseLarge_eta1 = 1200;
+    constexpr float maxNoiseLarge_eta2 = 3000;
+    constexpr float minNoiseLarge_eta2 = 2500;
 
 }  // namespace
 
@@ -169,37 +175,96 @@ StatusCode MM_DigitizationTool::initialize() {
     if (m_doSmearing) ATH_MSG_INFO("Running in smeared mode!");
 
     // get shortest and longest strip length for threshold scaling
-    //
-    Identifier tmpId;  // temporary identifier to work with ReadoutElement
-
-    // identifier for first gas gap in a large MM sector, layer is eta layer, station Eta = 1 to get shortest strip
-    tmpId = m_idHelperSvc->mmIdHelper().channelID("MML", 1, 1, 1, 1, 1);
-
+    Identifier tmpId{0};  // temporary identifier to work with ReadoutElement
     const MuonGM::MuonDetectorManager* muonGeoMgr{nullptr};
     ATH_CHECK(detStore()->retrieve(muonGeoMgr));
+    int stripNumberShortestStrip{-1}, stripNumberLongestStrip{-1};
+    Identifier tmpIdShortestStrip{0},tmpIdLongestStrip{0};
+    float shortestStripLength{FLT_MAX}, longestStripLength{0};   
+
+    //============================
+    //SMALL SECTORS - ETA 1 CHAMBERS
+    // identifier for first gas gap in a small MM sector, layer is eta layer
+    tmpId = m_idHelperSvc->mmIdHelper().channelID("MMS", 1, 1, 1, 1, 1);
     const MuonGM::MMReadoutElement* detectorReadoutElement = muonGeoMgr->getMMReadoutElement(tmpId);
-    int stripNumberShortestStrip = (detectorReadoutElement->getDesign(tmpId))->nMissedBottomEta + 1;
-    Identifier tmpIdShortestStrip = m_idHelperSvc->mmIdHelper().channelID(
-        "MML", 1, 1, 1, 1, stripNumberShortestStrip);  // identifier for first gas gap in a large MM sector, layer is eta layer
-    float shortestStripLength = detectorReadoutElement->stripLength(tmpIdShortestStrip);
+    stripNumberShortestStrip = (detectorReadoutElement->getDesign(tmpId))->nMissedBottomEta + 1;
+    tmpIdShortestStrip = m_idHelperSvc->mmIdHelper().channelID("MMS", 1, 1, 1, 1, stripNumberShortestStrip);  // identifier for the shortest strip
+    shortestStripLength = detectorReadoutElement->stripLength(tmpIdShortestStrip);
+
+    stripNumberLongestStrip = (detectorReadoutElement->getDesign(tmpId))->totalStrips - (detectorReadoutElement->getDesign(tmpId))->nMissedTopEta;
+    tmpIdLongestStrip = m_idHelperSvc->mmIdHelper().channelID("MMS", 1, 1, 1, 1, stripNumberLongestStrip);  // identifier for the longest strip
+    longestStripLength = detectorReadoutElement->stripLength(tmpIdLongestStrip);
+
+    // now get the slope and intercept for the threshold scaling
+    // function is m_noiseSlope * stripLength + m_noiseIntercept
+    /// Small wedges first eta station
+    NoiseCalibConstants& noise_smallEta1 = m_noiseParams[m_idHelperSvc->stationName(tmpId)];
+    noise_smallEta1.slope = (maxNoiseSmall_eta1 - minNoiseSmall_eta1) / (longestStripLength - shortestStripLength);
+    noise_smallEta1.intercept =  minNoiseSmall_eta1 - noise_smallEta1.slope* shortestStripLength;    
+    //============================
+
+    //============================
+    //SMALL SECTORS - ETA 2 CHAMBERS
+    // identifier for first gas gap in a small MM sector, layer is eta layer
+    tmpId = m_idHelperSvc->mmIdHelper().channelID("MMS", 2, 1, 1, 1, 1);
+    detectorReadoutElement = muonGeoMgr->getMMReadoutElement(tmpId);
+    stripNumberShortestStrip = (detectorReadoutElement->getDesign(tmpId))->nMissedBottomEta + 1;
+    tmpIdShortestStrip = m_idHelperSvc->mmIdHelper().channelID("MMS", 2, 1, 1, 1, stripNumberShortestStrip);  // identifier for the shortest strip
+    shortestStripLength = detectorReadoutElement->stripLength(tmpIdShortestStrip);
+
+    stripNumberLongestStrip = (detectorReadoutElement->getDesign(tmpId))->totalStrips - (detectorReadoutElement->getDesign(tmpId))->nMissedTopEta;
+    tmpIdLongestStrip = m_idHelperSvc->mmIdHelper().channelID("MMS", 2, 1, 1, 1, stripNumberLongestStrip);  // identifier for the longest strip
+    longestStripLength = detectorReadoutElement->stripLength(tmpIdLongestStrip);
 
-    // identifier for first gas gap in a large MM sector, layer is eta layer, station Eta = 2 to get the longest strip
-    Identifier tmpId2 = m_idHelperSvc->mmIdHelper().channelID("MML", 2, 1, 1, 1, 1);
+    // now get the slope and intercept for the threshold scaling
+    // function is m_noiseSlope * stripLength + m_noiseIntercept
+    /// Small wedges second eta station
+    NoiseCalibConstants& noise_smallEta2 = m_noiseParams[m_idHelperSvc->stationName(tmpId)*2];
+    noise_smallEta2.slope = (maxNoiseSmall_eta2 - minNoiseSmall_eta2) / (longestStripLength - shortestStripLength);
+    noise_smallEta2.intercept =  minNoiseSmall_eta2 - noise_smallEta2.slope * shortestStripLength;
+    //============================
+
+    //============================
+    //LARGE SECTORS - ETA 1 CHAMBERS
+    // identifier for first gas gap in a small MM sector, layer is eta layer
+    tmpId = m_idHelperSvc->mmIdHelper().channelID("MML", 1, 1, 1, 1, 1);
+    detectorReadoutElement = muonGeoMgr->getMMReadoutElement(tmpId);
+    stripNumberShortestStrip = (detectorReadoutElement->getDesign(tmpId))->nMissedBottomEta + 1;
+    tmpIdShortestStrip = m_idHelperSvc->mmIdHelper().channelID("MML", 1, 1, 1, 1, stripNumberShortestStrip);  // identifier for the shortest strip
+    shortestStripLength = detectorReadoutElement->stripLength(tmpIdShortestStrip);
 
-    detectorReadoutElement = muonGeoMgr->getMMReadoutElement(tmpId2);
-    int stripNumberLongestStrip =
-        (detectorReadoutElement->getDesign(tmpId))->totalStrips - (detectorReadoutElement->getDesign(tmpId))->nMissedTopEta;
-    Identifier tmpIdLongestStrip = m_idHelperSvc->mmIdHelper().channelID(
-        "MML", 2, 1, 1, 1, stripNumberLongestStrip);  // identifier for first gas gap in a large MM sector, layer is eta layer
-    float longestStripLength = detectorReadoutElement->stripLength(tmpIdLongestStrip);
+    stripNumberLongestStrip = (detectorReadoutElement->getDesign(tmpId))->totalStrips - (detectorReadoutElement->getDesign(tmpId))->nMissedTopEta;
+    tmpIdLongestStrip = m_idHelperSvc->mmIdHelper().channelID("MML", 1, 1, 1, 1, stripNumberLongestStrip);  // identifier for the longest strip
+    longestStripLength = detectorReadoutElement->stripLength(tmpIdLongestStrip);
 
     // now get the slope and intercept for the threshold scaling
     // function is m_noiseSlope * stripLength + m_noiseIntercept
-    m_noiseSlope = (maxNoise - minNoise) / (longestStripLength - shortestStripLength);
-    m_noiseIntercept = minNoise - m_noiseSlope * shortestStripLength;
+    /// Large wedges first eta station
+    NoiseCalibConstants& noise_largeEta1 = m_noiseParams[m_idHelperSvc->stationName(tmpId)];
+    noise_largeEta1.slope = (maxNoiseLarge_eta1 - minNoiseLarge_eta1) / (longestStripLength - shortestStripLength);
+    noise_largeEta1.intercept =  minNoiseLarge_eta1 - noise_largeEta1.slope * shortestStripLength;    
+    //============================
+
+    //============================
+    //LARGE SECTORS - ETA 2 CHAMBERS
+    // identifier for first gas gap in a small MM sector, layer is eta layer
+    tmpId = m_idHelperSvc->mmIdHelper().channelID("MML", 2, 1, 1, 1, 1);
+    detectorReadoutElement = muonGeoMgr->getMMReadoutElement(tmpId);
+    stripNumberShortestStrip = (detectorReadoutElement->getDesign(tmpId))->nMissedBottomEta + 1;
+    tmpIdShortestStrip = m_idHelperSvc->mmIdHelper().channelID("MML", 2, 1, 1, 1, stripNumberShortestStrip);  // identifier for the shortest strip
+    shortestStripLength = detectorReadoutElement->stripLength(tmpIdShortestStrip);
+
+    stripNumberLongestStrip = (detectorReadoutElement->getDesign(tmpId))->totalStrips - (detectorReadoutElement->getDesign(tmpId))->nMissedTopEta;
+    tmpIdLongestStrip = m_idHelperSvc->mmIdHelper().channelID("MML", 2, 1, 1, 1, stripNumberLongestStrip);  // identifier for the longest strip
+    longestStripLength = detectorReadoutElement->stripLength(tmpIdLongestStrip);
 
-    ATH_MSG_DEBUG(" Shortest strip number: " << stripNumberShortestStrip << " length: " << shortestStripLength
-                                             << " longest strip number: " << stripNumberLongestStrip << " length " << longestStripLength);
+    // now get the slope and intercept for the threshold scaling
+    // function is m_noiseSlope * stripLength + m_noiseIntercept
+    /// Large wedges first eta station
+    NoiseCalibConstants& noise_largeEta2 = m_noiseParams[m_idHelperSvc->stationName(tmpId)*2];
+    noise_largeEta2.slope = (maxNoiseLarge_eta2 - minNoiseLarge_eta2) / (longestStripLength - shortestStripLength);
+    noise_largeEta2.intercept =  minNoiseLarge_eta1 - noise_largeEta2.slope * shortestStripLength;    
+    //============================
 
     ATH_MSG_DEBUG("Configuration  MM_DigitizationTool ");
     ATH_MSG_INFO("RndmSvc                " << m_rndmSvc);
@@ -963,7 +1028,9 @@ MM_ElectronicsToolInput MM_DigitizationTool::combinedStripResponseAllHits(const
                                                                           m_idHelperSvc->mmIdHelper().gasGap(digitID), strip_id);
                     const MuonGM::MMReadoutElement* detectorReadoutElement = muonGeoMgr->getMMReadoutElement(id);
                     float stripLength = detectorReadoutElement->stripLength(id);
-                    float threshold = (m_noiseSlope * stripLength + m_noiseIntercept) * m_thresholdScaleFactor;
+                    const int noise_id = m_idHelperSvc->stationName(digitID) * std::abs(m_idHelperSvc->stationEta(digitID));
+                    const NoiseCalibConstants& noise = m_noiseParams.at(noise_id);
+                    float threshold = (noise.slope * stripLength + noise.intercept) * m_thresholdScaleFactor;
                     v_stripThresholdResponseAllHits.push_back(threshold);
                 } else {
                     v_stripThresholdResponseAllHits.push_back(m_electronicsThreshold);
diff --git a/Tools/PROCTools/data/q445_AOD_digest.ref b/Tools/PROCTools/data/q445_AOD_digest.ref
index 82010641af84..f6319da67585 100644
--- a/Tools/PROCTools/data/q445_AOD_digest.ref
+++ b/Tools/PROCTools/data/q445_AOD_digest.ref
@@ -14,7 +14,7 @@
       410000    43250013          57          83          30           4           0           6           1           5           8           6           2
       410000    43250014          48          52          20           5           1           7           2           5           6           5           1
       410000    43250015          68          82          36           5           1           8           2           6           4           4           0
-      410000    43250016          66          55          30           4           3          11           1          10           7           6           1
+      410000    43250016          66          55          30           4           4          11           1          10           7           6           1
       410000    43250017          48          44          36           7           1           4           1           3           8           6           2
       410000    43250018          54          55          21           5           2           3           0           3           5           4           1
       410000    43250019          92          94          69           9           2          12           0          12          11           4           7
diff --git a/Tools/WorkflowTestRunner/python/References.py b/Tools/WorkflowTestRunner/python/References.py
index eb8fe30a55ab..30f00210e067 100644
--- a/Tools/WorkflowTestRunner/python/References.py
+++ b/Tools/WorkflowTestRunner/python/References.py
@@ -22,10 +22,10 @@ references_map = {
     # Overlay
     "d1590": "v12",
     "d1726": "v8",
-    "d1759": "v14",
+    "d1759": "v15",
     # Reco
     "q442": "v1",
     "q443": "v1",
-    "q445": "v1",
+    "q445": "v2",
     "q449": "v2",
 }
-- 
GitLab