Update FEMC for new geometry:

src/algorithms/calorimetry/CalorimeterHitDigi.cc

@@ -209,9 +209,17 @@ void CalorimeterHitDigi::process(const CalorimeterHitDigi::Input& input,
 
     double ped = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC;
 
+    //SiPM Saturation
+    double saturation = 1.0;
+    if (m_cfg.totalPixel > 0) {
+      double nPhoton = edep * corrMeanScale_value * (1.0 + eResRel) * m_cfg.nPhotonPerGeV;
+      saturation     = (1.0 - exp(-nPhoton / m_cfg.totalPixel)) * m_cfg.totalPixel / nPhoton;
+      //printf("edep=%8.4f nPhoton=%8.2f totalPixel=%8.1f Saturation=%8.6f\n",edep,nPhoton,m_cfg.totalPixel,saturation);
+    }
+
     // Note: both adc and tdc values must be positive numbers to avoid integer wraparound
     unsigned long long adc =
-        std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) /
+        std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) * saturation /
                                         m_cfg.dyRangeADC * m_cfg.capADC),
                  0LL);
     unsigned long long tdc = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC);

src/algorithms/calorimetry/CalorimeterHitDigiConfig.h

@@ -28,6 +28,9 @@ struct CalorimeterHitDigiConfig {
   // signal sums
   std::string readout{""};
   std::vector<std::string> fields{};
-};
 
+  //SiPM Saturation
+  double totalPixel{0};
+  double nPhotonPerGeV{0};
+};
 } // namespace eicrecon

src/detectors/FEMC/FEMC.cc

@@ -1,3 +1,4 @@
+
 // SPDX-License-Identifier: LGPL-3.0-or-later
 // Copyright (C) 2021 - 2025, Chao Peng, Sylvester Joosten, Whitney Armstrong, David Lawrence, Friederike Bock, Wouter Deconinck, Kolja Kauder, Sebouh Paul
 
@@ -8,6 +9,7 @@
 #include <string>
 
 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
+#include "algorithms/calorimetry/CalorimeterHitRecoConfig.h"
 #include "extensions/jana/JOmniFactoryGeneratorT.h"
 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
