PedeSteerer.cc

Go to the documentation of this file.

#include "PedeSteerer.h"
#include "PedeSteererWeakModeConstraints.h"

#include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Alignment/CommonAlignment/interface/Alignable.h"
#include "Alignment/CommonAlignment/interface/Utilities.h"
#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterSelector.h"
#include "Alignment/CommonAlignmentAlgorithm/interface/SelectionUserVariables.h"
#include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
#include "Alignment/CommonAlignmentParametrization/interface/RigidBodyAlignmentParameters.h"

#include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
// for 'type identification' as Alignable
#include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
#include "Alignment/MuonAlignment/interface/AlignableMuon.h"
#include "Alignment/CommonAlignment/interface/AlignableExtras.h"
// GF doubts the need of these includes from include checker campaign:
#include <FWCore/Framework/interface/EventSetup.h>
#include <Geometry/CommonDetUnit/interface/GeomDet.h>
#include <Geometry/CommonDetUnit/interface/GeomDetType.h>
#include <DataFormats/GeometrySurface/interface/LocalError.h>
#include <Geometry/DTGeometry/interface/DTLayer.h>
// end of doubt

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include <atomic>
#include <fstream>
#include <sstream>
#include <algorithm>

// from ROOT
#include <TSystem.h>
#include <TMath.h>

#include <iostream>

//_________________________________________________________________________
PedeSteerer::PedeSteerer(AlignableTracker *aliTracker,
                         AlignableMuon *aliMuon,
                         AlignableExtras *aliExtras,
                         AlignmentParameterStore *store,
                         const PedeLabelerBase *labels,
                         const edm::ParameterSet &config,
                         const std::string &defaultDir,
                         bool noSteerFiles)
    : myParameterStore(store),
      myLabels(labels),
      alignableObjectId_{AlignableObjectId::commonObjectIdProvider(aliTracker, aliMuon)},
      myConfig(config),
      myDirectory(myConfig.getUntrackedParameter<std::string>("fileDir")),
      myRunDirectory(myConfig.getUntrackedParameter<std::string>("runDir")),
      myNoSteerFiles(noSteerFiles),
      myIsSteerFileDebug(myConfig.getUntrackedParameter<bool>("steerFileDebug")),
      myParameterSign(myConfig.getUntrackedParameter<int>("parameterSign")),
      theMinHieraConstrCoeff(myConfig.getParameter<double>("minHieraConstrCoeff")),
      theMinHieraParPerConstr(myConfig.getParameter<unsigned int>("minHieraParPerConstr")),
      theConstrPrecision(myConfig.getParameter<unsigned int>("constrPrecision")),
      theCoordMaster(nullptr) {
  if (myParameterSign != 1 && myParameterSign != -1) {
    cms::Exception("BadConfig") << "Expect PedeSteerer.parameterSign = +/-1, "
                                << "found " << myParameterSign << ".";
  }

  // Correct directory, needed before asking for fileName(..):
  if (myDirectory.empty())
    myDirectory = defaultDir;
  if (!myDirectory.empty() && myDirectory.find_last_of('/') != myDirectory.size() - 1) {
    myDirectory += '/';  // directory may need '/'
  }

  if (myRunDirectory.empty())
    myRunDirectory = defaultDir;
  if (!myRunDirectory.empty() && myRunDirectory.find_last_of('/') != myRunDirectory.size() - 1) {
    myRunDirectory += '/';  // directory may need '/'
    myDirectory = myRunDirectory;
  }

  const auto &alis = myParameterStore->alignables();
  if (!this->checkParameterChoices(alis)) {
  }  // anyway thrown exception

  // Coordinate system selection and correction before everything
  theCoordDefiners = this->selectCoordinateAlis(alis);
  if (!theCoordDefiners.empty()) {  // Create steering with constraints to define coordinate system:
    // OK, some hacks:
    // - we want a composite with global coordinates where tracker and muon are components
    //   (to call RigidBodyAl.Param.->globalParameters() in correctToReferenceSystem(..))
    // - so we create a AlignableComposite and add tracker and muon
    // - but the addComponent(..) method is so cute that it calculates position from
    //   daughters' deepComponents()
    // - so we want to move it back to (0,0,0), but ali->move(..) would move daughters as well
    //   => move only the surface back
    // - this master object does not have a label for its parameters
    //   => some warnings if debug output selected in pedeSteer files
    // - we must not delete our new master (little mem. leak...) since that would delete
    //   the daughters as well!
    if (aliTracker) {
      theCoordMaster = new AlignableComposite(aliTracker->id(), align::invalid);
      theCoordMaster->addComponent(aliTracker);
    } else if (aliMuon) {
      theCoordMaster = new AlignableComposite(aliMuon->id(), align::invalid);
    } else {
      throw cms::Exception("BadConfig") << "[PedeSteerer]"
                                        << "Cannot define global coordinate system "
                                        << "with neither tracker nor muon!";
    }
    if (aliMuon)
      theCoordMaster->addComponent(aliMuon);  // tracker is already added if existing
    if (aliExtras) {                          // tracker and/or muon are already added if existing
      align::Alignables allExtras = aliExtras->components();
      for (const auto &it : allExtras)
        theCoordMaster->addComponent(it);
    }

    theCoordMaster->recenterSurface();

    if (this->isCorrectToRefSystem(theCoordDefiners)) {  // defined by 's' (MC): 'correct' misalignment
      this->correctToReferenceSystem();                  // really before 'defineCoordinates'?
    }
  }
}

