Page MenuHomeHEPForge

buildInterpolatedCheck.cc
No OneTemporary

Size
3 KB
Referenced Files
None
Subscribers
None

buildInterpolatedCheck.cc

#include <cstdlib>
#include <iostream>
#include "test_utils.hh"
#include "geners/CPP11_array.hh"
#include "geners/BinaryFileArchive.hh"
#include "geners/Record.hh"
#include "npstat/rng/MersenneTwister.hh"
#include "npstat/nm/LinearMapper1d.hh"
#include "npstat/stat/Distributions1D.hh"
#include "npstat/stat/buildInterpolatedDistro1DNP.hh"
#include "npstat/stat/NonparametricDistro1DBuilder.hh"
#include "npstat/stat/LOrPEMarginalSmoother.hh"
#include "npstat/stat/RatioResponseIntervalBuilder.hh"
#include "npstat/stat/distro1DStats.hh"
#include "npstat/stat/InMemoryNtuple.hh"
using namespace npstat;
using namespace std;
int main(int argc, char const* argv[])
{
typedef CPP11_array<float,3> Point;
if (argc != 5)
{
cout << "\nUsage: " << parse_progname(argv[0])
<< " npoints nEffective nGridPoints archive\n" << endl;
exit(1);
}
const unsigned pointDim = PointDimensionality<Point>::dim_size;
unsigned npoints, nGridPoints;
double nEffective;
if (parse_unsigned(argv[1], &npoints)) exit(1);
if (parse_double(argv[2], &nEffective)) exit(1);
if (parse_unsigned(argv[3], &nGridPoints)) exit(1);
const char* archName = argv[4];
MersenneTwister rng;
InMemoryNtuple<float> nt(ntupleColumns("r", "pt", "sigma"), "Points");
assert(nt.nColumns() == pointDim);
std::vector<Point> data;
data.reserve(npoints);
Point pt;
LinearMapper1d meanmap(10.0, 0.5, 110.0, 1.0);
const double ptcut = 10.0;
for (unsigned i=0; i<npoints; )
{
const double ppt = 10.0 + rng()*100.0;
const double rmean = meanmap(ppt);
const double sigma = 0.1 + 0.1*rng();
Gauss1D g(rmean, sigma);
double random;
g.random(rng, &random);
const double ptmeas = ppt*random;
if (ptmeas > ptcut)
{
pt[0] = random;
pt[1] = ppt;
pt[2] = sigma;
data.push_back(pt);
nt.fill(&pt[0], pointDim);
++i;
}
}
unsigned dimPredictors[2] = {1U, 2U};
unsigned predictorNumBins[2];
predictorNumBins[0] = nGridPoints;
predictorNumBins[1] = nGridPoints;
const int predictorSymbetaPower = -1;
const unsigned responseDimToUse = 0U;
BoundaryHandling bm("BM_FOLD");
LOrPEMarginalSmoother smoother(1000, 0.0, 2.0, -1, 0.0, 1.0, bm, "r");
RatioResponseIntervalBuilder<Point> dib(1U, ptcut, 0.01);
NonparametricDistro1DBuilder<Point> builder(&smoother, &dib, false);
gs::BinaryFileArchive ar(archName, "w");
ar << gs::Record(nt, "Points", "");
builder.setArchive(&ar, "Histos");
CPP11_auto_ptr<InterpolatedDistro1DNP> idro = buildInterpolatedDistro1DNP(
data, dimPredictors, 2U, 0, predictorNumBins, predictorSymbetaPower,
nEffective, false, responseDimToUse, builder);
assert(idro->nAxes() == 2U);
assert(idro->getAxis(0).nCoords() == predictorNumBins[0]);
assert(idro->getAxis(1).nCoords() == predictorNumBins[1]);
unsigned cell[2];
unsigned count = 0;
for (unsigned i0=0; i0<predictorNumBins[0]; ++i0)
{
cell[0] = i0;
const double ppt = idro->getAxis(0).coordinate(i0);
const double mean = meanmap(ppt);
for (unsigned i1=0; i1<predictorNumBins[1]; ++i1, ++count)
{
cell[1] = i1;
const double sigma = idro->getAxis(1).coordinate(i1);
const AbsDistribution1D* distro = idro->getGridDistro(cell, 2U);
double fitMean, fitSigma;
distro1DStats(*distro, 0.0, 2.0, 2000, &fitMean, &fitSigma);
cout << count << ' ' << i0 << ' ' << i1 << ' ' << ppt << " "
<< fitMean << " (" << mean << ") "
<< fitSigma << " (" << sigma << ")"
<< endl;
}
}
return 0;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:48 AM (7 h, 14 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6558527
Default Alt Text
buildInterpolatedCheck.cc (3 KB)

Event Timeline