Irt 2.1c branch update

CMakeLists.txt

@@ -12,6 +12,7 @@ cmake_policy(SET CMP0079 NEW) # allow `target_link_libraries` from any dir
 option(EVALUATION  "Build evaluation programs"                          OFF)
 option(DELPHES     "Delphes card production"                            ON)
 option(IRT_ROOT_IO "Generate dictionary for ROOT I/O of libIRT objects" ON)
+option(JSON_EXPORT "Be able to export calibrations in JSON format"      OFF)
 
 # Sanitizer options (TSAN is mutually exclusive with ASAN/UBSAN)
 include(CMakeDependentOption)
@@ -84,9 +85,6 @@ include_directories(
   # For now assume that newly installed edm4eic::(CherenkovPID, ...) event structures are available 
   # in the same place where IRT is supposed to be installed;
   ${CMAKE_INSTALL_PREFIX}/include
-
-  # Help ePIC installation find json.hpp if needed;
-  #/opt/local/include
 )
 file(GLOB HEADERS ${PROJECT_SOURCE_DIR}/include/*.h)
 list(FILTER HEADERS EXCLUDE REGEX "LinkDef\\.h$")
@@ -97,6 +95,9 @@ list(FILTER HEADERS EXCLUDE REGEX "LinkDef\\.h$")
 
 # sources
 set(IRT_SRC
+  ${PROJECT_SOURCE_DIR}/source/CherenkovDetector.cc
+  ${PROJECT_SOURCE_DIR}/source/CherenkovDetectorCollection.cc
+
   ${PROJECT_SOURCE_DIR}/source/ParametricSurface.cc
   ${PROJECT_SOURCE_DIR}/source/SphericalSurface.cc
   ${PROJECT_SOURCE_DIR}/source/ToricSurface.cc
@@ -127,6 +128,12 @@ if(NOT IRT_ROOT_IO)
   message(STATUS "NOTE: disabling ROOT dictionary generation")
   target_compile_definitions(${CMAKE_PROJECT_NAME} PUBLIC DISABLE_ROOT_IO)
 endif()
+IF(JSON_EXPORT)
+  message(STATUS "NOTE: JSON_EXPORT enabled")
+  set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DJSON_EXPORT")
+  find_package(nlohmann_json 3.12.0 REQUIRED)
+  target_link_libraries(${CMAKE_PROJECT_NAME} PRIVATE nlohmann_json::nlohmann_json)
+ENDIF()
 target_include_directories(${CMAKE_PROJECT_NAME} PUBLIC
   $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
   $<INSTALL_INTERFACE:include>

include/CherenkovDetector.h

@@ -134,6 +134,8 @@ class CherenkovDetector: public TObject {
   // readout ID -> pixel position converter (for external usage)
   std::function<TVector3(long long int)> m_ReadoutIDToPosition; //!
 
+  void ExportJsonFormatCalibrations(const char *fname);
+
  private:  
   TString m_Name;
   // This is needed for dd4hep cell index decoding;

include/CherenkovDetectorCollection.h

@@ -20,19 +20,21 @@ namespace IRT2 {
 
 class CherenkovDetectorCollection: public BitMask {
  public:
-  CherenkovDetectorCollection() {};
+  CherenkovDetectorCollection();
   // FIXME: populate the dtor, for completeness;
   ~CherenkovDetectorCollection() {};
 
+  static CherenkovDetectorCollection *m_Instance;              //!
+  static CherenkovDetectorCollection *Instance() {
+    return m_Instance ? m_Instance : new CherenkovDetectorCollection();
+  };
+
   CherenkovDetector *AddNewDetector(const char *name) {
     auto det = new CherenkovDetector(name);
     m_Detectors[det->GetName()] = det;
 
     return det;
   };
-  //void AddNewDetector(CherenkovDetector *det) {
-  //m_Detectors[det->GetName()] = det;
-  //};
 
   CherenkovRadiator *FindOrAddRadiator(CherenkovDetector *det, const char *name, const G4LogicalVolume *volume, 
 				 const G4RadiatorMaterial *material) {

include/CherenkovRadiator.h

@@ -8,6 +8,7 @@
 #include <TVector3.h>
 #include <TString.h>
 class TH1D;
+class TCanvas;
 
 #include "ParametricSurface.h"
 class G4LogicalVolume;
@@ -73,8 +74,9 @@ class CherenkovRadiator: public TObject {
    m_ID(0), m_TrajectoryBinCount(1), m_Smearing(0.0), 
    m_GaussianSmearing(false), m_CalibrationPhotonCount(0), m_DetectedPhotonCount(0), m_YieldStat(0), 
    m_YieldCff(0.0), m_DetectedToCalibrationPhotonRatio(0.0), 
-   m_UsedInRingImaging(false), m_Plots(0) {
-    m_LogicalVolumes.push_back(volume);
+   m_UsedInRingImaging(false), m_Plots(0),
+   m_OutputPlotVisualizationEnabled(false), m_wtopx(0), m_wtopy(0), m_wx(0), m_wy(0) {
+   m_LogicalVolumes.push_back(volume);
   };
   ~CherenkovRadiator() {};
 
@@ -182,12 +184,20 @@ class CherenkovRadiator: public TObject {
     return this;
   };
   CherenkovRadiatorPlots *Plots( void ) const { return m_Plots; };
-  void DisplayStandardPlots(const char *cname, int wtopx, unsigned wtopy, unsigned wx, unsigned wy) const;
+  TCanvas *DisplayStandardPlots(const char *cname, const char *wname,
+				int wtopx, unsigned wtopy, unsigned wx, unsigned wy) const;
 
   CherenkovRadiatorPlots  *m_Plots;                         //!
 
   G4DataInterpolation *m_RefractiveIndex;                   //!
 
+  bool m_OutputPlotVisualizationEnabled;                    //!
+  int m_wtopx;                                              //!
+  unsigned m_wtopy;                                         //!
+  unsigned m_wx;                                            //!
+  unsigned m_wy;                                            //!
+  //TCanvas *m_cv;                                            //!
+
 #ifndef DISABLE_ROOT_IO
   ClassDef(CherenkovRadiator, 9);
 #endif

include/ReconstructionFactory.h

@@ -1,6 +1,8 @@
 #pragma once
 
+#include <TH1D.h>
 class TParticlePDG;
+class TCanvas;
 
 #include "Calibration.h"
 #include "Digitization.h"
@@ -9,7 +11,14 @@ namespace IRT2 {
 
 struct ReconstructionFactoryPlots {
   ReconstructionFactoryPlots();
-  ~ReconstructionFactoryPlots() {};
+  ~ReconstructionFactoryPlots() {
+    delete m_hnpe;
+    delete m_hnhits;
+    delete m_hccdftr;
+    delete m_hccdfev;
+    delete m_hdtph;
+    delete m_hmatch;
+  };
 
   // Monte-Carlo plot(s);
   TH1D *hnpe()              const { return m_hnpe;  };
@@ -29,7 +38,9 @@ class ReconstructionFactory : public Digitization, public Calibration {
  public:
   ReconstructionFactory(const char *dfname = 0, const char *cfname = 0, const char *dname = 0);
   ReconstructionFactory(CherenkovDetectorCollection *geometry, CherenkovDetector *cdet, CherenkovEvent *event);
-  virtual ~ReconstructionFactory() {};
+  virtual ~ReconstructionFactory() {
+    if (m_Plots) delete m_Plots;
+  };
 
   void IgnoreTimingInChiSquare( void )               { m_UseTimingInChiSquare = false; };
   void IgnorePoissonTermInChiSquare( void )          { m_UsePoissonTermInChiSquare = false; };
@@ -63,9 +74,11 @@ class ReconstructionFactory : public Digitization, public Calibration {
 
   void InitializePlots( void ) { m_Plots = new ReconstructionFactoryPlots(); };
   ReconstructionFactoryPlots *Plots( void ) const { return m_Plots; };
-  void DisplayStandardPlots(const char *cname, int wtopx, unsigned wtopy, unsigned wx, unsigned wy) const;
+  TCanvas *DisplayStandardPlots(const char *cname, int wtopx, unsigned wtopy, unsigned wx, unsigned wy) const;
 
   void SetHitCountCutoff(unsigned value) { m_HitCountCutoff = value; };
+
+  unsigned GetProcessedEventCount() const { return m_ProcessedEventCount; };
 
  private:
   bool m_VerboseMode;
@@ -91,6 +104,7 @@ class ReconstructionFactory : public Digitization, public Calibration {
   ReconstructionFactoryPlots  *m_Plots;      
 
   unsigned m_HitCountCutoff;
+  unsigned m_ProcessedEventCount;
 
   bool BeVerbose( void )    const { return m_VerboseMode; };
   void LaunchRingFinder(bool calibration);

source/Calibration.cc

@@ -23,7 +23,8 @@ Calibration::Calibration():
   m_UseActsTracking(true),//false),
   m_DefaultSinglePhotonThetaResolution(_SPE_THETA_RESOLUTION_DEFAULT_)
 {
-  m_DatabasePDG = new TDatabasePDG();
+  //m_DatabasePDG = new TDatabasePDG();
+  m_DatabasePDG = TDatabasePDG::Instance();
 
   memset(m_CalibrationBinStat, 0x00, sizeof(m_CalibrationBinStat));
 } // Calibration::Calibration()

source/CherenkovDetector.cc

@@ -0,0 +1,160 @@
+
+#include <iostream>
+#include <fstream>
+
+#ifdef JSON_EXPORT
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+#endif
+
+#include <CherenkovDetector.h>
+
+namespace IRT2 {
+
+// -------------------------------------------------------------------------------------
+//
+// Cut'n'paste from https://stackoverflow.com/questions/40029100/c-how-to-write-back-json-into-file;
+//// Source - https://stackoverflow.com
+// Posted by Danie
+// Retrieved 2025-12-21, License - CC BY-SA 4.0
+
+#ifdef JSON_EXPORT
+template<typename ITERATOR_TYPE>
+void print_element(ITERATOR_TYPE it, unsigned int depth, std::ostream& output_stream);
+void print_array(nlohmann::json::array_t arr, unsigned int depth, std::ostream& output_stream);
+void print_object(nlohmann::json obj, unsigned int depth, std::ostream& output_stream);
+
+void export_json(nlohmann::json json_file, std::ostream& output_stream) {
+    output_stream << "{\n";
+
+    for (nlohmann::json::iterator iter = json_file.begin();;) {
+        output_stream << '\t' << '"' << iter.key() << "\": ";
+
+        print_element(iter, 1, output_stream);
+
+        if (++iter != json_file.end()) {
+            output_stream << ',' << '\n';
+        } else {
+            break;
+        }
+    }
+
+    output_stream << "\n}";
+}
+
+// Adjust it in a way arrays are printed in a single line;
+void print_array(nlohmann::json::array_t arr, unsigned int depth, std::ostream& output_stream) {
+    if (!arr.size())
+        return;
+
+    //std::string indent = std::string(depth, '\t');
+    for (nlohmann::json::array_t::iterator iter = arr.begin();;) {
+      //output_stream << '\n' << indent << '\t';
+      output_stream << ' ';
+        print_element(iter, depth + 1, output_stream);
+
+        if (++iter != arr.end()) {
+            output_stream << ',';
+        } else {
+            break;
+        }
+    }
+    //output_stream << '\n' << indent;
+}
+
+void print_object(nlohmann::json obj, unsigned int depth, std::ostream& output_stream) {
+    if (!obj.size())
+        return;
+
+    std::string indent = std::string(depth, '\t');
+    for (nlohmann::json::iterator iter = obj.begin();;) {
+        output_stream << '\n' << indent << '\t' << '"' << iter.key() << "\": ";
+        print_element(iter, depth + 1, output_stream);
+
+        if (++iter != obj.end()) {
+            output_stream << ',';
+        } else {
+            break;
+        }
+    }
+    output_stream << '\n' << indent;
+}
+
+template<typename ITERATOR_TYPE>
+void print_element(ITERATOR_TYPE iter, unsigned int depth, std::ostream& output_stream) {
+    switch (iter->type())
+    {
+    case nlohmann::json::value_t::object:
+      output_stream << '{';
+      print_object(nlohmann::json(*iter), depth, output_stream);
+      output_stream << '}';
+      break;
+
+    case nlohmann::json::value_t::array:     
+      output_stream << '[';
+      print_array(nlohmann::json::array_t(*iter), depth, output_stream);
+      output_stream << ']';
+      break;
+
+    default:
+      output_stream << *iter;
+      break;
+    }
+}
+#endif
+
+// -------------------------------------------------------------------------------------
+
+void CherenkovDetector::ExportJsonFormatCalibrations(const char *fname)
+{
+#ifdef JSON_EXPORT
+  std::ofstream ofile(fname);
+
+  if (!ofile.is_open())  {
+    std::cout << "\n Failed to open output file";
+  } else {
+    nlohmann::json data;
+    auto &r1data = data["Radiators"] = nlohmann::json::object_t();
+    for(auto [name,radiator]: Radiators()) {
+      //auto radiator = rptr.second;
+      auto &r2data = r1data[/*rptr.first*/name] = nlohmann::json::object_t();
+
+      auto &r3data = r2data["theta-bins"] = nlohmann::json::object_t();
+      for(unsigned iq=0; iq<radiator->m_Calibrations.size(); iq++) {
+	auto &calib = radiator->m_Calibrations[iq];
+
+	if (calib.m_Stat) {
+	  std::vector<std::string> ri_strings;
+	  TString bin, stat, zvtx, offset, sigma, ristr;
+	  bin.Form("%02d", iq);
+	  stat.Form  (  "%6d",      calib.m_Stat);
+	  zvtx.Form  ("%7.1f",      calib.m_AverageZvtx);
+	  // '1000': prefer to have them in a readable [mrad] format;
+	  offset.Form("%7.2f", 1000.*calib.m_Coffset);
+	  sigma.Form ("%7.2f", 1000.*calib.m_Csigma);
+	  //r3data[bin.Data()] = {stat.Data(), zvtx.Data(), offset.Data(), sigma.Data()};
+	  ri_strings.push_back(stat.Data());
+	  ri_strings.push_back(zvtx.Data());
+	  ri_strings.push_back(offset.Data());
+	  ri_strings.push_back(sigma.Data());
+
+	  for(auto ri: calib.m_AverageRefractiveIndices) {
+	    ristr.Form("%8.6f", ri);
+	    ri_strings.push_back(ristr.Data());
+	  } //for ri
+
+	  //r2data["refractive-indices"] = ri_strings;
+	  r3data[bin.Data()] = ri_strings;
+	} //if
+      } //for calib
+    } //for rptr
+
+    export_json(data, ofile);
+    ofile.close();
+  } //if
+#endif
+} // CherenkovDetector::ExportJsonFormatCalibrations()
+
+} // namespace IRT2
+
+// -------------------------------------------------------------------------------------

