summaryrefslogtreecommitdiff
path: root/process_mags.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'process_mags.cpp')
-rw-r--r--process_mags.cpp155
1 files changed, 155 insertions, 0 deletions
diff --git a/process_mags.cpp b/process_mags.cpp
new file mode 100644
index 0000000..ab7112f
--- /dev/null
+++ b/process_mags.cpp
@@ -0,0 +1,155 @@
+
+#include <cinttypes>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <iostream>
+#include <iomanip>
+
+signed M(unsigned N, unsigned i) {
+ return -2 * ((signed)(N - 1) / 2 - (signed)i);
+}
+
+int main(int argc, char *argv[]) {
+ std::vector<std::string> all_args;
+
+ if (argc > 1) {
+ all_args.assign(argv + 1, argv + argc);
+ }
+
+ for (std::string filename : all_args) {
+ std::ifstream in_file;
+ in_file.open(filename);
+
+ if (!in_file.is_open()) {
+ std::cout << "Failed to open " << filename << ".\n";
+ } else {
+
+ std::vector<uint64_t> dat;
+
+ std::string line;
+ while (std::getline(in_file, line)) {
+ std::istringstream iss(line);
+ uint64_t a;
+ while (iss >> a) {
+ dat.push_back(a);
+ }
+ }
+
+ in_file.close();
+
+ double Mt = 0;
+ double δMt = 0;
+
+ double Ma = 0;
+ double δMa = 0;
+
+ double Me = 0;
+ double δMe = 0;
+
+ double Mm = 0;
+ double δMm = 0;
+
+ double χt = 0;
+ double δχt = 0;
+
+ double χa = 0;
+ double δχa = 0;
+
+ double χe = 0;
+ double δχe = 0;
+
+ double χm = 0;
+ double δχm = 0;
+
+ uint64_t normt = 0;
+ uint64_t norme = 0;
+ uint64_t normm = 0;
+
+ for (unsigned i = 0; i < dat.size(); i++) {
+ int64_t MM = M(dat.size(), i);
+ int64_t ni = dat[i];
+ normt += ni;
+
+ Mt += ni * MM;
+ δMt += ni * pow(MM, 2);
+
+ Ma += ni * std::abs(MM);
+ δMa += ni * pow(MM, 2);
+
+ if (i + 1 >= (dat.size()+1) / 2) {
+ norme += ni;
+ Me += ni * MM;
+ δMe += ni * pow(MM, 2);
+ }
+
+ if (i + 1 <= (dat.size()+1) / 2) {
+ normm += ni;
+ Mm += ni * MM;
+ δMm += ni * pow(MM, 2);
+ }
+ }
+
+ Mt /= normt;
+ Ma /= normt;
+ Me /= norme;
+ Mm /= normm;
+
+ δMt = sqrt(δMt) / normt;
+ δMa = sqrt(δMa) / normt;
+ δMe = sqrt(δMe) / norme;
+ δMm = sqrt(δMm) / normm;
+
+ for (unsigned i = 0; i < dat.size(); i++) {
+ int64_t MM = M(dat.size(), i);
+ int64_t ni = dat[i];
+
+ χt += ni * MM * (MM - Mt);
+ δχt += ni * pow(MM, 2) * (pow(MM - Mt, 2) + ni * pow(δMt, 2));
+
+ χa += ni * std::abs(MM) * (std::abs(MM) - Ma);
+ δχa += ni * pow(MM, 2) * (pow(std::abs(MM) - Ma, 2) + ni * pow(δMa, 2));
+
+ if (i + 1 >= (dat.size()+1) / 2) {
+ χe += ni * MM * (MM - Me);
+ δχe += ni * pow(MM, 2) * (pow(MM - Me, 2) + ni * pow(δMe, 2));
+ }
+
+ if (i + 1 <= (dat.size()+1) / 2) {
+ χm += ni * MM * (MM - Mm);
+ δχm += ni * pow(MM, 2) * (pow(MM - Mm, 2) + ni * pow(δMm, 2));
+ }
+ }
+
+ χt /= normt;
+ χa /= normt;
+ χe /= norme;
+ χm /= normm;
+
+ δχt = sqrt(δχt) / normt;
+ δχa = sqrt(δχa) / normt;
+ δχe = sqrt(δχe) / norme;
+ δχm = sqrt(δχm) / normm;
+
+ std::ofstream out_file;
+ out_file.open(filename + ".processed");
+
+ out_file << Me << " " << δMe << "\n";
+ out_file << Mm << " " << δMm << "\n";
+ out_file << Mt << " " << δMt << "\n";
+ out_file << Ma << " " << δMa << "\n";
+ out_file << χe << " " << δχe << "\n";
+ out_file << χm << " " << δχm << "\n";
+ out_file << χt << " " << δχt << "\n";
+ out_file << χa << " " << δχa << "\n";
+
+ out_file.close();
+
+ std::cout << "Finished processing " << filename << ".\n";
+ }
+ }
+
+ return 0;
+}