From adab7006cbba4fac5fa04ca1be36bb0f7d859e81 Mon Sep 17 00:00:00 2001 From: Raquel Alvarez Banos Date: Tue, 18 Nov 2014 10:58:57 +0000 Subject: [PATCH] Re #10498 Working on Phase Quad Muon algorithm --- .../inc/MantidAlgorithms/PhaseQuadMuon.h | 18 +- .../Algorithms/src/PhaseQuadMuon.cpp | 228 +++++++++++++++--- 2 files changed, 208 insertions(+), 38 deletions(-) diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PhaseQuadMuon.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PhaseQuadMuon.h index 0404732607e5..cdd2d9b8f183 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PhaseQuadMuon.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PhaseQuadMuon.h @@ -39,7 +39,7 @@ namespace Mantid { public: /// Default constructor - PhaseQuadMuon() : API::Algorithm(), m_muLife(2.19703E-06) {}; + PhaseQuadMuon() : API::Algorithm(), m_muLife(2.19703), m_bigNumber(1e10) {}; /// Destructor virtual ~PhaseQuadMuon() {}; /// Algorithm's name for identification overriding a virtual method @@ -72,12 +72,16 @@ namespace Mantid /// Load the PhaseTable void loadPhaseTable(const std::string& filename); /// TODO + void normaliseAlphas (std::vector& m_histData); + /// Remove exponential decay from input histograms void removeExponentialDecay (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs); - /// TODO - void squash(API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs); - /// Number of histograms (spectra) + /// Remove exponential decay from input histograms + void loseExponentialDecay (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs); + /// Create squashograms + void squash(API::MatrixWorkspace_sptr tempWs, API::MatrixWorkspace_sptr outputWs); + /// Number of input histograms size_t m_nHist; - /// Number of datapoints in each histogram + /// Number of datapoints per histogram size_t m_nData; /// TODO double m_res; @@ -93,8 +97,10 @@ namespace Mantid double m_tau; /// Freq (of silver run) double m_w; - /// Mu life + /// Muon lifetime double m_muLife; + /// Poisson limit + double m_poissonLim; /// TODO double m_bigNumber; /// Vector of histograms data diff --git a/Code/Mantid/Framework/Algorithms/src/PhaseQuadMuon.cpp b/Code/Mantid/Framework/Algorithms/src/PhaseQuadMuon.cpp index ede5cbc37ce2..8d9e1a022e3a 100644 --- a/Code/Mantid/Framework/Algorithms/src/PhaseQuadMuon.cpp +++ b/Code/Mantid/Framework/Algorithms/src/PhaseQuadMuon.cpp @@ -26,7 +26,13 @@ void PhaseQuadMuon::init() declareProperty(new API::WorkspaceProperty("InputWorkspace", "", Direction::Input), "Name of the input workspace"); - declareProperty(new API::WorkspaceProperty("OutputWorkspace","", Direction::Output), + declareProperty(new API::WorkspaceProperty("OutputWorkspace", "", Direction::Output), + "Name of the output workspace" ); + + declareProperty(new API::WorkspaceProperty("OutputWorkspace1","", Direction::Output), + "Name of the output workspace" ); + + declareProperty(new API::WorkspaceProperty("OutputWorkspace2","", Direction::Output), "Name of the output workspace" ); declareProperty(new API::FileProperty("PhaseTable", "", API::FileProperty::Load, ".INF"), @@ -50,16 +56,43 @@ void PhaseQuadMuon::exec() m_nData = inputWs->getSpectrum(0)->readY().size(); m_nHist = inputWs->getNumberHistograms(); + // Create temporary workspace to perform operations on it + API::MatrixWorkspace_sptr tempWs1 = boost::dynamic_pointer_cast( + API::WorkspaceFactory::Instance().create("Workspace2D", inputWs->getNumberHistograms(), inputWs->getSpectrum(0)->readX().size(), inputWs->getSpectrum(0)->readY().size())); + // Create temporary workspace to perform operations on it + API::MatrixWorkspace_sptr tempWs2 = boost::dynamic_pointer_cast( + API::WorkspaceFactory::Instance().create("Workspace2D", inputWs->getNumberHistograms(), inputWs->getSpectrum(0)->readX().size(), inputWs->getSpectrum(0)->readY().size())); + // Create output workspace with two spectra API::MatrixWorkspace_sptr outputWs = getProperty("OutputWorkspace"); outputWs = boost::dynamic_pointer_cast( - API::WorkspaceFactory::Instance().create("Workspace2D", inputWs->getNumberHistograms(), inputWs->getSpectrum(0)->readX().size(), inputWs->getSpectrum(0)->readY().size())); // TODO change number histograms + API::WorkspaceFactory::Instance().create("Workspace2D", 2, inputWs->getSpectrum(0)->readX().size(), inputWs->getSpectrum(0)->readY().size())); outputWs->getAxis(0)->unit() = inputWs->getAxis(0)->unit(); - removeExponentialDecay(inputWs,outputWs); - //squash (inputWs, outputWs); + normaliseAlphas(m_histData); + removeExponentialDecay(inputWs,tempWs1); + loseExponentialDecay(inputWs,tempWs2); + squash (tempWs2, outputWs); - setProperty("OutputWorkspace", outputWs); + setProperty("OutputWorkspace", outputWs); // TODO: change this + setProperty("OutputWorkspace1", tempWs1); // TODO: change this + setProperty("OutputWorkspace2", tempWs2); // TODO: change this +} + +void PhaseQuadMuon::normaliseAlphas (std::vector& histData) +{ + double max=0; + for (size_t h=0; hmax) + { + max = histData[h].alpha; + } + } + for (size_t h=0; hgetNumberHistograms(); h++) + if ( m_nHist != inputWs->getNumberHistograms() ) + { + throw std::runtime_error("InputWorkspace and PhaseTable do not have the same number of spectra"); + } + else { - auto specIn = inputWs->getSpectrum(h); - auto specOut = outputWs->getSpectrum(h); - for (int i=0; ireadY()[i]; - double oops = ( usey<= 0 ) || ( specIn->readE()[i] > m_bigNumber ); - specOut->dataY()[i] = oops ? 0 : log(usey); - specOut->dataE()[i] = oops ? m_bigNumber : specIn->readE()[i] / usey; - std::cout << - } - double s=0; - double sx=0; - double sy=0; + // Muon decay: N(t) = N0 exp(-t/tau) [ 1 + a P(t) ] + // N(t) - N0 exp(-t/tau) = 1 + a P(t) + // N0? + // Int N(t) dt = N0 Int exp(-t/tau) dt - for (int i=0; igetNumberHistograms(); h++) { - double sigma = specOut->dataE()[i] * specOut->dataE()[i]; - s+= 1./sigma; - sx+= specOut->dataX()[i]/sigma; - sy+= specOut->dataY()[i]/sigma; - } - double N0 = (sy+sx/m_muLife/1E9)/s; -// std::cout << specIn->readY()[40] << " " << sy << " " << sx << " " << s << " " << N0 << std::endl; - } + auto specIn = inputWs->getSpectrum(h); + auto specOut = outputWs->getSpectrum(h); -} + specOut->dataX()=specIn->readX(); + double sum1, sum2; + sum1=sum2=0; + for(int i=0; ireadX()[i]>=0 ) + { + sum1 += specIn->dataY()[i]; + sum2 += exp( -((specIn->dataX()[i+1]-specIn->dataX()[i])*0.5+specIn->dataX()[i])/m_muLife); + } + } + double N0 = sum1/sum2; -void PhaseQuadMuon::squash(API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs) + for (int i=0; ireadX()[i]>=0 ) + { + specOut->dataY()[i] = specIn->dataY()[i] - N0*exp(-((specIn->dataX()[i+1]-specIn->dataX()[i])*0.5+specIn->dataX()[i])/m_muLife); + specOut->dataE()[i] = ( specIn->readY()[i] > m_poissonLim) ? specIn->readE()[i] : sqrt(N0*exp(-specIn->readX()[i]/m_muLife)); + } + } + } // Histogram loop + + }// else + +} + +//---------------------------------------------------------------------------------------------- +/** Remove exponential decay from input histograms, i.e., calculate asymmetry +* @param inputWs :: input workspace containing the spectra +* @param outputWs :: output workspace to hold temporary results +*/ +void PhaseQuadMuon::loseExponentialDecay (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs) { if ( m_nHist != inputWs->getNumberHistograms() ) @@ -178,9 +230,121 @@ void PhaseQuadMuon::squash(API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspa } else { - // TODO + + for (size_t h=0; hgetNumberHistograms(); h++) + { + auto specIn = inputWs->getSpectrum(h); + auto specOut = outputWs->getSpectrum(h); + + specOut->dataX()=specIn->readX(); + + for(int i=0; ireadX()[i]>=0 ) + { + double usey = specIn->readY()[i]; + double oops = ( (usey<=0) || (specIn->readE()[i]>=m_bigNumber)); + specOut->dataY()[i] = oops ? 0 : log(usey); + specOut->dataE()[i] = oops ? m_bigNumber : specIn->readE()[i]/usey; + } + } + + double s, sx, sy, sig; + s = sx = sy =0; + for (int i=0; ireadX()[i]>=0 ) + { + sig = specOut->readE()[i]*specOut->readE()[i]; + s += 1./sig; + sx+= specOut->readX()[i]/sig; + sy+= specOut->readY()[i]/sig; + } + } + double N0 = (sy+sx/m_muLife)/s; + N0=exp(N0); + + for (int i=0; ireadX()[i]>=0 ) + { + specOut->dataY()[i] = specIn->dataY()[i] - N0 *exp(-specIn->readX()[i]/m_muLife); + specOut->dataE()[i] = ( specIn->readY()[i] > m_poissonLim) ? specIn->readE()[i] : sqrt(N0*exp(-specIn->readX()[i]/m_muLife)); + } + } + } // Histogram loop + + }// else + +} + + +void PhaseQuadMuon::squash(API::MatrixWorkspace_sptr tempWs, API::MatrixWorkspace_sptr outputWs) +{ + + double sxx=0; + double syy=0; + double sxy=0; + + for (size_t h=0; hgetNumberHistograms(); h++) + { + auto data = m_histData[h]; + double X = data.n0 * data.alpha * cos(data.phi); + double Y = data.n0 * data.alpha * sin(data.phi); + sxx += X*X; + syy += Y*Y; + sxy += X*Y; } + double lam1 = 2 * syy / (sxx*syy - sxy*sxy); + double mu1 = 2 * sxy / (sxy*sxy - sxx*syy); + double lam2 = 2 * sxy / (sxy*sxy - sxx*syy); + double mu2 = 2 * sxx / (sxx*syy - sxy*sxy); + std::vector aj, bj; + + for (size_t h=0; h data1(m_nData,0), data2(m_nData,0); + std::vector sigm1(m_nData,0), sigm2(m_nData,0); + for (size_t i=0; igetSpectrum(h); + data1[i] += aj[h] * spec->readY()[i]; + data2[i] += bj[h] * spec->readY()[i]; + sigm1[i] += aj[h]*aj[h] * spec->readE()[i] * spec->readE()[i]; + sigm2[i] += bj[h]*bj[h] * spec->readE()[i] * spec->readE()[i]; + } + sigm1[i] = sqrt(sigm1[i]); + sigm2[i] = sqrt(sigm2[i]); + } + + for (size_t i=0; igetSpectrum(0)->readX()[i]; + double xp= tempWs->getSpectrum(0)->readX()[i+1]; + double e = exp(-( (xp-x)*0.5+x )/m_muLife); + data1[i] /= e; + data2[i] /= e; + sigm1[i] /= e; + sigm2[i] /= e; + } + + outputWs->getSpectrum(0)->dataX() = tempWs->getSpectrum(0)->readX(); + outputWs->getSpectrum(0)->dataY() = data1; + outputWs->getSpectrum(0)->dataE() = sigm1; + outputWs->getSpectrum(1)->dataX() = tempWs->getSpectrum(1)->readX(); + outputWs->getSpectrum(1)->dataY() = data2; + outputWs->getSpectrum(1)->dataE() = sigm2; + }