PedeSteererWeakModeConstraints.cc

Go to the documentation of this file.

#include "Alignment/MillePedeAlignmentAlgorithm/src/PedeSteererWeakModeConstraints.h"
#include "Alignment/MillePedeAlignmentAlgorithm/src/PedeSteerer.h"

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

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

#include "Alignment/CommonAlignment/interface/Alignable.h"
#include <boost/cstdint.hpp>
#include "Alignment/CommonAlignment/interface/AlignableObjectId.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/CommonAlignmentParametrization/interface/BowedSurfaceAlignmentDerivatives.h"
#include "Alignment/CommonAlignmentParametrization/interface/TwoBowedSurfacesAlignmentParameters.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/GeomDetUnit.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 <fstream>
#include <sstream>
#include <algorithm>

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

#include <iostream>





GeometryConstraintConfigData::GeometryConstraintConfigData(const std::vector<double>& co,
                                                           const std::string& c,
                                                           const std::vector<std::pair<Alignable*,std::string> >& alisFile,
                                                           const int sd,
                                                           const std::vector<Alignable*>& ex,
                                                           const int instance,
                               const bool downToLowestLevel
                                                           ) :
  coefficients_(co),
  constraintName_(c),
  levelsFilenames_(alisFile),
  excludedAlignables_(ex),
  sysdeformation_(sd),
  instance_(instance),
  downToLowestLevel_(downToLowestLevel)
{
}

//_________________________________________________________________________
PedeSteererWeakModeConstraints::PedeSteererWeakModeConstraints(AlignableTracker *aliTracker,
                                                               const PedeLabelerBase *labels,
                                                               const std::vector<edm::ParameterSet> &config,
                                                               std::string sf
                                                               ) :
  myLabels_(labels),
  myConfig_(config),
  steerFile_(sf)
{
  unsigned int psetnr = 0;
  std::set<std::string> steerFilePrefixContainer;
  for(const auto& pset: myConfig_) {
    this->verifyParameterNames(pset, psetnr);
    psetnr++;

    const auto coefficients = pset.getParameter<std::vector<double> > ("coefficients");
    const auto dm = pset.exists("deadmodules") ?
      pset.getParameter<std::vector<unsigned int> >("deadmodules") : std::vector<unsigned int>();
    std::string name = pset.getParameter<std::string> ("constraint");
    std::transform(name.begin(), name.end(), name.begin(), ::tolower);
    const auto ignoredInstances = pset.exists("ignoredInstances") ?
      pset.getUntrackedParameter<std::vector<unsigned int> >("ignoredInstances"):
      std::vector<unsigned int>();

    const auto downToLowestLevel = pset.exists("downToLowestLevel") ?
      pset.getUntrackedParameter<bool>("downToLowestLevel"): false;

    AlignmentParameterSelector selector(aliTracker, nullptr, nullptr);
    selector.clear();
    selector.addSelections(pset.getParameter<edm::ParameterSet> ("levels"));

    const auto& alis = selector.selectedAlignables();

    AlignmentParameterSelector selector_excludedalignables(aliTracker, nullptr, nullptr);
    selector_excludedalignables.clear();
    if(pset.exists("excludedAlignables")) {
      selector_excludedalignables.addSelections(pset.getParameter<edm::ParameterSet> ("excludedAlignables"));
    }
    const auto& excluded_alis = selector_excludedalignables.selectedAlignables();

    //check that the name of the deformation is known and that the number
    //of provided parameter is right.
    auto sysdeformation = this->verifyDeformationName(name,coefficients);

    if(deadmodules_.size() == 0) { //fill the list of dead modules only once
      edm::LogInfo("Alignment") << "@SUB=PedeSteererWeakModeConstraints"
                                << "Load list of dead modules (size = " << dm.size()<< ").";
      for(const auto& it: dm) deadmodules_.push_back(it);
    }

    // loop over all IOVs/momentum ranges
    for (unsigned int instance = 0;
         instance < myLabels_->maxNumberOfParameterInstances();
         ++instance) {
      // check if this IOV/momentum range is to be ignored:
      if (std::find(ignoredInstances.begin(), ignoredInstances.end(), instance)
          != ignoredInstances.end()) {
        continue;
      }
      std::stringstream defaultsteerfileprefix;
      defaultsteerfileprefix << "autosteerFilePrefix_" << name << "_" << psetnr
                             << "_" << instance;

      const auto steerFilePrefix = pset.exists("steerFilePrefix") ?
        pset.getParameter<std::string> ("steerFilePrefix") + "_" + std::to_string(instance) :
        defaultsteerfileprefix.str();

      const auto levelsFilenames = this->makeLevelsFilenames(steerFilePrefixContainer,
                                                             alis,
                                                             steerFilePrefix);

      //Add the configuration data for this constraint to the container of config data
      ConstraintsConfigContainer_.emplace_back
    (GeometryConstraintConfigData(coefficients,
                      name,
                      levelsFilenames,
                      sysdeformation,
                      excluded_alis,
                      instance,
                      downToLowestLevel));
    }

  }
}