//___________________________________________________________________________
PedeSteerer::~PedeSteerer() {
  // delete theCoordMaster; NO, see above
}

//_________________________________________________________________________
bool PedeSteerer::isNoHiera(const Alignable *ali) const {
  return (myNoHieraCollection.find(ali) != myNoHieraCollection.end());
}

//_________________________________________________________________________
double PedeSteerer::cmsToPedeFactor(unsigned int parNum) const {
  return 1.;  // mmh, otherwise would need to FIXME hierarchyConstraint...

  switch (parNum) {
    case RigidBodyAlignmentParameters::dx:
    case RigidBodyAlignmentParameters::dy:
      return 1000.;  // cm to mum *1/10 to get smaller values
    case RigidBodyAlignmentParameters::dz:
      return 2500.;  // cm to mum *1/4
    case RigidBodyAlignmentParameters::dalpha:
    case RigidBodyAlignmentParameters::dbeta:
      return 1000.;  // rad to mrad (no first guess for sensitivity yet)
    case RigidBodyAlignmentParameters::dgamma:
      return 10000.;  // rad to mrad *10 to get larger values
    default:
      return 1.;
  }
}

//_________________________________________________________________________
unsigned int PedeSteerer::buildNoHierarchyCollection(const align::Alignables &alis) {
  myNoHieraCollection.clear();  // just in case of re-use...

  for (const auto &iAli : alis) {
    AlignmentParameters *params = iAli->alignmentParameters();
    SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(params->userVariables());
    if (!selVar)
      continue;
    // Now check whether taking out of hierarchy is selected - must be consistent!
    unsigned int numNoHieraPar = 0;
    unsigned int numHieraPar = 0;
    for (unsigned int iParam = 0; static_cast<int>(iParam) < params->size(); ++iParam) {
      const char selector = selVar->fullSelection()[iParam];
      if (selector == 'C' || selector == 'F' || selector == 'H') {
        ++numNoHieraPar;
      } else if (selector == 'c' || selector == 'f' || selector == '1' || selector == 'r' || selector == 's') {
        ++numHieraPar;
      }  // else ... accept '0' as undetermined
    }
    if (numNoHieraPar) {  // Selected to be taken out.
      if (numHieraPar) {  // Inconsistent: Some parameters still in hierarchy ==> exception!
        throw cms::Exception("BadConfig")
            << "[PedeSteerer::buildNoHierarchyCollection] All active parameters of alignables to be "
            << " taken out of the hierarchy must be marked with capital letters 'C', 'F' or 'H'!";
      }
      bool isInHiera = false;  // Check whether Alignable is really part of hierarchy:
      auto mother = iAli;
      while ((mother = mother->mother())) {
        if (mother->alignmentParameters())
          isInHiera = true;  // could 'break;', but loop is short
      }
      // Complain, but keep collection short if not in hierarchy:
      if (isInHiera)
        myNoHieraCollection.insert(iAli);
      else
        edm::LogWarning("Alignment") << "@SUB=PedeSteerer::buildNoHierarchyCollection"
                                     << "Alignable not in hierarchy, no need to remove it!";
    }
  }  // end loop on alignables

  return myNoHieraCollection.size();
}

