u1casc-ordinary/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[2];
  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
  31   static void prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
  32   static double phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
  33
  34   sim *Sim;
  35   obstat<double,double> oPhi2[2];
  36 };
  37
  38 obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2",
  39                                             _O815->paraQ->getParaNames() +
  40                                             "phi2one:phi2one_err:phi2onesus:phi2onesus_err"
  41                                             ":phi2two:phi2two_err:phi2twosus:phi2twosus_err",
  42                                             _O815, sizeof(obsmem) ) {
  43   OM = (obsmem*)obsMem;
  44   Sim = (sim*)O815->Sim;
  45 }
  46
  47 void obs_phi2::_start() {
  48   //*out << "OBS_test: start" << endl;
  49 };
  50
  51 void obs_phi2::_meas(bool loadedobs, const int& nthmeas) {
  52   if (!loadedobs)
  53     phi2Compute();
  54
  55   for (int iflav=0; iflav<2; iflav++)
  56     oPhi2[iflav].addMeas(OM->phi2[iflav]);
  57 };
  58
  59 void obs_phi2::_finish() {
  60   *out << O815->paraQ->getParaVals();
  61
  62   int compid_phi2[2], compid_phi2sus[2];
  63
  64   for (int iflav=0; iflav<2; iflav++) {
  65     compid_phi2[iflav] = oPhi2[iflav].computeEasy();
  66     compid_phi2sus[iflav] = oPhi2[iflav].computeJack(obs_phi2::prePhi2Sus, obs_phi2::phi2Sus, &(Sim->lsize4));
  67
  68     *out << "\t" << oPhi2[iflav].getMean(compid_phi2[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2[iflav]);
  69     *out << "\t" << oPhi2[iflav].getMean(compid_phi2sus[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2sus[iflav]);
  70
  71     oPhi2[iflav].reset();
  72   }
  73
  74   *out << endl;
  75 };
  76
  77 void obs_phi2::phi2Compute()
  78 {
  79   for (int iflav=0; iflav<2; iflav++) {
  80     OM->phi2[iflav] = 0;
  81
  82     for (int ix = 0; ix < Sim->lsize4; ix++)
  83       OM->phi2[iflav] += norm(Sim->phi[ iflav*Sim->lsize4 + ix ]);
  84
  85     OM->phi2[iflav] /= Sim->lsize4;
  86   }
  87 }
  88
  89 void obs_phi2::prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para) {
  90   preCalculated->push_back(0);
  91   preCalculated->push_back(0);
  92
  93   for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
  94     (*preCalculated)[0] += (*valIt)[0];
  95     (*preCalculated)[1] += pow((*valIt)[0],2);
  96   }
  97 }
  98
  99 double obs_phi2::phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
 100   return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) ) / (nmeas-1)
 101           - pow( ( (*preCalculated)[0] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
 102     * *(int*)para;
 103 }
 104
 105 #endif