Mille.cc

Go to the documentation of this file.

#include "Mille.h"

#include <iostream>
#include "FWCore/MessageLogger/interface/MessageLogger.h"

//___________________________________________________________________________

Mille::Mille(const char *outFileName, bool asBinary, bool writeZero) :
  fileMode_(asBinary ? (std::ios::binary | std::ios::out) : std::ios::out),
  fileName_(outFileName),
  outFile_(fileName_, fileMode_),
  asBinary_(asBinary), writeZero_(writeZero), bufferPos_(-1), hasSpecial_(false)
{
  // opens outFileName, by default as binary file

  // Instead bufferPos_(-1), hasSpecial_(false) and the following two lines
  // we could call newSet() and kill()...
  bufferInt_[0]   = 0;
  bufferFloat_[0] = 0.;

  if (!outFile_.is_open()) {
    edm::LogError("Alignment")
      << "Mille::Mille: Could not open " << fileName_ << " as output file.";
  }
}

//___________________________________________________________________________

Mille::~Mille()
{
  // closes file
  outFile_.close();
}

//___________________________________________________________________________

void Mille::mille(int NLC, const float *derLc,
                  int NGL, const float *derGl, const int *label,
                  float rMeas, float sigma)
{
  if (sigma <= 0.) return;
  if (bufferPos_ == -1) this->newSet(); // start, e.g. new track
  if (!this->checkBufferSize(NLC, NGL)) return;

  // first store measurement
  ++bufferPos_;
  bufferFloat_[bufferPos_] = rMeas;
  bufferInt_  [bufferPos_] = 0;

  // store local derivatives and local 'lables' 1,...,NLC
  for (int i = 0; i < NLC; ++i) {
    if (derLc[i] || writeZero_) { // by default store only non-zero derivatives
      ++bufferPos_;
      bufferFloat_[bufferPos_] = derLc[i]; // local derivatives
      bufferInt_  [bufferPos_] = i+1;      // index of local parameter
    }
  }

  // store uncertainty of measurement in between locals and globals
  ++bufferPos_;
  bufferFloat_[bufferPos_] = sigma;
  bufferInt_  [bufferPos_] = 0;

  // store global derivatives and their lables
  for (int i = 0; i < NGL; ++i) {
    if (derGl[i] || writeZero_) { // by default store only non-zero derivatives
      if ((label[i] > 0 || writeZero_) && label[i] <= maxLabel_) { // and for valid labels
        ++bufferPos_;
        bufferFloat_[bufferPos_] = derGl[i]; // global derivatives
        bufferInt_  [bufferPos_] = label[i]; // index of global parameter
      } else {
        edm::LogError("Alignment")
          << "Mille::mille: Invalid label " << label[i]
          << " <= 0 or > " << maxLabel_;
      }
    }
  }
}

//___________________________________________________________________________
void Mille::special(int nSpecial, const float *floatings, const int *integers)
{
  if (nSpecial == 0) return;
  if (bufferPos_ == -1) this->newSet(); // start, e.g. new track
  if (hasSpecial_) {
    edm::LogError("Alignment")
      << "Mille::special: Special values already stored for this record.";
    return;
  }
  if (!this->checkBufferSize(nSpecial, 0)) return;
  hasSpecial_ = true; // after newSet() (Note: MILLSP sets to buffer position...)

  //  bufferFloat_[.]   | bufferInt_[.]
  // ------------------------------------
  //      0.0           |      0
  //  -float(nSpecial)  |      0
  //  The above indicates special data, following are nSpecial floating and nSpecial integer data.

  ++bufferPos_; // zero pair
  bufferFloat_[bufferPos_] = 0.;
  bufferInt_  [bufferPos_] = 0;

  ++bufferPos_; // nSpecial and zero
  bufferFloat_[bufferPos_] = -nSpecial; // automatic conversion to float
  bufferInt_  [bufferPos_] = 0;

  for (int i = 0; i < nSpecial; ++i) {
    ++bufferPos_;
    bufferFloat_[bufferPos_] = floatings[i];
    bufferInt_  [bufferPos_] = integers[i];
  }
}

//___________________________________________________________________________

void Mille::kill()
{
  // reset buffers, i.e. kill derivatives accumulated for current set
  bufferPos_ = -1;
}

//___________________________________________________________________________


void Mille::flushOutputFile() {
  // flush output file
  outFile_.flush();
}

//___________________________________________________________________________


void Mille::resetOutputFile() {
  // flush output file
  outFile_.close();
  outFile_.open(fileName_, fileMode_);
  if (!outFile_.is_open()) {
    edm::LogError("Alignment")
      << "Mille::resetOutputFile: Could not reopen " << fileName_ << ".";
  }
}

//___________________________________________________________________________

void Mille::end()
{
  // write set of derivatives with same local parameters to file
  if (bufferPos_ > 0) { // only if anything stored...
    const int numWordsToWrite = (bufferPos_ + 1)*2;

    if (asBinary_) {
      outFile_.write(reinterpret_cast<const char*>(&numWordsToWrite),
                     sizeof(numWordsToWrite));
      outFile_.write(reinterpret_cast<char*>(bufferFloat_),
                     (bufferPos_+1) * sizeof(bufferFloat_[0]));
      outFile_.write(reinterpret_cast<char*>(bufferInt_),
                     (bufferPos_+1) * sizeof(bufferInt_[0]));
    } else {
      outFile_ << numWordsToWrite << "\n";
      for (int i = 0; i < bufferPos_+1; ++i) {
        outFile_ << bufferFloat_[i] << " ";
      }
      outFile_ << "\n";

      for (int i = 0; i < bufferPos_+1; ++i) {
        outFile_ << bufferInt_[i] << " ";
      }
      outFile_ << "\n";
    }
  }
  bufferPos_ = -1; // reset buffer for next set of derivatives
}

//___________________________________________________________________________

void Mille::newSet()
{
  // initilise for new set of locals, e.g. new track
  bufferPos_ = 0;
  hasSpecial_ = false;
  bufferFloat_[0] = 0.0;
  bufferInt_  [0] = 0;   // position 0 used as error counter
}

//___________________________________________________________________________

bool Mille::checkBufferSize(int nLocal, int nGlobal)
{
  // enough space for next nLocal + nGlobal derivatives incl. measurement?

  if (bufferPos_ + nLocal + nGlobal + 2 >= bufferSize_) {
    ++(bufferInt_[0]); // increase error count
    edm::LogError("Alignment")
      << "Mille::checkBufferSize: Buffer too short ("
      << bufferSize_ << "),"
      << "\n need space for nLocal (" << nLocal<< ")"
      << "/nGlobal (" << nGlobal << ") local/global derivatives, "
      << bufferPos_ + 1 << " already stored!";
    return false;
  } else {
    return true;
  }
}