Blatt 2 Aufg 2

exc2/Data.cc

@@ -37,9 +37,9 @@ Data::Data(const std::string& filename) {
   }
 
   for (int i = 0; i < size; i++) {
-    double uncertainty;
-    file >> uncertainty;
-    m_uncertainty.push_back(uncertainty);
+    double error;
+    file >> error;
+    m_error.push_back(error);
   }
 
   // done! close the file
@@ -51,29 +51,37 @@ int Data::checkCompatibility(const Data* in, int N) const {
   int fails = 0;
 
   for (int i = 0; i < size(); i++) {
-    if (abs(in->measurement(i) - measurement(i)) > N * (in->uncertainty(i)
-          + uncertainty(i))) {
+    if (abs(in->measurement(i) - measurement(i)) > N * (in->error(i)
+          + error(i))) {
       fails++;
     }
   }
   return fails;
 }
 
+double Data::chisquare(double (*f)(double)) const {
+  double res = 0;
+  for (unsigned int i = 0; i < size(); i++) {
+    res += pow(measurement(i) - f(binCenter(i)), 2) / pow(error(i), 2);
+  }
+  return res / 52;
+}
+
 Data operator+(const Data& a, const Data& b) {
 
   assert(a.checkCompatibility(&b, 3) == 0);
   vector<double> avg_y;
   vector<double> avg_w;
   for (unsigned int i = 0; i < a.size(); i++) {
-    double w1 = 1/(a.uncertainty(i) * a.uncertainty(i));
-    double w2 = 1/(b.uncertainty(i) * b.uncertainty(i));
+    double w1 = 1/(a.error(i) * a.error(i));
+    double w2 = 1/(b.error(i) * b.error(i));
     avg_y.push_back((w1 * a.measurement(i) + w2 * b.measurement(i)) / (w1 + w2));
     avg_w.push_back(sqrt(1/(w1 + w2)));
   }
   Data d;
   d.m_bins = a.m_bins;
   d.m_data = avg_y;
-  d.m_uncertainty = avg_w;
+  d.m_error = avg_w;
   return d;
 }
 

exc2/Data.hh

@@ -11,14 +11,14 @@ class Data {
 
   unsigned int size() const { return m_data.size(); }
   double measurement(int i) const { return m_data[i]; }
-  double uncertainty(int i) const { return m_uncertainty[i]; }
-  double binCenter(int i) const { return 0; }
-  double binLow(int i) const { return 0; }
-  double binHigh(int i) const { return 0; }
-  double error(int i) const { return 0; }
+  double binCenter(int i) const { return m_bins[i] + (m_bins[i+1] - m_bins[i]) / 2.; }
+  double binLow(int i) const { return m_bins[i]; }
+  double binHigh(int i) const { return m_bins[i+1]; }
+  double error(int i) const { return m_error[i]; }
   std::string name() const { return m_name; }
 
   int checkCompatibility(const Data* in, int N) const;
+  double chisquare(double (*f)(double)) const;
 
   friend Data operator+(const Data& a, const Data& b);
 
@@ -27,7 +27,7 @@ class Data {
   void assertSizes();
   std::string m_name;
   std::vector<double> m_data;
-  std::vector<double> m_uncertainty;
+  std::vector<double> m_error;
   std::vector<double> m_bins;
 };
 

exc2/run.cc

@@ -2,6 +2,7 @@
 #include <vector>
 #include <functional>
 #include <string>
+#include <math.h>
 
 #include "Data.hh"
 
@@ -59,6 +60,10 @@ void runTests() {
     std::cout << (test() ? " ok" : " FAILED!") << std::endl;
 }
 
+double background(double x) {
+  return 0.005 - 0.00001 * x + 0.08 * exp(-0.015 * x);
+}
+
 int main() {
   using namespace std;
 
@@ -83,10 +88,11 @@ int main() {
 
   for(auto data : datas) {
     cout << "measurement of " << data.name() << " in bin 27: "
-      << data.measurement(27) << " +- " << data.uncertainty(27) << endl;
+      << data.measurement(27) << " +- " << data.error(27) << endl;
+    cout << "Chisquare of " << data.name() << ": " << data.chisquare(background) << endl;
     for (auto data2 : datas) {
       if (abs(data.measurement(27) - data2.measurement(27)) >
-          data.uncertainty(27) + data2.uncertainty(27)) {
+          data.error(27) + data2.error(27)) {
         cout << "Data point doesnt agree" << endl;
       }
       if (data.name() != data2.name()) {
@@ -98,7 +104,8 @@ int main() {
   }
   Data dat_avg = datA + datB + datC + datD;
   cout << "Combined data: " << dat_avg.measurement(27) << endl
-    << "Combined uncertainty: " << dat_avg.uncertainty(27) << endl;
+    << "Combined error: " << dat_avg.error(27) << endl
+    << "Combined chisquare: " << dat_avg.chisquare(background) << endl;
 
   return 0;
 }