//_________________________________________________________________________
std::pair<align::GlobalPoint, align::GlobalPoint>
PedeSteererWeakModeConstraints::getDoubleSensorPosition(const Alignable *ali) const
{
  const auto aliPar =
    dynamic_cast<TwoBowedSurfacesAlignmentParameters*>(ali->alignmentParameters());
  if(aliPar) {
    const auto ySplit = aliPar->ySplit();
    const auto halfLength  = 0.5 * ali->surface().length();
    const auto yM1 = 0.5 * (ySplit - halfLength); // y_mean of surface 1
    const auto yM2 = yM1 + halfLength;            // y_mean of surface 2
    const auto pos_sensor0(ali->surface().toGlobal(align::LocalPoint(0.,yM1,0.)));
    const auto pos_sensor1(ali->surface().toGlobal(align::LocalPoint(0.,yM2,0.)));
    return std::make_pair(pos_sensor0, pos_sensor1);
  } else {
    throw cms::Exception("Alignment")
      << "[PedeSteererWeakModeConstraints::getDoubleSensorPosition]"
      << " Dynamic cast to double sensor parameters failed.";
    return std::make_pair( align::GlobalPoint(0.0, 0.0, 0.0), align::GlobalPoint(0.0, 0.0, 0.0));
  }
}

//_________________________________________________________________________
unsigned int
PedeSteererWeakModeConstraints::createAlignablesDataStructure()
{
  unsigned int nConstraints = 0;
  for(auto& iC: ConstraintsConfigContainer_) {
    //loop over all HLS for which the constraint is to be determined
    for(const auto& iHLS: iC.levelsFilenames_) {
      //determine next active sub-alignables for iHLS
      std::vector<Alignable*> aliDaughts;
      if (iC.downToLowestLevel_) {
    if (!iHLS.first->lastCompsWithParams(aliDaughts)) {
      edm::LogWarning("Alignment")
        << "@SUB=PedeSteererWeakModeConstraints::createAlignablesDataStructure"
        << "Some but not all component branches "
        << AlignableObjectId::idToString(iHLS.first->alignableObjectId())
        << " with params!";
    }
      } else {
    if (!iHLS.first->firstCompsWithParams(aliDaughts)) {
      edm::LogWarning("Alignment")
        << "@SUB=PedeSteererWeakModeConstraints::createAlignablesDataStructure"
        << "Some but not all daughters of "
        << AlignableObjectId::idToString(iHLS.first->alignableObjectId())
        << " with params!";
    }
      }
      ++nConstraints;

      std::list<Alignable*> usedinconstraint;
      for (const auto& iD: aliDaughts) {
        bool isNOTdead = true;
        for(const auto& iDeadmodules: deadmodules_) {
          if( (iD->alignableObjectId() == align::AlignableDetUnit
               || iD->alignableObjectId() == align::AlignableDet)
              && iD->geomDetId().rawId() == iDeadmodules) {
            isNOTdead = false;
            break;
          }
        }
        //check if the module is excluded
        for(const auto& iEx: iC.excludedAlignables_) {
          if(iD->id() == iEx->id() &&  iD->alignableObjectId() == iEx->alignableObjectId() ) {
            //if(iD->geomDetId().rawId() == (*iEx)->geomDetId().rawId()) {
            isNOTdead = false;
            break;
          }
        }
        const bool issubcomponent = this->checkMother(iD,iHLS.first);
        if(issubcomponent) {
          if(isNOTdead) {
            usedinconstraint.push_back(iD);
          }
        } else {
          //sanity check
          throw cms::Exception("Alignment")
            << "[PedeSteererWeakModeConstraints::createAlignablesDataStructure]"
            << " Sanity check failed. Alignable defined as active sub-component, "
            << " but in fact its not a daugther of " << AlignableObjectId::idToString(iHLS.first->alignableObjectId());
        }
      }

      if( usedinconstraint.size() > 0){
        iC.HLSsubdets_.push_back(std::make_pair(iHLS.first, usedinconstraint));
      } else {
        edm::LogInfo("Alignment") << "@SUB=PedeSteererWeakModeConstraints"
                                  << "No sub-components for "
                                  << AlignableObjectId::idToString(iHLS.first->alignableObjectId())
                                  << " at (" << iHLS.first->globalPosition().x()
                  << ","<< iHLS.first->globalPosition().y()
                  << "," << iHLS.first->globalPosition().z()
                  << ") selected. Skip constraint";
      }
      if(aliDaughts.size() == 0) {
        edm::LogWarning("Alignment") << "@SUB=PedeSteererWeakModeConstraints::createAlignablesDataStructure"
                                     << "No active sub-alignables found for "
                                     << AlignableObjectId::idToString(iHLS.first->alignableObjectId())
                                     << " at (" << iHLS.first->globalPosition().x() << ","<< iHLS.first->globalPosition().y() << "," << iHLS.first->globalPosition().z() << ").";
      }


    }
  }
  return nConstraints;
}

