u1casc-ordinary/obs_corrphiphimass.hpp

   1 #ifndef OBS_CORRPHIPHIMASS_HPP
   2 #define OBS_CORRPHIPHIMASS_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 #include <complex>
  12
  13 #include <cmath>
  14
  15 #define U1 0
  16 #define U2 1
  17 #define F1 2
  18 #define F2 3
  19
  20 using namespace std;
  21
  22 class obs_corrphiphimass : public o815::obs {
  23
  24 public:
  25   obs_corrphiphimass(o815 *_O815);
  26
  27 private:
  28   void _start();
  29   void _meas(bool loadedobs, const int& nthmeas);
  30   void _finish();
  31
  32   void corrCompute();
  33   static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
  34   static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
  35
  36   sim *Sim;
  37   obstat< complex<double>, complex<double>  > oC;
  38
  39   int spatialV;
  40
  41   complex<double> *OM[4];
  42 };
  43
  44 obs_corrphiphimass::obs_corrphiphimass(o815 *_O815) : o815::obs("corr",
  45                                                                   _O815->paraQ->getParaNames() +
  46                                                                   "tsep"
  47                                                                   ":u1_effmass:u1_effmass_err",
  48                                                                   _O815, 4*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2+1) ) {
  49
  50   for (int ivar = 0; ivar<4; ivar++)
  51     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 + 1 ) );
  52
  53   Sim = (sim*)O815->Sim;
  54   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
  55 }
  56
  57 void obs_corrphiphimass::_start() {
  58   //*out << "OBS_test: start" << endl;
  59 };
  60
  61 void obs_corrphiphimass::_meas(bool loadedobs, const int& nthmeas) {
  62   if (!loadedobs)
  63     corrCompute();
  64
  65   oC.addMeas( OM[U1], O815->comargs.lsize[1]/2+1 );
  66 };
  67
  68 void obs_corrphiphimass::_finish() {
  69   int compid_corr[O815->comargs.lsize[1]/2][4];
  70   int compid_effmass[O815->comargs.lsize[1]/2-1][4];
  71
  72   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++)
  73     compid_corr[itsep][U1] = oC.computeEasy(itsep);
  74
  75   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
  76     pair<int,int> effmasspass( itsep, O815->comargs.lsize[1]/2 );
  77     compid_effmass[itsep][U1] = oC.computeJack(obs_corrphiphimass::preEffMass, obs_corrphiphimass::effMass, &effmasspass);
  78   }
  79
  80   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
  81     *out << O815->paraQ->getParaVals();
  82     *out << "\t" << itsep;
  83
  84     *out << "\t" << real( oC.getMean(compid_effmass[itsep][U1]) )
  85          << "\t" << real( oC.getErr (compid_effmass[itsep][U1]) );
  86
  87     *out << endl;
  88   }
  89
  90   oC.reset();
  91 };
  92
  93 void obs_corrphiphimass::corrCompute()
  94 {
  95   complex<double> phislice[4][O815->comargs.lsize[1]];
  96
  97   OM[U1][O815->comargs.lsize[1]/2] = 0;
  98   OM[U2][O815->comargs.lsize[1]/2] = 0;
  99   OM[F1][O815->comargs.lsize[1]/2] = 0;
 100   OM[F2][O815->comargs.lsize[1]/2] = 0;
 101
 102   for (int it = 0; it < O815->comargs.lsize[1]; it++) {
 103     phislice[U1][it] = 0;
 104     phislice[U2][it] = 0;
 105     phislice[F1][it] = 0;
 106     phislice[F2][it] = 0;
 107
 108     for (int ix = 0; ix < spatialV; ix++) {
 109       phislice[U1][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 110       phislice[U2][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 111       phislice[F1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 112       phislice[F2][it] += Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 113     }
 114     phislice[U1][it] /= spatialV;
 115     phislice[U2][it] /= spatialV;
 116     phislice[F1][it] /= spatialV;
 117     phislice[F2][it] /= spatialV;
 118
 119     OM[U1][O815->comargs.lsize[1]/2] += phislice[U1][it];
 120     OM[U2][O815->comargs.lsize[1]/2] += phislice[U2][it];
 121     OM[F1][O815->comargs.lsize[1]/2] += phislice[F1][it];
 122     OM[F2][O815->comargs.lsize[1]/2] += phislice[F2][it];
 123   }
 124
 125   for (int icorr = 0; icorr < 4; icorr++) {
 126     for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 127       OM[icorr][itsep] = 0;
 128
 129       for (int it = 0; it < O815->comargs.lsize[1]; it++)
 130         OM[icorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[icorr][it] );
 131
 132       OM[icorr][itsep] /= O815->comargs.lsize[1];
 133     }
 134
 135     OM[icorr][O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
 136   }
 137 }
 138
 139 void obs_corrphiphimass::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
 140   pair<int,int> *myparas = (pair<int,int>*)para;
 141
 142   preCalculated->push_back(0);
 143   preCalculated->push_back(0);
 144   preCalculated->push_back(0);
 145
 146   for(vector< vector< complex<double> > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
 147     (*preCalculated)[0] += (*valIt)[myparas->first];
 148     (*preCalculated)[1] += (*valIt)[myparas->first+1];
 149     (*preCalculated)[2] += (*valIt)[myparas->second];
 150   }
 151 }
 152
 153 complex<double> obs_corrphiphimass::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
 154   pair<int,int> *myparas = (pair<int,int>*)para;
 155
 156   double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
 157
 158   return std::log( abs(
 159                        ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) /
 160                        ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected )
 161                         ) );
 162 }
 163
 164 #endif