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