Adding PCF spectrum writer

radixio/spectrum.cc

@@ -2,9 +2,25 @@
 
 namespace radix
 {
-short Spectrum::nrps() const { return mNrps; }
+short Spectrum::nrps() const
+{
+  // Calculate the number of records per spectrum required
+  // Used by PCF file to allocate space in binary file structure
+
+  // Go through each SpectrumData object and find the maximum number of channels
+  int maxChannels = 0;
+  for (size_t i = 0; i < mSpectrumData.size(); ++i)
+  {
+    if (mSpectrumData[i].numberOfChannels > maxChannels)
+    {
+      maxChannels = mSpectrumData[i].numberOfChannels;
+    }
+  }
 
-void Spectrum::setNrps(short nrps) { mNrps = nrps; }
+  short nrps = short((maxChannels / 64) + 1);
+
+  return nrps;
+}
 
 bool Spectrum::dhsVersion() const { return mDhsVersion; }
 

radixio/spectrum.hh

@@ -26,7 +26,6 @@ class RADIX_PUBLIC Spectrum
   };  // struct SpectrumData
 
   short nrps() const;
-  void setNrps(short nrps);
   bool dhsVersion() const;
   void setDhsVersion(bool dhsVersion);
   std::string energyCalibrationLabel() const;
@@ -109,8 +108,6 @@ class RADIX_PUBLIC Spectrum
   // PCF file objects -
   // TODO work out what/how to change to support SPE and PCF
 
-  // Number of records per spectrum
-  short mNrps;
   // Is the file a DHS version?
   bool mDhsVersion = false;
   // Non-DHS header objects

radixio/spectrumpcfstream.i.hh

@@ -38,7 +38,6 @@ bool SpectrumPCFStream<data_type>::read_from(const std::string &file)
   radix_line("  Reading header...");
   // Number of records per spectrum (NRPS)
   short nrps = stream.readShort();
-  mData->setNrps(nrps);
   radix_line("    nrps: " << nrps);
   // Version (3-character string)
   std::string version = stream.readString(3);
