Table.cpp
No OneTemporary
Actions

Size

26 KB

Subscribers

None

Table.cpp
View Options

	//==============================================================================
	// Table.cpp
	//
	// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
	//
	// This file is part of Sartre version: 1.1
	//
	// This program is free software: you can redistribute it and/or modify
	// it under the terms of the GNU General Public License as published by
	// the Free Software Foundation.
	// This program is distributed in the hope that it will be useful,
	// but without any warranty; without even the implied warranty of
	// merchantability or fitness for a particular purpose. See the
	// GNU General Public License for more details.
	// You should have received a copy of the GNU General Public License
	// along with this program. If not, see <http://www.gnu.org/licenses/>.
	//
	// Author: Thomas Ullrich
	// Last update:
	// $Date: 2013-05-29 21:26:19 +0100 (Wed, 29 May 2013) $
	// $Author: thomas.ullrich@bnl.gov $
	//==============================================================================
	//
	// Table data is stored in a TH3F.
	//
	// All information is stored in the histogram title.
	//
	// Table ID is encoded as follows:
	//
	// histogram title is a number that is to be interpreted
	// as an uint64_t with the bits set as follows:
	//
	// bit 0: content type: 0 for <A>, 1 for <A2>
	// bit 1: polarization: T=1, L=0
	// bit 2: t encoding: 0 for \|t\|, 1 for log(\|t\|)
	// bit 3: W2 encoding: 0 for lin, 1 for log
	// bit 4: Q2 encoding: 0 for lin, 1 for log
	// bit 5-7: dipole model type (Enumerations.h)
	// bit 8-15: mass number A
	// bit 16-31: vector meson ID (PDG)
	// bit 32: content encoding: 0 in lin, 1 in log
	// bit 33: content type is lambda_<A> (bit 0 == 0 in this case)
	// bit 34-63: not used
	//
	// log implies ln
	// bit 0 is the LSB
	//
	// Internal histogram: x = Q2, y = W2, z=t (or the logs)
	//
	// Actual value is taken at the center of the bin.
	// min/max functions (e.g. minQ2()) return the value of the first/last
	// bin center.
	//==============================================================================
	#include "Table.h"
	#include "TFile.h"
	#include "TH3F.h"
	#include <cstdio>
	#include <cmath>
	#include <ctime>
	#include <algorithm>
	#include <string>
	#include <sstream>
	#include <fstream>
	#include <iomanip>
	#include <limits>
	#include <unistd.h>
	#include <sys/types.h>
	#include "GridSpline.h"

	#define PR(x) cout << #x << " = " << (x) << endl;

	Table::Table()
	{
	mID = 0;
	m3DHisto = 0;
	mFillCounter = 0;
	mBackupFrequence = 0;
	}

	Table& Table::operator=(const Table& tab)
	{
	if (this != &tab) {
	delete m3DHisto;

	mID = tab.mID;
	m3DHisto = new TH3F(*(tab.m3DHisto));
	m3DHisto->SetDirectory(0);
	mFilename = tab.mID;;
	mFillCounter = tab.mID;;
	mBackupFrequence = tab.mID;;
	mBackupPrefix = tab.mID;;
	mLastBackupFilename = tab.mID;;
	}
	return *this;
	}

	Table::Table(const Table& tab)
	{
	mID = tab.mID;
	m3DHisto = new TH3F(*(tab.m3DHisto));
	m3DHisto->SetDirectory(0);
	mFilename = tab.mID;;
	mFillCounter = tab.mID;;
	mBackupFrequence = tab.mID;;
	mBackupPrefix = tab.mID;;
	mLastBackupFilename = tab.mID;;
	}

	Table::~Table()
	{
	delete m3DHisto;
	}

	unsigned int Table::create(int nbinsQ2, double Q2min, double Q2max,
	int nbinsW2, double W2min, double W2max,
	int nbinsT, double tmin, double tmax,
	bool logQ2, bool logW2, bool logt, bool logC,
	AmplitudeMoment mom, GammaPolarization pol,
	unsigned int A, int vm,
	DipoleModelType model, const char* filename)
	{
	if (m3DHisto != 0) {
	cout << "Table:create(): cannot create, table already exists." << endl;
	return 0;
	}

	if (filename) mFilename = string(filename); // name of file where table is written to

	//
	// Create table ID first
	//
	mID = (vm << 16);
	mID \|= (A << 8);
	mID \|= (model << 5);
	if (mom == mean_A2) mID \|= 1;
	if (pol == transverse) mID \|= (1 << 1);
	if (logt) mID \|= (1 << 2);
	if (logW2) mID \|= (1 << 3);
	if (logQ2) mID \|= (1 << 4);
	if (logC) mID \|= (static_cast<uint64_t>(1) << 32);
	if (mom == lambda_A) mID \|= (static_cast<uint64_t>(1) << 33);
	ostringstream titlestream;
	titlestream << mID;
	string title = titlestream.str();

	//
	// Binning
	//
	// Interpolate() only works up to the bin center
	// of the first and last bin. To guarantee that
	// the full range is accessible we shift the min
	// and max of each axis.
	//
	tmin = fabs(tmin);
	tmax = fabs(tmax);
	if (tmin > tmax) swap(tmin, tmax);
	if (logQ2) Q2min = log(Q2min);
	if (logQ2) Q2max = log(Q2max);
	if (logW2) W2min = log(W2min);
	if (logW2) W2max = log(W2max);
	if (logt) tmin = log(tmin);
	if (logt) tmax = log(tmax);
	double binWidthQ2 = (Q2max - Q2min)/(nbinsQ2-1);
	double binWidthW2 = (W2max - W2min)/(nbinsW2-1);
	double binWidthT = (tmax - tmin)/(nbinsT-1);

	//
	// Book 3D histo to hold table
	//
	m3DHisto = new TH3F("table", title.c_str(),
	nbinsQ2, Q2min-binWidthQ2/2., Q2max+binWidthQ2/2.,
	nbinsW2, W2min-binWidthW2/2., W2max+binWidthW2/2.,
	nbinsT, tmin-binWidthT/2., tmax+binWidthT/2.);
	m3DHisto->SetDirectory(0);

	return nbinsQ2nbinsW2nbinsT; // return total number of bins
	}

	int Table::numberOfEntries() const
	{
	if (m3DHisto) {
	int nx = m3DHisto->GetXaxis()->GetNbins();
	int ny = m3DHisto->GetYaxis()->GetNbins();
	int nz = m3DHisto->GetZaxis()->GetNbins();
	return nxnynz;
	}
	else
	return 0;
	}

	void Table::binXYZ(int globalBin, int& binx, int& biny, int& binz) const
	{
	//
	// Find correct bins for each axis.
	// Don't use ROOT global bins here,
	// they are a mess since they cross
	// over underflow and overflow bins.
	//
	// All bins returned count starting
	// at 1. The global bin starts at
	// 0 as usual in C++.
	//
	int nx = m3DHisto->GetXaxis()->GetNbins();
	int ny = m3DHisto->GetYaxis()->GetNbins();

	binx = globalBin%nx;
	biny = ((globalBin-binx)/nx)%ny;
	binz = ((globalBin-binx)/nx -biny)/ny;
	binx++; biny++; binz++;
	}

	void Table::binCenter(int i, double& Q2, double& W2, double& t) const
	{
	int binx, biny, binz;
	binXYZ(i, binx, biny, binz);
	// cout << i << '\t' << binx << '\t' << biny << '\t' << binz << endl;

	double x = m3DHisto->GetXaxis()->GetBinCenter(binx);
	double y = m3DHisto->GetYaxis()->GetBinCenter(biny);
	double z = m3DHisto->GetZaxis()->GetBinCenter(binz);
	Q2 = isLogQ2() ? exp(x) : x;
	W2 = isLogW2() ? exp(y) : y;
	t = isLogT() ? -exp(z) : -z;
	if (t > 0) t = 0; // avoid numeric glitch
	}

	void Table::fill(int i, double val, double err)
	{
	int binx, biny, binz;
	binXYZ(i, binx, biny, binz);

	if (isLogContent()) {
	if (val == 0) {
	cout << "Table::fill(): warning, log scale requested but value = 0." << endl;
	val = numeric_limits<float>::min();
	}
	val = log(fabs(val));
	}
	m3DHisto->SetBinContent(binx, biny, binz, val);
	if(err!=0.)
	m3DHisto->SetBinError(binx, biny, binz, err);
	mFillCounter++;

	//
	// Check if backup is due
	//
	if (mBackupFrequence)
	if (mFillCounter%mBackupFrequence == 0) backup(i);
	}

	bool Table::fexists(const char* filename) const
	{
	ifstream ifs(filename);
	return ifs;
	}

	void Table::write(const char* filename)
	{
	//
	// 'filename' is optional argument. Null value implies
	// that one was already passed via create(). Check
	// that this is the case. If one is given, it is used
	// no matter if it already was defined or not.
	//
	if (filename)
	mFilename = string(filename);
	else {
	if (mFilename.empty()) { // should be defined but is not
	cout << "Table::write(): Warning, no filename supplied. Will use 'sartre_table.root'." << endl;
	mFilename = string("sartre_table.root");
	}
	}

	//
	// Check if file exist. If so alter
	// the filename. Creating tables is
	// a time consuming business so we do
	// not want to lose precious data if
	// something goes wrong here and we do
	// want to prevent accidents.
	//
	while (fexists(mFilename.c_str())) {
	mFilename += string(".save");
	}

	TFile hfile(mFilename.c_str(),"RECREATE");
	m3DHisto->Write();
	hfile.Close();

	cout << "Table::write(): table stored in file '" << mFilename.c_str()<< "'." << endl;

	if (mBackupFrequence && !mLastBackupFilename.empty())
	remove(mLastBackupFilename.c_str());
	}

	bool Table::read(const char* filename)
	{
	//
	// Read histogram into memory
	//
	TFile *file = TFile::Open(filename,"READ");
	if (file == 0) {
	cout << "Table::read(): failed opening file '" << filename << "'." << endl;
	return false;
	}

	m3DHisto = reinterpret_cast<TH3F*>(file->Get("table"));
	if (m3DHisto == 0) {
	cout << "Table::read(): failed retrieving table from file '" << filename << "'." << endl;
	return false;
	}
	m3DHisto->SetDirectory(0);
	file->Close();

	//
	// Get table ID
	//
	mID = atoll(m3DHisto->GetTitle());

	mFilename = string(filename);
	return true;
	}

	int Table::vectorMesonId() const {return ((mID >> 16) & 0xFFFF);}

	unsigned int Table::dipoleModelType() const {return ((mID >> 5) & 0x7);}

	unsigned int Table::A() const {return ((mID >> 8) & 0xFF);}

	bool Table::isLongitudinal() const {return !isTransverse();}

	bool Table::isTransverse() const {return (mID & (1 << 1));}

	GammaPolarization Table::polarization() const {return (isTransverse() ? transverse : longitudinal);}

	AmplitudeMoment Table::moment() const {
	if (isLambdaA())
	return lambda_A;
	else if (isMeanA())
	return mean_A;
	else
	return mean_A2;
	}

	bool Table::isMeanA() const {return !(mID & 1) && !(mID & (static_cast<uint64_t>(1) << 33));}

	bool Table::isMeanA2() const {return (mID & 1) && !(mID & (static_cast<uint64_t>(1) << 33));}

	bool Table::isLogQ2() const {return (mID & (1 << 4));}

	bool Table::isLogW2() const {return (mID & (1 << 3));}

	bool Table::isLogT() const {return (mID & (1 << 2));}

	bool Table::isLogContent() const {return (mID & (static_cast<uint64_t>(1) << 32));}

	bool Table::isLambdaA() const {return !(mID & 1) && (mID & (static_cast<uint64_t>(1) << 33));}

	uint64_t Table::id() const {return mID;}

	double Table::binWidthW2() const
	{
	return m3DHisto->GetYaxis()->GetBinWidth(1);
	}

	double Table::get(double Q2, double W2, double t) const
	{
	if (m3DHisto == 0) return 0;

	//
	// Transform variables to how they are stored in the table
	//
	double x = isLogQ2() ? log(Q2) : Q2;
	double y = isLogW2() ? log(W2) : W2;
	t = fabs(t);
	double z = isLogT() ? log(t) : t;

	//
	// Tiny numerical glitches will cause TH3F::Interpolate() to fail
	// since it requires that the variables lie between the first
	// and last bin center, excluding the centers.
	// Here we enforce that this is the case. The downside of this
	// is that all values of Q2, W2, t will be forced to lie within.
	// Hence the user has to make sure that the values are within
	// the boundaries (see minQ2(), maxQ2(), minW2() etc.) before
	// calling get().
	// In this case the corrections below are tiny and are only
	// applied to avoid minor glitches that do not affect the results.
	//

	TAxis *axis = m3DHisto->GetXaxis();
	x = max(x, axis->GetBinCenter(1)+numeric_limits<float>::epsilon());
	x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::epsilon());
	axis = m3DHisto->GetYaxis();
	y = max(y, axis->GetBinCenter(1)+numeric_limits<float>::epsilon());
	y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::epsilon());
	axis = m3DHisto->GetZaxis();
	z = max(z, axis->GetBinCenter(1)+numeric_limits<float>::epsilon());
	z = min(z, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::epsilon());

	// double result = InterpolateGridSpline(x, y, z); // tmp uncommented until 0's in tables are cleared
	double result = m3DHisto->Interpolate(x, y, z);

	if (result == 0 && isLogContent()) {
	cout << "Table::get(): warning, 0 is not a valid table content when working in log scale." << endl;
	}
	if (isLogContent()) result = exp(result);
	return result;
	}

	double Table::minQ2() const
	{
	if (m3DHisto == 0) return 0;
	TAxis *axis = m3DHisto->GetXaxis();
	double val = axis->GetBinCenter(1);
	return isLogQ2() ? exp(val) : val;
	}

	double Table::maxQ2() const
	{
	if (m3DHisto == 0) return 0;
	TAxis *axis = m3DHisto->GetXaxis();
	double val = axis->GetBinCenter(axis->GetNbins());
	return isLogQ2() ? exp(val) : val;
	}

	double Table::minW2() const
	{
	if (m3DHisto == 0) return 0;
	TAxis *axis = m3DHisto->GetYaxis();
	double val = axis->GetBinCenter(1);
	return isLogW2() ? exp(val) : val;
	}

	double Table::maxW2() const
	{
	if (m3DHisto == 0) return 0;
	TAxis *axis = m3DHisto->GetYaxis();
	double val = axis->GetBinCenter(axis->GetNbins());
	return isLogW2() ? exp(val) : val;
	}

	double Table::minT() const
	{
	if (m3DHisto == 0) return 0;
	TAxis *axis = m3DHisto->GetZaxis();
	double val = axis->GetBinCenter(axis->GetNbins()); // t always as \|t\| in table
	return isLogT() ? -exp(val) : -val;
	}

	double Table::maxT() const
	{
	if (m3DHisto == 0) return 0;
	TAxis *axis = m3DHisto->GetZaxis();
	double val = axis->GetBinCenter(1); // t always as \|t\| in table
	double result = isLogT() ? -exp(val) : -val;
	if (result > 0) {
	// cout << "Table::maxT(): warning, t > 0, t (" << result << ") set to 0 now." << endl;
	result = 0;
	}
	return result;
	}

	string Table::filename() const {return mFilename;}

	void Table::list(ostream& os, bool printContent, bool printStatistics) const
	{
	ios::fmtflags fmt = os.flags(); // store current I/O flags

	if (!m3DHisto) {
	os << "Table::list(): table is undefined." << endl;
	return;
	}

	os << setw(11) << "Table ID = " << setw(14) << left << mID;
	os << setw(6) << right << "A = " << setw(4) << left << A();
	os << setw(18) << right << "content = " << (isLogContent() ? "log" : "");
	if (isMeanA())
	os << "<A>" << endl;
	else if (isMeanA2())
	os << "<A^2>" << endl;
	else if (isLambdaA())
	os << "lambda_<A>" << endl;
	else
	os << "unknown" << endl;
	os << setw(6+11+14) << right << "vmId = " << setw(4) <<left << vectorMesonId();
	os << setw(18) << right << "polarization = " << (isTransverse() ? 'T' : 'L') << endl;
	os << setw(6+11+14) << right << "model = ";
	if (dipoleModelType() == bSat)
	os << "bSat" << endl;
	else if (dipoleModelType() == bNonSat)
	os << "bNonSat" << endl;
	else
	os << "bCGC" << endl;
	os << setw(6+11+14+1) << right << "Q2 = [" << minQ2() << ", " << maxQ2()
	<< "; bins = " << m3DHisto->GetXaxis()->GetNbins() << (isLogQ2() ? "; log]" : "; lin]") << endl;
	os << setw(6+11+14+1) << right << "W2 = [" << minW2() << ", " << maxW2()
	<< "; bins = " << m3DHisto->GetYaxis()->GetNbins() << (isLogW2() ? "; log]" : "; lin]") << endl;
	os << setw(6+11+14+1) << right << "t = [" << minT() << ", " << maxT()
	<< "; bins = " << m3DHisto->GetZaxis()->GetNbins() << (isLogT() ? "; log]" : "; lin]") << endl;
	os << setw(6+11+14) << right << "file = " << mFilename.c_str() << endl;
	os << endl;

	double Q2, W2, t;
	int binx, biny, binz;

	if (printContent) {
	for (int i=0; i<numberOfEntries(); i++) {
	binXYZ(i, binx, biny, binz);
	binCenter(i, Q2, W2, t);
	double value = get(Q2, W2, t);
	double rawContent = m3DHisto->GetBinContent(binx, biny, binz);
	double error = m3DHisto->GetBinError(binx, biny, binz);
	os << "bin = " << setw(8) << left << i;
	os << "Q2 = " << setw(10) << left << fixed << setprecision(3) << Q2;
	os << "W2 = " << setw(11) << left << fixed << W2;
	os << "t = " << setw(14) << left << scientific << t;
	os << "value = " << setw(12) << left << scientific << value;
	os << "(binx = " << setw(5) << left << binx;
	os << "biny = " << setw(5) << left << biny;
	os << "binz = " << setw(5) << left << binz;
	os << "content = " << setw(12) << left << scientific << rawContent;
	os << "error = " << left << error<< ')';
	if ( (isLogContent() && rawContent <= log(numeric_limits<float>::min()*2)) \|\| (rawContent == 0)) cout << " Z";
	cout << endl;
	}
	os << endl;
	}

	if (printStatistics) {
	int nEmpty = 0;
	int nInvalid = 0;
	int nNegative = 0;
	double sum = 0;
	double maxContent = 0;
	double minContent = numeric_limits<float>::max();
	int minBin, maxBin;
	minBin = maxBin = 0;
	for (int i=0; i<numberOfEntries(); i++) {
	binXYZ(i, binx, biny, binz);
	double c = m3DHisto->GetBinContent(binx, biny, binz);
	if (isLogContent()) {
	if (c <= log(numeric_limits<float>::min()*2)) nEmpty++;
	}
	else {
	if (c == 0) nEmpty++;
	}
	if ( !isLogContent() && c < 0) nNegative++;
	if (isnan(c) \|\| isinf(c)) nInvalid++;
	if (isLogContent()) c = exp(c);
	if (c > maxContent) {maxContent = c; maxBin = i;}
	if (c > 0 && c < minContent) {minContent = c; minBin = i;}
	if (c > 0) sum += c;
	}

	os << setw(31) << right << "total number of cells = " << numberOfEntries() << endl;
	os << setw(31) << right << "cells with negative content = " << nNegative << endl;
	os << setw(31) << right << "cells with no (0) content = " << nEmpty << endl;
	os << setw(31) << right << "cells with invalid content = " << nInvalid << endl;
	binCenter(minBin, Q2, W2, t);
	os << setw(31) << right << "minimum content = " << minContent << " at Q2 = " << Q2 << ", W2 = " << W2 << ", t = " << t << endl;
	binCenter(maxBin, Q2, W2, t);
	os << setw(31) << right << "maximum content = " << maxContent << " at Q2 = " << Q2 << ", W2 = " << W2 << ", t = " << t << endl;
	os << setw(31) << right << "sum = " << sum << endl;
	os << endl;
	}

	os.flags(fmt); // restore I/O flags
	}

	double Table::InterpolateGridSpline(double x, double y, double z) const // Q2, W2, t
	{
	//
	// The algorithm was taken from Cristian Constantin Lalescu.
	// See http://arxiv.org/abs/0905.3564 for details.
	//
	// Grid points: 4, 6, or 8
	// Spline order: 3,5 for order 4
	// 3,5,7,9 for order 6
	// 3,5,7,9,11,13 for order 8
	//

	//
	// Find cell that refer to the position
	//
	int nx = m3DHisto->GetXaxis()->FindBin(x);
	int ny = m3DHisto->GetYaxis()->FindBin(y);
	int nz = m3DHisto->GetZaxis()->FindBin(z);

	//
	// Define the # of grid points depending on how far we are away
	// from the edge. If at the edge we fall back to linear interpolation
	// as provided by TH3.
	// There must be a smarter way doing this than the code below
	//
	int gridPoints = 0;

	if (nx-3 >= 1 && nx+4 <= m3DHisto->GetXaxis()->GetNbins() &&
	ny-3 >= 1 && ny+4 <= m3DHisto->GetYaxis()->GetNbins() &&
	nz-3 >= 1 && nz+4 <= m3DHisto->GetZaxis()->GetNbins()) {
	gridPoints = 8;
	}
	else if (nx-2 >= 1 && nx+3 <= m3DHisto->GetXaxis()->GetNbins() &&
	ny-2 >= 1 && ny+3 <= m3DHisto->GetYaxis()->GetNbins() &&
	nz-2 >= 1 && nz+3 <= m3DHisto->GetZaxis()->GetNbins()) {
	gridPoints = 6;
	}
	else if (nx-1 >= 1 && nx+2 <= m3DHisto->GetXaxis()->GetNbins() &&
	ny-1 >= 1 && ny+2 <= m3DHisto->GetYaxis()->GetNbins() &&
	nz-1 >= 1 && nz+2 <= m3DHisto->GetZaxis()->GetNbins()) {
	gridPoints = 4;
	}
	else {
	return m3DHisto->Interpolate(x,y,z);
	}

	//
	// Find bin centers
	//
	double xc = m3DHisto->GetXaxis()->GetBinCenter(nx);
	double yc = m3DHisto->GetYaxis()->GetBinCenter(ny);
	double zc = m3DHisto->GetZaxis()->GetBinCenter(nz);

	//
	// Define the scale for coordinate transformations
	// grid_spline expects x,y,z in bin units
	//
	double xscale = m3DHisto->GetXaxis()->GetBinWidth(1);
	double yscale = m3DHisto->GetYaxis()->GetBinWidth(1);
	double zscale = m3DHisto->GetZaxis()->GetBinWidth(1);

	//
	// Prepare grid spline alogorithm
	//
	double result, splineOrder;
	if (gridPoints == 8) {
	local_scal_3D<8> lf;
	for (int i=-3;i<=4;i++) {
	for(int j=-3;j<=4;j++) {
	lf(-3,i,j) = m3DHisto->GetBinContent(nx-3, ny+i, nz+j);
	lf(-2,i,j) = m3DHisto->GetBinContent(nx-2, ny+i, nz+j);
	lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
	lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
	lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
	lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
	lf( 3,i,j) = m3DHisto->GetBinContent(nx+3, ny+i, nz+j);
	lf( 4,i,j) = m3DHisto->GetBinContent(nx+4, ny+i, nz+j);
	}
	}
	splineOrder = 7;
	result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
	}
	else if (gridPoints == 6) {
	local_scal_3D<6> lf;
	for (int i=-2;i<=3;i++) {
	for(int j=-2;j<=3;j++) {
	lf(-2,i,j) = m3DHisto->GetBinContent(nx-2, ny+i, nz+j);
	lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
	lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
	lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
	lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
	lf( 3,i,j) = m3DHisto->GetBinContent(nx+3, ny+i, nz+j);
	}
	}
	splineOrder = 5;
	result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
	}
	else if (gridPoints == 4) {
	local_scal_3D<4> lf;
	for (int i=-1;i<=2;i++) {
	for(int j=-1;j<=2;j++) {
	lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
	lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
	lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
	lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
	}
	}
	splineOrder = 3;
	result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
	}
	else {
	cout << "Table::InterpolateGridSpline(): Error, illegal number of grid points." << endl;
	result = 0;
	}

	return result;
	}

	void Table::setAutobackup(const char* prefix, int freq)
	{
	mBackupPrefix = string(prefix);
	mBackupFrequence = freq;
	}

	void Table::backup(int backupBin)
	{
	ostringstream backupFilenameStream;
	backupFilenameStream << mBackupPrefix.c_str() << "_backup." << static_cast<int>(getpid())
	<< '.' << mID << '.' << backupBin << ".root";

	string filename = backupFilenameStream.str();
	TFile file(filename.c_str(),"RECREATE");
	m3DHisto->Write();
	file.Close();

	if (filename != mLastBackupFilename)
	remove(mLastBackupFilename.c_str());
	mLastBackupFilename = filename;
	time_t now = time(0);
	cout << "Table::backup(): autobackup performed, file = '" << mLastBackupFilename.c_str() << "', time = " << ctime(&now);
	}

	int Table::globalBin(int binx, int biny, int binz) const
	{
	int nbinx = m3DHisto->GetXaxis()->GetNbins();
	int nbiny = m3DHisto->GetYaxis()->GetNbins();

	return (binz-1)nbinxnbiny+(biny-1)*nbinx+(binx-1);
	}

File Metadata

Mime Type: text/x-c
Expires: Sat, Dec 21, 1:25 PM (1 d, 55 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4015342
Default Alt Text: Table.cpp (26 KB)

Table.cppNo OneTemporaryActions

Table.cppView Options

File Metadata

Event Timeline

Table.cpp
No OneTemporary
Actions

Table.cpp
View Options