u1casc-ordinary/obs_corrtetra.hpp

   1 #ifndef OBS_CORRTETRA_HPP
   2 #define OBS_CORRTETRA_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_corrtetra : public o815::obs {
  18
  19 public:
  20   obs_corrtetra(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_corrtetra::obs_corrtetra(o815 *_O815) : o815::obs("corrtetra",
  40                                                       _O815->paraQ->getParaNames() +
  41                                                       "tsep"
  42                                                       ":tetra_real:tetra_imag:tetra_abs:tetra_mass:tetra_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_corrtetra::_start() {
  52   //*out << "OBS_test: start" << endl;
  53 };
  54
  55 void obs_corrtetra::_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_corrtetra::_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_corrtetra::preEffMass, obs_corrtetra::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_corrtetra::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[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]
 106         * conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 107
 108     phislice[it] /= sqrt(2)*spatialV;
 109
 110     OM[O815->comargs.lsize[1]/2] += phislice[it];
 111   }
 112
 113   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 114     OM[itsep] = 0;
 115
 116     for (int it = 0; it < O815->comargs.lsize[1]; it++)
 117       OM[itsep] += phislice[ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[it] );
 118
 119     OM[itsep] /= O815->comargs.lsize[1];
 120   }
 121
 122   OM[O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
 123 }
 124
 125 void obs_corrtetra::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
 126   pair<int,int> *myparas = (pair<int,int>*)para;
 127
 128   preCalculated->push_back(0);
 129   preCalculated->push_back(0);
 130   preCalculated->push_back(0);
 131
 132   for(vector< vector< complex<double> > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
 133     (*preCalculated)[0] += (*valIt)[myparas->first];
 134     (*preCalculated)[1] += (*valIt)[myparas->first+1];
 135     (*preCalculated)[2] += (*valIt)[myparas->second];
 136   }
 137 }
 138
 139 complex<double> obs_corrtetra::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
 140   pair<int,int> *myparas = (pair<int,int>*)para;
 141
 142   double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
 143
 144   return std::log( abs(
 145                        ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) /
 146                        ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected )
 147                         ) );
 148 }
 149
 150 #endif