detectors/RICH-IRT/RICH-IRT.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024, 2025, 2026, Alexander Kiselev
0003
0004 #if IRT2_VERSION_MAJOR > 2 || (IRT2_VERSION_MAJOR == 2 && IRT2_VERSION_MINOR >= 2)
0005
0006 #include <DD4hep/DetElement.h>
0007 #include <DD4hep/Detector.h>
0008 #include <DD4hep/Volumes.h>
0009 #include <Evaluator/DD4hepUnits.h>
0010 #include <IRT2/CherenkovDetector.h>
0011 #include <IRT2/CherenkovDetectorCollection.h>
0012 #include <IRT2/CherenkovRadiator.h>
0013 #include <IRT2/ParametricSurface.h>
0014 #include <JANA/JApplication.h>
0015 #include <JANA/JApplicationFwd.h>
0016 #include <JANA/JException.h>
0017 #include <JANA/Utils/JTypeInfo.h>
0018 #include <TGeoVolume.h>
0019 #include <TString.h>
0020 #include <TVector3.h>
0021 #include <cctype>
0022 #include <edm4eic/TrackPoint.h>
0023 #include <nlohmann/json.hpp>
0024 #include <nlohmann/json_fwd.hpp>
0025 #include <string.h>
0026 #include <algorithm>
0027 #include <cmath>
0028 #include <format>
0029 #include <fstream>
0030 #include <functional>
0031 #include <gsl/pointers>
0032 #include <iterator>
0033 #include <map>
0034 #include <memory>
0035 #include <string>
0036 #include <utility>
0037 #include <variant>
0038 #include <vector>
0039
0040 #include "algorithms/pid/IrtInterfaceConfig.h"
0041 // ACTS projection algorithm configuration;
0042 #include "algorithms/tracking/TrackPropagationConfig.h"
0043 // JANA;
0044 #include "extensions/jana/JOmniFactoryGeneratorT.h"
0045 // Factories;
0046 #include "factories/tracking/TrackPropagation_factory.h"
0047 #include "global/pid/IrtInterface_factory.h"
0048 // DD4HEP geometry services;
0049 #include "services/geometry/dd4hep/DD4hep_service.h"
0050
0051 using json = nlohmann::json;
0052
0053 // Have to start hardcoding known RICH detectors in some way, before trying
0054 // to pull their configuration data out; for the time being the list is:
0055 //
0056 // PFRICH: ePIC backward proximity focusing RICH
0057 //  DRICH: ePIC forward dual radiator RICH
0058 //
0059 static const std::vector<std::string> RICHes = {"DRICH", "PFRICH"};
0060
0061 using namespace eicrecon;
0062
0063 extern "C" {
0064 void InitPlugin(JApplication* app) {
0065   InitJANAPlugin(app);
0066
0067   auto dd4hep_service = app->GetService<DD4hep_service>();
0068   auto det            = dd4hep_service->detector();
0069
0070   // Loop through all known RICH detectors handled by IRT 2.0 algorithm through this plugin;
0071   for (const auto& RICH : RICHes) {
0072     // First a sanity cross-check: detector with this name should be present in the geometry;
0073     {
0074       bool exists = false;
0075
0076       // FIXME: is there a way to poll a detector existence without a segfault?;
0077       // For now resort to using a loop; assume RICH detector category is "tracker"
0078       // (it presently is); FIXME: may want to loop through all types?;
0079       for (const auto& detElem : det->detectors("tracker"))
0080         if (RICH == detElem.volume()->GetName()) {
0081           exists = true;
0082           break;
0083         } //for detElem
0084
0085       // Detector was not present in the geometry -> throw error
0086       if (!exists)
0087         continue;
0088     }
0089
0090     // Check command line options; expect a "config" key with a JSON configuration file name;
0091     {
0092       std::string rich_lowercase = RICH;
0093       transform(rich_lowercase.begin(), rich_lowercase.end(), rich_lowercase.begin(), ::tolower);
0094       static std::map<std::string, std::string> config_files;
0095       std::string& config_file = config_files[RICH];
0096       config_file              = std::format("calibrations/irt2/{}-reco_v2.json", rich_lowercase);
0097       app->SetDefaultParameter(std::format("{}:config", RICH), config_file,
0098                                std::format("Path to config file for {}", RICH));
0099
0100       if (config_file == "") {
0101         throw JException("RICH detector '%s' requires a config file", RICH.c_str());
0102       }
0103
0104       IrtConfig config;
0105       config.m_detector_name         = RICH;
0106       config.m_json_config_file_name = config_file;
0107
0108       // Import JSON configuration file; sanity checks for several keys which are supposed
0109       // to be present; FIXME: add warning / error printouts;
0110       {
0111         std::ifstream fcfg(config_file);
0112         if (!fcfg)
0113           throw JException("RICH detector '%s' cannot open config file '%s'", RICH.c_str(),
0114                            config_file.c_str());
0115         auto json_config = json::parse(fcfg);
0116         // For less typing;
0117         json* jptr = &json_config;
0118
0119         // An entry describing a nominal acceptance should be present;
0120         if (jptr->find("Acceptance") == jptr->end())
0121           throw JException("RICH detector '%s' config missing 'Acceptance' section", RICH.c_str());
0122         {
0123           const auto& aconfig = (*jptr)["Acceptance"];
0124
0125           if (aconfig.find("eta-min") == aconfig.end() || aconfig.find("eta-max") == aconfig.end())
0126             throw JException(
0127                 "RICH detector '%s' config missing 'eta-min' or 'eta-max' in Acceptance section",
0128                 RICH.c_str());
0129
0130           config.m_eta_min = aconfig["eta-min"].template get<double>();
0131           config.m_eta_max = aconfig["eta-max"].template get<double>();
0132           // FIXME: it will for now (and perhaps forever) be assumed that these min
0133           // and max values have the same sign (and are therefore usable to identify
0134           // what endcap the detector belongs to);
0135           if (config.m_eta_min * config.m_eta_max < 0.0)
0136             throw JException(
0137                 "RICH detector '%s' eta-min and eta-max must have same sign, got %f and %f",
0138                 RICH.c_str(), config.m_eta_min, config.m_eta_max);
0139         }
0140
0141         // And a group of entries describing various radiator parameters;
0142         if (jptr->find("Radiators") == jptr->end())
0143           throw JException("RICH detector '%s' config missing 'Radiators' section", RICH.c_str());
0144         const auto& rconfig = (*jptr)["Radiators"];
0145
0146         // Obtain a handle to a Cherenkov detector optics configuration;
0147         {
0148           auto irt_geometry = IRT2::CherenkovDetectorCollection::Instance();
0149           if (!irt_geometry)
0150             throw JException(
0151                 "RICH detector '%s' failed to get CherenkovDetectorCollection instance",
0152                 RICH.c_str());
0153
0154           auto cdet = irt_geometry->GetDetector(RICH.c_str());
0155           if (!cdet)
0156             throw JException("RICH detector '%s' not found in IRT geometry collection",
0157                              RICH.c_str());
0158
0159           //
0160           // Everything is fine, proceed with the essential part;
0161           //
0162
0163           // Track propagation config;
0164           TrackPropagationConfig track_cfg;
0165
0166           // Eventually do all this in a proper way: consider all sectors, and make use
0167           // of eta boundaries to define a Z-range for a particular radiator; reorder radiators in Z if needed;
0168           std::map<double, std::pair<double, double>> disks;
0169
0170           for (auto [name, radiator] : cdet->Radiators()) {
0171             // There should be an entry in JSON file; skip otherwise;
0172             if (rconfig.find(name.Data()) == rconfig.end())
0173               continue;
0174
0175             const auto& rrconfig = rconfig[name.Data()];
0176             // FIXME: per default assume that all defined radiators participate in ring imaging?;
0177             if (rrconfig.find("imaging") != rrconfig.end() &&
0178                 !strcmp(rrconfig["imaging"].template get<std::string>().c_str(), "no"))
0179               continue;
0180
0181             // Well, this option which instructs ACTS where to build track projections must be present;
0182             if (rrconfig.find("acts-planes") == rrconfig.end())
0183               continue;
0184
0185             unsigned numPlanes = rrconfig["acts-planes"].template get<int>();
0186             if (!numPlanes)
0187               throw JException("RICH detector '%s' radiator '%s' has zero acts-planes",
0188                                RICH.c_str(), name.Data());
0189
0190             double theta_min = 2 * std::atan(exp(-config.m_eta_min));
0191             double theta_max = 2 * std::atan(exp(-config.m_eta_max));
0192
0193             // Estimate a required Z-range;
0194             double zmin = 0.0, zmax = 0.0;
0195             for (unsigned iq = 0; iq < 2; iq++) {
0196               double theta = iq ? theta_max : theta_min;
0197               TVector3 x0(0, 0, 0), n0(0.0, sin(theta), cos(theta)), from, to;
0198               unsigned isec = cdet->GetSector(n0);
0199
0200               auto sf = radiator->GetFrontSide(isec);
0201               auto sr = radiator->GetRearSide(isec);
0202
0203               // FIXME: may want to check return codes?;
0204               [[maybe_unused]] bool bf = sf->GetCrossing(x0, n0, &from, false);
0205               [[maybe_unused]] bool br = sr->GetCrossing(x0, n0, &to, false);
0206
0207               double zf = from.Z(), zr = to.Z();
0208               if (!iq || fabs(zf) < fabs(zmin))
0209                 zmin = zf;
0210               if (!iq || fabs(zr) > fabs(zmax))
0211                 zmax = zr;
0212             } //for iq
0213
0214             // "+1": avoid a "coordinate at the boundary condition"; essentially make 'numPlanes'
0215             // bins and use bin centers;
0216             double zmid = (zmin + zmax) / 2 * dd4hep::mm,
0217                    step = fabs(zmax - zmin) / (numPlanes + 1) * dd4hep::mm;
0218
0219             for (unsigned i = 0; i < numPlanes; i++) {
0220               auto zCoord =
0221                   zmid + step * (static_cast<int>(i) - (static_cast<int>(numPlanes) - 1) / 2.);
0222               double rmin = fabs(zCoord * tan(theta_min)), rmax = fabs(zCoord * tan(theta_max));
0223               if (rmax < rmin)
0224                 std::swap(rmin, rmax);
0225
0226               // Yes, prefer to order in ascending fabs(z) order in both endcaps; FIXME: implicitly assume
0227               // all these coordinates are >0 in the hadron-going endcap and <0 in the electron-going one;
0228               disks[fabs(zCoord)] = std::make_pair(rmin, rmax);
0229             } //for i
0230           } //for radiators
0231
0232           const char* tag = config.m_eta_min > 0.0 ? "ForwardRICH_ID" : "BackwardRICH_ID";
0233           // Will be ordered along the track direction (assuming IP is the origin),
0234           // no matter ACTS performs reordering internally or not;
0235           for (auto disk : disks) {
0236             double z    = (config.m_eta_min > 0.0 ? 1.0 : -1.0) * disk.first;
0237             double rmin = disk.second.first, rmax = disk.second.second;
0238
0239             track_cfg.target_surfaces.push_back(eicrecon::DiscSurfaceConfig{tag, z, rmin, rmax});
0240           } //for disk
0241
0242           // FIXME: this should be fine?;
0243           track_cfg.filter_surfaces.push_back(track_cfg.target_surfaces.back());
0244
0245           // FIXME: may not be a good idea for dRICH;
0246           track_cfg.track_point_cut =
0247               std::function<bool(edm4eic::TrackPoint)>([](edm4eic::TrackPoint) { return true; });
0248
0249           // Eventually define the factories and collections to be used;
0250           {
0251             TString RICHstr(RICH), RICHtracks = RICHstr + "Tracks";
0252
0253             // Charged particle track projections;
0254             app->Add(new JOmniFactoryGeneratorT<TrackPropagation_factory>(
0255                 RICHtracks.Data(),
0256                 {"CentralCKFTracks", "CentralCKFActsTrackStates", "CentralCKFActsTracks"},
0257                 {RICHtracks.Data()}, track_cfg, app));
0258
0259             // A unified IRT 2.1 algorithm; FIXME: split digitization step off later;
0260             app->Add(new JOmniFactoryGeneratorT<IrtInterface_factory>(
0261                 (RICHstr + "IrtInterface").Data(),
0262                 {"EventHeader", "MCParticles", "CentralCKFTracks", "CentralCKFTrackLinks",
0263                  RICHtracks.Data(), (RICHstr + "Hits").Data()},
0264                 {(RICHstr + "IrtRadiatorInfo").Data(), (RICHstr + "IrtParticles").Data()}, config,
0265                 app));
0266           }
0267         }
0268       }
0269     }
0270   } //for RICH
0271 }
0272 }
0273
0274 #endif