19#include "boost/math/special_functions/round.hpp"
43 "An input PeaksWorkspace.");
45 "Select the directory and base name for the output files.");
46 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
47 mustBePositive->setLower(0.0);
48 declareProperty(
"ScalePeaks", 1.0, mustBePositive,
"Multiply FSQ and sig(FSQ) by scaleFactor");
49 declareProperty(
"MinDSpacing", 0.0,
"Minimum d-spacing (Angstroms)");
50 declareProperty(
"MinWavelength", 0.0,
"Minimum wavelength (Angstroms)");
51 declareProperty(
"MaxWavelength",
EMPTY_DBL(),
"Maximum wavelength (Angstroms)");
52 std::vector<std::string> histoTypes{
"Bank",
"RunNumber",
"Both Bank and RunNumber"};
53 declareProperty(
"SortFilesBy", histoTypes[0], std::make_shared<StringListValidator>(histoTypes),
54 "Sort into files by bank(default), run number or both.");
55 declareProperty(
"MinIsigI",
EMPTY_DBL(), mustBePositive,
"The minimum I/sig(I) ratio");
56 declareProperty(
"WidthBorder",
EMPTY_INT(),
"Width of border of detectors");
57 declareProperty(
"MinIntensity",
EMPTY_DBL(), mustBePositive,
"The minimum Intensity");
58 declareProperty(
"UseDetScale",
false,
59 "Scale intensity and sigI by scale "
60 "factor of detector if set in "
62 "If false, no change (default).");
64 "Comma deliminated string of bank numbers to exclude for example 1,2,5");
65 declareProperty(
"LaueScaleFormat",
false,
"New format for Lauescale");
67 declareProperty(
"CrystalSystem", m_typeList[0], std::make_shared<Kernel::StringListValidator>(m_typeList),
68 "The conventional cell type to use");
70 declareProperty(
"Centering", m_centeringList[0], std::make_shared<Kernel::StringListValidator>(m_centeringList),
71 "The centering for the conventional cell");
80 Poco::Path path(filename);
81 std::string basename = path.getBaseName();
92 std::string cellType =
getProperty(
"CrystalSystem");
94 auto cellNo = 1 + std::distance(
m_typeList.begin(), iter);
100 int oldSequence = -1;
103 std::ostringstream ss;
106 std::vector<std::pair<std::string, bool>> criteria;
107 if (type.compare(0, 2,
"Ba") == 0)
108 criteria.emplace_back(
"BankName",
true);
109 else if (type.compare(0, 2,
"Ru") == 0)
110 criteria.emplace_back(
"RunNumber",
true);
112 criteria.emplace_back(
"RunNumber",
true);
113 criteria.emplace_back(
"BankName",
true);
115 criteria.emplace_back(
"h",
true);
116 criteria.emplace_back(
"k",
true);
117 criteria.emplace_back(
"l",
true);
120 std::vector<Peak> peaks =
ws->getPeaks();
128 if (convention ==
"Crystallography")
132 auto inst =
ws->getInstrument();
135 type =
"Both Bank and RunNumber";
136 if (!
ws->sample().hasOrientedLattice()) {
138 const std::string fft(
"FindUBUsingIndexedPeaks");
140 findUB->initialize();
142 findUB->executeAsChildAlg();
144 if (!
ws->sample().hasOrientedLattice()) {
145 g_log.
notice(std::string(
"Could not find UB for ") + std::string(
ws->getName()));
146 throw std::invalid_argument(std::string(
"Could not find UB for ") + std::string(
ws->getName()));
149 lattice =
ws->sample().getOrientedLattice();
152 std::vector<int> numPeaks;
154 std::vector<double> maxLamVec;
155 std::vector<double> minLamVec;
156 std::vector<double> sumLamVec;
157 std::vector<double> minDVec;
162 for (
int wi = 0; wi <
ws->getNumberPeaks(); wi++) {
167 if (intensity == 0.0 || !(std::isfinite(sigI)))
169 if (minIsigI !=
EMPTY_DBL() && intensity < std::abs(minIsigI * sigI))
176 p.
getCol() > (nCols - widthBorder) || p.
getRow() > (nRows - widthBorder)))
179 bankName.erase(remove_if(bankName.begin(), bankName.end(), std::not_fn(::isdigit)), bankName.end());
180 if (type.compare(0, 2,
"Ba") == 0) {
184 std::vector<std::string> notBanks =
getProperty(
"EliminateBankNumbers");
185 if (std::find(notBanks.begin(), notBanks.end(), bankName) != notBanks.end())
188 if (inst->hasParameter(
"detScale" + bankName)) {
189 double correc =
static_cast<double>(inst->getNumberParameter(
"detScale" + bankName)[0]);
194 if (minIntensity !=
EMPTY_DBL() && intensity < minIntensity)
203 if (sequence != oldSequence) {
204 oldSequence = sequence;
205 numPeaks.emplace_back(
count);
206 maxLamVec.emplace_back(maxLam);
207 minLamVec.emplace_back(minLam);
208 sumLamVec.emplace_back(sumLam);
209 minDVec.emplace_back(minD);
225 numPeaks.emplace_back(
count);
226 maxLamVec.emplace_back(maxLam);
227 minLamVec.emplace_back(minLam);
228 sumLamVec.emplace_back(sumLam);
229 minDVec.emplace_back(minD);
232 for (
int wi = 0; wi <
ws->getNumberPeaks(); wi++) {
237 if (intensity == 0.0 || !(std::isfinite(sigI)))
239 if (minIsigI !=
EMPTY_DBL() && intensity < std::abs(minIsigI * sigI))
246 p.
getCol() > (nCols - widthBorder) || p.
getRow() > (nRows - widthBorder)))
249 bankName.erase(remove_if(bankName.begin(), bankName.end(), std::not_fn(::isdigit)), bankName.end());
250 if (type.compare(0, 2,
"Ba") == 0) {
254 std::vector<std::string> notBanks =
getProperty(
"EliminateBankNumbers");
255 if (std::find(notBanks.begin(), notBanks.end(), bankName) != notBanks.end())
258 if (inst->hasParameter(
"detScale" + bankName)) {
259 double correc =
static_cast<double>(inst->getNumberParameter(
"detScale" + bankName)[0]);
264 if (minIntensity !=
EMPTY_DBL() && intensity < minIntensity)
274 if (sequence != oldSequence) {
275 oldSequence = sequence;
277 out <<
"END-OF-REFLECTION-DATA\n";
278 out <<
"HARMONICS DATA 0 REFLECTIONS\n";
279 out <<
"END-OF-FILE\n";
286 ss << std::setw(3) << std::setfill(
'0') << sequenceNo;
290 path.append(basename + ss.str());
292 path.setExtension(
"geasc");
293 Poco::File fileobj(path);
294 out.open(path.toString().c_str(), std::ios::out);
297 out << basename <<
"\n";
298 out <<
"CRYS " << basename.substr(0, 6) <<
"\n";
299 out <<
"FIDX 1.00000 1.00000 1.00000 1.00000 "
301 out <<
"FIDY 1.00000 1.00000 1.00000 1.00000 "
303 out <<
"OMEG 1.00000 1.00000 1.00000 1.00000 "
305 out <<
"CELL " << std::setw(11) << std::setprecision(4) << 1.0 / lattice.
a() << std::setw(12)
306 << std::setprecision(4) << 1.0 / lattice.
b() << std::setw(12) << std::setprecision(4) << 1.0 / lattice.
c()
307 << std::setw(9) << boost::math::iround(lattice.
alpha()) << std::setw(9)
308 << boost::math::iround(lattice.
beta()) << std::setw(9) << boost::math::iround(lattice.
gamma()) <<
'\n';
310 out <<
"SYST " << cellNo <<
" " << centerNo <<
" 1 3"
314 out <<
"IBOX 1 1 1 1 1"
319 double phi = angles[2];
320 double chi = angles[1];
321 double omega = angles[0];
323 out <<
"PHIS " << std::setw(11) << std::setprecision(4) << phi << std::setw(12) << std::setprecision(4) << chi
324 << std::setw(12) << std::setprecision(4) << omega <<
"\n";
327 out << std::setprecision(1) << std::fixed << sumLamVec[sequenceNo] / numPeaks[sequenceNo] <<
" "
328 << minLamVec[sequenceNo] <<
" " << maxLamVec[sequenceNo] <<
"\n";
331 out << std::setprecision(2) << std::fixed << minDVec[sequenceNo] <<
"\n";
335 out <<
"RADI " << std::setprecision(0) << std::fixed << L2 <<
"\n";
338 out <<
"XC_S 0.00000 0.00000 0.00000 0.00000 "
340 out <<
"YC_S 0.00000 0.00000 0.00000 0.00000 "
342 out <<
"WC_S 0.00000 0.00000 0.00000 0.00000 "
346 out <<
"TWIS 0.00000 0.00000 0.00000 0.00000 "
347 "0.00000 0.00000 \n";
348 out <<
"TILT 0.00000 0.00000 0.00000 0.00000 "
349 "0.00000 0.00000 \n";
350 out <<
"BULG 0.00000 0.00000 0.00000 0.00000 "
351 "0.00000 0.00000 \n";
352 out <<
"CTOF " << L2 <<
"\n";
353 out <<
"YSCA 1.00000 1.00000 1.00000 1.00000 "
355 out <<
"CRAT 1.00000 1.00000 1.00000 1.00000 "
359 out << minIntensity <<
"\n";
364 out << numPeaks[sequenceNo]
367 out <<
" 0 0 0 0 0 0 0 0 "
370 out <<
"LAMH " << numPeaks[sequenceNo]
371 <<
" 0 0 0 0 0 0 0 0 "
373 out <<
" 0 0 0 0 0 0\n";
378 out <<
"NSPT " << numPeaks[sequenceNo] <<
" 0 0 0 0"
380 out <<
"NODH " << numPeaks[sequenceNo] <<
" 0 0 0 0 0 0 0 0 0\n"
384 out <<
"REFLECTION DATA " << numPeaks[sequenceNo] <<
" REFLECTIONS"
398 out << std::setw(10) << std::fixed << std::setprecision(5) << (p.
getCol() - 127.5) * 150.0 / 256.0;
399 out << std::setw(10) << std::fixed << std::setprecision(5) << (p.
getRow() - 127.5) * 150.0 / 256.0 <<
"\n";
401 out << std::setw(10) << std::fixed << std::setprecision(5) <<
lambda;
404 out <<
" 1 0 0 0 0 0 0 0 0.0 0 ";
406 out << std::setw(10) << std::fixed << std::setprecision(5) << dsp * dsp * 0.25 <<
"\n";
409 out << std::setw(10) << std::fixed << std::setprecision(5) << 0.5 * scattering;
414 double ckIntensity = scaleFactor * intensity;
415 if (ckIntensity > 999999999.985)
416 g_log.
warning() <<
"Scaled intensity, " << ckIntensity <<
" is too large for format. Decrease ScalePeaks.\n";
420 out <<
" -9999 -9999 -9999 -9999 -9999 \n";
423 out << std::setw(10) <<
Utils::round(scaleFactor * sigI);
426 out <<
" -9999 -9999 -9999 -9999 -9999 \n";
428 out <<
" -9999 -9999 -9999 -9999 -9999 \n";
429 out << std::setw(10) <<
Utils::round(scaleFactor * sigI);
430 out <<
" -9999 -9999 -9999 -9999 -9999 * ";
437 out <<
" -9999 -9999 -9999 -9999 -9999 \n";
442 out <<
" -9999 -9999 -9999 -9999 -9999 \n";
444 out <<
" -9999 -9999 -9999 -9999 -9999 \n";
446 out <<
" -9999 -9999 -9999 -9999 -9999 * ";
453 out <<
"END-OF-REFLECTION-DATA\n";
454 out <<
"HARMONICS DATA 0 REFLECTIONS\n";
455 out <<
"END-OF-FILE\n";
461 if (bankName ==
"None")
463 std::shared_ptr<const IComponent> parent =
ws->getInstrument()->getComponentByName(bankName);
466 if (parent->type() ==
"RectangularDetector") {
467 std::shared_ptr<const RectangularDetector> RDet = std::dynamic_pointer_cast<const RectangularDetector>(parent);
469 nCols = RDet->xpixels();
470 nRows = RDet->ypixels();
472 std::vector<Geometry::IComponent_const_sptr> children;
473 std::shared_ptr<const Geometry::ICompAssembly> asmb =
474 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(parent);
475 asmb->getChildren(children,
false);
476 std::shared_ptr<const Geometry::ICompAssembly> asmb2 =
477 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[0]);
478 std::vector<Geometry::IComponent_const_sptr> grandchildren;
479 asmb2->getChildren(grandchildren,
false);
480 nRows =
static_cast<int>(grandchildren.size());
481 nCols =
static_cast<int>(children.size());
#define DECLARE_ALGORITHM(classname)
const std::vector< double > * lambda
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
@ Save
to specify a file to write to, the file may or may not exist
A property class for workspaces.
Save a PeaksWorkspace to a lauenorm format http://www.ccp4.ac.uk/cvs/viewvc.cgi/laue/doc/lauenorm....
void sizeBanks(const std::string &bankName, int &nCols, int &nRows)
const std::vector< std::string > m_centeringList
const std::vector< std::string > m_typeList
DataObjects::PeaksWorkspace_sptr ws
void exec() override
Run the algorithm.
double getL() const override
Get the L index of the peak.
double getIntensity() const override
Return the integrated peak intensity.
double getSigmaIntensity() const override
Return the error on the integrated peak intensity.
double getK() const override
Get the K index of the peak.
int getRunNumber() const override
Return the run number this peak was measured at.
double getH() const override
Get the H index of the peak.
Mantid::Kernel::Matrix< double > getGoniometerMatrix() const override
Get the goniometer rotation matrix at which this peak was measured.
Structure describing a single-crystal peak.
int getCol() const override
For RectangularDetectors only, returns the column (x) of the pixel of the detector or -1 if not found...
double getDSpacing() const override
Calculate the d-spacing of the peak, in 1/Angstroms
std::string getBankName() const
Find the name of the bank that is the parent of the detector.
int getRow() const override
For RectangularDetectors only, returns the row (y) of the pixel of the detector or -1 if not found.
double getWavelength() const override
Calculate the neutron wavelength (in angstroms) at the peak (Note for inelastic scattering - it is th...
int getDetectorID() const
Get the ID of the detector at the center of the peak
double getScattering() const override
Calculate the scattering angle of the peak
Class to represent a particular goniometer setting, which is described by the rotation matrix.
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
Class to implement UB matrix.
double alpha() const
Get lattice parameter.
double a(int nd) const
Get lattice parameter a1-a3 as function of index (0-2)
double c() const
Get lattice parameter.
double beta() const
Get lattice parameter.
double b() const
Get lattice parameter.
double gamma() const
Get lattice parameter.
Support for a property that holds an array of values.
void notice(const std::string &msg)
Logs at notice level.
void warning(const std::string &msg)
Logs at warning level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
int convert(const std::string &A, T &out)
Convert a string into a number.
long round(double x)
Custom rounding method for a double->long because none is portable in C++ (!)
A namespace containing physical constants that are required by algorithms and unit routines.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.