//_________________________________________________________________________
double
PedeSteererWeakModeConstraints::getX(const int sysdeformation,
                     const align::GlobalPoint &pos,
                     const double phase) const
{
  double x = 0.0;

  const double r = TMath::Sqrt(pos.x() * pos.x() + pos.y() * pos.y());

  switch(sysdeformation) {
  case SystematicDeformations::kTwist:
  case SystematicDeformations::kZexpansion:
    x = pos.z();
    break;
  case SystematicDeformations::kSagitta:
  case SystematicDeformations::kRadial:
  case SystematicDeformations::kTelescope:
  case SystematicDeformations::kLayerRotation:
    x = r;
    break;
  case SystematicDeformations::kBowing:
    x = pos.z() * pos.z(); //TMath::Abs(pos.z());
    break;
  case SystematicDeformations::kElliptical:
    x = r * TMath::Cos(2.0 * pos.phi() + phase);
    break;
  case SystematicDeformations::kSkew:
    x = TMath::Cos(pos.phi() + phase);
    break;
  };

  return x;
}

//_________________________________________________________________________
double
PedeSteererWeakModeConstraints::getCoefficient(const int sysdeformation,
                           const align::GlobalPoint &pos,
                           const GlobalPoint gUDirection,
                           const GlobalPoint gVDirection,
                           const GlobalPoint gWDirection,
                           const int iParameter, const double &x0,
                           const std::vector<double> &constraintparameters) const
{


  if(iParameter < 0 || iParameter > 2) {
    throw cms::Exception("Alignment")
      << "[PedeSteererWeakModeConstraints::getCoefficient]" << " iParameter has to be in the range [0,2] but"
      << " it is equal to " << iParameter << ".";
  }


  //global vectors pointing in u,v,w direction
  const std::vector<double> vec_u = {pos.x() - gUDirection.x(), pos.y() - gUDirection.y(), pos.z() - gUDirection.z()};
  const std::vector<double> vec_v = {pos.x() - gVDirection.x(), pos.y() - gVDirection.y(), pos.z() - gVDirection.z()};
  const std::vector<double> vec_w = {pos.x() - gWDirection.x(), pos.y() - gWDirection.y(), pos.z() - gWDirection.z()};

  //FIXME: how to make inner vectors const?
  const std::vector<std::vector<double> > global_vecs = {vec_u, vec_v, vec_w};

  const double n = TMath::Sqrt( global_vecs.at(iParameter).at(0) * global_vecs.at(iParameter).at(0)
                                + global_vecs.at(iParameter).at(1) * global_vecs.at(iParameter).at(1)
                                + global_vecs.at(iParameter).at(2) * global_vecs.at(iParameter).at(2) );
  const double r = TMath::Sqrt( pos.x() * pos.x() + pos.y() * pos.y() );

  const double phase = this->getPhase(constraintparameters);
  //const double radial_direction[3] = {TMath::Sin(phase), TMath::Cos(phase), 0.0};
  const std::vector<double> radial_direction = {TMath::Sin(phase), TMath::Cos(phase), 0.0};
  //is equal to unity by construction ...
  const double norm_radial_direction =  TMath::Sqrt(radial_direction.at(0) * radial_direction.at(0)
                                                    + radial_direction.at(1) * radial_direction.at(1)
                                                    + radial_direction.at(2) * radial_direction.at(2));

  //const double phi_direction[3] = { -1.0 * pos.y(), pos.x(), 0.0};
  const std::vector<double> phi_direction = {-1.0 * pos.y(), pos.x(), 0.0};
  const double norm_phi_direction = TMath::Sqrt(phi_direction.at(0) * phi_direction.at(0)
                                                + phi_direction.at(1) * phi_direction.at(1)
                                                + phi_direction.at(2) * phi_direction.at(2));

  //const double z_direction[3] = {0.0, 0.0, 1.0};
  static const std::vector<double> z_direction = {0.0, 0.0, 1.0};
  const double norm_z_direction = TMath::Sqrt(z_direction.at(0)*z_direction.at(0)
                                              + z_direction.at(1)*z_direction.at(1)
                                              + z_direction.at(2)*z_direction.at(2));

  //unit vector pointing from the origin to the module position in the transverse plane
  const std::vector<double> rDirection = {pos.x(), pos.y(), 0.0};
  const double norm_rDirection =  TMath::Sqrt(rDirection.at(0) * rDirection.at(0)
                                              + rDirection.at(1) * rDirection.at(1)
                                              + rDirection.at(2) * rDirection.at(2));

  double coeff = 0.0;
  double dot_product = 0.0;
  double normalisation_factor = 1.0;

  //see https://indico.cern.ch/getFile.py/access?contribId=15&sessionId=1&resId=0&materialId=slides&confId=127126
  switch(sysdeformation) {
  case SystematicDeformations::kTwist:
  case SystematicDeformations::kLayerRotation:
    dot_product = phi_direction.at(0) * global_vecs.at(iParameter).at(0)
      + phi_direction.at(1) * global_vecs.at(iParameter).at(1)
      + phi_direction.at(2) * global_vecs.at(iParameter).at(2);
    normalisation_factor = r * n * norm_phi_direction;
    break;
  case SystematicDeformations::kZexpansion :
  case SystematicDeformations::kTelescope:
  case SystematicDeformations::kSkew:
    dot_product = global_vecs.at(iParameter).at(0) * z_direction.at(0)
      + global_vecs.at(iParameter).at(1) * z_direction.at(1)
      + global_vecs.at(iParameter).at(2) * z_direction.at(2);
    normalisation_factor = ( n * norm_z_direction);
    break;
  case SystematicDeformations::kRadial:
  case SystematicDeformations::kBowing:
  case SystematicDeformations::kElliptical:
    dot_product = global_vecs.at(iParameter).at(0) * rDirection.at(0)
      + global_vecs.at(iParameter).at(1) * rDirection.at(1)
      + global_vecs.at(iParameter).at(2) * rDirection.at(2);
    normalisation_factor = ( n * norm_rDirection);
    break;
  case SystematicDeformations::kSagitta:
    dot_product = global_vecs.at(iParameter).at(0) * radial_direction.at(0)
      + global_vecs.at(iParameter).at(1) * radial_direction.at(1)
      + global_vecs.at(iParameter).at(2) * radial_direction.at(2);
    normalisation_factor = ( n * norm_radial_direction);
    break;
  default: break;
  }

  if(TMath::Abs(normalisation_factor) > 0.0) {
    coeff = dot_product * ( this->getX(sysdeformation,pos,phase) - x0 ) / normalisation_factor;
  } else {
    throw cms::Exception("Alignment")
      << "[PedeSteererWeakModeConstraints::getCoefficient]" << " Normalisation factor"
      << "for coefficient calculation equal to zero! Misconfiguration?";
  }
  return coeff;
}

