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