//_________________________________________________________________________
bool PedeSteerer::checkParameterChoices(const align::Alignables &alis) const {
  for (const auto &iAli : alis) {
    AlignmentParameters *paras = iAli->alignmentParameters();
    SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(paras->userVariables());
    if (!selVar)
      continue;
    for (unsigned int iParam = 0; static_cast<int>(iParam) < paras->size(); ++iParam) {
      const char sel = selVar->fullSelection()[iParam];
      if (sel != 'f' && sel != 'F' && sel != 'c' && sel != 'C' && sel != '0' && sel != '1' && sel != 'H' &&
          sel != 'r' && sel != 's') {
        throw cms::Exception("BadConfig")
            << "[PedeSteerer::unknownParameterChoices] "
            << "Unexpected parameter selector '" << sel
            << "', use \n'f/F' (fix),\n'c/C' (fix at correct pos.),\n'1/H' (free),\n"
            << "'r/s' (free, but defining reference system, trying to correct misalignment if 's')"
            << " or \n'0' (ignore).\n"
            << "Capital letters mean that the Alignable is taken out of a possible hierarchy,\n"
            << "but must be used consistently for all its parameters.";
        return false;  // unreached
      }
    }
  }

  return true;
}

//_________________________________________________________________________
std::pair<unsigned int, unsigned int> PedeSteerer::fixParameters(const align::Alignables &alis,
                                                                 const std::string &fileName) {
  // return number of parameters fixed at 0. and fixed at original position
  std::pair<unsigned int, unsigned int> numFixNumFixCor(0, 0);

  std::ofstream *filePtr = nullptr;

  for (const auto &iAli : alis) {
    AlignmentParameters *params = iAli->alignmentParameters();
    SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(params->userVariables());
    if (!selVar)
      continue;

    for (unsigned int iParam = 0; static_cast<int>(iParam) < params->size(); ++iParam) {
      const unsigned int nInstances = myLabels->numberOfParameterInstances(iAli, iParam);
      for (unsigned int iInstance = 0; iInstance < nInstances; ++iInstance) {
        int whichFix = this->fixParameter(iAli, iInstance, iParam, selVar->fullSelection()[iParam], filePtr, fileName);
        if (whichFix == 1) {
          ++(numFixNumFixCor.first);
        } else if (whichFix == -1) {
          ++(numFixNumFixCor.second);
        }
      }
    }
  }

  delete filePtr;  // automatically flushes, no problem if NULL ptr.

  return numFixNumFixCor;
}

//_________________________________________________________________________
int PedeSteerer::fixParameter(Alignable *ali,
                              unsigned int iInstance,
                              unsigned int iParam,
                              char selector,
                              std::ofstream *&filePtr,
                              const std::string &fileName) {
  int result = 0;
  float fixAt = 0.;
  if (selector == 'c' || selector == 'C') {
    if (ali->alignmentParameters()->type() != AlignmentParametersFactory::kRigidBody) {
      throw cms::Exception("BadConfig")
          << "PedeSteerer::fixParameter: correction (c/C) possible only for RigidBodyParameters";
    }
    fixAt = -this->parameterSign() * RigidBodyAlignmentParameters(ali, true).parameters()[iParam];
    result = -1;
  } else if (selector == 'f' || selector == 'F') {
    result = 1;
  }

  if (result) {
    if (!filePtr) {
      filePtr = this->createSteerFile(fileName, true);
      (*filePtr) << "Parameter\n";
    }
    std::ofstream &file = *filePtr;

    const unsigned int aliLabel = myLabels->alignableLabelFromParamAndInstance(ali, iParam, iInstance);
    file << myLabels->parameterLabel(aliLabel, iParam) << "  " << fixAt * this->cmsToPedeFactor(iParam) << " -1.0";
    if (myIsSteerFileDebug) {  // debug
      const GlobalPoint position(ali->globalPosition());
      file << "  * id " << ali->id() << ", eta " << position.eta() << ", z " << position.z() << ", r "
           << position.perp() << ", phi " << position.phi();
    }
    file << "\n";
  }

  return result;
}

//_________________________________________________________________________
align::Alignables PedeSteerer::selectCoordinateAlis(const align::Alignables &alis) const {
  align::Alignables coordAlis;

  for (const auto &iAli : alis) {
    AlignmentParameters *params = iAli->alignmentParameters();
    SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(params->userVariables());
    if (!selVar)
      continue;
    unsigned int refParam = 0;
    unsigned int nonRefParam = 0;
    for (unsigned int iParam = 0; static_cast<int>(iParam) < params->size(); ++iParam) {
      const char selector = selVar->fullSelection()[iParam];
      if (selector == 'r' || selector == 's') {
        ++refParam;
      } else if (selector != '0' && selector != 'f') {  // allow also 'c'?
        ++nonRefParam;
      }
    }
    // Check whether some 'r/s' selection string. If yes and selection makes sense, add to result:
    if (refParam) {
      if (nonRefParam) {
        throw cms::Exception("BadConfig")
            << "[PedeSteerer::selectCoordinateAlis] All active parameters of alignables defining "
            << "the coordinate system must be marked with 'r/s' (or fixed, 'f')!";
      } else {
        auto mother = iAli;
        while ((mother = mother->mother())) {
          if (mother->alignmentParameters()) {
            throw cms::Exception("BadConfig") << "[PedeSteerer::selectCoordinateAlis] "
                                              << "Alignables defining the coordinate system must "
                                              << "be highest level!";
          }
        }
        coordAlis.push_back(iAli);
      }
    }
  }  // end loop on alignables

  return coordAlis;
}