//_________________________________________________________________________
bool
PedeSteererWeakModeConstraints::checkSelectionShiftParameter(const Alignable *ali,
                                 unsigned int iParameter) const
{
  bool isselected = false;
  const std::vector<bool> &aliSel= ali->alignmentParameters()->selector();
  //exclude non-shift parameters
  if((iParameter <= 2)
     || (iParameter >= 9 && iParameter <=11)) {
    if(!aliSel.at(iParameter)) {
      isselected = false;
    } else {
      auto params = ali->alignmentParameters();
      auto selVar = dynamic_cast<SelectionUserVariables*>(params->userVariables());
      if (selVar) {
        if(selVar->fullSelection().size() <= (iParameter+1)) {
          throw cms::Exception("Alignment")
            << "[PedeSteererWeakModeConstraints::checkSelectionShiftParameter]"
            << " Can not access selected alignment variables of alignable "
            <<  AlignableObjectId::idToString(ali->alignableObjectId())
            << "at (" << ali->globalPosition().x() << ","<< ali->globalPosition().y() << "," << ali->globalPosition().z()<< ") "
            << "for parameter number " << (iParameter+1) << ".";
        }
      }
      const char selChar = (selVar ? selVar->fullSelection().at(iParameter) : '1');
      // if(selChar == '1') { //FIXME??? what about 'r'?
      if(selChar == '1' || selChar == 'r') {
        isselected = true;
      } else {
        isselected = false;
      }
    }
  }
  return isselected;
}

