Checkpointing work. Refs #6737.

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CreateLogTimeCorrection.h

@@ -4,7 +4,7 @@
 #include "MantidKernel/System.h"
 #include "MantidAPI/Algorithm.h"
 #include "MantidAPI/MatrixWorkspace.h"
-#include "MantidDataObjects/Workspace2D.h"
+#include "MantidDataObjects/TableWorkspace.h"
 #include "MantidGeometry/Instrument.h"
 
 namespace Mantid
@@ -56,7 +56,23 @@ namespace Algorithms
     /// Implement abstract Algorithm methods
     void exec();
 
+    /// Get instrument geometry setup including L2 for each detector and L1
+    void getInstrumentSetup();
+
+    /// Calculate the log time correction for each pixel, i.e., correcton from event time at detector to time at sample
+    void calculateCorrection();
+
+    /// Write L2 map and correction map to a TableWorkspace
+    DataObjects::TableWorkspace_sptr generateCorrectionTable();
+
+    /// Write correction map to a text file
+    void writeCorrectionToFile(std::string filename);
+
     API::MatrixWorkspace_sptr m_dataWS;
+    Geometry::Instrument_const_sptr m_instrument;
+    std::map<int, double> m_l2map;
+    std::map<int, double> m_correctionMap;
+    double m_L1;
   };
 
 

Code/Mantid/Framework/Algorithms/src/CreateLogTimeCorrection.cpp

@@ -1,7 +1,10 @@
 #include "MantidAlgorithms/CreateLogTimeCorrection.h"
 #include "MantidAPI/WorkspaceProperty.h"
 #include "MantidAPI/FileProperty.h"
+#include "MantidAPI/TableRow.h"
 
+#include <fstream>
+#include <iomanip>
 
 using namespace Mantid;
 using namespace Mantid::API;
@@ -53,7 +56,7 @@ namespace Algorithms
     auto inpwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input);
     declareProperty(inpwsprop, "Name of the input workspace to generate log correct from.");
 
-    auto outwsprop = new WorkspaceProperty<Workspace2D>("Outputworkspace", "AnonymousOut", Direction::Output);
+    auto outwsprop = new WorkspaceProperty<TableWorkspace>("Outputworkspace", "AnonymousOut", Direction::Output);
     declareProperty(outwsprop, "Name of the output workspace containing the corrections.");
 
     auto fileprop = new FileProperty("OutputFilename", "", FileProperty::Save);
