diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d62376..7e95aa7 100644 --- a/CMakeLists.txt +++ b/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 $ $ diff --git a/include/CherenkovDetector.h b/include/CherenkovDetector.h index 382ff41..0d9e583 100644 --- a/include/CherenkovDetector.h +++ b/include/CherenkovDetector.h @@ -134,6 +134,8 @@ class CherenkovDetector: public TObject { // readout ID -> pixel position converter (for external usage) std::function m_ReadoutIDToPosition; //! + void ExportJsonFormatCalibrations(const char *fname); + private: TString m_Name; // This is needed for dd4hep cell index decoding; diff --git a/include/CherenkovDetectorCollection.h b/include/CherenkovDetectorCollection.h index 2692ff0..6dd458c 100644 --- a/include/CherenkovDetectorCollection.h +++ b/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) { diff --git a/include/CherenkovRadiator.h b/include/CherenkovRadiator.h index 816813a..af0609c 100644 --- a/include/CherenkovRadiator.h +++ b/include/CherenkovRadiator.h @@ -8,6 +8,7 @@ #include #include 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 diff --git a/include/ReconstructionFactory.h b/include/ReconstructionFactory.h index 0c3f733..d590c81 100644 --- a/include/ReconstructionFactory.h +++ b/include/ReconstructionFactory.h @@ -1,6 +1,8 @@ #pragma once +#include 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); diff --git a/source/Calibration.cc b/source/Calibration.cc index 583aca8..5085eaa 100644 --- a/source/Calibration.cc +++ b/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() diff --git a/source/CherenkovDetector.cc b/source/CherenkovDetector.cc new file mode 100644 index 0000000..ad63f82 --- /dev/null +++ b/source/CherenkovDetector.cc @@ -0,0 +1,160 @@ + +#include +#include + +#ifdef JSON_EXPORT +#include +using json = nlohmann::json; +#endif + +#include + +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 +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 +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; iqm_Calibrations.size(); iq++) { + auto &calib = radiator->m_Calibrations[iq]; + + if (calib.m_Stat) { + std::vector 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 + +// ------------------------------------------------------------------------------------- diff --git a/source/CherenkovDetectorCollection.cc b/source/CherenkovDetectorCollection.cc new file mode 100644 index 0000000..d719638 --- /dev/null +++ b/source/CherenkovDetectorCollection.cc @@ -0,0 +1,22 @@ + +#include + +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 + +// ------------------------------------------------------------------------------------- diff --git a/source/CherenkovRadiator.cc b/source/CherenkovRadiator.cc index 0b3d423..9f76644 100644 --- a/source/CherenkovRadiator.cc +++ b/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() // ------------------------------------------------------------------------------------- diff --git a/source/ReconstructionFactory.cc b/source/ReconstructionFactory.cc index 44665c2..2f2474c 100644 --- a/source/ReconstructionFactory.cc +++ b/source/ReconstructionFactory.cc @@ -35,7 +35,8 @@ ReconstructionFactory::ReconstructionFactory(const char *dfname, const char *cfn m_UseMcTruthPhotonDirectionSeed(true), m_Plots(0), // Require at least one associated hit per default; - m_HitCountCutoff(1) + m_HitCountCutoff(1), + m_ProcessedEventCount(0) { } // ReconstructionFactory::ReconstructionFactory() // ------------------------------------------------------------------------------------- @@ -54,7 +55,8 @@ ReconstructionFactory::ReconstructionFactory(CherenkovDetectorCollection *geomet m_UseMcTruthPhotonDirectionSeed(true), m_Plots(0), // Require at least one associated hit per default; - m_HitCountCutoff(1) + m_HitCountCutoff(1), + m_ProcessedEventCount(0) { } // ReconstructionFactory::ReconstructionFactory() @@ -347,9 +349,11 @@ void ReconstructionFactory::LaunchRingFinder(bool calibration) // And eventually calculate the overall event chi^2; { +#ifdef _EXPECTED_BG_PHOTON_COUNT_ unsigned nbg = 0; for(auto &hit: Hits()) if (hit.m_BackgroundCandidate) nbg++; +#endif unsigned ndfev = 0, idx = 0; double chi2ev = 0; @@ -580,6 +584,8 @@ CherenkovEvent *ReconstructionFactory::GetEvent(unsigned ev, bool calibration) } //for mcparticle if (BeVerbose() && !(ev%100)) printf("Event %5d ...\n", ev); + + m_ProcessedEventCount++; return Event(); } // ReconstructionFactory::GetEvent() @@ -602,13 +608,13 @@ ReconstructionFactoryPlots::ReconstructionFactoryPlots( void ) // ------------------------------------------------------------------------------------- -void ReconstructionFactory::DisplayStandardPlots(const char *cname, int wtopx, - unsigned wtopy, unsigned wx, unsigned wy) const +TCanvas *ReconstructionFactory::DisplayStandardPlots(const char *cname, int wtopx, + unsigned wtopy, unsigned wx, unsigned wy) const { - if (!Plots()) return; - - auto cv = new TCanvas("", cname, wtopx, wtopy, wx, wy); - cv->Divide(2, 4); + if (!Plots()) return 0; + + auto cv = new TCanvas("cx", cname, wtopx, wtopy, wx, wy); + cv->Divide(2, 3); cv->cd(1); Plots()->hnpe() ->Draw(); cv->cd(2); Plots()->hnhits() ->Draw(); @@ -616,6 +622,8 @@ void ReconstructionFactory::DisplayStandardPlots(const char *cname, int wtopx, cv->cd(4); Plots()->hccdfev()->Draw(); cv->cd(5); Plots()->hdtph() ->Fit("gaus"); cv->cd(6); Plots()->hmatch() ->Draw(); + + return cv; } // ReconstructionFactory::DisplayStandardPlots() // -------------------------------------------------------------------------------------