Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
/**
 * \file Mille.cc
 *  \author    : Gero Flucke
 *  date       : October 2006
 *  $Revision: 1.3 $
 *  $Date: 2007/04/16 17:47:38 $
 *  (last update by $Author: flucke $)
 */

#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;
  }
}