@@ -68,45 +71,154 @@ namespace Algorithms
   {
     // 1. Process input
     m_dataWS = getProperty("InputWorkspace");
+    m_instrument = m_dataWS->getInstrument();
+    if (!m_instrument)
+    {
+      stringstream errss;
+      errss << "Input matrix workspace " << m_dataWS->name() << " does not have instrument. ";
+      g_log.error(errss.str());
+      throw runtime_error(errss.str());
+    }
 
     // 2. Explore geometry
-    Instrument_const_sptr m_instrument = m_dataWS->getInstrument();
+    getInstrumentSetup();
 
-    std::vector<detid_t> detids = m_instrument->getDetectorIDs(true);
+    // 3. Calculate log time correction
+    calculateCorrection();
+
+    // 4. Output
+    TableWorkspace_sptr outWS = generateCorrectionTable();
+    setProperty("OutputWorkspace", outWS);
+
+    string filename = getProperty("OutputFilename");
+    if (filename.size() > 0)
+    {
+      writeCorrectionToFile(filename);
+    }
 
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Get instrument geometry setup including L2 for each detector and L1
+    */
+  void CreateLogTimeCorrection::getInstrumentSetup()
+  {
+    // 1. Get sample position and source position
     IObjComponent_const_sptr sample = m_instrument->getSample();
+    if (!sample)
+    {
+      throw runtime_error("No sample has been set.");
+    }
     V3D samplepos = sample->getPos();
 
     IObjComponent_const_sptr source = m_instrument->getSource();
+    if (!source)
+    {
+      throw runtime_error("No source has been set.");
+    }
     V3D sourcepos = source->getPos();
-    double l1 = sourcepos.distance(samplepos);
-
-    g_log.notice() << "Sample position = " << samplepos << ".\n";
-    g_log.notice() << "Source position = " << sourcepos << ", L1 = " << l1 << ".\n";
-    g_log.notice() << "Number of detector/pixels = " << detids.size() << ".\n";
-
-    // 3. Output
-    Workspace2D_sptr m_outWS = boost::dynamic_pointer_cast<Workspace2D>
-        (WorkspaceFactory::Instance().create("Workspace2D", 1, detids.size(), detids.size()));
+    m_L1 = sourcepos.distance(samplepos);
 
-    MantidVec& vecX = m_outWS->dataX(0);
-    MantidVec& vecY = m_outWS->dataY(0);
+    // 2. Get detector IDs
+    std::vector<detid_t> detids = m_instrument->getDetectorIDs(true);
     for (size_t i = 0; i < detids.size(); ++i)
     {
-      vecX[i] = static_cast<double>(detids[i]);
-
       IDetector_const_sptr detector = m_instrument->getDetector(detids[i]);
       V3D detpos = detector->getPos();
-      vecY[i] = detpos.distance(samplepos);
+      double l2 = detpos.distance(samplepos);
+      m_l2map.insert(make_pair(detids[i], l2));
     }
 
-    setProperty("OutputWorkspace", m_outWS);
+    // 3. Output information
+    g_log.information() << "Sample position = " << samplepos << "; "
+                        << "Source position = " << sourcepos << ", L1 = " << m_L1 << "; "
+                        << "Number of detector/pixels = " << detids.size() << ".\n";
+  }
 
+  //----------------------------------------------------------------------------------------------
+  /** Calculate the log time correction for each pixel, i.e., correcton from event time at detector
+    * to time at sample
+    */
+  void CreateLogTimeCorrection::calculateCorrection()
+  {
+    map<int, double>::iterator miter;
+    for (miter = m_l2map.begin(); miter != m_l2map.end(); ++miter)
+    {
+      int detid = miter->first;
+      double l2 = miter->second;
+      double corrfactor = m_L1/(m_L1+l2);
+      m_correctionMap.insert(make_pair(detid, corrfactor));
+    }
   }
 
+  //----------------------------------------------------------------------------------------------
+  /** Write L2 map and correction map to a TableWorkspace
+    */
+  TableWorkspace_sptr CreateLogTimeCorrection::generateCorrectionTable()
+  {
+    TableWorkspace_sptr tablews(new TableWorkspace());
 
+    tablews->addColumn("int", "DetectorID");
+    tablews->addColumn("double", "Correction");
+    tablews->addColumn("double", "L2");
 
+    if (m_l2map.size() != m_correctionMap.size())
+      throw runtime_error("Program logic error!");
+
+    map<int, double>::iterator l2iter, citer;
+    for (l2iter = m_l2map.begin(); l2iter != m_l2map.end(); ++l2iter)
+    {
+      int detid = l2iter->first;
+      double l2 = l2iter->second;
+
+      citer = m_correctionMap.find(detid);
+      if (citer == m_correctionMap.end())
+      {
+        throw runtime_error("Program logic error (B)!");
+      }
+      double correction = citer->second;
+
+      TableRow newrow = tablews->appendRow();
+      newrow << detid << correction << l2;
+    }
+
+    return tablews;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Write correction map to a text file
+    */
+  void CreateLogTimeCorrection::writeCorrectionToFile(string filename)
+  {
+    ofstream ofile;
+    ofile.open(filename.c_str());
+
+    if (ofile.is_open())
+    {
+      map<int, double>::iterator miter;
+      for (miter = m_correctionMap.begin(); miter != m_correctionMap.end(); ++miter)
+      {
+        int detid = miter->first;
+        double corr = miter->second;
+        ofile << detid << "\t" << setw(20) << setprecision(5) << corr << "\n";
+      }
+      ofile.close();
+    }
+    else
+    {
+      g_log.error() << "Unable to open file " << filename << " to write!\n";
+    }
 
+    return;
+  }
 
 } // namespace Algorithms
 } // namespace Mantid
+
+
+
+
+
+
+
+

Code/Mantid/Framework/Algorithms/test/CreateLogTimeCorrectionTest.h

@@ -4,9 +4,19 @@
 #include <cxxtest/TestSuite.h>
 
 #include "MantidAlgorithms/CreateLogTimeCorrection.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidDataHandling/LoadInstrument.h"
