Commit 6280f59a authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

cpp,py: binary files with _double_ values implemented

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@109 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 83a7762d
...@@ -52,7 +52,7 @@ public: ...@@ -52,7 +52,7 @@ public:
void getLocalData(double &aValue, double &aWeight, void getLocalData(double &aValue, double &aWeight,
std::vector<unsigned int>* &indLocal, std::vector<unsigned int>* &indLocal,
std::vector<double>* &derLocal); std::vector<double>* &derLocal);
void getAllData(float &fValue, float &fErr, void getAllData(double &aValue, double &aErr,
std::vector<unsigned int>* &indLocal, std::vector<unsigned int>* &indLocal,
std::vector<double>* &derLocal, std::vector<int>* &labGlobal, std::vector<double>* &derLocal, std::vector<int>* &labGlobal,
std::vector<double>* &derGlobal); std::vector<double>* &derGlobal);
......
...@@ -18,8 +18,10 @@ namespace gbl { ...@@ -18,8 +18,10 @@ namespace gbl {
/** /**
* Containing information for local (track) and global fit. * Containing information for local (track) and global fit.
* *
* The data blocks are collected in two arrays, a real array and * The data blocks are collected in two arrays, a real array
* an integer array, of same length. The content of the record is: * (containing float or double values) and integer array, of same length.
* A positive record length indicate _float_ and a negative one _double_ values.
* The content of the record is:
*\verbatim *\verbatim
* real array integer array * real array integer array
* 0 0.0 error count (this record) * 0 0.0 error count (this record)
...@@ -44,9 +46,9 @@ namespace gbl { ...@@ -44,9 +46,9 @@ namespace gbl {
class MilleBinary { class MilleBinary {
public: public:
MilleBinary(const std::string fileName = "milleBinaryISN.dat", MilleBinary(const std::string fileName = "milleBinaryISN.dat",
unsigned int aSize = 2000); bool doublePrec = false, unsigned int aSize = 2000);
virtual ~MilleBinary(); virtual ~MilleBinary();
void addData(float aMeas, float aPrec, void addData(double aMeas, double aPrec,
const std::vector<unsigned int> &indLocal, const std::vector<unsigned int> &indLocal,
const std::vector<double> &derLocal, const std::vector<double> &derLocal,
const std::vector<int> &labGlobal, const std::vector<int> &labGlobal,
...@@ -57,6 +59,8 @@ private: ...@@ -57,6 +59,8 @@ private:
std::ofstream binaryFile; ///< Binary File std::ofstream binaryFile; ///< Binary File
std::vector<int> intBuffer; ///< Integer buffer std::vector<int> intBuffer; ///< Integer buffer
std::vector<float> floatBuffer; ///< Float buffer std::vector<float> floatBuffer; ///< Float buffer
std::vector<double> doubleBuffer; ///< Double buffer
bool doublePrecision; ///< Flag for storage in as *double* values
}; };
} }
#endif /* MILLEBINARY_H_ */ #endif /* MILLEBINARY_H_ */
...@@ -210,18 +210,18 @@ void GblData::getLocalData(double &aValue, double &aWeight, ...@@ -210,18 +210,18 @@ void GblData::getLocalData(double &aValue, double &aWeight,
/// Get all Data for MP-II binary record. /// Get all Data for MP-II binary record.
/** /**
* \param [out] fValue Value * \param [out] aValue Value
* \param [out] fErr Error * \param [out] aErr Error
* \param [out] indLocal List of labels of local parameters * \param [out] indLocal List of labels of local parameters
* \param [out] derLocal List of derivatives for local parameters * \param [out] derLocal List of derivatives for local parameters
* \param [out] labGlobal List of labels of global parameters * \param [out] labGlobal List of labels of global parameters
* \param [out] derGlobal List of derivatives for global parameters * \param [out] derGlobal List of derivatives for global parameters
*/ */
void GblData::getAllData(float &fValue, float &fErr, void GblData::getAllData(double &aValue, double &aErr,
std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal, std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal,
std::vector<int>* &labGlobal, std::vector<double>* &derGlobal) { std::vector<int>* &labGlobal, std::vector<double>* &derGlobal) {
fValue = theValue; aValue = theValue;
fErr = 1.0 / sqrt(thePrecision); aErr = 1.0 / sqrt(thePrecision);
indLocal = &theParameters; indLocal = &theParameters;
derLocal = &theDerivatives; derLocal = &theDerivatives;
labGlobal = &globalLabels; labGlobal = &globalLabels;
......
...@@ -1001,8 +1001,8 @@ unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight, ...@@ -1001,8 +1001,8 @@ unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
/// Write valid trajectory to Millepede-II binary file. /// Write valid trajectory to Millepede-II binary file.
void GblTrajectory::milleOut(MilleBinary &aMille) { void GblTrajectory::milleOut(MilleBinary &aMille) {
float fValue; double aValue;
float fErr; double aErr;
std::vector<unsigned int>* indLocal; std::vector<unsigned int>* indLocal;
std::vector<double>* derLocal; std::vector<double>* derLocal;
std::vector<int>* labGlobal; std::vector<int>* labGlobal;
...@@ -1014,9 +1014,9 @@ void GblTrajectory::milleOut(MilleBinary &aMille) { ...@@ -1014,9 +1014,9 @@ void GblTrajectory::milleOut(MilleBinary &aMille) {
// data: measurements, kinks and external seed // data: measurements, kinks and external seed
std::vector<GblData>::iterator itData; std::vector<GblData>::iterator itData;
for (itData = theData.begin(); itData != theData.end(); ++itData) { for (itData = theData.begin(); itData != theData.end(); ++itData) {
itData->getAllData(fValue, fErr, indLocal, derLocal, labGlobal, itData->getAllData(aValue, aErr, indLocal, derLocal, labGlobal,
derGlobal); derGlobal);
aMille.addData(fValue, fErr, *indLocal, *derLocal, *labGlobal, aMille.addData(aValue, aErr, *indLocal, *derLocal, *labGlobal,
*derGlobal); *derGlobal);
} }
aMille.writeRecord(); aMille.writeRecord();
......
...@@ -13,14 +13,23 @@ namespace gbl { ...@@ -13,14 +13,23 @@ namespace gbl {
/// Create binary file. /// Create binary file.
/** /**
* \param [in] fileName File name * \param [in] fileName File name
* \param [in] doublePrec Flag for storage as double values
* \param [in] aSize Buffer size * \param [in] aSize Buffer size
*/ */
MilleBinary::MilleBinary(const std::string fileName, unsigned int aSize) : MilleBinary::MilleBinary(const std::string fileName, bool doublePrec,
binaryFile(fileName.c_str(), std::ios::binary | std::ios::out), intBuffer(), floatBuffer() { unsigned int aSize) :
binaryFile(fileName.c_str(), std::ios::binary | std::ios::out), intBuffer(), floatBuffer(), doubleBuffer(), doublePrecision(
doublePrec) {
intBuffer.reserve(aSize); intBuffer.reserve(aSize);
floatBuffer.reserve(aSize);
intBuffer.push_back(0); // first word is error counter intBuffer.push_back(0); // first word is error counter
floatBuffer.push_back(0.); if (doublePrecision) {
doubleBuffer.reserve(aSize);
doubleBuffer.push_back(0.);
} else {
floatBuffer.reserve(aSize);
floatBuffer.push_back(0.);
}
} }
MilleBinary::~MilleBinary() { MilleBinary::~MilleBinary() {
...@@ -36,23 +45,42 @@ MilleBinary::~MilleBinary() { ...@@ -36,23 +45,42 @@ MilleBinary::~MilleBinary() {
* \param [in] labGlobal List of labels of global parameters * \param [in] labGlobal List of labels of global parameters
* \param [in] derGlobal List of derivatives for global parameters * \param [in] derGlobal List of derivatives for global parameters
*/ */
void MilleBinary::addData(float aMeas, float aErr, void MilleBinary::addData(double aMeas, double aErr,
const std::vector<unsigned int> &indLocal, const std::vector<unsigned int> &indLocal,
const std::vector<double> &derLocal, const std::vector<int> &labGlobal, const std::vector<double> &derLocal, const std::vector<int> &labGlobal,
const std::vector<double> &derGlobal) { const std::vector<double> &derGlobal) {
intBuffer.push_back(0); if (doublePrecision) {
floatBuffer.push_back(aMeas); // double values
for (unsigned int i = 0; i < indLocal.size(); ++i) { intBuffer.push_back(0);
intBuffer.push_back(indLocal[i]); doubleBuffer.push_back(aMeas);
floatBuffer.push_back(derLocal[i]); for (unsigned int i = 0; i < indLocal.size(); ++i) {
} intBuffer.push_back(indLocal[i]);
intBuffer.push_back(0); doubleBuffer.push_back(derLocal[i]);
floatBuffer.push_back(aErr); }
for (unsigned int i = 0; i < labGlobal.size(); ++i) { intBuffer.push_back(0);
if (derGlobal[i]) { doubleBuffer.push_back(aErr);
intBuffer.push_back(labGlobal[i]); for (unsigned int i = 0; i < labGlobal.size(); ++i) {
floatBuffer.push_back(derGlobal[i]); if (derGlobal[i]) {
intBuffer.push_back(labGlobal[i]);
doubleBuffer.push_back(derGlobal[i]);
}
}
} else {
// float values
intBuffer.push_back(0);
floatBuffer.push_back(aMeas);
for (unsigned int i = 0; i < indLocal.size(); ++i) {
intBuffer.push_back(indLocal[i]);
floatBuffer.push_back(derLocal[i]);
}
intBuffer.push_back(0);
floatBuffer.push_back(aErr);
for (unsigned int i = 0; i < labGlobal.size(); ++i) {
if (derGlobal[i]) {
intBuffer.push_back(labGlobal[i]);
floatBuffer.push_back(derGlobal[i]);
}
} }
} }
} }
...@@ -60,15 +88,23 @@ void MilleBinary::addData(float aMeas, float aErr, ...@@ -60,15 +88,23 @@ void MilleBinary::addData(float aMeas, float aErr,
/// Write record to file. /// Write record to file.
void MilleBinary::writeRecord() { void MilleBinary::writeRecord() {
const int recordLength = intBuffer.size() * 2; const int recordLength =
(doublePrecision) ? -intBuffer.size() * 2 : intBuffer.size() * 2;
binaryFile.write(reinterpret_cast<const char*>(&recordLength), binaryFile.write(reinterpret_cast<const char*>(&recordLength),
sizeof(recordLength)); sizeof(recordLength));
binaryFile.write(reinterpret_cast<char*>(&floatBuffer[0]), if (doublePrecision)
floatBuffer.size() * sizeof(floatBuffer[0])); binaryFile.write(reinterpret_cast<char*>(&doubleBuffer[0]),
doubleBuffer.size() * sizeof(doubleBuffer[0]));
else
binaryFile.write(reinterpret_cast<char*>(&floatBuffer[0]),
floatBuffer.size() * sizeof(floatBuffer[0]));
binaryFile.write(reinterpret_cast<char*>(&intBuffer[0]), binaryFile.write(reinterpret_cast<char*>(&intBuffer[0]),
intBuffer.size() * sizeof(intBuffer[0])); intBuffer.size() * sizeof(intBuffer[0]));
// start with new record // start with new record
intBuffer.resize(1); intBuffer.resize(1);
floatBuffer.resize(1); if (doublePrecision)
doubleBuffer.resize(1);
else
floatBuffer.resize(1);
} }
} }
...@@ -585,12 +585,13 @@ class GblTrajectory(object): ...@@ -585,12 +585,13 @@ class GblTrajectory(object):
for d in self.__data: for d in self.__data:
d.printData() d.printData()
## Write (data blocks of) trajectory to MP (binary) file. ## Write (data blocks of) trajectory to MP (binary) file (as *float* or *double* values).
# #
# @param aFile MP file # @param aFile MP file
# @param doublePrec flag for storage in as *double* values
# #
def milleOut(self, aFile): def milleOut(self, aFile, doublePrec=False):
rec = MilleRecord() rec = MilleRecord(doublePrec)
# data measurements and kinks # data measurements and kinks
for aData in self.__data: for aData in self.__data:
rec.addData(aData.toRecord()) rec.addData(aData.toRecord())
......
...@@ -15,9 +15,11 @@ import array, math ...@@ -15,9 +15,11 @@ import array, math
# #
# Containing information for local (track) and global fit. # Containing information for local (track) and global fit.
# #
# The data blocks are collected in two arrays, a real array and # The data blocks are collected in two arrays, a real array
# an integer array, of same length. The content of the arrays:: # (containing float or double values) and integer array, of same length.
# # A positive record length indicate _float_ and a negative one _double_ values.
# The content of the record is:
#
#\verbatim #\verbatim
# real array integer array # real array integer array
# 0 0.0 error count (this record) # 0 0.0 error count (this record)
...@@ -43,7 +45,9 @@ class MilleRecord(object): ...@@ -43,7 +45,9 @@ class MilleRecord(object):
## Create MP-II binary record. ## Create MP-II binary record.
# #
def __init__(self): def __init__(self, doublePrec=False):
## flag for storage in as *double* values
self.__doublePrecision = doublePrec
## position in record, usually start of next data block; int ## position in record, usually start of next data block; int
self.__position = 1 self.__position = 1
## number of data blocks in record; int ## number of data blocks in record; int
...@@ -56,15 +60,15 @@ class MilleRecord(object): ...@@ -56,15 +60,15 @@ class MilleRecord(object):
self.__iErr = 0 self.__iErr = 0
## array with markers (0) and labels; array(int32) ## array with markers (0) and labels; array(int32)
self.__inder = array.array('i') self.__inder = array.array('i')
## array with values, errors and derivatives; (float32) ## array with values, errors and derivatives; (float32 or float64)
self.__glder = array.array('f') self.__glder = array.array('d' if doublePrec else 'f')
## Add data block to (end of) record. ## Add data block to (end of) record.
# #
# @param dataList list with measurement, error, labels and derivatives; list # @param dataList list with measurement, error, labels and derivatives; list
# #
def addData(self, dataList): def addData(self, dataList):
if (self.__numData == 0): # first word is error counter if (self.__numData == 0): # first word is error counter
self.__inder.append(0) self.__inder.append(0)
self.__glder.append(0.) self.__glder.append(0.)
self.__numData += 1 self.__numData += 1
...@@ -75,7 +79,7 @@ class MilleRecord(object): ...@@ -75,7 +79,7 @@ class MilleRecord(object):
self.__inder.fromlist(indLocal) self.__inder.fromlist(indLocal)
self.__glder.fromlist(derLocal) self.__glder.fromlist(derLocal)
self.__inder.append(0) self.__inder.append(0)
self.__glder.append(1.0 / math.sqrt(aPrec)) # convert to error self.__glder.append(1.0 / math.sqrt(aPrec)) # convert to error
self.__inder.fromlist(labGlobal) self.__inder.fromlist(labGlobal)
self.__glder.fromlist(derGlobal) self.__glder.fromlist(derGlobal)
...@@ -90,7 +94,7 @@ class MilleRecord(object): ...@@ -90,7 +94,7 @@ class MilleRecord(object):
for i in range(self.__iMeas + 1, self.__iErr): for i in range(self.__iMeas + 1, self.__iErr):
indLocal.append(self.__inder[i]) indLocal.append(self.__inder[i])
derLocal.append(self.__glder[i]) derLocal.append(self.__glder[i])
aPrec = 1.0 / self.__glder[self.__iErr] ** 2 # convert to precision aPrec = 1.0 / self.__glder[self.__iErr] ** 2 # convert to precision
indGlobal = [] indGlobal = []
derGlobal = [] derGlobal = []
for i in range(self.__iErr + 1, self.__position): for i in range(self.__iErr + 1, self.__position):
...@@ -110,8 +114,8 @@ class MilleRecord(object): ...@@ -110,8 +114,8 @@ class MilleRecord(object):
# @param aFile (binary) file # @param aFile (binary) file
# #
def writeRecord(self, aFile): def writeRecord(self, aFile):
header = array.array('i') # header with number of words header = array.array('i') # header with number of words
header.append(len(self.__inder) * 2) header.append(-len(self.__inder) * 2 if self.__doublePrecision else len(self.__inder) * 2)
header.tofile(aFile) header.tofile(aFile)
self.__glder.tofile(aFile) self.__glder.tofile(aFile)
self.__inder.tofile(aFile) self.__inder.tofile(aFile)
...@@ -121,9 +125,11 @@ class MilleRecord(object): ...@@ -121,9 +125,11 @@ class MilleRecord(object):
# @param aFile (binary) file # @param aFile (binary) file
# #
def readRecord(self, aFile): def readRecord(self, aFile):
header = array.array('i') # header with number of words header = array.array('i') # header with number of words
header.fromfile(aFile, 1) header.fromfile(aFile, 1)
self.__recLen = header[0] / 2 self.__recLen = abs(header[0] / 2)
if header[0] < 0:
self.__glder = array.array('d')
self.__glder.fromfile(aFile, self.__recLen) self.__glder.fromfile(aFile, self.__recLen)
self.__inder.fromfile(aFile, self.__recLen) self.__inder.fromfile(aFile, self.__recLen)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment