Modify EditTOFDiffractomerGeometry. Refs #3599.

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

+/*WIKI*
+
+
+
+This algorithm edits the diffractomer detectors' parameters
+
+
+
+
+*WIKI*/
+
 #include "MantidAlgorithms/EditTOFPowderDiffractomerGeometry.h"
+#include "MantidAPI/IAlgorithm.h"
+#include "MantidGeometry/IDetector.h"
+#include "MantidAPI/ISpectrum.h"
 #include "MantidKernel/System.h"
 #include "MantidKernel/ArrayProperty.h"
 
@@ -10,7 +24,7 @@ namespace Mantid
 namespace Algorithms
 {
 
-  //DECLARE_ALGORITHM(EditTOFPowderDiffractomerGeometry);
+  DECLARE_ALGORITHM(EditTOFPowderDiffractomerGeometry);
 
   //---------------------------------------------- ------------------------------------------------
   /** Constructor
@@ -41,15 +55,19 @@ namespace Algorithms
    */
   void EditTOFPowderDiffractomerGeometry::init()
   {
-    declareProperty(new WorkspaceProperty<>("Workspace","",Direction::InOut));
-    declareProperty("Name", "");
+    declareProperty(new WorkspaceProperty<>("Workspace", "", Direction::InOut),
+        "Workspace to edit the detector information");
+    declareProperty("InstrumentName", "GenericPowder", "Name of the instrument. ");
     declareProperty("PrimaryFlightPath", -1.0);
-    declareProperty(new ArrayProperty<int32_t>("DetectorIDs", new MandatoryValidator<std::vector<int32_t> >),
-        "Detector IDs.");
-    declareProperty(new ArrayProperty<double>("SecondaryFlightPaths", new MandatoryValidator<std::vector<double> >),
-        "Seconary flight paths for each detector");
-    declareProperty(new ArrayProperty<double>("TwoThetas", new MandatoryValidator<std::vector<double> >),
-        "Two thetas for all detectors");
+    declareProperty(new ArrayProperty<int32_t>("SpectrumIDs", new MandatoryValidator<std::vector<int32_t> >),
+        "Spectrum IDs (note that it is not detector ID or workspace indices).");
+    declareProperty(new ArrayProperty<double>("L2", new MandatoryValidator<std::vector<double> >),
+        "Seconary flight (L2) paths for each detector");
+    declareProperty(new ArrayProperty<double>("Polar", new MandatoryValidator<std::vector<double> >),
+        "Polar angles (two thetas) for detectors");
+    declareProperty(new ArrayProperty<double>("Azimuthal", new MandatoryValidator<std::vector<double> >),
+        "Azimuthal angles (out-of-plain) for detectors");
+    declareProperty("NewInstrument", false, "Add detectors information from scratch");
   }
 
   //----------------------------------------------------------------------------------------------
@@ -60,25 +78,28 @@ namespace Algorithms
 
     // 1. Get Input
     MatrixWorkspace_sptr workspace = getProperty("Workspace");
-    std::string name = getProperty("Name");
     double l1 = this->getProperty("PrimaryFlightPath");
-    const std::vector<int32_t> detectorids = this->getProperty("DetectorIDs");
-    const std::vector<double> l2s = this->getProperty("SecondaryFlightPaths");
-    const std::vector<double> tths = this->getProperty("TwoThetas");
+    const std::vector<int32_t> specids = this->getProperty("SpectrumIDs");
+    const std::vector<double> l2s = this->getProperty("L2");
+    const std::vector<double> tths = this->getProperty("Polar");
+    const std::vector<double> phis = this->getProperty("Azimuthal");
+    const bool newinstrument = getProperty("NewInstrument");
 
     // 2. Check validity
-    g_log.notice() << "L1 = " << l1 << "  # Detector = " << detectorids.size() << std::endl;
+    bool userL1 = true;
     if (l1 <= 0){
-      throw std::invalid_argument("Primary flight path cannot be less or equal to 0");
+      userL1 = false;
+    } else {
+      g_log.information() << "L1 = " << l1 << "  # Detector = " << std::endl;
     }
 
-    if (detectorids.size() != l2s.size() || l2s.size() != tths.size()){
-      throw std::invalid_argument("Input Detector IDs, Secondary Flight Paths, and Two Thetas have different items number");
+    if (specids.size() != l2s.size() || l2s.size() != tths.size() || phis.size() != l2s.size()){
+      throw std::invalid_argument("Input Spectrum, L2 (Secondary Flight Paths), Polar angles (Two Thetas) and Azimuthal have different items number");
     }
 
     for (size_t ib = 0; ib < l2s.size(); ib ++){
-      g_log.information() << "Detector " << detectorids[ib] << "  L2 = " << l2s[ib] << "  2Theta = " << tths[ib] << std::endl;
-      if (detectorids[ib] < 0){
+      g_log.information() << "Detector " << specids[ib] << "  L2 = " << l2s[ib] << "  2Theta = " << tths[ib] << std::endl;
+      if (specids[ib] < 0){
         throw std::invalid_argument("Detector ID cannot be less than 0");
       }
       if (l2s[ib] <= 0.0){
@@ -86,22 +107,43 @@ namespace Algorithms
       }
     }
 
-    // 3. Generate a new instrument or use the old one
-    // a) Create a new instrument
-    if (name == ""){
-      name = "Generic";
+    // 3. Generate a new instrument
+    std::string name = "";
+    if (newinstrument){
+      std::string tempname = getProperty("InstrumentName");
+      name += tempname;
+    } else {
+      name += workspace->getInstrument()->getName();
     }
-    Geometry::Instrument_sptr instrument(new Geometry::Instrument(name));
-    workspace->setInstrument(instrument);
 
-    // getBaseInstrument();
+    //  Create a new instrument from scratch
+    Geometry::Instrument_sptr instrument(new Geometry::Instrument(name));
     if (instrument.get() == 0)
     {
-        g_log.error("Trying to use a Parametrized Instrument as an Instrument.");
-        throw std::runtime_error("Trying to use a Parametrized Instrument as an Instrument.");
+      g_log.error("Trying to use a Parametrized Instrument as an Instrument.");
+      throw std::runtime_error("Trying to use a Parametrized Instrument as an Instrument.");
+    }
+
+    // 4. Get information from original: L1
+    double mL1;
+    if (newinstrument){
+      // use input L1
+      mL1 = l1;
+      if (l1 <= 0){
+        g_log.error() << "Input L1 = " << l1 << "  < 0 Is Not Allowed!" << std::endl;
+        throw std::invalid_argument("Input L1 is not allowed.");
+      }
+    } else {
+      if (l1 > 0){
+        mL1 = l1;
+      } else {
+        Kernel::V3D sourcepos = workspace->getInstrument()->getSource()->getPos();
+        mL1 = sourcepos.Z()*-1;
+        g_log.notice() << "source position = " << sourcepos.X() << ", " << sourcepos.Y() << ", " << sourcepos.Z() << std::endl;
+      }
     }
 
-    // 3. Source and sample information
+    // 5. Source and sample information
     Geometry::ObjComponent *samplepos = new Geometry::ObjComponent("Sample", instrument.get());
     instrument->add(samplepos);
     instrument->markAsSamplePos(samplepos);
@@ -110,23 +152,58 @@ namespace Algorithms
     Geometry::ObjComponent *source = new Geometry::ObjComponent("Source", instrument.get());
     instrument->add(source);
     instrument->markAsSource(source);
-    source->setPos(0.0, 0.0, -1.0 * l1);
+    source->setPos(0.0, 0.0, -1.0 * mL1);
 
-    // 4. Add detector information
-    // Create a new detector. Instrument will take ownership of pointer so no need to delete.
-    for (size_t ib = 0; ib < detectorids.size(); ib ++){
-      Geometry::Detector *detector = new Geometry::Detector("det", detectorids[ib], samplepos);
-      Kernel::V3D pos;
+    // 6. Keep original instrument and set the new instrument
+    Geometry::Instrument_const_sptr oldinstrument = workspace->getInstrument();
+    workspace->setInstrument(instrument);
+
+    // 6. Add detector information
+    spec2index_map * spec2indexmap = workspace->getSpectrumToWorkspaceIndexMap();
 
-      double l2 = l2s[ib];
-      double tth = tths[ib];
-      pos.spherical(l2, tth, 0.0 );
+    for (size_t i = 0; i < specids.size(); i ++){
+      // Create a new detector. Instrument will take ownership of pointer so no need to delete.
+      Geometry::Detector *detector = new Geometry::Detector("det", specids[i]+100, samplepos);
+      spec2index_map::iterator it = spec2indexmap->find(specids[i]);
+      if (it == spec2indexmap->end()){
+        g_log.error() << "Spectrum ID " << specids[i] << " is not found" << std::endl;
+        continue;
+      }
+      size_t workspaceindex = it->second;
+      g_log.notice() << "workspace index = " << workspaceindex << " is for Spectrum " << specids[i] << std::endl;
+
+      Kernel::V3D pos;
+      double l2 = l2s[i];
+      double tth = tths[i];
+      double phi = phis[i];
+      pos.spherical(l2, tth, phi);
       detector->setPos(pos);
+      double tl2, t2t, tph;
+      detector->getPos().getSpherical(tl2, t2t, tph);
+      g_log.notice() << "Detector Position (Input) = " << l2 << ", " << tth << ", " << phi << std::endl;
+      g_log.notice() << "New Detector " << specids[i] << ":  L2, 2theta, phi = " << tl2 << ", " << t2t << ", " << tph << std::endl;
 
       // add copy to instrument and mark it
+      API::ISpectrum *spectrum = workspace->getSpectrum(workspaceindex);
+      if (!spectrum){
+        g_log.error() << g_log.error() << "Spectrum ID " << specids[i] << " does not exist!  Skip setting detector parameters to this spectrum" << std::endl;
+        continue;
+      } else {
+        g_log.notice() << "Raw: Spectrum " << specids[i] << ": # Detectors =  " << spectrum->getDetectorIDs().size() << std::endl;
+      }
+
+      spectrum->clearDetectorIDs();
+      spectrum->addDetectorID(specids[i]+100);
       instrument->add(detector);
       instrument->markAsDetector(detector);
-    }
+
+      // Check
+      g_log.notice() << "New: Spectrum: " << specids[i] << ": # Detectros = " << spectrum->getDetectorIDs().size() << std::endl;
+
+
+    } // for i
+
+    delete spec2indexmap;
 
     return;
   }

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

@@ -2,12 +2,17 @@
 #define MANTID_ALGORITHMS_EDITTOFPOWDERDIFFRACTOMERGEOMETRYTEST_H_
 
 #include <cxxtest/TestSuite.h>
+#include "MantidAlgorithms/EditTOFPowderDiffractomerGeometry.h"
+#include "MantidDataHandling/LoadNexusProcessed.h"
+#include "MantidGeometry/Instrument.h"
+#include "MantidAPI/ISpectrum.h"
+#include "MantidGeometry/IDetector.h"
+
 #include "MantidKernel/Timer.h"
 #include "MantidKernel/System.h"
 #include <iostream>
 #include <iomanip>
 
-#include "MantidAlgorithms/EditTOFPowderDiffractomerGeometry.h"
 
 using namespace Mantid;
 using namespace Mantid::Algorithms;
@@ -17,8 +22,62 @@ class EditTOFPowderDiffractomerGeometryTest : public CxxTest::TestSuite
 {
 public:
 
-  void test_Something()
+  void test_Initialize(){
+
+    EditTOFPowderDiffractomerGeometry editdetector;
+    TS_ASSERT_THROWS_NOTHING(editdetector.initialize());
+    TS_ASSERT(editdetector.isInitialized());
+
+  }
+
+  void test_SingleSpectrum()
   {
+    // 1. Init
+    EditTOFPowderDiffractomerGeometry editdetector;
+    editdetector.initialize();
+
+    // 2. Load Data
+    Mantid::DataHandling::LoadNexusProcessed loader;
+    loader.initialize();
+    loader.setProperty("Filename","PG3_2583.nxs");
+    const std::string inputWS = "inputWS";
+    loader.setPropertyValue("OutputWorkspace",inputWS);
+    loader.execute();
+
+    // 3. Set Property
+    TS_ASSERT_THROWS_NOTHING( editdetector.setPropertyValue("Workspace", inputWS) );
+    TS_ASSERT_THROWS_NOTHING( editdetector.setPropertyValue("SpectrumIDs","1") );
+    TS_ASSERT_THROWS_NOTHING( editdetector.setPropertyValue("L2","3.45") );
+    TS_ASSERT_THROWS_NOTHING( editdetector.setPropertyValue("Polar","90.09") );
+    TS_ASSERT_THROWS_NOTHING( editdetector.setPropertyValue("Azimuthal","1.84") );
+    TS_ASSERT_THROWS_NOTHING( editdetector.setProperty("NewInstrument", false) );
+
+    // 4. Run
+    TS_ASSERT_THROWS_NOTHING( editdetector.execute() );
+    TS_ASSERT( editdetector.isExecuted() );
+
+    // 5. Check result
+    Mantid::API::MatrixWorkspace_sptr workspace;
+    TS_ASSERT_THROWS_NOTHING( workspace = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+      (Mantid::API::AnalysisDataService::Instance().retrieve(inputWS)) );
+
+    API::ISpectrum* spectrum1 = workspace->getSpectrum(0);
+    Geometry::Instrument_const_sptr instrument = workspace->getInstrument();
+
+    std::set<detid_t> detids = spectrum1->getDetectorIDs();
+    TS_ASSERT_EQUALS(detids.size(), 1);
+    detid_t detid = 0;
+    std::set<detid_t>::iterator it;
+    for (it = detids.begin(); it != detids.end(); ++it){
+      detid = *it;
+    }
+    Geometry::IDetector_const_sptr detector = instrument->getDetector(detid);
+    double r, tth, phi;
+    detector->getPos().getSpherical(r, tth, phi);
+    TS_ASSERT_DELTA(r, 3.45, 0.000001);
+    TS_ASSERT_DELTA(tth, 90.09, 0.000001);
+    TS_ASSERT_DELTA(phi, 1.84, 0.000001);
+
   }