diff --git a/include/MFITSWriterL1a.h b/include/MFITSWriterL1a.h new file mode 100644 index 00000000..0b826d3c --- /dev/null +++ b/include/MFITSWriterL1a.h @@ -0,0 +1,99 @@ +/* + * MFITSWriterL1a.h + * + * Copyright (C) by Andreas Zoglauer, WingYeung Ma. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MFITSWriterL1a__ +#define __MFITSWriterL1a__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MString.h" + +// Nuclearizer libs: +#include "MReadOutAssembly.h" + +// CCfits libs: +#include + + +//////////////////////////////////////////////////////////////////////////////// + + +class MFITSWriterL1a +{ + // public interface: + public: + //! Default constructor + MFITSWriterL1a(); + //! Default destructor + virtual ~MFITSWriterL1a(); + + //! Create the FITS file + GED_L1A extension. Returns false on error. + bool Create(const MString& FileName); + //! Append one event's strip hits as a row. Returns false on error. + bool Write(MReadOutAssembly* Event); + //! Flush buffered rows to the table. Returns false on error. + bool FlushBatch(); + + void Close(); + + + private: + //! The FITS file object + CCfits::FITS* m_FITSFile; + //! The GED_L1A science table extension + CCfits::ExtHDU* m_Table; + + long m_BatchStartRow; + long m_BatchEventCount; + //! Batch size + static const long m_BatchSize = 100; + + bool m_HasEvents; + + //! First / last event time seen, RTS (mission seconds since 2025-01-01) + double m_FirstEventTime_RTS; + double m_LastEventTime_RTS; + + std::vector m_BatchTIME; + std::vector m_BatchEVENTID; + std::vector m_BatchEVENTTYPE; + std::vector m_BatchNUMSTRIPHIT; + + std::vector> m_BatchHITTYPE; // PB + std::vector> m_BatchDETID; // PB + std::vector> m_BatchSTRIPID; // PB + std::vector> m_BatchSIDEID; // PX (0/1) + std::vector> m_BatchFASTTIME; // PX (0/1) + std::vector> m_BatchPHA; // PI + std::vector> m_BatchTAC; // PI + + +#ifdef ___CLING___ + public: + ClassDef(MFITSWriterL1a, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MGUIOptionsLoaderMeasurementsL0.h b/include/MGUIOptionsLoaderMeasurementsL0.h new file mode 100644 index 00000000..5c96f52c --- /dev/null +++ b/include/MGUIOptionsLoaderMeasurementsL0.h @@ -0,0 +1,90 @@ +/* + * MGUIOptionsLoaderMeasurementsL0.h + * + * Copyright (C) by Andreas Zoglauer, WingYeung Ma. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MGUIOptionsLoaderMeasurementsL0__ +#define __MGUIOptionsLoaderMeasurementsL0__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// 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 L0 measurements loader +class MGUIOptionsLoaderMeasurementsL0 : public MGUIOptions +{ + // public Session: + public: + //! Default constructor + MGUIOptionsLoaderMeasurementsL0(MModule* Module); + //! Default destructor + virtual ~MGUIOptionsLoaderMeasurementsL0(); + + //! 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 the L0 binary file to load (.bin / .dat) + MGUIEFileSelector* m_FileSelectorL0; + + //! Select the strip map file + MGUIEFileSelector* m_FileSelectorStripMap; + + + +#ifdef ___CLING___ + public: + ClassDef(MGUIOptionsLoaderMeasurementsL0, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MGUIOptionsSaverMeasurementsL0.h b/include/MGUIOptionsSaverMeasurementsL0.h index 7dc9bd71..2539e5bd 100644 --- a/include/MGUIOptionsSaverMeasurementsL0.h +++ b/include/MGUIOptionsSaverMeasurementsL0.h @@ -72,6 +72,8 @@ class MGUIOptionsSaverMeasurementsL0 : public MGUIOptions //! Select which output file to save to MGUIEFileSelector* m_FileSelectorOutput; + MGUIEFileSelector* m_FileSelectorStripMap; + #ifdef ___CLING___ diff --git a/include/MModuleLoaderMeasurementsFITS.h b/include/MModuleLoaderMeasurementsFITS.h index e5a2caeb..22e3b963 100644 --- a/include/MModuleLoaderMeasurementsFITS.h +++ b/include/MModuleLoaderMeasurementsFITS.h @@ -113,11 +113,12 @@ class MModuleLoaderMeasurementsFITS : public MModuleLoaderMeasurements //! Batch data storage for scalar columns std::vector m_BatchTIME; + std::vector m_BatchEVENTID; std::vector m_BatchEVENTTYPE; std::vector m_BatchNUMSTRIPHIT; //! Batch data storage for variable-length array columns - std::vector> m_BatchTYPEHIT; + std::vector> m_BatchHITTYPE; std::vector> m_BatchDETID; std::vector> m_BatchSTRIPID; std::vector> m_BatchSIDEID; diff --git a/include/MModuleLoaderMeasurementsL0.h b/include/MModuleLoaderMeasurementsL0.h new file mode 100644 index 00000000..15ef61a9 --- /dev/null +++ b/include/MModuleLoaderMeasurementsL0.h @@ -0,0 +1,122 @@ +/* + * MModuleLoaderMeasurementsL0.h + * + * Copyright (C) by Andreas Zoglauer, WingYeung Ma. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleLoaderMeasurementsL0__ +#define __MModuleLoaderMeasurementsL0__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MString.h" + +// Nuclearizer libs: +#include "MModuleLoaderMeasurements.h" +#include "MStripMap.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +//! A module to load L0 binary files (CCSDS DD packets) into MReadOutAssembly events. +class MModuleLoaderMeasurementsL0 : public MModuleLoaderMeasurements +{ + // public interface: + public: + //! Default constructor + MModuleLoaderMeasurementsL0(); + //! Default destructor + virtual ~MModuleLoaderMeasurementsL0(); + + //! Create a new object of this class + virtual MModuleLoaderMeasurementsL0* Clone() { return new MModuleLoaderMeasurementsL0(); } + + //! 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 L0 binary file name + void SetFileName(const MString& FileName) { m_FileName = FileName; } + //! Get the L0 binary file name + MString GetFileName() const { return m_FileName; } + + //! Set the strip map file name + void SetFileNameStripMap(const MString& FileName) { m_FileNameStripMap = FileName; } + //! Get the strip map file name + MString GetFileNameStripMap() const { return m_FileNameStripMap; } + + + // protected methods: + protected: + //! Open the L0 file. Auto-detects whether it has a 20-byte L0 header (from MOC simulator) or is raw CCSDS packets. + bool OpenL0File(MString FileName); + + //! Read the next CCSDS packet from the file and populate Event. + //! Returns false on EOF or unrecoverable read error. + bool ReadNextPacket(MReadOutAssembly* Event); + + //! Decode the bit-packed HIT_DATA payload into MStripHit objects on Event. + //! Returns false on parsing error. + bool DecodeHitData(const std::vector& hitData, unsigned int expectedHits, MReadOutAssembly* Event); + + + // private members: + private: + //! Input L0 binary file name (.dat or .bin) + MString m_FileName; + + //! Strip map file name for read-out ID -> (detector, side, strip) lookup + MString m_FileNameStripMap; + + //! Strip map (loaded at Initialize) + MStripMap m_StripMap; + + //! Whether the strip map was successfully loaded + bool m_StripMapLoaded; + + //! The opened input file stream + std::ifstream m_InFile; + + //! Total number of packets read so far + unsigned long m_PacketsRead; + + +#ifdef ___CLING___ + public: + ClassDef(MModuleLoaderMeasurementsL0, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MModuleSaverMeasurementsFITS.h b/include/MModuleSaverMeasurementsFITS.h index 589579d6..bf860c9c 100644 --- a/include/MModuleSaverMeasurementsFITS.h +++ b/include/MModuleSaverMeasurementsFITS.h @@ -30,6 +30,7 @@ // Nuclearizer libs: #include "MModule.h" +#include "MFITSWriterL1a.h" // CCfits libs #include @@ -101,9 +102,12 @@ class MModuleSaverMeasurementsFITS : public MModule //! Output file name MString m_FileName; - //! Output data level: 1 = L1b (all events, with QUALITY_FLAG), 2 = L2 (screened, no QUALITY_FLAG) + //! Output data level: 0 = L1a (raw hits), 1 = L1b (all events, with QUALITY_FLAG), 2 = L2 (screened, no QUALITY_FLAG) int m_OutputDataLevel; + //! L1a writer, used only when m_OutputDataLevel == 0 + MFITSWriterL1a m_L1aWriter; + //! The FITS file object pointer FITS* m_FITSFile; diff --git a/include/MModuleSaverMeasurementsL0.h b/include/MModuleSaverMeasurementsL0.h index 1684da22..ebd5dc1e 100644 --- a/include/MModuleSaverMeasurementsL0.h +++ b/include/MModuleSaverMeasurementsL0.h @@ -29,6 +29,7 @@ // Nuclearizer libs: #include "MModule.h" #include "MReadOutAssembly.h" +#include "MStripMap.h" //////////////////////////////////////////////////////////////////////////////// @@ -69,6 +70,11 @@ class MModuleSaverMeasurementsL0 : public MModule //! Get the output file name MString GetFileName() const { return m_FileName; } + //! Set the strip map file name (used for reverse (det,side,strip) → ReadOutID lookup) + void SetFileNameStripMap(const MString& FileName) { m_FileNameStripMap = FileName; } + //! Get the strip map file name + MString GetFileNameStripMap() const { return m_FileNameStripMap; } + // protected methods: protected: //! Write CCSDS Primary Header @@ -104,6 +110,14 @@ class MModuleSaverMeasurementsL0 : public MModule //! Output file name MString m_FileName; + //! Strip map file name. Used for reverse lookup (det,side,strip) → read-out ID + MString m_FileNameStripMap; + + //! The loaded strip map + MStripMap m_StripMap; + + bool m_StripMapLoaded; + //! Output file stream std::ofstream m_OutFile; diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index 9d5005f5..8e46ee70 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -279,6 +279,10 @@ class MReadOutAssembly : public MReadOutSequence MTime ComputeRTSfromUTCTime(MTime UTCTime); //! Compute the UTC time from known RTS MTime ComputeUTCfromRTSTime(MTime RTSTime); + //! Compute the RTS time from GPS time + MTime ComputeRTSfromGPSTime(MTime GPSTime); + //! Compute GPS time from known RTS + MTime ComputeGPSfromRTSTime(MTime RTSTime); // protected methods: protected: diff --git a/include/MStripMap.h b/include/MStripMap.h index 821def58..b7e3f8f0 100644 --- a/include/MStripMap.h +++ b/include/MStripMap.h @@ -19,6 +19,7 @@ // Standard libs: #include #include +#include using namespace std; // ROOT libs: @@ -60,6 +61,12 @@ class MStripMap //! Get strip ID by read out ID - check with HasReadOutID(ROI) first unsigned int GetStripNumber(unsigned int ROI) const; + //! Check if a (detector, side, strip) tuple is mapped to a read-out ID + bool HasROIDetSideStrip(unsigned int DetectorID, bool IsLowVoltage, unsigned int StripNumber) const; + + //! Reverse lookup: get read-out ID for a (detector, side, strip) tuple. + unsigned int GetReadOutID(unsigned int DetectorID, bool IsLowVoltage, unsigned int StripNumber) const; + // protected methods: protected: @@ -94,6 +101,10 @@ class MStripMap //! The strip mapping data vector m_StripMappings; + //! Reverse index: (det<<8) | (side<<7) | strip → read-out ID. 12-bit key + //! det 0-15 (4 bits), side 0/1 (1 bit), strip 0-64 (7 bits - guard-ring strip 64) + unordered_map m_DetSideStripToROI; + #ifdef ___CLING___ public: diff --git a/src/MAssembly.cxx b/src/MAssembly.cxx index 70340d41..29af9d41 100644 --- a/src/MAssembly.cxx +++ b/src/MAssembly.cxx @@ -65,6 +65,7 @@ using namespace std; #include "MModuleLoaderMeasurementsROA.h" #include "MModuleLoaderMeasurementsHDF.h" #include "MModuleLoaderMeasurementsFITS.h" +#include "MModuleLoaderMeasurementsL0.h" #include "MModuleEnergyCalibration.h" #include "MModuleDepthCalibration.h" #include "MModuleStripPairingMultiRoundChiSquare.h" @@ -128,6 +129,7 @@ MAssembly::MAssembly() m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsROA()); m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsHDF()); m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsFITS()); + m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsL0()); m_Supervisor->AddAvailableModule(new MModuleDEESMEX()); diff --git a/src/MFITSWriterL1a.cxx b/src/MFITSWriterL1a.cxx new file mode 100644 index 00000000..ee0e448d --- /dev/null +++ b/src/MFITSWriterL1a.cxx @@ -0,0 +1,339 @@ +/* + * MFITSWriterL1a.cxx + * + * + * Copyright (C) by Andreas Zoglauer, WingYeung Ma. + * 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. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MFITSWriterL1a +// +// Writes a GED_L1A (Level 1a, raw strip hits) FITS file. Owned by +// MModuleSaverMeasurementsFITS; not a pipeline module itself. +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MFITSWriterL1a.h" + +// Standard libs: + +// MEGAlib libs: +#include "MStripHit.h" +#include "MStreams.h" + +using namespace std; +using namespace CCfits; + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MFITSWriterL1a) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MFITSWriterL1a::MFITSWriterL1a() +{ + // Construct an instance of MFITSWriterL1a + + m_FITSFile = nullptr; + m_Table = nullptr; + m_BatchStartRow = 1; + m_BatchEventCount = 0; + m_HasEvents = false; + m_FirstEventTime_RTS = 0; + m_LastEventTime_RTS = 0; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MFITSWriterL1a::~MFITSWriterL1a() +{ + delete m_FITSFile; + m_FITSFile = nullptr; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MFITSWriterL1a::Create(const MString& FileName) +{ + // Create the FITS file with an empty primary HDU and the GED_L1A extension. + try { + m_FITSFile = new FITS(string(FileName), RWmode::Write); + + // Primary header (Table 6.2d primary) + PHDU& primary = m_FITSFile->pHDU(); + primary.addKey("TELESCOP", "COSI", "Telescope mission name"); + primary.addKey("INSTRUME", "GED", "Instrument name"); + primary.addKey("ORIGIN", "SSL", "Origin of the FITS file"); + primary.addKey("CREATOR", "TBD", "Software that create 1st the file"); + + // File creation date (UTC, computed at write time) + time_t now = time(nullptr); + char dateBuffer[32]; + strftime(dateBuffer, sizeof(dateBuffer), "%Y-%m-%dT%H:%M:%S", gmtime(&now)); + primary.addKey("DATE", string(dateBuffer), "File creation date (UTC)"); + + // add placeholder at creation, write at the end + primary.addKey("OBS_ID", "YYYYMMDD", "Observation ID"); + primary.addKey("DATE-OBS", "yyyy-mm-ddThh:mm:ss", "Start Date"); + primary.addKey("DATE-END", "yyyy-mm-ddThh:mm:ss", "Stop Date"); + + // GED_L1A binary table: 11 columns. P* = variable-length arrays (max 2080). + struct ColSpec { string name, format, unit, comment; }; + const std::vector colsGedL1a = { + {"TIME", "1D", "s", "Mission Time in sec since 01 Jan 2025 00:00:00"}, + {"EVENTID", "1J", "", "Event ID, unique number that starts at 1 each day"}, + {"EVENTTYPE", "1B", "", "Type of event"}, + {"NUMSTRIPHIT", "1I", "", "Number of strips involved in the hit"}, + {"HITTYPE", "PB(2080)", "", "Type of Hits"}, + {"DETID", "PB(2080)", "", "DETECTOR ID 0-15"}, + {"STRIPID", "PB(2080)", "", "Strip ID Number 0-64"}, + {"SIDEID", "PX(2080)", "", "Strip Side ID 0-1"}, + {"FASTTIME", "PX(2080)", "", "Interaction fast timing value 0 or 1"}, + {"PHA", "PI(2080)", "chan", "PHA ADC value for each strip hit"}, + {"TAC", "PI(2080)", "chan", "Timing analog converter for each strip hit"} + }; + + std::vector colNames, colFormats, colUnits; + colNames.reserve(colsGedL1a.size()); + colFormats.reserve(colsGedL1a.size()); + colUnits.reserve(colsGedL1a.size()); + for (const auto& c : colsGedL1a) { + colNames.push_back(c.name); + colFormats.push_back(c.format); + colUnits.push_back(c.unit); + } + + m_Table = m_FITSFile->addTable("GED_L1A", 0, colNames, colFormats, colUnits); + + { + m_Table->makeThisCurrent(); + int status = 0; + for (size_t i = 0; i < colsGedL1a.size(); ++i) { + char keyName[16]; + snprintf(keyName, sizeof(keyName), "TTYPE%zu", i + 1); + fits_modify_comment(m_FITSFile->fitsPointer(), + keyName, + const_cast(colsGedL1a[i].comment.c_str()), + &status); + } + if (status != 0 && g_Verbosity >= c_Warning) { + cout << "MFITSWriterL1a: fits_modify_comment status=" << status << endl; + } + } + + // Extension header keywords (Table 6.2d) + m_Table->addKey("EXTNAME", "GED_L1A", "name of this HDU"); + m_Table->addKey("TELESCOP", "COSI", "Telescope mission name"); + m_Table->addKey("INSTRUME", "GED", "Instrument name"); + m_Table->addKey("DATAMODE", "SYNC", "Instrument datamode"); + m_Table->addKey("OBSMODE", "SCANNING", "Spacecraft observing mode"); + m_Table->addKey("OBSERVER", "John Tomsick", "Principal Investigator"); + m_Table->addKey("OBJECT", "ALLSKY", "Object/target name of ALLSKY"); + m_Table->addKey("MJDREFI", 60676, "MJD reference day 01 Jan 2025 00:00:00"); + m_Table->addKey("MJDREFF", 8.007407407407E-04, "MJD reference (fraction of day)"); + m_Table->addKey("TIMEREF", "LOCAL", "Reference Frame"); + m_Table->addKey("TASSIGN", "SATELLITE", "Time assigned"); + m_Table->addKey("TIMESYS", "TT", "Time System"); + m_Table->addKey("TIMEUNIT", "s", "Time unit for timing header keywords"); + m_Table->addKey("CLOCKAPP", "F", "If clock corrections are applied (T/F)"); + m_Table->addKey("HDUCLASS", "OGIP", "format conforms to OGIP standard"); + m_Table->addKey("HDUCLAS1", "EVENTS", "hduclass1"); + m_Table->addKey("HDUCLAS2", "ALL", "hduclas2"); + m_Table->addKey("CREATOR", "TBD", "Software that create 1st the file"); + m_Table->addKey("PROCVER", "MM.XX.NN.YY", "Processing version"); + m_Table->addKey("CALDBVER", "csYYYYMMDD", "CALDB version"); + m_Table->addKey("SEQPNUM", "nn", "Times the dataset has been processed"); + m_Table->addKey("ORIGIN", "SSL", "Origin of the FITS files"); + + m_Table->addKey("DATE", string(dateBuffer), "File creation date (UTC)"); + m_Table->addKey("OBS_ID", "YYYYMMDD", "Observation ID"); + m_Table->addKey("DATE-OBS", "yyyy-mm-ddThh:mm:ss", "Start Date"); + m_Table->addKey("DATE-END", "yyyy-mm-ddThh:mm:ss", "Stop Date"); + + m_BatchStartRow = 1; + m_BatchEventCount = 0; + m_HasEvents = false; + return true; + } catch (const FitsException& e) { + if (g_Verbosity >= c_Error) cout << "MFITSWriterL1a: Create error: " << e.message() << endl; + return false; + } +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MFITSWriterL1a::Write(MReadOutAssembly* Event) +{ + + double time = Event->GetTimeRTS().GetAsDouble(); + if (!m_HasEvents) { + m_FirstEventTime_RTS = time; + m_LastEventTime_RTS = time; + m_HasEvents = true; + } else { + if (time < m_FirstEventTime_RTS) m_FirstEventTime_RTS = time; + if (time > m_LastEventTime_RTS) m_LastEventTime_RTS = time; + } + + unsigned int nStripHits = Event->GetNStripHits(); + + std::valarray hittype(nStripHits), detid(nStripHits), stripid(nStripHits), sideid(nStripHits), fasttime(nStripHits); + std::valarray pha(nStripHits), tac(nStripHits); + + for (unsigned int i = 0; i < nStripHits; ++i) { + MStripHit* hit = Event->GetStripHit(i); + uint8_t hitType = 0; + if (hit->IsGuardRing()) hitType = 2; + else if (hit->IsNearestNeighbor()) hitType = 1; + hittype[i] = hitType; + detid[i] = hit->GetDetectorID(); + stripid[i] = hit->GetStripID(); + sideid[i] = hit->IsLowVoltageStrip() ? 0 : 1; + fasttime[i] = hit->HasFastTiming() ? 1 : 0; + pha[i] = hit->GetADCUnits(); + tac[i] = hit->GetTAC(); + } + + m_BatchTIME.push_back(time); + m_BatchEVENTID.push_back((uint32_t)Event->GetID()); + m_BatchEVENTTYPE.push_back(0); // TODO + m_BatchNUMSTRIPHIT.push_back(nStripHits); + m_BatchHITTYPE.push_back(hittype); + m_BatchDETID.push_back(detid); + m_BatchSTRIPID.push_back(stripid); + m_BatchSIDEID.push_back(sideid); + m_BatchFASTTIME.push_back(fasttime); + m_BatchPHA.push_back(pha); + m_BatchTAC.push_back(tac); + m_BatchEventCount++; + + if (m_BatchEventCount >= m_BatchSize) { + return FlushBatch(); + } + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MFITSWriterL1a::FlushBatch() +{ + if (m_BatchEventCount == 0) return true; + try { + long firstRow = m_BatchStartRow; + m_Table->column("TIME").write(m_BatchTIME, firstRow); + m_Table->column("EVENTID").write(m_BatchEVENTID, firstRow); + m_Table->column("EVENTTYPE").write(m_BatchEVENTTYPE, firstRow); + m_Table->column("NUMSTRIPHIT").write(m_BatchNUMSTRIPHIT, firstRow); + + // Variable-length array columns + m_Table->column("HITTYPE").writeArrays(m_BatchHITTYPE, firstRow); + m_Table->column("DETID").writeArrays(m_BatchDETID, firstRow); + m_Table->column("STRIPID").writeArrays(m_BatchSTRIPID, firstRow); + m_Table->column("SIDEID").writeArrays(m_BatchSIDEID, firstRow); + m_Table->column("FASTTIME").writeArrays(m_BatchFASTTIME, firstRow); + m_Table->column("PHA").writeArrays(m_BatchPHA, firstRow); + m_Table->column("TAC").writeArrays(m_BatchTAC, firstRow); + + m_BatchStartRow += m_BatchEventCount; + + m_BatchTIME.clear(); + m_BatchEVENTID.clear(); + m_BatchEVENTTYPE.clear(); + m_BatchNUMSTRIPHIT.clear(); + m_BatchHITTYPE.clear(); + m_BatchDETID.clear(); + m_BatchSTRIPID.clear(); + m_BatchSIDEID.clear(); + m_BatchFASTTIME.clear(); + m_BatchPHA.clear(); + m_BatchTAC.clear(); + m_BatchEventCount = 0; + return true; + } catch (const FitsException& e) { + if (g_Verbosity >= c_Error) cout << "MFITSWriterL1a: FlushBatch error: " << e.message() << endl; + return false; + } +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MFITSWriterL1a::Close() +{ + if (m_BatchEventCount > 0) FlushBatch(); + + if (m_HasEvents && m_Table != nullptr) { + try { + m_Table->addKey("TSTART", m_FirstEventTime_RTS, "Start time [s] since MJDREFI"); + m_Table->addKey("TSTOP", m_LastEventTime_RTS, "Stop time [s] since MJDREFI"); + + constexpr time_t MISSION_EPOCH_UNIX = 1735689600; // 2025-01-01T00:00:00 UTC + time_t startUnix = MISSION_EPOCH_UNIX + (time_t)m_FirstEventTime_RTS; + time_t stopUnix = MISSION_EPOCH_UNIX + (time_t)m_LastEventTime_RTS; + char startBuf[32], stopBuf[32]; + strftime(startBuf, sizeof(startBuf), "%Y-%m-%dT%H:%M:%S", gmtime(&startUnix)); + strftime(stopBuf, sizeof(stopBuf), "%Y-%m-%dT%H:%M:%S", gmtime(&stopUnix)); + + PHDU& primary = m_FITSFile->pHDU(); + primary.addKey("DATE-OBS", string(startBuf), "Start Date"); + primary.addKey("DATE-END", string(stopBuf), "Stop Date"); + m_Table->addKey("DATE-OBS", string(startBuf), "Start Date"); + m_Table->addKey("DATE-END", string(stopBuf), "Stop Date"); + + char obsIdBuf[16]; + strftime(obsIdBuf, sizeof(obsIdBuf), "%Y%m%d", gmtime(&startUnix)); + primary.addKey("OBS_ID", string(obsIdBuf), "Observation ID"); + m_Table->addKey("OBS_ID", string(obsIdBuf), "Observation ID"); + } catch (const FitsException& e) { + if (g_Verbosity >= c_Error) cout << "MFITSWriterL1a: Close keyword error: " << e.message() << endl; + } + } + + if (m_FITSFile != nullptr) { + try { + m_FITSFile->pHDU().writeChecksum(); + if (m_Table != nullptr) { + m_Table->writeChecksum(); + } + } catch (const FitsException& e) { + if (g_Verbosity >= c_Error) cout << "MFITSWriterL1a: writeChecksum error: " << e.message() << endl; + } + } + + delete m_FITSFile; // CCfits flushes + closes on destruction + m_FITSFile = nullptr; + m_Table = nullptr; +} + + +// MFITSWriterL1a.cxx: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MGUIOptionsLoaderMeasurementsL0.cxx b/src/MGUIOptionsLoaderMeasurementsL0.cxx new file mode 100644 index 00000000..b624a6c3 --- /dev/null +++ b/src/MGUIOptionsLoaderMeasurementsL0.cxx @@ -0,0 +1,132 @@ +/* + * MGUIOptionsLoaderMeasurementsL0.cxx + * + * + * Copyright (C) by Andreas Zoglauer, WingYeung Ma. + * 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 "MGUIOptionsLoaderMeasurementsL0.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" +#include "MModuleLoaderMeasurementsL0.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIOptionsLoaderMeasurementsL0) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsLoaderMeasurementsL0::MGUIOptionsLoaderMeasurementsL0(MModule* Module) + : MGUIOptions(Module) +{ + // standard constructor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsLoaderMeasurementsL0::~MGUIOptionsLoaderMeasurementsL0() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIOptionsLoaderMeasurementsL0::Create() +{ + PreCreate(); + + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + + m_FileSelectorL0 = new MGUIEFileSelector(m_OptionsFrame, "Please select an L0 binary file:", + dynamic_cast(m_Module)->GetFileName()); + m_FileSelectorL0->SetFileType("L0 binary file", "*.bin"); + m_FileSelectorL0->SetFileType("L0 binary file", "*.dat"); + m_OptionsFrame->AddFrame(m_FileSelectorL0, LabelLayout); + + m_FileSelectorStripMap = new MGUIEFileSelector(m_OptionsFrame, "Please select a strip map file:", + dynamic_cast(m_Module)->GetFileNameStripMap()); + m_FileSelectorStripMap->SetFileType("Strip map file", "*.map"); + m_OptionsFrame->AddFrame(m_FileSelectorStripMap, LabelLayout); + + PostCreate(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsLoaderMeasurementsL0::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 MGUIOptionsLoaderMeasurementsL0::OnApply() +{ + // Modify this to store the data in the module! + + dynamic_cast(m_Module)->SetFileName(m_FileSelectorL0->GetFileName()); + dynamic_cast(m_Module)->SetFileNameStripMap(m_FileSelectorStripMap->GetFileName()); + + return true; +} + + +// MGUIOptionsLoaderMeasurementsL0: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MGUIOptionsSaverMeasurementsFITS.cxx b/src/MGUIOptionsSaverMeasurementsFITS.cxx index 7b096bb5..2160ccb6 100644 --- a/src/MGUIOptionsSaverMeasurementsFITS.cxx +++ b/src/MGUIOptionsSaverMeasurementsFITS.cxx @@ -79,6 +79,7 @@ void MGUIOptionsSaverMeasurementsFITS::Create() m_OptionsFrame->AddFrame(LevelLabel, LabelLayout); m_OutputDataLevelCombo = new TGComboBox(m_OptionsFrame); + m_OutputDataLevelCombo->AddEntry("L1a (raw hits, no calibration)", 0); m_OutputDataLevelCombo->AddEntry("L1b (all events, with QUALITY_FLAG)", 1); m_OutputDataLevelCombo->AddEntry("L2 (screened, no QUALITY_FLAG)", 2); m_OutputDataLevelCombo->Select(dynamic_cast(m_Module)->GetOutputDataLevel()); diff --git a/src/MGUIOptionsSaverMeasurementsL0.cxx b/src/MGUIOptionsSaverMeasurementsL0.cxx index 3b389d4c..ab2a9222 100644 --- a/src/MGUIOptionsSaverMeasurementsL0.cxx +++ b/src/MGUIOptionsSaverMeasurementsL0.cxx @@ -74,6 +74,11 @@ void MGUIOptionsSaverMeasurementsL0::Create() m_FileSelectorOutput->SetFileType("Binary file", "*.bin"); m_OptionsFrame->AddFrame(m_FileSelectorOutput, LabelLayout); + m_FileSelectorStripMap = new MGUIEFileSelector(m_OptionsFrame, "Please select strip map file (required):", + dynamic_cast(m_Module)->GetFileNameStripMap()); + m_FileSelectorStripMap->SetFileType("Strip map file", "*.map"); + m_OptionsFrame->AddFrame(m_FileSelectorStripMap, LabelLayout); + PostCreate(); } @@ -117,6 +122,7 @@ bool MGUIOptionsSaverMeasurementsL0::OnApply() // Modify this to store the data in the module! dynamic_cast(m_Module)->SetFileName(m_FileSelectorOutput->GetFileName()); + dynamic_cast(m_Module)->SetFileNameStripMap(m_FileSelectorStripMap->GetFileName()); return true; } diff --git a/src/MModuleLoaderMeasurementsFITS.cxx b/src/MModuleLoaderMeasurementsFITS.cxx index 2c2fa969..284c9fc7 100644 --- a/src/MModuleLoaderMeasurementsFITS.cxx +++ b/src/MModuleLoaderMeasurementsFITS.cxx @@ -167,7 +167,7 @@ bool MModuleLoaderMeasurementsFITS::OpenFITSFile(MString FileName) // Validate required columns exist const vector requiredColumns = { - "TIME", "EVENTTYPE", "NUMSTRIPHIT", "TYPEHIT", "DETID", + "TIME", "EVENTID", "EVENTTYPE", "NUMSTRIPHIT", "HITTYPE", "DETID", "STRIPID", "SIDEID", "FASTTIME", "PHA", "TAC" }; @@ -242,9 +242,10 @@ bool MModuleLoaderMeasurementsFITS::ReadBatch() // Resize vectors to read this batch m_BatchTIME.resize(rowsToRead); + m_BatchEVENTID.resize(rowsToRead); m_BatchEVENTTYPE.resize(rowsToRead); m_BatchNUMSTRIPHIT.resize(rowsToRead); - m_BatchTYPEHIT.resize(rowsToRead); + m_BatchHITTYPE.resize(rowsToRead); m_BatchDETID.resize(rowsToRead); m_BatchSTRIPID.resize(rowsToRead); m_BatchSIDEID.resize(rowsToRead); @@ -254,11 +255,12 @@ bool MModuleLoaderMeasurementsFITS::ReadBatch() // Read scalar columns - third parameter is LAST row index, NOT count m_ComptonTable->column("TIME").read(m_BatchTIME, m_BatchStartRow, lastRow); + m_ComptonTable->column("EVENTID").read(m_BatchEVENTID, 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("HITTYPE").readArrays(m_BatchHITTYPE, 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); @@ -278,7 +280,7 @@ bool MModuleLoaderMeasurementsFITS::ReadBatch() // <<"TIME="< +using namespace std; + +// ROOT libs: +#include "TGClient.h" + +// MEGAlib libs: +#include "MFile.h" +#include "MStripHit.h" +#include "MTime.h" + +// Nuclearizer libs: +#include "MGUIOptionsLoaderMeasurementsL0.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleLoaderMeasurementsL0) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +// CCSDS / L0 file constants +static const uint32_t COSI_IDENTIFIER = 0x434F5349; // "COSI" in ASCII +static const size_t L0_FILE_HEADER_SIZE = 20; // 20 byte: 4 bytes ID, 2 bytes SC ID, 2 bytes PKT APID, 4 bytes UTC first pkt, 4 bytes UTC time last pkt, 4 bytes pkt count +static const size_t L0_CRC_TRAILER_SIZE = 2; //16-bit CRC-CCITT checksum of the file +static const size_t CCSDS_PRIMARY_HEADER_SIZE = 6; +static const size_t CCSDS_SECONDARY_HEADER_SIZE = 8; +static const uint16_t APID_DD = 0x0DD; + + +//////////////////////////////////////////////////////////////////////////////// + + +// bit-level reader for the variable-length HIT_DATA payload. +namespace { +class BitReader { +public: + BitReader(const std::vector& data) : m_Data(data), m_BytePos(0), m_BitPos(0) {} + + // Read up to 32 bits and return them as an unsigned integer + uint32_t ReadBits(unsigned int n) { + uint32_t result = 0; + for (unsigned int i = 0; i < n; ++i) { + if (m_BytePos >= m_Data.size()) { + throw std::out_of_range("BitReader: ran out of bits"); + } + // read each bit and accumulate into result + uint8_t bit = (m_Data[m_BytePos] >> (7 - m_BitPos)) & 0x1; + result = (result << 1) | bit; + + // Advance the cursor. Move to the next bit + ++m_BitPos; + if (m_BitPos == 8) { + ++m_BytePos; + m_BitPos = 0; + } + } + return result; + } + + size_t BitsRemaining() const { + if (m_BytePos >= m_Data.size()) return 0; + return (m_Data.size() - m_BytePos) * 8 - m_BitPos; + } + +private: + const std::vector& m_Data; + size_t m_BytePos; + unsigned int m_BitPos; +}; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleLoaderMeasurementsL0::MModuleLoaderMeasurementsL0() : MModuleLoaderMeasurements() +{ + // Construct an instance of MModuleLoaderMeasurementsL0 + + m_Name = "Measurement loader for L0 binary files (DD packets)"; + m_XmlTag = "XmlTagMeasurementLoaderL0"; + + m_IsStartModule = true; + + m_AllowMultiThreading = true; + m_AllowMultipleInstances = false; + + m_StripMapLoaded = false; + m_PacketsRead = 0; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleLoaderMeasurementsL0::~MModuleLoaderMeasurementsL0() +{ + // Delete this instance of MModuleLoaderMeasurementsL0 +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleLoaderMeasurementsL0::Initialize() +{ + m_PacketsRead = 0; + + if (m_FileName == "") { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout<(magicBuf), 4); + if (m_InFile.gcount() != 4) { + if (g_Verbosity >= c_Error) cout<= c_Info) cout<= c_Info) cout<= c_Info) cout<SetAnalysisProgress(MAssembly::c_EventLoader | MAssembly::c_EventLoaderMeasurement); + + ++m_PacketsRead; + + // set the current packet # as the event ID + Event->SetID(m_PacketsRead); + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleLoaderMeasurementsL0::ReadNextPacket(MReadOutAssembly* Event) +{ + // Read the 6-byte CCSDS primary header + uint8_t priHdr[CCSDS_PRIMARY_HEADER_SIZE]; + m_InFile.read(reinterpret_cast(priHdr), CCSDS_PRIMARY_HEADER_SIZE); + if (m_InFile.gcount() < (streamsize)CCSDS_PRIMARY_HEADER_SIZE) { + return false; // EOF or truncated + } + + uint16_t word0 = (uint16_t(priHdr[0]) << 8) | priHdr[1]; + uint16_t pktDataLen = (uint16_t(priHdr[4]) << 8) | priHdr[5]; + uint16_t apid = word0 & 0x07FF; + bool secHdrFlag = (word0 >> 11) & 0x1; + + // Total bytes after primary header = pktDataLen + 1 + size_t payloadSize = pktDataLen + 1; + + if (apid != APID_DD) { + if (g_Verbosity >= c_Error) { + cout << m_XmlTag << ": Unexpected APID 0x" << hex << apid << dec + << " in file (expected 0x" << hex << APID_DD << dec << "). Stopping." << endl; + } + m_IsFinished = true; + return false; + } + + // Read the rest of the packet (secondary header + payload) + std::vector rest(payloadSize); + m_InFile.read(reinterpret_cast(rest.data()), payloadSize); + if (m_InFile.gcount() < (streamsize)payloadSize) { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<> 2; + unsigned int hlen = ((uint16_t(byte0 & 0x3) << 8) | byte1); + cursor += 3; + + if (rest.size() < cursor + hlen) { + if (g_Verbosity >= c_Error) cout< hitData(rest.begin() + cursor, rest.begin() + cursor + hlen); + + // Build GPS MTime from seconds and subseconds, then convert and set to RTS + // Subseconds are 40 MHz DCB ticks → ns = ticks * 25 + long int pktNanoseconds = (long int)pktSubseconds * 25; + MTime gpsTime((long int)pktSeconds, pktNanoseconds); + Event->SetTimeRTS(Event->ComputeRTSfromGPSTime(gpsTime)); + + // Decode the bit-packed HIT_DATA into individual MStripHit objects on the event + if (DecodeHitData(hitData, hits, Event) == false) { + return false; + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleLoaderMeasurementsL0::DecodeHitData(const std::vector& hitData, + unsigned int expectedHits, + MReadOutAssembly* Event) +{ + BitReader br(hitData); + unsigned int parsed = 0; + + while (br.BitsRemaining() >= 4 && parsed < expectedHits) { + uint32_t hitType = 0; + try { + hitType = br.ReadBits(4); + } catch (const std::out_of_range&) { + break; + } + + unsigned int stripID = 0; + int adc = 0; + int tac = 0; + bool fastTiming = false; + bool isGuardRing = false; + bool isNeighbor = false; + + try { + if (hitType == 0x0) { + // Normal: 11(strip) + 1(fastTimingFlag) + 14(adc) + 14(tac) = 40 more bits + if (br.BitsRemaining() < 40) break; + stripID = br.ReadBits(11); + fastTiming = br.ReadBits(1) == 1; + adc = br.ReadBits(14); + tac = br.ReadBits(14); + } else if (hitType == 0x1) { + // Neighbor: 11 + 1 + 10 + 10 = 32 more bits + if (br.BitsRemaining() < 32) break; + stripID = br.ReadBits(11); + fastTiming = br.ReadBits(1) == 1; + adc = br.ReadBits(10); + tac = br.ReadBits(10); + isNeighbor = true; + } else if (hitType == 0x2) { + // Guard ring: 5 + 14 + 1(pad) = 20 more bits + if (br.BitsRemaining() < 20) break; + stripID = br.ReadBits(5); + adc = br.ReadBits(14); + (void)br.ReadBits(1); // pad + isGuardRing = true; + tac = 0; + fastTiming = false; + } else if (hitType == 0x3) { + // Test pulser: 11(strip) + 1(fastTimingFlag) + 14(adc) + 14(tac) = 40 more bits + if (br.BitsRemaining() < 40) break; + stripID = br.ReadBits(11); + fastTiming = br.ReadBits(1) == 1; + adc = br.ReadBits(14); + tac = br.ReadBits(14); + } else { + // Unknown — abort this packet + if (g_Verbosity >= c_Warning) { + cout<= c_Warning) { + cout<SetDetectorID(detectorID); + SH->SetStripID(stripNumber); + SH->IsLowVoltageStrip(isLowVoltage); + SH->SetADCUnits(adc); + SH->SetTAC(tac); + SH->HasFastTiming(fastTiming); + SH->IsNearestNeighbor(isNeighbor); + SH->IsGuardRing(isGuardRing); + + if (isGuardRing) { + Event->SetGuardRingVeto(true); + } + + Event->AddStripHit(SH); + ++parsed; + } + + if (parsed != expectedHits && g_Verbosity >= c_Warning) { + cout<= c_Info) { + cout<GetNode("FileName"); + if (FileNameNode != nullptr) { + m_FileName = FileNameNode->GetValue(); + } + + MXmlNode* StripMapNode = Node->GetNode("FileNameStripMap"); + if (StripMapNode != nullptr) { + m_FileNameStripMap = StripMapNode->GetValue(); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleLoaderMeasurementsL0::CreateXmlConfiguration() +{ + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + new MXmlNode(Node, "FileName", m_FileName); + new MXmlNode(Node, "FileNameStripMap", m_FileNameStripMap); + return Node; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleLoaderMeasurementsL0::ShowOptionsGUI() +{ + //! Show the options GUI + + MGUIOptionsLoaderMeasurementsL0* Options = new MGUIOptionsLoaderMeasurementsL0(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +// MModuleLoaderMeasurementsL0.cxx: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleSaverMeasurementsFITS.cxx b/src/MModuleSaverMeasurementsFITS.cxx index d016579b..ff196997 100644 --- a/src/MModuleSaverMeasurementsFITS.cxx +++ b/src/MModuleSaverMeasurementsFITS.cxx @@ -116,6 +116,20 @@ bool MModuleSaverMeasurementsFITS::Initialize() return false; } + // If output data level is 0, call the external m_L1aWriter + if (m_OutputDataLevel == 0) { + // L1b need strip pairing which is not a requirement for L1a. + m_PreceedingModules.clear(); + m_PreceedingModulesHardRequirement.clear(); + AddPreceedingModuleType(MAssembly::c_EventLoader); + + if (m_L1aWriter.Create(m_FileName) == false) { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<addKey("TELESCOP", "COSI", "Mission name"); m_PrimaryHDU->addKey("INSTRUME", "GeD", "Instrument name"); - m_PrimaryHDU->addKey("OBS_ID", "YYMMDD", "Observation ID"); //OBS_ID should have the same YYMMDD as the filename - m_PrimaryHDU->addKey("DATE-OBS", "yyyy-mm-ddThh:mm:ss", "Start Date"); //DATE-OBS should have the start date and time of the data, and this should match the YYMMDD in the filename + m_PrimaryHDU->addKey("OBS_ID", "YYYYMMDD", "Observation ID"); //OBS_ID should have the same YYYYMMDD as the filename + m_PrimaryHDU->addKey("DATE-OBS", "yyyy-mm-ddThh:mm:ss", "Start Date"); //DATE-OBS should have the start date and time of the data, and this should match the YYYYMMDD in the filename m_PrimaryHDU->addKey("DATE-END", "yyyy-mm-ddThh:mm:ss", "Stop Date"); //DATE-END should have the stop time of the data, i.e. the last timestamp m_PrimaryHDU->addKey("ORIGIN", "SSL", "Organization"); @@ -165,79 +179,81 @@ bool MModuleSaverMeasurementsFITS::CreateFITSFile(MString FileName) string seqHitFormat = isL1b ? "PB(50)" : "10B"; string hitFormat = isL1b ? "PE(50)" : "10E"; - std::vector colNames = { - "TIME", "EVENTID", "EVENTCLASS", "NUMHIT", "SEQHIT", - "X", "Y", "Z", - "X_ERR", "Y_ERR", "Z_ERR", - "ENERGY", "ENERGY_ERR", - "RECOILDIR", "RECOILDIR_ERR" - }; - - std::vector colFormats = { - "1D", // TIME - "1J", // EVENTID - "1B", // EVENTCLASS - "1B", // NUMHIT - seqHitFormat, // SEQHIT - hitFormat, // X - hitFormat, // Y - hitFormat, // Z - hitFormat, // X_ERR - hitFormat, // Y_ERR - hitFormat, // Z_ERR - hitFormat, // ENERGY - hitFormat, // ENERGY_ERR - "3E", // RECOILDIR - "3E", // RECOILDIR_ERR - }; - - std::vector colUnits = { - "s", "", "", "", "", - "cm", "cm", "cm", - "cm", "cm", "cm", - "keV", "keV", - "", "" + struct ColSpec { string name, format, unit, comment; bool isL1bOnly; }; + const std::vector colsGedL1b = { + {"TIME", "1D", "s", "Mission Time in sec since 01 Jan 2025 00:00:00", false}, + {"EVENTID", "1J", "", "Event ID, unique number that starts at 1 each day", false}, + {"EVENTTYPE", "1B", "", "Type of event", true }, + {"EVENTCLASS", "1B", "", "0=Compton, 1=photoabsor, 2=charge particle, 4=pair, 5=unknown", false}, + {"NUMHIT", "1B", "", "Number of Hits", false}, + {"STATTEST", "8E", "", "Statistical Test", true }, + {"VETO", "1B", "", "Veto flag: 0=none, 1=hard ACD veto, 2=soft ACD veto, 3=guard ring veto", true }, + {"SEQHIT", seqHitFormat, "", "Sequence of the Hits", false}, + {"X", hitFormat, "cm", "X Location of the HITS", false}, + {"Y", hitFormat, "cm", "Y Location of the HITS", false}, + {"Z", hitFormat, "cm", "Z Location of the HITS", false}, + {"X_ERR", hitFormat, "cm", "Error of X Location of the HITS", false}, + {"Y_ERR", hitFormat, "cm", "Error on Y Location of the HITS", false}, + {"Z_ERR", hitFormat, "cm", "Error on Z Location of the HITS", false}, + {"ENERGY", hitFormat, "keV", "ENERGY of the HITS", false}, + {"ENERGY_ERR", hitFormat, "keV", "Error on ENERGY of the HITS", false}, + {"RECOILDIR", "3E", "", "Recoil electron direction", false}, + {"RECOILDIR_ERR", "3E", "", "Recoil electron direction error", false}, + {"QUALITY_FLAG", "64A", "", "String to specify quality of event calibration and reconstruction", true }, }; - // L2 drops EVENTTYPE, STATTEST, VETO, and QUALITY_FLAG. - // VETO: 0=none, 1=hard ACD veto, 2=soft ACD veto, 3=guard ring veto - if (m_OutputDataLevel == 1) { - colNames.push_back("EVENTTYPE"); - colFormats.push_back("1B"); - colUnits.push_back(""); - - colNames.push_back("STATTEST"); - colFormats.push_back("8E"); - colUnits.push_back(""); - - colNames.push_back("VETO"); - colFormats.push_back("1B"); - colUnits.push_back(""); - - colNames.push_back("QUALITY_FLAG"); - colFormats.push_back("64A"); - colUnits.push_back(""); + std::vector colNames, colFormats, colUnits; + std::vector colComments; + colNames.reserve(colsGedL1b.size()); + colFormats.reserve(colsGedL1b.size()); + colUnits.reserve(colsGedL1b.size()); + colComments.reserve(colsGedL1b.size()); + + for (const auto& col : colsGedL1b) { + // if we in L2, and col is L1b only, skip it + if (col.isL1bOnly && !isL1b) continue; + colNames.push_back(col.name); + colFormats.push_back(col.format); + colUnits.push_back(col.unit); + colComments.push_back(col.comment); } // Create binary table extension string extName = (m_OutputDataLevel == 2) ? "GED_L2" : "GED_L1B"; m_ScienceTable = m_FITSFile->addTable(extName, 0, colNames, colFormats, colUnits); + { + m_ScienceTable->makeThisCurrent(); + int status = 0; + for (size_t i = 0; i < colComments.size(); ++i) { + char keyName[16]; + snprintf(keyName, sizeof(keyName), "TTYPE%zu", i + 1); + fits_modify_comment(m_FITSFile->fitsPointer(), + keyName, + const_cast(colComments[i].c_str()), + &status); + } + if (status != 0 && g_Verbosity >= c_Warning) { + cout << m_XmlTag << ": fits_modify_comment status=" << status << endl; + } + } + // Add keywords to science table m_ScienceTable->addKey("EXTNAME", extName, "name of this HDU"); m_ScienceTable->addKey("TELESCOP", "COSI", "Telescope mission name"); m_ScienceTable->addKey("INSTRUME", "GED", "Instrument name"); m_ScienceTable->addKey("DATAMODE", "SYNC", "Instrument datamode: SYNC or ASYNC"); + m_ScienceTable->addKey("OBSMODE", "SCANNING", "Spacecraft observing mode"); m_ScienceTable->addKey("OBSERVER", "John Tomsick", "Principal Investigator"); - m_ScienceTable->addKey("OBS_ID", "YYMMDD", "Observation ID"); //should match the YYMMDD of the filename - m_ScienceTable->addKey("OBJECT", "ALL SKY", "Object/Target name or ALL SKY"); + m_ScienceTable->addKey("OBS_ID", "YYYYMMDD", "Observation ID"); //should match the YYYYMMDD of the filename + m_ScienceTable->addKey("OBJECT", "ALLSKY", "Object/Target name or ALLSKY"); 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("CLOCKAPP", false, "If clock corrections are applied (T/F)"); + m_ScienceTable->addKey("CLOCKAPP", "F", "If clock corrections are applied (T/F)"); m_ScienceTable->addKey("DATE-OBS", "yyyy-mm-ddThh:mm:ss", "Start Date"); //placeholder, this will be written after we read through all the events m_ScienceTable->addKey("DATE-END", "yyyy-mm-ddThh:mm:ss", "Stop Date"); // m_ScienceTable->addKey("TSTART", 0.0, "Start time"); //placeholder, this will be written after we read through all the events @@ -246,13 +262,12 @@ bool MModuleSaverMeasurementsFITS::CreateFITSFile(MString FileName) m_ScienceTable->addKey("HDUCLAS1", "ARRAY", "hduclass1"); m_ScienceTable->addKey("HDUCLAS2", "TOTAL", "hduclas2"); m_ScienceTable->addKey("CREATOR", "TBD", "Software that create 1st the file"); - m_ScienceTable->addKey("PROCVER", "TBD", "Processing Version"); - m_ScienceTable->addKey("CALDBVER", "TBD", "CALDB version"); - m_ScienceTable->addKey("SEQPNUM", "TBD", "Times the dataset has been processed"); + m_ScienceTable->addKey("PROCVER", "MM.XX.NN.YY", "Processing Version"); + m_ScienceTable->addKey("CALDBVER", "csYYYYMMDD", "CALDB version"); + m_ScienceTable->addKey("SEQPNUM", "nn", "Times the dataset has been processed"); m_ScienceTable->addKey("ORIGIN", "SSL", "Origin of the FITS files"); - m_ScienceTable->addKey("DATE", "TOTAL", "File creation date"); //DATE should have the date of the file creation (same as primary header) - //CHECKSUM - //DATESUM + m_ScienceTable->addKey("DATE", string(dateBuffer), "File creation date (UTC)"); //DATE should have the date of the file creation (same as primary header) + // CHECKSUM and DATASUM are written at file close in Finalize() if (g_Verbosity >= c_Info) cout<SetAnalysisProgress(MAssembly::c_EventSaver); + return ok; + } + // L2 mode: skip bad events (screening) if (m_OutputDataLevel == 2 && Event->IsBad()) { m_TotalEventsSkipped++; @@ -552,6 +574,12 @@ void MModuleSaverMeasurementsFITS::Finalize() { // Finalize the module + if (m_OutputDataLevel == 0) { + m_L1aWriter.Close(); + MModule::Finalize(); + return; + } + // Write any remaining events in the batch if (m_BatchEventCount > 0) { FlushBatch(); @@ -578,12 +606,12 @@ void MModuleSaverMeasurementsFITS::Finalize() // Update science table HDU m_ScienceTable->addKey("DATE-OBS", string(startBuf), "Start Date"); m_ScienceTable->addKey("DATE-END", string(stopBuf), "Stop Date"); - m_ScienceTable->addKey("TSTART", m_FirstEventTime_RTS, "Start time"); - m_ScienceTable->addKey("TSTOP", m_LastEventTime_RTS, "Stop time"); + m_ScienceTable->addKey("TSTART", m_FirstEventTime_RTS, "Start time (RTS)"); + m_ScienceTable->addKey("TSTOP", m_LastEventTime_RTS, "Stop time (RTS)"); - // Also update OBS_ID to match the start date (YYMMDD format) - char obsIdBuf[8]; - strftime(obsIdBuf, sizeof(obsIdBuf), "%y%m%d", gmtime(&startUnix)); + // Also update OBS_ID to match the start date (YYYYMMDD format) + char obsIdBuf[16]; + strftime(obsIdBuf, sizeof(obsIdBuf), "%Y%m%d", gmtime(&startUnix)); m_PrimaryHDU->addKey("OBS_ID", string(obsIdBuf), "Observation ID"); m_ScienceTable->addKey("OBS_ID", string(obsIdBuf), "Observation ID"); @@ -599,6 +627,20 @@ void MModuleSaverMeasurementsFITS::Finalize() } } + // CHECKSUM and DATASUM keywords on every HDU + if (m_FITSFile != nullptr) { + try { + m_FITSFile->pHDU().writeChecksum(); + if (m_ScienceTable != nullptr) { + m_ScienceTable->writeChecksum(); + } + } catch (const CCfits::FitsException& e) { + if (g_Verbosity >= c_Error) { + cout<= c_Info) { @@ -635,6 +677,8 @@ bool MModuleSaverMeasurementsFITS::ReadXmlConfiguration(MXmlNode* Node) MString Level = OutputLevelNode->GetValue(); if (Level == "L2" || Level == "l2") { m_OutputDataLevel = 2; + } else if (Level == "L1a" || Level == "l1a") { + m_OutputDataLevel = 0; } else { m_OutputDataLevel = 1; } @@ -653,7 +697,14 @@ MXmlNode* MModuleSaverMeasurementsFITS::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0, m_XmlTag); new MXmlNode(Node, "FileName", m_FileName); - new MXmlNode(Node, "OutputLevel", (m_OutputDataLevel == 2) ? "L2" : "L1b"); + + const char* lvl = "L1b"; + switch (m_OutputDataLevel) { + case 0: lvl = "L1a"; break; + case 2: lvl = "L2"; break; + default: lvl = "L1b"; break; // 1 = L1b + } + new MXmlNode(Node, "OutputLevel", lvl); return Node; } diff --git a/src/MModuleSaverMeasurementsL0.cxx b/src/MModuleSaverMeasurementsL0.cxx index dbb6fc6f..d263aa62 100644 --- a/src/MModuleSaverMeasurementsL0.cxx +++ b/src/MModuleSaverMeasurementsL0.cxx @@ -80,6 +80,7 @@ MModuleSaverMeasurementsL0::MModuleSaverMeasurementsL0() : MModule() m_SequenceCount = 0; m_TotalEventsWritten = 0; + m_StripMapLoaded = false; } @@ -111,6 +112,16 @@ bool MModuleSaverMeasurementsL0::Initialize() return false; } + if (m_FileNameStripMap == "") { + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": Strip map file name is required but not set." << endl; + return false; + } + if (m_StripMap.Open(m_FileNameStripMap) == false) { + if (g_Verbosity >= c_Error) cout << m_XmlTag << ": Failed to open strip map file: " << m_FileNameStripMap << endl; + return false; + } + m_StripMapLoaded = true; + m_SequenceCount = 0; m_TotalEventsWritten = 0; @@ -394,8 +405,14 @@ bool MModuleSaverMeasurementsL0::AnalyzeEvent(MReadOutAssembly* Event) } else { eventTime = Event->GetTimeRTS(); } - uint32_t eventSeconds = (uint32_t)eventTime.GetAsSeconds(); - uint32_t eventNanoseconds = eventTime.GetNanoSeconds(); + + // Save time in GPS format + MTime gpsTime = Event->ComputeGPSfromRTSTime(eventTime); + uint32_t eventSeconds = (uint32_t)gpsTime.GetAsSeconds(); + uint32_t eventNanoseconds = (uint32_t)gpsTime.GetNanoSeconds(); + + // CCSDS secondary header subseconds are 40 MHz DCB ticks (25 ns each) + uint32_t subsec40MHz = eventNanoseconds / 25; // Convert to TRUNC_TIME format: 10 bits seconds + 22 bits 4MHz subseconds uint32_t seconds10bit = eventSeconds & 0x3FF; @@ -412,13 +429,13 @@ bool MModuleSaverMeasurementsL0::AnalyzeEvent(MReadOutAssembly* Event) MStripHit* hit = Event->GetStripHit(i); // Get strip info - // custom encoding 11-bit strip ID: [detectorID:4][side:1][stripNumber:6] - // GetStripID() dont actually return 0-2079 in simulation... Cannot map through the mapping file - // TODO: figure out exactly what to do in the future, but now this is a workaround int detID = hit->GetDetectorID() & 0xF; // 4 bits (0-15) int side = hit->IsLowVoltageStrip() ? 0 : 1; // 1 bit int stripNum = hit->GetStripID() & 0x3F; // 6 bits (0-63) - int stripID = (detID << 7) | (side << 6) | stripNum; + int stripID = 0; + if (m_StripMap.HasROIDetSideStrip(detID, hit->IsLowVoltageStrip(), stripNum)) { + stripID = (int)m_StripMap.GetReadOutID(detID, hit->IsLowVoltageStrip(), stripNum); + } bool fastTiming = hit->HasFastTiming(); int adc = (int)hit->GetADCUnits(); int tac = (int)hit->GetTAC(); @@ -456,7 +473,7 @@ bool MModuleSaverMeasurementsL0::AnalyzeEvent(MReadOutAssembly* Event) m_SequenceCount = (m_SequenceCount + 1) & 0x3FFF; // Write Secondary Header - WriteSecondaryHeader(eventSeconds, eventNanoseconds); + WriteSecondaryHeader(eventSeconds, subsec40MHz); // Write TRUNC_TIME (4 bytes) WriteUInt32BE(truncTime); @@ -518,6 +535,11 @@ bool MModuleSaverMeasurementsL0::ReadXmlConfiguration(MXmlNode* Node) m_FileName = FileNameNode->GetValue(); } + MXmlNode* StripMapNode = Node->GetNode("FileNameStripMap"); + if (StripMapNode != nullptr) { + m_FileNameStripMap = StripMapNode->GetValue(); + } + return true; } @@ -531,6 +553,7 @@ MXmlNode* MModuleSaverMeasurementsL0::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0, m_XmlTag); new MXmlNode(Node, "FileName", m_FileName); + new MXmlNode(Node, "FileNameStripMap", m_FileNameStripMap); return Node; } diff --git a/src/MReadOutAssembly.cxx b/src/MReadOutAssembly.cxx index d4084d83..adfe48bc 100644 --- a/src/MReadOutAssembly.cxx +++ b/src/MReadOutAssembly.cxx @@ -449,6 +449,32 @@ MTime MReadOutAssembly::ComputeUTCfromRTSTime(MTime RTSTime) //////////////////////////////////////////////////////////////////////////////// +MTime MReadOutAssembly::ComputeRTSfromGPSTime(MTime GPSTime) +{ + //! Compute RTS time from GPS by converting to UTC, then call ComputeRTSfromUTCTime + //! UTC = GPS - 18 + MTime GPS_Unix = MTime(1980,1,6,0,0,0,0); + MTime UTCTime = GPS_Unix + GPSTime - 18; + return ComputeRTSfromUTCTime(UTCTime); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MTime MReadOutAssembly::ComputeGPSfromRTSTime(MTime RTSTime) +{ + //! Compute GPS time from RTS by calling ComputeUTCfromRTSTime then converting UTC to GPS + //! GPS = UTC + 18 + MTime GPS_Unix = MTime(1980,1,6,0,0,0,0); + MTime UTCTime = ComputeUTCfromRTSTime(RTSTime); + return UTCTime - GPS_Unix + 18; +} + + +//////////////////////////////////////////////////////////////////////////////// + + bool MReadOutAssembly::Parse(MString& Line, int Version) { // Handles SE, TI, RO, IA diff --git a/src/MStripMap.cxx b/src/MStripMap.cxx index 8302a8e7..c8a4f765 100644 --- a/src/MStripMap.cxx +++ b/src/MStripMap.cxx @@ -115,6 +115,14 @@ bool MStripMap::Open(MString FileName) // Sort by m_ReadOutID: sort(m_StripMappings.begin(), m_StripMappings.end(), [](const MSingleStripMapping& A, const MSingleStripMapping& B) { return A.m_ReadOutID < B.m_ReadOutID; }); + // build the map + m_DetSideStripToROI.clear(); + m_DetSideStripToROI.reserve(m_StripMappings.size()); + for (const MSingleStripMapping& SM : m_StripMappings) { + unsigned int Key = (SM.m_DetectorID << 8) | ((SM.m_IsLowVoltage ? 0u : 1u) << 7) | SM.m_StripNumber; + m_DetSideStripToROI[Key] = SM.m_ReadOutID; + } + return true; } @@ -126,6 +134,14 @@ bool MStripMap::UpdateASICPolarities(vector>> ASICPolarit for (MSingleStripMapping& S : m_StripMappings) { S.m_IsLowVoltage = ASICPolarities[S.m_DetectorID][S.m_IsPrimary][S.m_ASICID]; } + + // rebuild the map since the key (low voltage) may have changed + m_DetSideStripToROI.clear(); + m_DetSideStripToROI.reserve(m_StripMappings.size()); + for (const MSingleStripMapping& SM : m_StripMappings) { + unsigned int Key = (SM.m_DetectorID << 8) | ((SM.m_IsLowVoltage ? 0u : 1u) << 7) | SM.m_StripNumber; + m_DetSideStripToROI[Key] = SM.m_ReadOutID; + } } return true; } @@ -188,5 +204,29 @@ unsigned int MStripMap::GetStripNumber(unsigned int ROI) const } +//////////////////////////////////////////////////////////////////////////////// + + +//! Check whether (detector, side, strip) maps to a known read-out ID +bool MStripMap::HasROIDetSideStrip(unsigned int DetectorID, bool IsLowVoltage, unsigned int StripNumber) const +{ + unsigned int Key = (DetectorID << 8) | ((IsLowVoltage ? 0u : 1u) << 7) | StripNumber; + return m_DetSideStripToROI.find(Key) != m_DetSideStripToROI.end(); +} + + +//////////////////////////////////////////////////////////////////////////////// + +unsigned int MStripMap::GetReadOutID(unsigned int DetectorID, bool IsLowVoltage, unsigned int StripNumber) const +{ + unsigned int Key = (DetectorID << 8) | ((IsLowVoltage ? 0u : 1u) << 7) | StripNumber; + auto Iter = m_DetSideStripToROI.find(Key); + if (Iter == m_DetSideStripToROI.end()) { + throw MExceptionValueNotFound(Key, "(det/side/strip) tuple in strip map"); + } + return Iter->second; +} + + // MStripMap.cxx: the end... ////////////////////////////////////////////////////////////////////////////////