//_________________________________________________________________________
void
PedeSteererWeakModeConstraints::closeOutputfiles()
{
   //'delete' output files which means: close them
  for(auto& it: ConstraintsConfigContainer_) {
    for(auto& iFile:it.mapFileName_) {
      if(iFile.second) {
        delete iFile.second;
        iFile.second = nullptr;
      } else {
        throw cms::Exception("FileCloseProblem")
          << "[PedeSteererWeakModeConstraints]" << " can not close file " << iFile.first << ".";
      }
    }
  }
}

//_________________________________________________________________________
void
PedeSteererWeakModeConstraints::writeOutput(const std::list<std::pair<unsigned int,double> > &output,
                        const GeometryConstraintConfigData &it,
                        const Alignable* iHLS,
                        double sum_xi_x0)
{
  std::ofstream* ofile = getFile(it, iHLS);

  if(ofile == nullptr) {
    throw cms::Exception("FileFindError")
      << "[PedeSteererWeakModeConstraints] Cannot find output file.";
  } else {
    if(output.size() > 0) {
      const double constr = sum_xi_x0 * it.coefficients_.front();
      (*ofile) << "Constraint " << std::scientific << constr << std::endl;
      for(const auto& ioutput: output) {
        (*ofile) << std::fixed << ioutput.first << " " << std::scientific << ioutput.second << std::endl;
      }
    }
  }
}

//_________________________________________________________________________
std::ofstream*
PedeSteererWeakModeConstraints::getFile(const GeometryConstraintConfigData &it,
                    const Alignable* iHLS) const
{
  std::ofstream* file = nullptr;

  for(const auto& ilevelsFilename: it.levelsFilenames_) {
    if(ilevelsFilename.first->id() == iHLS->id() &&
       ilevelsFilename.first->alignableObjectId() == iHLS->alignableObjectId()) {

      const auto iFile = it.mapFileName_.find(ilevelsFilename.second);
      if(iFile != it.mapFileName_.end()) {
        file = iFile->second;
      }
    }
  }

  return file;
}

