Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-11-01 07:57:32

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 Whitney Armstrong, Sylvester Joosten, Chao Peng, David Lawrence, Thomas Britton, Wouter Deconinck, Maria Zurek, Akshaya Vijay, Nathan Brei, Dmitry Kalinkin, Derek Anderson, Minho Kim
0003 
0004 #include <Evaluator/DD4hepUnits.h>
0005 #include <JANA/JApplication.h>
0006 #include <JANA/JApplicationFwd.h>
0007 #include <JANA/Utils/JTypeInfo.h>
0008 #include <edm4eic/unit_system.h>
0009 #include <edm4hep/SimCalorimeterHit.h>
0010 #include <cmath>
0011 #include <map>
0012 #include <memory>
0013 #include <string>
0014 #include <variant>
0015 #include <vector>
0016 
0017 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
0018 #include "algorithms/calorimetry/ImagingTopoClusterConfig.h"
0019 #include "algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h"
0020 #include "algorithms/digi/PulseGenerationConfig.h"
0021 #include "algorithms/digi/PulseCombinerConfig.h"
0022 #include "algorithms/digi/PulseNoiseConfig.h"
0023 #include "extensions/jana/JOmniFactoryGeneratorT.h"
0024 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
0025 #include "factories/calorimetry/CalorimeterClusterShape_factory.h"
0026 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
0027 #include "factories/calorimetry/CalorimeterHitReco_factory.h"
0028 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
0029 #include "factories/calorimetry/EnergyPositionClusterMerger_factory.h"
0030 #include "factories/calorimetry/ImagingClusterReco_factory.h"
0031 #include "factories/calorimetry/ImagingTopoCluster_factory.h"
0032 #include "factories/calorimetry/SimCalorimeterHitProcessor_factory.h"
0033 #include "factories/calorimetry/TruthEnergyPositionClusterMerger_factory.h"
0034 #include "factories/digi/PulseGeneration_factory.h"
0035 #include "factories/digi/PulseCombiner_factory.h"
0036 #include "factories/digi/PulseNoise_factory.h"
0037 
0038 extern "C" {
0039 void InitPlugin(JApplication* app) {
0040 
0041   using namespace eicrecon;
0042 
0043   InitJANAPlugin(app);
0044 
0045   // Make sure left and right use the same value
0046   decltype(SimCalorimeterHitProcessorConfig::attenuationParameters) EcalBarrelScFi_attPars = {
0047       0.416212, 747.39875 * edm4eic::unit::mm, 7521.88383 * edm4eic::unit::mm};
0048   decltype(SimCalorimeterHitProcessorConfig::hitMergeFields) EcalBarrelScFi_hitMergeFields = {
0049       "fiber", "z"};
0050   decltype(SimCalorimeterHitProcessorConfig::contributionMergeFields)
0051       EcalBarrelScFi_contributionMergeFields = {"fiber"};
0052   decltype(SimCalorimeterHitProcessorConfig::inversePropagationSpeed)
0053       EcalBarrelScFi_inversePropagationSpeed = {(1. / 160) * edm4eic::unit::ns / edm4eic::unit::mm};
0054   decltype(SimCalorimeterHitProcessorConfig::fixedTimeDelay) EcalBarrelScFi_fixedTimeDelay = {
0055       2 * edm4eic::unit::ns};
0056   decltype(SimCalorimeterHitProcessorConfig::timeWindow) EcalBarrelScFi_timeWindow = {
0057       100 * edm4eic::unit::ns};
0058 
0059   decltype(PulseGenerationConfig::pulse_shape_function) EcalBarrelScFi_pulse_shape_function = {
0060       "LandauPulse"};
0061   decltype(PulseGenerationConfig::pulse_shape_params) EcalBarrelScFi_pulse_shape_params = {
0062       1.0, 2 * edm4eic::unit::ns};
0063   decltype(PulseGenerationConfig::ignore_thres) EcalBarrelScFi_ignore_thres = {5.0e-5};
0064   decltype(PulseGenerationConfig::timestep) EcalBarrelScFi_timestep = {0.5 * edm4eic::unit::ns};
0065 
0066   decltype(PulseCombinerConfig::combine_field) EcalBarrelScFi_combine_field           = {"grid"};
0067   decltype(PulseCombinerConfig::minimum_separation) EcalBarrelScFi_minimum_separation = {
0068       100 * edm4eic::unit::ns};
0069   decltype(PulseNoiseConfig::poles) EcalBarrelScFi_poles       = {2};
0070   decltype(PulseNoiseConfig::variance) EcalBarrelScFi_variance = {0.5};
0071   decltype(PulseNoiseConfig::alpha) EcalBarrelScFi_alpha       = {0};
0072   decltype(PulseNoiseConfig::scale) EcalBarrelScFi_scale       = {5.4e-5};
0073   decltype(PulseNoiseConfig::pedestal) EcalBarrelScFi_pedestal = {1.6e-4};
0074 
0075   // Make sure digi and reco use the same value
0076   decltype(CalorimeterHitDigiConfig::capADC) EcalBarrelScFi_capADC = 16384; //16384,  14bit ADC
0077   decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalBarrelScFi_dyRangeADC   = 1500 * dd4hep::MeV;
0078   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalBarrelScFi_pedMeanADC   = 100;
0079   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalBarrelScFi_pedSigmaADC = 1;
0080   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalBarrelScFi_resolutionTDC =
0081       10 * dd4hep::picosecond;
0082   app->Add(new JOmniFactoryGeneratorT<SimCalorimeterHitProcessor_factory>(
0083       "EcalBarrelScFiPAttenuatedHits", {"EcalBarrelScFiHits"},
0084       {"EcalBarrelScFiPAttenuatedHits", "EcalBarrelScFiPAttenuatedHitContributions"},
0085       {
0086           .attenuationParameters            = EcalBarrelScFi_attPars,
0087           .readout                          = "EcalBarrelScFiHits",
0088           .attenuationReferencePositionName = "EcalBarrel_LightGuide_PositivePosZ",
0089           .hitMergeFields                   = EcalBarrelScFi_hitMergeFields,
0090           .contributionMergeFields          = EcalBarrelScFi_contributionMergeFields,
0091           .inversePropagationSpeed          = EcalBarrelScFi_inversePropagationSpeed,
0092           .fixedTimeDelay                   = EcalBarrelScFi_fixedTimeDelay,
0093           .timeWindow                       = EcalBarrelScFi_timeWindow,
0094       },
0095       app // TODO: Remove me once fixed
0096       ));
0097   app->Add(new JOmniFactoryGeneratorT<SimCalorimeterHitProcessor_factory>(
0098       "EcalBarrelScFiNAttenuatedHits", {"EcalBarrelScFiHits"},
0099       {"EcalBarrelScFiNAttenuatedHits", "EcalBarrelScFiNAttenuatedHitContributions"},
0100       {
0101           .attenuationParameters            = EcalBarrelScFi_attPars,
0102           .readout                          = "EcalBarrelScFiHits",
0103           .attenuationReferencePositionName = "EcalBarrel_LightGuide_NegativePosZ",
0104           .hitMergeFields                   = EcalBarrelScFi_hitMergeFields,
0105           .contributionMergeFields          = EcalBarrelScFi_contributionMergeFields,
0106           .inversePropagationSpeed          = EcalBarrelScFi_inversePropagationSpeed,
0107           .fixedTimeDelay                   = EcalBarrelScFi_fixedTimeDelay,
0108           .timeWindow                       = EcalBarrelScFi_timeWindow,
0109       },
0110       app // TODO: Remove me once fixed
0111       ));
0112   app->Add(new JOmniFactoryGeneratorT<PulseGeneration_factory<edm4hep::SimCalorimeterHit>>(
0113       "EcalBarrelScFiPPulses", {"EcalBarrelScFiPAttenuatedHits"}, {"EcalBarrelScFiPPulses"},
0114       {
0115           .pulse_shape_function = EcalBarrelScFi_pulse_shape_function,
0116           .pulse_shape_params   = EcalBarrelScFi_pulse_shape_params,
0117           .ignore_thres         = EcalBarrelScFi_ignore_thres,
0118           .timestep             = EcalBarrelScFi_timestep,
0119       },
0120       app // TODO: Remove me once fixed
0121       ));
0122   app->Add(new JOmniFactoryGeneratorT<PulseGeneration_factory<edm4hep::SimCalorimeterHit>>(
0123       "EcalBarrelScFiNPulses", {"EcalBarrelScFiNAttenuatedHits"}, {"EcalBarrelScFiNPulses"},
0124       {
0125           .pulse_shape_function = EcalBarrelScFi_pulse_shape_function,
0126           .pulse_shape_params   = EcalBarrelScFi_pulse_shape_params,
0127           .ignore_thres         = EcalBarrelScFi_ignore_thres,
0128           .timestep             = EcalBarrelScFi_timestep,
0129       },
0130       app // TODO: Remove me once fixed
0131       ));
0132   app->Add(new JOmniFactoryGeneratorT<PulseCombiner_factory>(
0133       "EcalBarrelScFiPCombinedPulses", {"EcalBarrelScFiPPulses"}, {"EcalBarrelScFiPCombinedPulses"},
0134       {
0135           .minimum_separation = EcalBarrelScFi_minimum_separation,
0136           .readout            = "EcalBarrelScFiHits",
0137           .combine_field      = EcalBarrelScFi_combine_field,
0138       },
0139       app // TODO: Remove me once fixed
0140       ));
0141   app->Add(new JOmniFactoryGeneratorT<PulseCombiner_factory>(
0142       "EcalBarrelScFiNCombinedPulses", {"EcalBarrelScFiNPulses"}, {"EcalBarrelScFiNCombinedPulses"},
0143       {
0144           .minimum_separation = EcalBarrelScFi_minimum_separation,
0145           .readout            = "EcalBarrelScFiHits",
0146           .combine_field      = EcalBarrelScFi_combine_field,
0147       },
0148       app // TODO: Remove me once fixed
0149       ));
0150   app->Add(new JOmniFactoryGeneratorT<PulseNoise_factory>(
0151       "EcalBarrelScFiPCombinedPulsesWithNoise", {"EventHeader", "EcalBarrelScFiPCombinedPulses"},
0152       {"EcalBarrelScFiPCombinedPulsesWithNoise"},
0153       {
0154           .poles    = EcalBarrelScFi_poles,
0155           .variance = EcalBarrelScFi_variance,
0156           .alpha    = EcalBarrelScFi_alpha,
0157           .scale    = EcalBarrelScFi_scale,
0158           .pedestal = EcalBarrelScFi_pedestal,
0159       },
0160       app // TODO: Remove me once fixed
0161       ));
0162   app->Add(new JOmniFactoryGeneratorT<PulseNoise_factory>(
0163       "EcalBarrelScFiNCombinedPulsesWithNoise", {"EventHeader", "EcalBarrelScFiNCombinedPulses"},
0164       {"EcalBarrelScFiNCombinedPulsesWithNoise"},
0165       {
0166           .poles    = EcalBarrelScFi_poles,
0167           .variance = EcalBarrelScFi_variance,
0168           .alpha    = EcalBarrelScFi_alpha,
0169           .scale    = EcalBarrelScFi_scale,
0170           .pedestal = EcalBarrelScFi_pedestal,
0171       },
0172       app // TODO: Remove me once fixed
0173       ));
0174   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0175       "EcalBarrelScFiRawHits", {"EventHeader", "EcalBarrelScFiHits"},
0176       {"EcalBarrelScFiRawHits", "EcalBarrelScFiRawHitAssociations"},
0177       {
0178           .eRes          = {0.0 * sqrt(dd4hep::GeV), 0.0, 0.0 * dd4hep::GeV},
0179           .tRes          = 0.0 * dd4hep::ns,
0180           .threshold     = 0.0 * dd4hep::keV, // threshold is set in ADC in reco
0181           .capADC        = EcalBarrelScFi_capADC,
0182           .dyRangeADC    = EcalBarrelScFi_dyRangeADC,
0183           .pedMeanADC    = EcalBarrelScFi_pedMeanADC,
0184           .pedSigmaADC   = EcalBarrelScFi_pedSigmaADC,
0185           .resolutionTDC = EcalBarrelScFi_resolutionTDC,
0186           .corrMeanScale = "1.0",
0187           .readout       = "EcalBarrelScFiHits",
0188           .fields        = {"fiber", "z"},
0189       },
0190       app // TODO: Remove me once fixed
0191       ));
0192   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0193       "EcalBarrelScFiRecHits", {"EcalBarrelScFiRawHits"}, {"EcalBarrelScFiRecHits"},
0194       {
0195           .capADC          = EcalBarrelScFi_capADC,
0196           .dyRangeADC      = EcalBarrelScFi_dyRangeADC,
0197           .pedMeanADC      = EcalBarrelScFi_pedMeanADC,
0198           .pedSigmaADC     = EcalBarrelScFi_pedSigmaADC, // not needed; use only thresholdValue
0199           .resolutionTDC   = EcalBarrelScFi_resolutionTDC,
0200           .thresholdFactor = 0.0, // use only thresholdValue
0201           .thresholdValue  = 5.0, // 16384 ADC counts/1500 MeV * 0.5 MeV (desired threshold) = 5.46
0202           .sampFrac        = "0.09285755",
0203           .readout         = "EcalBarrelScFiHits",
0204           .layerField      = "layer",
0205           .sectorField     = "sector",
0206           .localDetFields  = {"system", "sector"},
0207           // here we want to use grid center position (XY) but keeps the z information from fiber-segment
0208           // TODO: a more realistic way to get z is to reconstruct it from timing
0209           .maskPos       = "xy",
0210           .maskPosFields = {"fiber", "z"},
0211       },
0212       app // TODO: Remove me once fixed
0213       ));
0214   app->Add(new JOmniFactoryGeneratorT<CalorimeterIslandCluster_factory>(
0215       "EcalBarrelScFiProtoClusters", {"EcalBarrelScFiRecHits"}, {"EcalBarrelScFiProtoClusters"},
0216       {
0217           .adjacencyMatrix{},
0218           .peakNeighbourhoodMatrix{},
0219           .readout{},
0220           .sectorDist = 50. * dd4hep::mm,
0221           .localDistXY{},
0222           .localDistXZ = {80 * dd4hep::mm, 80 * dd4hep::mm},
0223           .localDistYZ{},
0224           .globalDistRPhi{},
0225           .globalDistEtaPhi{},
0226           .dimScaledLocalDistXY{},
0227           .splitCluster         = false,
0228           .minClusterHitEdep    = 5.0 * dd4hep::MeV,
0229           .minClusterCenterEdep = 100.0 * dd4hep::MeV,
0230           .transverseEnergyProfileMetric{},
0231           .transverseEnergyProfileScale{},
0232           .transverseEnergyProfileScaleUnits{},
0233       },
0234       app // TODO: Remove me once fixed
0235       ));
0236   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
0237       "EcalBarrelScFiClustersWithoutShapes",
0238       {"EcalBarrelScFiProtoClusters",         // edm4eic::ProtoClusterCollection
0239        "EcalBarrelScFiRawHitAssociations"},   // edm4eic::MCRecoCalorimeterHitAssociation
0240       {"EcalBarrelScFiClustersWithoutShapes", // edm4eic::Cluster
0241        "EcalBarrelScFiClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
0242       {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, .enableEtaBounds = false},
0243       app // TODO: Remove me once fixed
0244       ));
0245   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
0246       "EcalBarrelScFiClusters",
0247       {"EcalBarrelScFiClustersWithoutShapes", "EcalBarrelScFiClusterAssociationsWithoutShapes"},
0248       {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations"},
0249       {.longitudinalShowerInfoAvailable = true, .energyWeight = "log", .logWeightBase = 6.2}, app));
0250 
0251   // Make sure digi and reco use the same value
0252   decltype(CalorimeterHitDigiConfig::capADC) EcalBarrelImaging_capADC = 8192; //8192,  13bit ADC
0253   decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalBarrelImaging_dyRangeADC = 3 * dd4hep::MeV;
0254   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalBarrelImaging_pedMeanADC =
0255       14; // Noise floor at 5 keV: 8192 / 3 * 0.005
0256   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalBarrelImaging_pedSigmaADC =
0257       5; // Upper limit for sigma for AstroPix
0258   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalBarrelImaging_resolutionTDC =
0259       3.25 * dd4hep::nanosecond;
0260   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
0261       "EcalBarrelImagingRawHits", {"EventHeader", "EcalBarrelImagingHits"},
0262       {"EcalBarrelImagingRawHits", "EcalBarrelImagingRawHitAssociations"},
0263       {
0264           .eRes          = {0.0 * sqrt(dd4hep::GeV), 0.02, 0.0 * dd4hep::GeV},
0265           .tRes          = 0.0 * dd4hep::ns,
0266           .capADC        = EcalBarrelImaging_capADC,
0267           .dyRangeADC    = EcalBarrelImaging_dyRangeADC,
0268           .pedMeanADC    = EcalBarrelImaging_pedMeanADC,
0269           .pedSigmaADC   = EcalBarrelImaging_pedSigmaADC,
0270           .resolutionTDC = EcalBarrelImaging_resolutionTDC,
0271           .corrMeanScale = "1.0",
0272           .readout       = "EcalBarrelImagingHits",
0273       },
0274       app // TODO: Remove me once fixed
0275       ));
0276   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
0277       "EcalBarrelImagingRecHits", {"EcalBarrelImagingRawHits"}, {"EcalBarrelImagingRecHits"},
0278       {
0279           .capADC          = EcalBarrelImaging_capADC,
0280           .dyRangeADC      = EcalBarrelImaging_dyRangeADC,
0281           .pedMeanADC      = EcalBarrelImaging_pedMeanADC,
0282           .pedSigmaADC     = EcalBarrelImaging_pedSigmaADC, // not needed; use only thresholdValue
0283           .resolutionTDC   = EcalBarrelImaging_resolutionTDC,
0284           .thresholdFactor = 0.0, // use only thresholdValue
0285           .thresholdValue  = 41,  // 8192 ADC counts/3 MeV * 0.015 MeV (desired threshold) = 41
0286           .sampFrac        = "0.00429453",
0287           .readout         = "EcalBarrelImagingHits",
0288           .layerField      = "layer",
0289           .sectorField     = "sector",
0290       },
0291       app // TODO: Remove me once fixed
0292       ));
0293   app->Add(new JOmniFactoryGeneratorT<ImagingTopoCluster_factory>(
0294       "EcalBarrelImagingProtoClusters", {"EcalBarrelImagingRecHits"},
0295       {"EcalBarrelImagingProtoClusters"},
0296       {
0297           .neighbourLayersRange = 2, //  # id diff for adjacent layer
0298           .sameLayerDistTZ      = {2.0 * dd4hep::mm, 2 * dd4hep::mm},     //  # same layer
0299           .diffLayerDistEtaPhi  = {10 * dd4hep::mrad, 10 * dd4hep::mrad}, //  # adjacent layer
0300           .sameLayerMode        = eicrecon::ImagingTopoClusterConfig::ELayerMode::tz,
0301           .diffLayerMode        = eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi,
0302           .sectorDist           = 3.0 * dd4hep::cm,
0303           .minClusterHitEdep    = 0,
0304           .minClusterCenterEdep = 0,
0305           .minClusterEdep       = 100 * dd4hep::MeV,
0306           .minClusterNhits      = 10,
0307       },
0308       app // TODO: Remove me once fixed
0309       ));
0310 
0311   app->Add(new JOmniFactoryGeneratorT<ImagingClusterReco_factory>(
0312       "EcalBarrelImagingClustersWithoutShapes",
0313       {"EcalBarrelImagingProtoClusters", "EcalBarrelImagingRawHitAssociations"},
0314       {"EcalBarrelImagingClustersWithoutShapes",
0315        "EcalBarrelImagingClusterAssociationsWithoutShapes", "EcalBarrelImagingLayers"},
0316       {
0317           .trackStopLayer = 6,
0318       },
0319       app // TODO: Remove me once fixed
0320       ));
0321   app->Add(new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
0322       "EcalBarrelImagingClusters",
0323       {"EcalBarrelImagingClustersWithoutShapes",
0324        "EcalBarrelImagingClusterAssociationsWithoutShapes"},
0325       {"EcalBarrelImagingClusters", "EcalBarrelImagingClusterAssociations"},
0326       {.longitudinalShowerInfoAvailable = false, .energyWeight = "log", .logWeightBase = 6.2},
0327       app));
0328   app->Add(new JOmniFactoryGeneratorT<EnergyPositionClusterMerger_factory>(
0329       "EcalBarrelClusters",
0330       {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations", "EcalBarrelImagingClusters",
0331        "EcalBarrelImagingClusterAssociations"},
0332       {"EcalBarrelClusters", "EcalBarrelClusterAssociations"},
0333       {
0334           .energyRelTolerance = 0.5,
0335           .phiTolerance       = 0.1,
0336           .etaTolerance       = 0.2,
0337       },
0338       app // TODO: Remove me once fixed
0339       ));
0340   app->Add(new JOmniFactoryGeneratorT<TruthEnergyPositionClusterMerger_factory>(
0341       "EcalBarrelTruthClusters",
0342       {"MCParticles", "EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations",
0343        "EcalBarrelImagingClusters", "EcalBarrelImagingClusterAssociations"},
0344       {"EcalBarrelTruthClusters", "EcalBarrelTruthClusterAssociations"},
0345       app // TODO: Remove me once fixed
0346       ));
0347 }
0348 }