Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:18

0001 /**
0002  * \file PedeSteererWeakModeConstraints.cc
0003  *
0004  *  \author    : Joerg Behr
0005  *  date       : February 2013
0006  */
0007 
0008 #include "Alignment/MillePedeAlignmentAlgorithm/src/PedeSteererWeakModeConstraints.h"
0009 #include "Alignment/MillePedeAlignmentAlgorithm/src/PedeSteerer.h"
0010 
0011 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
0012 
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 
0015 #include "Alignment/CommonAlignment/interface/Alignable.h"
0016 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
0017 #include "Alignment/CommonAlignment/interface/Utilities.h"
0018 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0019 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterSelector.h"
0020 #include "Alignment/CommonAlignmentAlgorithm/interface/SelectionUserVariables.h"
0021 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
0022 #include "Alignment/CommonAlignmentParametrization/interface/RigidBodyAlignmentParameters.h"
0023 #include "Alignment/CommonAlignmentParametrization/interface/BowedSurfaceAlignmentDerivatives.h"
0024 #include "Alignment/CommonAlignmentParametrization/interface/TwoBowedSurfacesAlignmentParameters.h"
0025 
0026 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
0027 // for 'type identification' as Alignable
0028 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0029 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0030 #include "Alignment/CommonAlignment/interface/AlignableExtras.h"
0031 // GF doubts the need of these includes from include checker campaign:
0032 #include <FWCore/Framework/interface/EventSetup.h>
0033 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0034 #include <Geometry/CommonDetUnit/interface/GeomDetType.h>
0035 #include <DataFormats/GeometrySurface/interface/LocalError.h>
0036 #include <Geometry/DTGeometry/interface/DTLayer.h>
0037 // end of doubt
0038 
0039 #include <DataFormats/GeometryVector/interface/GlobalPoint.h>
0040 
0041 #include <fstream>
0042 #include <sstream>
0043 #include <algorithm>
0044 
0045 // from ROOT
0046 #include <TSystem.h>
0047 #include <TMath.h>
0048 
0049 #include <iostream>
0050 
0051 GeometryConstraintConfigData::GeometryConstraintConfigData(
0052     const std::vector<double>& co,
0053     const std::string& c,
0054     const std::vector<std::pair<Alignable*, std::string> >& alisFile,
0055     const int sd,
0056     const align::Alignables& ex,
0057     const int instance,
0058     const bool downToLowestLevel)
0059     : coefficients_(co),
0060       constraintName_(c),
0061       levelsFilenames_(alisFile),
0062       excludedAlignables_(ex),
0063       sysdeformation_(sd),
0064       instance_(instance),
0065       downToLowestLevel_(downToLowestLevel) {}
0066 
0067 //_________________________________________________________________________
0068 PedeSteererWeakModeConstraints::PedeSteererWeakModeConstraints(AlignableTracker* aliTracker,
0069                                                                const PedeLabelerBase* labels,
0070                                                                const std::vector<edm::ParameterSet>& config,
0071                                                                std::string sf)
0072     : myLabels_(labels),
0073       myConfig_(config),
0074       steerFile_(sf),
0075       alignableObjectId_{AlignableObjectId::commonObjectIdProvider(aliTracker, nullptr)} {
0076   unsigned int psetnr = 0;
0077   std::set<std::string> steerFilePrefixContainer;
0078   for (const auto& pset : myConfig_) {
0079     this->verifyParameterNames(pset, psetnr);
0080     psetnr++;
0081 
0082     const auto coefficients = pset.getParameter<std::vector<double> >("coefficients");
0083     const auto dm = pset.exists("deadmodules") ? pset.getParameter<std::vector<unsigned int> >("deadmodules")
0084                                                : std::vector<unsigned int>();
0085     std::string name = pset.getParameter<std::string>("constraint");
0086     std::transform(name.begin(), name.end(), name.begin(), ::tolower);
0087     const auto ignoredInstances = pset.exists("ignoredInstances")
0088                                       ? pset.getUntrackedParameter<std::vector<unsigned int> >("ignoredInstances")
0089                                       : std::vector<unsigned int>();
0090 
0091     const auto downToLowestLevel =
0092         pset.exists("downToLowestLevel") ? pset.getUntrackedParameter<bool>("downToLowestLevel") : false;
0093 
0094     AlignmentParameterSelector selector(aliTracker, nullptr, nullptr);
0095     selector.clear();
0096     selector.addSelections(pset.getParameter<edm::ParameterSet>("levels"));
0097 
0098     const auto& alis = selector.selectedAlignables();
0099 
0100     AlignmentParameterSelector selector_excludedalignables(aliTracker, nullptr, nullptr);
0101     selector_excludedalignables.clear();
0102     if (pset.exists("excludedAlignables")) {
0103       selector_excludedalignables.addSelections(pset.getParameter<edm::ParameterSet>("excludedAlignables"));
0104     }
0105     const auto& excluded_alis = selector_excludedalignables.selectedAlignables();
0106 
0107     //check that the name of the deformation is known and that the number
0108     //of provided parameter is right.
0109     auto sysdeformation = this->verifyDeformationName(name, coefficients);
0110 
0111     if (deadmodules_.empty()) {  //fill the list of dead modules only once
0112       edm::LogInfo("Alignment") << "@SUB=PedeSteererWeakModeConstraints"
0113                                 << "Load list of dead modules (size = " << dm.size() << ").";
0114       for (const auto& it : dm)
0115         deadmodules_.push_back(it);
0116     }
0117 
0118     // loop over all IOVs/momentum ranges
0119     for (unsigned int instance = 0; instance < myLabels_->maxNumberOfParameterInstances(); ++instance) {
0120       // check if this IOV/momentum range is to be ignored:
0121       if (std::find(ignoredInstances.begin(), ignoredInstances.end(), instance) != ignoredInstances.end()) {
0122         continue;
0123       }
0124       std::stringstream defaultsteerfileprefix;
0125       defaultsteerfileprefix << "autosteerFilePrefix_" << name << "_" << psetnr << "_" << instance;
0126 
0127       const auto steerFilePrefix = pset.exists("steerFilePrefix") ? pset.getParameter<std::string>("steerFilePrefix") +
0128                                                                         "_" + std::to_string(instance)
0129                                                                   : defaultsteerfileprefix.str();
0130 
0131       const auto levelsFilenames = this->makeLevelsFilenames(steerFilePrefixContainer, alis, steerFilePrefix);
0132 
0133       //Add the configuration data for this constraint to the container of config data
0134       ConstraintsConfigContainer_.emplace_back(GeometryConstraintConfigData(
0135           coefficients, name, levelsFilenames, sysdeformation, excluded_alis, instance, downToLowestLevel));
0136     }
0137   }
0138 }
0139 
0140 //_________________________________________________________________________
0141 std::pair<align::GlobalPoint, align::GlobalPoint> PedeSteererWeakModeConstraints::getDoubleSensorPosition(
0142     const Alignable* ali) const {
0143   const auto aliPar = dynamic_cast<TwoBowedSurfacesAlignmentParameters*>(ali->alignmentParameters());
0144   if (aliPar) {
0145     const auto ySplit = aliPar->ySplit();
0146     const auto halfLength = 0.5 * ali->surface().length();
0147     const auto yM1 = 0.5 * (ySplit - halfLength);  // y_mean of surface 1
0148     const auto yM2 = yM1 + halfLength;             // y_mean of surface 2
0149     const auto pos_sensor0(ali->surface().toGlobal(align::LocalPoint(0., yM1, 0.)));
0150     const auto pos_sensor1(ali->surface().toGlobal(align::LocalPoint(0., yM2, 0.)));
0151     return std::make_pair(pos_sensor0, pos_sensor1);
0152   } else {
0153     throw cms::Exception("Alignment") << "[PedeSteererWeakModeConstraints::getDoubleSensorPosition]"
0154                                       << " Dynamic cast to double sensor parameters failed.";
0155     return std::make_pair(align::GlobalPoint(0.0, 0.0, 0.0), align::GlobalPoint(0.0, 0.0, 0.0));
0156   }
0157 }
0158 
0159 //_________________________________________________________________________
0160 unsigned int PedeSteererWeakModeConstraints::createAlignablesDataStructure() {
0161   unsigned int nConstraints = 0;
0162   for (auto& iC : ConstraintsConfigContainer_) {
0163     //loop over all HLS for which the constraint is to be determined
0164     for (const auto& iHLS : iC.levelsFilenames_) {
0165       //determine next active sub-alignables for iHLS
0166       align::Alignables aliDaughts;
0167       if (iC.downToLowestLevel_) {
0168         if (!iHLS.first->lastCompsWithParams(aliDaughts)) {
0169           edm::LogWarning("Alignment") << "@SUB=PedeSteererWeakModeConstraints::createAlignablesDataStructure"
0170                                        << "Some but not all component branches "
0171                                        << alignableObjectId_.idToString(iHLS.first->alignableObjectId())
0172                                        << " with params!";
0173         }
0174       } else {
0175         if (!iHLS.first->firstCompsWithParams(aliDaughts)) {
0176           edm::LogWarning("Alignment") << "@SUB=PedeSteererWeakModeConstraints::createAlignablesDataStructure"
0177                                        << "Some but not all daughters of "
0178                                        << alignableObjectId_.idToString(iHLS.first->alignableObjectId())
0179                                        << " with params!";
0180         }
0181       }
0182       ++nConstraints;
0183 
0184       std::list<Alignable*> usedinconstraint;
0185       for (const auto& iD : aliDaughts) {
0186         bool isNOTdead = true;
0187         for (const auto& iDeadmodules : deadmodules_) {
0188           if ((iD->alignableObjectId() == align::AlignableDetUnit || iD->alignableObjectId() == align::AlignableDet) &&
0189               iD->geomDetId().rawId() == iDeadmodules) {
0190             isNOTdead = false;
0191             break;
0192           }
0193         }
0194         //check if the module is excluded
0195         for (const auto& iEx : iC.excludedAlignables_) {
0196           if (iD->id() == iEx->id() && iD->alignableObjectId() == iEx->alignableObjectId()) {
0197             //if(iD->geomDetId().rawId() == (*iEx)->geomDetId().rawId()) {
0198             isNOTdead = false;
0199             break;
0200           }
0201         }
0202         const bool issubcomponent = this->checkMother(iD, iHLS.first);
0203         if (issubcomponent) {
0204           if (isNOTdead) {
0205             usedinconstraint.push_back(iD);
0206           }
0207         } else {
0208           //sanity check
0209           throw cms::Exception("Alignment") << "[PedeSteererWeakModeConstraints::createAlignablesDataStructure]"
0210                                             << " Sanity check failed. Alignable defined as active sub-component, "
0211                                             << " but in fact its not a daugther of "
0212                                             << alignableObjectId_.idToString(iHLS.first->alignableObjectId());
0213         }
0214       }
0215 
0216       if (!usedinconstraint.empty()) {
0217         iC.HLSsubdets_.push_back(std::make_pair(iHLS.first, usedinconstraint));
0218       } else {
0219         edm::LogInfo("Alignment") << "@SUB=PedeSteererWeakModeConstraints"
0220                                   << "No sub-components for "
0221                                   << alignableObjectId_.idToString(iHLS.first->alignableObjectId()) << " at ("
0222                                   << iHLS.first->globalPosition().x() << "," << iHLS.first->globalPosition().y() << ","
0223                                   << iHLS.first->globalPosition().z() << ") selected. Skip constraint";
0224       }
0225       if (aliDaughts.empty()) {
0226         edm::LogWarning("Alignment") << "@SUB=PedeSteererWeakModeConstraints::createAlignablesDataStructure"
0227                                      << "No active sub-alignables found for "
0228                                      << alignableObjectId_.idToString(iHLS.first->alignableObjectId()) << " at ("
0229                                      << iHLS.first->globalPosition().x() << "," << iHLS.first->globalPosition().y()
0230                                      << "," << iHLS.first->globalPosition().z() << ").";
0231       }
0232     }
0233   }
0234   return nConstraints;
0235 }
0236 
0237 //_________________________________________________________________________
0238 double PedeSteererWeakModeConstraints::getX(const int sysdeformation,
0239                                             const align::GlobalPoint& pos,
0240                                             const double phase) const {
0241   double x = 0.0;
0242 
0243   const double r = TMath::Sqrt(pos.x() * pos.x() + pos.y() * pos.y());
0244 
0245   switch (sysdeformation) {
0246     case SystematicDeformations::kTwist:
0247     case SystematicDeformations::kZexpansion:
0248       x = pos.z();
0249       break;
0250     case SystematicDeformations::kSagitta:
0251     case SystematicDeformations::kRadial:
0252     case SystematicDeformations::kTelescope:
0253     case SystematicDeformations::kLayerRotation:
0254       x = r;
0255       break;
0256     case SystematicDeformations::kBowing:
0257       x = pos.z() * pos.z();  //TMath::Abs(pos.z());
0258       break;
0259     case SystematicDeformations::kElliptical:
0260       x = r * TMath::Cos(2.0 * pos.phi() + phase);
0261       break;
0262     case SystematicDeformations::kSkew:
0263       x = TMath::Cos(pos.phi() + phase);
0264       break;
0265   };
0266 
0267   return x;
0268 }
0269 
0270 //_________________________________________________________________________
0271 double PedeSteererWeakModeConstraints::getCoefficient(const int sysdeformation,
0272                                                       const align::GlobalPoint& pos,
0273                                                       const GlobalPoint gUDirection,
0274                                                       const GlobalPoint gVDirection,
0275                                                       const GlobalPoint gWDirection,
0276                                                       const int iParameter,
0277                                                       const double& x0,
0278                                                       const std::vector<double>& constraintparameters) const {
0279   if (iParameter < 0 || iParameter > 2) {
0280     throw cms::Exception("Alignment") << "[PedeSteererWeakModeConstraints::getCoefficient]"
0281                                       << " iParameter has to be in the range [0,2] but"
0282                                       << " it is equal to " << iParameter << ".";
0283   }
0284 
0285   //global vectors pointing in u,v,w direction
0286   const std::vector<double> vec_u = {pos.x() - gUDirection.x(), pos.y() - gUDirection.y(), pos.z() - gUDirection.z()};
0287   const std::vector<double> vec_v = {pos.x() - gVDirection.x(), pos.y() - gVDirection.y(), pos.z() - gVDirection.z()};
0288   const std::vector<double> vec_w = {pos.x() - gWDirection.x(), pos.y() - gWDirection.y(), pos.z() - gWDirection.z()};
0289 
0290   //FIXME: how to make inner vectors const?
0291   const std::vector<std::vector<double> > global_vecs = {vec_u, vec_v, vec_w};
0292 
0293   const double n = TMath::Sqrt(global_vecs.at(iParameter).at(0) * global_vecs.at(iParameter).at(0) +
0294                                global_vecs.at(iParameter).at(1) * global_vecs.at(iParameter).at(1) +
0295                                global_vecs.at(iParameter).at(2) * global_vecs.at(iParameter).at(2));
0296   const double r = TMath::Sqrt(pos.x() * pos.x() + pos.y() * pos.y());
0297 
0298   const double phase = this->getPhase(constraintparameters);
0299   //const double radial_direction[3] = {TMath::Sin(phase), TMath::Cos(phase), 0.0};
0300   const std::vector<double> radial_direction = {TMath::Sin(phase), TMath::Cos(phase), 0.0};
0301   //is equal to unity by construction ...
0302   const double norm_radial_direction =
0303       TMath::Sqrt(radial_direction.at(0) * radial_direction.at(0) + radial_direction.at(1) * radial_direction.at(1) +
0304                   radial_direction.at(2) * radial_direction.at(2));
0305 
0306   //const double phi_direction[3] = { -1.0 * pos.y(), pos.x(), 0.0};
0307   const std::vector<double> phi_direction = {-1.0 * pos.y(), pos.x(), 0.0};
0308   const double norm_phi_direction =
0309       TMath::Sqrt(phi_direction.at(0) * phi_direction.at(0) + phi_direction.at(1) * phi_direction.at(1) +
0310                   phi_direction.at(2) * phi_direction.at(2));
0311 
0312   //const double z_direction[3] = {0.0, 0.0, 1.0};
0313   static const std::vector<double> z_direction = {0.0, 0.0, 1.0};
0314   const double norm_z_direction =
0315       TMath::Sqrt(z_direction.at(0) * z_direction.at(0) + z_direction.at(1) * z_direction.at(1) +
0316                   z_direction.at(2) * z_direction.at(2));
0317 
0318   //unit vector pointing from the origin to the module position in the transverse plane
0319   const std::vector<double> rDirection = {pos.x(), pos.y(), 0.0};
0320   const double norm_rDirection = TMath::Sqrt(rDirection.at(0) * rDirection.at(0) + rDirection.at(1) * rDirection.at(1) +
0321                                              rDirection.at(2) * rDirection.at(2));
0322 
0323   double coeff = 0.0;
0324   double dot_product = 0.0;
0325   double normalisation_factor = 1.0;
0326 
0327   //see https://indico.cern.ch/getFile.py/access?contribId=15&sessionId=1&resId=0&materialId=slides&confId=127126
0328   switch (sysdeformation) {
0329     case SystematicDeformations::kTwist:
0330     case SystematicDeformations::kLayerRotation:
0331       dot_product = phi_direction.at(0) * global_vecs.at(iParameter).at(0) +
0332                     phi_direction.at(1) * global_vecs.at(iParameter).at(1) +
0333                     phi_direction.at(2) * global_vecs.at(iParameter).at(2);
0334       normalisation_factor = r * n * norm_phi_direction;
0335       break;
0336     case SystematicDeformations::kZexpansion:
0337     case SystematicDeformations::kTelescope:
0338     case SystematicDeformations::kSkew:
0339       dot_product = global_vecs.at(iParameter).at(0) * z_direction.at(0) +
0340                     global_vecs.at(iParameter).at(1) * z_direction.at(1) +
0341                     global_vecs.at(iParameter).at(2) * z_direction.at(2);
0342       normalisation_factor = (n * norm_z_direction);
0343       break;
0344     case SystematicDeformations::kRadial:
0345     case SystematicDeformations::kBowing:
0346     case SystematicDeformations::kElliptical:
0347       dot_product = global_vecs.at(iParameter).at(0) * rDirection.at(0) +
0348                     global_vecs.at(iParameter).at(1) * rDirection.at(1) +
0349                     global_vecs.at(iParameter).at(2) * rDirection.at(2);
0350       normalisation_factor = (n * norm_rDirection);
0351       break;
0352     case SystematicDeformations::kSagitta:
0353       dot_product = global_vecs.at(iParameter).at(0) * radial_direction.at(0) +
0354                     global_vecs.at(iParameter).at(1) * radial_direction.at(1) +
0355                     global_vecs.at(iParameter).at(2) * radial_direction.at(2);
0356       normalisation_factor = (n * norm_radial_direction);
0357       break;
0358     default:
0359       break;
0360   }
0361 
0362   if (TMath::Abs(normalisation_factor) > 0.0) {
0363     coeff = dot_product * (this->getX(sysdeformation, pos, phase) - x0) / normalisation_factor;
0364   } else {
0365     throw cms::Exception("Alignment") << "[PedeSteererWeakModeConstraints::getCoefficient]"
0366                                       << " Normalisation factor"
0367                                       << "for coefficient calculation equal to zero! Misconfiguration?";
0368   }
0369   return coeff;
0370 }
0371 
0372 //_________________________________________________________________________
0373 bool PedeSteererWeakModeConstraints::checkSelectionShiftParameter(const Alignable* ali, unsigned int iParameter) const {
0374   bool isselected = false;
0375   const std::vector<bool>& aliSel = ali->alignmentParameters()->selector();
0376   //exclude non-shift parameters
0377   if ((iParameter <= 2) || (iParameter >= 9 && iParameter <= 11)) {
0378     if (!aliSel.at(iParameter)) {
0379       isselected = false;
0380     } else {
0381       auto params = ali->alignmentParameters();
0382       auto selVar = dynamic_cast<SelectionUserVariables*>(params->userVariables());
0383       if (selVar) {
0384         if (selVar->fullSelection().size() <= (iParameter + 1)) {
0385           throw cms::Exception("Alignment")
0386               << "[PedeSteererWeakModeConstraints::checkSelectionShiftParameter]"
0387               << " Can not access selected alignment variables of alignable "
0388               << alignableObjectId_.idToString(ali->alignableObjectId()) << "at (" << ali->globalPosition().x() << ","
0389               << ali->globalPosition().y() << "," << ali->globalPosition().z() << ") "
0390               << "for parameter number " << (iParameter + 1) << ".";
0391         }
0392       }
0393       const char selChar = (selVar ? selVar->fullSelection().at(iParameter) : '1');
0394       // if(selChar == '1') { //FIXME??? what about 'r'?
0395       if (selChar == '1' || selChar == 'r') {
0396         isselected = true;
0397       } else {
0398         isselected = false;
0399       }
0400     }
0401   }
0402   return isselected;
0403 }
0404 
0405 //_________________________________________________________________________
0406 void PedeSteererWeakModeConstraints::closeOutputfiles() {
0407   //'delete' output files which means: close them
0408   for (auto& it : ConstraintsConfigContainer_) {
0409     for (auto& iFile : it.mapFileName_) {
0410       if (iFile.second) {
0411         delete iFile.second;
0412         iFile.second = nullptr;
0413       } else {
0414         throw cms::Exception("FileCloseProblem") << "[PedeSteererWeakModeConstraints]"
0415                                                  << " can not close file " << iFile.first << ".";
0416       }
0417     }
0418   }
0419 }
0420 
0421 //_________________________________________________________________________
0422 void PedeSteererWeakModeConstraints::writeOutput(const std::list<std::pair<unsigned int, double> >& output,
0423                                                  const GeometryConstraintConfigData& it,
0424                                                  const Alignable* iHLS,
0425                                                  double sum_xi_x0) {
0426   std::ofstream* ofile = getFile(it, iHLS);
0427 
0428   if (ofile == nullptr) {
0429     throw cms::Exception("FileFindError") << "[PedeSteererWeakModeConstraints] Cannot find output file.";
0430   } else {
0431     if (!output.empty()) {
0432       const double constr = sum_xi_x0 * it.coefficients_.front();
0433       (*ofile) << "Constraint " << std::scientific << constr << std::endl;
0434       for (const auto& ioutput : output) {
0435         (*ofile) << std::fixed << ioutput.first << " " << std::scientific << ioutput.second << std::endl;
0436       }
0437     }
0438   }
0439 }
0440 
0441 //_________________________________________________________________________
0442 std::ofstream* PedeSteererWeakModeConstraints::getFile(const GeometryConstraintConfigData& it,
0443                                                        const Alignable* iHLS) const {
0444   std::ofstream* file = nullptr;
0445 
0446   for (const auto& ilevelsFilename : it.levelsFilenames_) {
0447     if (ilevelsFilename.first->id() == iHLS->id() &&
0448         ilevelsFilename.first->alignableObjectId() == iHLS->alignableObjectId()) {
0449       const auto iFile = it.mapFileName_.find(ilevelsFilename.second);
0450       if (iFile != it.mapFileName_.end()) {
0451         file = iFile->second;
0452       }
0453     }
0454   }
0455 
0456   return file;
0457 }
0458 
0459 //_________________________________________________________________________
0460 double PedeSteererWeakModeConstraints::getX0(const std::pair<Alignable*, std::list<Alignable*> >& iHLS,
0461                                              const GeometryConstraintConfigData& it) const {
0462   double nmodules = 0.0;
0463   double x0 = 0.0;
0464 
0465   for (const auto& ali : iHLS.second) {
0466     align::PositionType pos = ali->globalPosition();
0467     bool alignableIsFloating = false;  //means: true=alignable is able to move in at least one direction
0468 
0469     //test whether at least one variable has been selected in the configuration
0470     for (unsigned int iParameter = 0; static_cast<int>(iParameter) < ali->alignmentParameters()->size(); iParameter++) {
0471       if (this->checkSelectionShiftParameter(ali, iParameter)) {
0472         alignableIsFloating = true;
0473         // //verify that alignable has just one label -- meaning no IOV-dependence etc
0474         // const unsigned int nInstances = myLabels_->numberOfParameterInstances(ali, iParameter);
0475         // if(nInstances > 1) {
0476         //   throw cms::Exception("PedeSteererWeakModeConstraints")
0477         //     << "@SUB=PedeSteererWeakModeConstraints::ConstructConstraints"
0478         //     << " Weak mode constraints are only supported for alignables which have"
0479         //     << " just one label. However, e.g. alignable"
0480         //     << " " << alignableObjectId_.idToString(ali->alignableObjectId())
0481         //     << "at (" << ali->globalPosition().x() << ","<< ali->globalPosition().y() << "," << ali->globalPosition().z()<< "), "
0482         //     << " was configured to have >1 label. Remove e.g. IOV-dependence for this (and other) alignables which are used in the constraint.";
0483         // }
0484         break;
0485       }
0486     }
0487     //at least one parameter of the alignable can be changed in the alignment
0488     if (alignableIsFloating) {
0489       const auto phase = this->getPhase(it.coefficients_);
0490       if (ali->alignmentParameters()->type() != AlignmentParametersFactory::kTwoBowedSurfaces) {
0491         x0 += this->getX(it.sysdeformation_, pos, phase);
0492         nmodules++;
0493       } else {
0494         std::pair<align::GlobalPoint, align::GlobalPoint> sensorpositions = this->getDoubleSensorPosition(ali);
0495         x0 += this->getX(it.sysdeformation_, sensorpositions.first, phase) +
0496               this->getX(it.sysdeformation_, sensorpositions.second, phase);
0497         nmodules++;
0498         nmodules++;
0499       }
0500     }
0501   }
0502   if (nmodules > 0) {
0503     x0 = x0 / nmodules;
0504   } else {
0505     throw cms::Exception("Alignment") << "@SUB=PedeSteererWeakModeConstraints::ConstructConstraints"
0506                                       << " Number of selected modules equal to zero. Check configuration!";
0507     x0 = 1.0;
0508   }
0509   return x0;
0510 }
0511 
0512 //_________________________________________________________________________
0513 unsigned int PedeSteererWeakModeConstraints::constructConstraints(const align::Alignables& alis) {
0514   //FIXME: split the code of the method into smaller pieces/submethods
0515 
0516   //create the data structures that store the alignables
0517   //for which the constraints need to be calculated and
0518   //their association to high-level structures
0519   const auto nConstraints = this->createAlignablesDataStructure();
0520 
0521   std::vector<std::list<std::pair<unsigned int, double> > > createdConstraints;
0522 
0523   //calculate constraints
0524   //loop over all constraints
0525   for (const auto& it : ConstraintsConfigContainer_) {
0526     //loop over all subdets for which constraints are determined
0527     for (const auto& iHLS : it.HLSsubdets_) {
0528       double sum_xi_x0 = 0.0;
0529       std::list<std::pair<unsigned int, double> > output;
0530 
0531       const double x0 = this->getX0(iHLS, it);
0532 
0533       for (std::list<Alignable*>::const_iterator iAlignables = iHLS.second.begin(); iAlignables != iHLS.second.end();
0534            iAlignables++) {
0535         const Alignable* ali = (*iAlignables);
0536         const auto aliLabel =
0537             myLabels_->alignableLabelFromParamAndInstance(const_cast<Alignable*>(ali), 0, it.instance_);
0538         const AlignableSurface& surface = ali->surface();
0539 
0540         const LocalPoint lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
0541 
0542         GlobalPoint gUDirection = surface.toGlobal(lUDirection), gVDirection = surface.toGlobal(lVDirection),
0543                     gWDirection = surface.toGlobal(lWDirection);
0544 
0545         const bool isDoubleSensor = ali->alignmentParameters()->type() == AlignmentParametersFactory::kTwoBowedSurfaces;
0546 
0547         const auto sensorpositions = isDoubleSensor ? this->getDoubleSensorPosition(ali)
0548                                                     : std::make_pair(ali->globalPosition(), align::PositionType());
0549 
0550         const auto& pos_sensor0 = sensorpositions.first;
0551         const auto& pos_sensor1 = sensorpositions.second;
0552         const auto phase = this->getPhase(it.coefficients_);
0553         const auto x_sensor0 = this->getX(it.sysdeformation_, pos_sensor0, phase);
0554         const auto x_sensor1 = isDoubleSensor ? this->getX(it.sysdeformation_, pos_sensor1, phase) : 0.0;
0555 
0556         sum_xi_x0 += (x_sensor0 - x0) * (x_sensor0 - x0);
0557         if (isDoubleSensor) {
0558           sum_xi_x0 += (x_sensor1 - x0) * (x_sensor1 - x0);
0559         }
0560         const int numparameterlimit = ali->alignmentParameters()->size();  //isDoubleSensor ? 18 : 3;
0561 
0562         for (int iParameter = 0; iParameter < numparameterlimit; iParameter++) {
0563           int localindex = 0;
0564           if (iParameter == 0 || iParameter == 9)
0565             localindex = 0;
0566           if (iParameter == 1 || iParameter == 10)
0567             localindex = 1;
0568           if (iParameter == 2 || iParameter == 11)
0569             localindex = 2;
0570 
0571           if ((iParameter >= 0 && iParameter <= 2) || (iParameter >= 9 && iParameter <= 11)) {
0572           } else {
0573             continue;
0574           }
0575           if (!this->checkSelectionShiftParameter(ali, iParameter)) {
0576             continue;
0577           }
0578           //do it for each 'instance' separately? -> IOV-dependence, no
0579           const auto paramLabel = myLabels_->parameterLabel(aliLabel, iParameter);
0580 
0581           const auto& pos = (iParameter <= 2) ? pos_sensor0 : pos_sensor1;
0582           //select only u,v,w
0583           if (iParameter == 0 || iParameter == 1 || iParameter == 2 || iParameter == 9 || iParameter == 10 ||
0584               iParameter == 11) {
0585             const double coeff = this->getCoefficient(
0586                 it.sysdeformation_, pos, gUDirection, gVDirection, gWDirection, localindex, x0, it.coefficients_);
0587             if (TMath::Abs(coeff) > 0.0) {
0588               //nothing
0589             } else {
0590               edm::LogWarning("PedeSteererWeakModeConstraints")
0591                   << "@SUB=PedeSteererWeakModeConstraints::getCoefficient"
0592                   << "Coefficient of alignable " << alignableObjectId_.idToString(ali->alignableObjectId()) << " at ("
0593                   << ali->globalPosition().x() << "," << ali->globalPosition().y() << "," << ali->globalPosition().z()
0594                   << ") "
0595                   << " in subdet " << alignableObjectId_.idToString(iHLS.first->alignableObjectId())
0596                   << " for parameter " << localindex << " equal to zero. This alignable is used in the constraint"
0597                   << " '" << it.constraintName_
0598                   << "'. The id is: alignable->geomDetId().rawId() = " << ali->geomDetId().rawId() << ".";
0599             }
0600             output.push_back(std::make_pair(paramLabel, coeff));
0601           }
0602         }
0603       }
0604 
0605       if (std::find(createdConstraints.begin(), createdConstraints.end(), output) != createdConstraints.end()) {
0606         // check if linearly dependent constraint exists already:
0607         auto outFile = getFile(it, iHLS.first);
0608         if (outFile == nullptr) {
0609           throw cms::Exception("FileFindError") << "[PedeSteererWeakModeConstraints] Cannot find output file.";
0610         } else {
0611           *outFile << "! The constraint for this IOV/momentum range" << std::endl
0612                    << "! has been removed because the used parameters" << std::endl
0613                    << "! are not IOV or momentum-range dependent." << std::endl;
0614         }
0615         continue;
0616       }
0617       this->writeOutput(output, it, iHLS.first, sum_xi_x0);
0618       createdConstraints.push_back(output);
0619     }
0620   }
0621   this->closeOutputfiles();
0622 
0623   return nConstraints;
0624 }
0625 
0626 //_________________________________________________________________________
0627 bool PedeSteererWeakModeConstraints::checkMother(const Alignable* const lowleveldet, const Alignable* const HLS) const {
0628   if (lowleveldet->id() == HLS->id() && lowleveldet->alignableObjectId() == HLS->alignableObjectId()) {
0629     return true;
0630   } else {
0631     if (lowleveldet->mother() == nullptr)
0632       return false;
0633     else
0634       return this->checkMother(lowleveldet->mother(), HLS);
0635   }
0636 }
0637 
0638 //_________________________________________________________________________
0639 void PedeSteererWeakModeConstraints::verifyParameterNames(const edm::ParameterSet& pset, unsigned int psetnr) const {
0640   const auto parameterNames = pset.getParameterNames();
0641   for (const auto& name : parameterNames) {
0642     if (name != "coefficients" && name != "deadmodules" && name != "constraint" && name != "steerFilePrefix" &&
0643         name != "levels" && name != "excludedAlignables" && name != "ignoredInstances" && name != "downToLowestLevel") {
0644       throw cms::Exception("BadConfig") << "@SUB=PedeSteererWeakModeConstraints::verifyParameterNames:"
0645                                         << " Unknown parameter name '" << name << "' in PSet number " << psetnr
0646                                         << ". Maybe a typo?";
0647     }
0648   }
0649 }
0650 
0651 //_________________________________________________________________________
0652 const std::vector<std::pair<Alignable*, std::string> > PedeSteererWeakModeConstraints::makeLevelsFilenames(
0653     std::set<std::string>& steerFilePrefixContainer,
0654     const align::Alignables& alis,
0655     const std::string& steerFilePrefix) const {
0656   //check whether the prefix is unique
0657   if (steerFilePrefixContainer.find(steerFilePrefix) != steerFilePrefixContainer.end()) {
0658     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints] Steering file"
0659                                       << " prefix '" << steerFilePrefix << "' already exists. Specify unique names!";
0660   } else {
0661     steerFilePrefixContainer.insert(steerFilePrefix);
0662   }
0663 
0664   std::vector<std::pair<Alignable*, std::string> > levelsFilenames;
0665   for (const auto& ali : alis) {
0666     std::stringstream n;
0667     n << steerFile_ << "_" << steerFilePrefix  //<< "_" << name
0668       << "_" << alignableObjectId_.idToString(ali->alignableObjectId()) << "_" << ali->id() << "_"
0669       << ali->alignableObjectId() << ".txt";
0670 
0671     levelsFilenames.push_back(std::make_pair(ali, n.str()));
0672   }
0673   return levelsFilenames;
0674 }
0675 
0676 //_________________________________________________________________________
0677 int PedeSteererWeakModeConstraints::verifyDeformationName(const std::string& name,
0678                                                           const std::vector<double>& coefficients) const {
0679   int sysdeformation = SystematicDeformations::kUnknown;
0680 
0681   if (name == "twist") {
0682     sysdeformation = SystematicDeformations::kTwist;
0683   } else if (name == "zexpansion") {
0684     sysdeformation = SystematicDeformations::kZexpansion;
0685   } else if (name == "sagitta") {
0686     sysdeformation = SystematicDeformations::kSagitta;
0687   } else if (name == "radial") {
0688     sysdeformation = SystematicDeformations::kRadial;
0689   } else if (name == "telescope") {
0690     sysdeformation = SystematicDeformations::kTelescope;
0691   } else if (name == "layerrotation") {
0692     sysdeformation = SystematicDeformations::kLayerRotation;
0693   } else if (name == "bowing") {
0694     sysdeformation = SystematicDeformations::kBowing;
0695   } else if (name == "skew") {
0696     sysdeformation = SystematicDeformations::kSkew;
0697   } else if (name == "elliptical") {
0698     sysdeformation = SystematicDeformations::kElliptical;
0699   }
0700 
0701   if (sysdeformation == SystematicDeformations::kUnknown) {
0702     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0703                                       << " specified configuration option '" << name << "' not known.";
0704   }
0705   if ((sysdeformation == SystematicDeformations::kSagitta || sysdeformation == SystematicDeformations::kElliptical ||
0706        sysdeformation == SystematicDeformations::kSkew) &&
0707       coefficients.size() != 2) {
0708     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0709                                       << " Excactly two parameters using the coefficient"
0710                                       << " variable have to be provided for the " << name << " constraint.";
0711   }
0712   if ((sysdeformation == SystematicDeformations::kTwist || sysdeformation == SystematicDeformations::kZexpansion ||
0713        sysdeformation == SystematicDeformations::kTelescope ||
0714        sysdeformation == SystematicDeformations::kLayerRotation || sysdeformation == SystematicDeformations::kRadial ||
0715        sysdeformation == SystematicDeformations::kBowing) &&
0716       coefficients.size() != 1) {
0717     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0718                                       << " Excactly ONE parameter using the coefficient"
0719                                       << " variable have to be provided for the " << name << " constraint.";
0720   }
0721 
0722   if (coefficients.empty()) {
0723     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0724                                       << " At least one coefficient has to be specified.";
0725   }
0726   return sysdeformation;
0727 }
0728 
0729 //_________________________________________________________________________
0730 double PedeSteererWeakModeConstraints::getPhase(const std::vector<double>& coefficients) const {
0731   return coefficients.size() == 2 ? coefficients.at(1) : 0.0;  //treat second parameter as phase otherwise return 0
0732 }
0733 
0734 //_________________________________________________________________________
0735 PedeSteererWeakModeConstraints::~PedeSteererWeakModeConstraints() = default;