+#ifndef OBS_PHI2_HPP
+#define OBS_PHI2_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+
+using namespace std;
+
+class obs_phi2 : public o815::obs {
+
+public:
+ struct obsmem {
+ double A[2];
+ double B[2];
+ double C[2];
+ };
+ obs_phi2(o815 *_O815);
+
+private:
+ void _start();
+ void _meas(bool loadedobs, const int& nthmeas);
+ void _finish();
+
+ obsmem* OM;
+
+ void phi2Compute();
+
+ static void prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
+ static double phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
+
+ sim *Sim;
+ obstat<double,double> oPhi2[2];
+};
+
+obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2",
+ _O815->paraQ->getParaNames() +
+ "phi2one:phi2one_err:phi2onesus:phi2onesus_err"
+ ":phi2two:phi2two_err:phi2twosus:phi2twosus_err",
+ _O815, sizeof(obsmem) ) {
+ OM = (obsmem*)obsMem;
+ Sim = (sim*)O815->Sim;
+}
+
+void obs_phi2::_start() {
+ //*out << "OBS_test: start" << endl;
+};
+
+void obs_phi2::_meas(bool loadedobs, const int& nthmeas) {
+ if (!loadedobs)
+ phi2Compute();
+
+ for (int iflav=0; iflav<2; iflav++) {
+ double tmpMeas[3] = {OM->A[iflav], OM->B[iflav], OM->C[iflav]};
+ oPhi2[iflav].addMeas(tmpMeas, 3);
+ }
+};
+
+void obs_phi2::_finish() {
+ *out << O815->paraQ->getParaVals();
+
+ int compid_phi2[2], compid_phi2sus[2];
+
+ for (int iflav=0; iflav<2; iflav++) {
+ compid_phi2[iflav] = oPhi2[iflav].computeEasy();
+ compid_phi2sus[iflav] = oPhi2[iflav].computeJack(obs_phi2::prePhi2Sus, obs_phi2::phi2Sus, &(Sim->lsize4));
+
+ *out << "\t" << oPhi2[iflav].getMean(compid_phi2[iflav]) / Sim->lsize4 << "\t" << oPhi2[iflav].getErr(compid_phi2[iflav]) / Sim->lsize4;
+ *out << "\t" << oPhi2[iflav].getMean(compid_phi2sus[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2sus[iflav]);
+
+ oPhi2[iflav].reset();
+ }
+
+ *out << endl;
+};
+
+void obs_phi2::phi2Compute()
+{
+ for(int i=0; i<2; i++) {
+ OM->A[i] = 0;
+ OM->B[i] = 0;
+ OM->C[i] = 0;
+
+ for (int x = 0; x < Sim->lsize4; x++) {
+ double wp2w;
+
+ wp2w = Sim->WF( Sim->wArg(x,i) + 2, Sim->wArg(x,i), i );
+
+ OM->A[i] += wp2w;
+ OM->B[i] += Sim->WF( Sim->wArg(x,i) + 4, Sim->wArg(x,i), i );
+ OM->C[i] += pow(wp2w, 2);
+ }
+ }
+}
+
+void obs_phi2::prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para)
+{
+ preCalculated->push_back(0);
+ preCalculated->push_back(0);
+
+ for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt)
+ {
+ (*preCalculated)[0] += pow( (*valIt)[0], 2 ) + (*valIt)[1] - (*valIt)[2];
+ (*preCalculated)[1] += (*valIt)[0];
+ }
+}
+
+double obs_phi2::phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
+ return (( (*preCalculated)[0] - pow((*excludedMeas)[0], 2) - (*excludedMeas)[1] + (*excludedMeas)[2] ) / (nmeas-1)
+ - pow( ( (*preCalculated)[1] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
+ / *(int*)para;
+}
+
+#endif