diff --git a/Code/Mantid/Framework/Algorithms/src/ModeratorTzero2.cpp b/Code/Mantid/Framework/Algorithms/src/ModeratorTzero2.cpp index fb44bd195b46..1259e9e398a5 100644 --- a/Code/Mantid/Framework/Algorithms/src/ModeratorTzero2.cpp +++ b/Code/Mantid/Framework/Algorithms/src/ModeratorTzero2.cpp @@ -84,12 +84,12 @@ void ModeratorTzero2::exec() //deltaE-mode (should be "indirect") std::vector Emode=m_instrument->getStringParameter("deltaE-mode"); - if(!Emode.size()) throw Exception::InstrumentDefinitionError("Unable to retrieve instrument geometry (direct or indirect) parameter", inputWS->getTitle()); + if(Emode.empty()) throw Exception::InstrumentDefinitionError("Unable to retrieve instrument geometry (direct or indirect) parameter", inputWS->getTitle()); if(Emode[0]!= "indirect") throw Exception::InstrumentDefinitionError("Instrument geometry must be of type indirect."); // extract formula from instrument parameters std::vector t0_formula=m_instrument->getStringParameter("t0_formula"); - if(!t0_formula.size()) throw Exception::InstrumentDefinitionError("Unable to retrieve t0_formula among instrument parameters"); + if(t0_formula.empty()) throw Exception::InstrumentDefinitionError("Unable to retrieve t0_formula among instrument parameters"); m_formula=t0_formula[0]; //Run execEvent if eventWorkSpace @@ -125,14 +125,19 @@ void ModeratorTzero2::exec() mu::Parser parser; parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser parser.SetExpr(m_formula); + E1=m_convfactor*(L1/m_t1min)*(L1/m_t1min); + double min_t0_next=parser.Eval(); // fast neutrons are shifted by min_t0_next, irrespective of tof MantidVec &inbins = inputWS->dataX(i); MantidVec &outbins = outputWS->dataX(i); // iterate over the time-of-flight values for(unsigned int ibin=0; ibin < inbins.size(); ibin++) { tof=inbins[ibin]; // current time-of-flight - if(tof<0.0) tof=0.0; // impossible substraction to negative time-of-flight values - outbins[ibin] = inbins[ibin]-CalculateT0(tof, L1, t2, E1, parser); + if(tof #include "MantidAlgorithms/ModeratorTzero2.h" +#include "MantidAlgorithms/Rebin.h" #include "MantidTestHelpers/WorkspaceCreationHelper.h" #include "MantidKernel/UnitFactory.h" #include "MantidDataObjects/EventList.h" @@ -13,7 +14,6 @@ #include "MantidGeometry/muParser_Silent.h" #include - using namespace Mantid::Kernel; using namespace Mantid::API; using namespace Mantid::Algorithms; @@ -36,7 +36,6 @@ void testExecThrowsDeltaEmode() { const std::string inputName("testWS"); Workspace2D_sptr testWS = createTestHistoWorkspace(); - AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service ModeratorTzero2 alg; initializeAlgorithm(alg,inputName,inputName); TS_ASSERT_THROWS(alg.execute(), Mantid::Kernel::Exception::InstrumentDefinitionError); @@ -48,10 +47,7 @@ void testExecThrowsTzeroFormula() { Workspace2D_sptr testWS = createTestHistoWorkspace(); const std::string inputName("testWS"); - AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service - Mantid::Geometry::Instrument_const_sptr instrument = testWS->getInstrument(); // pointer to the instrument - testWS->instrumentParameters().addString(instrument->getChild(0).get(),"deltaE-mode", "indirect"); - std::vector Emode=instrument->getStringParameter("deltaE-mode"); + testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),"deltaE-mode", "indirect"); ModeratorTzero2 alg; initializeAlgorithm(alg,inputName,inputName); TS_ASSERT_THROWS(alg.execute(), Mantid::Kernel::Exception::InstrumentDefinitionError); @@ -62,86 +58,59 @@ void testExecThrowsTzeroFormula() void testExecHistogramWorkspace() { Workspace2D_sptr testWS = createTestHistoWorkspace(); - const std::string inputName("testWS"); - AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service - const std::string formula("0.1*incidentEnergy^(-0.25)"); + const std::string formula("10*incidentEnergy^(-0.05)"); // formula for delay emission time from moderator versus incident energy addFormula(formula, boost::dynamic_pointer_cast(testWS)); // insert t0-formula as an instrument parameter ModeratorTzero2 alg; - initializeAlgorithm(alg,inputName,inputName); + initializeAlgorithm(alg,"testWS","testWS"); TS_ASSERT_THROWS_NOTHING(alg.execute()); - const int numHists(2); - const int numBins(2000); - double L1,L2,E2,v2,t2,tof,t1,v1,E1; - mu::Parser parser; - parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser - parser.SetExpr(formula); - Mantid::Geometry::IObjComponent_const_sptr sample = testWS->getInstrument()->getSample(); - L1=testWS->getInstrument()->getSource()->getDistance(*sample); - for(int iHist=0; iHistgetNumberHistograms(); iHist++) { Mantid::MantidVec &xdata=testWS->dataX(iHist); - Mantid::Geometry::IDetector_const_sptr det=testWS->getDetector(iHist); - L2=det->getDistance(*sample); - std::vector wsProp=det->getNumberParameter("Efixed"); - E2=wsProp.at(0); //[E2]=meV - v2=m_convFact*sqrt(E2); - t2=v2/L2; - for(int iBin=0; iBin<=numBins; ++iBin) - { - tof = 5.0 + 5.5*iBin; // original time of flight - t1=tof-t2; - v1=L1/t1; - E1=m_convFact*v1*v1; // Energy in meV if v1 in meter/microsecond - TS_ASSERT_DELTA(xdata[iBin], tof-parser.Eval(), 1e-08); // evaluate the formula on the X-axis - } + for(unsigned int iBin=0; iBin(testWS)); // insert t0-formula as an instrument parameter ModeratorTzero2 alg; - const std::string outputName("outWS"); - initializeAlgorithm(alg,inputName,outputName); + initializeAlgorithm(alg,"testWS","outWS"); TS_ASSERT_THROWS_NOTHING(alg.execute()); - MatrixWorkspace_sptr outws = AnalysisDataService::Instance().retrieveWS(outputName); + // check a few values + double tofH[17]={-6.22573, 993.774, 1993.77, 2993.77, 3993.77, 4993.77, 5993.77, 6993.77, 7993.77, 8992.63, 9992.13, 10991.8170, 11991.5835, 12991.4, 13991.2408, 14991.1, 15990.9881}; + double tofL[10]={143.774, 593.774, 2993.77, 6993.77, 7993.77, 8992.63, 9992.13, 10991.8, 12991.4, 14991.1}; + MatrixWorkspace_sptr outws = AnalysisDataService::Instance().retrieveWS("outWS"); Mantid::DataObjects::EventWorkspace_sptr outWS(boost::dynamic_pointer_cast(outws)); - Mantid::Geometry::IObjComponent_const_sptr sample = testWS->getInstrument()->getSample(); - double L1,L2,E2,v2,t2,t1,v1,E1; - mu::Parser parser; - parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser - parser.SetExpr(formula); - L1=testWS->getInstrument()->getSource()->getDistance(*sample); - const size_t numHists=testWS->getNumberHistograms(); - for(size_t iHist=0; iHistgetNumberHistograms(); ++iHist) { - Mantid::Geometry::IDetector_const_sptr det=testWS->getDetector(iHist); - L2=det->getDistance(*sample); - std::vector< double > wsProp=det->getNumberParameter("Efixed"); - E2=wsProp.at(0); //[E2]=meV - v2=m_convFact*sqrt(E2); - t2=v2/L2; - Mantid::DataObjects::EventList &inevlist=testWS->getEventList(iHist); - Mantid::MantidVec intofs=inevlist.getTofs(); - Mantid::DataObjects::EventList &outevlist=outWS->getEventList(iHist); - Mantid::MantidVec outtofs=outevlist.getTofs(); - for(unsigned int itof=0; itofgetEventList(iHist); + Mantid::MantidVec & x = evlist.dataX(); + for(unsigned int iBin=0; iBinsetX(0, xdata); testWS->setX(1, xdata); + const std::string inputName("testWS"); + AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service return testWS; } -/// Create an event workspace with time-of-flights values normally distributed +/// Create an event workspace containing two detectors and two monitors Mantid::DataObjects::EventWorkspace_sptr createTestEventWorkspace() { - const int numBanks(1); // one bank per instrument - const int numPixels(2); // Two pixels per bank - const int numEvents(200); // number of events per pixel + const int numBanks(1); // one squared bank per instrument + const int numPixels(2); // 2x2=4 pixels per bank + const int numEvents(10); // number of events per pixel + double tofs[10]={150.0, 600.0, 3000.0, 7000.0, 8000.0, 9000.0, 10000.0, 11000.0, 13000.0, 15000.0}; // use this list of TOF at each detector pixel Mantid::DataObjects::EventWorkspace_sptr testWS=WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(numBanks, numPixels); testWS->getAxis(0)->unit()=Mantid::Kernel::UnitFactory::Instance().create("TOF"); - boost::mt19937 rng; - boost::normal_distribution<> nd(6493.0, 724.0); //normal distribution with mean=6493 and standard distribution=724 - boost::variate_generator > var_nor(rng, nd); // random number generator const size_t numHists=testWS->getNumberHistograms(); - ModeratorTzero2 alg; - const double t1min(alg.gett1min()); - double tof; for(size_t iHist=0; iHistgetEventList(iHist); for(int iEvent=0; iEvent params {0.,1000.,16000.}; + alg.setProperty("Params",params); + alg.setProperty("PreserveEvents",true); + TS_ASSERT_THROWS_NOTHING(alg.execute()); return testWS; } -/// Add formula for the emission time from moderator to instrument parameters +/* Add formula for the emission time from moderator to instrument parameters + * Also define instrument as of type "indirect" + * Also define "Efixed" for each detector (but nor for the monitors) + */ void addFormula(const std::string &formula, MatrixWorkspace_sptr workspace) { - //deltaE-mode (should be "indirect") - Mantid::Geometry::Instrument_const_sptr instrument = workspace->getInstrument(); // pointer to the instrument - workspace->instrumentParameters().addString(instrument->getChild(0).get(),"deltaE-mode", "indirect"); - workspace->instrumentParameters().addString(instrument->getChild(0).get(),"t0_formula",formula); + const double efixed=2.082; // energy for analyzers, in meV + workspace->instrumentParameters().addString(workspace->getInstrument()->getComponentID(),"deltaE-mode", "indirect"); + workspace->instrumentParameters().addString(workspace->getInstrument()->getComponentID(),"t0_formula", formula); + for (size_t i=0; igetNumberHistograms(); ++i) + { + Mantid::Geometry::IDetector_const_sptr detector=workspace->getDetector(i); + if( !detector->isMonitor() ) + workspace->instrumentParameters().addDouble(detector->getComponentID(),"Efixed", efixed); + } } /// Avoids writing these over and over