obs_phi2.hpp

   1 #ifndef OBS_PHI2_HPP
   2 #define OBS_PHI2_HPP
   3
   4 #include "latlib/o815/o815.h"
   5
   6 #include "latlib/writeout.h"
   7
   8 #include "latlib/obstat.hpp"
   9
  10 #include <iostream>
  11
  12 using namespace std;
  13
  14 class obs_phi2 : public o815::obs {
  15
  16 public:
  17   struct obsmem {
  18     double phi2;
  19   };
  20   obs_phi2(o815 *_O815);
  21
  22 private:
  23   void _start();
  24   void _meas(bool loadedobs, const int& nthmeas);
  25   void _finish();
  26
  27   obsmem* OM;
  28
  29   void phi2Compute();
  30   static double phi2Sus(vector< vector<double> > *vals, void *para);
  31
  32   sim *Sim;
  33   obstat<double,double> oPhi2;
  34 };
  35
  36 obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2",
  37                                             _O815->paraQ->getParaNames() + "phi2:phi2_err:phi2sus:phi2sus_err",
  38                                             _O815, sizeof(obsmem) ) {
  39   OM = (obsmem*)obsMem;
  40   Sim = (sim*)O815->Sim;
  41 }
  42
  43 void obs_phi2::_start() {
  44   //*out << "OBS_test: start" << endl;
  45 };
  46
  47 void obs_phi2::_meas(bool loadedobs, const int& nthmeas) {
  48   if (!loadedobs)
  49     phi2Compute();
  50
  51   oPhi2.addMeas(OM->phi2);
  52 };
  53
  54 void obs_phi2::_finish() {
  55   *out << O815->paraQ->getParaVals();
  56
  57   int compid_phi2, compid_phi2sus;
  58
  59   compid_phi2 = oPhi2.computeEasy();
  60   compid_phi2sus = oPhi2.computeJack(obs_phi2::phi2Sus, &(Sim->LSIZE2));
  61
  62   *out << "\t" << oPhi2.getMean(compid_phi2) << "\t" << oPhi2.getErr(compid_phi2);
  63   *out << "\t" << oPhi2.getMean(compid_phi2sus) << "\t" << oPhi2.getErr(compid_phi2sus);
  64
  65   *out << endl;
  66
  67   oPhi2.reset();
  68 };
  69
  70 void obs_phi2::phi2Compute()
  71 {
  72   OM->phi2 = 0;
  73
  74   for (int ix = 0; ix < Sim->LSIZE2; ix++)
  75     OM->phi2 += norm( Sim->conf[ix].phi );
  76
  77   OM->phi2 /= Sim->LSIZE2;
  78 }
  79
  80 double obs_phi2::phi2Sus(vector< vector<double> > *vals, void *para) {
  81   double mean=0, mean2=0;
  82   for(vector< vector<double> >::iterator valIt = vals->begin(); valIt != vals->end(); ++valIt)
  83     {
  84       mean += (*valIt)[0];
  85       mean2 += pow((*valIt)[0],2);
  86     }
  87   mean /= vals->size();
  88   mean2 /= vals->size();
  89
  90   return ( mean2 - pow(mean,2) ) * *(int*)para;
  91 }
  92
  93 #endif