//_________________________________________________________________________
double
PedeSteererWeakModeConstraints::getX0(const std::pair<Alignable*, std::list<Alignable*> > &iHLS,
                      const GeometryConstraintConfigData &it) const
{
  double nmodules = 0.0;
  double x0 = 0.0;

  for(const auto& ali: iHLS.second) {
    align::PositionType pos = ali->globalPosition();
    bool alignableIsFloating = false; //means: true=alignable is able to move in at least one direction

    //test whether at least one variable has been selected in the configuration
    for(unsigned int iParameter = 0;
        static_cast<int>(iParameter) < ali->alignmentParameters()->size(); iParameter++) {
      if(this->checkSelectionShiftParameter(ali,iParameter) ) {
        alignableIsFloating = true;
        // //verify that alignable has just one label -- meaning no IOV-dependence etc
        // const unsigned int nInstances = myLabels_->numberOfParameterInstances(ali, iParameter);
        // if(nInstances > 1) {
        //   throw cms::Exception("PedeSteererWeakModeConstraints")
        //     << "@SUB=PedeSteererWeakModeConstraints::ConstructConstraints"
        //     << " Weak mode constraints are only supported for alignables which have"
        //     << " just one label. However, e.g. alignable"
        //     << " " << AlignableObjectId::idToString(ali->alignableObjectId())
        //     << "at (" << ali->globalPosition().x() << ","<< ali->globalPosition().y() << "," << ali->globalPosition().z()<< "), "
        //     << " was configured to have >1 label. Remove e.g. IOV-dependence for this (and other) alignables which are used in the constraint.";
        // }
        break;
      }
    }
    //at least one parameter of the alignable can be changed in the alignment
    if(alignableIsFloating) {
      const auto phase = this->getPhase(it.coefficients_);
      if(ali->alignmentParameters()->type() != AlignmentParametersFactory::kTwoBowedSurfaces ) {
        x0 += this->getX(it.sysdeformation_,pos,phase);
        nmodules++;
      } else {
        std::pair<align::GlobalPoint, align::GlobalPoint> sensorpositions = this->getDoubleSensorPosition(ali);
        x0 += this->getX(it.sysdeformation_,sensorpositions.first,phase) + this->getX(it.sysdeformation_,sensorpositions.second,phase);
        nmodules++;
        nmodules++;
      }
    }
  }
  if(nmodules>0) {
    x0 = x0 / nmodules;
  } else {
    throw cms::Exception("Alignment")
      << "@SUB=PedeSteererWeakModeConstraints::ConstructConstraints"
      << " Number of selected modules equal to zero. Check configuration!";
    x0 = 1.0;
  }
  return x0;
}