//_________________________________________________________________________
void PedeSteerer::defineCoordinates(const align::Alignables &alis, Alignable *aliMaster, const std::string &fileName) {
  std::ofstream *filePtr = this->createSteerFile(fileName, true);
  (*filePtr) << "* Constraints to define coordinate system:\n";
  if (!aliMaster || aliMaster->alignmentParameters()) {
    throw cms::Exception("BadConfig") << "[PedeSteerer::defineCoordinates] "
                                      << "No master alignable or it has parameters!";
  }
  if (myIsSteerFileDebug) {  // See constructor comments about hack:
    edm::LogError("Alignment") << "@SUB=PedeSteerer::defineCoordinates"
                               << "Ignore following LogicErrors from PedeLabeler.";
  }
  AlignmentParameters *par = new RigidBodyAlignmentParameters(aliMaster, false);
  aliMaster->setAlignmentParameters(par);  // hierarchyConstraint needs parameters
  this->hierarchyConstraint(aliMaster, alis, *filePtr);
  aliMaster->setAlignmentParameters(nullptr);  // erase dummy parameters

  delete filePtr;  // automatically flushes, no problem if NULL ptr.
}

//_________________________________________________________________________
bool PedeSteerer::isCorrectToRefSystem(const align::Alignables &coordDefiners) const {
  bool doCorrect = false;
  bool doNotCorrect = false;
  for (const auto &it : coordDefiners) {
    SelectionUserVariables *selVar =
        (it->alignmentParameters() ? dynamic_cast<SelectionUserVariables *>(it->alignmentParameters()->userVariables())
                                   : nullptr);
    if (!selVar)
      continue;  // is an error!?

    for (unsigned int i = 0; i < selVar->fullSelection().size(); ++i) {
      if (selVar->fullSelection()[i] == 'r')
        doNotCorrect = true;
      else if (selVar->fullSelection()[i] == 's')
        doCorrect = true;
    }
  }

  if (doCorrect && doNotCorrect) {
    throw cms::Exception("BadConfig")
        << "[PedeSteerer::doCorrectToRefSystem]: Parameter selection 's' and 'r' must not coexist!";
  }

  return doCorrect;
}

//_________________________________________________________________________
void PedeSteerer::correctToReferenceSystem() {
  typedef RigidBodyAlignmentParameters RbPars;
  if (!theCoordMaster || theCoordDefiners.empty())
    return;  // nothing was defined

  align::Alignables definerDets;  // or ...DetUnits
  for (const auto &it :
       theCoordDefiners) {  // find lowest level objects of alignables that define the coordinate system
    const auto &comp = it->deepComponents();
    definerDets.insert(definerDets.end(), comp.begin(), comp.end());
  }

  for (unsigned int iLoop = 0;; ++iLoop) {  // iterate: shifts and rotations are not independent
    AlgebraicVector meanPars(RbPars::N_PARAM);
    for (const auto &it : definerDets) {                // sum up mean displacements/misrotations:
      meanPars += RbPars(it, true).globalParameters();  // requires theCoordMaster has global frame
    }
    meanPars /= definerDets.size();
    const align::Scalar squareSum = meanPars.normsq();

    if (squareSum < 1.e-20)
      break;  // sqrt(1.e-20)=1.e-10: close enough to stop iterating
    if (iLoop == 0) {
      edm::LogInfo("Alignment") << "@SUB=PedeSteerer::correctToReferenceSystem"
                                << "Loop " << iLoop << " "
                                << "Mean misalignment of dets of defined coordinate system"
                                << (squareSum < 1.e-20 ? ":" : " (will be iteratively corrected to < 1.e-10):")
                                << meanPars;
    }
    if (iLoop >= 5) {  // 3 iterations should be safe, use 5 for 'more' safety...
      edm::LogError("Alignment") << "@SUB=PedeSteerer::correctToReferenceSystem"
                                 << "No convergence in " << iLoop << " iterations, "
                                 << "remaining misalignment: " << meanPars;
      break;
    }

    const GlobalVector globalShift(meanPars[RbPars::dx], meanPars[RbPars::dy], meanPars[RbPars::dz]);
    theCoordMaster->move(-globalShift);  // sign to revert
    align::EulerAngles globalAngles(3);
    globalAngles[0] = meanPars[RbPars::dalpha];
    globalAngles[1] = meanPars[RbPars::dbeta];
    globalAngles[2] = meanPars[RbPars::dgamma];
    theCoordMaster->rotateInGlobalFrame(align::toMatrix(-globalAngles));  // sign to revert
  }
}

