Update FEMC for new geometry:

src/detectors/FEMC/FEMC.cc

@@ -1,10 +1,16 @@
 // 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
+// Copyright (C) 2021 - 2025, Chao Peng, Sylvester Joosten, Whitney Armstrong, David Lawrence, Friederike Bock, Wouter Deconinck, Kolja Kauder, Sebouh Paul, Akio Ogawa
 
+#include <DD4hep/Detector.h>
 #include <Evaluator/DD4hepUnits.h>
+#include <JANA/JApplication.h>
 #include <JANA/JApplicationFwd.h>
 #include <JANA/Utils/JTypeInfo.h>
+#include <fmt/core.h>
+#include <spdlog/logger.h>
 #include <cmath>
+#include <gsl/pointers>
+#include <memory>
 #include <string>
 #include <variant>
 #include <vector>
@@ -18,41 +24,99 @@
 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
 #include "factories/calorimetry/CalorimeterTruthClustering_factory.h"
 #include "factories/calorimetry/TrackClusterMergeSplitter_factory.h"
+#include "services/geometry/dd4hep/DD4hep_service.h"
+#include "services/log/Log_service.h"
 
 extern "C" {
 void InitPlugin(JApplication* app) {
 
   using namespace eicrecon;
 
   InitJANAPlugin(app);
+
+  auto log_service = app->GetService<Log_service>();
+  auto mLog        = log_service->logger("FEMC");
+
   // 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", {"EventHeader", "EcalEndcapPHits"},
-      {"EcalEndcapPRawHits", "EcalEndcapPRawHitAssociations"},
-      {
-          .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
-      ));
+  const double EcalEndcapP_sampFrac = 0.029043; // updated with ratio to ScFi model
+  decltype(CalorimeterHitDigiConfig::corrMeanScale) EcalEndcapP_corrMeanScale =
+      fmt::format("{}", 1.0 / EcalEndcapP_sampFrac);
+  const double EcalEndcapP_nPhotonPerGeV          = 1500;
+  const double EcalEndcapP_PhotonCollectionEff    = 0.5;
+  const unsigned long long EcalEndcapP_totalPixel = 4 * 159565ULL;
+
+  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 (...) {
+    FEMCHomoScFi = 0;
+  };
+
+  if (FEMCHomoScFi <= 1) {
+    app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
+        "EcalEndcapPRawHits", {"EcalEndcapPHits"},
+        {"EcalEndcapPRawHits", "EcalEndcapPRawHitAssociations"},
+        {
+            .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
+            .readoutType = "sipm",
+            .lightYield  = EcalEndcapP_nPhotonPerGeV / EcalEndcapP_PhotonCollectionEff,
+            .photonDetectionEfficiency = EcalEndcapP_PhotonCollectionEff,
+            .numEffectiveSipmPixels    = EcalEndcapP_totalPixel,
+            .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",
+        },
+        app // TODO: Remove me once fixed
+        ));
+  } else if (FEMCHomoScFi == 2) {
+    app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
+        "EcalEndcapPRawHits", {"EcalEndcapPHits"},
+        {"EcalEndcapPRawHits", "EcalEndcapPRawHitAssociations"},
+        {
+            .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
+            .readoutType = "sipm",
+            .lightYield  = EcalEndcapP_nPhotonPerGeV / EcalEndcapP_PhotonCollectionEff,
+            .photonDetectionEfficiency = EcalEndcapP_PhotonCollectionEff,
+            .numEffectiveSipmPixels    = EcalEndcapP_totalPixel,
+            .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"},
+        },
+        app // TODO: Remove me once fixed
+        ));
+  }
+
   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
       "EcalEndcapPRecHits", {"EcalEndcapPRawHits"}, {"EcalEndcapPRecHits"},
       {
@@ -63,8 +127,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