detectors/FEMC/FEMC.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2021 - 2025, Chao Peng, Sylvester Joosten, Whitney Armstrong, David Lawrence, Friederike Bock, Wouter Deconinck, Kolja Kauder, Sebouh Paul, Akio Ogawa
0003
0004 #include <DD4hep/Detector.h>
0005 #include <Evaluator/DD4hepUnits.h>
0006 #include <JANA/JApplication.h>
0007 #include <JANA/JApplicationFwd.h>
0008 #include <JANA/Utils/JTypeInfo.h>
0009 #include <edm4eic/EDM4eicVersion.h>
0010 #include <fmt/format.h>
0011 #include <spdlog/logger.h>
0012 #include <cmath>
0013 #include <gsl/pointers>
0014 #include <memory>
0015 #include <string>
0016 #include <variant>
0017 #include <vector>
0018
0019 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
0020 #include "extensions/jana/JOmniFactoryGeneratorT.h"
0021 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
0022 #include "factories/calorimetry/CalorimeterClusterShape_factory.h"
0023 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
0024 #include "factories/calorimetry/CalorimeterHitReco_factory.h"
0025 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
0026 #include "factories/calorimetry/CalorimeterTruthClustering_factory.h"
0027 #include "factories/calorimetry/TrackClusterMergeSplitter_factory.h"
0028 #include "services/geometry/dd4hep/DD4hep_service.h"
0029 #include "services/log/Log_service.h"
0030
0031 extern "C" {
0032 void InitPlugin(JApplication* app) {
0033
0034   using namespace eicrecon;
0035
0036   InitJANAPlugin(app);
0037
0038   auto log_service = app->GetService<Log_service>();
0039   auto mLog        = log_service->logger("FEMC");
0040
0041   // Make sure digi and reco use the same value
0042   decltype(CalorimeterHitDigiConfig::capADC) EcalEndcapP_capADC =
0043       16384; //16384, assuming 14 bits. For approximate HGCROC resolution use 65536
0044   decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalEndcapP_dyRangeADC   = 100 * dd4hep::GeV;
0045   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalEndcapP_pedMeanADC   = 200;
0046   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalEndcapP_pedSigmaADC = 2.4576;
0047   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalEndcapP_resolutionTDC =
0048       10 * dd4hep::picosecond;
0049   const double EcalEndcapP_sampFrac = 0.029043; // updated with ratio to ScFi model
0050   decltype(CalorimeterHitDigiConfig::corrMeanScale) EcalEndcapP_corrMeanScale =
0051       fmt::format("{}", 1.0 / EcalEndcapP_sampFrac); //only used for ScFi model
0052   const double EcalEndcapP_nPhotonPerGeV          = 1500;
0053   const double EcalEndcapP_PhotonCollectionEff    = 0.285;
0054   const unsigned long long EcalEndcapP_totalPixel = 4 * 159565ULL;
0055
0056   int EcalEndcapP_homogeneousFlag = 0;
0057   try {
0058     auto detector               = app->GetService<DD4hep_service>()->detector();
0059     EcalEndcapP_homogeneousFlag = detector->constant<int>("EcalEndcapP_Homogeneous_ScFi");
0060     if (EcalEndcapP_homogeneousFlag <= 1) {
0061       mLog->info("Homogeneous geometry loaded");
0062     } else {
0063       mLog->info("ScFi geometry loaded");
0064     }
0065   } catch (...) {
0066     // Variable not present apply legacy homogeneous geometry implementation
0067     EcalEndcapP_homogeneousFlag = 0;
0068   };
0069
0070   if (EcalEndcapP_homogeneousFlag <= 1) {
0071     app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0072         "EcalEndcapPRawHits", {"EventHeader", "EcalEndcapPHits"},
0073         {"EcalEndcapPRawHits",
0074 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0075          "EcalEndcapPRawHitLinks",
0076 #endif
0077          "EcalEndcapPRawHitAssociations"},
0078         {
0079             .eRes = {0.11333 * sqrt(dd4hep::GeV), 0.03,
0080                      0.0 * dd4hep::GeV}, // (11.333% / sqrt(E)) \oplus 3%
0081             .tRes = 0.0,
0082             .threshold =
0083                 0.0, // 15MeV threshold for a single tower will be applied on ADC at Reco below
0084             .readoutType = "sipm",
0085             .lightYield  = EcalEndcapP_nPhotonPerGeV / EcalEndcapP_PhotonCollectionEff,
0086             .photonDetectionEfficiency = EcalEndcapP_PhotonCollectionEff,
0087             .numEffectiveSipmPixels    = EcalEndcapP_totalPixel,
0088             .capADC                    = EcalEndcapP_capADC,
0089             .capTime                   = 100, // given in ns, 4 samples in HGCROC
0090             .dyRangeADC                = EcalEndcapP_dyRangeADC,
0091             .pedMeanADC                = EcalEndcapP_pedMeanADC,
0092             .pedSigmaADC               = EcalEndcapP_pedSigmaADC,
0093             .resolutionTDC             = EcalEndcapP_resolutionTDC,
0094             .corrMeanScale             = "1.0",
0095             .readout                   = "EcalEndcapPHits",
0096         },
0097         app // TODO: Remove me once fixed
0098         ));
0099   } else if (EcalEndcapP_homogeneousFlag == 2) {
0100     app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0101         "EcalEndcapPRawHits", {"EventHeader", "EcalEndcapPHits"},
0102         {"EcalEndcapPRawHits",
0103 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0104          "EcalEndcapPRawHitLinks",
0105 #endif
0106          "EcalEndcapPRawHitAssociations"},
0107         {
0108             .eRes = {0.0, 0.022, 0.0}, // just constant term 2.2% based on MC data comparison
0109             .tRes = 0.0,
0110             .threshold =
0111                 0.0, // 15MeV threshold for a single tower will be applied on ADC at Reco below
0112             .readoutType = "sipm",
0113             .lightYield =
0114                 EcalEndcapP_nPhotonPerGeV / EcalEndcapP_PhotonCollectionEff / EcalEndcapP_sampFrac,
0115             .photonDetectionEfficiency = EcalEndcapP_PhotonCollectionEff,
0116             .numEffectiveSipmPixels    = EcalEndcapP_totalPixel,
0117             .capADC                    = EcalEndcapP_capADC,
0118             .capTime                   = 100, // given in ns, 4 samples in HGCROC
0119             .dyRangeADC                = EcalEndcapP_dyRangeADC,
0120             .pedMeanADC                = EcalEndcapP_pedMeanADC,
0121             .pedSigmaADC               = EcalEndcapP_pedSigmaADC,
0122             .resolutionTDC             = EcalEndcapP_resolutionTDC,
0123             .corrMeanScale             = EcalEndcapP_corrMeanScale,
0124             .readout                   = "EcalEndcapPHits",
0125             .fields                    = {"fiber_x", "fiber_y"},
0126         },
0127         app // TODO: Remove me once fixed
0128         ));
0129   }
0130
0131   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0132       "EcalEndcapPRecHits", {"EcalEndcapPRawHits"}, {"EcalEndcapPRecHits"},
0133       {
0134           .capADC          = EcalEndcapP_capADC,
0135           .dyRangeADC      = EcalEndcapP_dyRangeADC,
0136           .pedMeanADC      = EcalEndcapP_pedMeanADC,
0137           .pedSigmaADC     = EcalEndcapP_pedSigmaADC,
0138           .resolutionTDC   = EcalEndcapP_resolutionTDC,
0139           .thresholdFactor = 0.0,
0140           .thresholdValue =
0141               3, // The ADC of a 15 MeV particle is adc = 200 + 15 * 0.03 * ( 1.0 + 0) / 3000 * 16384 = 200 + 2.4576
0142           // 15 MeV = 2.4576, but adc=llround(dE) and cut off is "<". So 3 here = 15.25MeV
0143           .sampFrac = "1.00", // already taken care in DIGI code above
0144           .readout  = "EcalEndcapPHits",
0145       },
0146       app // TODO: Remove me once fixed
0147       ));
0148   app->Add(new JOmniFactoryGeneratorT<CalorimeterTruthClustering_factory>(
0149       "EcalEndcapPTruthProtoClusters", {"EcalEndcapPRecHits", "EcalEndcapPHits"},
0150       {"EcalEndcapPTruthProtoClusters"},
0151       app // TODO: Remove me once fixed
0152       ));
0153   app->Add(new JOmniFactoryGeneratorT<CalorimeterIslandCluster_factory>(
0154       "EcalEndcapPIslandProtoClusters", {"EcalEndcapPRecHits"}, {"EcalEndcapPIslandProtoClusters"},
0155       {.adjacencyMatrix{},
0156        .peakNeighbourhoodMatrix{},
0157        .readout{},
0158        .sectorDist = 5.0 * dd4hep::cm,
0159        .localDistXY{},
0160        .localDistXZ{},
0161        .localDistYZ{},
0162        .globalDistRPhi{},
0163        .globalDistEtaPhi{},
0164        .dimScaledLocalDistXY          = {1.5, 1.5},
0165        .splitCluster                  = false,
0166        .minClusterHitEdep             = 0.0 * dd4hep::MeV,
0167        .minClusterCenterEdep          = 60.0 * dd4hep::MeV,
0168        .transverseEnergyProfileMetric = "dimScaledLocalDistXY",
0169        .transverseEnergyProfileScale  = 1.,
0170        .transverseEnergyProfileScaleUnits{}},
0171       app // TODO: Remove me once fixed
0172       ));
0173
0174   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0175       "EcalEndcapPTruthClustersWithoutShapes",
0176       {
0177           "EcalEndcapPTruthProtoClusters", // edm4eic::ProtoClusterCollection
0178 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0179           "EcalEndcapPRawHitLinks", // edm4eic::MCRecoCalorimeterHitLink
0180 #endif
0181           "EcalEndcapPRawHitAssociations" // edm4eic::MCRecoCalorimeterHitAssociationCollection
0182       },
0183       {"EcalEndcapPTruthClustersWithoutShapes",
0184 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0185        "EcalEndcapPTruthClusterLinksWithoutShapes",
0186 #endif
0187        "EcalEndcapPTruthClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
0188       {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, .enableEtaBounds = true},
0189       app // TODO: Remove me once fixed
0190       ));
0191
0192   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
0193       "EcalEndcapPTruthClusters",
0194       {"EcalEndcapPTruthClustersWithoutShapes", "EcalEndcapPTruthClusterAssociationsWithoutShapes"},
0195       {"EcalEndcapPTruthClusters",
0196 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0197        "EcalEndcapPTruthClusterLinks",
0198 #endif
0199        "EcalEndcapPTruthClusterAssociations"},
0200       {.energyWeight = "log", .logWeightBase = 6.2}, app));
0201
0202   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0203       "EcalEndcapPClustersWithoutShapes",
0204       {
0205           "EcalEndcapPIslandProtoClusters", // edm4eic::ProtoClusterCollection
0206 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0207           "EcalEndcapPRawHitLinks", // edm4eic::MCRecoCalorimeterHitLink
0208 #endif
0209           "EcalEndcapPRawHitAssociations" // edm4eic::MCRecoCalorimeterHitAssociationCollection
0210       },
0211       {"EcalEndcapPClustersWithoutShapes",
0212 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0213        "EcalEndcapPClusterLinksWithoutShapes",
0214 #endif
0215        "EcalEndcapPClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
0216       {
0217           .energyWeight    = "log",
0218           .sampFrac        = 1.0,
0219           .logWeightBase   = 3.6,
0220           .enableEtaBounds = false,
0221       },
0222       app // TODO: Remove me once fixed
0223       ));
0224
0225   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
0226       "EcalEndcapPClusters",
0227       {"EcalEndcapPClustersWithoutShapes", "EcalEndcapPClusterAssociationsWithoutShapes"},
0228       {"EcalEndcapPClusters",
0229 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0230        "EcalEndcapPClusterLinks",
0231 #endif
0232        "EcalEndcapPClusterAssociations"},
0233       {.energyWeight = "log", .logWeightBase = 3.6}, app));
0234
0235   app->Add(new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
0236       "EcalEndcapPSplitMergeProtoClusters",
0237       {"EcalEndcapPIslandProtoClusters", "CalorimeterTrackProjections"},
0238       {"EcalEndcapPSplitMergeProtoClusters"},
0239       {.idCalo                       = "EcalEndcapP_ID",
0240        .minSigCut                    = -2.0,
0241        .avgEP                        = 1.0,
0242        .sigEP                        = 0.10,
0243        .drAdd                        = 0.30,
0244        .sampFrac                     = 1.0,
0245        .transverseEnergyProfileScale = 1.0},
0246       app // TODO: remove me once fixed
0247       ));
0248
0249   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0250       "EcalEndcapPSplitMergeClustersWithoutShapes",
0251       {
0252           "EcalEndcapPSplitMergeProtoClusters", // edm4eic::ProtoClusterCollection
0253 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0254           "EcalEndcapPRawHitLinks", // edm4eic::MCRecoCalorimeterHitLink
0255 #endif
0256           "EcalEndcapPRawHitAssociations" // edm4hep::MCRecoCalorimeterHitAssociationCollection
0257       },
0258       {"EcalEndcapPSplitMergeClustersWithoutShapes",
0259 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0260        "EcalEndcapPSplitMergeClusterLinksWithoutShapes",
0261 #endif
0262        "EcalEndcapPSplitMergeClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
0263       {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 3.6, .enableEtaBounds = false},
0264       app // TODO: Remove me once fixed
0265       ));
0266
0267   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
0268       "EcalEndcapPSplitMergeClusters",
0269       {"EcalEndcapPSplitMergeClustersWithoutShapes",
0270        "EcalEndcapPSplitMergeClusterAssociationsWithoutShapes"},
0271       {"EcalEndcapPSplitMergeClusters",
0272 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0273        "EcalEndcapPSplitMergeClusterLinks",
0274 #endif
0275        "EcalEndcapPSplitMergeClusterAssociations"},
0276       {.energyWeight = "log", .logWeightBase = 3.6}, app));
0277 }
0278 }