//_________________________________________________________________________
unsigned int PedeSteerer::hierarchyConstraints(const align::Alignables &alis, const std::string &fileName) {
  std::ofstream *filePtr = nullptr;

  unsigned int nConstraints = 0;
  align::Alignables aliDaughts;
  for (const auto &iA : alis) {
    aliDaughts.clear();
    if (!(iA->firstCompsWithParams(aliDaughts))) {
      edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraints"
                                   << "Some but not all daughters of "
                                   << alignableObjectId_.idToString(iA->alignableObjectId()) << " with params!";
    }
    //     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::hierarchyConstraints"
    //                << aliDaughts.size() << " ali param components";
    if (aliDaughts.empty())
      continue;
    //     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::hierarchyConstraints"
    //                << aliDaughts.size() << " alignable components ("
    //                << iA->size() << " in total) for "
    //                << aliId.alignableTypeName(iA)
    //                << ", layer " << aliId.typeAndLayerFromAlignable(iA).second
    //                << ", position " << (iA)->globalPosition()
    //                << ", r = " << (iA)->globalPosition().perp();
    if (!filePtr)
      filePtr = this->createSteerFile(fileName, true);
    ++nConstraints;
    this->hierarchyConstraint(iA, aliDaughts, *filePtr);
  }

  delete filePtr;  // automatically flushes, no problem if NULL ptr.

  return nConstraints;
}
//_________________________________________________________________________
void PedeSteerer::hierarchyConstraint(const Alignable *ali,
                                      const align::Alignables &components,
                                      std::ofstream &file) const {
  typedef AlignmentParameterStore::ParameterId ParameterId;

  std::vector<std::vector<ParameterId> > paramIdsVec;
  std::vector<std::vector<double> > factorsVec;
  const bool allConstr = false;  // true; // make configurable?
  static std::atomic<bool> allConstrWarning{false};
  bool expected{false};
  if (allConstr && allConstrWarning.compare_exchange_strong(expected, true)) {
    edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraint"
                                 << "changed to use all 6 constraints";
  }
  if (!myParameterStore->hierarchyConstraints(
          ali, components, paramIdsVec, factorsVec, allConstr, theMinHieraConstrCoeff)) {
    edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraint"
                                 << "Problems from store.";
  }

  for (unsigned int iConstr = 0; iConstr < paramIdsVec.size(); ++iConstr) {
    std::ostringstream aConstr;

    const std::vector<ParameterId> &parIds = paramIdsVec[iConstr];
    const std::vector<double> &factors = factorsVec[iConstr];
    unsigned int nParPerConstr = 0;  // keep track of used factor/parId pair
    // parIds.size() == factors.size() granted by myParameterStore->hierarchyConstraints
    for (unsigned int iParam = 0; iParam < parIds.size(); ++iParam) {
      Alignable *aliSubComp = parIds[iParam].first;
      const unsigned int compParNum = parIds[iParam].second;
      if (this->isNoHiera(aliSubComp)) {
        if (myIsSteerFileDebug)
          aConstr << "* Taken out of hierarchy: ";
        continue;
      }
      const unsigned int aliLabel = myLabels->alignableLabel(aliSubComp);
      const unsigned int paramLabel = myLabels->parameterLabel(aliLabel, compParNum);
      // FIXME: multiply by cmsToPedeFactor(subcomponent)/cmsToPedeFactor(mother) (or vice a versa?)
      if (theConstrPrecision > 0)
        aConstr << paramLabel << "    " << std::setprecision(theConstrPrecision) << factors[iParam];
      else
        aConstr << paramLabel << "    " << factors[iParam];
      if (myIsSteerFileDebug) {  // debug
        aConstr << "   ! for param " << compParNum << " of a "
                << alignableObjectId_.idToString(aliSubComp->alignableObjectId()) << " at "
                << aliSubComp->globalPosition() << ", r=" << aliSubComp->globalPosition().perp();
      }
      aConstr << "\n";
      ++nParPerConstr;  // OK, we used one.
    }                   // end loop on params

    //
    if (nParPerConstr && nParPerConstr >= theMinHieraParPerConstr) {  // Enough to make sense?
      if (myIsSteerFileDebug) {                                       //debug
        file << "\n* Nr. " << iConstr << " of a '" << alignableObjectId_.idToString(ali->alignableObjectId())
             << "' (label " << myLabels->alignableLabel(const_cast<Alignable *>(ali))  // ugly cast: FIXME!
             << "), position " << ali->globalPosition() << ", r = " << ali->globalPosition().perp();
      }
      file << "\nConstraint   0.\n" << aConstr.str();  // in future 'Wconstraint'?
    } else if (nParPerConstr > 0) {                    // no warning for trivial case...
      edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraint"
                                   << "Skip constraint on " << nParPerConstr << " parameter(s):\n"
                                   << aConstr.str();
    }
  }  // end loop on constraints
}

