packages/TrackerMillepedeAlignment/Mille.cc

0001
0002 /** \file
0003  *  Create Millepede-II C-binary record.
0004  *
0005  *  \author Gero Flucke, University Hamburg, 2006
0006  *
0007  *  \copyright
0008  *  Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton,
0009  *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
0010  *  This library is free software; you can redistribute it and/or modify
0011  *  it under the terms of the GNU Library General Public License as
0012  *  published by the Free Software Foundation; either version 2 of the
0013  *  License, or (at your option) any later version. \n\n
0014  *  This library is distributed in the hope that it will be useful,
0015  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
0016  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0017  *  GNU Library General Public License for more details. \n\n
0018  *  You should have received a copy of the GNU Library General Public
0019  *  License along with this program (see the file COPYING.LIB for more
0020  *  details); if not, write to the Free Software Foundation, Inc.,
0021  *  675 Mass Ave, Cambridge, MA 02139, USA.
0022  */
0023
0024 /**
0025  *  author    : Gero Flucke, University Hamburg, 2006
0026  *  date      : October 2006
0027  *  $Revision: 1.3 $
0028  *  $Date: 2007/04/16 17:47:38 $
0029  *  (last update by $Author: flucke $)
0030  */
0031
0032 #include "Mille.h"
0033
0034 #include <fstream>
0035 #include <iostream>
0036
0037 //___________________________________________________________________________
0038
0039 /// Opens outFileName (by default as binary file).
0040 /**
0041  * \param[in] outFileName  file name
0042  * \param[in] asBinary     flag for binary
0043  * \param[in] writeZero    flag for keeping of zeros
0044  */
0045 Mille::Mille(const char *outFileName, bool asBinary, bool writeZero)
0046   : myOutFile(outFileName, (asBinary ? (std::ios::binary | std::ios::out) : std::ios::out))
0047   , myAsBinary(asBinary)
0048   , myWriteZero(writeZero)
0049   , myBufferPos(-1)
0050   , myHasSpecial(false)
0051 {
0052   // Instead myBufferPos(-1), myHasSpecial(false) and the following two lines
0053   // we could call newSet() and kill()...
0054   myBufferInt[0] = 0;
0055   myBufferFloat[0] = 0.;
0056
0057   if (!myOutFile.is_open())
0058   {
0059     std::cerr << "Mille::Mille: Could not open " << outFileName
0060               << " as output file." << std::endl;
0061   }
0062 }
0063
0064 //___________________________________________________________________________
0065 /// Closes file.
0066 Mille::~Mille()
0067 {
0068   myOutFile.close();
0069 }
0070
0071 //___________________________________________________________________________
0072 /// Add measurement to buffer.
0073 /**
0074  * \param[in]    NLC    number of local derivatives
0075  * \param[in]    derLc  local derivatives
0076  * \param[in]    NGL    number of global derivatives
0077  * \param[in]    derGl  global derivatives
0078  * \param[in]    label  global labels
0079  * \param[in]    rMeas  measurement (residuum)
0080  * \param[in]    sigma  error
0081  */
0082 void Mille::mille(int NLC, const float *derLc,
0083                   int NGL, const float *derGl, const int *label,
0084                   float rMeas, float sigma)
0085 {
0086   if (sigma <= 0.)
0087   {
0088     return;
0089   }
0090   if (myBufferPos == -1)
0091   {
0092     this->newSet();  // start, e.g. new track
0093   }
0094   if (!this->checkBufferSize(NLC, NGL))
0095   {
0096     return;
0097   }
0098
0099   // first store measurement
0100   ++myBufferPos;
0101   myBufferFloat[myBufferPos] = rMeas;
0102   myBufferInt[myBufferPos] = 0;
0103
0104   // store local derivatives and local 'lables' 1,...,NLC
0105   for (int i = 0; i < NLC; ++i)
0106   {
0107     if (derLc[i] || myWriteZero)
0108     {  // by default store only non-zero derivatives
0109       ++myBufferPos;
0110       myBufferFloat[myBufferPos] = derLc[i];  // local derivatives
0111       myBufferInt[myBufferPos] = i + 1;       // index of local parameter
0112     }
0113   }
0114
0115   // store uncertainty of measurement in between locals and globals
0116   ++myBufferPos;
0117   myBufferFloat[myBufferPos] = sigma;
0118   myBufferInt[myBufferPos] = 0;
0119
0120   // store global derivatives and their labels
0121   for (int i = 0; i < NGL; ++i)
0122   {
0123     if (derGl[i] || myWriteZero)
0124     {  // by default store only non-zero derivatives
0125       if ((label[i] > 0 || myWriteZero) && label[i] <= myMaxLabel)
0126       {  // and for valid labels
0127         ++myBufferPos;
0128         myBufferFloat[myBufferPos] = derGl[i];  // global derivatives
0129         myBufferInt[myBufferPos] = label[i];    // index of global parameter
0130       }
0131       else
0132       {
0133         std::cerr << "Mille::mille: Invalid label " << label[i]
0134                   << " <= 0 or > " << myMaxLabel << std::endl;
0135       }
0136     }
0137   }
0138 }
0139
0140 //___________________________________________________________________________
0141 /// Add special data to buffer.
0142 /**
0143  * \param[in]    nSpecial   number of floats/ints
0144  * \param[in]    floatings  floats
0145  * \param[in]    integers   ints
0146  */
0147 void Mille::special(int nSpecial, const float *floatings, const int *integers)
0148 {
0149   if (nSpecial == 0)
0150   {
0151     return;
0152   }
0153   if (myBufferPos == -1)
0154   {
0155     this->newSet();  // start, e.g. new track
0156   }
0157   if (myHasSpecial)
0158   {
0159     std::cerr << "Mille::special: Special values already stored for this record."
0160               << std::endl;
0161     return;
0162   }
0163   if (!this->checkBufferSize(nSpecial, 0))
0164   {
0165     return;
0166   }
0167   myHasSpecial = true;  // after newSet() (Note: MILLSP sets to buffer position...)
0168
0169   //  myBufferFloat[.]  | myBufferInt[.]
0170   // ------------------------------------
0171   //      0.0           |      0
0172   //  -float(nSpecial)  |      0
0173   //  The above indicates special data, following are nSpecial floating and nSpecial integer data.
0174
0175   ++myBufferPos;  // zero pair
0176   myBufferFloat[myBufferPos] = 0.;
0177   myBufferInt[myBufferPos] = 0;
0178
0179   ++myBufferPos;                           // nSpecial and zero
0180   myBufferFloat[myBufferPos] = -nSpecial;  // automatic conversion to float
0181   myBufferInt[myBufferPos] = 0;
0182
0183   for (int i = 0; i < nSpecial; ++i)
0184   {
0185     ++myBufferPos;
0186     myBufferFloat[myBufferPos] = floatings[i];
0187     myBufferInt[myBufferPos] = integers[i];
0188   }
0189 }
0190
0191 //___________________________________________________________________________
0192 /// Reset buffers, i.e. kill derivatives accumulated for current set.
0193 void Mille::kill()
0194 {
0195   myBufferPos = -1;
0196 }
0197
0198 //___________________________________________________________________________
0199 /// Write buffer (set of derivatives with same local parameters) to file.
0200 void Mille::end()
0201 {
0202   // std:: cout << " Mille::end() called with myBufferPos " << myBufferPos << std::endl;
0203
0204   if (myBufferPos > 0)
0205   {  // only if anything stored...
0206     const int numWordsToWrite = (myBufferPos + 1) * 2;
0207
0208     if (myAsBinary)
0209     {
0210       myOutFile.write(reinterpret_cast<const char *>(&numWordsToWrite),
0211                       sizeof(numWordsToWrite));
0212       myOutFile.write(reinterpret_cast<char *>(myBufferFloat),
0213                       (myBufferPos + 1) * sizeof(myBufferFloat[0]));
0214       myOutFile.write(reinterpret_cast<char *>(myBufferInt),
0215                       (myBufferPos + 1) * sizeof(myBufferInt[0]));
0216     }
0217     else
0218     {
0219       myOutFile << numWordsToWrite << "\n";
0220       for (int i = 0; i < myBufferPos + 1; ++i)
0221       {
0222         myOutFile << myBufferFloat[i] << " ";
0223       }
0224       myOutFile << "\n";
0225
0226       for (int i = 0; i < myBufferPos + 1; ++i)
0227       {
0228         myOutFile << myBufferInt[i] << " ";
0229       }
0230       myOutFile << "\n";
0231     }
0232   }
0233   myBufferPos = -1;  // reset buffer for next set of derivatives
0234
0235   //  std:: cout << " Mille::end() finished with myBufferPos " << myBufferPos << std::endl;
0236 }
0237
0238 //___________________________________________________________________________
0239 /// Initialize for new set of locals, e.g. new track.
0240 void Mille::newSet()
0241 {
0242   myBufferPos = 0;
0243   myHasSpecial = false;
0244   myBufferFloat[0] = 0.0;
0245   myBufferInt[0] = 0;  // position 0 used as error counter
0246 }
0247
0248 //___________________________________________________________________________
0249 /// Enough space for next nLocal + nGlobal derivatives incl. measurement?
0250 /**
0251  * \param[in]   nLocal  number of local derivatives
0252  * \param[in]   nGlobal number of global derivatives
0253  * \return      true if sufficient space available (else false)
0254  */
0255 bool Mille::checkBufferSize(int nLocal, int nGlobal)
0256 {
0257   if (myBufferPos + nLocal + nGlobal + 2 >= myBufferSize)
0258   {
0259     ++(myBufferInt[0]);  // increase error count
0260     std::cerr << "Mille::checkBufferSize: Buffer too short ("
0261               << myBufferSize << "),"
0262               << "\n need space for nLocal (" << nLocal << ")"
0263               << "/nGlobal (" << nGlobal << ") local/global derivatives, "
0264               << myBufferPos + 1 << " already stored!"
0265               << std::endl;
0266     return false;
0267   }
0268   else
0269   {
0270     return true;
0271   }
0272 }