diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h
index 042371ba4e150994917722a0d77e52e9df69bd1d..d79bac83a162b8193dc033187ac6e2aca684f356 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h
@@ -87,7 +87,9 @@ protected: // for testing, private otherwise.
/// if log, which identifies that instrument is running is available on workspace.
/// The log should be positive when instrument is running and negative or 0 otherwise.
- bool m_useFilterLog;
+ bool m_useFilterLog;
+ /// if true, take derivate of the filter log to identify interval when instrument is running.
+ bool m_FilterWithDerivative;
/// maximal relative peak width to consider acceptable. Defined by minimal instrument resolution
/// and does not exceed 0.08
double m_min_Eresolution;
diff --git a/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp b/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp
index c5783578fc0caecc79780b4d73e628ff74abadbc..d29e58432dbb793042221bf4ca4832e8861d667b 100644
--- a/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/GetAllEi.cpp
@@ -19,7 +19,9 @@ namespace Mantid {
/// Empty default constructor
GetAllEi::GetAllEi() : Algorithm(),
- m_useFilterLog(true),m_min_Eresolution(0.08),
+ m_useFilterLog(true),m_FilterWithDerivative(true),
+ // minimal resolution for all instruments
+ m_min_Eresolution(0.08),
// half maximal resolution for LET
m_max_Eresolution(0.5e-3),
m_peakEnergyRatio2reject(0.1){}
@@ -52,6 +54,13 @@ namespace Mantid {
" chopper speed and chopper delay which matter.\n"
"If such log is not present, log values are calculated "
"from experiment start&end times.");
+ declareProperty("FilterWithDerivative", true, "Use derivative of 'FilterBaseLog' "
+ "rather then log values itself to filter invalid time intervals.\n"
+ "Invalid values are then the "
+ "values where the derivative of the log turns zero.\n"
+ "E.g. would be 'proton_chage' log which grows for each frame "
+ "when instrument is counting and constant otherwise.");
+
auto maxInRange = boost::make_shared<Kernel::BoundedValidator<double>>();
maxInRange->setLower(1.e-6);
@@ -121,11 +130,21 @@ namespace Mantid {
}
};
+ struct peakKeeper2{
+ double left_rng;
+ double right_rng;
+ peakKeeper2(){};
+ peakKeeper2(double left,double right):
+ left_rng(left),right_rng(right)
+ {}
+ };
+
+
/** Executes the algorithm -- found all existing monitor peaks. */
void GetAllEi::exec() {
// Get pointers to the workspace, parameter map and table
API::MatrixWorkspace_sptr inputWS = getProperty("Workspace");
-
+ m_FilterWithDerivative = getProperty("FilterWithDerivative");
m_min_Eresolution = getProperty("MinInstrResolution");
m_max_Eresolution = getProperty("MaxInstrResolution");
m_peakEnergyRatio2reject = getProperty("PeaksRatioToReject");
@@ -333,7 +352,7 @@ namespace Mantid {
double SmoothRange = 2*maxSigma;
std::vector<double> SAvrg,binsAvrg;
- Kernel::VectorHelper::smoothAtNPoints(S,SAvrg,SmoothRange,&X,
+ Kernel::VectorHelper::smoothInRange(S,SAvrg,SmoothRange,&X,
ind_min,ind_max,&binsAvrg);
double realPeakPos; // this position is less shifted due to the skew in averaging formula
@@ -349,7 +368,7 @@ namespace Mantid {
size_t ic(0);
while((nPeaks>1 || nHills>2)&&(nPrevHills!=nHills) && ic<50){
- Kernel::VectorHelper::smoothAtNPoints(SAvrg,SAvrg1,SmoothRange,&binsAvrg,
+ Kernel::VectorHelper::smoothInRange(SAvrg,SAvrg1,SmoothRange,&binsAvrg,
0,ind_max-ind_min,&binsAvrg1);
nPrevHills=nHills;
@@ -497,11 +516,12 @@ namespace Mantid {
std::vector<size_t> & irangeMax, std::vector<bool> &guessValid){
- size_t index_min,index_max;
size_t nBins = eBins.size();
guessValid.resize(guess_energy.size());
- // get bin range corresponding to energy range
- auto getBinRange=[&](double eMin,double eMax){
+ // get bin range corresponding to the energy range
+ auto getBinRange=[&](double eMin,double eMax,
+ size_t &index_min,size_t &index_max)
+ {
if(eMin<=eBins[0]){
index_min =0;
}else{
@@ -509,18 +529,18 @@ namespace Mantid {
}
if(eMax>=eBins[nBins-1]){
- index_max =nBins+1;
+ index_max =nBins-1;
}else{
index_max = Kernel::VectorHelper::getBinIndex(eBins,eMax)+1;
- if(index_max>=nBins)index_max=nBins-1;
+ if(index_max>=nBins)index_max=nBins-1; // last bin range anyway. Should not happen
}
-
};
// refine bin range. May need better procedure for this.
- auto refineEGuess =[&](double &eGuess){
+ auto refineEGuess =[&](double &eGuess,size_t index_min,size_t index_max)
+ {
size_t ind_Emax=index_min;
double SMax(0);
- for(size_t i=index_min;i<index_max-1;i++){
+ for(size_t i=index_min;i<index_max;i++){
double dX = eBins[i+1]-eBins[i];
double sig = signal[i]/dX;
if(sig>SMax){
@@ -538,16 +558,38 @@ namespace Mantid {
};
// Do the job
+ size_t ind_min,ind_max;
+ irangeMin.resize(0);
+ irangeMax.resize(0);
+
+ // identify guess bin ranges
+ std::vector<peakKeeper2> guess_peak(guess_energy.size());
+ for(size_t nGuess = 0;nGuess<guess_energy.size();nGuess++){
+ double eGuess = guess_energy[nGuess];
+ getBinRange(eGuess*(1-4*Eresolution),eGuess*(1+4*Eresolution),ind_min,ind_max);
+ guess_peak[nGuess] = peakKeeper2(eBins[ind_min],eBins[ind_max]);
+ }
+ // verify that the ranges not intercept and refine interceptions
+ for(size_t i = 1;i<guess_energy.size();i++){
+ if(guess_peak[i-1].right_rng>guess_peak[i].left_rng){
+ double mid_pnt = 0.5*(guess_peak[i-1].right_rng+guess_peak[i].left_rng);
+ guess_peak[i-1].right_rng = mid_pnt;
+ guess_peak[i].left_rng = mid_pnt;
+ }
+ }
+ // identify final bin ranges
for(size_t nGuess = 0;nGuess<guess_energy.size();nGuess++){
double eGuess = guess_energy[nGuess];
- getBinRange(eGuess*(1-4*Eresolution),eGuess*(1+4*Eresolution));
+ getBinRange(guess_peak[nGuess].left_rng,guess_peak[nGuess].right_rng,ind_min,ind_max);
- if(refineEGuess(eGuess)){
+ if(refineEGuess(eGuess,ind_min,ind_max)){
- getBinRange(eGuess*(1-3*Eresolution),eGuess*(1+3*Eresolution));
- irangeMin.push_back(index_min);
- irangeMax.push_back(index_max);
+ getBinRange(std::max(guess_peak[nGuess].left_rng,eGuess*(1-3*Eresolution)),
+ std::max(guess_peak[nGuess].right_rng,eGuess*(1+3*Eresolution)),
+ ind_min,ind_max);
+ irangeMin.push_back(ind_min);
+ irangeMax.push_back(ind_max);
guessValid[nGuess]=true;
}else{
guessValid[nGuess]=false;
diff --git a/Code/Mantid/Framework/Algorithms/test/GetAllEiTest.h b/Code/Mantid/Framework/Algorithms/test/GetAllEiTest.h
index fbe0b79628937974d99737fa0304d6a27fc3b04d..264d2cc8197e9192a89834d50a1358ce8498ded7 100644
--- a/Code/Mantid/Framework/Algorithms/test/GetAllEiTest.h
+++ b/Code/Mantid/Framework/Algorithms/test/GetAllEiTest.h
@@ -33,17 +33,19 @@ public:
API::MatrixWorkspace_sptr buildWorkspaceToFit(const API::MatrixWorkspace_sptr &inputWS, size_t &wsIndex0){
return GetAllEi::buildWorkspaceToFit(inputWS,wsIndex0);
}
- void findBinRanges(const MantidVec & eBins,
- const std::vector<double> & guess_energy,double Eresolution,std::vector<size_t> & irangeMin,
- std::vector<size_t> & irangeMax){
- GetAllEi::findBinRanges(eBins,guess_energy,Eresolution,irangeMin,irangeMax);
+ void findBinRanges(const MantidVec & eBins,const MantidVec & signal,
+ const std::vector<double> & guess_energies,
+ double Eresolution,std::vector<size_t> & irangeMin,
+ std::vector<size_t> & irangeMax, std::vector<bool> &guessValid){
+ GetAllEi::findBinRanges(eBins,signal,guess_energies,Eresolution,irangeMin,
+ irangeMax,guessValid);
}
void setResolution(double newResolution){
this->m_max_Eresolution = newResolution;
}
size_t calcDerivativeAndCountZeros(const std::vector<double> &bins,const std::vector<double> &signal,
- std::vector<double> &deriv){
- return GetAllEi::calcDerivativeAndCountZeros(bins,signal,deriv);
+ std::vector<double> &deriv,std::vector<double> &zeros){
+ return GetAllEi::calcDerivativeAndCountZeros(bins,signal,deriv,zeros);
}
};
@@ -258,36 +260,15 @@ public:
// TS_ASSERT_DELTA(Xsp1[i],Xsp2[i],1.e-6);
//}
}
- void test_binRanges(){
- std::vector<size_t> bin_min,bin_max;
- //Index 0 1 2 3 4 5 6 7 8 9 10 11 12 13
- double debin[] = {1,2,3,4,5,6,7,8,9,10,11,12,13,15};
- std::vector<double> ebin(debin,debin+sizeof(debin)/sizeof(double));
-
- double dguess[]= {1,2,10,15};
- std::vector<double> guess(dguess,dguess+sizeof(dguess)/sizeof(double));
- m_getAllEi.findBinRanges(ebin,guess,0.01,bin_min,bin_max);
-
- TS_ASSERT_EQUALS(bin_min.size(),4)
- TS_ASSERT_EQUALS(bin_max.size(),4)
-
- TS_ASSERT_EQUALS(bin_min[0],0);
- TS_ASSERT_EQUALS(bin_max[0],1);
- TS_ASSERT_EQUALS(bin_min[1],0);
- TS_ASSERT_EQUALS(bin_max[1],2);
- TS_ASSERT_EQUALS(bin_min[2],8);
- TS_ASSERT_EQUALS(bin_max[2],10);
- TS_ASSERT_EQUALS(bin_min[3],12);
- TS_ASSERT_EQUALS(bin_max[3],13);
- }
void test_calcDerivative(){
double sig[]={1,2,3,4,5,6};
std::vector<double> signal(sig,sig+sizeof(sig)/sizeof(double));
double bin[]={2,3,4,5,6,7,8};
std::vector<double> bins(bin,bin+sizeof(bin)/sizeof(double));
+ std::vector<double> zeros;
std::vector<double> deriv;
- size_t nZer = m_getAllEi.calcDerivativeAndCountZeros(bins,signal,deriv);
+ size_t nZer = m_getAllEi.calcDerivativeAndCountZeros(bins,signal,deriv,zeros);
TS_ASSERT_EQUALS(nZer,0);
TS_ASSERT_DELTA(deriv[0],deriv[1],1.e-9);
TS_ASSERT_DELTA(deriv[0],deriv[5],1.e-9);
@@ -296,7 +277,7 @@ public:
double bin1[]={0,1,3,6,10,15,21};
std::vector<double> bins1(bin1,bin1+sizeof(bin1)/sizeof(double));
- nZer = m_getAllEi.calcDerivativeAndCountZeros(bins1,signal,deriv);
+ nZer = m_getAllEi.calcDerivativeAndCountZeros(bins1,signal,deriv,zeros);
TS_ASSERT_EQUALS(nZer,0);
TS_ASSERT_DELTA(deriv[0],deriv[1],1.e-9);
TS_ASSERT_DELTA(deriv[0],deriv[5],1.e-9);
@@ -311,16 +292,70 @@ public:
for(size_t i=0;i<100;i++){
signal[i]=std::sin(0.5*(bins[i]+bins[i+1]));
}
- nZer = m_getAllEi.calcDerivativeAndCountZeros(bins,signal,deriv);
+ nZer = m_getAllEi.calcDerivativeAndCountZeros(bins,signal,deriv,zeros);
TS_ASSERT_EQUALS(nZer,3);
for(size_t i=0;i<99;i++){ // intentionally left boundary point -- its accuracy is much lower
TSM_ASSERT_DELTA("At i="+std::to_string(i),deriv[i],10.*std::cos(0.5*(bins[i]+bins[i+1])),1.e-1);
}
+ TS_ASSERT_DELTA(zeros[0],1.65,1.e-3);
+ TS_ASSERT_DELTA(zeros[1],4.75,1.e-3);
+ TS_ASSERT_DELTA(zeros[2],7.95,1.e-3);
+ }
+ void test_binRanges(){
+ std::vector<size_t> bin_min,bin_max,zeros;
+ //Index 0 1 2 3 4 5 6 7 8 9 10 11 12 13
+ double debin[] = {1,2,3,4,5,6,7,8,9,10,11,12,13,15};
+ std::vector<double> ebin(debin,debin+sizeof(debin)/sizeof(double));
+ // Not yet supported in VC 2012
+ //std::vector<double> ebin={1,2,3,4,5,6,7,8,9,10,11,12,13,15};
+ //Ind 0 1 2 3 4 5 6 7 8 9 10 11 12 13
+ double sig[] = {0,0,0,3,0,0,4,0,0,0,11,0,0};
+ std::vector<double> signal(sig,sig+sizeof(sig)/sizeof(double));
+
+ double dguess[]= {1,6,10,12};
+ std::vector<double> guess(dguess,dguess+sizeof(dguess)/sizeof(double));
+ std::vector<bool> guessValid;
+
+ m_getAllEi.findBinRanges(ebin,signal,guess,0.1,bin_min,bin_max,guessValid);
+
+ TS_ASSERT_EQUALS(bin_min.size(),2)
+ TS_ASSERT_EQUALS(bin_max.size(),2)
+ TS_ASSERT_EQUALS(guessValid.size(),4)
+ TS_ASSERT_EQUALS(bin_min[0],4)
+ TS_ASSERT_EQUALS(bin_max[0],9)
+ TS_ASSERT_EQUALS(bin_min[1],7)
+ TS_ASSERT_EQUALS(bin_max[1],13)
+
+
+ signal[10]=0;
+ signal[11]=11;
+ guess[1] = 3;
+ guess[2]=6;
+ guess[3]=11;
+ m_getAllEi.findBinRanges(ebin,signal,guess,0.01,bin_min,bin_max,guessValid);
+ TS_ASSERT_EQUALS(bin_min.size(),3)
+ TS_ASSERT_EQUALS(bin_max.size(),3)
+ TS_ASSERT_EQUALS(guessValid.size(),4)
+
+ TS_ASSERT_EQUALS(bin_min[0],3);
+ TS_ASSERT_EQUALS(bin_max[0],4);
+ TS_ASSERT(guessValid[1]);
+
+ TS_ASSERT_EQUALS(bin_min[1],6);
+ TS_ASSERT_EQUALS(bin_max[1],7);
+ TS_ASSERT(guessValid[2]);
+
+ TS_ASSERT_EQUALS(bin_min[2],11);
+ TS_ASSERT_EQUALS(bin_max[2],12);
+ TS_ASSERT(guessValid[3]);
+
+ TS_ASSERT(!guessValid[0]);
}
+
private:
GetAllEiTester m_getAllEi;
diff --git a/Code/Mantid/Framework/Kernel/inc/MantidKernel/VectorHelper.h b/Code/Mantid/Framework/Kernel/inc/MantidKernel/VectorHelper.h
index c791d5057ce2ff62f8e12d7fc40f3d8325f62927..e455cb5edeca1ae304641c0021d90edcb2cb3366 100644
--- a/Code/Mantid/Framework/Kernel/inc/MantidKernel/VectorHelper.h
+++ b/Code/Mantid/Framework/Kernel/inc/MantidKernel/VectorHelper.h
@@ -82,13 +82,14 @@ MANTID_KERNEL_DLL int getBinIndex(const std::vector<double> &bins,
MANTID_KERNEL_DLL void linearlyInterpolateY(const std::vector<double> &x,
std::vector<double> &y,
const double stepSize);
-// Do n-point running average of input vector considering bin-boundaries if provided
-MANTID_KERNEL_DLL void smoothAtNPoints(const std::vector<double> &input,
+// Do running average of input vector within specified range, considering heterogeneous bin-boundaries
+// if such boundaries are provided
+MANTID_KERNEL_DLL void smoothInRange(const std::vector<double> &input,
std::vector<double> &output,
double avrgInterval,
- std::vector<double> const * const binBoundaris =NULL,
+ std::vector<double> const * const binBoundaris =nullptr,
size_t startIndex = 0,size_t endIndex = 0,
- std::vector<double> * const outputBinBoundaries=NULL);
+ std::vector<double> * const outputBinBoundaries=nullptr);
//-------------------------------------------------------------------------------------
/** Return the length of the vector (in the physical sense),
diff --git a/Code/Mantid/Framework/Kernel/src/VectorHelper.cpp b/Code/Mantid/Framework/Kernel/src/VectorHelper.cpp
index 88b212d8dc83115af253c37175cfb8a8f3972dec..3be67dd244a68aa734905ec950fd5cfb30807f8f 100644
--- a/Code/Mantid/Framework/Kernel/src/VectorHelper.cpp
+++ b/Code/Mantid/Framework/Kernel/src/VectorHelper.cpp
@@ -534,8 +534,10 @@ void linearlyInterpolateY(const std::vector<double> &x, std::vector<double> &y,
step++;
}
}
-/**Basic running averaging procedure with possibility to account for
- * variable bin ranges.
+/** Basic running average of input vector within specified range, considering variable bin-boundaries
+ * if such boundaries are provided.
+ * The algorithm is trapezium integration, so some peak shift related to first derivative of the
+ * integrated function does happen.
*
* @param input:: input vector to smooth
* @param output:: resulting vector (can not coincide with input)
@@ -554,7 +556,7 @@ void linearlyInterpolateY(const std::vector<double> &x, std::vector<double> &y,
* @param outBins :: if present, pointer to a vector to return
* bin boundaries for output array.
*/
-void smoothAtNPoints(const std::vector<double> &input,
+void smoothInRange(const std::vector<double> &input,
std::vector<double> &output,
double avrgInterval,
std::vector<double> const * const binBndrs,
diff --git a/Code/Mantid/Framework/Kernel/test/VectorHelperTest.h b/Code/Mantid/Framework/Kernel/test/VectorHelperTest.h
index 498ebc5cd2ee0f6748febd495f8ce2c0f08aae47..f985d51d31e3ab1787d6564cd0a848a6e332eeb3 100644
--- a/Code/Mantid/Framework/Kernel/test/VectorHelperTest.h
+++ b/Code/Mantid/Framework/Kernel/test/VectorHelperTest.h
@@ -275,9 +275,9 @@ public:
std::vector<double> inputBoundaries(ib,ib+sizeof(ib)/sizeof(double));
std::vector<double> output;
- TS_ASSERT_THROWS(VectorHelper::smoothAtNPoints(inputData,output,6,&inputBoundaries),std::invalid_argument);
+ TS_ASSERT_THROWS(VectorHelper::smoothInRange(inputData,output,6,&inputBoundaries),std::invalid_argument);
inputBoundaries.push_back(6);
- VectorHelper::smoothAtNPoints(inputData,output,6,&inputBoundaries);
+ VectorHelper::smoothInRange(inputData,output,6,&inputBoundaries);
TS_ASSERT_DELTA(output[1]-output[0],0.492,1.e-3);
TS_ASSERT_DELTA(output[3]-output[2],0.4545,1.e-3);
@@ -288,14 +288,14 @@ public:
inputBoundaries[4]=10;
inputBoundaries[5]=15;
inputBoundaries[6]=21;
- VectorHelper::smoothAtNPoints(inputData,output,6,&inputBoundaries);
+ VectorHelper::smoothInRange(inputData,output,6,&inputBoundaries);
TS_ASSERT_DELTA(output[2],3.,1.e-8);
TS_ASSERT_DELTA(output[0],1.,1.e-8);
TS_ASSERT_DELTA(output[5],6.,1.e-8);
std::vector<double> out_bins;
- VectorHelper::smoothAtNPoints(inputData,output,3,&inputBoundaries,1,5,&out_bins);
+ VectorHelper::smoothInRange(inputData,output,3,&inputBoundaries,1,5,&out_bins);
TS_ASSERT_EQUALS(output.size(),4);
TS_ASSERT_DELTA(output[1],3.,1.e-8);
}
@@ -328,7 +328,7 @@ public:
TS_ASSERT(iLeft<fMax);
TS_ASSERT(iRight<fMax);
- VectorHelper::smoothAtNPoints(inputData,output,10,&inputBoundaries);
+ VectorHelper::smoothInRange(inputData,output,10,&inputBoundaries);
fMax = output[indOfMax]/(inputBoundaries[indOfMax+1]-inputBoundaries[indOfMax]);
iLeft = inputData[indOfMax-1]/(inputBoundaries[indOfMax]-inputBoundaries[indOfMax-1]);
iRight = inputData[indOfMax+1]/(inputBoundaries[indOfMax+2]-inputBoundaries[indOfMax+1]);
@@ -337,7 +337,7 @@ public:
//TS_ASSERT(iRight<fMax);
output.swap(inputData);
- VectorHelper::smoothAtNPoints(inputData,output,10,&inputBoundaries);
+ VectorHelper::smoothInRange(inputData,output,10,&inputBoundaries);
fMax = output[indOfMax]/(inputBoundaries[indOfMax+1]-inputBoundaries[indOfMax]);
iLeft = inputData[indOfMax-1]/(inputBoundaries[indOfMax]-inputBoundaries[indOfMax-1]);
@@ -348,7 +348,7 @@ public:
output.swap(inputData);
- VectorHelper::smoothAtNPoints(inputData,output,10,&inputBoundaries);
+ VectorHelper::smoothInRange(inputData,output,10,&inputBoundaries);
fMax = output[indOfMax]/(inputBoundaries[indOfMax+1]-inputBoundaries[indOfMax]);
iLeft = inputData[indOfMax-1]/(inputBoundaries[indOfMax]-inputBoundaries[indOfMax-1]);
@@ -358,7 +358,7 @@ public:
//TS_ASSERT(iRight<fMax);
output.swap(inputData);
- VectorHelper::smoothAtNPoints(inputData,output,10,&inputBoundaries);
+ VectorHelper::smoothInRange(inputData,output,10,&inputBoundaries);
fMax = output[indOfMax]/(inputBoundaries[indOfMax+1]-inputBoundaries[indOfMax]);
iLeft = inputData[indOfMax-1]/(inputBoundaries[indOfMax]-inputBoundaries[indOfMax-1]);