//_________________________________________________________________________
unsigned int PedeSteerer::presigmas(const std::vector<edm::ParameterSet> &cffPresi,
                                    const std::string &fileName,
                                    const align::Alignables &alis,
                                    AlignableTracker *aliTracker,
                                    AlignableMuon *aliMuon,
                                    AlignableExtras *aliExtras) {
  // We loop on given PSet's, each containing a parameter selection and the presigma value
  // The resulting presigmas are stored in a map with Alignable* as key.
  // This map, 'fileName' and 'alis' are passed further to create the steering file.

  AlignmentParameterSelector selector(aliTracker, aliMuon, aliExtras);
  AlignablePresigmasMap aliPresiMap;  // map to store alis with presigmas of their parameters
  for (std::vector<edm::ParameterSet>::const_iterator iSet = cffPresi.begin(), iE = cffPresi.end(); iSet != iE;
       ++iSet) {  // loop on individual PSets defining ali-params with their presigma
    selector.clear();
    selector.addSelections((*iSet).getParameter<edm::ParameterSet>("Selector"));
    const auto &alis = selector.selectedAlignables();
    const std::vector<std::vector<char> > &sels = selector.selectedParameters();
    const float presigma = (*iSet).getParameter<double>("presigma");
    if (presigma <= 0.) {  // given presigma > 0., 0. later used if not (yet) chosen for parameter
      throw cms::Exception("BadConfig") << "[PedeSteerer::presigmas]: Pre-sigma must be > 0., but is " << presigma
                                        << ".";
    }
    // now loop on alis of present selection
    for (unsigned int iAli = 0; iAli < alis.size(); ++iAli) {
      std::vector<float> &presigmas = aliPresiMap[alis[iAli]];  // existing or empty, so ensure length:
      if (presigmas.size() < sels[iAli].size())
        presigmas.resize(sels[iAli].size(), 0.);
      for (unsigned int iParam = 0; iParam < sels[iAli].size(); ++iParam) {  // loop on parameters
        if (sels[iAli][iParam] != '0') {  // all but '0' means to apply the chosen presigma
          if (presigmas[iParam] != 0.) {  // reset forbidden (would make it order dependent!)
            throw cms::Exception("BadConfig")
                << "[PedeSteerer::presigmas]: Try to set pre-sigma " << presigma << ", but already "
                << "set " << presigmas[iParam] << " (for a "
                << alignableObjectId_.idToString(alis[iAli]->alignableObjectId()) << ").";
          }
          presigmas[iParam] = presigma;
        }  // end if selected for presigma
      }    // end loop on params
    }      // end loop on alignables for given selection and presigma
  }        // end loop on PSets

  if (aliPresiMap.empty())
    return 0;
  else
    return this->presigmasFile(fileName, alis, aliPresiMap);
}

