diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h index 7133038e53899bba21368405ee159de44c85f252..0f319b63fc365c54aaabe44964d79ead76449b4e 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h @@ -78,6 +78,7 @@ private: double noiseScale); void replaceAll(std::string &str, const std::string &from, const std::string &to); + void addChopperParameters(API::MatrixWorkspace_sptr &ws); /// A pointer to the random number generator Kernel::PseudoRandomNumberGenerator *m_randGen; diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h index 21e97a89ed4834e2adf04305545c529bf157b406..4ccdffcf5fbf2fae40d1b45f7a30ed294d95db74 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h @@ -102,6 +102,8 @@ protected: // for testing, private otherwise. // the value of constant phase shift on the chopper used to calculate // tof at chopper from recorded delay. double m_phase; + // internal poĆnter to accelerate access to chopper + boost::shared_ptr<const Geometry::IComponent> m_chopper; }; } // namespace Algorithms diff --git a/Code/Mantid/Framework/Algorithms/src/CreateSampleWorkspace.cpp b/Code/Mantid/Framework/Algorithms/src/CreateSampleWorkspace.cpp index 96e9ffe3304873bfd7a06e7ad7e4d61e7a79e4fa..26a5736ee7fc5f83e53b5075e5b5ac62800a247a 100644 --- a/Code/Mantid/Framework/Algorithms/src/CreateSampleWorkspace.cpp +++ b/Code/Mantid/Framework/Algorithms/src/CreateSampleWorkspace.cpp @@ -203,7 +203,6 @@ void CreateSampleWorkspace::exec() { } m_randGen = new Kernel::MersenneTwister(seedValue); } - // Create an instrument with one or more rectangular banks. Instrument_sptr inst = createTestInstrumentRectangular( numBanks, bankPixelWidth, pixelSpacing, bankDistanceFromSample, @@ -221,6 +220,8 @@ void CreateSampleWorkspace::exec() { numBanks * bankPixelWidth * bankPixelWidth, num_bins, xMin, binWidth, bankPixelWidth * bankPixelWidth, inst, functionString, isRandom); } + // add chopper + this->addChopperParameters(ws); // Set the Unit of the X Axis try { @@ -249,6 +250,26 @@ void CreateSampleWorkspace::exec() { setProperty("OutputWorkspace", ws); ; } +/** Add chopper to the existing matrix workspace +@param ws -- shared pointer to existing matrix workspace which has instrument and chopper + +@returns workspace modified to have Fermi chopper added to it. +*/ +void CreateSampleWorkspace::addChopperParameters(API::MatrixWorkspace_sptr &ws){ + + auto testInst = ws->getInstrument(); + auto chopper = testInst->getComponentByName("chopper-position"); + + // add chopper parameters + auto ¶mMap = ws->instrumentParameters(); + const std::string description("The initial rotation phase of the disk used to calculate the time" + " for neutrons arriving at the chopper according to the formula time = delay + initial_phase/Speed"); + paramMap.add<double>("double", chopper.get(),"initial_phase",-3000.,&description); + paramMap.add<std::string>("string",chopper.get(),"ChopperDelayLog","fermi_delay"); + paramMap.add<std::string>("string",chopper.get(),"ChopperSpeedLog","fermi_speed"); + paramMap.add<std::string>("string",chopper.get(),"FilterBaseLog","is_running"); + paramMap.add<bool>("bool",chopper.get(),"filter_with_derivative",false); +} /** Create histogram workspace */ @@ -478,13 +499,13 @@ Instrument_sptr CreateSampleWorkspace::createTestInstrumentRectangular( source->setPos(V3D(0.0, 0.0, -sourceSampleDistance)); testInst->add(source); testInst->markAsSource(source); - // add chopper - ObjComponent *chopper = + + // Add chopper + ObjComponent *chopper = new ObjComponent("chopper-position", Object_sptr(new Object), testInst.get()); chopper->setPos(V3D(0.0, 0.0, -0.25*sourceSampleDistance)); testInst->add(chopper); - auto ID = chopper->getComponentID(); - + // Define a sample as a simple sphere diff --git a/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp b/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp index d6165820765f88d824b9ad0efc4958d69f5eae3d..457492c509d66fbb274d4963669d845c2277e22f 100644 --- a/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp +++ b/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp @@ -27,7 +27,7 @@ namespace Mantid { // half maximal resolution for LET m_max_Eresolution(0.5e-3), m_peakEnergyRatio2reject(0.1), - m_phase(0){} + m_phase(0),m_chopper(){} /// Initialization method. void GetAllEi::init() { @@ -165,11 +165,11 @@ namespace Mantid { auto pInstrument = inputWS->getInstrument(); //auto lastChopPositionComponent = pInstrument->getComponentByName("chopper-position"); //auto chopPoint1 = pInstrument->getChopperPoint(0); ->TODO: BUG! this operation loses parameters map. - auto chopPoint = pInstrument->getComponentByName("chopper-position"); - if(!chopPoint)throw std::runtime_error("Instrument "+pInstrument->getName()+ + m_chopper = pInstrument->getComponentByName("chopper-position"); + if(!m_chopper)throw std::runtime_error("Instrument "+pInstrument->getName()+ " does not have 'chopper-position' component"); - auto phase = chopPoint->getNumberParameter("initial_phase"); + auto phase = m_chopper->getNumberParameter("initial_phase"); if(phase.size()==0){ throw std::runtime_error("Can not find initial_phase parameter" @@ -182,12 +182,22 @@ namespace Mantid { m_phase = phase[0]; std::string filterBase = getProperty("FilterBaseLog"); if(boost::iequals(filterBase, "Defined in IDF")){ - m_FilterLogName = chopPoint->getStringParameter("FilterBaseLog")[0]; - m_FilterWithDerivative = chopPoint->getBoolParameter("filter_with_derivative")[0]; + m_FilterLogName = m_chopper->getStringParameter("FilterBaseLog")[0]; + m_FilterWithDerivative = m_chopper->getBoolParameter("filter_with_derivative")[0]; }else{ m_FilterLogName = filterBase; m_FilterWithDerivative = getProperty("FilterWithDerivative"); } + try{ + auto pTimeSeries = dynamic_cast<Kernel::TimeSeriesProperty<double> *> + (inputWS->run().getProperty(m_FilterLogName)); + }catch(std::runtime_error &){ + g_log.warning()<<" Can not retrieve (double) filtering log: "+m_FilterLogName+" from current workspace\n" + " Using total experiment range to find logs averages for chopper parameters\n"; + m_FilterLogName=""; + m_FilterWithDerivative=false; + } + //auto chopPoint1 = pInstrument->getComponentByName("fermi-chopper"); //auto par = chopPoint1->getDoubleParameter("Delay (us)"); double chopSpeed,chopDelay; @@ -199,8 +209,8 @@ namespace Mantid { auto moderator = pInstrument->getSource(); - double chopDistance = chopPoint->getDistance(*moderator); //location[0].distance(moderator->getPos()); - double velocity = chopDistance/chopDelay; + double chopDistance = m_chopper->getDistance(*moderator); //location[0].distance(moderator->getPos()); + double velocity = chopDistance/chopDelay; // build workspace to find monitor's peaks size_t det1WSIndex; @@ -801,16 +811,14 @@ namespace Mantid { *@param -- propertyName name of the property to find log for * *@return -- pointer to property which contain the log requested or nullptr if - * no log found or othr errors identified. + * no log found or other errors identified. */ Kernel::Property * GetAllEi::getPLogForProperty(const API::MatrixWorkspace_sptr &inputWS, const std::string &propertyName){ std::string LogName = this->getProperty(propertyName); if(boost::iequals(LogName, "Defined in IDF")){ - auto chopper = inputWS->getInstrument()->getComponentByName("chopper-position"); - if(!chopper)return nullptr; - auto AllNames = chopper->getStringParameter(propertyName); + auto AllNames = m_chopper->getStringParameter(propertyName); if(AllNames.size()!=1)return nullptr; LogName = AllNames[0]; } diff --git a/Code/Mantid/docs/source/algorithms/GetAllEi-v1.rst b/Code/Mantid/docs/source/algorithms/GetAllEi-v1.rst index 0d33b3319bd6a38e2c08d00cd68fe0f9e7d62310..67f7f9518fcd0ac22816807bba59a7680a03f9c6 100644 --- a/Code/Mantid/docs/source/algorithms/GetAllEi-v1.rst +++ b/Code/Mantid/docs/source/algorithms/GetAllEi-v1.rst @@ -10,60 +10,84 @@ Description ----------- Algorithm finds the estimate for all incident energies allowed by chopper system of an inelastic instrument and returns a workspace, -with the estimates for positions, heights and width of incident energies provided by the choppers. These estimates can be used as guess value for -:ref:`algm-GetEi` algorithm or as inputs for a peak fitting procedure. +with the estimates for positions, heights and width of incident energies provided by the choppers and registered on monitors. +These estimates can be used as guess value for :ref:`algm-GetEi` algorithm or as inputs for a peak fitting procedure. Algorithm performs number of steps to identify the values requested: -* It takes appropriate log names from insrument definition file (IDF), namely chopper-position component and calculates last chopper speed and delay as average -of the filtered log values. Guess chopper opening times are calculated from chopper speed and delay time. The "chopper-position" component with appopriate properties -has to be present in IDF for this algorithm to work. See ISIS MARI or MAPS instrument definition files for example of "chopper-position" component. +#. It takes appropriate log names from instrument definition file (IDF), namely chopper-position component and calculates last chopper speed and delay as average +of the filtered log values. Guess chopper opening times are calculated from chopper speed and delay time. The "chopper-position" component with appropriate properties +has to be present in IDF for this algorithm to work. See ISIS MARI or MAPS instrument definition files for example of "chopper-position" component. -* Algorithm takes mininal instrument resolution estimate and searches for real peaks around guess values above within 4 sigma of this resolution interval. +#. Algorithm uses estimate for the minimal energy resolution of an instrument and searches for real peaks around guess values obtained +earlier within 4 sigma of this resolution interval. -* If peaks are found, the algorithm performs running averages over signal in the appropriate time interval until first derivative -of the signal has only one zero. This value is accepted as the guess energy postion and the the distance between closest to -the guess energy zeros of the second derivative are accpted as the guess values for the peak width. The peak amplitude -is estimated from the integral intensity of the signal, assuming that the peak shape is Gaussian. +#. If peaks are found, the algorithm performs running averages over signal in the appropriate time interval until first derivative +of the signal has only one zero. This value is returned as the guess energy and the distance between closest to +the guess energy zeros of the second derivative are returned as the guess values for the peak width. The peak amplitude +is estimated from the total intensity of the signal within the search interval, assuming that the peak shape is Gaussian. -* The similar procedure is performed for second monitors. The peak is accepted as real only if the peak width is within the limits of the instrument resolution and -the distance between peaks positions on two monitors (on energy scalse) is smaller then two sigma. +#. Similar procedure is performed for second monitors. The peak is accepted only if the peak width is between the minimal and maximal instrument resolution and +the distance between peaks positions on two monitors (on energy scale) is smaller then two sigma. -Algorithm returns matrix workspace contatingn single sepctra, with x-value representing peak positions, y-values: peak heigths and the error: peak width. +Algorithm returns matrix workspace containing single spectrum, with x-values representing peak positions, y-values: peak heights and the error: peak width. Used Subalgorithms ------------------ The algorithm uses :ref:`Unit Factory <Unit Factory>` and :ref:`algm-ConvertUnits` algorithm to convert units from TOF to energy. -Usage ------ -.. include:: ../usagedata-note.txt - -**Example: Fixing the Ei** - -.. testcode:: fixEi +**Example: Find all incident energies for test workspace** + +.. testcode:: foundAllEi + + # BUILD SAMPLE WORKSPACE + # Build sample workspace with chopper and in energy units to + # have defined peaks in defined energy positions + wsEn=CreateSampleWorkspace(Function='Multiple Peaks', NumBanks=1, BankPixelWidth=2, NumEvents=10000, XUnit='Energy', XMin=10, XMax=200, BinWidth=0.1) + # convert units to TOF to simulate real workspace obtained from experiment + ws = ConvertUnits(InputWorkspace=wsEn, Target='TOF') + # find chopper log values would be present in real workspace + l_chop = 7.5 # chopper position build into test workspace + l_mon1 = 15. # monitor 1 position (detector 1), build into test workspace + t_mon1 = 3100. # the time of flight defined by incident energy of the peak generated by CreateSampelpWorkspace algorithm. + t_chop = (l_chop/l_mon1)*t_mon1 + # Add these log values to simulated workspace to represent real sample logs + AddTimeSeriesLog(ws, Name="fermi_delay", Time="2010-01-01T00:00:00", Value=t_chop ,DeleteExisting=True) + AddTimeSeriesLog(ws, Name="fermi_delay", Time="2010-01-01T00:30:00", Value=t_chop ) + AddTimeSeriesLog(ws, Name="fermi_speed", Time="2010-01-01T00:00:00", Value=900 ,DeleteExisting=True) + AddTimeSeriesLog(ws, Name="fermi_speed", Time="2010-01-01T00:30:00", Value=900) + #------------------------------------------------------------- + + # FIND GUESS PEAKS + allEiWs=GetAllEi(ws,Monitor1SpecID=1,Monitor2SpecID=2) + # Analyze results + allEi = allEiWs.readX(0); + peakHeight = allEiWs.readY(0); + peakWidth = allEiWs.readE(0); + + # Check if peaks positions are indeed correct: + #------------------------------------------------------------- + resEi=[] + for ei_guess in allEi: + nop,t_peak,monIndex,tZero=GetEi(InputWorkspace=ws, Monitor1Spec=1, Monitor2Spec=2, EnergyEstimate=ei_guess) + resEi.append((nop,t_peak)); + print "! Guess Ei ! peak TOF ! peak height ! peak width !" + for ind,val in enumerate(resEi): + print "! {0: >6.2f} ! {1: >6.2f} ! {2: >6.2f} ! {3: >6.2f} !".format(allEi[ind],val[1],peakHeight[ind],peakWidth[ind]) + # + # NOTE: incident energy of GetEi is calculated from distance between monitor 1 and 2, and this distance is not correct in + # the test workspace. The important pint is that getEi can find energies from guess values and TOF for peaks is correct. - ws = CreateSampleWorkspace(bankPixelWidth=1,binWidth=10) - - (ei, firstMonitorPeak, FirstMonitorIndex, tzero) = GetEi(ws,Monitor1Spec=1,Monitor2Spec=2,EnergyEstimate=15.0,FixEi=True) - - print "ei: %.2f" % ei - print "firstMonitorPeak: %.2f" % firstMonitorPeak - print "FirstMonitorIndex: %i" % FirstMonitorIndex - print "tzero: %.2f" % tzero - Output: -.. testoutput:: fixEi +.. testoutput:: foundAllEi :options: +NORMALIZE_WHITESPACE - ei: 15.00 - firstMonitorPeak: 8854.69 - FirstMonitorIndex: 0 - tzero: 0.00 - + ! Guess Ei ! peak TOF ! peak height ! peak width ! + ! 67.05 ! 4188.16 ! 34.68 ! 2.35 ! + ! 124.15 ! 3079.10 ! 14.01 ! 4.35 ! .. categories::