Mantid
Loading...
Searching...
No Matches
OptimizeLatticeForCellType.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
12#include "MantidAPI/Sample.h"
25#include <fstream>
26
27using namespace Mantid::Geometry;
28
29namespace Mantid::Crystal {
30// Register the class into the algorithm factory
31DECLARE_ALGORITHM(OptimizeLatticeForCellType)
32
33using namespace Kernel;
34using namespace API;
35using std::size_t;
36using namespace DataObjects;
37
41
42 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
43 "An input PeaksWorkspace with an instrument.");
44 std::vector<std::string> cellTypes;
45 cellTypes.emplace_back(ReducedCell::CUBIC());
46 cellTypes.emplace_back(ReducedCell::TETRAGONAL());
47 cellTypes.emplace_back(ReducedCell::ORTHORHOMBIC());
48 cellTypes.emplace_back(ReducedCell::HEXAGONAL());
49 cellTypes.emplace_back(ReducedCell::RHOMBOHEDRAL());
50 cellTypes.emplace_back(ReducedCell::MONOCLINIC());
51 cellTypes.emplace_back(ReducedCell::TRICLINIC());
52 declareProperty("CellType", cellTypes[0], std::make_shared<StringListValidator>(cellTypes), "Select the cell type.");
53 declareProperty("Apply", false, "Re-index the peaks");
54 declareProperty("PerRun", false, "Make per run orientation matrices");
55 declareProperty("Tolerance", 0.12, "Indexing Tolerance");
56 declareProperty("EdgePixels", 0, "Remove peaks that are at pixels this close to edge. ");
57 declareProperty(std::make_unique<PropertyWithValue<double>>("OutputChi2", 0.0, Direction::Output),
58 "Returns the goodness of the fit");
59 declareProperty(std::make_unique<FileProperty>("OutputDirectory", ".", FileProperty::Directory),
60 "The directory where the per run peaks files and orientation matrices "
61 "will be written.");
62
63 // Disable default gsl error handler (which is to call abort!)
64 gsl_set_error_handler_off();
65}
66
73 bool apply = this->getProperty("Apply");
74 bool perRun = this->getProperty("PerRun");
75 double tolerance = this->getProperty("Tolerance");
76 int edge = this->getProperty("EdgePixels");
77 std::string cell_type = getProperty("CellType");
78 IPeaksWorkspace_sptr ws = getProperty("PeaksWorkspace");
79
80 std::vector<int> badPeaks;
81 std::vector<IPeaksWorkspace_sptr> runWS;
82 if (edge > 0)
83 if (auto pw = std::dynamic_pointer_cast<PeaksWorkspace>(ws)) {
84 Geometry::Instrument_const_sptr inst = ws->getInstrument();
85 for (int i = int(pw->getNumberPeaks()) - 1; i >= 0; --i) {
86 const std::vector<Peak> &peaks = pw->getPeaks();
87 if (edgePixel(inst, peaks[i].getBankName(), peaks[i].getCol(), peaks[i].getRow(), edge)) {
88 badPeaks.emplace_back(i);
89 }
90 }
91 pw->removePeaks(std::move(badPeaks));
92 }
93 runWS.emplace_back(ws);
94
95 if (perRun) {
96 std::vector<std::pair<std::string, bool>> criteria;
97 // Sort by run number
98 criteria.emplace_back("runnumber", true);
99 ws->sort(criteria);
100 int run = 0;
101 int count = 0;
102 for (int i = 0; i < ws->getNumberPeaks(); i++) {
103 IPeak &peak = ws->getPeak(i);
104 if (peak.getRunNumber() != run) {
105 count++; // first entry in runWS is input workspace
106 auto cloneWS = std::dynamic_pointer_cast<IPeaksWorkspace>(WorkspaceFactory::Instance().createPeaks(ws->id()));
107 cloneWS->copyExperimentInfoFrom(ws.get());
108 runWS.emplace_back(cloneWS);
109 runWS[count]->addPeak(peak);
110 run = peak.getRunNumber();
111 AnalysisDataService::Instance().addOrReplace(std::to_string(run) + ws->getName(), runWS[count]);
112 } else {
113 runWS[count]->addPeak(peak);
114 }
115 }
116 }
117 // finally do the optimization
118 for (auto &i_run : runWS) {
119 IPeaksWorkspace_sptr peakWS(i_run->clone());
120 AnalysisDataService::Instance().addOrReplace("_peaks", peakWS);
121 const DblMatrix UB = peakWS->sample().getOrientedLattice().getUB();
122 auto ol = peakWS->sample().getOrientedLattice();
123 DblMatrix modUB = peakWS->mutableSample().getOrientedLattice().getModUB();
124 int maxOrder = peakWS->mutableSample().getOrientedLattice().getMaxOrder();
125 bool crossTerms = peakWS->mutableSample().getOrientedLattice().getCrossTerm();
126 std::vector<double> lat(6);
128
129 API::ILatticeFunction_sptr latticeFunction = getLatticeFunction(cell_type, peakWS->sample().getOrientedLattice());
130
131 IAlgorithm_sptr fit_alg;
132 try {
133 fit_alg = createChildAlgorithm("Fit", -1, -1, false);
134 } catch (Exception::NotFoundError &) {
135 g_log.error("Can't locate Fit algorithm");
136 throw;
137 }
138
139 fit_alg->setProperty("Function", std::static_pointer_cast<IFunction>(latticeFunction));
140 fit_alg->setProperty("Ties", "ZeroShift=0.0");
141 fit_alg->setProperty("InputWorkspace", peakWS);
142 fit_alg->setProperty("CostFunction", "Unweighted least squares");
143 fit_alg->setProperty("CreateOutput", true);
144 fit_alg->executeAsChildAlg();
145
146 double chisq = fit_alg->getProperty("OutputChi2overDoF");
147 Geometry::UnitCell refinedCell = latticeFunction->getUnitCell();
148
149 IAlgorithm_sptr ub_alg;
150 try {
151 ub_alg = createChildAlgorithm("CalculateUMatrix", -1, -1, false);
152 } catch (Exception::NotFoundError &) {
153 g_log.error("Can't locate CalculateUMatrix algorithm");
154 throw;
155 }
156
157 ub_alg->setProperty("PeaksWorkspace", peakWS);
158 ub_alg->setProperty("a", refinedCell.a());
159 ub_alg->setProperty("b", refinedCell.b());
160 ub_alg->setProperty("c", refinedCell.c());
161 ub_alg->setProperty("alpha", refinedCell.alpha());
162 ub_alg->setProperty("beta", refinedCell.beta());
163 ub_alg->setProperty("gamma", refinedCell.gamma());
164 ub_alg->executeAsChildAlg();
165 DblMatrix UBnew = peakWS->mutableSample().getOrientedLattice().getUB();
166 auto o_lattice = std::make_unique<OrientedLattice>();
167 o_lattice->setUB(UBnew);
168 if (maxOrder > 0) {
169 o_lattice->setModUB(modUB);
170 o_lattice->setMaxOrder(maxOrder);
171 o_lattice->setCrossTerm(crossTerms);
172 }
173 o_lattice->set(refinedCell.a(), refinedCell.b(), refinedCell.c(), refinedCell.alpha(), refinedCell.beta(),
174 refinedCell.gamma());
175 o_lattice->setError(refinedCell.errora(), refinedCell.errorb(), refinedCell.errorc(), refinedCell.erroralpha(),
176 refinedCell.errorbeta(), refinedCell.errorgamma());
177
178 // Show the modified lattice parameters
179 g_log.notice() << i_run->getName() << " " << *o_lattice << "\n";
180
181 i_run->mutableSample().setOrientedLattice(std::move(o_lattice));
182
183 setProperty("OutputChi2", chisq);
184
185 if (apply) {
186 // Reindex peaks with new UB
187 auto alg = createChildAlgorithm("IndexPeaks");
188 alg->setPropertyValue("PeaksWorkspace", i_run->getName());
189 alg->setProperty("Tolerance", tolerance);
190 alg->executeAsChildAlg();
191 }
192 AnalysisDataService::Instance().remove("_peaks");
193 if (perRun) {
194 std::string outputdir = getProperty("OutputDirectory");
195 if (outputdir.back() != '/')
196 outputdir += "/";
197 // Save Peaks
198 auto savePks_alg = createChildAlgorithm("SaveIsawPeaks");
199 savePks_alg->setPropertyValue("InputWorkspace", i_run->getName());
200 savePks_alg->setProperty("Filename", outputdir + "ls" + i_run->getName() + ".integrate");
201 savePks_alg->executeAsChildAlg();
202 g_log.notice() << "See output file: " << outputdir + "ls" + i_run->getName() + ".integrate"
203 << "\n";
204 // Save UB
205 auto saveUB_alg = createChildAlgorithm("SaveIsawUB");
206 saveUB_alg->setPropertyValue("InputWorkspace", i_run->getName());
207 saveUB_alg->setProperty("Filename", outputdir + "ls" + i_run->getName() + ".mat");
208 saveUB_alg->executeAsChildAlg();
209 // Show the names of files written
210 g_log.notice() << "See output file: " << outputdir + "ls" + i_run->getName() + ".mat"
211 << "\n";
212 }
213 }
214}
215//-----------------------------------------------------------------------------------------
222 const UnitCell &cell) const {
223 std::ostringstream fun_str;
224 fun_str << "name=LatticeFunction,LatticeSystem=" << cellType;
225
226 API::IFunction_sptr rawFunction = API::FunctionFactory::Instance().createInitialized(fun_str.str());
227 API::ILatticeFunction_sptr latticeFunction = std::dynamic_pointer_cast<API::ILatticeFunction>(rawFunction);
228 if (latticeFunction) {
229 latticeFunction->setUnitCell(cell);
230 }
231
232 return latticeFunction;
233}
234
235} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const bool crossTerms
Definition: IndexPeaks.cpp:57
const int maxOrder
Definition: IndexPeaks.cpp:55
int count
counter
Definition: Matrix.cpp:37
double tolerance
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
@ Directory
to specify a directory that must exist
Definition: FileProperty.h:54
A property class for workspaces.
void exec() override
Executes the algorithm.
API::ILatticeFunction_sptr getLatticeFunction(const std::string &cellType, const Geometry::UnitCell &cell) const
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual int getRunNumber() const =0
static bool GetLatticeParameters(const Kernel::DblMatrix &UB, std::vector< double > &lattice_par)
Get the lattice parameters for the specified orientation matrix.
static const std::string HEXAGONAL()
Definition: ReducedCell.h:53
static const std::string MONOCLINIC()
Definition: ReducedCell.h:57
static const std::string RHOMBOHEDRAL()
Definition: ReducedCell.h:54
static const std::string CUBIC()
Definition: ReducedCell.h:52
static const std::string TRICLINIC()
Definition: ReducedCell.h:58
static const std::string TETRAGONAL()
Definition: ReducedCell.h:55
static const std::string ORTHORHOMBIC()
Definition: ReducedCell.h:56
Class to implement unit cell of crystals.
Definition: UnitCell.h:44
double alpha() const
Get lattice parameter.
Definition: UnitCell.cpp:133
double erroralpha(const int angleunit=angDegrees) const
Get lattice parameter error.
Definition: UnitCell.cpp:225
double a(int nd) const
Get lattice parameter a1-a3 as function of index (0-2)
Definition: UnitCell.cpp:94
double c() const
Get lattice parameter.
Definition: UnitCell.cpp:128
double errorgamma(const int angleunit=angDegrees) const
Get lattice parameter error.
Definition: UnitCell.cpp:252
double errorbeta(const int angleunit=angDegrees) const
Get lattice parameter error.
Definition: UnitCell.cpp:239
double beta() const
Get lattice parameter.
Definition: UnitCell.cpp:138
double b() const
Get lattice parameter.
Definition: UnitCell.cpp:123
double errorc() const
Get lattice parameter error.
Definition: UnitCell.cpp:218
double gamma() const
Get lattice parameter.
Definition: UnitCell.cpp:143
double errora() const
Get lattice parameter error.
Definition: UnitCell.cpp:208
double errorb() const
Get lattice parameter error.
Definition: UnitCell.cpp:213
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< ILatticeFunction > ILatticeFunction_sptr
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
MANTID_GEOMETRY_DLL bool edgePixel(const Geometry::Instrument_const_sptr &inst, const std::string &bankName, int col, int row, int Edge)
Function to find peaks near detector edge.
Definition: EdgePixel.cpp:22
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Output
An output workspace.
Definition: Property.h:54