+#include "MantidAPI/TableRow.h"
 
 using Mantid::Algorithms::CreateLogTimeCorrection;
 
+using namespace Mantid;
+using namespace Mantid::API;
+using namespace Mantid::DataObjects;
+using namespace Mantid::Kernel;
+
+using namespace std;
+
 class CreateLogTimeCorrectionTest : public CxxTest::TestSuite
 {
 public:
@@ -15,14 +25,96 @@ class CreateLogTimeCorrectionTest : public CxxTest::TestSuite
   static CreateLogTimeCorrectionTest *createSuite() { return new CreateLogTimeCorrectionTest(); }
   static void destroySuite( CreateLogTimeCorrectionTest *suite ) { delete suite; }
 
+  /** Test against a Vulcan run
+    */
+  void test_Vulcan()
+  {
+    MatrixWorkspace_sptr inpws = createEmptyWorkspace("VULCAN");
+    AnalysisDataService::Instance().addOrReplace("Vulcan_Fake", inpws);
+
+    CreateLogTimeCorrection alg;
+    alg.initialize();
+    TS_ASSERT(alg.isInitialized());
+
+    alg.setProperty("InputWorkspace", inpws);
+    alg.setProperty("OutputWorkspace", "CorrectionTable");
+    alg.setProperty("OutputFilename", "VucanCorrection.dat");
+
+    TS_ASSERT_THROWS_NOTHING(alg.execute());
+    TS_ASSERT(alg.isExecuted());
+
+    TableWorkspace_sptr outws = boost::dynamic_pointer_cast<TableWorkspace>(
+          AnalysisDataService::Instance().retrieve("CorrectionTable"));
+    TS_ASSERT(outws);
+
+    int numrows = static_cast<int>(outws->rowCount());
+    TS_ASSERT_EQUALS(numrows, 7392);
+
+    // get some value to check
+    double l1 = 43.754;
+
+    vector<size_t> checkrows;
+    checkrows.push_back(0);
+    checkrows.push_back(100);
+    checkrows.push_back(1000);
+    checkrows.push_back(5000);
+
+    for (size_t i = 0; i < checkrows.size(); ++i)
+    {
+      TableRow row = outws->getRow(i);
+      int detid;
+      double correction, l2;
+      row >> detid >> correction >> l2;
+      TS_ASSERT(detid > 0);
+
+      TS_ASSERT_DELTA(correction*(l1+l2)/l1, 1.0, 0.0001);
+    }
 
-  void test_Something()
+    return;
+  }
+
+  /** Test if there is no instrument in given workspace
+    */
+  void test_NoInstrument()
   {
-    TSM_ASSERT( "You forgot to write a test!", 0);
+    MatrixWorkspace_sptr inpws = createEmptyWorkspace("");
+    AnalysisDataService::Instance().addOrReplace("Vulcan_Fake", inpws);
+
+    CreateLogTimeCorrection alg;
+    alg.initialize();
+    TS_ASSERT(alg.isInitialized());
+
+    alg.setProperty("InputWorkspace", inpws);
+    alg.setProperty("OutputWorkspace", "CorrectionTable");
+    alg.setProperty("OutputFilename", "VucanCorrection.dat");
+
+    TS_ASSERT_THROWS_NOTHING(alg.execute());
+    TS_ASSERT(!alg.isExecuted());
+
+    return;
+  }
+
+  /** Generate an empty Vulcan workspace
+    */
+  API::MatrixWorkspace_sptr createEmptyWorkspace(string instrument)
+  {
+    MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(
+        WorkspaceFactory::Instance().create("Workspace2D", 1, 1, 1));
+
+    if (instrument.size() > 0)
+    {
+      DataHandling::LoadInstrument load;
+      load.initialize();
+      load.setProperty("Workspace", ws);
+      load.setProperty("InstrumentName", instrument);
+      load.execute();
+    }
+
+    return ws;
   }
 
 
 };
 
 
-#endif /* MANTID_ALGORITHMS_CREATELOGTIMECORRECTIONTEST_H_ */
+#endif /* MANTID_ALGORITHMS_CREATELOGTIMECORRECTIONTEST_H_ */