@@ -111,6 +110,7 @@ bool SpectrumPCFStream<data_type>::read_from(const std::string &file)
   std::string deviationPairPresence = radix::trim_string(stream.readString(30));
   if (deviationPairPresence.compare("DeviationPairsInFile") == 0)
   {
+    mData->setDeviationPairPresence(deviationPairPresence);
     radix_line("    Deviation pairs present; reading...");
     // Set SRSI to correct value
     mData->setSrsi(83);
@@ -144,6 +144,7 @@ bool SpectrumPCFStream<data_type>::read_from(const std::string &file)
   }
   else if (deviationPairPresence.compare("DeviationPairsInFileCompressed") == 0)
   {
+    mData->setDeviationPairPresence(deviationPairPresence);
     radix_line("    Compressed deviation pairs present; reading...");
     // Set SRSI to correct value
     mData->setSrsi(83);
@@ -161,7 +162,7 @@ bool SpectrumPCFStream<data_type>::read_from(const std::string &file)
           for (size_t pair = 0; pair < 20; ++pair)
           {
             float energy = float(stream.readShort()) / 10.f;
-            short offset = float(stream.readShort()) / 10.f;
+            float offset = float(stream.readShort()) / 10.f;
             deviationPairEnergies[column][panel][mca][pair] = energy;
             deviationPairOffsets[column][panel][mca][pair]  = offset;
           }
@@ -174,6 +175,7 @@ bool SpectrumPCFStream<data_type>::read_from(const std::string &file)
   }
   else
   {
+    mData->setDeviationPairPresence("");
     radix_line("    No deviation pairs present");
     // Set SRSI to correct value
     mData->setSrsi(2);
@@ -311,6 +313,215 @@ bool SpectrumPCFStream<data_type>::write_to(const std::string &file) const
 {
   bool result = false;
 
+  radix_line("Writing data to PCF file: " << file);
+
+  // Open a stream of the file
+  eafstream stream(file.c_str(), std::ios::out | std::ios::binary);
+  if (stream.is_open() == false)
+  {
+    return result;
+  }
+
+  // File should be little endian
+  stream.setReverseBytes(false);
+
+  // Write the header
+  radix_line("  Writing header...");
+  // Number of records per spectrum (NRPS)
+  stream.writeShort(mData->nrps());
+  // Rest of header depends on version type
+  if (mData->dhsVersion())
+  {
+    radix_line("    File has DHS-type header");
+    stream.write("DHS", 3);
+    stream.write(mData->lastModifiedHash().c_str(), 7);
+    stream.write(mData->uuid().c_str(), 36);
+    stream.write(mData->inspection().c_str(), 16);
+    stream.writeShort(mData->laneNumber());
+    stream.write(mData->measurementRemark().c_str(), 26);
+    stream.write(mData->instrumentType().c_str(), 28);
+    stream.write(mData->manufacturer().c_str(), 28);
+    stream.write(mData->instrumentModel().c_str(), 18);
+    stream.write(mData->instrumentID().c_str(), 18);
+    stream.write(mData->itemDescription().c_str(), 20);
+    stream.write(mData->measurementLocationName().c_str(), 16);
+    stream.write(mData->measurementLocationCoords().c_str(), 16);
+    stream.writeShort(mData->itemDetectorDistance());
+    stream.writeShort(mData->occupancyNumber());
+    stream.write(mData->cargoType().c_str(), 16);
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(stream.bytesWritten() + 248);
+  }
+  else
+  {
+    radix_line("    File has non-DHS-type header");
+    stream.write("   ", 3);
+    stream.write(mData->energyCalibrationLabel().c_str(), 4);
+    stream.writeFloat(mData->energyCalibrationOffset());
+    stream.writeFloat(mData->energyCalibrationGain());
+    stream.writeFloat(mData->energyCalibrationQuadraticTerm());
+    stream.writeFloat(mData->energyCalibrationCubicTerm());
+    stream.writeFloat(mData->energyCalibrationLowEnergy());
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(stream.bytesWritten() + 7);
+    // Skip to 256 bytes (end of header)
+    stream.write(" ", 256 - stream.bytesWritten());
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(256);
+  }
+  radix_line("  End of header: " << stream.bytesWritten()
+                                 << " bytes written so far");
+
+  // Write deviation pairs
+  radix_line("  Writing deviation pairs...");
+  if (mData->deviationPairPresence().compare("DeviationPairsInFile") == 0)
+  {
+    radix_line("    Regular deviation pairs found:");
+    // Write regular deviation pairs
+    stream.write(mData->deviationPairPresence().c_str(), 20);
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(stream.bytesWritten() + 20);
+    // Skip to 512 bytes (start of pairs)
+    stream.write(" ", 512 - stream.bytesWritten());
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(512);
+
+    // Write the pairs themselves
+    for (size_t column = 0; column < 2; ++column)
+    {
+      radix("      ");
+      for (size_t panel = 0; panel < 8; ++panel)
+      {
+        for (size_t mca = 0; mca < 8; ++mca)
+        {
+          for (size_t pair = 0; pair < 20; ++pair)
+          {
+            stream.writeFloat(
+                mData->deviationPairEnergies()[column][panel][mca][pair]);
+            stream.writeFloat(
+                mData->deviationPairOffsets()[column][panel][mca][pair]);
+          }
+          radix(".");
+        }
+        radix("-");
+      }
+      radix_line("|");
+    }
+  }
+  else if (mData->deviationPairPresence().compare(
+               "DeviationPairsInFileCompressed") == 0)
+  {
+    radix_line("    Compressed deviation pairs found:");
+    // Write compressed deviation pairs
+    stream.write(mData->deviationPairPresence().c_str(), 30);
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(stream.bytesWritten() + 30);
+    // Skip to 512 bytes (start of pairs)
+    stream.write(" ", 512 - stream.bytesWritten());
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(512);
+
+    // Write the pairs themselves
+    for (size_t column = 0; column < 4; ++column)
+    {
+      radix("      ");
+      for (size_t panel = 0; panel < 8; ++panel)
+      {
+        for (size_t mca = 0; mca < 8; ++mca)
+        {
+          for (size_t pair = 0; pair < 20; ++pair)
+          {
+            // Calculate the pair values
+            short energy =
+                short(mData->deviationPairEnergies()[column][panel][mca][pair] *
+                      10.f);
+            short offset = short(
+                mData->deviationPairOffsets()[column][panel][mca][pair] * 10.f);
+            stream.writeShort(energy);
+            stream.writeShort(offset);
+          }
+          radix(".");
+        }
+        radix("-");
+      }
+      radix_line("|");
+    }
+  }
+  else
+  {
+    // Don't write deviation pairs
+    radix_line("    No deviation pairs found:");
+  }
+
+  // Write spectrum data
+  // Spectrum data values
+  std::string title, source, description, dateTime, tag;
+  float liveTime, totalTime, energyCalOffset, energyCalGain,
+      energyCalQuadraticTerm, energyCalCubicTerm, energyCalLowEnergyTerm,
+      occupancyFlag, totalNeutronCount;
+  int numberOfChannels;
+  std::vector<float> countsByChannel;
+  for (size_t i = 0; i < mData->spectrumDataCount(); ++i)
+  {
+    radix_line("  Writing spectrum data for spectrum "
+               << i + 1 << " of " << mData->spectrumDataCount());
+
+    // Get the spectrum
+    mData->spectrumData(i, title, source, description, dateTime, tag, liveTime,
+                        totalTime, energyCalOffset, energyCalGain,
+                        energyCalQuadraticTerm, energyCalCubicTerm,
+                        energyCalLowEnergyTerm, occupancyFlag,
+                        totalNeutronCount, numberOfChannels, countsByChannel);
+
+    // Write spectrum header
+    radix_line("    Writing spectrum header:");
+    std::stringstream ss;
+    char delim = 255;
+    ss << delim << trim_string(title) << delim << trim_string(description)
+       << delim << trim_string(source);
+    stream.write(ss.str().c_str(), 180);
+    stream.write(dateTime.c_str(), 23);
+    stream.write(tag.c_str(), 1);
+    // Add bytes written directly to eafstream count
+    stream.setBytesWritten(stream.bytesWritten() + 204);
+    stream.writeFloat(liveTime);
+    stream.writeFloat(totalTime);
+    // Three unused floats in data structure
+    stream.writeFloat(0.f);
+    stream.writeFloat(0.f);
+    stream.writeFloat(0.f);
+    stream.writeFloat(energyCalOffset);
+    stream.writeFloat(energyCalGain);
+    stream.writeFloat(energyCalQuadraticTerm);
+    stream.writeFloat(energyCalCubicTerm);
+    stream.writeFloat(energyCalLowEnergyTerm);
+    stream.writeFloat(occupancyFlag);
+    stream.writeFloat(totalNeutronCount);
+    stream.writeInt(numberOfChannels);
+
+    // Write spectrum data
+    radix_line("    Writing spectrum data:");
+    for (size_t channel = 0; channel < countsByChannel.size(); ++channel)
+    {
+      stream.writeFloat(countsByChannel[channel]);
+    }
+    int maxChannels = 64 * (mData->nrps() - 1);
+    radix_line("      Written " << countsByChannel.size()
+                                << " data points; max size is " << maxChannels);
+    if (countsByChannel.size() < maxChannels)
+    {
+      radix_line("      Writing " << maxChannels - countsByChannel.size()
+                                  << " zeroes to pad to max");
+      for (int i = countsByChannel.size(); i < maxChannels; ++i)
+      {
+        stream.writeFloat(0.f);
+      }
+    }
+
+    radix_line("  Spectrum write complete");
+  }
+
+  result = true;
   return result;
 }
 

radixio/tests/tstSpectrum.cc

@@ -20,16 +20,19 @@ TEST(RadixIO, SpectrumFromPCF)
   EXPECT_EQ(30, testSpectrum->spectrumDataCount());
 }
 
