detectors/FHCAL/FHCAL.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Friederike Bock, Wouter Deconinck
0003
0004 #include <DD4hep/Detector.h>
0005 #include <edm4eic/EDM4eicVersion.h>
0006 #include <Evaluator/DD4hepUnits.h>
0007 #include <JANA/JApplication.h>
0008 #include <TString.h>
0009 #include <math.h>
0010 #include <algorithm>
0011 #include <gsl/pointers>
0012 #include <memory>
0013 #include <stdexcept>
0014 #include <string>
0015
0016 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
0017 #include "algorithms/calorimetry/ImagingTopoClusterConfig.h"
0018 #include "extensions/jana/JOmniFactoryGeneratorT.h"
0019 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
0020 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
0021 #include "factories/calorimetry/CalorimeterHitReco_factory.h"
0022 #include "factories/calorimetry/CalorimeterHitsMerger_factory.h"
0023 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
0024 #include "factories/calorimetry/CalorimeterTruthClustering_factory.h"
0025 #include "factories/calorimetry/HEXPLIT_factory.h"
0026 #include "factories/calorimetry/ImagingTopoCluster_factory.h"
0027 #include "factories/calorimetry/TrackClusterMergeSplitter_factory.h"
0028 #include "services/geometry/dd4hep/DD4hep_service.h"
0029
0030 extern "C" {
0031     void InitPlugin(JApplication *app) {
0032
0033         using namespace eicrecon;
0034
0035         InitJANAPlugin(app);
0036
0037         // Make sure digi and reco use the same value
0038         decltype(CalorimeterHitDigiConfig::capADC)        HcalEndcapPInsert_capADC = 32768;
0039         decltype(CalorimeterHitDigiConfig::dyRangeADC)    HcalEndcapPInsert_dyRangeADC = 200 * dd4hep::MeV;
0040         decltype(CalorimeterHitDigiConfig::pedMeanADC)    HcalEndcapPInsert_pedMeanADC = 10;
0041         decltype(CalorimeterHitDigiConfig::pedSigmaADC)   HcalEndcapPInsert_pedSigmaADC = 2;
0042         decltype(CalorimeterHitDigiConfig::resolutionTDC) HcalEndcapPInsert_resolutionTDC = 10 * dd4hep::picosecond;
0043
0044         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0045            "HcalEndcapPInsertRawHits",
0046            {"HcalEndcapPInsertHits"},
0047 #if EDM4EIC_VERSION_MAJOR >= 7
0048            {"HcalEndcapPInsertRawHits", "HcalEndcapPInsertRawHitAssociations"},
0049 #else
0050            {"HcalEndcapPInsertRawHits"},
0051 #endif
0052            {
0053              .eRes = {},
0054              .tRes = 0.0 * dd4hep::ns,
0055              .capADC = HcalEndcapPInsert_capADC,
0056              .dyRangeADC = HcalEndcapPInsert_dyRangeADC,
0057              .pedMeanADC = HcalEndcapPInsert_pedMeanADC,
0058              .pedSigmaADC = HcalEndcapPInsert_pedSigmaADC,
0059              .resolutionTDC = HcalEndcapPInsert_resolutionTDC,
0060              .corrMeanScale = "1.0",
0061              .readout = "HcalEndcapPInsertHits",
0062            },
0063           app   // TODO: Remove me once fixed
0064         ));
0065         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0066           "HcalEndcapPInsertRecHits", {"HcalEndcapPInsertRawHits"}, {"HcalEndcapPInsertRecHits"},
0067           {
0068             .capADC = HcalEndcapPInsert_capADC,
0069             .dyRangeADC = HcalEndcapPInsert_dyRangeADC,
0070             .pedMeanADC = HcalEndcapPInsert_pedMeanADC,
0071             .pedSigmaADC = HcalEndcapPInsert_pedSigmaADC,
0072             .resolutionTDC = HcalEndcapPInsert_resolutionTDC,
0073             .thresholdFactor = 0.,
0074             .thresholdValue = 41.0, // 0.25 MeV --> 0.25 / 200 * 32768 = 41
0075
0076             .sampFrac = "1.0",
0077             .readout = "HcalEndcapPInsertHits",
0078             .layerField="layer",
0079           },
0080           app   // TODO: Remove me once fixed
0081         ));
0082         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitsMerger_factory>(
0083           "HcalEndcapPInsertMergedHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertMergedHits"},
0084           {
0085             .readout = "HcalEndcapPInsertHits",
0086             .fields = {"layer", "slice"},
0087             .refs = {1, 0},
0088           },
0089           app   // TODO: Remove me once fixed
0090         ));
0091         app->Add(new JOmniFactoryGeneratorT<CalorimeterTruthClustering_factory>(
0092           "HcalEndcapPInsertTruthProtoClusters", {"HcalEndcapPInsertMergedHits", "HcalEndcapPInsertHits"}, {"HcalEndcapPInsertTruthProtoClusters"},
0093           app   // TODO: Remove me once fixed
0094         ));
0095
0096       app->Add(new JOmniFactoryGeneratorT<HEXPLIT_factory>(
0097         "HcalEndcapPInsertSubcellHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertSubcellHits"},
0098         {
0099           .MIP = 480. * dd4hep::keV,
0100           .Emin_in_MIPs=0.5,
0101           .tmax=162 * dd4hep::ns, //150 ns + (z at front face)/(speed of light)
0102         },
0103         app   // TODO: Remove me once fixed
0104       ));
0105
0106       // define the distance between neighbors in terms of the largest possible distance between subcell hits
0107       auto detector = app->GetService<DD4hep_service>()->detector();
0108       double side_length;
0109       try {
0110         side_length = std::max({detector->constant<double>("HcalEndcapPInsertCellSizeLGRight"), detector->constant<double>("HcalEndcapPInsertCellSizeLGLeft")});
0111       } catch (std::runtime_error&) {
0112         side_length = 31. * dd4hep::mm;
0113       }
0114       app->Add(new JOmniFactoryGeneratorT<ImagingTopoCluster_factory>(
0115           "HcalEndcapPInsertImagingProtoClusters", {"HcalEndcapPInsertSubcellHits"}, {"HcalEndcapPInsertImagingProtoClusters"},
0116           {
0117               .neighbourLayersRange = 1,
0118               .localDistXY = {0.5*side_length, 0.5*side_length*sin(M_PI/3)},
0119               .layerDistXY = {0.25*side_length, 0.25*side_length*sin(M_PI/3)},
0120               .layerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xy,
0121               .sectorDist = 10.0 * dd4hep::cm,
0122               .minClusterHitEdep = 5.0 * dd4hep::keV,
0123               .minClusterCenterEdep = 3.0 * dd4hep::MeV,
0124               .minClusterEdep = 11.0 * dd4hep::MeV,
0125               .minClusterNhits = 100,
0126           },
0127           app   // TODO: Remove me once fixed
0128       ));
0129
0130         app->Add(
0131           new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0132              "HcalEndcapPInsertTruthClusters",
0133             {"HcalEndcapPInsertTruthProtoClusters",        // edm4eic::ProtoClusterCollection
0134 #if EDM4EIC_VERSION_MAJOR >= 7
0135              "HcalEndcapPInsertRawHitAssociations"},       // edm4eic::MCRecoCalorimeterHitAssociationCollection
0136 #else
0137              "HcalEndcapPInsertHits"},                     // edm4hep::SimCalorimeterHitCollection
0138 #endif
0139             {"HcalEndcapPInsertTruthClusters",             // edm4eic::Cluster
0140              "HcalEndcapPInsertTruthClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation
0141             {
0142               .energyWeight = "log",
0143               .sampFrac = 0.0257,
0144               .logWeightBase = 3.6,
0145               .longitudinalShowerInfoAvailable = true,
0146               .enableEtaBounds = true
0147             },
0148             app   // TODO: Remove me once fixed
0149           )
0150         );
0151
0152         app->Add(
0153           new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0154              "HcalEndcapPInsertClusters",
0155             {"HcalEndcapPInsertImagingProtoClusters",  // edm4eic::ProtoClusterCollection
0156 #if EDM4EIC_VERSION_MAJOR >= 7
0157              "HcalEndcapPInsertRawHitAssociations"},  // edm4eic::MCRecoCalorimeterHitAssociationCollection
0158 #else
0159              "HcalEndcapPInsertHits"},                // edm4hep::SimCalorimeterHitCollection
0160 #endif
0161             {"HcalEndcapPInsertClusters",             // edm4eic::Cluster
0162              "HcalEndcapPInsertClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation
0163             {
0164               .energyWeight = "log",
0165               .sampFrac = 0.0257,
0166               .logWeightBase = 6.2,
0167               .longitudinalShowerInfoAvailable = true,
0168               .enableEtaBounds = false,
0169             },
0170             app   // TODO: Remove me once fixed
0171           )
0172         );
0173
0174         // Make sure digi and reco use the same value
0175         decltype(CalorimeterHitDigiConfig::capADC)        LFHCAL_capADC = 65536;
0176         decltype(CalorimeterHitDigiConfig::dyRangeADC)    LFHCAL_dyRangeADC = 1 * dd4hep::GeV;
0177         decltype(CalorimeterHitDigiConfig::pedMeanADC)    LFHCAL_pedMeanADC = 50;
0178         decltype(CalorimeterHitDigiConfig::pedSigmaADC)   LFHCAL_pedSigmaADC = 10;
0179         decltype(CalorimeterHitDigiConfig::resolutionTDC) LFHCAL_resolutionTDC = 10 * dd4hep::picosecond;
0180
0181         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0182           "LFHCALRawHits",
0183           {"LFHCALHits"},
0184 #if EDM4EIC_VERSION_MAJOR >= 7
0185           {"LFHCALRawHits", "LFHCALRawHitAssociations"},
0186 #else
0187           {"LFHCALRawHits"},
0188 #endif
0189           {
0190             .eRes = {},
0191             .tRes = 0.0 * dd4hep::ns,
0192             .capADC = LFHCAL_capADC,
0193             .capTime = 100,
0194             .dyRangeADC = LFHCAL_dyRangeADC,
0195             .pedMeanADC = LFHCAL_pedMeanADC,
0196             .pedSigmaADC = LFHCAL_pedSigmaADC,
0197             .resolutionTDC = LFHCAL_resolutionTDC,
0198             .corrMeanScale = "1.0",
0199             .readout = "LFHCALHits",
0200             .fields = {"layerz"},
0201           },
0202           app   // TODO: Remove me once fixed
0203         ));
0204         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0205           "LFHCALRecHits", {"LFHCALRawHits"}, {"LFHCALRecHits"},
0206           {
0207             .capADC = LFHCAL_capADC,
0208             .dyRangeADC = LFHCAL_dyRangeADC,
0209             .pedMeanADC = LFHCAL_pedMeanADC,
0210             .pedSigmaADC = LFHCAL_pedSigmaADC,
0211             .resolutionTDC = LFHCAL_resolutionTDC,
0212             .thresholdFactor = 0.0,
0213             .thresholdValue = 20, // 0.3 MeV deposition --> adc = 50 + 0.3 / 1000 * 65536 == 70
0214             .sampFrac = "(rlayerz == 0) ? 0.019 : 0.037", // 0.019 only in the 0-th tile
0215             .readout = "LFHCALHits",
0216             .layerField = "rlayerz",
0217           },
0218           app   // TODO: Remove me once fixed
0219         ));
0220         app->Add(new JOmniFactoryGeneratorT<CalorimeterTruthClustering_factory>(
0221           "LFHCALTruthProtoClusters", {"LFHCALRecHits", "LFHCALHits"}, {"LFHCALTruthProtoClusters"},
0222           app   // TODO: Remove me once fixed
0223         ));
0224
0225         // Magic constants:
0226         //  54 - number of modules in a row/column
0227         //  2  - number of towers in a module
0228         // sign for towerx and towery are *negative* to maintain linearity with global X and Y
0229         std::string cellIdx_1  = "(54*2-moduleIDx_1*2-towerx_1)";
0230         std::string cellIdx_2  = "(54*2-moduleIDx_2*2-towerx_2)";
0231         std::string cellIdy_1  = "(54*2-moduleIDy_1*2-towery_1)";
0232         std::string cellIdy_2  = "(54*2-moduleIDy_2*2-towery_2)";
0233         std::string cellIdz_1  = "rlayerz_1";
0234         std::string cellIdz_2  = "rlayerz_2";
0235         std::string deltaX     = Form("abs(%s-%s)", cellIdx_2.data(), cellIdx_1.data());
0236         std::string deltaY     = Form("abs(%s-%s)", cellIdy_2.data(), cellIdy_1.data());
0237         std::string deltaZ     = Form("abs(%s-%s)", cellIdz_2.data(), cellIdz_1.data());
0238         std::string neighbor   = Form("(%s+%s+%s==1)", deltaX.data(), deltaY.data(), deltaZ.data());
0239         std::string corner2D   = Form("((%s==0&&%s==1&&%s==1)||(%s==1&&%s==0&&%s==1)||(%s==1&&%s==1&&%s==0))",
0240                                   deltaZ.data(), deltaX.data(), deltaY.data(),
0241                                   deltaZ.data(), deltaX.data(), deltaY.data(),
0242                                   deltaZ.data(), deltaX.data(), deltaY.data());
0243
0244         app->Add(new JOmniFactoryGeneratorT<CalorimeterIslandCluster_factory>(
0245           "LFHCALIslandProtoClusters", {"LFHCALRecHits"}, {"LFHCALIslandProtoClusters"},
0246           {
0247             .adjacencyMatrix = Form("%s||%s", neighbor.data(), corner2D.data()),
0248             .readout = "LFHCALHits",
0249             .sectorDist = 0 * dd4hep::cm,
0250             .splitCluster = false,
0251             .minClusterHitEdep = 1 * dd4hep::MeV,
0252             .minClusterCenterEdep = 100.0 * dd4hep::MeV,
0253             .transverseEnergyProfileMetric = "globalDistEtaPhi",
0254             .transverseEnergyProfileScale = 1.,
0255           },
0256           app   // TODO: Remove me once fixed
0257         ));
0258
0259         app->Add(
0260           new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0261              "LFHCALTruthClusters",
0262             {"LFHCALTruthProtoClusters",        // edm4eic::ProtoClusterCollection
0263 #if EDM4EIC_VERSION_MAJOR >= 7
0264              "LFHCALRawHitAssociations"},       // edm4eic::MCRecoCalorimeterHitAssociationCollection
0265 #else
0266              "LFHCALHits"},                     // edm4hep::SimCalorimeterHitCollection
0267 #endif
0268             {"LFHCALTruthClusters",             // edm4eic::Cluster
0269              "LFHCALTruthClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation
0270             {
0271               .energyWeight = "log",
0272               .sampFrac = 1.0,
0273               .logWeightBase = 4.5,
0274               .longitudinalShowerInfoAvailable = true,
0275               .enableEtaBounds = false
0276             },
0277             app   // TODO: Remove me once fixed
0278           )
0279         );
0280
0281         app->Add(
0282           new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0283              "LFHCALClusters",
0284             {"LFHCALIslandProtoClusters",  // edm4eic::ProtoClusterCollection
0285 #if EDM4EIC_VERSION_MAJOR >= 7
0286              "LFHCALRawHitAssociations"},  // edm4eic::MCRecoCalorimeterHitAssociationCollection
0287 #else
0288              "LFHCALHits"},                // edm4hep::SimCalorimeterHitCollection
0289 #endif
0290             {"LFHCALClusters",             // edm4eic::Cluster
0291              "LFHCALClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation
0292             {
0293               .energyWeight = "log",
0294               .sampFrac = 1.0,
0295               .logWeightBase = 4.5,
0296               .longitudinalShowerInfoAvailable = true,
0297               .enableEtaBounds = false,
0298             },
0299             app   // TODO: Remove me once fixed
0300           )
0301         );
0302
0303         app->Add(
0304           new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
0305             "LFHCALSplitMergeProtoClusters",
0306             {"LFHCALIslandProtoClusters",
0307              "CalorimeterTrackProjections"},
0308             {"LFHCALSplitMergeProtoClusters"},
0309             {
0310               .idCalo = "LFHCAL_ID",
0311               .minSigCut = -2.0,
0312               .avgEP = 0.50,
0313               .sigEP = 0.25,
0314               .drAdd = 0.30,
0315               .sampFrac = 1.0,
0316               .transverseEnergyProfileScale = 1.0
0317             },
0318             app   // TODO: remove me once fixed
0319           )
0320         );
0321
0322         app->Add(
0323           new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0324              "LFHCALSplitMergeClusters",
0325             {"LFHCALSplitMergeProtoClusters",        // edm4eic::ProtoClusterCollection
0326              "LFHCALHits"},                          // edm4hep::SimCalorimeterHitCollection
0327             {"LFHCALSplitMergeClusters",             // edm4eic::Cluster
0328              "LFHCALSplitMergeClusterAssociations"}, // edm4eic::MCRecoClusterParticleAssociation
0329             {
0330               .energyWeight = "log",
0331               .sampFrac = 1.0,
0332               .logWeightBase = 4.5,
0333               .enableEtaBounds = false
0334             },
0335             app   // TODO: Remove me once fixed
0336           )
0337         );
0338     }
0339 }