//_________________________________________________________________________
unsigned int PedeSteerer::presigmasFile(const std::string &fileName,
                                        const align::Alignables &alis,
                                        const AlignablePresigmasMap &aliPresiMap) {
  // Check if 'alis' are in aliPresiMap,
  // if yes apply presigma - but NOT if parameter is fixed!
  std::ofstream *filePtr = nullptr;

  unsigned int nPresiParam = 0;
  for (const auto &iAli : alis) {
    // Any presigma chosen for alignable?
    AlignablePresigmasMap::const_iterator presigmasIt = aliPresiMap.find(iAli);
    if (presigmasIt == aliPresiMap.end())
      continue;  // no presigma chosen for alignable

    // Why does the following not work? It does with CMSSW_1_3_X on SLC3...
    // const AlignablePresigmasMap::data_type &presigmas = presigmasIt->second;
    const std::vector<float> &presigmas = presigmasIt->second;  // I want to hide float or double...
    for (unsigned int iParam = 0; iParam < presigmas.size(); ++iParam) {
      // Now check whether a presigma value > 0. chosen:
      if (presigmas[iParam] <= 0.)
        continue;  // must be positive, '<' checked above
      // Do not apply presigma to inactive or fixed values.
      if (!(iAli->alignmentParameters()->selector()[iParam]))
        continue;
      SelectionUserVariables *selVar =
          dynamic_cast<SelectionUserVariables *>(iAli->alignmentParameters()->userVariables());
      const char selChar = (selVar ? selVar->fullSelection()[iParam] : '1');
      if (selChar == 'f' || selChar == 'F' || selChar == 'c' || selChar == 'C')
        continue;
      // Finally create and write steering file:
      if (!filePtr) {
        filePtr = this->createSteerFile(fileName, true);
        (*filePtr) << "* Presigma values for active parameters: \nParameter\n";
      }
      const unsigned int aliLabel = myLabels->alignableLabel(iAli);
      (*filePtr) << myLabels->parameterLabel(aliLabel, iParam) << "   0.   "
                 << presigmas[iParam] * fabs(this->cmsToPedeFactor(iParam));
      if (myIsSteerFileDebug) {
        (*filePtr) << "  ! for a " << alignableObjectId_.idToString((iAli)->alignableObjectId());
      }
      (*filePtr) << '\n';

      ++nPresiParam;
    }  // end loop on parameters for alignables with chosen presigmas
  }    // end loop on alignables

  delete filePtr;  // close properly file
  return nPresiParam;
}

//_________________________________________________________________________
std::ofstream *PedeSteerer::createSteerFile(const std::string &name, bool addToList) {
  const std::string realName(myNoSteerFiles ? "/dev/null" : name.c_str());

  std::ofstream *result = new std::ofstream(realName.c_str(), std::ios::out);
  if (!result || !result->is_open()) {
    delete result;  // needed before exception in case just open failed
    throw cms::Exception("FileOpenProblem") << "[PedeSteerer::createSteerFile]"
                                            << "Could not open " << realName << " as output file.";
  } else if (addToList) {
    mySteeringFiles.push_back(realName);  // keep track
  }

  return result;
}

//_________________________________________________________________________
std::string PedeSteerer::fileName(const std::string &addendum) const {
  std::string name(myDirectory);
  name += myConfig.getParameter<std::string>("steerFile");
  name += addendum;
  name += ".txt";

  return name;
}

//___________________________________________________________________________
void PedeSteerer::buildSubSteer(AlignableTracker *aliTracker, AlignableMuon *aliMuon, AlignableExtras *aliExtras) {
  const auto &alis = myParameterStore->alignables();

  if (theCoordMaster && !theCoordDefiners.empty()) {
    const std::string nameCoordFile(this->fileName("Coord"));
    this->defineCoordinates(theCoordDefiners, theCoordMaster, nameCoordFile);
    edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer" << theCoordDefiners.size()
                              << " highest level objects define the "
                              << "coordinate system, steering file " << nameCoordFile << ".";
  }

  const std::string nameFixFile(this->fileName("FixPara"));
  const std::pair<unsigned int, unsigned int> nFixFixCor(this->fixParameters(alis, nameFixFile));
  if (nFixFixCor.first != 0 || nFixFixCor.second != 0) {
    edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer" << nFixFixCor.first << " parameters fixed at 0. and "
                              << nFixFixCor.second << " at 'original' position, "
                              << "steering file " << nameFixFile << ".";
  }

  if (this->buildNoHierarchyCollection(alis)) {  // before hierarchyConstraints(..)
    edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer" << myNoHieraCollection.size()
                              << " alignables taken out of hierarchy.";
  }

  const std::string nameHierarchyFile(this->fileName("Hierarchy"));
  unsigned int nConstraint = this->hierarchyConstraints(alis, nameHierarchyFile);
  if (nConstraint) {
    edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer"
                              << "Hierarchy constraints for " << nConstraint << " alignables, "
                              << "steering file " << nameHierarchyFile << ".";
  }

  //construct the systematic geometry deformations
  if (!(myConfig.getParameter<std::vector<edm::ParameterSet> >("constraints")).empty()) {
    PedeSteererWeakModeConstraints GeometryConstraints(
        aliTracker,
        myLabels,
        myConfig.getParameter<std::vector<edm::ParameterSet> >("constraints"),
        myConfig.getParameter<std::string>("steerFile"));

    //prepare the output files
    //Get the data structure in which the configuration data are stored.
    //The relation between the ostream* and the corresponding file name needs to be filled
    auto &ConstraintsConfigContainer = GeometryConstraints.getConfigData();

    //loop over all configured constraints
    for (auto &it : ConstraintsConfigContainer) {
      //each level has its own constraint which means the output is stored in a separate file
      for (const auto &ilevelsFilename : it.levelsFilenames_) {
        it.mapFileName_.insert(
            std::make_pair(ilevelsFilename.second, this->createSteerFile(ilevelsFilename.second, true)));
      }
    }

    unsigned int nGeometryConstraint = GeometryConstraints.constructConstraints(alis);
    if (nGeometryConstraint) {
      edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer"
                                << "Geometry constraints for " << nGeometryConstraint << " alignables.";
    }
  }

  const std::string namePresigmaFile(this->fileName("Presigma"));
  unsigned int nPresigma = this->presigmas(myConfig.getParameter<std::vector<edm::ParameterSet> >("Presigmas"),
                                           namePresigmaFile,
                                           alis,
                                           aliTracker,
                                           aliMuon,
                                           aliExtras);
  if (nPresigma) {
    edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer"
                              << "Presigma values set for " << nPresigma << " parameters, "
                              << "steering file " << namePresigmaFile << ".";
  }

  // Delete all SelectionUserVariables now? They will anyway be overwritten by MillePedeVariables...
}