//_________________________________________________________________________
unsigned int
PedeSteererWeakModeConstraints::constructConstraints(const std::vector<Alignable*> &alis)
{
  //FIXME: split the code of the method into smaller pieces/submethods

  //create the data structures that store the alignables
  //for which the constraints need to be calculated and
  //their association to high-level structures
  const auto nConstraints = this->createAlignablesDataStructure();

  std::vector<std::list<std::pair<unsigned int,double> > > createdConstraints;

  //calculate constraints
  //loop over all constraints
  for(const auto& it: ConstraintsConfigContainer_) {

    //loop over all subdets for which constraints are determined
    for(const auto& iHLS : it.HLSsubdets_) {
      double sum_xi_x0 = 0.0;
      std::list<std::pair<unsigned int,double> > output;

      const double x0 = this->getX0(iHLS, it);

      for(std::list<Alignable*>::const_iterator iAlignables = iHLS.second.begin();
          iAlignables != iHLS.second.end(); iAlignables++) {
        const Alignable *ali = (*iAlignables);
        const auto aliLabel = myLabels_->alignableLabelFromParamAndInstance(const_cast<Alignable*>(ali), 0, it.instance_);
        const AlignableSurface &surface = ali->surface();

        const LocalPoint  lUDirection(1.,0.,0.),
          lVDirection(0.,1.,0.),
          lWDirection(0.,0.,1.);

        GlobalPoint gUDirection = surface.toGlobal(lUDirection),
          gVDirection = surface.toGlobal(lVDirection),
          gWDirection = surface.toGlobal(lWDirection);

        const bool isDoubleSensor = ali->alignmentParameters()->type() == AlignmentParametersFactory::kTwoBowedSurfaces;

        const auto sensorpositions =
          isDoubleSensor ? this->getDoubleSensorPosition(ali) : std::make_pair(ali->globalPosition(), align::PositionType ());

        const auto& pos_sensor0 = sensorpositions.first;
        const auto& pos_sensor1 = sensorpositions.second;
        const auto phase = this->getPhase(it.coefficients_);
        const auto x_sensor0 = this->getX(it.sysdeformation_,pos_sensor0,phase);
        const auto x_sensor1 = isDoubleSensor ? this->getX(it.sysdeformation_,pos_sensor1,phase) : 0.0;

        sum_xi_x0 += ( x_sensor0 - x0 ) * ( x_sensor0 - x0 );
        if(isDoubleSensor) {
          sum_xi_x0 += ( x_sensor1 - x0 ) * ( x_sensor1 - x0 );
        }
        const int numparameterlimit = ali->alignmentParameters()->size(); //isDoubleSensor ? 18 : 3;

        for(int iParameter = 0; iParameter < numparameterlimit; iParameter++) {
          int localindex = 0;
          if(iParameter == 0 || iParameter == 9)
            localindex = 0;
          if(iParameter == 1 || iParameter == 10)
            localindex = 1;
          if(iParameter == 2 || iParameter == 11)
            localindex = 2;

          if((iParameter >= 0 && iParameter <= 2)
             || (iParameter >= 9 && iParameter <=11)) {
          } else {
            continue;
          }
          if(! this->checkSelectionShiftParameter(ali,iParameter) ) {
            continue;
          }
          //do it for each 'instance' separately? -> IOV-dependence, no
          const auto paramLabel = myLabels_->parameterLabel(aliLabel, iParameter);

          const auto& pos = (iParameter <= 2) ? pos_sensor0 : pos_sensor1;
          //select only u,v,w
          if(iParameter == 0 || iParameter == 1 || iParameter == 2
             || iParameter == 9 || iParameter == 10 || iParameter == 11) {
            const double coeff = this->getCoefficient(it.sysdeformation_,
                                                      pos,
                                                      gUDirection,
                                                      gVDirection,
                                                      gWDirection,
                                                      localindex,
                                                      x0,
                                                      it.coefficients_);
            if(TMath::Abs(coeff) > 0.0) {
              //nothing
            } else {
              edm::LogWarning("PedeSteererWeakModeConstraints")
                << "@SUB=PedeSteererWeakModeConstraints::getCoefficient"
                << "Coefficient of alignable "
                <<  AlignableObjectId::idToString(ali->alignableObjectId())
                << " at (" << ali->globalPosition().x() << ","<< ali->globalPosition().y() << "," << ali->globalPosition().z()<< ") "
                << " in subdet " << AlignableObjectId::idToString(iHLS.first->alignableObjectId())
                << " for parameter " << localindex << " equal to zero. This alignable is used in the constraint"
                << " '" << it.constraintName_ << "'. The id is: alignable->geomDetId().rawId() = "
                << ali->geomDetId().rawId() << ".";
            }
            output.push_back(std::make_pair (paramLabel, coeff));
          }
        }


      }

      if (std::find(createdConstraints.begin(), createdConstraints.end(), output)
      != createdConstraints.end()) {
    // check if linearly dependent constraint exists already:
    auto outFile = getFile(it, iHLS.first);
    if (outFile == nullptr) {
      throw cms::Exception("FileFindError")
        << "[PedeSteererWeakModeConstraints] Cannot find output file.";
    } else {
      *outFile << "! The constraint for this IOV/momentum range" << std::endl
           << "! has been removed because the used parameters" << std::endl
           << "! are not IOV or momentum-range dependent." << std::endl;
    }
    continue;
      }
      this->writeOutput(output, it, iHLS.first, sum_xi_x0);
      createdConstraints.push_back(output);
    }
  }
  this->closeOutputfiles();

  return nConstraints;
}

//_________________________________________________________________________
bool
PedeSteererWeakModeConstraints::checkMother(const Alignable * const lowleveldet,
                        const Alignable * const HLS) const
{
  if(lowleveldet->id() == HLS->id() && lowleveldet->alignableObjectId() == HLS->alignableObjectId()) {
    return true;
  } else {
    if(lowleveldet->mother() == nullptr)
      return false;
    else
      return this->checkMother(lowleveldet->mother(),HLS);
  }
}

//_________________________________________________________________________
void
PedeSteererWeakModeConstraints::verifyParameterNames(const edm::ParameterSet &pset,
                             unsigned int psetnr) const
{
  const auto parameterNames = pset.getParameterNames();
  for (const auto& name: parameterNames) {
    if(name != "coefficients"
       && name != "deadmodules" && name != "constraint"
       && name != "steerFilePrefix" && name != "levels"
       && name != "excludedAlignables" && name != "ignoredInstances"
       && name != "downToLowestLevel"
       ) {
      throw cms::Exception("BadConfig")
        << "@SUB=PedeSteererWeakModeConstraints::verifyParameterNames:"
        << " Unknown parameter name '" << name << "' in PSet number " << psetnr << ". Maybe a typo?";
    }
  }
}