-// TEST(RadixIO, SpectrumToPCF)
-//{
-//  std::string testPCFFile = "test.pcf";
-
-//  Spectrum::SP testSpectrum = std::make_shared<Spectrum>();
+TEST(RadixIO, SpectrumToPCF)
+{
+  // Read in a test spectrum
+  std::string testReadPCFFile =
+      radix::to_native_path(std::string(dirname(__FILE__) + "/data/235F.pcf"));
+  Spectrum::SP testSpectrum1 = std::make_shared<Spectrum>();
+  SpectrumPCFStream<Spectrum::SP> testStream1(testSpectrum1);
+  ASSERT_TRUE(testStream1.read_from(testReadPCFFile));
 
-//  // Load the spectrum
-//  SpectrumPCFStream<Spectrum::SP> stream(testSpectrum);
-//  ASSERT_TRUE(stream.write_to(testPCFFile));
-//}
+  // Write the spectrum back out to a file
+  std::string testWritePCFFile = "test.pcf";
+  ASSERT_TRUE(testStream1.write_to(testWritePCFFile));
+}
 
 TEST(RadixIO, SpectrumFromSPE)
 {

[Purves, Murray]

@jap I think it is the first time addSpectrumData is called that is the issue: the PCF read/write also spits out an empty first SpectrumData object