diff --git a/src/detectors/FEMC/FEMC.cc b/src/detectors/FEMC/FEMC.cc index a8d4a3b4aa..5b9551fc13 100644 --- a/src/detectors/FEMC/FEMC.cc +++ b/src/detectors/FEMC/FEMC.cc @@ -1,5 +1,5 @@ // 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 #include @@ -8,8 +8,10 @@ #include #include #include +#include #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h" +#include "algorithms/calorimetry/CalorimeterHitRecoConfig.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h" #include "factories/calorimetry/CalorimeterClusterShape_factory.h" @@ -25,34 +27,85 @@ void InitPlugin(JApplication* app) { using namespace eicrecon; InitJANAPlugin(app); + + auto log_service = app->GetService(); + 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( - "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 sampFrac = 0.029043; // updated with ratio to ScFi model + decltype(CalorimeterHitDigiConfig::corrMeanScale) EcalEndcapP_corrMeanScale = + fmt::format("{}", 1.0 / sampFrac); + decltype(CalorimeterHitRecoConfig::sampFrac) EcalEndcapP_sampFrac = fmt::format("{}", sampFrac); + const double EcalEndcapP_nPhotonPerGeV = 1500; + const double EcalEndcapP_PhotonCollectionEff = 0.5; + const unsigned long long EcalEndcapP_totalPixel = 4 * 159565; + + int FEMCHomoScfi = 0; + try { + auto detector = app->GetService()->detector(); + FEMCHomoScfi = detector->constant("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( + "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 + .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", + .readoutType = "sipm", + .lightYield = EcalEndcapP_nPhotonPerGeV / EcalEndcapP_PhotonCollectionEff, + .photonDetectionEfficiency = EcalEndcapP_PhotonCollectionEff, + .numEffectiveSipmPixels = EcalEndcapP_totalPixel}, + app // TODO: Remove me once fixed + )); + } else if (FEMCHomoScfi == 2) { + app->Add(new JOmniFactoryGeneratorT( + "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 + .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", + .readoutType = "sipm", + .fields = {"fiberx", "fibery", "x", "y"}, + .lightYield = EcalEndcapP_nPhotonPerGeV / EcalEndcapP_PhotonCollectionEff, + .photonDetectionEfficiency = EcalEndcapP_PhotonCollectionEff, + .numEffectiveSipmPixels = EcalEndcapP_totalPixel}, + app // TODO: Remove me once fixed + )); + } + app->Add(new JOmniFactoryGeneratorT( "EcalEndcapPRecHits", {"EcalEndcapPRawHits"}, {"EcalEndcapPRecHits"}, { @@ -63,8 +116,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