u1casc-ordinary/obs_corrchichi.hpp

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