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