MillePedeAlignmentAlgorithm/src/Mille.cc

0001 /**
0002  * \file Mille.cc
0003  *  \author    : Gero Flucke
0004  *  date       : October 2006
0005  *  $Revision: 1.3 $
0006  *  $Date: 2007/04/16 17:47:38 $
0007  *  (last update by $Author: flucke $)
0008  */
0009
0010 #include "Mille.h"
0011
0012 #include <iostream>
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 //___________________________________________________________________________
0016
0017 Mille::Mille(const char *outFileName, bool asBinary, bool writeZero)
0018     : fileMode_(asBinary ? (std::ios::binary | std::ios::out) : std::ios::out),
0019       fileName_(outFileName),
0020       outFile_(fileName_, fileMode_),
0021       asBinary_(asBinary),
0022       writeZero_(writeZero),
0023       bufferPos_(-1),
0024       hasSpecial_(false) {
0025   // opens outFileName, by default as binary file
0026
0027   // Instead bufferPos_(-1), hasSpecial_(false) and the following two lines
0028   // we could call newSet() and kill()...
0029   bufferInt_[0] = 0;
0030   bufferFloat_[0] = 0.;
0031
0032   if (!outFile_.is_open()) {
0033     edm::LogError("Alignment") << "Mille::Mille: Could not open " << fileName_ << " as output file.";
0034   }
0035 }
0036
0037 //___________________________________________________________________________
0038
0039 Mille::~Mille() {
0040   // closes file
0041   outFile_.close();
0042 }
0043
0044 //___________________________________________________________________________
0045
0046 void Mille::mille(int NLC, const float *derLc, int NGL, const float *derGl, const int *label, float rMeas, float sigma) {
0047   if (sigma <= 0.)
0048     return;
0049   if (bufferPos_ == -1)
0050     this->newSet();  // start, e.g. new track
0051   if (!this->checkBufferSize(NLC, NGL))
0052     return;
0053
0054   // first store measurement
0055   ++bufferPos_;
0056   bufferFloat_[bufferPos_] = rMeas;
0057   bufferInt_[bufferPos_] = 0;
0058
0059   // store local derivatives and local 'lables' 1,...,NLC
0060   for (int i = 0; i < NLC; ++i) {
0061     if (derLc[i] || writeZero_) {  // by default store only non-zero derivatives
0062       ++bufferPos_;
0063       bufferFloat_[bufferPos_] = derLc[i];  // local derivatives
0064       bufferInt_[bufferPos_] = i + 1;       // index of local parameter
0065     }
0066   }
0067
0068   // store uncertainty of measurement in between locals and globals
0069   ++bufferPos_;
0070   bufferFloat_[bufferPos_] = sigma;
0071   bufferInt_[bufferPos_] = 0;
0072
0073   // store global derivatives and their lables
0074   for (int i = 0; i < NGL; ++i) {
0075     if (derGl[i] || writeZero_) {                                   // by default store only non-zero derivatives
0076       if ((label[i] > 0 || writeZero_) && label[i] <= maxLabel_) {  // and for valid labels
0077         ++bufferPos_;
0078         bufferFloat_[bufferPos_] = derGl[i];  // global derivatives
0079         bufferInt_[bufferPos_] = label[i];    // index of global parameter
0080       } else {
0081         edm::LogError("Alignment") << "Mille::mille: Invalid label " << label[i] << " <= 0 or > " << maxLabel_;
0082       }
0083     }
0084   }
0085 }
0086
0087 //___________________________________________________________________________
0088 void Mille::special(int nSpecial, const float *floatings, const int *integers) {
0089   if (nSpecial == 0)
0090     return;
0091   if (bufferPos_ == -1)
0092     this->newSet();  // start, e.g. new track
0093   if (hasSpecial_) {
0094     edm::LogError("Alignment") << "Mille::special: Special values already stored for this record.";
0095     return;
0096   }
0097   if (!this->checkBufferSize(nSpecial, 0))
0098     return;
0099   hasSpecial_ = true;  // after newSet() (Note: MILLSP sets to buffer position...)
0100
0101   //  bufferFloat_[.]   | bufferInt_[.]
0102   // ------------------------------------
0103   //      0.0           |      0
0104   //  -float(nSpecial)  |      0
0105   //  The above indicates special data, following are nSpecial floating and nSpecial integer data.
0106
0107   ++bufferPos_;  // zero pair
0108   bufferFloat_[bufferPos_] = 0.;
0109   bufferInt_[bufferPos_] = 0;
0110
0111   ++bufferPos_;                          // nSpecial and zero
0112   bufferFloat_[bufferPos_] = -nSpecial;  // automatic conversion to float
0113   bufferInt_[bufferPos_] = 0;
0114
0115   for (int i = 0; i < nSpecial; ++i) {
0116     ++bufferPos_;
0117     bufferFloat_[bufferPos_] = floatings[i];
0118     bufferInt_[bufferPos_] = integers[i];
0119   }
0120 }
0121
0122 //___________________________________________________________________________
0123
0124 void Mille::kill() {
0125   // reset buffers, i.e. kill derivatives accumulated for current set
0126   bufferPos_ = -1;
0127 }
0128
0129 //___________________________________________________________________________
0130
0131 void Mille::flushOutputFile() {
0132   // flush output file
0133   outFile_.flush();
0134 }
0135
0136 //___________________________________________________________________________
0137
0138 void Mille::resetOutputFile() {
0139   // flush output file
0140   outFile_.close();
0141   outFile_.open(fileName_, fileMode_);
0142   if (!outFile_.is_open()) {
0143     edm::LogError("Alignment") << "Mille::resetOutputFile: Could not reopen " << fileName_ << ".";
0144   }
0145 }
0146
0147 //___________________________________________________________________________
0148
0149 void Mille::end() {
0150   // write set of derivatives with same local parameters to file
0151   if (bufferPos_ > 0) {  // only if anything stored...
0152     const int numWordsToWrite = (bufferPos_ + 1) * 2;
0153
0154     if (asBinary_) {
0155       outFile_.write(reinterpret_cast<const char *>(&numWordsToWrite), sizeof(numWordsToWrite));
0156       outFile_.write(reinterpret_cast<char *>(bufferFloat_), (bufferPos_ + 1) * sizeof(bufferFloat_[0]));
0157       outFile_.write(reinterpret_cast<char *>(bufferInt_), (bufferPos_ + 1) * sizeof(bufferInt_[0]));
0158     } else {
0159       outFile_ << numWordsToWrite << "\n";
0160       for (int i = 0; i < bufferPos_ + 1; ++i) {
0161         outFile_ << bufferFloat_[i] << " ";
0162       }
0163       outFile_ << "\n";
0164
0165       for (int i = 0; i < bufferPos_ + 1; ++i) {
0166         outFile_ << bufferInt_[i] << " ";
0167       }
0168       outFile_ << "\n";
0169     }
0170   }
0171   bufferPos_ = -1;  // reset buffer for next set of derivatives
0172 }
0173
0174 //___________________________________________________________________________
0175
0176 void Mille::newSet() {
0177   // initilise for new set of locals, e.g. new track
0178   bufferPos_ = 0;
0179   hasSpecial_ = false;
0180   bufferFloat_[0] = 0.0;
0181   bufferInt_[0] = 0;  // position 0 used as error counter
0182 }
0183
0184 //___________________________________________________________________________
0185
0186 bool Mille::checkBufferSize(int nLocal, int nGlobal) {
0187   // enough space for next nLocal + nGlobal derivatives incl. measurement?
0188
0189   if (bufferPos_ + nLocal + nGlobal + 2 >= bufferSize_) {
0190     ++(bufferInt_[0]);  // increase error count
0191     edm::LogError("Alignment") << "Mille::checkBufferSize: Buffer too short (" << bufferSize_ << "),"
0192                                << "\n need space for nLocal (" << nLocal << ")"
0193                                << "/nGlobal (" << nGlobal << ") local/global derivatives, " << bufferPos_ + 1
0194                                << " already stored!";
0195     return false;
0196   } else {
0197     return true;
0198   }
0199 }