From 388217efc8ebf0ad9eecf6e4a0aee0c9cf498526 Mon Sep 17 00:00:00 2001 From: Wenduo Zhou Date: Mon, 1 Jun 2015 17:28:12 -0400 Subject: [PATCH] Refs #11847. Checkpointing progress on solving issues during test. --- .../ConvertCWSDExpToMomentum.h | 17 +-- .../src/ConvertCWSDExpToMomentum.cpp | 105 +++++++----------- .../algorithms/CollectHB3AExperimentInfo.py | 93 ++++++++-------- 3 files changed, 99 insertions(+), 116 deletions(-) diff --git a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertCWSDExpToMomentum.h b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertCWSDExpToMomentum.h index e7748f0a4869..f1494ac05257 100644 --- a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertCWSDExpToMomentum.h +++ b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertCWSDExpToMomentum.h @@ -6,6 +6,7 @@ #include "MantidAPI/ITableWorkspace.h" #include "MantidAPI/IMDEventWorkspace.h" #include "MantidDataObjects/MDEventInserter.h" +#include "MantidKernel/Matrix.h" namespace Mantid { namespace MDAlgorithms { @@ -39,9 +40,7 @@ class DLLExport ConvertCWSDExpToMomentum : public API::Algorithm { virtual ~ConvertCWSDExpToMomentum(); /// Algorithm's name - virtual const std::string name() const { - return "ConvertCWSDExpToMomentum"; - } + virtual const std::string name() const { return "ConvertCWSDExpToMomentum"; } /// Summary of algorithms purpose virtual const std::string summary() const { @@ -66,8 +65,10 @@ class DLLExport ConvertCWSDExpToMomentum : public API::Algorithm { const detid_t &startdetid, const int runnumber); - void convertToMomentum(const std::vector &detPos, const double &wavelength, - std::vector &qSample); + void convertToMomentum(const std::vector &detPos, + const double &wavelength, + std::vector &qSample, + const Kernel::DblMatrix &rotationMatrix); API::IMDEventWorkspace_sptr createExperimentMDWorkspace(); @@ -80,9 +81,9 @@ class DLLExport ConvertCWSDExpToMomentum : public API::Algorithm { std::vector &vec_detid); void setupTransferMatrix(API::MatrixWorkspace_sptr dataws, - const double &omega, - const double &phi, - const double &chi); + const double &omega, const double &phi, + const double &chi, + Kernel::DblMatrix &rotationMatrix); API::ITableWorkspace_sptr m_expDataTableWS; API::ITableWorkspace_sptr m_detectorListTableWS; diff --git a/Code/Mantid/Framework/MDAlgorithms/src/ConvertCWSDExpToMomentum.cpp b/Code/Mantid/Framework/MDAlgorithms/src/ConvertCWSDExpToMomentum.cpp index 700d51b69ffa..7a5859d7970c 100644 --- a/Code/Mantid/Framework/MDAlgorithms/src/ConvertCWSDExpToMomentum.cpp +++ b/Code/Mantid/Framework/MDAlgorithms/src/ConvertCWSDExpToMomentum.cpp @@ -16,7 +16,7 @@ using namespace Mantid::Geometry; namespace Mantid { namespace MDAlgorithms { - DECLARE_ALGORITHM(ConvertCWSDExpToMomentum) +DECLARE_ALGORITHM(ConvertCWSDExpToMomentum) //---------------------------------------------------------------------------------------------- /** Constructor @@ -42,7 +42,7 @@ void ConvertCWSDExpToMomentum::init() { "Name of table workspace containing all the detectors."); declareProperty(new WorkspaceProperty( - "OutputWorkspace", "", Direction::Output), + "OutputWorkspace", "", Direction::Output), "Name of MDEventWorkspace containing all experimental data."); declareProperty(new ArrayProperty("SourcePosition"), @@ -53,7 +53,6 @@ void ConvertCWSDExpToMomentum::init() { declareProperty(new ArrayProperty("PixelDimension"), "A vector of 8 doubles to determine a cubic pixel's size."); - } /** @@ -87,15 +86,14 @@ ConvertCWSDExpToMomentum::createExperimentMDWorkspace() { parseDetectorTable(vec_detpos, vec_detid); // FIXME - Should set create a solid instrument - m_virtualInstrument = Geometry::ComponentHelper::createVirtualInstrument(m_sourcePos, m_samplePos, - vec_detpos, - vec_detid); - if (!m_virtualInstrument) - throw std::runtime_error("Failed to create virtual instrument."); - g_log.notice() << "[DB] Virtual Instrument has " << m_virtualInstrument->getDetectorIDs().size() + m_virtualInstrument = Geometry::ComponentHelper::createVirtualInstrument( + m_sourcePos, m_samplePos, vec_detpos, vec_detid); + if (!m_virtualInstrument) + throw std::runtime_error("Failed to create virtual instrument."); + g_log.notice() << "[DB] Virtual Instrument has " + << m_virtualInstrument->getDetectorIDs().size() << "Detectors\n"; - // Create workspace // FIXME - Should I give out a better value??? m_extentMins.resize(3, -10.0); @@ -118,15 +116,18 @@ ConvertCWSDExpToMomentum::createExperimentMDWorkspace() { dimensionNames[0] = "Q_sample_x"; dimensionNames[1] = "Q_sample_y"; dimensionNames[2] = "Q_sample_z"; - Mantid::Kernel::SpecialCoordinateSystem coordinateSystem = Mantid::Kernel::QSample; + Mantid::Kernel::SpecialCoordinateSystem coordinateSystem = + Mantid::Kernel::QSample; // Add dimensions - if (m_extentMins.size() != 3 || m_extentMaxs.size() != 3 || m_numBins.size() != 3) - throw std::runtime_error("The range and number of bins are not set up correctly to create MDEventWorkspace."); + if (m_extentMins.size() != 3 || m_extentMaxs.size() != 3 || + m_numBins.size() != 3) + throw std::runtime_error("The range and number of bins are not set up " + "correctly to create MDEventWorkspace."); else for (size_t d = 0; d < 3; ++d) g_log.debug() << "Direction " << d << ", Range = " << m_extentMins[d] - << ", " << m_extentMaxs[d] << "\n"; + << ", " << m_extentMaxs[d] << "\n"; for (size_t i = 0; i < m_nDimensions; ++i) { std::string id = vec_ID[i]; @@ -187,46 +188,19 @@ void ConvertCWSDExpToMomentum::addMDEvents() { return; } -void ConvertCWSDExpToMomentum::setupTransferMatrix(API::MatrixWorkspace_sptr dataws, - const double &omega, - const double &phi, - const double &chi) -{ - // Set up MDWSDescription from MatrixWorkspace - MDWSDescription targetwsdescription; - std::vector minVal(4,-3),maxVal(4,3); - targetwsdescription.setMinMax(minVal,maxVal); - - targetwsdescription.buildFromMatrixWS(dataws, "Q3D", "Direct"); - - // Now it is in LabFrame - MDWSTransform transf; - - std::string frameid = "Q_sample"; - CnvrtToMD::TargetFrame targetframe = transf.getTargetFrame(frameid); - /* - bool foundtarget = transf.findTargetFrame(targetwsdescription); - if (!foundtarget) - throw std::runtime_error("Unable to convert to ..."); - */ - - // SetGoniometer(fakews2d,0,0,0); +void ConvertCWSDExpToMomentum::setupTransferMatrix( + API::MatrixWorkspace_sptr dataws, const double &omega, const double &phi, + const double &chi, Kernel::DblMatrix &rotationMatrix) { // set axis-0 dataws->mutableRun().mutableGoniometer().setRotationAngle(0, omega); // set axis-1 dataws->mutableRun().mutableGoniometer().setRotationAngle(1, phi); // set axis-2 dataws->mutableRun().mutableGoniometer().setRotationAngle(2, chi); - // TS_ASSERT_EQUALS(CnvrtToMD::SampleFrame,Transf.findTargetFrame(TargWSDescription)); - // FIXME - How to define *pLattice? - Geometry::OrientedLattice *pLattice; - dataws->mutableSample().setOrientedLattice(pLattice); - // It should be in sample frame - // FIXME : TS_ASSERT_EQUALS(CnvrtToMD::HKLFrame,Transf.findTargetFrame(TargWSDescription)); - - CnvrtToMD::CoordScaling scaling(CnvrtToMD::NoScaling); - std::vector matrix = - transf.getTransfMatrix(targetwsdescription, "Sample", "NoScaling"); + // Get rotation matrix from Q-lab to Q-sample + rotationMatrix = dataws->run().getGoniometer().getR(); + g_log.notice() << "[DB] Rotation Matrix Rows = " << rotationMatrix.numRows() + << "\n"; return; } @@ -234,12 +208,17 @@ void ConvertCWSDExpToMomentum::setupTransferMatrix(API::MatrixWorkspace_sptr dat void ConvertCWSDExpToMomentum::convertSpiceMatrixToMomentumMDEvents( MatrixWorkspace_sptr dataws, const detid_t &startdetid, const int runnumber) { + // Create transformation matrix from + double omega = dataws->run().getPropertyAsSingleValue("_omega"); + double phi = dataws->run().getPropertyAsSingleValue("_phi"); + double chi = dataws->run().getPropertyAsSingleValue("_chi"); + Kernel::DblMatrix rotationMatrix; + setupTransferMatrix(dataws, omega, phi, chi, rotationMatrix); // Creates a new instance of the MDEventInserter. MDEventWorkspace, 3>::sptr mdws_mdevt_3 = boost::dynamic_pointer_cast, 3> >(m_outputWS); - MDEventInserter, 3>::sptr> inserter( - mdws_mdevt_3); + MDEventInserter, 3>::sptr> inserter(mdws_mdevt_3); size_t numspec = dataws->getNumberHistograms(); for (size_t iws = 0; iws < numspec; ++iws) { @@ -250,12 +229,12 @@ void ConvertCWSDExpToMomentum::convertSpiceMatrixToMomentumMDEvents( // Create the MDEvent Kernel::V3D detpos = dataws->getDetector(iws)->getPos(); std::vector momentum(3); - convertToMomentum(detpos, wavelength, momentum); + convertToMomentum(detpos, wavelength, momentum, rotationMatrix); detid_t detid = dataws->getDetector(iws)->getID() + startdetid; // Insert - inserter.insertMDEvent(static_cast(signal), static_cast(error*error), - static_cast(runnumber), detid, - momentum.data()); + inserter.insertMDEvent( + static_cast(signal), static_cast(error * error), + static_cast(runnumber), detid, momentum.data()); } ExperimentInfo_sptr expinfo = boost::make_shared(); @@ -320,19 +299,22 @@ bool ConvertCWSDExpToMomentum::getInputs(std::string &errmsg) { return (errmsg.size() > 0); } -void ConvertCWSDExpToMomentum::convertToMomentum(const std::vector &detPos, const double &wavelength, std::vector &qSample) -{ - // FIXME - At this stage it is only for test on compiling purpose - return; - +void ConvertCWSDExpToMomentum::convertToMomentum( + const std::vector &detPos, const double &wavelength, + std::vector &qSample, const Kernel::DblMatrix &rotationMatrix) { + // TODO - Use detector position and wavelength/Q to calcualte Q_lab + Kernel::V3D q_lab; + // TODO - Use matrix workspace + Kernel::V3D q_sample = rotationMatrix * q_lab; return; } -API::MatrixWorkspace_sptr ConvertCWSDExpToMomentum::loadSpiceData(const std::string &filename, bool &loaded, std::string &errmsg) -{ +API::MatrixWorkspace_sptr +ConvertCWSDExpToMomentum::loadSpiceData(const std::string &filename, + bool &loaded, std::string &errmsg) { IAlgorithm_sptr loader = createChildAlgorithm("LoadSpiceXML2DDet"); loader->initialize(); loader->setProperty("Filename", filename); @@ -353,7 +335,6 @@ API::MatrixWorkspace_sptr ConvertCWSDExpToMomentum::loadSpiceData(const std::str return dataws; } - void ConvertCWSDExpToMomentum::parseDetectorTable( std::vector &vec_detpos, std::vector &vec_detid) { // Set vectors' sizes diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CollectHB3AExperimentInfo.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CollectHB3AExperimentInfo.py index a5a6510a0f27..c6baa2741eff 100644 --- a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CollectHB3AExperimentInfo.py +++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CollectHB3AExperimentInfo.py @@ -5,11 +5,14 @@ from mantid.kernel import * import os +iColPtNumber = 0 +iCol2Theta = 10 + + class CollectHB3AExperimentInfo(PythonAlgorithm): """ Python algorithm to export sample logs to spread sheet file for VULCAN """ - def category(self): """ Category """ @@ -30,48 +33,26 @@ def PyInit(self): # Format of the input should be self.declareProperty("ExperimentNumber", -1, "Integer for experiment number.") - self.declareProperty(IntArrayProperty("ScanList", [], Direction.Input), - "List of scans of the experiment to be loaded. It cannot be left blank.") + scanlistprop = mantid.kernel.IntArrayProperty("ScanList", []) + self.declareProperty(scanlistprop, "List of scans of the experiment to be loaded. It cannot be left blank") + #self.declareProperty(IntArrayProperty("ScanList", [], Direction.Input), + # "List of scans of the experiment to be loaded. It cannot be left blank.") - self.declareProperty(InArrayProperty("PtLists", [], Direction.Input), + self.declareProperty(mantid.kernel.IntArrayProperty("PtLists", []), "List of Pts to be loaded for scans specified by 'ScanList'. Each scan's Pt. must be started with -1 as a flag. If no Pt. given, then all Pt. are all loaded.") self.declareProperty(FileProperty(name="DataDirectory",defaultValue="",action=FileAction.Directory)) self.declareProperty("GetFileFromServer", False, "Obtain data file directly from neutrons.ornl.gov.") - self.declarerPoperty("Detector2ThetaTolerance", 0.01, "Tolerance of 2 detector's 2theta to consider as being at same position.") - - - ## Input workspace - #self.declareProperty(MatrixWorkspaceProperty("InputWorkspace", "", Direction.Input), - # "Name of data workspace containing sample logs to be exported. ") - - ## Output file name - #self.declareProperty(FileProperty("OutputFilename", "", FileAction.Save, [".txt"]), - # "Name of the output sample environment log file name.") - - ## Sample log names - #self.declareProperty(StringArrayProperty("SampleLogNames", values=[], direction=Direction.Input), - # "Names of sample logs to be exported in a same file.") - - ## Header - #self.declareProperty("WriteHeaderFile", False, "Flag to generate a sample log header file.") + self.declareProperty("Detector2ThetaTolerance", 0.01, "Tolerance of 2 detector's 2theta to consider as being at same position.") - #self.declareProperty("Header", "", "String in the header file.") + tableprop = mantid.api.ITableWorkspaceProperty("OutputWorkspace", "", mantid.kernel.Direction.Output) + self.declareProperty(tableprop, "TableWorkspace for experiment number, scan, file name and starting detector IDs.") - ## Time zone - #timezones = ["UTC", "America/New_York", "Asia/Shanghai", "Australia/Sydney", "Europe/London", "GMT+0",\ - # "Europe/Paris", "Europe/Copenhagen"] + tableprop2 = mantid.api.ITableWorkspaceProperty("DetectorTableWorkspace", "", mantid.kernel.Direction.Output) + self.declareProperty(tableprop2, "TableWorkspace for detector Id and information.") - #description = "Sample logs recorded in NeXus files (in SNS) are in UTC time. TimeZone " + \ - # "can allow the algorithm to output the log with local time." - #self.declareProperty("TimeZone", "America/New_York", StringListValidator(timezones), description) - - ## Log time tolerance - #self.declareProperty("TimeTolerance", 0.01, - # "If any 2 log entries with log times within the time tolerance, " + \ - # "they will be recorded in one line. Unit is second. ") return @@ -83,7 +64,7 @@ def PyExec(self): self._getProperties() # Get scan list - self._getExpScanPtLDict() + self._getExpScanPtDict() # Get 2theta-(scan-pt) dictionary self._getDetectorPositionScanPtDict() @@ -100,14 +81,15 @@ def PyExec(self): def _getExpScanPtDict(self): """ Get the scan-Pt. dictionary - """ - for iscan in self._scanList: + """ + for iscan in xrange(len(self._scanList)): # Loop over scan number scan = self._scanList[iscan] # Load data file - spicefilename = 'HB3A_Exp%04d_Scan%4d.dat' % (self._expNo, scan) + spicefilename = os.path.join(self._dataDir, 'HB3A_exp%04d_scan%04d.dat' % (self._expNumber, scan)) spicetablews = self._loadSpiceFile(spicefilename) + self._spiceTableDict[scan] = spicetablews if spicetablews is None: self.glog.warning("Unable to access Exp %d Scan %d's SPICE file %s." % (self._expNo, scan, spicefilename)) @@ -128,10 +110,11 @@ def _getExpScanPtDict(self): def _getProperties(self): """ Get properties from user input """ - self._expNumber = getProperty("ExperimentNumer").value() - self._scanList = getProperty("ScanList").value() - rawptlist = self.getProperty("PtLists").value() - self._tol2Theta = getProperty("2ThetaTolerance").value() + self._expNumber = self.getProperty("ExperimentNumber").value + self._scanList = self.getProperty("ScanList").value + rawptlist = self.getProperty("PtLists").value + self._tol2Theta = self.getProperty("Detector2ThetaTolerance").value + self._dataDir = self.getProperty("DataDirectory").value # process Pt number self._ptListList = [] @@ -148,9 +131,14 @@ def _getProperties(self): # ENDFOR # check - if len(self._ptListList) != self._scanList: + if len(self._ptListList) != len(self._scanList): raise RuntimeError("Number of sub Pt. list is not equal to number of scans.") + # Define class variable Scan-Pt. dictionary: Key scan number, Value list of pt numbers + self._scanPtDict = {} + + self._spiceTableDict = {} + return @@ -164,7 +152,15 @@ def _getDetectorPositionScanPtDict(self): # Scan every row of every SPICE table to get a dictionary with key as 2theta value for scannumber in self._spiceTableDict.keys(): spicetable = self._spiceTableDict[scannumber] - for irow in xrange(len(spicetable.rowCount())): + + # check column names + colnames = spicetable.getColumnNames() + if colnames[iColPtNumber].lower().startswith('pt') is False or \ + colnames[iCol2Theta].lower().startswith('2theta') is False: + raise NotImplementedError("Colulmn %d is not Pt. but %s; OR Column %d is not 2theta, but %s." % ( + iColPtNumber, colnames[iColPtNumber], iCol2Theta, colnames[iCol2Theta])) + + for irow in xrange(spicetable.rowCount()): ptnumber = spicetable.cell(irow, iColPtNumber) twotheta = spicetable.cell(irow, iCol2Theta) if self._2thetaScanPtDict.has_key(twotheta) is False: @@ -180,7 +176,7 @@ def _getDetectorPositionScanPtDict(self): twotheta_prev = twothetalist[0] for itt in xrange(1, len(twothetalist)): twotheta_curr = twothetalist[itt] - if twotheta_curr - twotheta_prev < self._2thetaTol: + if twotheta_curr - twotheta_prev < self._tol2Theta: # two keys (2theta) are close enough, combine then self._2thetaScanPtDict[twotheta_prev].extend(self._2thetaScanPtDict[twotheta_curr][:]) del self._2thetaScanPtDict[twotehta] @@ -198,9 +194,11 @@ def _collectPixelsPositions(self): for twotheta in sorted(self._2thetaScanPtDict.keys()): if len(self._2thetaScanPtDict[twotheta]) == 0: raise RuntimeError("Logic error to have empty list.") + else: + self.log().notice("[DB] Process pixels of detector centered at 2theta = %.5f." % (twotheta)) # Load detector counts file (.xml) - scannumber, ptnumber = self._2thetaScanPtDict[twotehta] + scannumber, ptnumber = self._2thetaScanPtDict[twotheta] dataws = self._loadHB3ADetCountFile(scannumber, ptnumber) self._detStartID[twotheta] = self._newDetID @@ -237,7 +235,10 @@ def _getAllPtFromTable(self, spicetablews): def _loadSpiceFile(self, spicefilename): """ Load SPICE file """ - spicetablews, spicematrixws = api.LoadSpiceAscii(Filename=spicefilename) + outwsname = os.path.basename(spicefilename).split(".")[0] + self.log().notice("Loading SPICE file %s." % (spicefilename)) + spicetablews, spicematrixws = api.LoadSpiceAscii(Filename=spicefilename, OutputWorkspace=outwsname) + self.log().notice("[DB] SPICE table workspace has %d rows."%(spicetablews.rowCount())) return spicetablews