diff --git a/Makefile b/Makefile index 39638f1..c394a6c 100755 --- a/Makefile +++ b/Makefile @@ -54,17 +54,26 @@ else $(error "Unable to find HDF5 headers and libraries") endif +CCFITSCXXFLAGS = +CCFITSLIBS = +ifeq ("$(shell pkg-config --exists CCfits 1>&2 2> /dev/null; echo $$?)", "0") +CCFITSCXXFLAGS += $(shell pkg-config --cflags CCfits) +CCFITSLIBS += $(shell pkg-config --libs CCfits) +else +$(error "Unable to find CCfits headers and libraries") +endif + #---------------------------------------------------------------- # Definitions based on architecture and user options # CMD="" -CXXFLAGS += -I$(IN) -I$(MEGALIB)/include -I/opt/local/include $(H5CXXFLAGS) +CXXFLAGS += -I$(IN) -I$(MEGALIB)/include -I/opt/local/include $(H5CXXFLAGS) $(CCFITSCXXFLAGS) # Comment this line out if you want to accept warnings #CXXFLAGS += -Werror -Wno-unused-variable -LIBS += $(H5LIBS) +LIBS += $(H5LIBS) $(CCFITSLIBS) # Definitions NUCLEARIZER_DIR := $(NUCLEARIZER) @@ -162,7 +171,7 @@ $(NUCLEARIZER_DICT): $(FRETALON_H_FILES) $(NUCLEARIZER_H_FILES) @echo "Generating LinkDef ..." @$(BN)/generatelinkdef -o $(NUCLEARIZER_LINKDEF) -i $(NUCLEARIZER_H_FILES) $(FRETALON_H_FILES) @echo "Generating dictionary ..." - @rootcling -f $@ -I$(IN) -I$(MEGALIB)/include $(H5CXXFLAGS) -D___CLING___ -rmf $(NUCLEARIZER_ROOTMAP) -s libNuclearizer -c $(NUCLEARIZER_H_FILES) $(FRETALON_H_FILES) $(NUCLEARIZER_LINKDEF) + @rootcling -f $@ -I$(IN) -I$(MEGALIB)/include $(H5CXXFLAGS) $(CCFITSCXXFLAGS) -D___CLING___ -rmf $(NUCLEARIZER_ROOTMAP) -s libNuclearizer -c $(NUCLEARIZER_H_FILES) $(FRETALON_H_FILES) $(NUCLEARIZER_LINKDEF) @mv $(NUCLEARIZER_ROOTPCM) $(LB) $(NUCLEARIZER_DICT_LIB): $(NUCLEARIZER_DICT) diff --git a/include/MGUIOptionsLoaderMeasurementsFITS.h b/include/MGUIOptionsLoaderMeasurementsFITS.h new file mode 100644 index 0000000..1d2d9aa --- /dev/null +++ b/include/MGUIOptionsLoaderMeasurementsFITS.h @@ -0,0 +1,87 @@ +/* + * MGUIOptionsLoaderMeasurementsFITS.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MGUIOptionsLoaderMeasurementsFITS__ +#define __MGUIOptionsLoaderMeasurementsFITS__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIEFileSelector.h" +#include "MGUIOptions.h" + +// Nuclearizer libs: +#include "MModule.h" + + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +//! UI settings for the FITS measurements loader +class MGUIOptionsLoaderMeasurementsFITS : public MGUIOptions +{ + // public Session: + public: + //! Default constructor + MGUIOptionsLoaderMeasurementsFITS(MModule* Module); + //! Default destructor + virtual ~MGUIOptionsLoaderMeasurementsFITS(); + + //! Process all button, etc. messages + virtual bool ProcessMessage(long Message, long Parameter1, long Parameter2); + + //! The creation part which gets overwritten + virtual void Create(); + + // protected methods: + protected: + + //! Actions after the Apply or OK button has been pressed + virtual bool OnApply(); + + + // protected members: + protected: + + // private members: + private: + //! Select which FITS file to load + MGUIEFileSelector* m_FileSelectorFITS; + + + +#ifdef ___CLING___ + public: + ClassDef(MGUIOptionsLoaderMeasurementsFITS, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MGUIOptionsSaverMeasurementsFITS.h b/include/MGUIOptionsSaverMeasurementsFITS.h new file mode 100644 index 0000000..31da9b2 --- /dev/null +++ b/include/MGUIOptionsSaverMeasurementsFITS.h @@ -0,0 +1,87 @@ +/* + * MGUIOptionsSaverMeasurementsFITS.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MGUIOptionsSaverMeasurementsFITS__ +#define __MGUIOptionsSaverMeasurementsFITS__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIEFileSelector.h" +#include "MGUIOptions.h" + +// Nuclearizer libs: +#include "MModule.h" + + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +//! UI settings for the FITS measurements saver +class MGUIOptionsSaverMeasurementsFITS : public MGUIOptions +{ + // public Session: + public: + //! Default constructor + MGUIOptionsSaverMeasurementsFITS(MModule* Module); + //! Default destructor + virtual ~MGUIOptionsSaverMeasurementsFITS(); + + //! Process all button, etc. messages + virtual bool ProcessMessage(long Message, long Parameter1, long Parameter2); + + //! The creation part which gets overwritten + virtual void Create(); + + // protected methods: + protected: + + //! Actions after the Apply or OK button has been pressed + virtual bool OnApply(); + + + // protected members: + protected: + + // private members: + private: + //! Select which FITS file to save to + MGUIEFileSelector* m_FileSelectorFITS; + + + +#ifdef ___CLING___ + public: + ClassDef(MGUIOptionsSaverMeasurementsFITS, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MModuleLoaderMeasurementsFITS.h b/include/MModuleLoaderMeasurementsFITS.h new file mode 100644 index 0000000..a701c84 --- /dev/null +++ b/include/MModuleLoaderMeasurementsFITS.h @@ -0,0 +1,139 @@ +/* + * MModuleLoaderMeasurementsFITS.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleLoaderMeasurementsFITS__ +#define __MModuleLoaderMeasurementsFITS__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include +#include + +// ROOT libs: + +// MEGAlib libs: +#include "MGlobal.h" +#include "MFileReadOuts.h" + +// Nuclearizer libs: +#include "MModuleLoaderMeasurements.h" + +// CCfits libs +#include +using namespace CCfits; + + +//////////////////////////////////////////////////////////////////////////////// + +//! A module to load FITS data files +class MModuleLoaderMeasurementsFITS : public MModuleLoaderMeasurements +{ + // public interface: + public: + //! Default constructor + MModuleLoaderMeasurementsFITS(); + //! Default destructor + virtual ~MModuleLoaderMeasurementsFITS(); + + //! Create a new object of this class + virtual MModuleLoaderMeasurementsFITS* Clone() { return new MModuleLoaderMeasurementsFITS(); } + + //! Initialize the module + virtual bool Initialize(); + + //! Finalize the module + virtual void Finalize(); + + //! Main data analysis routine, which updates the event to a new level + virtual bool AnalyzeEvent(MReadOutAssembly* Event); + + //! Show the options GUI + virtual void ShowOptionsGUI(); + + //! Read the configuration data from an XML node + virtual bool ReadXmlConfiguration(MXmlNode* Node); + //! Create an XML node tree from the configuration + virtual MXmlNode* CreateXmlConfiguration(); + + + // protected methods: + protected: + //! Open the FITS file + bool OpenFITSFile(MString FileName); + //! Read next event (loads new batch if needed) + bool ReadBatch(); + + // private methods: + private: + + + // protected members: + protected: + + + // private members: + private: + //! The FITS file object pointer + FITS* m_FITSFile; + + //! Primary HDU (HDU 0) - header/metadata + const PHDU* m_PrimaryHDU; + + //! Compton L1a table extension (extension 1) + const ExtHDU* m_ComptonTable; + + //! Current row number in the FITS table + long m_CurrentRow; + + //! Total number of rows in the FITS table + long m_TotalRows; + + //! Batch size for reading FITS data + static const long m_BatchSize = 100; + + //! Current batch start row + long m_BatchStartRow = 1; + + //! Current event index within the batch + long m_CurrentEventInBatch = 0; + + //! Number of events in current batch + long m_BatchEventCount = 0; + + //! Batch data storage for scalar columns + std::vector m_BatchTIME; + std::vector m_BatchEVENTTYPE; + std::vector m_BatchNUMSTRIPHIT; + + //! Batch data storage for variable-length array columns + std::vector> m_BatchTYPEHIT; + std::vector> m_BatchDETID; + std::vector> m_BatchSTRIPID; + std::vector> m_BatchSIDEID; + std::vector> m_BatchFASTTIME; + std::vector> m_BatchPHA; + std::vector> m_BatchTAC; + + +#ifdef ___CLING___ + public: + ClassDef(MModuleLoaderMeasurementsFITS, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MModuleSaverMeasurementsFITS.h b/include/MModuleSaverMeasurementsFITS.h new file mode 100644 index 0000000..0e2d906 --- /dev/null +++ b/include/MModuleSaverMeasurementsFITS.h @@ -0,0 +1,152 @@ +/* + * MModuleSaverMeasurementsFITS.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleSaverMeasurementsFITS__ +#define __MModuleSaverMeasurementsFITS__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include +#include + +// ROOT libs: + +// MEGAlib libs: +#include "MGlobal.h" +#include "MString.h" + +// Nuclearizer libs: +#include "MModule.h" + +// CCfits libs +#include +using namespace CCfits; + + +//////////////////////////////////////////////////////////////////////////////// + + +//! A module to save hit-level events to FITS data files +class MModuleSaverMeasurementsFITS : public MModule +{ + // public interface: + public: + //! Default constructor + MModuleSaverMeasurementsFITS(); + //! Default destructor + virtual ~MModuleSaverMeasurementsFITS(); + + //! Create a new object of this class + virtual MModuleSaverMeasurementsFITS* Clone() { return new MModuleSaverMeasurementsFITS(); } + + //! Initialize the module + virtual bool Initialize(); + + //! Finalize the module + virtual void Finalize(); + + //! Main data analysis routine, which updates the event to a new level + virtual bool AnalyzeEvent(MReadOutAssembly* Event); + + //! Show the options GUI + virtual void ShowOptionsGUI(); + + //! Read the configuration data from an XML node + virtual bool ReadXmlConfiguration(MXmlNode* Node); + //! Create an XML node tree from the configuration + virtual MXmlNode* CreateXmlConfiguration(); + + //! Set the output file name + void SetFileName(const MString& FileName) { m_FileName = FileName; } + //! Get the output file name + MString GetFileName() const { return m_FileName; } + + + // protected methods: + protected: + //! Create the FITS file and extensions + bool CreateFITSFile(MString FileName); + //! Flush current batch to FITS file + bool FlushBatch(); + + + // private methods: + private: + + + // protected members: + protected: + + + // private members: + private: + //! Output file name + MString m_FileName; + + //! The FITS file object pointer + FITS* m_FITSFile; + + //! Primary HDU (HDU 0) - header/metadata + PHDU* m_PrimaryHDU; + + //! Science data table extension + ExtHDU* m_ScienceTable; + + //! Total number of events written + long m_TotalEventsWritten; + + //! Batch size for writing FITS data + static const long m_BatchSize = 100; + + //! Current row where next batch will be written + long m_BatchStartRow; + + //! Number of events in current batch + long m_BatchEventCount; + + //! Batch data storage for scalar columns + std::vector m_BatchTIME; + std::vector m_BatchEVENTTYPE; + std::vector m_BatchEVENTCLASS; + std::vector m_BatchNUMHIT; + std::vector m_BatchSEQHIT; + + //! Batch data storage for fixed-length array columns (event-level) + std::vector> m_BatchSTATTEST; + std::vector> m_BatchRECOILDIR; + std::vector> m_BatchRECOILDIR_ERR; + + //! Batch data storage for variable-length array columns (hit-level data) + std::vector> m_BatchX; + std::vector> m_BatchY; + std::vector> m_BatchZ; + std::vector> m_BatchX_ERR; + std::vector> m_BatchY_ERR; + std::vector> m_BatchZ_ERR; + std::vector> m_BatchENERGY; + std::vector> m_BatchENERGY_ERR; + std::vector> m_BatchBAD_FLAG; + + +#ifdef ___CLING___ + public: + ClassDef(MModuleSaverMeasurementsFITS, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MAssembly.cxx b/src/MAssembly.cxx index f1e3721..d96e52c 100644 --- a/src/MAssembly.cxx +++ b/src/MAssembly.cxx @@ -65,6 +65,7 @@ using namespace std; #include "MModuleLoaderSimulationsCosima.h" #include "MModuleLoaderMeasurementsROA.h" #include "MModuleLoaderMeasurementsHDF.h" +#include "MModuleLoaderMeasurementsFITS.h" #include "MModuleEnergyCalibration.h" #include "MModuleEnergyCalibrationUniversal.h" #include "MModuleDepthCalibration.h" @@ -75,6 +76,7 @@ using namespace std; #include "MModuleStripPairingChiSquare.h" #include "MModuleEventFilter.h" #include "MModuleEventSaver.h" +#include "MModuleSaverMeasurementsFITS.h" #include "MModuleResponseGenerator.h" #include "MModuleTACcut.h" #include "MModuleNearestNeighbor.h" @@ -129,7 +131,8 @@ MAssembly::MAssembly() m_Supervisor->AddAvailableModule(new MModuleLoaderSimulationsCosima()); m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsROA()); m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsHDF()); - + m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsFITS()); + m_Supervisor->AddAvailableModule(new MModuleDEESMEX()); m_Supervisor->AddAvailableModule(new MModuleEventFilter()); @@ -143,6 +146,7 @@ MAssembly::MAssembly() m_Supervisor->AddAvailableModule(new MModuleDepthCalibration2024()); m_Supervisor->AddAvailableModule(new MModuleEventSaver()); + m_Supervisor->AddAvailableModule(new MModuleSaverMeasurementsFITS()); m_Supervisor->AddAvailableModule(new MModuleTransmitterRealta()); m_Supervisor->AddAvailableModule(new MModuleResponseGenerator()); m_Supervisor->AddAvailableModule(new MModuleTACcut()); diff --git a/src/MGUIOptionsLoaderMeasurementsFITS.cxx b/src/MGUIOptionsLoaderMeasurementsFITS.cxx new file mode 100644 index 0000000..be9497e --- /dev/null +++ b/src/MGUIOptionsLoaderMeasurementsFITS.cxx @@ -0,0 +1,126 @@ +/* + * MGUIOptionsLoaderMeasurementsFITS.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIOptionsLoaderMeasurementsFITS.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" +#include "MModuleLoaderMeasurementsFITS.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIOptionsLoaderMeasurementsFITS) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsLoaderMeasurementsFITS::MGUIOptionsLoaderMeasurementsFITS(MModule* Module) + : MGUIOptions(Module) +{ + // standard constructor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsLoaderMeasurementsFITS::~MGUIOptionsLoaderMeasurementsFITS() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIOptionsLoaderMeasurementsFITS::Create() +{ + PreCreate(); + + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + + m_FileSelectorFITS = new MGUIEFileSelector(m_OptionsFrame, "Please select a FITS file:", + dynamic_cast(m_Module)->GetFileName()); + m_FileSelectorFITS->SetFileType("FITS file", "*.fits"); + m_FileSelectorFITS->SetFileType("FITS file", "*.fit"); + m_OptionsFrame->AddFrame(m_FileSelectorFITS, LabelLayout); + + PostCreate(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsLoaderMeasurementsFITS::ProcessMessage(long Message, long Parameter1, long Parameter2) +{ + // Modify here if you have more buttons + + bool Status = true; + + switch (GET_MSG(Message)) { + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + break; + default: + break; + } + break; + default: + break; + } + + if (Status == false) { + return false; + } + + // Call also base class + return MGUIOptions::ProcessMessage(Message, Parameter1, Parameter2); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsLoaderMeasurementsFITS::OnApply() +{ + // Modify this to store the data in the module! + + dynamic_cast(m_Module)->SetFileName(m_FileSelectorFITS->GetFileName()); + + return true; +} + + +// MGUIOptionsLoaderMeasurementsFITS: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MGUIOptionsSaverMeasurementsFITS.cxx b/src/MGUIOptionsSaverMeasurementsFITS.cxx new file mode 100644 index 0000000..c070a6b --- /dev/null +++ b/src/MGUIOptionsSaverMeasurementsFITS.cxx @@ -0,0 +1,126 @@ +/* + * MGUIOptionsSaverMeasurementsFITS.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIOptionsSaverMeasurementsFITS.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" +#include "MModuleSaverMeasurementsFITS.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIOptionsSaverMeasurementsFITS) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsSaverMeasurementsFITS::MGUIOptionsSaverMeasurementsFITS(MModule* Module) + : MGUIOptions(Module) +{ + // standard constructor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsSaverMeasurementsFITS::~MGUIOptionsSaverMeasurementsFITS() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIOptionsSaverMeasurementsFITS::Create() +{ + PreCreate(); + + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + + m_FileSelectorFITS = new MGUIEFileSelector(m_OptionsFrame, "Please select output FITS file:", + dynamic_cast(m_Module)->GetFileName()); + m_FileSelectorFITS->SetFileType("FITS file", "*.fits"); + m_FileSelectorFITS->SetFileType("FITS file", "*.fit"); + m_OptionsFrame->AddFrame(m_FileSelectorFITS, LabelLayout); + + PostCreate(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsSaverMeasurementsFITS::ProcessMessage(long Message, long Parameter1, long Parameter2) +{ + // Modify here if you have more buttons + + bool Status = true; + + switch (GET_MSG(Message)) { + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + break; + default: + break; + } + break; + default: + break; + } + + if (Status == false) { + return false; + } + + // Call also base class + return MGUIOptions::ProcessMessage(Message, Parameter1, Parameter2); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsSaverMeasurementsFITS::OnApply() +{ + // Modify this to store the data in the module! + + dynamic_cast(m_Module)->SetFileName(m_FileSelectorFITS->GetFileName()); + + return true; +} + + +// MGUIOptionsSaverMeasurementsFITS: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleLoaderMeasurementsFITS.cxx b/src/MModuleLoaderMeasurementsFITS.cxx new file mode 100644 index 0000000..85a072c --- /dev/null +++ b/src/MModuleLoaderMeasurementsFITS.cxx @@ -0,0 +1,409 @@ +/* + * MModuleLoaderMeasurementsFITS.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MModuleLoaderMeasurementsFITS +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MModuleLoaderMeasurementsFITS.h" + +// Standard libs: +#include +using namespace std; + +// ROOT libs: +#include "TGClient.h" + +// MEGAlib libs: +#include "MGUIOptionsLoaderMeasurementsFITS.h" +#include "MReadOut.h" +#include "MReadOutSequence.h" +#include "MReadOutElementDoubleStrip.h" +#include "MReadOutDataADCValue.h" +#include "MReadOutDataTiming.h" +#include "MReadOutDataOrigins.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleLoaderMeasurementsFITS) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleLoaderMeasurementsFITS::MModuleLoaderMeasurementsFITS() : MModuleLoaderMeasurements() +{ + // Construct an instance of MModuleLoaderMeasurementsFITS + + // Set all module relevant information + + // Set the module name --- has to be unique + m_Name = "Measurement loader for FITS files"; + + // Set the XML tag --- has to be unique --- no spaces allowed + m_XmlTag = "XmlTagMeasurementLoaderFITS"; + + // This is a special start module which can generate its own events + m_IsStartModule = true; + + // Allow the use of multiple threads and instances + m_AllowMultiThreading = true; + m_AllowMultipleInstances = false; + + m_FITSFile = nullptr; + m_PrimaryHDU = nullptr; + m_ComptonTable = nullptr; + m_CurrentRow = 1; // FITS uses 1-based row indexing + m_TotalRows = 0; +} + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleLoaderMeasurementsFITS::~MModuleLoaderMeasurementsFITS() +{ + // Delete this instance of MModuleLoaderMeasurementsFITS +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleLoaderMeasurementsFITS::Initialize() +{ + // Initialize the module + + // Clean: + m_FileType = "Unknown"; + + if (MFile::Exists(m_FileName) == false) { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Error) cout<pHDU(); + + // Try to get Compton extension by index 1 (first extension after primary) + // Note: FITS files have HDU 0 = Primary, HDU 1 = first extension + try { + m_ComptonTable = &m_FITSFile->extension(1); + if (g_Verbosity >= c_Info) { + cout<<"Opened extension at index 1: "<name()<extension("Science Data Table 1st Extension"); + } + + //Nothing in the GTI Table now. + //m_GtiTable = &m_FITSFile->extension(2); + + // Get table metadata + m_TotalRows = m_ComptonTable->rows(); + int nCols = m_ComptonTable->numCols(); + + cout<<"FITS table metadata:"<column(); + + // Print column information to verify structure, can comment out + cout<<" Column details:"<type() + <<", width: "<width() + <<", repeat: "<repeat() + <<")"<= m_BatchEventCount) { + + // Check if we reach the end of file + if (m_BatchStartRow > m_TotalRows) { + return false; // No more data + } + + try { + // Calculate how many rows to read in this batch + long rowsToRead = std::min(m_BatchSize, m_TotalRows - m_BatchStartRow + 1); + long lastRow = m_BatchStartRow + rowsToRead - 1; + + // Resize vectors to read this batch + m_BatchTIME.resize(rowsToRead); + m_BatchEVENTTYPE.resize(rowsToRead); + m_BatchNUMSTRIPHIT.resize(rowsToRead); + m_BatchTYPEHIT.resize(rowsToRead); + m_BatchDETID.resize(rowsToRead); + m_BatchSTRIPID.resize(rowsToRead); + m_BatchSIDEID.resize(rowsToRead); + m_BatchFASTTIME.resize(rowsToRead); + m_BatchPHA.resize(rowsToRead); + m_BatchTAC.resize(rowsToRead); + + // Read scalar columns - third parameter is LAST row index, NOT count + m_ComptonTable->column("TIME").read(m_BatchTIME, m_BatchStartRow, lastRow); + m_ComptonTable->column("EVENTTYPE").read(m_BatchEVENTTYPE, m_BatchStartRow, lastRow); + m_ComptonTable->column("NUMSTRIPHIT").read(m_BatchNUMSTRIPHIT, m_BatchStartRow, lastRow); + + // Read variable-length array columns + m_ComptonTable->column("TYPEHIT").readArrays(m_BatchTYPEHIT, m_BatchStartRow, lastRow); + m_ComptonTable->column("DETID").readArrays(m_BatchDETID, m_BatchStartRow, lastRow); + m_ComptonTable->column("STRIPID").readArrays(m_BatchSTRIPID, m_BatchStartRow, lastRow); + m_ComptonTable->column("SIDEID").readArrays(m_BatchSIDEID, m_BatchStartRow, lastRow); + m_ComptonTable->column("FASTTIME").readArrays(m_BatchFASTTIME, m_BatchStartRow, lastRow); + m_ComptonTable->column("PHA").readArrays(m_BatchPHA, m_BatchStartRow, lastRow); + m_ComptonTable->column("TAC").readArrays(m_BatchTAC, m_BatchStartRow, lastRow); + + // Update batch tracking + m_BatchEventCount = rowsToRead; + m_CurrentEventInBatch = 0; + m_BatchStartRow += rowsToRead; + + // Display first 3 events in this batch to verify event, can comment out + // cout<<" First "<SetID(); // TODO: No EventID + Event->SetCL(eventTime); // Mission time in seconds + + // Loop through strip hits and create MStripHit objects + for (uint8_t hitIdx = 0; hitIdx < numStripHit; ++hitIdx) { + // Extract hit data from arrays + uint8_t typeHit = m_BatchTYPEHIT[idx][hitIdx]; + int detID = m_BatchDETID[idx][hitIdx]; + int stripID = m_BatchSTRIPID[idx][hitIdx]; + int sideID = m_BatchSIDEID[idx][hitIdx]; + uint8_t fastTime = m_BatchFASTTIME[idx][hitIdx]; + int pha = m_BatchPHA[idx][hitIdx]; + int tac = m_BatchTAC[idx][hitIdx]; + + // Create new strip hit + MStripHit* H = new MStripHit(); + + // Set detector and strip information + H->SetDetectorID(detID); + H->SetStripID(stripID); + H->IsLowVoltageStrip(sideID == 0); + + // Set measured data + H->SetADCUnits(pha); // I think pha is the ADCunits? + H->SetTAC(tac); // Timing + + // Set boolean flags based on hit type + H->IsGuardRing(typeHit == 2); + if (H->IsGuardRing() == true) { + Event->SetGuardRingVeto(true); + } + H->IsNearestNeighbor(typeHit == 1); + H->HasFastTiming(fastTime == 1); + + // Add strip hit to event + Event->AddStripHit(H); + + // DEBUG + // if (m_CurrentRow > m_TotalRows - 10) { + // cout<<"Event "<GetNode("FileNameFITS"); + if (FileNameFITSNode != nullptr) { + m_FileName = FileNameFITSNode->GetValue(); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleLoaderMeasurementsFITS::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + new MXmlNode(Node, "FileNameFITS", m_FileName); + + return Node; +} + + +/////////////////////////////////////////////////////////////////////////////// + + +void MModuleLoaderMeasurementsFITS::ShowOptionsGUI() +{ + //! Show the options GUI + + MGUIOptionsLoaderMeasurementsFITS* Options = new MGUIOptionsLoaderMeasurementsFITS(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +// MModuleLoaderMeasurementsFITS.cxx: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleSaverMeasurementsFITS.cxx b/src/MModuleSaverMeasurementsFITS.cxx new file mode 100644 index 0000000..3848169 --- /dev/null +++ b/src/MModuleSaverMeasurementsFITS.cxx @@ -0,0 +1,455 @@ +/* + * MModuleSaverMeasurementsFITS.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MModuleSaverMeasurementsFITS +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MModuleSaverMeasurementsFITS.h" + +// Standard libs: +#include +using namespace std; + +// ROOT libs: +#include "TGClient.h" + +// MEGAlib libs: +#include "MGUIOptionsSaverMeasurementsFITS.h" +#include "MHit.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleSaverMeasurementsFITS) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleSaverMeasurementsFITS::MModuleSaverMeasurementsFITS() : MModule() +{ + // Construct an instance of MModuleSaverMeasurementsFITS + + // Set all module relevant information + + // Set the module name --- has to be unique + m_Name = "Save events to FITS files (L1b)"; + + // Set the XML tag --- has to be unique --- no spaces allowed + m_XmlTag = "XmlTagSaverMeasurementsFITS"; + + // Set all modules, which have to be done before this module + AddPreceedingModuleType(MAssembly::c_EventLoader); + AddPreceedingModuleType(MAssembly::c_StripPairing); + + // Set all types this modules handles + AddModuleType(MAssembly::c_EventSaver); + + // Set if this module has an options GUI + m_HasOptionsGUI = true; + + // Allow the use of multiple threads and instances + m_AllowMultiThreading = true; + m_AllowMultipleInstances = false; + + m_FITSFile = nullptr; + m_PrimaryHDU = nullptr; + m_ScienceTable = nullptr; + m_TotalEventsWritten = 0; + m_BatchStartRow = 1; + m_BatchEventCount = 0; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleSaverMeasurementsFITS::~MModuleSaverMeasurementsFITS() +{ + // Delete this instance of MModuleSaverMeasurementsFITS +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleSaverMeasurementsFITS::Initialize() +{ + // Initialize the module + + if (m_FileName == "") { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Info) cout<pHDU(); + + // Add some keywords to primary HDU + m_PrimaryHDU->addKey("CREATOR", "Nuclearizer", "Software that created this file"); + m_PrimaryHDU->addKey("ORIGIN", "UC Berkeley SSL", "Organization"); + m_PrimaryHDU->addKey("TELESCOP", "COSI", "Mission name"); + m_PrimaryHDU->addKey("INSTRUME", "GeD", "Instrument name"); + + // Define columns for science data table per specification + // PE(100) = variable-length single-precision float array (max 100) + // 4E = fixed-length array of 4 single-precision floats + // 3E = fixed-length array of 3 single-precision floats + std::vector colNames = { + "TIME", "EVENTTYPE", "EVENTCLASS", "NUMHIT", "SEQHIT", + "STATTEST", "RECOILDIR", "RECOILDIR_ERR", + "X", "Y", "Z", + "X_ERR", "Y_ERR", "Z_ERR", + "ENERGY", "ENERGY_ERR", "BAD_FLAG" + }; + + std::vector colFormats = { + "1D", // TIME - scalar double + "1B", // EVENTTYPE - scalar byte + "1B", // EVENTCLASS - scalar byte + "1B", // NUMHIT - scalar byte + "1B", // SEQHIT - scalar byte + "4E", // STATTEST + "3E", // RECOILDIR + "3E", // RECOILDIR_ERR + "PE(100)", // X - variable-length float array + "PE(100)", // Y + "PE(100)", // Z + "PE(100)", // X_ERR + "PE(100)", // Y_ERR + "PE(100)", // Z_ERR + "PE(100)", // ENERGY + "PE(100)", // ENERGY_ERR + "PE(100)" // BAD_FLAG + }; + + std::vector colUnits = { + "s", "", "", "", "", + "", "unit", "unit", + "cm", "unit", "unit", + "unit", "unit", "unit", + "keV", "unit", "" + }; + + // Create binary table extension + m_ScienceTable = m_FITSFile->addTable("Compton_L1b_1st_Ext", 0, colNames, colFormats, colUnits); + + // Add keywords to science table + m_ScienceTable->addKey("EXTNAME", "COMPTON_L1B", "name of this HDU"); + m_ScienceTable->addKey("TELESCOP", "COSI", "Telescope mission name"); + m_ScienceTable->addKey("INSTRUME", "GED", "Instrument name"); + m_ScienceTable->addKey("DATAMODE", "string", "Instrument datamode"); + m_ScienceTable->addKey("OBSERVER", "string", "Principal Investigator"); + m_ScienceTable->addKey("OBS_ID", "YYMMDD", "Observation ID"); + m_ScienceTable->addKey("OBJECT", "string", "Object/Target name"); + m_ScienceTable->addKey("MJDREFI", 60676, "MJD reference day 01 Jan 2025 00:00:00"); + m_ScienceTable->addKey("MJDREFF", 8.007407407407E-04, "MJD reference (fraction of day)"); + m_ScienceTable->addKey("TIMEREF", "LOCAL", "Reference Frame"); + m_ScienceTable->addKey("TASSIGN", "SATELLITE", "Time assigned"); + m_ScienceTable->addKey("TIMESYS", "TT", "Time System"); + m_ScienceTable->addKey("TIMEUNIT", "s", "Time unit for timing header keywords"); + m_ScienceTable->addKey("TIMEDEL", 0.0, "Integration time"); + m_ScienceTable->addKey("CLOCKAPP", false, "If clock corrections are applied (T/F)"); + m_ScienceTable->addKey("DATE-OBS", "yyyy-mm-ddThh:mm:ss", "Start Date"); + m_ScienceTable->addKey("DATE-END", "yyyy-mm-ddThh:mm:ss", "Stop Date"); + m_ScienceTable->addKey("TSTART", 0.0, "Start time"); + m_ScienceTable->addKey("TSTOP", 0.0, "Stop time"); + m_ScienceTable->addKey("HDUCLASS", "OGIP", "format conforms to OGIP standard"); + m_ScienceTable->addKey("HDUCLAS1", "ARRAY", "hduclass1"); + m_ScienceTable->addKey("HDUCLAS2", "TOTAL", "hduclas2"); + + cout<<"FITS file created successfully"<GetCL(); + unsigned int numHits = Event->GetNHits(); + + // Event-level metadata (placeholders for now - can be filled in later) + uint8_t eventType = 0; // 0 = unknown/default + uint8_t eventClass = 0; // 0 = unknown (can check for Compton/photoabsorption later) + uint8_t seqHit = 0; // 0 = first/only sequence + + // Fixed-length arrays for event-level data (initialize to zeros) + std::valarray statTest(0.0f, 4); // 4 statistical test values + std::valarray recoilDir(0.0f, 3); // Recoil electron direction (x,y,z) + std::valarray recoilDirErr(0.0f, 3); // Recoil direction error + + // Resize arrays for this event's hits (using float to match PE format) + std::valarray x(numHits); + std::valarray y(numHits); + std::valarray z(numHits); + std::valarray x_err(numHits); + std::valarray y_err(numHits); + std::valarray z_err(numHits); + std::valarray energy(numHits); + std::valarray energy_err(numHits); + std::valarray bad_flag(0.0f, numHits); // Initialize flags to 0 + + // Extract hit-level data + for (unsigned int i = 0; i < numHits; ++i) { + MHit* hit = Event->GetHit(i); + + // Get position + MVector position = hit->GetPosition(); + x[i] = (float)position.X(); + y[i] = (float)position.Y(); + z[i] = (float)position.Z(); + + // Get position errors + MVector positionResolution = hit->GetPositionResolution(); + x_err[i] = (float)positionResolution.X(); + y_err[i] = (float)positionResolution.Y(); + z_err[i] = (float)positionResolution.Z(); + + // Get energy + energy[i] = (float)hit->GetEnergy(); + energy_err[i] = (float)hit->GetEnergyResolution(); + } + + // Add to batch + m_BatchTIME.push_back(time); + m_BatchEVENTTYPE.push_back(eventType); + m_BatchEVENTCLASS.push_back(eventClass); + m_BatchNUMHIT.push_back((uint8_t)numHits); + m_BatchSEQHIT.push_back(seqHit); + m_BatchSTATTEST.push_back(statTest); + m_BatchRECOILDIR.push_back(recoilDir); + m_BatchRECOILDIR_ERR.push_back(recoilDirErr); + m_BatchX.push_back(x); + m_BatchY.push_back(y); + m_BatchZ.push_back(z); + m_BatchX_ERR.push_back(x_err); + m_BatchY_ERR.push_back(y_err); + m_BatchZ_ERR.push_back(z_err); + m_BatchENERGY.push_back(energy); + m_BatchENERGY_ERR.push_back(energy_err); + m_BatchBAD_FLAG.push_back(bad_flag); + + m_BatchEventCount++; + + // Write batch if full + if (m_BatchEventCount >= m_BatchSize) { + if (FlushBatch() == false) { + m_IsOK = false; + return false; + } + } + + Event->SetAnalysisProgress(MAssembly::c_EventSaver); + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleSaverMeasurementsFITS::FlushBatch() +{ + // Write the current batch to the FITS file + + if (m_BatchEventCount == 0) { + return true; // Nothing to write + } + + try { + // Calculate last row for this batch + long lastRow = m_BatchStartRow + m_BatchEventCount - 1; + + if (g_Verbosity >= c_Info) { + cout<<"Writing batch: "<column("TIME").write(m_BatchTIME, m_BatchStartRow); + m_ScienceTable->column("EVENTTYPE").write(m_BatchEVENTTYPE, m_BatchStartRow); + m_ScienceTable->column("EVENTCLASS").write(m_BatchEVENTCLASS, m_BatchStartRow); + m_ScienceTable->column("NUMHIT").write(m_BatchNUMHIT, m_BatchStartRow); + m_ScienceTable->column("SEQHIT").write(m_BatchSEQHIT, m_BatchStartRow); + + // Write fixed-length array columns (event-level) + m_ScienceTable->column("STATTEST").writeArrays(m_BatchSTATTEST, m_BatchStartRow); + m_ScienceTable->column("RECOILDIR").writeArrays(m_BatchRECOILDIR, m_BatchStartRow); + m_ScienceTable->column("RECOILDIR_ERR").writeArrays(m_BatchRECOILDIR_ERR, m_BatchStartRow); + + // Write variable-length array columns (hit-level) + m_ScienceTable->column("X").writeArrays(m_BatchX, m_BatchStartRow); + m_ScienceTable->column("Y").writeArrays(m_BatchY, m_BatchStartRow); + m_ScienceTable->column("Z").writeArrays(m_BatchZ, m_BatchStartRow); + m_ScienceTable->column("X_ERR").writeArrays(m_BatchX_ERR, m_BatchStartRow); + m_ScienceTable->column("Y_ERR").writeArrays(m_BatchY_ERR, m_BatchStartRow); + m_ScienceTable->column("Z_ERR").writeArrays(m_BatchZ_ERR, m_BatchStartRow); + m_ScienceTable->column("ENERGY").writeArrays(m_BatchENERGY, m_BatchStartRow); + m_ScienceTable->column("ENERGY_ERR").writeArrays(m_BatchENERGY_ERR, m_BatchStartRow); + m_ScienceTable->column("BAD_FLAG").writeArrays(m_BatchBAD_FLAG, m_BatchStartRow); + + // Update tracking + m_TotalEventsWritten += m_BatchEventCount; + m_BatchStartRow += m_BatchEventCount; + + // Clear batch vectors - scalar columns + m_BatchTIME.clear(); + m_BatchEVENTTYPE.clear(); + m_BatchEVENTCLASS.clear(); + m_BatchNUMHIT.clear(); + m_BatchSEQHIT.clear(); + + // Clear batch vectors - fixed-length arrays + m_BatchSTATTEST.clear(); + m_BatchRECOILDIR.clear(); + m_BatchRECOILDIR_ERR.clear(); + + // Clear batch vectors - variable-length arrays + m_BatchX.clear(); + m_BatchY.clear(); + m_BatchZ.clear(); + m_BatchX_ERR.clear(); + m_BatchY_ERR.clear(); + m_BatchZ_ERR.clear(); + m_BatchENERGY.clear(); + m_BatchENERGY_ERR.clear(); + m_BatchBAD_FLAG.clear(); + + m_BatchEventCount = 0; + + return true; + + } catch (FitsException& e) { + cout<<"Error writing FITS batch: "< 0) { + FlushBatch(); + } + + MModule::Finalize(); + + cout<<"MModuleSaverMeasurementsFITS: "<GetNode("FileName"); + if (FileNameNode != nullptr) { + m_FileName = FileNameNode->GetValue(); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleSaverMeasurementsFITS::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + new MXmlNode(Node, "FileName", m_FileName); + + return Node; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleSaverMeasurementsFITS::ShowOptionsGUI() +{ + //! Show the options GUI + + MGUIOptionsSaverMeasurementsFITS* Options = new MGUIOptionsSaverMeasurementsFITS(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +// MModuleSaverMeasurementsFITS.cxx: the end... +////////////////////////////////////////////////////////////////////////////////