source/CherenkovDetectorCollection.cc

@@ -0,0 +1,22 @@
+
+#include <CherenkovDetectorCollection.h>
+
+namespace IRT2 {
+
+CherenkovDetectorCollection *CherenkovDetectorCollection::m_Instance = 0;
+
+// -------------------------------------------------------------------------------------
+
+CherenkovDetectorCollection::CherenkovDetectorCollection()
+{
+  if (m_Instance) {
+    printf("Attempt to instantiate CherenkovDetectorCollection singleton twice!\n");
+    exit(0);
+  } //if
+
+  m_Instance = this;
+} // CherenkovDetectorCollection::CherenkovDetectorCollection()
+
+} // namespace IRT2
+
+// -------------------------------------------------------------------------------------

source/CherenkovRadiator.cc

@@ -68,12 +68,12 @@ void CherenkovRadiatorPlots::SetCherenkovAngleRange(double min, double max)
 
 // -------------------------------------------------------------------------------------
 
-void CherenkovRadiator::DisplayStandardPlots(const char *cname, int wtopx,
-					     unsigned wtopy, unsigned wx, unsigned wy) const
+TCanvas *CherenkovRadiator::DisplayStandardPlots(const char *cname, const char *wname, int wtopx,
+						 unsigned wtopy, unsigned wx, unsigned wy) const
 {
-  if (!Plots()) return;
-  
-  auto cv = new TCanvas("", cname, wtopx, wtopy, wx, wy);
+  if (!Plots()) return 0;
+
+  auto cv = new TCanvas(cname, wname, wtopx, wtopy, wx, wy);
   cv->Divide(4, 2);
 
   cv->cd(1); if (Plots()->hvtx()) Plots()->hvtx()   ->Draw();
@@ -84,6 +84,8 @@ void CherenkovRadiator::DisplayStandardPlots(const char *cname, int wtopx,
   cv->cd(6);                      Plots()->hnhits() ->Draw();
   cv->cd(7);                      Plots()->hccdfph()->Draw();
   cv->cd(8);                      Plots()->hthtr()  ->Fit("gaus");
+
+  return cv;
 } // CherenkovRadiator::DisplayStandardPlots()
 
 // -------------------------------------------------------------------------------------