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