@@ -23,38 +25,107 @@ void InitPlugin(JApplication* app) {
   using namespace eicrecon;
 
   InitJANAPlugin(app);
+
+  auto log_service = app->GetService<Log_service>();
+  auto mLog        = log_service->logger("FEMC");
+
+  int SiPMSaturation = 1;
+  std::string opt("ON");
+  app->SetDefaultParameter("FEMC:SiPMSaturation", opt, "Turn ON(default) or OFF SiPM Saturation");
+  if (opt.find("OFF") != std::string::npos || opt.find("off") != std::string::npos ||
+      opt.find("Off") != std::string::npos) {
+    mLog->info("SiPM Saturation OFF");
+    SiPMSaturation = 0;
+  } else {
+    mLog->info("SiPM Saturation ON");
+  }
+
   // Make sure digi and reco use the same value
   decltype(CalorimeterHitDigiConfig::capADC) EcalEndcapP_capADC =
       16384; //16384, assuming 14 bits. For approximate HGCROC resolution use 65536
-  decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalEndcapP_dyRangeADC   = 3 * dd4hep::GeV;
+  decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalEndcapP_dyRangeADC   = 100 * dd4hep::GeV;
   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalEndcapP_pedMeanADC   = 200;
   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalEndcapP_pedSigmaADC = 2.4576;
   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalEndcapP_resolutionTDC =
       10 * dd4hep::picosecond;
-  app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
-      "EcalEndcapPRawHits", {"EcalEndcapPHits"},
+  //const double sampFrac =  0.03;
+  const double sampFrac = 0.029043; //updated with ratio to ScFi model
+  decltype(CalorimeterHitDigiConfig::corrMeanScale) EcalEndcapP_corrMeanScale =
+      Form("%f", 1.0 / sampFrac);
+  decltype(CalorimeterHitRecoConfig::sampFrac) EcalEndcapP_sampFrac = Form("%f", sampFrac);
+  decltype(CalorimeterHitDigiConfig::nPhotonPerGeV) EcalEndcapP_nPhotonPerGeV = 1500;
+  decltype(CalorimeterHitDigiConfig::totalPixel) EcalEndcapP_totalPixel       = 4 * 159565;
+  if (SiPMSaturation == 0)
+    EcalEndcapP_totalPixel = 0;
+
+  int FEMCHomoScfi = 0;
+  try {
+    auto detector = app->GetService<DD4hep_service>()->detector();
+    FEMCHomoScfi  = detector->constant<int>("ForwardEcal_Homogeneous_Scfi");
+    if (FEMCHomoScfi <= 1)
+      mLog->info("Homogeneous geometry loaded");
+    else
+      mLog->info("ScFi geometry loaded");
+  } catch (...) {
+  };
+
+  if (FEMCHomoScfi <= 1) {
+    app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
+        "EcalEndcapPRawHits", {"EcalEndcapPHits"},
 #if EDM4EIC_VERSION_MAJOR >= 7
-      {"EcalEndcapPRawHits", "EcalEndcapPRawHitAssociations"},
+        {"EcalEndcapPRawHits", "EcalEndcapPRawHitAssociations"},
 #else
-      {"EcalEndcapPRawHits"},
+        {"EcalEndcapPRawHits"},
 #endif
-      {
-          .eRes      = {0.11333 * sqrt(dd4hep::GeV), 0.03,
-                        0.0 * dd4hep::GeV}, // (11.333% / sqrt(E)) \oplus 3%
-          .tRes      = 0.0,
-          .threshold = 0.0,
-          // .threshold = 15 * dd4hep::MeV for a single tower, applied on ADC level
-          .capADC        = EcalEndcapP_capADC,
-          .capTime       = 100, // given in ns, 4 samples in HGCROC
-          .dyRangeADC    = EcalEndcapP_dyRangeADC,
-          .pedMeanADC    = EcalEndcapP_pedMeanADC,
-          .pedSigmaADC   = EcalEndcapP_pedSigmaADC,
-          .resolutionTDC = EcalEndcapP_resolutionTDC,
-          .corrMeanScale = "0.03",
-          .readout       = "EcalEndcapPHits",
-      },
-      app // TODO: Remove me once fixed
-      ));
+        {
+            .eRes = {0.11333 * sqrt(dd4hep::GeV), 0.03,
+                     0.0 * dd4hep::GeV}, // (11.333% / sqrt(E)) \oplus 3%
+            .tRes = 0.0,
+            .threshold =
+                0.0, // 15MeV threshold for a single tower will be applied on ADC at Reco below
+            .capADC        = EcalEndcapP_capADC,
+            .capTime       = 100, // given in ns, 4 samples in HGCROC
+            .dyRangeADC    = EcalEndcapP_dyRangeADC,
+            .pedMeanADC    = EcalEndcapP_pedMeanADC,
+            .pedSigmaADC   = EcalEndcapP_pedSigmaADC,
+            .resolutionTDC = EcalEndcapP_resolutionTDC,
+            .corrMeanScale = "1.0",
+            .readout       = "EcalEndcapPHits",
+            .totalPixel    = EcalEndcapP_totalPixel,
+            .nPhotonPerGeV = EcalEndcapP_nPhotonPerGeV,
+        },
+        app // TODO: Remove me once fixed
+        ));
+  } else if (FEMCHomoScfi == 2) {
+    app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
+        "EcalEndcapPRawHits", {"EcalEndcapPHits"},
+#if EDM4EIC_VERSION_MAJOR >= 7
+        {"EcalEndcapPRawHits", "EcalEndcapPRawHitAssociations"},
+#else
+        {"EcalEndcapPRawHits"},
+#endif
+        {
+            //.eRes = {0.0, 0.0, 0.0}, // No smearing for ScFi
+            .eRes = {0.0, 0.022, 0.0}, // just constant term 2.2% based on MC data comparison
+            .tRes = 0.0,
+            .threshold =
+                0.0, // 15MeV threshold for a single tower will be applied on ADC at Reco below
+            .capADC        = EcalEndcapP_capADC,
+            .capTime       = 100, // given in ns, 4 samples in HGCROC
+            .dyRangeADC    = EcalEndcapP_dyRangeADC,
+            .pedMeanADC    = EcalEndcapP_pedMeanADC,
+            .pedSigmaADC   = EcalEndcapP_pedSigmaADC,
+            .resolutionTDC = EcalEndcapP_resolutionTDC,
+            .corrMeanScale = EcalEndcapP_corrMeanScale,
+            .readout       = "EcalEndcapPHits",
+            .fields        = {"fiberx", "fibery", "x", "y"},
+            .totalPixel    = EcalEndcapP_totalPixel,
+            .nPhotonPerGeV = EcalEndcapP_nPhotonPerGeV,
+        },
+        app // TODO: Remove me once fixed
+        ));
+  }
+
   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
       "EcalEndcapPRecHits", {"EcalEndcapPRawHits"}, {"EcalEndcapPRecHits"},
       {
@@ -65,8 +136,9 @@ void InitPlugin(JApplication* app) {
           .resolutionTDC   = EcalEndcapP_resolutionTDC,
           .thresholdFactor = 0.0,
           .thresholdValue =
-              2, // The ADC of a 15 MeV particle is adc = 200 + 15 * 0.03 * ( 1.0 + 0) / 3000 * 16384 = 200 + 2.4576
-          .sampFrac = "0.03",
+              //    2, // The ADC of a 15 MeV particle is adc = 200 + 15 * 0.03 * ( 1.0 + 0) / 3000 * 16384 = 200 + 2.4576
+          3, // 15 MeV = 2.4576, but adc=llround(dE) and cut off is "<". So 3 here = 15.25MeV
+          .sampFrac = "1.00", // already taken care in DIGI code above
           .readout  = "EcalEndcapPHits",
       },
       app // TODO: Remove me once fixed

src/factories/calorimetry/CalorimeterHitDigi_factory.h

@@ -36,6 +36,8 @@ class CalorimeterHitDigi_factory
   ParameterRef<std::string> m_corrMeanScale{this, "scaleResponse", config().corrMeanScale};
   ParameterRef<std::vector<std::string>> m_fields{this, "signalSumFields", config().fields};
   ParameterRef<std::string> m_readout{this, "readoutClass", config().readout};
+  ParameterRef<double> m_nPhotonPerGeV{this, "numberOfPhotonPerGeV", config().nPhotonPerGeV};
+  ParameterRef<double> m_totalPixel{this, "totalNnumberOfPixel", config().totalPixel};
 
   Service<AlgorithmsInit_service> m_algorithmsInit{this};