u1casc-flux/obs_n.hpp

   1 #ifndef OBS_N_HPP
   2 #define OBS_N_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_n : public o815::obs {
  15
  16 public:
  17   struct obsmem {
  18     int lx4sum[2];
  19   };
  20   obs_n(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   static void preNSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
  30   static double nSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
  31
  32   sim *Sim;
  33   obstat<double,double> oN[2];
  34
  35   int lx4sum(int i);
  36 };
  37
  38 obs_n::obs_n(o815 *_O815) : o815::obs("n",
  39                                             _O815->paraQ->getParaNames() +
  40                                             "none:none_err:nonesus:nonesus_err"
  41                                             ":ntwo:ntwo_err:ntwosus:ntwosus_err",
  42                                             _O815, sizeof(obsmem) )
  43 {
  44   OM = (obsmem*)obsMem;
  45   Sim = (sim*)O815->Sim;
  46 }
  47
  48 void obs_n::_start() {
  49   //*out << "OBS_test: start" << endl;
  50 };
  51
  52 void obs_n::_meas(bool loadedobs, const int& nthmeas) {
  53   if (!loadedobs) {
  54     OM->lx4sum[0] = lx4sum(0);
  55     OM->lx4sum[1] = lx4sum(1);
  56   }
  57
  58   oN[0].addMeas(OM->lx4sum[0]);
  59   oN[1].addMeas(OM->lx4sum[1]);
  60 };
  61
  62 void obs_n::_finish() {
  63   *out << O815->paraQ->getParaVals();
  64
  65   int compid_n[2], compid_nsus[2];
  66
  67   for (int iflav=0; iflav<2; iflav++) {
  68     compid_n[iflav] = oN[iflav].computeEasy();
  69     compid_nsus[iflav] = oN[iflav].computeJack(obs_n::preNSus, obs_n::nSus, &(Sim->lsize4));
  70
  71     *out << "\t" << oN[iflav].getMean(compid_n[iflav]) / Sim->lsize4 << "\t" << oN[iflav].getErr(compid_n[iflav]) / Sim->lsize4;
  72     *out << "\t" << oN[iflav].getMean(compid_nsus[iflav]) << "\t" << oN[iflav].getErr(compid_nsus[iflav]);
  73
  74     oN[iflav].reset();
  75   }
  76
  77   *out << endl;
  78 };
  79
  80 int obs_n::lx4sum(int i)
  81 {
  82   int lx4sum = 0;
  83
  84   for( int ix = 0; ix < Sim->lsize4; ix++ )
  85     lx4sum += Sim->l[ i*Sim->lsize4*4 + ix*4 + 3 ];
  86
  87   return lx4sum;
  88 }
  89
  90 void obs_n::preNSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para)
  91 {
  92   preCalculated->push_back(0);
  93   preCalculated->push_back(0);
  94
  95   for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt)
  96     {
  97       (*preCalculated)[0] += (*valIt)[0];
  98       (*preCalculated)[1] += pow( (*valIt)[0], 2 );
  99     }
 100 }
 101
 102 double obs_n::nSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
 103   return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) ) / (nmeas-1)
 104           - pow( ( (*preCalculated)[0] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
 105     / *(int*)para;
 106 }
 107
 108 #endif