detectors/BEMC/BEMC.cc

0001 // Copyright 2022, David Lawrence
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004 //
0005
0006 #include <Evaluator/DD4hepUnits.h>
0007 #include <JANA/JApplicationFwd.h>
0008 #include <edm4eic/EDM4eicVersion.h>
0009 #include <edm4eic/unit_system.h>
0010 #include <cmath>
0011 #include <string>
0012 #include <variant>
0013
0014 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
0015 #include "algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h"
0016 #include "extensions/jana/JOmniFactoryGeneratorT.h"
0017 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
0018 #include "factories/calorimetry/CalorimeterClusterShape_factory.h"
0019 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
0020 #include "factories/calorimetry/CalorimeterHitReco_factory.h"
0021 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
0022 #include "factories/calorimetry/EnergyPositionClusterMerger_factory.h"
0023 #include "factories/calorimetry/ImagingClusterReco_factory.h"
0024 #include "factories/calorimetry/ImagingTopoCluster_factory.h"
0025 #include "factories/calorimetry/SimCalorimeterHitProcessor_factory.h"
0026 #include "factories/calorimetry/TruthEnergyPositionClusterMerger_factory.h"
0027
0028 extern "C" {
0029 void InitPlugin(JApplication* app) {
0030
0031   using namespace eicrecon;
0032
0033   InitJANAPlugin(app);
0034
0035   // Make sure left and right use the same value
0036   decltype(SimCalorimeterHitProcessorConfig::attenuationParameters) EcalBarrelScFi_attPars = {
0037       0.416212, 747.39875 * edm4eic::unit::mm, 7521.88383 * edm4eic::unit::mm};
0038   decltype(SimCalorimeterHitProcessorConfig::hitMergeFields) EcalBarrelScFi_hitMergeFields = {
0039       "fiber", "z"};
0040   decltype(SimCalorimeterHitProcessorConfig::contributionMergeFields)
0041       EcalBarrelScFi_contributionMergeFields = {"fiber"};
0042   decltype(SimCalorimeterHitProcessorConfig::inversePropagationSpeed)
0043       EcalBarrelScFi_inversePropagationSpeed = {(1. / 160) * edm4eic::unit::ns / edm4eic::unit::mm};
0044   decltype(SimCalorimeterHitProcessorConfig::fixedTimeDelay) EcalBarrelScFi_fixedTimeDelay = {
0045       2 * edm4eic::unit::ns};
0046
0047   // Make sure digi and reco use the same value
0048   decltype(CalorimeterHitDigiConfig::capADC) EcalBarrelScFi_capADC = 16384; //16384,  14bit ADC
0049   decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalBarrelScFi_dyRangeADC   = 1500 * dd4hep::MeV;
0050   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalBarrelScFi_pedMeanADC   = 100;
0051   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalBarrelScFi_pedSigmaADC = 1;
0052   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalBarrelScFi_resolutionTDC =
0053       10 * dd4hep::picosecond;
0054   app->Add(new JOmniFactoryGeneratorT<SimCalorimeterHitProcessor_factory>(
0055       "EcalBarrelScFiPAttenuatedHits", {"EcalBarrelScFiHits"},
0056       {"EcalBarrelScFiPAttenuatedHits", "EcalBarrelScFiPAttenuatedHitContributions"},
0057       {
0058           .attenuationParameters            = EcalBarrelScFi_attPars,
0059           .readout                          = "EcalBarrelScFiHits",
0060           .attenuationReferencePositionName = "EcalBarrel_Readout_zmax",
0061           .hitMergeFields                   = EcalBarrelScFi_hitMergeFields,
0062           .contributionMergeFields          = EcalBarrelScFi_contributionMergeFields,
0063           .inversePropagationSpeed          = EcalBarrelScFi_inversePropagationSpeed,
0064           .fixedTimeDelay                   = EcalBarrelScFi_fixedTimeDelay,
0065       },
0066       app // TODO: Remove me once fixed
0067       ));
0068   app->Add(new JOmniFactoryGeneratorT<SimCalorimeterHitProcessor_factory>(
0069       "EcalBarrelScFiNAttenuatedHits", {"EcalBarrelScFiHits"},
0070       {"EcalBarrelScFiNAttenuatedHits", "EcalBarrelScFiNAttenuatedHitContributions"},
0071       {
0072           .attenuationParameters            = EcalBarrelScFi_attPars,
0073           .readout                          = "EcalBarrelScFiHits",
0074           .attenuationReferencePositionName = "EcalBarrel_Readout_zmin",
0075           .hitMergeFields                   = EcalBarrelScFi_hitMergeFields,
0076           .contributionMergeFields          = EcalBarrelScFi_contributionMergeFields,
0077           .inversePropagationSpeed          = EcalBarrelScFi_inversePropagationSpeed,
0078           .fixedTimeDelay                   = EcalBarrelScFi_fixedTimeDelay,
0079       },
0080       app // TODO: Remove me once fixed
0081       ));
0082   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0083       "EcalBarrelScFiRawHits", {"EcalBarrelScFiHits"},
0084 #if EDM4EIC_VERSION_MAJOR >= 7
0085       {"EcalBarrelScFiRawHits", "EcalBarrelScFiRawHitAssociations"},
0086 #else
0087       {"EcalBarrelScFiRawHits"},
0088 #endif
0089       {
0090           .eRes          = {0.0 * sqrt(dd4hep::GeV), 0.0, 0.0 * dd4hep::GeV},
0091           .tRes          = 0.0 * dd4hep::ns,
0092           .threshold     = 0.0 * dd4hep::keV, // threshold is set in ADC in reco
0093           .capADC        = EcalBarrelScFi_capADC,
0094           .dyRangeADC    = EcalBarrelScFi_dyRangeADC,
0095           .pedMeanADC    = EcalBarrelScFi_pedMeanADC,
0096           .pedSigmaADC   = EcalBarrelScFi_pedSigmaADC,
0097           .resolutionTDC = EcalBarrelScFi_resolutionTDC,
0098           .corrMeanScale = "1.0",
0099           .readout       = "EcalBarrelScFiHits",
0100           .fields        = {"fiber", "z"},
0101       },
0102       app // TODO: Remove me once fixed
0103       ));
0104   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0105       "EcalBarrelScFiRecHits", {"EcalBarrelScFiRawHits"}, {"EcalBarrelScFiRecHits"},
0106       {
0107           .capADC          = EcalBarrelScFi_capADC,
0108           .dyRangeADC      = EcalBarrelScFi_dyRangeADC,
0109           .pedMeanADC      = EcalBarrelScFi_pedMeanADC,
0110           .pedSigmaADC     = EcalBarrelScFi_pedSigmaADC, // not needed; use only thresholdValue
0111           .resolutionTDC   = EcalBarrelScFi_resolutionTDC,
0112           .thresholdFactor = 0.0, // use only thresholdValue
0113           .thresholdValue  = 5.0, // 16384 ADC counts/1500 MeV * 0.5 MeV (desired threshold) = 5.46
0114           .sampFrac        = "0.09285755",
0115           .readout         = "EcalBarrelScFiHits",
0116           .layerField      = "layer",
0117           .sectorField     = "sector",
0118           .localDetFields  = {"system"},
0119           // here we want to use grid center position (XY) but keeps the z information from fiber-segment
0120           // TODO: a more realistic way to get z is to reconstruct it from timing
0121           .maskPos       = "xy",
0122           .maskPosFields = {"fiber", "z"},
0123       },
0124       app // TODO: Remove me once fixed
0125       ));
0126   app->Add(new JOmniFactoryGeneratorT<CalorimeterIslandCluster_factory>(
0127       "EcalBarrelScFiProtoClusters", {"EcalBarrelScFiRecHits"}, {"EcalBarrelScFiProtoClusters"},
0128       {
0129           .adjacencyMatrix{},
0130           .peakNeighbourhoodMatrix{},
0131           .readout{},
0132           .sectorDist = 50. * dd4hep::mm,
0133           .localDistXY{},
0134           .localDistXZ = {80 * dd4hep::mm, 80 * dd4hep::mm},
0135           .localDistYZ{},
0136           .globalDistRPhi{},
0137           .globalDistEtaPhi{},
0138           .dimScaledLocalDistXY{},
0139           .splitCluster         = false,
0140           .minClusterHitEdep    = 5.0 * dd4hep::MeV,
0141           .minClusterCenterEdep = 100.0 * dd4hep::MeV,
0142           .transverseEnergyProfileMetric{},
0143           .transverseEnergyProfileScale{},
0144           .transverseEnergyProfileScaleUnits{},
0145       },
0146       app // TODO: Remove me once fixed
0147       ));
0148   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0149       "EcalBarrelScFiClustersWithoutShapes",
0150       {
0151         "EcalBarrelScFiProtoClusters", // edm4eic::ProtoClusterCollection
0152 #if EDM4EIC_VERSION_MAJOR >= 7
0153             "EcalBarrelScFiRawHitAssociations"
0154       }, // edm4eic::MCRecoCalorimeterHitAssociation
0155 #else
0156             "EcalBarrelScFiHits"
0157       }, // edm4hep::SimCalorimeterHitCollection
0158 #endif
0159       {"EcalBarrelScFiClustersWithoutShapes",             // edm4eic::Cluster
0160        "EcalBarrelScFiClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
0161       {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, .enableEtaBounds = false},
0162       app // TODO: Remove me once fixed
0163       ));
0164   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
0165       "EcalBarrelScFiClusters",
0166       {"EcalBarrelScFiClustersWithoutShapes", "EcalBarrelScFiClusterAssociationsWithoutShapes"},
0167       {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations"},
0168       {.longitudinalShowerInfoAvailable = true, .energyWeight = "log", .logWeightBase = 6.2}, app));
0169
0170   // Make sure digi and reco use the same value
0171   decltype(CalorimeterHitDigiConfig::capADC) EcalBarrelImaging_capADC = 8192; //8192,  13bit ADC
0172   decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalBarrelImaging_dyRangeADC = 3 * dd4hep::MeV;
0173   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalBarrelImaging_pedMeanADC =
0174       14; // Noise floor at 5 keV: 8192 / 3 * 0.005
0175   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalBarrelImaging_pedSigmaADC =
0176       5; // Upper limit for sigma for AstroPix
0177   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalBarrelImaging_resolutionTDC =
0178       3.25 * dd4hep::nanosecond;
0179   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0180       "EcalBarrelImagingRawHits", {"EcalBarrelImagingHits"},
0181 #if EDM4EIC_VERSION_MAJOR >= 7
0182       {"EcalBarrelImagingRawHits", "EcalBarrelImagingRawHitAssociations"},
0183 #else
0184       {"EcalBarrelImagingRawHits"},
0185 #endif
0186       {
0187           .eRes          = {0.0 * sqrt(dd4hep::GeV), 0.02, 0.0 * dd4hep::GeV},
0188           .tRes          = 0.0 * dd4hep::ns,
0189           .capADC        = EcalBarrelImaging_capADC,
0190           .dyRangeADC    = EcalBarrelImaging_dyRangeADC,
0191           .pedMeanADC    = EcalBarrelImaging_pedMeanADC,
0192           .pedSigmaADC   = EcalBarrelImaging_pedSigmaADC,
0193           .resolutionTDC = EcalBarrelImaging_resolutionTDC,
0194           .corrMeanScale = "1.0",
0195           .readout       = "EcalBarrelImagingHits",
0196       },
0197       app // TODO: Remove me once fixed
0198       ));
0199   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0200       "EcalBarrelImagingRecHits", {"EcalBarrelImagingRawHits"}, {"EcalBarrelImagingRecHits"},
0201       {
0202           .capADC          = EcalBarrelImaging_capADC,
0203           .dyRangeADC      = EcalBarrelImaging_dyRangeADC,
0204           .pedMeanADC      = EcalBarrelImaging_pedMeanADC,
0205           .pedSigmaADC     = EcalBarrelImaging_pedSigmaADC, // not needed; use only thresholdValue
0206           .resolutionTDC   = EcalBarrelImaging_resolutionTDC,
0207           .thresholdFactor = 0.0, // use only thresholdValue
0208           .thresholdValue  = 41,  // 8192 ADC counts/3 MeV * 0.015 MeV (desired threshold) = 41
0209           .sampFrac        = "0.00429453",
0210           .readout         = "EcalBarrelImagingHits",
0211           .layerField      = "layer",
0212           .sectorField     = "sector",
0213       },
0214       app // TODO: Remove me once fixed
0215       ));
0216   app->Add(new JOmniFactoryGeneratorT<ImagingTopoCluster_factory>(
0217       "EcalBarrelImagingProtoClusters", {"EcalBarrelImagingRecHits"},
0218       {"EcalBarrelImagingProtoClusters"},
0219       {
0220           .neighbourLayersRange = 2, //  # id diff for adjacent layer
0221           .localDistXY          = {2.0 * dd4hep::mm, 2 * dd4hep::mm},     //  # same layer
0222           .layerDistEtaPhi      = {10 * dd4hep::mrad, 10 * dd4hep::mrad}, //  # adjacent layer
0223           .sectorDist           = 3.0 * dd4hep::cm,
0224           .minClusterHitEdep    = 0,
0225           .minClusterCenterEdep = 0,
0226           .minClusterEdep       = 100 * dd4hep::MeV,
0227           .minClusterNhits      = 10,
0228       },
0229       app // TODO: Remove me once fixed
0230       ));
0231
0232   app->Add(new JOmniFactoryGeneratorT<ImagingClusterReco_factory>(
0233       "EcalBarrelImagingClusters",
0234       {
0235         "EcalBarrelImagingProtoClusters",
0236 #if EDM4EIC_VERSION_MAJOR >= 7
0237             "EcalBarrelImagingRawHitAssociations"
0238       },
0239 #else
0240             "EcalBarrelImagingHits"
0241       },
0242 #endif
0243       {"EcalBarrelImagingClusters", "EcalBarrelImagingClusterAssociations",
0244        "EcalBarrelImagingLayers"},
0245       {
0246           .trackStopLayer = 6,
0247       },
0248       app // TODO: Remove me once fixed
0249       ));
0250   app->Add(new JOmniFactoryGeneratorT<EnergyPositionClusterMerger_factory>(
0251       "EcalBarrelClusters",
0252       {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations", "EcalBarrelImagingClusters",
0253        "EcalBarrelImagingClusterAssociations"},
0254       {"EcalBarrelClusters", "EcalBarrelClusterAssociations"},
0255       {
0256           .energyRelTolerance = 0.5,
0257           .phiTolerance       = 0.1,
0258           .etaTolerance       = 0.2,
0259       },
0260       app // TODO: Remove me once fixed
0261       ));
0262   app->Add(new JOmniFactoryGeneratorT<TruthEnergyPositionClusterMerger_factory>(
0263       "EcalBarrelTruthClusters",
0264       {"MCParticles", "EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations",
0265        "EcalBarrelImagingClusters", "EcalBarrelImagingClusterAssociations"},
0266       {"EcalBarrelTruthClusters", "EcalBarrelTruthClusterAssociations"},
0267       app // TODO: Remove me once fixed
0268       ));
0269 }
0270 }