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