Commit 23894fee authored by Ruth Mikkelson's avatar Ruth Mikkelson
Browse files

Refs #5487 LoadIsawPeaks now incorporates Calibration info

Values in test programs were slighly off.  Also, changes were needed
when the instrument was parameterized from the start.
parent c1ce8b01
......@@ -37,7 +37,13 @@ namespace Crystal
/// Run the algorithm
void exec();
std::string readHeader( Mantid::DataObjects::PeaksWorkspace_sptr outWS, std::ifstream& in );
std::string ApplyCalibInfo(std::ifstream & in,
std::string startChar,
Geometry::Instrument_const_sptr instr_old,
Geometry::Instrument_const_sptr instr,
double &T0);
std::string readHeader( Mantid::DataObjects::PeaksWorkspace_sptr outWS, std::ifstream& in,double &T0 );
void appendFile( Mantid::DataObjects::PeaksWorkspace_sptr outWS, std::string filename);
......
......@@ -7,7 +7,6 @@ Reads an ISAW-style .peaks or .integrate file into a PeaksWorkspace. Any detecto
NOTE: The instrument used is determined by reading the 'Instrument:' and 'Date:' tags at the start of the file.If the date is not present, the latest [[Instrument Definition File]] is used.
*WIKI*/
#include "MantidAPI/FileProperty.h"
#include "MantidCrystal/LoadIsawPeaks.h"
......@@ -18,6 +17,7 @@ NOTE: The instrument used is determined by reading the 'Instrument:' and 'Date:'
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/System.h"
#include "MantidCrystal/SCDCalibratePanels.h"
#include <algorithm>
#include <boost/shared_ptr.hpp>
#include <exception>
......@@ -28,9 +28,11 @@ NOTE: The instrument used is determined by reading the 'Instrument:' and 'Date:'
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include "MantidKernel/Unit.h"
using Mantid::Kernel::Strings::readToEndOfLine;
using Mantid::Kernel::Strings::getWord;
using Mantid::Kernel::Units::Wavelength;
namespace Mantid
{
......@@ -84,6 +86,137 @@ namespace Crystal
}
std::string LoadIsawPeaks::ApplyCalibInfo(std::ifstream & in, std::string startChar,Geometry::Instrument_const_sptr instr_old, Geometry::Instrument_const_sptr instr,
double &T0)
{
ParameterMap_sptr parMap1= instr_old->getParameterMap();
ParameterMap_sptr parMap= instr->getParameterMap();
std::string lastString ="";
while( in.good() && (startChar.size() <1 || startChar !="7") )
{
readToEndOfLine( in, true);
startChar = getWord(in, false);
}
if( !(in.good()))
{
//g_log.error()<<"Peaks file has no time shift and L0 info"<<std::endl;
throw std::invalid_argument("Peaks file has no time shift and L0 info");
}
std::string L1s= getWord(in,false);
std::string T0s =getWord(in, false);
if( L1s.length() < 1 || T0s.length() < 1)
{
g_log.error()<<"Missing L1 or Time offset"<<std::endl;
throw std::invalid_argument("Missing L1 or Time offset");
}
double L1;
try
{
std::istringstream iss( L1s+" "+T0s, std::istringstream::in);
iss>>L1;
iss>>T0;
SCDCalibratePanels::FixUpSourceParameterMap(instr, L1/100, parMap1);
}catch(...)
{
g_log.error()<<"Invalid L1 or Time offset"<<std::endl;
throw std::invalid_argument("Invalid L1 or Time offset");
}
readToEndOfLine( in, true);
startChar = getWord(in , false);
while( in.good() && (startChar.size() <1 || startChar !="5") )
{
readToEndOfLine( in, true);
startChar = getWord(in, false);
}
if( !(in.good()))
{
g_log.error()<<"Peaks file has no detector panel info"<<std::endl;
throw std::invalid_argument("Peaks file has no detector panel info");
}
while( startChar =="5")
{
std::string line;
for( int i=0; i<16;i++)
{
std::string s= getWord(in, false);
if( s.size() < 1)
{
g_log.error()<<"Not enough info to describe panel "<<std::endl;
throw std::length_error("Not enough info to describe panel ");
}
line +=" "+s;;
}
readToEndOfLine(in, true);
startChar = getWord( in, false);// blank lines ?? and # lines ignore
std::istringstream iss( line, std::istringstream::in);
int bankNum,nrows,ncols;
double width,height,depth,detD,Centx,Centy,Centz,Basex,Basey,Basez,
Upx,Upy,Upz;
try
{
iss>>bankNum>>nrows>>ncols>>width>>height>>depth>>detD
>>Centx>>Centy>>Centz>>Basex>>Basey>>Basez
>>Upx>>Upy>>Upz;
}catch(...)
{
g_log.error()<<"incorrect type of data for panel "<<std::endl;
throw std::length_error("incorrect type of data for panel ");
}
std::string SbankNum = boost::lexical_cast<std::string>(bankNum);
std::string bankName = "bank"+SbankNum;
boost::shared_ptr<const Geometry::IComponent> bank =instr_old->getComponentByName( bankName );
if( !bank)
{
g_log.error()<<"There is no bank "<< bankName<<" in the instrument"<<std::endl;
throw std::length_error("There is no bank "+ bankName+" in the instrument");
}
V3D dPos= V3D(Centx,Centy,Centz)/100.0- bank->getPos();
V3D Base(Basex,Basey,Basez), Up(Upx,Upy,Upz);
V3D ToSamp =Base.cross_prod(Up);
Base.normalize();
Up.normalize();
ToSamp.normalize();
Quat thisRot(Base,Up,ToSamp);
Quat bankRot(bank->getRotation());
bankRot.inverse();
Quat dRot = thisRot*bankRot;
boost::shared_ptr< const Geometry::RectangularDetector>bankR= boost::dynamic_pointer_cast
<const Geometry::RectangularDetector>( bank);
double DetWScale = 1, DetHtScale = 1;
if( bank)
{
DetWScale = width/bankR->xsize()/100;
DetHtScale = height/bankR->ysize()/100;
}
std::vector<std::string> bankNames;
bankNames.push_back(bankName);
SCDCalibratePanels::FixUpBankParameterMap(bankNames,instr, dPos,
dRot,DetWScale,DetHtScale , parMap1);
}
return startChar;
}
//-----------------------------------------------------------------------------------------------
/** Reads the header of a .peaks file
......@@ -91,7 +224,7 @@ namespace Crystal
* @param in :: stream of the input file
* @return the first word on the next line
*/
std::string LoadIsawPeaks::readHeader( PeaksWorkspace_sptr outWS, std::ifstream& in )
std::string LoadIsawPeaks::readHeader( PeaksWorkspace_sptr outWS, std::ifstream& in, double &T0 )
{
std::string tag;
std::string r = getWord( in , false );
......@@ -139,18 +272,25 @@ namespace Crystal
// Populate the instrument parameters in this workspace - this works around a bug
tempWS->populateInstrumentParameters();
outWS->setInstrument( tempWS->getInstrument() );
Geometry::Instrument_const_sptr instr_old = tempWS->getInstrument() ;
boost::shared_ptr< ParameterMap > map(new ParameterMap());
Geometry::Instrument_const_sptr instr ( new Geometry::Instrument(instr_old->baseInstrument(), map ));
//std::string s;
std::string s = ApplyCalibInfo(in, "", instr_old, instr, T0);
outWS->setInstrument( instr);
// Now skip all lines on L1, detector banks, etc. until we get to a block of peaks. They start with 0.
readToEndOfLine( in , true );
readToEndOfLine( in , true );
std::string s = getWord(in, false);
while (s != "0")
// readToEndOfLine( in , true );
// readToEndOfLine( in , true );
// s = getWord(in, false);
while (s != "0" && in.good())
{
readToEndOfLine( in , true );
s = getWord(in, false);
}
return s;
}
......@@ -243,7 +383,6 @@ namespace Crystal
peak.setIntensity(Inti);
peak.setSigmaIntensity(SigI);
peak.setBinCount(IPK);
// Return the peak
return peak;
}
......@@ -303,14 +442,16 @@ namespace Crystal
* @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)
void LoadIsawPeaks::appendFile( PeaksWorkspace_sptr outWS, std::string filename )
{
// Open the file
std::ifstream in( filename.c_str() );
// Read the header, load the instrument
std::string s = readHeader( outWS, in );
double T0;
std::string s = readHeader( outWS, in , T0);
if( !in.good() || s.length() < 1 )
throw std::runtime_error( "End of Peaks file before peaks" );
......@@ -359,6 +500,13 @@ namespace Crystal
peak.setGoniometerMatrix(gonMat);
peak.setRunNumber(run);
double tof = peak.getTOF()+T0;
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);
}
......
......@@ -31,7 +31,7 @@ Some features:
*/
//TODO
// Rotate should also rotate around centroid of Group
// 1. Change xxx --> SCDCalibratePanelsx where x=0,1,2,3,...
......@@ -393,6 +393,7 @@ namespace Crystal
boost::shared_ptr<const Instrument> newInstrument = ws->getInstrument();
newInstrument->getInstrumentParameters(L0,beamline,norm,samplePos);
timeOffset = LoadDetCal->getProperty("TimeOffset");
AnalysisDataService::Instance().remove("fff");
return newInstrument;
}
......@@ -484,7 +485,10 @@ namespace Crystal
Quat ChgRot = rotPre*rotI;
Quat2RotxRotyRotz(ChgRot,Xrot0,Yrot0,Zrot0);
std::cout<<"new Initial Values="<<std::endl;
std::cout<<" "<< detWidthScale0<<","<<detHeightScale0
<<","<<Xoffset0<<","<<Yoffset0<<","<<Zoffset0
<<","<<Xrot0<<","<<Yrot0<<","<<Zrot0<<std::endl;
}
......@@ -558,8 +562,13 @@ namespace Crystal
(string) getProperty("PreProcFilename"),
T0, L0, banksVec);
g_log.debug()<<"Initial L0,T0="<<L0<<","<<T0<<std::endl;
AnalysisDataService::Instance().addOrReplace("xxx",peaksWs );
string PeakWSName = getPropertyValue( "PeakWorkspace");
if(PeakWSName.length()<1)
{
PeakWSName = "xxx";
AnalysisDataService::Instance().addOrReplace("xxx",peaksWs );
}
string FunctionArgument;
string Constraints("");
......@@ -600,7 +609,7 @@ namespace Crystal
}
// if( i > 0 ) oss << ";";
oss << "name=SCDPanelErrors, PeakWorkspaceName=\"xxx\",";
oss << "name=SCDPanelErrors, PeakWorkspaceName=\""<<PeakWSName<<"\",";
oss << "a=" << fixed << a << "," << "b=" << fixed << b << "," << "c=" << fixed << c << "," << "alpha=" << fixed << alpha << "," << "beta=" << fixed << beta
<< "," << "gamma=" << fixed << gamma << ","<< "NGroups="<<NGroups<<",BankNames ="<<BankNameString<<","
<<"startX=-1,endX=-1,";
......@@ -920,6 +929,7 @@ namespace Crystal
string instName = instrument->getName();
i = -1;
std::cout<<"Ere set new values into instrument"<<std::endl;
for( vector<vector< string > >::iterator itv = Groups.begin(); itv != Groups.end(); ++itv )
{
i++;
......@@ -958,7 +968,7 @@ namespace Crystal
{
boost::shared_ptr<Algorithm>SaveDetCal = createSubAlgorithm( "SaveIsawDetCal");
//clone this
peaksWs->setInstrument( NewInstrument);
SaveDetCal->initialize();
SaveDetCal->setProperty( "InputWorkspace",peaksWs);
......@@ -971,6 +981,7 @@ namespace Crystal
g_log.notice()<<"Saved DetCal file in "<<DetCalFileName<< std::endl;
}
SaveXmlFile( XmlFileName, Groups,NewInstrument);
......@@ -1161,7 +1172,10 @@ namespace Crystal
fit->setAttribute("alpha",IFunction::Attribute((double)getProperty("alpha")));
fit->setAttribute("beta",IFunction::Attribute((double)getProperty("beta")));
fit->setAttribute("gamma",IFunction::Attribute((double)getProperty("gamma")));
fit->setAttribute("PeakWorkspaceName",IFunction::Attribute("xxx"));
string PeakWSName = getPropertyValue("PeakWorkspace");
if( PeakWSName.length()<1)
PeakWSName="xxx";
fit->setAttribute("PeakWorkspaceName",IFunction::Attribute(PeakWSName));
fit->setAttribute("startX",IFunction::Attribute(-1));
fit->setAttribute("endX",IFunction::Attribute(-1));
fit->setAttribute("NGroups",IFunction::Attribute(NGroups));
......@@ -1306,8 +1320,10 @@ namespace Crystal
const string bankName = (*it1);
boost::shared_ptr<const IComponent> bank = NewInstrument->getComponentByName(bankName);
boost::shared_ptr<const IComponent> bank1 = NewInstrument->getComponentByName(bankName);
boost::shared_ptr<const Geometry::RectangularDetector> bank = boost::dynamic_pointer_cast<const RectangularDetector>(bank1);//Component
updateBankParams(bank, pmap, pmapOld);
Quat RelRot = bank->getRelativeRot();
Quat newRelRot = rot * RelRot;
double rotx, roty, rotz;
......@@ -1330,12 +1346,12 @@ namespace Crystal
double scalex, scaley;
if (!oldScalex.empty())
scalex = oldScalex[0] + DetWScale;
scalex = oldScalex[0] * DetWScale;
else
scalex = DetWScale;
if (!oldScaley.empty())
scaley = oldScaley[0] + DetHtScale;
scaley = oldScaley[0] * DetHtScale;
else
scaley = DetHtScale;
......
......@@ -102,7 +102,6 @@ public:
TS_ASSERT_DELTA( p.getIntensity(), 221.83, 0.01);
TS_ASSERT_DELTA( p.getSigmaIntensity(), 15.02, 0.01);
TS_ASSERT_DELTA( p.getBinCount(), 8, 0.01);
TS_ASSERT_DELTA( p.getWavelength(), 0.761095, 0.001);
TS_ASSERT_DELTA( p.getL1(), 18.0, 1e-3);
TS_ASSERT_DELTA( p.getL2(), 0.461, 1e-3);
......
......@@ -61,10 +61,11 @@ public:
AnalysisDataService::Instance().remove("TOPAZ_3007");
alg= AlgorithmFactory::Instance().create("SCDCalibratePanels", 1);
std::cout<<"A"<<std::endl;
alg->initialize();
Peakws->setName("PeaksWsp");
alg->setProperty("PeakWorkspace", Peakws );
alg->setProperty("a",14.0);
alg->setProperty("b",19.3);
alg->setProperty("c", 8.6);
......@@ -77,10 +78,11 @@ public:
//alg->setProperty("DetCalFilename","abc.DetCal");
alg->setProperty("PanelGroups","SpecifyGroups");
alg->setProperty("Grouping","26");
alg->setPropertyValue("ResultWorkspace","Result");
alg->setPropertyValue("QErrorWorkspace","QErrorResult");
alg->execute();
alg->setPropertyValue("ResultWorkspace","Result");
boost::shared_ptr<TableWorkspace> Results = alg->getProperty("ResultWorkspace");
......@@ -91,7 +93,7 @@ public:
TS_ASSERT_DELTA(-0.416479,Results->cell<double>(9,1),.01);
TS_ASSERT_DELTA(-0.298888,Results->cell<double>(8,1),.01);
TS_ASSERT_DELTA(0.212273,Results->cell<double>(17,1),.01);
std::cout<<"G"<<std::endl;
/* for( int i=0; i<(int)Results->rowCount(); i++)
{
......
......@@ -121,18 +121,18 @@ public:
//calib.setParameter("Yrot",90);
calib.function1D(out.data(), xVals.data(), (size_t) N);
std::cout<<out[0]<<","<<out[4]<<","<<out[8]<<","<<out[10]<<std::endl;
double d = .0001;
TS_ASSERT_DELTA(out[0], -0.0038554, d);
double d = .00001;
TS_ASSERT_DELTA(out[0], -0.00396681, d);
TS_ASSERT_DELTA(out[4], 0.00756805, d);
TS_ASSERT_DELTA(out[4], 0.00734371, d);
TS_ASSERT_DELTA(out[8], 0.0267926, d);
TS_ASSERT_DELTA(out[8], 0.0268435, d);
TS_ASSERT_DELTA(out[10], 0.00858511, d);
TS_ASSERT_DELTA(out[10], 0.00880572, d);
//-------------------------Test the derivative --------------------------------
boost::shared_ptr<Jacob> Jac(new Jacob(10, N));
boost::shared_ptr<Jacob> Jac(new Jacob(10, N));
calib.functionDeriv1D(Jac.get(), xVals.data(), (size_t) N);
......
......@@ -79,7 +79,7 @@ public:
int num_indexed = alg.getProperty("NumIndexed");
TS_ASSERT_EQUALS( num_indexed, 43 );
double average_error = alg.getProperty("AverageError");
TS_ASSERT_DELTA( average_error, 0.0119856, 1e-5 );
TS_ASSERT_DELTA( average_error, 0.0119856, .0001 );
AnalysisDataService::Instance().remove(WSName);
}
......
......@@ -78,7 +78,7 @@ public:
int num_indexed = alg.getProperty("NumIndexed");
TS_ASSERT_EQUALS( num_indexed, 43 );
double average_error = alg.getProperty("AverageError");
TS_ASSERT_DELTA( average_error, 0.0119856, 1e-5 );
TS_ASSERT_DELTA( average_error, 0.0119856, .0001 );
AnalysisDataService::Instance().remove(WSName);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment