Commit before breaking everything

[matches/honours.git] / research / transmission_spectroscopy / TOF / analyse.cpp

diff --git a/research/transmission_spectroscopy/TOF/analyse.cpp b/research/transmission_spectroscopy/TOF/analyse.cpp
new file mode 100644 (file)
index 0000000..6976501
--- /dev/null
+++ b/research/transmission_spectroscopy/TOF/analyse.cpp
@@ -0,0 +1,105 @@
+#include "analyse.h"\r
+\r
+#include <cmath>\r
+\r
+using namespace std;\r
+\r
+\r
+\r
+int Analysis::FindMaximum(vector<int> & data)\r
+{\r
+       int maxIndex = 0;\r
+       for (int ii=1; ii < data.size(); ++ii)\r
+       {\r
+               if (data[ii] > data[maxIndex])\r
+                       maxIndex = ii;\r
+       }\r
+       return maxIndex;\r
+}\r
+\r
+void Analysis::AnalyseData(vector<int> & input, map<float, vector<int> > & output, float peakEnergy, float timeRes, float m, float L, float peakCutOff, int spectraNumber)\r
+{\r
+       //Create vector for each spectra\r
+       vector<vector<int> > inputs;\r
+       for (int ii=0; ii < spectraNumber; ++ii)\r
+               inputs.push_back(vector<int>());\r
+       \r
+\r
+       int index = 0;\r
+       int zeroCount = 0;\r
+       bool nextReady = false;\r
+       for (int ii=0; ii < input.size() && index < inputs.size(); ++ii)\r
+       {\r
+               if (input[ii] <= 0)\r
+               {\r
+                       zeroCount++;\r
+                       if (zeroCount > 50)\r
+                       {\r
+                               zeroCount = 0;\r
+                               nextReady = true;\r
+                       }\r
+               }\r
+               else if (nextReady)\r
+               {\r
+                       index++;\r
+                       nextReady = false;\r
+               }\r
+               inputs[index].push_back(input[ii]);\r
+\r
+       }\r
+\r
+       for (int ii=0; ii < inputs.size(); ++ii)\r
+       {\r
+               AnalyseSpectrum(inputs[ii], output, peakEnergy, timeRes, m, L, peakCutOff);\r
+       }\r
+       \r
+}\r
+\r
+void Analysis::AnalyseSpectrum(vector<int> & input, map<float, vector<int> > & output, float peakEnergy, float timeRes, float m, float L, float peakCutOff)\r
+{\r
+       int peakBin = FindMaximum(input);\r
+\r
+       if (peakCutOff != 100)\r
+       {\r
+               peakCutOff /= 100; //Convert from percentage to a scaling factor\r
+\r
+               //Sum all valid values\r
+               int sum = 0;\r
+               for (int ii=0; ii < input.size(); ++ii)\r
+               {\r
+                       if (input[ii] >= input[peakBin]*peakCutOff)\r
+                               sum += input[ii];\r
+               }\r
+\r
+               //Now find the middle\r
+               int halfSum = 0;\r
+               int ii=0;\r
+               for (ii=0; ii < input.size() && halfSum < sum/2; ++ii)\r
+               {\r
+                       if (input[ii] >= input[peakBin]*peakCutOff)\r
+                               halfSum += input[ii];\r
+               }       \r
+               peakBin = ii;   \r
+       }\r
+       \r
+       timeRes *= pow(10.0f,-12.0f); //Convert resolution to s\r
+       peakEnergy *= (1.602f * pow(10.0f, -19.0f)); //Convert energy to J\r
+       m = m * 1.783f * pow(10.0f, -30.0f); //Convert mass to kg\r
+       L = L * pow(10.0f, -3.0f); //Convert length to m\r
+       \r
+\r
+       for (int ii=0; ii < input.size(); ++ii)\r
+       {\r
+\r
+               float dT = (peakBin - ii) * timeRes;\r
+               \r
+               \r
+\r
+               float E = 0.5f*m * pow(L/(L*pow(m/(2.0f*peakEnergy), 0.5f) + dT), 2.0f);\r
+               E = E / (1.602f * pow(10.0f, -19.0f)); //Convert J to eV\r
+\r
+               if (E <= 10*peakEnergy || output[E].size() > 0 || input[ii] != 0) //Don't add zero count values above 10*Epeak\r
+                       output[E].push_back(input[ii]);\r
+       }\r
+       \r
+}\r