本文整理汇总了C++中PeaksWorkspace_sptr类的典型用法代码示例。如果您正苦于以下问题:C++ PeaksWorkspace_sptr类的具体用法?C++ PeaksWorkspace_sptr怎么用?C++ PeaksWorkspace_sptr使用的例子?那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了PeaksWorkspace_sptr类的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: getParameter
/**
* Updates the map from run number to GoniometerMatrix
*
* @param Peaks The PeaksWorkspace whose peaks contain the run numbers
* along with the corresponding GoniometerMatrix
*
* @param OptRuns A '/' separated "list" of run numbers to include in the
* map. This string must also start and end with a '/'
*
* @param Res The resultant map.
*/
void PeakHKLErrors::getRun2MatMap(
PeaksWorkspace_sptr &Peaks, const std::string &OptRuns,
std::map<int, Mantid::Kernel::Matrix<double>> &Res) const {
for (int i = 0; i < Peaks->getNumberPeaks(); ++i) {
Geometry::IPeak &peak_old = Peaks->getPeak(i);
int runNum = peak_old.getRunNumber();
std::string runNumStr = std::to_string(runNum);
size_t N = OptRuns.find("/" + runNumStr + "/");
if (N < OptRuns.size()) {
double chi =
getParameter("chi" + boost::lexical_cast<std::string>(runNumStr));
double phi =
getParameter("phi" + boost::lexical_cast<std::string>(runNumStr));
double omega =
getParameter("omega" + boost::lexical_cast<std::string>(runNumStr));
Mantid::Geometry::Goniometer uniGonio;
uniGonio.makeUniversalGoniometer();
uniGonio.setRotationAngle("phi", phi);
uniGonio.setRotationAngle("chi", chi);
uniGonio.setRotationAngle("omega", omega);
Res[runNum] = uniGonio.getR();
}
}
}
示例2: fillPossibleHKLsUsingPeaksWorkspace
/// Fills possibleHKLs with all HKLs from the supplied PeaksWorkspace.
void PredictPeaks::fillPossibleHKLsUsingPeaksWorkspace(
const PeaksWorkspace_sptr &peaksWorkspace,
std::vector<V3D> &possibleHKLs) const {
possibleHKLs.clear();
possibleHKLs.reserve(peaksWorkspace->getNumberPeaks());
bool roundHKL = getProperty("RoundHKL");
/* Q is at the end multiplied with the factor determined in the
* constructor (-1 for crystallography, 1 otherwise). So to avoid
* "flippling HKLs" when it's not required, the HKLs of the input
* workspace are also multiplied by the factor that is appropriate
* for the convention stored in the workspace.
*/
double peaks_q_convention_factor =
get_factor_for_q_convention(peaksWorkspace->getConvention());
for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
IPeak &p = peaksWorkspace->getPeak(i);
// Get HKL from that peak
V3D hkl = p.getHKL() * peaks_q_convention_factor;
if (roundHKL)
hkl.round();
possibleHKLs.push_back(hkl);
} // for each hkl in the workspace
}
示例3: getProperty
/// Returns a PeaksWorkspace which is either the input workspace or a clone.
PeaksWorkspace_sptr SortHKL::getOutputPeaksWorkspace(
const PeaksWorkspace_sptr &inputPeaksWorkspace) const {
PeaksWorkspace_sptr outputPeaksWorkspace = getProperty("OutputWorkspace");
if (outputPeaksWorkspace != inputPeaksWorkspace) {
outputPeaksWorkspace.reset(inputPeaksWorkspace->clone().release());
}
return outputPeaksWorkspace;
}
示例4: optLattice
/**
@param inname Name of workspace containing peaks
@param params optimized cell parameters
@param out residuals from optimization
*/
void OptimizeLatticeForCellType::optLattice(std::string inname,
std::vector<double> ¶ms,
double *out) {
PeaksWorkspace_sptr ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
AnalysisDataService::Instance().retrieve(inname));
const std::vector<Peak> &peaks = ws->getPeaks();
size_t n_peaks = ws->getNumberPeaks();
std::vector<V3D> q_vector;
std::vector<V3D> hkl_vector;
for (size_t i = 0; i < params.size(); i++)
params[i] = std::abs(params[i]);
for (size_t i = 0; i < n_peaks; i++) {
q_vector.push_back(peaks[i].getQSampleFrame());
hkl_vector.push_back(peaks[i].getHKL());
}
Mantid::API::IAlgorithm_sptr alg = createChildAlgorithm("CalculateUMatrix");
alg->setPropertyValue("PeaksWorkspace", inname);
alg->setProperty("a", params[0]);
alg->setProperty("b", params[1]);
alg->setProperty("c", params[2]);
alg->setProperty("alpha", params[3]);
alg->setProperty("beta", params[4]);
alg->setProperty("gamma", params[5]);
alg->executeAsChildAlg();
ws = alg->getProperty("PeaksWorkspace");
OrientedLattice latt = ws->mutableSample().getOrientedLattice();
DblMatrix UB = latt.getUB();
DblMatrix A = aMatrix(params);
DblMatrix Bc = A;
Bc.Invert();
DblMatrix U1_B1 = UB * A;
OrientedLattice o_lattice;
o_lattice.setUB(U1_B1);
DblMatrix U1 = o_lattice.getU();
DblMatrix U1_Bc = U1 * Bc;
for (size_t i = 0; i < hkl_vector.size(); i++) {
V3D error = U1_Bc * hkl_vector[i] - q_vector[i] / (2.0 * M_PI);
out[i] = error.norm2();
}
return;
}
示例5: checkNumberPeaks
/** Count the peaks from a .peaks file and compare with the workspace
* @param outWS :: the workspace in which to place the information
* @param filename :: path to the .peaks file
*/
void LoadIsawPeaks::checkNumberPeaks( PeaksWorkspace_sptr outWS, std::string filename )
{
// Open the file
std::ifstream in( filename.c_str() );
std::string first;
int NumberPeaks = 0;
while (getline(in,first))
{
if (first[0] == '3')NumberPeaks++;
}
if(NumberPeaks != outWS->getNumberPeaks())
{
g_log.error()<<"Number of peaks in file is " << NumberPeaks << " but only read "
<<outWS->getNumberPeaks() << std::endl;
throw std::length_error("Wrong number of peaks read");
}
}
示例6: fillPossibleHKLsUsingPeaksWorkspace
/// Fills possibleHKLs with all HKLs from the supplied PeaksWorkspace.
void PredictPeaks::fillPossibleHKLsUsingPeaksWorkspace(
const PeaksWorkspace_sptr &peaksWorkspace,
std::vector<V3D> &possibleHKLs) const {
possibleHKLs.clear();
possibleHKLs.reserve(peaksWorkspace->getNumberPeaks());
bool roundHKL = getProperty("RoundHKL");
for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
IPeak &p = peaksWorkspace->getPeak(i);
// Get HKL from that peak
V3D hkl = p.getHKL();
if (roundHKL)
hkl.round();
possibleHKLs.push_back(hkl);
} // for each hkl in the workspace
}
示例7: CheckForOneRun
/**
* Checks that a PeaksWorkspace has only one run.
*
* @param Peaks The PeaksWorkspace
* @param GoniometerMatrix the goniometer matrix for the run
*/
bool GoniometerAnglesFromPhiRotation::CheckForOneRun(
const PeaksWorkspace_sptr &Peaks,
Kernel::Matrix<double> &GoniometerMatrix) const {
int RunNumber = -1;
for (int peak = 0; peak < Peaks->getNumberPeaks(); peak++) {
int thisRunNum = Peaks->getPeak(peak).getRunNumber();
GoniometerMatrix = Peaks->getPeak(peak).getGoniometerMatrix();
if (RunNumber < 0)
RunNumber = thisRunNum;
else if (thisRunNum != RunNumber)
return false;
}
return true;
}
示例8: invalid_argument
/**
* Creates a new parameterized instrument for which the parameter values can be
*changed
*
* @param Peaks - a PeaksWorkspace used to get the original instrument. The
*instrument from the 0th peak is
* the one that is used.
*
* NOTE: All the peaks in the PeaksWorkspace must use the same instrument.
*/
boost::shared_ptr<Geometry::Instrument>
PeakHKLErrors::getNewInstrument(PeaksWorkspace_sptr Peaks) const {
Geometry::Instrument_const_sptr instSave = Peaks->getPeak(0).getInstrument();
auto pmap = boost::make_shared<Geometry::ParameterMap>();
if (!instSave) {
g_log.error(" Peaks workspace does not have an instrument");
throw std::invalid_argument(" Not all peaks have an instrument");
}
if (!hasParameterMap) {
pmapSv = instSave->getParameterMap();
hasParameterMap = true;
if (!instSave->isParametrized()) {
boost::shared_ptr<Geometry::Instrument> instClone(instSave->clone());
auto Pinsta = boost::make_shared<Geometry::Instrument>(instSave, pmap);
instChange = Pinsta;
IComponent_const_sptr sample = instChange->getSample();
sampPos = sample->getRelativePos();
} else // catch(... )
{
auto P1 = boost::make_shared<Geometry::Instrument>(
instSave->baseInstrument(), instSave->makeLegacyParameterMap());
instChange = P1;
IComponent_const_sptr sample = instChange->getSample();
sampPos = sample->getRelativePos();
}
}
if (!instChange) {
g_log.error("Cannot 'clone' instrument");
throw std::logic_error("Cannot clone instrument");
}
//------------------"clone" orig instruments pmap -------------------
cLone(pmap, instSave, pmapSv);
V3D sampOffsets(getParameter("SampleXOffset"), getParameter("SampleYOffset"),
getParameter("SampleZOffset"));
IComponent_const_sptr sample = instChange->getSample();
pmap->addPositionCoordinate(sample.get(), std::string("x"),
sampPos.X() + sampOffsets.X());
pmap->addPositionCoordinate(sample.get(), std::string("y"),
sampPos.Y() + sampOffsets.Y());
pmap->addPositionCoordinate(sample.get(), std::string("z"),
sampPos.Z() + sampOffsets.Z());
return instChange;
}
示例9: edgePixel
/**
@param ws Name of workspace containing peaks
@param bankName Name of bank containing peak
@param col Column number containing peak
@param row Row number containing peak
@param Edge Number of edge points for each bank
@return True if peak is on edge
*/
bool OptimizeLatticeForCellType::edgePixel(PeaksWorkspace_sptr ws,
std::string bankName, int col,
int row, int Edge) {
if (bankName.compare("None") == 0)
return false;
Geometry::Instrument_const_sptr Iptr = ws->getInstrument();
boost::shared_ptr<const IComponent> parent =
Iptr->getComponentByName(bankName);
if (parent->type().compare("RectangularDetector") == 0) {
boost::shared_ptr<const RectangularDetector> RDet =
boost::dynamic_pointer_cast<const RectangularDetector>(parent);
return col < Edge || col >= (RDet->xpixels() - Edge) || row < Edge ||
row >= (RDet->ypixels() - Edge);
} else {
std::vector<Geometry::IComponent_const_sptr> children;
boost::shared_ptr<const Geometry::ICompAssembly> asmb =
boost::dynamic_pointer_cast<const Geometry::ICompAssembly>(parent);
asmb->getChildren(children, false);
int startI = 1;
if (children[0]->getName() == "sixteenpack") {
startI = 0;
parent = children[0];
children.clear();
boost::shared_ptr<const Geometry::ICompAssembly> asmb =
boost::dynamic_pointer_cast<const Geometry::ICompAssembly>(parent);
asmb->getChildren(children, false);
}
boost::shared_ptr<const Geometry::ICompAssembly> asmb2 =
boost::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[0]);
std::vector<Geometry::IComponent_const_sptr> grandchildren;
asmb2->getChildren(grandchildren, false);
int NROWS = static_cast<int>(grandchildren.size());
int NCOLS = static_cast<int>(children.size());
// Wish pixels and tubes start at 1 not 0
return col - startI < Edge || col - startI >= (NCOLS - Edge) ||
row - startI < Edge || row - startI >= (NROWS - Edge);
}
return false;
}
示例10: UBinv
void PeakHKLErrors::functionDeriv1D(Jacobian *out, const double *xValues,
const size_t nData) {
PeaksWorkspace_sptr Peaks =
AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(
PeakWorkspaceName);
boost::shared_ptr<Geometry::Instrument> instNew = getNewInstrument(Peaks);
const DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
DblMatrix UBinv(UB);
UBinv.Invert();
UBinv /= 2 * M_PI;
double GonRotx = getParameter("GonRotx");
double GonRoty = getParameter("GonRoty");
double GonRotz = getParameter("GonRotz");
Matrix<double> InvGonRotxMat = RotationMatrixAboutRegAxis(GonRotx, 'x');
Matrix<double> InvGonRotyMat = RotationMatrixAboutRegAxis(GonRoty, 'y');
Matrix<double> InvGonRotzMat = RotationMatrixAboutRegAxis(GonRotz, 'z');
Matrix<double> GonRot = InvGonRotxMat * InvGonRotyMat * InvGonRotzMat;
InvGonRotxMat.Invert();
InvGonRotyMat.Invert();
InvGonRotzMat.Invert();
std::map<int, Kernel::Matrix<double>> RunNums2GonMatrix;
getRun2MatMap(Peaks, OptRuns, RunNums2GonMatrix);
g_log.debug()
<< "----------------------------Derivative------------------------\n";
V3D samplePosition = instNew->getSample()->getPos();
IPeak &ppeak = Peaks->getPeak(0);
double L0 = ppeak.getL1();
double velocity = (L0 + ppeak.getL2()) / ppeak.getTOF();
double K =
2 * M_PI / ppeak.getWavelength() / velocity; // 2pi/lambda = K* velocity
V3D beamDir = instNew->getBeamDirection();
size_t paramNums[] = {parameterIndex(std::string("SampleXOffset")),
parameterIndex(std::string("SampleYOffset")),
parameterIndex(std::string("SampleZOffset"))};
for (size_t i = 0; i < nData; i += 3) {
int peakNum = boost::math::iround(xValues[i]);
IPeak &peak_old = Peaks->getPeak(peakNum);
Peak peak = createNewPeak(peak_old, instNew, 0, peak_old.getL1());
int runNum = peak_old.getRunNumber();
std::string runNumStr = std::to_string(runNum);
for (int kk = 0; kk < static_cast<int>(nParams()); kk++) {
out->set(i, kk, 0.0);
out->set(i + 1, kk, 0.0);
out->set(i + 2, kk, 0.0);
}
double chi, phi, omega;
size_t chiParamNum, phiParamNum, omegaParamNum;
size_t N = OptRuns.find("/" + runNumStr);
if (N < OptRuns.size()) {
chi = getParameter("chi" + (runNumStr));
phi = getParameter("phi" + (runNumStr));
omega = getParameter("omega" + (runNumStr));
peak.setGoniometerMatrix(GonRot * RunNums2GonMatrix[runNum]);
chiParamNum = parameterIndex("chi" + (runNumStr));
phiParamNum = parameterIndex("phi" + (runNumStr));
omegaParamNum = parameterIndex("omega" + (runNumStr));
} else {
Geometry::Goniometer Gon(peak.getGoniometerMatrix());
std::vector<double> phichiOmega = Gon.getEulerAngles("YZY");
chi = phichiOmega[1];
phi = phichiOmega[2];
omega = phichiOmega[0];
// peak.setGoniometerMatrix( GonRot*Gon.getR());
chiParamNum = phiParamNum = omegaParamNum = nParams() + 10;
peak.setGoniometerMatrix(GonRot * peak.getGoniometerMatrix());
}
V3D sampOffsets(getParameter("SampleXOffset"),
getParameter("SampleYOffset"),
getParameter("SampleZOffset"));
peak.setSamplePos(peak.getSamplePos() + sampOffsets);
// NOTE:Use getQLabFrame except for below.
// For parameters the getGoniometerMatrix should remove GonRot, for derivs
// wrt GonRot*, wrt chi*,phi*,etc.
// Deriv wrt chi phi and omega
if (phiParamNum < nParams()) {
Matrix<double> chiMatrix = RotationMatrixAboutRegAxis(chi, 'z');
Matrix<double> phiMatrix = RotationMatrixAboutRegAxis(phi, 'y');
Matrix<double> omegaMatrix = RotationMatrixAboutRegAxis(omega, 'y');
Matrix<double> dchiMatrix = DerivRotationMatrixAboutRegAxis(chi, 'z');
Matrix<double> dphiMatrix = DerivRotationMatrixAboutRegAxis(phi, 'y');
Matrix<double> domegaMatrix = DerivRotationMatrixAboutRegAxis(omega, 'y');
//.........这里部分代码省略.........
示例11: o
/** Execute the algorithm.
*/
void CalculateUMatrix::exec()
{
double a=this->getProperty("a");
double b=this->getProperty("b");
double c=this->getProperty("c");
double alpha=this->getProperty("alpha");
double beta=this->getProperty("beta");
double gamma=this->getProperty("gamma");
OrientedLattice o(a,b,c,alpha,beta,gamma);
Matrix<double> B=o.getB();
double H,K,L;
PeaksWorkspace_sptr ws;
ws = AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(this->getProperty("PeaksWorkspace") );
if (!ws) throw std::runtime_error("Problems reading the peaks workspace");
size_t nIndexedpeaks=0;
bool found2nc=false;
V3D old(0,0,0);
Matrix<double> Hi(4,4),Si(4,4),HS(4,4),zero(4,4);
for (int i=0;i<ws->getNumberPeaks();i++)
{
Peak p=ws->getPeaks()[i];
H=p.getH();
K=p.getK();
L=p.getL();
if(H*H+K*K+L*L>0)
{
nIndexedpeaks++;
if (!found2nc)
{
if (nIndexedpeaks==1)
{
old=V3D(H,K,L);
}
else
{
if (!old.coLinear(V3D(0,0,0),V3D(H,K,L))) found2nc=true;
}
}
V3D Qhkl=B*V3D(H,K,L);
Hi[0][0]=0.; Hi[0][1]=-Qhkl.X(); Hi[0][2]=-Qhkl.Y(); Hi[0][3]=-Qhkl.Z();
Hi[1][0]=Qhkl.X(); Hi[1][1]=0.; Hi[1][2]=Qhkl.Z(); Hi[1][3]=-Qhkl.Y();
Hi[2][0]=Qhkl.Y(); Hi[2][1]=-Qhkl.Z(); Hi[2][2]=0.; Hi[2][3]=Qhkl.X();
Hi[3][0]=Qhkl.Z(); Hi[3][1]=Qhkl.Y(); Hi[3][2]=-Qhkl.X(); Hi[3][3]=0.;
V3D Qgon=p.getQSampleFrame();
Si[0][0]=0.; Si[0][1]=-Qgon.X(); Si[0][2]=-Qgon.Y(); Si[0][3]=-Qgon.Z();
Si[1][0]=Qgon.X(); Si[1][1]=0.; Si[1][2]=-Qgon.Z(); Si[1][3]=Qgon.Y();
Si[2][0]=Qgon.Y(); Si[2][1]=Qgon.Z(); Si[2][2]=0.; Si[2][3]=-Qgon.X();
Si[3][0]=Qgon.Z(); Si[3][1]=-Qgon.Y(); Si[3][2]=Qgon.X(); Si[3][3]=0.;
HS+=(Hi*Si);
}
}
//check if enough peaks are indexed or if HS is 0
if ((nIndexedpeaks<2) || (found2nc==false)) throw std::invalid_argument("Less then two non-colinear peaks indexed");
if (HS==zero) throw std::invalid_argument("Something really bad happened");
Matrix<double> Eval;
Matrix<double> Diag;
HS.Diagonalise(Eval,Diag);
Eval.sortEigen(Diag);
Mantid::Kernel::Quat qR(Eval[0][0],Eval[1][0],Eval[2][0],Eval[3][0]);//the first column corresponds to the highest eigenvalue
DblMatrix U(qR.getRotation());
o.setU(U);
ws->mutableSample().setOrientedLattice(new OrientedLattice(o));
}
示例12: UNUSED_ARG
int PeakIntegration::fitneighbours(int ipeak, std::string det_name, int x0,
int y0, int idet, double qspan,
PeaksWorkspace_sptr &Peaks,
const detid2index_map &pixel_to_wi) {
UNUSED_ARG(ipeak);
UNUSED_ARG(det_name);
UNUSED_ARG(x0);
UNUSED_ARG(y0);
Geometry::IPeak &peak = Peaks->getPeak(ipeak);
// Number of slices
int TOFmax = 0;
IAlgorithm_sptr slice_alg = createChildAlgorithm("IntegratePeakTimeSlices");
slice_alg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputW);
std::ostringstream tab_str;
tab_str << "LogTable" << ipeak;
slice_alg->setPropertyValue("OutputWorkspace", tab_str.str());
slice_alg->setProperty<PeaksWorkspace_sptr>("Peaks", Peaks);
slice_alg->setProperty("PeakIndex", ipeak);
slice_alg->setProperty("PeakQspan", qspan);
int nPixels = std::max<int>(0, getProperty("NBadEdgePixels"));
slice_alg->setProperty("NBadEdgePixels", nPixels);
slice_alg->executeAsChildAlg();
Mantid::API::MemoryManager::Instance().releaseFreeMemory();
MantidVec &Xout = outputW->dataX(idet);
MantidVec &Yout = outputW->dataY(idet);
MantidVec &Eout = outputW->dataE(idet);
TableWorkspace_sptr logtable = slice_alg->getProperty("OutputWorkspace");
peak.setIntensity(slice_alg->getProperty("Intensity"));
peak.setSigmaIntensity(slice_alg->getProperty("SigmaIntensity"));
TOFmax = static_cast<int>(logtable->rowCount());
for (int iTOF = 0; iTOF < TOFmax; iTOF++) {
Xout[iTOF] = logtable->getRef<double>(std::string("Time"), iTOF);
if (m_IC) // Ikeda-Carpenter fit
{
Yout[iTOF] = logtable->getRef<double>(std::string("TotIntensity"), iTOF);
Eout[iTOF] =
logtable->getRef<double>(std::string("TotIntensityError"), iTOF);
} else {
Yout[iTOF] = logtable->getRef<double>(std::string("ISAWIntensity"), iTOF);
Eout[iTOF] =
logtable->getRef<double>(std::string("ISAWIntensityError"), iTOF);
}
}
outputW->getSpectrum(idet)->clearDetectorIDs();
// Find the pixel ID at that XY position on the rectangular detector
int pixelID = peak.getDetectorID(); // det->getAtXY(x0,y0)->getID();
// Find the corresponding workspace index, if any
auto wiEntry = pixel_to_wi.find(pixelID);
if (wiEntry != pixel_to_wi.end()) {
size_t wi = wiEntry->second;
// Set detectorIDs
outputW->getSpectrum(idet)
->addDetectorIDs(inputW->getSpectrum(wi)->getDetectorIDs());
}
return TOFmax - 1;
}
示例13: if
/**
@param inname Name of Filename containing peaks
@param cell_type cell type to optimize
@param params optimized cell parameters
@return chisq of optimization
*/
double OptimizeLatticeForCellType::optLatticeSum(std::string inname,
std::string cell_type,
std::vector<double> ¶ms) {
std::vector<double> lattice_parameters;
lattice_parameters.assign(6, 0);
if (cell_type == ReducedCell::CUBIC()) {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[0];
lattice_parameters[2] = params[0];
lattice_parameters[3] = 90;
lattice_parameters[4] = 90;
lattice_parameters[5] = 90;
} else if (cell_type == ReducedCell::TETRAGONAL()) {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[0];
lattice_parameters[2] = params[1];
lattice_parameters[3] = 90;
lattice_parameters[4] = 90;
lattice_parameters[5] = 90;
} else if (cell_type == ReducedCell::ORTHORHOMBIC()) {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[1];
lattice_parameters[2] = params[2];
lattice_parameters[3] = 90;
lattice_parameters[4] = 90;
lattice_parameters[5] = 90;
} else if (cell_type == ReducedCell::RHOMBOHEDRAL()) {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[0];
lattice_parameters[2] = params[0];
lattice_parameters[3] = params[1];
lattice_parameters[4] = params[1];
lattice_parameters[5] = params[1];
} else if (cell_type == ReducedCell::HEXAGONAL()) {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[0];
lattice_parameters[2] = params[1];
lattice_parameters[3] = 90;
lattice_parameters[4] = 90;
lattice_parameters[5] = 120;
} else if (cell_type == "Monoclinic ( a unique )") {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[1];
lattice_parameters[2] = params[2];
lattice_parameters[3] = params[3];
lattice_parameters[4] = 90;
lattice_parameters[5] = 90;
} else if (cell_type == ReducedCell::MONOCLINIC() ||
cell_type == "Monoclinic ( b unique )") {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[1];
lattice_parameters[2] = params[2];
lattice_parameters[3] = 90;
lattice_parameters[4] = params[3];
lattice_parameters[5] = 90;
} else if (cell_type == "Monoclinic ( c unique )") {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[1];
lattice_parameters[2] = params[2];
lattice_parameters[3] = 90;
lattice_parameters[4] = 90;
lattice_parameters[5] = params[3];
} else if (cell_type == ReducedCell::TRICLINIC()) {
lattice_parameters[0] = params[0];
lattice_parameters[1] = params[1];
lattice_parameters[2] = params[2];
lattice_parameters[3] = params[3];
lattice_parameters[4] = params[4];
lattice_parameters[5] = params[5];
}
PeaksWorkspace_sptr ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
AnalysisDataService::Instance().retrieve(inname));
size_t n_peaks = ws->getNumberPeaks();
double *out = new double[n_peaks];
optLattice(inname, lattice_parameters, out);
double ChiSqTot = 0;
for (size_t i = 0; i < n_peaks; i++)
ChiSqTot += out[i];
delete[] out;
return ChiSqTot;
}
示例14: appendFile
/** Append the peaks from a .peaks file into the workspace
* @param outWS :: the workspace in which to place the information
* @param filename :: path to the .peaks file
*/
void LoadIsawPeaks::appendFile( PeaksWorkspace_sptr outWS, std::string filename )
{
// Open the file
std::ifstream in( filename.c_str() );
// Read the header, load the instrument
double T0;
std::string s = readHeader( outWS, in , T0);
// set T0 in the run parameters
API::Run & m_run = outWS->mutableRun();
m_run.addProperty<double>("T0", T0, true);
if( !in.good() || s.length() < 1 )
throw std::runtime_error( "End of Peaks file before peaks" );
if( s.compare( std::string( "0" ) ) != 0 )
throw std::logic_error( "No header for Peak segments" );
readToEndOfLine( in , true );
s = getWord( in , false );
int run, bankNum;
double chi , phi , omega , monCount;
// Build the universal goniometer that will build the rotation matrix.
Mantid::Geometry::Goniometer uniGonio;
uniGonio.makeUniversalGoniometer();
// TODO: Can we find the number of peaks to get better progress reporting?
Progress prog(this, 0.0, 1.0, 100);
while( in.good() )
{
// Read the header if necessary
s = readPeakBlockHeader( s , in , run , bankNum , chi , phi ,
omega , monCount );
// Build the Rotation matrix using phi,chi,omega
uniGonio.setRotationAngle("phi", phi);
uniGonio.setRotationAngle("chi", chi);
uniGonio.setRotationAngle("omega", omega);
//Put goniometer into peaks workspace
outWS->mutableRun().setGoniometer(uniGonio, false);
std::ostringstream oss;
std::string bankString = "bank";
if (outWS->getInstrument()->getName() == "WISH") bankString = "WISHpanel0";
oss << bankString << bankNum;
std::string bankName = oss.str();
int seqNum = -1;
try
{
// Read the peak
Peak peak = readPeak(outWS, s, in, seqNum, bankName);
// Get the calculated goniometer matrix
Matrix<double> gonMat = uniGonio.getR();
peak.setGoniometerMatrix(gonMat);
peak.setRunNumber(run);
peak.setMonitorCount( monCount );
double tof = peak.getTOF();
Kernel::Units::Wavelength wl;
wl.initialize(peak.getL1(), peak.getL2(), peak.getScattering(), 0,
peak.getInitialEnergy(), 0.0);
peak.setWavelength(wl.singleFromTOF( tof));
// Add the peak to workspace
outWS->addPeak(peak);
}
catch (std::runtime_error & e)
{
g_log.warning() << "Error reading peak SEQN " << seqNum << " : " << e.what() << std::endl;
}
prog.report();
}
}
示例15: readPeak
/** Read one peak in a line of an ISAW peaks file.
*
* @param outWS :: workspace to add peaks to
* @param lastStr [in,out] :: last word (the one at the start of the line)
* @param in :: input stream
* @param seqNum [out] :: the sequence number of the peak
* @param bankName :: the bank number from the ISAW file.
* @return the Peak the Peak object created
*/
Mantid::DataObjects::Peak readPeak( PeaksWorkspace_sptr outWS, std::string & lastStr, std::ifstream& in, int & seqNum, std::string bankName)
{
double h ; double k ; double l ; double col ;
double row ; double wl ;
double IPK ; double Inti ;
double SigI ;
seqNum = -1;
std::string s = lastStr;
if( s.length() < 1 && in.good() )//blank line
{
readToEndOfLine( in , true );
s = getWord( in , false );;
}
if( s.length() < 1 )
throw std::runtime_error("Empty peak line encountered.");
if( s.compare( "2" ) == 0 )
{
readToEndOfLine( in , true );
for( s = getWord( in , false ) ; s.length() < 1 && in.good() ;
s = getWord( in , true ) )
{
s = getWord( in , false );
}
}
if( s.length() < 1 )
throw std::runtime_error("Empty peak line encountered.");
if( s.compare( "3" ) != 0 )
throw std::runtime_error("Empty peak line encountered.");
seqNum = atoi( getWord( in , false ).c_str() );
h = strtod( getWord( in , false ).c_str() , 0 ) ;
k = strtod( getWord( in , false ).c_str() , 0 ) ;
l = strtod( getWord( in , false ).c_str() , 0 ) ;
col = strtod( getWord( in , false ).c_str() , 0 ) ;
row = strtod( getWord( in , false ).c_str() , 0 ) ;
strtod( getWord( in , false ).c_str() , 0 ) ; //chan
strtod( getWord( in , false ).c_str() , 0 ) ; //L2
strtod( getWord( in , false ).c_str() , 0 ) ; //ScatAng
strtod( getWord( in , false ).c_str() , 0 ) ; //Az
wl = strtod( getWord( in , false ).c_str() , 0 ) ;
strtod( getWord( in , false ).c_str() , 0 ) ; //D
IPK = strtod( getWord( in , false ).c_str() , 0 ) ;
Inti = strtod( getWord( in , false ).c_str() , 0 ) ;
SigI = strtod( getWord( in , false ).c_str() , 0 ) ;
atoi( getWord( in , false ).c_str() ) ; // iReflag
// Finish the line and get the first word of next line
readToEndOfLine( in , true );
lastStr = getWord( in , false );
// Find the detector ID from row/col
Instrument_const_sptr inst = outWS->getInstrument();
if (!inst) throw std::runtime_error("No instrument in PeaksWorkspace!");
LoadIsawPeaks u;
int pixelID = u.findPixelID(inst, bankName, static_cast<int>(col), static_cast<int>(row));
//Create the peak object
Peak peak(outWS->getInstrument(), pixelID, wl);
// HKL's are flipped by -1 because of the internal Q convention
peak.setHKL(-h,-k,-l);
peak.setIntensity(Inti);
peak.setSigmaIntensity(SigI);
peak.setBinCount(IPK);
// Return the peak
return peak;
}