//_________________________________________________________________________
const std::vector<std::pair<Alignable*, std::string> >
PedeSteererWeakModeConstraints::makeLevelsFilenames(std::set<std::string> &steerFilePrefixContainer,
                                                    const std::vector<Alignable*> &alis,
                                                    const std::string &steerFilePrefix) const
{
  //check whether the prefix is unique
  if(steerFilePrefixContainer.find(steerFilePrefix) != steerFilePrefixContainer.end()) {
    throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints] Steering file"
                                      << " prefix '" << steerFilePrefix << "' already exists. Specify unique names!";
  } else {
    steerFilePrefixContainer.insert(steerFilePrefix);
  }

  std::vector<std::pair<Alignable*, std::string> > levelsFilenames;
  for(const auto& ali: alis) {
    std::stringstream n;
    n << steerFile_ << "_" << steerFilePrefix //<< "_" << name
      << "_" << AlignableObjectId::idToString(ali->alignableObjectId())
      << "_" << ali->id() << "_" << ali->alignableObjectId() << ".txt";

    levelsFilenames.push_back(std::make_pair(ali, n.str()));
  }
  return levelsFilenames;
}

//_________________________________________________________________________
int
PedeSteererWeakModeConstraints::verifyDeformationName(const std::string &name,
                              const std::vector<double> &coefficients) const
{
  int sysdeformation = SystematicDeformations::kUnknown;

  if(name == "twist") {
    sysdeformation = SystematicDeformations::kTwist;
  }  else if(name == "zexpansion") {
    sysdeformation = SystematicDeformations::kZexpansion;
  } else if(name == "sagitta") {
    sysdeformation = SystematicDeformations::kSagitta;
  } else if(name == "radial") {
    sysdeformation = SystematicDeformations::kRadial;
  } else if(name == "telescope") {
    sysdeformation = SystematicDeformations::kTelescope;
  } else if(name == "layerrotation") {
    sysdeformation = SystematicDeformations::kLayerRotation;
  } else if(name == "bowing") {
    sysdeformation = SystematicDeformations::kBowing;
  } else if(name == "skew") {
    sysdeformation = SystematicDeformations::kSkew;
  } else if(name == "elliptical") {
    sysdeformation = SystematicDeformations::kElliptical;
  }

  if(sysdeformation == SystematicDeformations::kUnknown) {
    throw cms::Exception("BadConfig")
      << "[PedeSteererWeakModeConstraints]" << " specified configuration option '"
      << name << "' not known.";
  }
  if((sysdeformation == SystematicDeformations::kSagitta
      || sysdeformation == SystematicDeformations::kElliptical
      || sysdeformation == SystematicDeformations::kSkew) && coefficients.size() != 2) {
    throw cms::Exception("BadConfig")
      << "[PedeSteererWeakModeConstraints]" << " Excactly two parameters using the coefficient"
      << " variable have to be provided for the " << name << " constraint.";
  }
  if((sysdeformation == SystematicDeformations::kTwist
      || sysdeformation == SystematicDeformations::kZexpansion
      || sysdeformation == SystematicDeformations::kTelescope
      || sysdeformation == SystematicDeformations::kLayerRotation
      || sysdeformation == SystematicDeformations::kRadial
      || sysdeformation == SystematicDeformations::kBowing) && coefficients.size() != 1) {
    throw cms::Exception("BadConfig")
      << "[PedeSteererWeakModeConstraints]" << " Excactly ONE parameter using the coefficient"
      << " variable have to be provided for the " << name << " constraint.";
  }

  if(coefficients.size() == 0) {
    throw cms::Exception("BadConfig")
      << "[PedeSteererWeakModeConstraints]" << " At least one coefficient has to be specified.";
  }
  return sysdeformation;
}


//_________________________________________________________________________
double
PedeSteererWeakModeConstraints::getPhase(const std::vector<double> &coefficients) const
{
  return coefficients.size() == 2 ? coefficients.at(1) : 0.0; //treat second parameter as phase otherwise return 0
}

//_________________________________________________________________________
PedeSteererWeakModeConstraints::~PedeSteererWeakModeConstraints() = default;