detectors/DRICH/DRICH.cc

0001 // Copyright (C) 2022, 2023, Christopher Dilks, Luigi Dello Stritto
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003
0004 // SPDX-License-Identifier: LGPL-3.0-or-later
0005 // Copyright (C) 2024, Dmitry Kalinkin
0006
0007 #include <DD4hep/DetElement.h>
0008 #include <DD4hep/Detector.h>
0009 #include <Evaluator/DD4hepUnits.h>
0010 #include <JANA/JApplication.h>
0011 #include <math.h>
0012 #include <functional>
0013 #include <gsl/pointers>
0014 #include <map>
0015 #include <memory>
0016 #include <stdexcept>
0017 #include <string>
0018 #include <utility>
0019 #include <vector>
0020
0021 // algorithm configurations
0022 #include "algorithms/digi/PhotoMultiplierHitDigiConfig.h"
0023 #include "algorithms/pid/IrtCherenkovParticleIDConfig.h"
0024 #include "algorithms/pid/MergeParticleIDConfig.h"
0025 #include "algorithms/pid_lut/PIDLookupConfig.h"
0026 #include "algorithms/tracking/TrackPropagationConfig.h"
0027 #include "extensions/jana/JOmniFactoryGeneratorT.h"
0028 // factories
0029 #include "global/digi/PhotoMultiplierHitDigi_factory.h"
0030 #include "global/pid/IrtCherenkovParticleID_factory.h"
0031 #include "global/pid/MergeCherenkovParticleID_factory.h"
0032 #include "global/pid/MergeTrack_factory.h"
0033 #include "global/pid/RichTrack_factory.h"
0034 #include "global/pid_lut/PIDLookup_factory.h"
0035 #include "services/geometry/dd4hep/DD4hep_service.h"
0036 #include "services/geometry/richgeo/ActsGeo.h"
0037 #include "services/geometry/richgeo/RichGeo.h"
0038 #include "services/geometry/richgeo/RichGeo_service.h"
0039
0040 extern "C" {
0041   void InitPlugin(JApplication *app) {
0042     InitJANAPlugin(app);
0043
0044     using namespace eicrecon;
0045
0046     // configuration parameters ///////////////////////////////////////////////
0047
0048     // digitization
0049     PhotoMultiplierHitDigiConfig digi_cfg;
0050     digi_cfg.seed            = 5; // FIXME: set to 0 for a 'unique' seed, but
0051                                   // that seems to delay the RNG from actually randomizing
0052     digi_cfg.hitTimeWindow   = 20.0; // [ns]
0053     digi_cfg.timeResolution  = 1/16.0; // [ns]
0054     digi_cfg.speMean         = 80.0;
0055     digi_cfg.speError        = 16.0;
0056     digi_cfg.pedMean         = 200.0;
0057     digi_cfg.pedError        = 3.0;
0058     digi_cfg.enablePixelGaps = true;
0059     digi_cfg.safetyFactor    = 0.7;
0060     digi_cfg.enableNoise     = false;
0061     digi_cfg.noiseRate       = 20000; // [Hz]
0062     digi_cfg.noiseTimeWindow = 20.0 * dd4hep::ns; // [ns]
0063     digi_cfg.quantumEfficiency.clear();
0064     digi_cfg.quantumEfficiency = { // wavelength units are [nm]
0065       {315,  0.00},
0066       {325,  0.04},
0067       {340,  0.10},
0068       {350,  0.20},
0069       {370,  0.30},
0070       {400,  0.35},
0071       {450,  0.40},
0072       {500,  0.38},
0073       {550,  0.35},
0074       {600,  0.27},
0075       {650,  0.20},
0076       {700,  0.15},
0077       {750,  0.12},
0078       {800,  0.08},
0079       {850,  0.06},
0080       {900,  0.04},
0081       {1000, 0.00}
0082     };
0083
0084     // Track propagation
0085     TrackPropagationConfig aerogel_track_cfg;
0086     TrackPropagationConfig gas_track_cfg;
0087
0088     // get RICH geo service
0089     auto richGeoSvc = app->GetService<RichGeo_service>();
0090     auto dd4hepGeo = richGeoSvc->GetDD4hepGeo();
0091     if (dd4hepGeo->world().children().contains("DRICH")) {
0092       auto actsGeo = richGeoSvc->GetActsGeo("DRICH");
0093       auto aerogel_tracking_planes = actsGeo->TrackingPlanes(richgeo::kAerogel, 5);
0094       auto aerogel_track_point_cut = actsGeo->TrackPointCut(richgeo::kAerogel);
0095       auto gas_tracking_planes = actsGeo->TrackingPlanes(richgeo::kGas, 10);
0096       auto gas_track_point_cut = actsGeo->TrackPointCut(richgeo::kGas);
0097       auto filter_surface = gas_tracking_planes.back();
0098       // track propagation to each radiator
0099       aerogel_track_cfg.filter_surfaces.push_back(filter_surface);
0100       aerogel_track_cfg.target_surfaces = aerogel_tracking_planes;
0101       aerogel_track_cfg.track_point_cut = aerogel_track_point_cut;
0102       gas_track_cfg.filter_surfaces.push_back(filter_surface);
0103       gas_track_cfg.target_surfaces = gas_tracking_planes;
0104       gas_track_cfg.track_point_cut = gas_track_point_cut;
0105     }
0106
0107     // IRT PID
0108     IrtCherenkovParticleIDConfig irt_cfg;
0109     // - refractive index interpolation
0110     irt_cfg.numRIndexBins = 100;
0111     // - aerogel
0112     irt_cfg.radiators.insert({"Aerogel", RadiatorConfig{}});
0113     irt_cfg.radiators.at("Aerogel").referenceRIndex = 1.0190;
0114     irt_cfg.radiators.at("Aerogel").attenuation     = 48; // [mm]
0115     irt_cfg.radiators.at("Aerogel").smearingMode    = "gaussian";
0116     irt_cfg.radiators.at("Aerogel").smearing        = 2e-3; // [radians]
0117     // - gas
0118     irt_cfg.radiators.insert({"Gas", RadiatorConfig{}});
0119     irt_cfg.radiators.at("Gas").referenceRIndex = 1.00076;
0120     irt_cfg.radiators.at("Gas").attenuation     = 0; // [mm]
0121     irt_cfg.radiators.at("Gas").smearingMode    = "gaussian";
0122     irt_cfg.radiators.at("Gas").smearing        = 5e-3; // [radians]
0123     // - PDG list
0124     irt_cfg.pdgList.insert(irt_cfg.pdgList.end(), { 11, 211, 321, 2212 });
0125     // - cheat modes
0126     irt_cfg.cheatPhotonVertex  = false;
0127     irt_cfg.cheatTrueRadiator  = false;
0128
0129     // Merge PID from radiators
0130     MergeParticleIDConfig merge_cfg;
0131     merge_cfg.mergeMode = MergeParticleIDConfig::kAddWeights;
0132
0133     // wiring between factories and data ///////////////////////////////////////
0134     // clang-format off
0135
0136     // digitization
0137     app->Add(new JOmniFactoryGeneratorT<PhotoMultiplierHitDigi_factory>(
0138           "DRICHRawHits",
0139           {"DRICHHits"},
0140           {"DRICHRawHits", "DRICHRawHitsAssociations"},
0141           digi_cfg,
0142           app
0143           ));
0144
0145     // charged particle tracks
0146     app->Add(new JOmniFactoryGeneratorT<RichTrack_factory>(
0147           "DRICHAerogelTracks",
0148           {"CentralCKFTracks", "CentralCKFActsTrajectories", "CentralCKFActsTracks"},
0149           {"DRICHAerogelTracks"},
0150           aerogel_track_cfg,
0151           app
0152           ));
0153     app->Add(new JOmniFactoryGeneratorT<RichTrack_factory>(
0154           "DRICHGasTracks",
0155           {"CentralCKFTracks", "CentralCKFActsTrajectories", "CentralCKFActsTracks"},
0156           {"DRICHGasTracks"},
0157           gas_track_cfg,
0158           app
0159           ));
0160
0161     app->Add(new JOmniFactoryGeneratorT<MergeTrack_factory>(
0162           "DRICHMergedTracks",
0163           {"DRICHAerogelTracks", "DRICHGasTracks"},
0164           {"DRICHMergedTracks"},
0165           {},
0166           app
0167           ));
0168
0169     // PID algorithm
0170     app->Add(new JOmniFactoryGeneratorT<IrtCherenkovParticleID_factory>(
0171           "DRICHIrtCherenkovParticleID",
0172           {
0173             "DRICHAerogelTracks", "DRICHGasTracks", "DRICHMergedTracks",
0174             "DRICHRawHits",
0175             "DRICHRawHitsAssociations"
0176           },
0177           {"DRICHAerogelIrtCherenkovParticleID", "DRICHGasIrtCherenkovParticleID"},
0178           irt_cfg,
0179           app
0180           ));
0181
0182     // merge aerogel and gas PID results
0183     app->Add(new JOmniFactoryGeneratorT<MergeCherenkovParticleID_factory>(
0184           "DRICHMergedIrtCherenkovParticleID",
0185           {"DRICHAerogelIrtCherenkovParticleID", "DRICHGasIrtCherenkovParticleID"},
0186           {"DRICHMergedIrtCherenkovParticleID"},
0187           merge_cfg,
0188           app
0189           ));
0190
0191     int ForwardRICH_ID = 0;
0192     try {
0193         auto detector = app->GetService<DD4hep_service>()->detector();
0194         ForwardRICH_ID = detector->constant<int>("ForwardRICH_ID");
0195     } catch(const std::runtime_error&) {
0196         // Nothing
0197     }
0198     PIDLookupConfig pid_cfg {
0199       .filename="calibrations/drich.lut",
0200       .system=ForwardRICH_ID,
0201       .pdg_values={211, 321, 2212},
0202       .charge_values={1},
0203       .momentum_edges={0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75, 10.25, 10.75, 11.25, 11.75, 12.25, 12.75, 13.25, 13.75, 14.25, 14.75, 15.25, 15.75, 16.25, 16.75, 17.25, 17.75, 18.25, 18.75, 19.25, 19.75, 20.50, 21.50, 22.50, 23.50, 24.50, 25.50, 26.50, 27.50, 28.50, 29.50, 30.50},
0204       .polar_edges={0.060, 0.164, 0.269, 0.439},
0205       .azimuthal_binning={0., 2 * M_PI, 2 * M_PI}, // lower, upper, step
0206       .polar_bin_centers_in_lut=true,
0207       .use_radians=true,
0208       .missing_electron_prob=true,
0209     };
0210
0211     app->Add(new JOmniFactoryGeneratorT<PIDLookup_factory>(
0212           "DRICHTruthSeededLUTPID",
0213           {
0214           "ReconstructedTruthSeededChargedWithPFRICHTOFDIRCPIDParticles",
0215           "ReconstructedTruthSeededChargedWithPFRICHTOFDIRCPIDParticleAssociations",
0216           },
0217           {
0218           "ReconstructedTruthSeededChargedParticles",
0219           "ReconstructedTruthSeededChargedParticleAssociations",
0220           "DRICHTruthSeededParticleIDs",
0221           },
0222           pid_cfg,
0223           app
0224           ));
0225
0226     app->Add(new JOmniFactoryGeneratorT<PIDLookup_factory>(
0227           "DRICHLUTPID",
0228           {
0229           "ReconstructedChargedWithPFRICHTOFDIRCPIDParticles",
0230           "ReconstructedChargedWithPFRICHTOFDIRCPIDParticleAssociations",
0231           },
0232           {
0233           "ReconstructedChargedParticles",
0234           "ReconstructedChargedParticleAssociations",
0235           "DRICHParticleIDs",
0236           },
0237           pid_cfg,
0238           app
0239           ));
0240     // clang-format on
0241   }
0242 }