research/transmission_spectroscopy/TOF/analyse.cpp

   1 #include "analyse.h"\r
   2 \r
   3 #include <cmath>\r
   4 \r
   5 using namespace std;\r
   6 \r
   7 \r
   8 \r
   9 int Analysis::FindMaximum(vector<int> & data)\r
  10 {\r
  11         int maxIndex = 0;\r
  12         for (int ii=1; ii < data.size(); ++ii)\r
  13         {\r
  14                 if (data[ii] > data[maxIndex])\r
  15                         maxIndex = ii;\r
  16         }\r
  17         return maxIndex;\r
  18 }\r
  19 \r
  20 void Analysis::AnalyseData(vector<int> & input, map<float, vector<int> > & output, float peakEnergy, float timeRes, float m, float L, float peakCutOff, int spectraNumber)\r
  21 {\r
  22         //Create vector for each spectra\r
  23         vector<vector<int> > inputs;\r
  24         for (int ii=0; ii < spectraNumber; ++ii)\r
  25                 inputs.push_back(vector<int>());\r
  26         \r
  27 \r
  28         int index = 0;\r
  29         int zeroCount = 0;\r
  30         bool nextReady = false;\r
  31         for (int ii=0; ii < input.size() && index < inputs.size(); ++ii)\r
  32         {\r
  33                 if (input[ii] <= 0)\r
  34                 {\r
  35                         zeroCount++;\r
  36                         if (zeroCount > 50)\r
  37                         {\r
  38                                 zeroCount = 0;\r
  39                                 nextReady = true;\r
  40                         }\r
  41                 }\r
  42                 else if (nextReady)\r
  43                 {\r
  44                         index++;\r
  45                         nextReady = false;\r
  46                 }\r
  47                 inputs[index].push_back(input[ii]);\r
  48 \r
  49         }\r
  50 \r
  51         for (int ii=0; ii < inputs.size(); ++ii)\r
  52         {\r
  53                 AnalyseSpectrum(inputs[ii], output, peakEnergy, timeRes, m, L, peakCutOff);\r
  54         }\r
  55         \r
  56 }\r
  57 \r
  58 void Analysis::AnalyseSpectrum(vector<int> & input, map<float, vector<int> > & output, float peakEnergy, float timeRes, float m, float L, float peakCutOff)\r
  59 {\r
  60         int peakBin = FindMaximum(input);\r
  61 \r
  62         if (peakCutOff != 100)\r
  63         {\r
  64                 peakCutOff /= 100; //Convert from percentage to a scaling factor\r
  65 \r
  66                 //Sum all valid values\r
  67                 int sum = 0;\r
  68                 for (int ii=0; ii < input.size(); ++ii)\r
  69                 {\r
  70                         if (input[ii] >= input[peakBin]*peakCutOff)\r
  71                                 sum += input[ii];\r
  72                 }\r
  73 \r
  74                 //Now find the middle\r
  75                 int halfSum = 0;\r
  76                 int ii=0;\r
  77                 for (ii=0; ii < input.size() && halfSum < sum/2; ++ii)\r
  78                 {\r
  79                         if (input[ii] >= input[peakBin]*peakCutOff)\r
  80                                 halfSum += input[ii];\r
  81                 }       \r
  82                 peakBin = ii;   \r
  83         }\r
  84         \r
  85         timeRes *= pow(10.0f,-12.0f); //Convert resolution to s\r
  86         peakEnergy *= (1.602f * pow(10.0f, -19.0f)); //Convert energy to J\r
  87         m = m * 1.783f * pow(10.0f, -30.0f); //Convert mass to kg\r
  88         L = L * pow(10.0f, -3.0f); //Convert length to m\r
  89         \r
  90 \r
  91         for (int ii=0; ii < input.size(); ++ii)\r
  92         {\r
  93 \r
  94                 float dT = (peakBin - ii) * timeRes;\r
  95                 \r
  96                 \r
  97 \r
  98                 float E = 0.5f*m * pow(L/(L*pow(m/(2.0f*peakEnergy), 0.5f) + dT), 2.0f);\r
  99                 E = E / (1.602f * pow(10.0f, -19.0f)); //Convert J to eV\r
 100 \r
 101                 if (E <= 10*peakEnergy || output[E].size() > 0 || input[ii] != 0) //Don't add zero count values above 10*Epeak\r
 102                         output[E].push_back(input[ii]);\r
 103         }\r
 104         \r
 105 }\r