MillePedeAlignmentAlgorithm/src/PedeReader.cc

0001 /**
0002  * \file PedeReader.cc
0003  *
0004  *  \author    : Gero Flucke
0005  *  date       : November 2006
0006  *  $Revision: 1.14 $
0007  *  $Date: 2012/08/10 09:01:11 $
0008  *  (last update by $Author: flucke $)
0009  */
0010
0011 #include "PedeReader.h"
0012 #include "PedeSteerer.h"
0013
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "Alignment/CommonAlignment/interface/Alignable.h"
0018 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
0019
0020 #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
0021
0022 #include "Alignment/CommonAlignmentParametrization/interface/RigidBodyAlignmentParameters.h"
0023
0024 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeVariables.h"
0025
0026 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0027
0028 #include <set>
0029 #include <string>
0030 #include <sstream>
0031
0032 const unsigned int PedeReader::myMaxNumValPerParam = 5;
0033
0034 //__________________________________________________________________________________________________
0035 PedeReader::PedeReader(const edm::ParameterSet &config,
0036                        const PedeSteerer &steerer,
0037                        const PedeLabelerBase &labels,
0038                        const RunRange &runrange)
0039     : mySteerer(steerer), myLabels(labels), myRunRange(runrange) {
0040   std::string pedeResultFile(config.getUntrackedParameter<std::string>("fileDir"));
0041   if (pedeResultFile.empty())
0042     pedeResultFile = steerer.directory();  // includes final '/'
0043   else if (pedeResultFile.find_last_of('/') != pedeResultFile.size() - 1) {
0044     pedeResultFile += '/';  // directory may need '/'
0045   }
0046
0047   pedeResultFile += config.getParameter<std::string>("readFile");
0048   myPedeResult.open(pedeResultFile.c_str(), std::ios::in);
0049   if (!myPedeResult.is_open()) {
0050     edm::LogError("Alignment") << "@SUB=PedeReader"
0051                                << "Problem opening pede output file " << pedeResultFile;
0052   }
0053 }
0054
0055 //__________________________________________________________________________________________________
0056 bool PedeReader::read(align::Alignables &alignables, bool setUserVars) {
0057   alignables.clear();
0058   myPedeResult.seekg(0, std::ios::beg);  // back to start
0059   bool isAllOk = true;
0060
0061   std::set<Alignable *> uniqueList;
0062
0063   edm::LogInfo("Alignment") << "@SUB=PedeReader::read"
0064                             << "will read parameters for run range " << myRunRange.first << " - " << myRunRange.second;
0065
0066   // loop on lines of text file
0067   unsigned int nParam = 0, nParamCalib = 0, nParamUnknown = 0;
0068   while (myPedeResult.good() && !myPedeResult.eof()) {
0069     // read label
0070     unsigned int paramLabel = 0;
0071     if (!this->readIfSameLine<unsigned int>(myPedeResult, paramLabel))
0072       continue;  // empty line?
0073
0074     // read up to maximal number of pede result per parameter
0075     float buffer[myMaxNumValPerParam] = {0.};
0076     unsigned int bufferPos = 0;
0077     for (; bufferPos < myMaxNumValPerParam; ++bufferPos) {
0078       if (!this->readIfSameLine<float>(myPedeResult, buffer[bufferPos]))
0079         break;
0080     }
0081
0082     // First check whether parameter is from any calibration (to be done before RunRange check:
0083     // run dependence for them is not handled here, but probably inside the calibration).
0084     // Double setting by calling read(..) twice for different RunRanges shouldn't harm.
0085     std::pair<IntegratedCalibrationBase *, unsigned int> calibParam = myLabels.calibrationParamFromLabel(paramLabel);
0086     if (calibParam.first) {  // label belongs to a calibration
0087       if (this->setCalibrationParameter(calibParam.first, calibParam.second, bufferPos, buffer)) {
0088         ++nParamCalib;
0089       } else {
0090         edm::LogError("Alignment") << "@SUB=PedeReader::read"
0091                                    << "Problems setting results of "
0092                                    << "parameter " << calibParam.second << " to calibration '"
0093                                    << calibParam.first->name() << "' (" << calibParam.first << ").";
0094         isAllOk = false;
0095       }
0096       continue;  // does not belong to any Alignable, so go to next line of file
0097     }
0098     // Now treat Alignables if paramLabel fits to run range, otherwise skip line:
0099     const RunRange &runRange = myLabels.runRangeFromLabel(paramLabel);
0100     if (!(runRange.first <= myRunRange.first && myRunRange.second <= runRange.second))
0101       continue;
0102
0103     Alignable *alignable = this->setParameter(paramLabel, bufferPos, buffer, setUserVars);
0104     if (!alignable) {
0105       // OK, e.g. for PedeReaderInputs: calibration parameters, but not yet known to labeler
0106       ++nParamUnknown;
0107       continue;
0108     }
0109     uniqueList.insert(alignable);
0110     ++nParam;
0111   }
0112
0113   // add Alignables to output
0114   alignables.insert(alignables.end(), uniqueList.begin(), uniqueList.end());
0115
0116   std::stringstream out;  // "@SUB=PedeReader::read" cannot yet go to 'out' for proper format
0117   out << nParam << " parameters for " << alignables.size() << " alignables and " << nParamCalib << " for calibrations, "
0118       << nParamUnknown << " parameters are unknown.\n";
0119   if (nParamUnknown) {
0120     edm::LogWarning("Alignment") << "@SUB=PedeReader::read" << out.str();
0121   } else {
0122     edm::LogInfo("Alignment") << "@SUB=PedeReader::read" << out.str();
0123   }
0124
0125   return isAllOk && (nParam + nParamCalib);  // nParam+nParamCalib == 0: empty or bad file
0126 }
0127
0128 //__________________________________________________________________________________________________
0129 template <class T>
0130 bool PedeReader::readIfSameLine(std::ifstream &aStream, T &outValue) const {
0131   while (true) {
0132     const int aChar = aStream.get();
0133     if (!aStream.good())
0134       return false;
0135
0136     switch (aChar) {
0137       case ' ':
0138       case '\t':
0139         continue;  // to next character
0140       case '\n':
0141         return false;  // end of line
0142       default:
0143         aStream.unget();
0144         aStream >> outValue;
0145         if (aStream.fail()) {  // not correct type 'T' (!aStream.good() is true also in case of EOF)
0146           aStream.clear();
0147           while (aStream.good() && aStream.get() != '\n')
0148             ;  // forward to end of line
0149           return false;
0150         } else {
0151           return true;
0152         }
0153     }  // switch
0154   }    // while
0155
0156   edm::LogError("Alignment") << "@SUB=PedeReader::readIfSameLine"
0157                              << "Should never come here!";
0158   return false;
0159 }
0160
0161 //__________________________________________________________________________________________________
0162 Alignable *PedeReader::setParameter(unsigned int paramLabel,
0163                                     unsigned int bufLength,
0164                                     const float *buf,
0165                                     bool setUserVars) const {
0166   Alignable *alignable = myLabels.alignableFromLabel(paramLabel);
0167   const unsigned int paramNum = myLabels.paramNumFromLabel(paramLabel);
0168   const double cmsToPede = mySteerer.cmsToPedeFactor(paramNum);
0169   if (alignable) {
0170     AlignmentParameters *params = this->checkAliParams(alignable, setUserVars);
0171     MillePedeVariables *userParams =  // static cast ensured by previous checkAliParams
0172         (setUserVars ? static_cast<MillePedeVariables *>(params->userVariables()) : nullptr);
0173     // if (userParams && userParams->label() != myLabels.alignableLabelFromLabel(paramLabel)) {
0174     if (userParams && userParams->label() != myLabels.alignableLabel(alignable)) {
0175       edm::LogError("Alignment") << "@SUB=PedeReader::setParameter"
0176                                  << "Label mismatch: paramLabel " << paramLabel << " for alignableLabel "
0177                                  << userParams->label();
0178     }
0179
0180     AlgebraicVector parVec(params->parameters());
0181     AlgebraicSymMatrix covMat(params->covariance());
0182
0183     if (userParams)
0184       userParams->setAllDefault(paramNum);
0185
0186     switch (bufLength) {
0187       case 5:  // global correlation
0188         if (userParams)
0189           userParams->globalCor()[paramNum] = buf[4];
0190         [[fallthrough]];
0191       case 4:  // uncertainty
0192         if (userParams)
0193           userParams->sigma()[paramNum] = buf[3] / cmsToPede;
0194         covMat[paramNum][paramNum] = buf[3] * buf[3] / (cmsToPede * cmsToPede);
0195         [[fallthrough]];
0196       case 3:  // difference to start value
0197         if (userParams)
0198           userParams->diffBefore()[paramNum] = buf[2] / cmsToPede;
0199         [[fallthrough]];
0200       case 2:
0201         params->setValid(true);
0202         parVec[paramNum] = buf[0] / cmsToPede * mySteerer.parameterSign();  // parameter
0203         if (userParams) {
0204           userParams->parameter()[paramNum] = parVec[paramNum];  // duplicate in millepede parameters
0205           userParams->preSigma()[paramNum] = buf[1];             // presigma given, probably means fixed
0206           if (!userParams->isFixed(paramNum)) {
0207             userParams->preSigma()[paramNum] /= cmsToPede;
0208             if (bufLength == 2) {
0209               edm::LogWarning("Alignment")
0210                   << "@SUB=PedeReader::setParameter"
0211                   << "Param " << paramLabel << " (from " << typeid(*alignable).name() << ") without result!";
0212               userParams->isValid()[paramNum] = false;
0213               params->setValid(false);
0214             }
0215           }
0216         }
0217         break;
0218       case 0:
0219       case 1:
0220       default:
0221         edm::LogError("Alignment") << "@SUB=PedeReader::setParameter"
0222                                    << "Expect 2 to 5 values, got " << bufLength << " for label " << paramLabel;
0223         break;
0224     }
0225     alignable->setAlignmentParameters(params->clone(parVec, covMat));  //transferred mem. responsib.
0226   } else {
0227     unsigned int lasBeamId = myLabels.lasBeamIdFromLabel(paramLabel);
0228     edm::LogError("Alignment") << "@SUB=PedeReader::setParameter"
0229                                << "No alignable for paramLabel " << paramLabel << ", probably LasBeam with Id "
0230                                << lasBeamId << ",\nparam " << paramNum << ": "
0231                                << buf[0] / cmsToPede * mySteerer.parameterSign()
0232                                << " += " << (bufLength >= 4 ? buf[3] / cmsToPede : -99.);
0233   }
0234
0235   return alignable;
0236 }
0237
0238 //__________________________________________________________________________________________________
0239 bool PedeReader::setCalibrationParameter(IntegratedCalibrationBase *calib,
0240                                          unsigned int paramNum,
0241                                          unsigned int bufLength,
0242                                          const float *buf) const {
0243   if (!calib || !buf)
0244     return false;
0245
0246   // FIXME: Should we attach MillePedeVariables to IntegratedCalibrationBase to store
0247   //        'other' results beyond value and error?
0248   switch (bufLength) {
0249     case 5:  // buf[4]: global correlation - not treated yet FIXME // no break;
0250     case 4:  // uncertainty
0251       calib->setParameterError(paramNum, buf[3]);
0252       [[fallthrough]];
0253     case 3:  // buf[2]: difference to start value - not treated yet
0254     case 2:
0255       if (bufLength == 2 && buf[1] >= 0.) {  // buf[1]: pre-sigma, < 0 means fixed
0256         edm::LogWarning("Alignment") << "@SUB=PedeReader::setCalibrationParameter"
0257                                      << "Param " << paramNum << " of calibration '" << calib->name()
0258                                      << "' without result!";
0259       }
0260       return calib->setParameter(paramNum, buf[0] * mySteerer.parameterSign());
0261     case 0:
0262     case 1:
0263     default:
0264       edm::LogError("Alignment") << "@SUB=PedeReader::setCalibrationParameter"
0265                                  << "Expect 2 to 5 values, got " << bufLength << ".";
0266       return false;
0267   }
0268 }
0269
0270 //__________________________________________________________________________________________________
0271 AlignmentParameters *PedeReader::checkAliParams(Alignable *alignable, bool createUserVars) const {
0272   // first check that we have parameters
0273   AlignmentParameters *params = alignable->alignmentParameters();
0274   if (!params) {
0275     throw cms::Exception("BadConfig") << "PedeReader::checkAliParams: "
0276                                       << "Alignable without parameters.";
0277   }
0278
0279   // now check that we have user parameters of correct type if requested:
0280   if (createUserVars && !dynamic_cast<MillePedeVariables *>(params->userVariables())) {
0281     edm::LogInfo("Alignment") << "@SUB=PedeReader::checkAliParams"
0282                               << "Add user variables for alignable with label " << myLabels.alignableLabel(alignable);
0283     params->setUserVariables(new MillePedeVariables(
0284         params->size(),
0285         myLabels.alignableLabel(alignable),
0286         myLabels.alignableTracker()->objectIdProvider().typeToName(alignable->alignableObjectId())));
0287   }
0288
0289   return params;
0290 }