//_________________________________________________________________________
std::string PedeSteerer::buildMasterSteer(const std::vector<std::string> &binaryFiles) {
  const std::string nameMasterSteer(this->fileName("Master"));
  std::ofstream *mainSteerPtr = this->createSteerFile(nameMasterSteer, false);
  if (!mainSteerPtr)
    return "";

  // add external steering files, if any
  std::vector<std::string> addfiles = myConfig.getParameter<std::vector<std::string> >("additionalSteerFiles");
  mySteeringFiles.insert(mySteeringFiles.end(), addfiles.begin(), addfiles.end());

  // add steering files to master steering file
  std::ofstream &mainSteerRef = *mainSteerPtr;
  for (unsigned int iFile = 0; iFile < mySteeringFiles.size(); ++iFile) {
    mainSteerRef << mySteeringFiles[iFile] << "\n";
  }

  // add binary files to master steering file
  mainSteerRef << "\nCfiles\n";
  for (unsigned int iFile = 0; iFile < binaryFiles.size(); ++iFile) {
    if (myRunDirectory.empty())
      mainSteerRef << binaryFiles[iFile] << "\n";
    else
      mainSteerRef << "../" + binaryFiles[iFile] << "\n";
  }

  // add method
  mainSteerRef << "\nmethod  " << myConfig.getParameter<std::string>("method") << "\n";

  // add further options
  const std::vector<std::string> opt(myConfig.getParameter<std::vector<std::string> >("options"));
  mainSteerRef << "\n* Outlier treatment and other options \n";
  for (unsigned int i = 0; i < opt.size(); ++i) {
    mainSteerRef << opt[i] << "\n";
  }

  delete mainSteerPtr;  // close (and flush) again

  return nameMasterSteer;
}

//_________________________________________________________________________
int PedeSteerer::runPede(const std::string &masterSteer) const {
  if (masterSteer.empty()) {
    edm::LogError("Alignment") << "@SUB=PedeSteerer::runPede"
                               << "Empty master steer file, stop";
    return 0;  //false;
  }

  // Change pede command if running in different directory
  std::string command;
  if (myRunDirectory.empty()) {
    command = myConfig.getUntrackedParameter<std::string>("pedeCommand");
    (command += " ") += masterSteer;
  } else {
    command = "(cd " + myRunDirectory + " && ";
    command += myConfig.getUntrackedParameter<std::string>("pedeCommand");
    (command += " ") += "../" + masterSteer;
    command += ")";
  }

  const std::string dump(myConfig.getUntrackedParameter<std::string>("pedeDump"));
  if (!dump.empty()) {
    command += " > ";
    (command += myDirectory) += dump;
  }

  edm::LogInfo("Alignment") << "@SUB=PedeSteerer::runPede"
                            << "Start running " << command;
  // FIXME: Recommended interface to system commands?
  int shellReturn = gSystem->Exec(command.c_str());
  edm::LogInfo("Alignment") << "@SUB=PedeSteerer::runPede"
                            << "Command returns " << shellReturn;

  return shellReturn;
}