portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2024-11-12 17:35:35 +00:00
Code Releases Activity

Now includes functioning const fit; tidying up now

Browse Source
This commit is contained in:
Andrew Zhen Ning Yong 2019-01-24 12:20:13 +00:00
parent 741762f17b
commit 82d41d0b5d
1 changed files with 93 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

140
physics/2pt-fit.cpp
View File

@ -117,7 +117,7 @@ int main(int argc, char *argv[])
// make model ////////////////////////////////////////////////////////////// // make model //////////////////////////////////////////////////////////////
DoubleModel mod; DoubleModel mod;
bool coshModel = false, linearModel = false; bool coshModel = false, linearModel = false, constModel = false;
if ((model == "exp") or (model == "exp1")) if ((model == "exp") or (model == "exp1"))
{ {
@ -185,11 +185,12 @@ int main(int argc, char *argv[])
} }
else if (model == "const") else if (model == "const")
{ {
nPar = 1; constModel = true;
mod.setFunction([](const double *x, const double *p) nPar = 1;
mod.setFunction([](const double *x __dumb, const double *p)
{ {
return p[0]; return p[0];
}, 0, nPar); }, 1, nPar);
} }
else else
{ {
@ -222,21 +223,35 @@ int main(int argc, char *argv[])
data.addXDim(nt, "t/a", true); data.addXDim(nt, "t/a", true);
data.addYDim("C(t)"); data.addYDim("C(t)");
data.setUnidimData(tvec, corr); data.setUnidimData(tvec, corr);
for (Index p = 0; p < nPar; p += 2) // set parameter name /////////////
if(constModel)
{ {
mod.parName().setName(p, "E_" + strFrom(p/2)); mod.parName().setName(0, "delta_m");
mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
} }
else
{
for (Index p = 0; p < nPar; p += 2)
{
mod.parName().setName(p, "E_" + strFrom(p/2));
mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
}
}
//set initial values ////////////////
if (linearModel) if (linearModel)
{ {
init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central]; init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central];
init(1) = data.y(nt/4, 0)[central] + nt/4*init(0); init(1) = data.y(nt/4, 0)[central] + nt/4*init(0);
} }
else if(constModel)
{
init(0) = data.y(nt/4, 0)[central];
}
else else
{ {
init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]); init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]);
init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4)); init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4));
} }
for (Index p = 2; p < nPar; p += 2) for (Index p = 2; p < nPar; p += 2)
{ {
init(p) = 2*init(p - 2); init(p) = 2*init(p - 2);
@ -249,14 +264,24 @@ int main(int argc, char *argv[])
globMin.setLowLimit(p, -10.*fabs(init(p))); globMin.setLowLimit(p, -10.*fabs(init(p)));
globMin.setHighLimit(p, 10.*fabs(init(p))); globMin.setHighLimit(p, 10.*fabs(init(p)));
} }
else if(constModel)
{
globMin.setLowLimit(p, -1.1*fabs(init(0)));
locMin.setLowLimit(p, 0.);
globMin.setHighLimit(p, 1.1*fabs(init(0)));
}
else else
{ {
globMin.setLowLimit(p, 0.); globMin.setLowLimit(p, 0.);
locMin.setLowLimit(p, 0.); locMin.setLowLimit(p, 0.);
globMin.setHighLimit(p, 10.*init(p)); globMin.setHighLimit(p, 10.*init(p));
} }
globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1))); if(!constModel)
globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1))); {
globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1)));
globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1)));
}
} }
globMin.setPrecision(0.001); globMin.setPrecision(0.001);
globMin.setMaxIteration(100000); globMin.setMaxIteration(100000);
@ -286,69 +311,89 @@ int main(int argc, char *argv[])
fit = data.fit(locMin, init, mod); fit = data.fit(locMin, init, mod);
fit.print(); fit.print();
} }
// plots /////////////////////////////////////////////////////////////////// // plots ///////////////////////////////////////////////////////////////////
if (doPlot) if (doPlot)
{ {
Plot p; Plot p;
DMatSample effMass(nSample);
DVec effMassT, fitErr;
Index maxT = (coshModel) ? (nt - 2) : (nt - 1);
double e0, e0Err;
// if (!constModel)
// {
// p << PlotRange(Axis::x, 0, nt - 1);
// if (!linearModel)
// {
// p << LogScale(Axis::y);
// }
// p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1); //<<<< problematic line for const fit
// p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1);
// p << Color("rgb 'red'") << PlotData(data.getData());
// p.display();
// }
p << PlotRange(Axis::x, 0, nt - 1); p << PlotRange(Axis::x, 0, nt - 1);
if (!linearModel) if (!linearModel and !constModel)
{ {
p << LogScale(Axis::y); p << LogScale(Axis::y);
} }
p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1); p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1);
p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1); p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1); //<<<< problematic line for const fit
p << Color("rgb 'red'") << PlotData(data.getData()); p << Color("rgb 'red'") << PlotData(data.getData());
p.display(); p.display();
effMass.resizeMat(maxT, 1);
effMassT.setLinSpaced(maxT, 1, maxT); // effective mass plot //////////////////////////////////////////////////////
fitErr = fit.variance().cwiseSqrt(); if (!constModel)
e0 = fit[central](0);
e0Err = fitErr(0);
if (coshModel)
{ {
FOR_STAT_ARRAY(effMass, s) DMatSample effMass(nSample);
DVec effMassT, fitErr;
Index maxT = (coshModel) ? (nt - 2) : (nt - 1);
double e0, e0Err;
effMass.resizeMat(maxT, 1);
effMassT.setLinSpaced(maxT, 0, maxT-1);
fitErr = fit.variance().cwiseSqrt();
e0 = fit[central](0);
e0Err = fitErr(0);
if (coshModel)
{ {
for (Index t = 1; t < nt - 1; ++t) FOR_STAT_ARRAY(effMass, s)
{ {
effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1)) for (Index t = 1; t < nt - 1; ++t)
/(2.*corr[s](t))); {
effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1))
/(2.*corr[s](t)));
}
} }
} }
} else if (linearModel)
else if (linearModel)
{
FOR_STAT_ARRAY(effMass, s)
{ {
for (Index t = 0; t < nt - 1; ++t) FOR_STAT_ARRAY(effMass, s)
{ {
effMass[s](t) = corr[s](t) - corr[s](t+1); for (Index t = 0; t < nt - 1; ++t)
{
effMass[s](t) = corr[s](t) - corr[s](t+1);
}
} }
} }
} else
else
{
FOR_STAT_ARRAY(effMass, s)
{ {
for (Index t = 1; t < nt; ++t) FOR_STAT_ARRAY(effMass, s)
{ {
effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t)); for (Index t = 1; t < nt; ++t)
{
effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t));
}
} }
} }
p.reset();
p << PlotRange(Axis::x, 0, maxT);
p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
p << Color("rgb 'blue'") << PlotHLine(e0);
p << Color("rgb 'red'") << PlotData(effMassT, effMass);
p << Caption("Effective Mass");
p.display();
} }
p.reset(); }
p << PlotRange(Axis::x, 1, maxT);
p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
p << Color("rgb 'blue'") << PlotHLine(e0);
p << Color("rgb 'red'") << PlotData(effMassT, effMass);
p.display();
}
if (doHeatmap) if (doHeatmap)
{ {
Plot p; Plot p;
@ -367,6 +412,7 @@ int main(int argc, char *argv[])
} }
} }
// output ////////////////////////////////////////////////////////////////// // output //////////////////////////////////////////////////////////////////
if (!outFileName.empty()) if (!outFileName.empty())
{ {