Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-07-31 22:36:59

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 = myLabels_->alignableLabelFromParamAndInstance(ali, 0, it.instance_);
0537         const AlignableSurface& surface = ali->surface();
0538 
0539         const LocalPoint lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
0540 
0541         GlobalPoint gUDirection = surface.toGlobal(lUDirection), gVDirection = surface.toGlobal(lVDirection),
0542                     gWDirection = surface.toGlobal(lWDirection);
0543 
0544         const bool isDoubleSensor = ali->alignmentParameters()->type() == AlignmentParametersFactory::kTwoBowedSurfaces;
0545 
0546         const auto sensorpositions = isDoubleSensor ? this->getDoubleSensorPosition(ali)
0547                                                     : std::make_pair(ali->globalPosition(), align::PositionType());
0548 
0549         const auto& pos_sensor0 = sensorpositions.first;
0550         const auto& pos_sensor1 = sensorpositions.second;
0551         const auto phase = this->getPhase(it.coefficients_);
0552         const auto x_sensor0 = this->getX(it.sysdeformation_, pos_sensor0, phase);
0553         const auto x_sensor1 = isDoubleSensor ? this->getX(it.sysdeformation_, pos_sensor1, phase) : 0.0;
0554 
0555         sum_xi_x0 += (x_sensor0 - x0) * (x_sensor0 - x0);
0556         if (isDoubleSensor) {
0557           sum_xi_x0 += (x_sensor1 - x0) * (x_sensor1 - x0);
0558         }
0559         const int numparameterlimit = ali->alignmentParameters()->size();  //isDoubleSensor ? 18 : 3;
0560 
0561         for (int iParameter = 0; iParameter < numparameterlimit; iParameter++) {
0562           int localindex = 0;
0563           if (iParameter == 0 || iParameter == 9)
0564             localindex = 0;
0565           if (iParameter == 1 || iParameter == 10)
0566             localindex = 1;
0567           if (iParameter == 2 || iParameter == 11)
0568             localindex = 2;
0569 
0570           if ((iParameter >= 0 && iParameter <= 2) || (iParameter >= 9 && iParameter <= 11)) {
0571           } else {
0572             continue;
0573           }
0574           if (!this->checkSelectionShiftParameter(ali, iParameter)) {
0575             continue;
0576           }
0577           //do it for each 'instance' separately? -> IOV-dependence, no
0578           const auto paramLabel = myLabels_->parameterLabel(aliLabel, iParameter);
0579 
0580           const auto& pos = (iParameter <= 2) ? pos_sensor0 : pos_sensor1;
0581           //select only u,v,w
0582           if (iParameter == 0 || iParameter == 1 || iParameter == 2 || iParameter == 9 || iParameter == 10 ||
0583               iParameter == 11) {
0584             const double coeff = this->getCoefficient(
0585                 it.sysdeformation_, pos, gUDirection, gVDirection, gWDirection, localindex, x0, it.coefficients_);
0586             if (TMath::Abs(coeff) > 0.0) {
0587               //nothing
0588             } else {
0589               edm::LogWarning("PedeSteererWeakModeConstraints")
0590                   << "@SUB=PedeSteererWeakModeConstraints::getCoefficient"
0591                   << "Coefficient of alignable " << alignableObjectId_.idToString(ali->alignableObjectId()) << " at ("
0592                   << ali->globalPosition().x() << "," << ali->globalPosition().y() << "," << ali->globalPosition().z()
0593                   << ") "
0594                   << " in subdet " << alignableObjectId_.idToString(iHLS.first->alignableObjectId())
0595                   << " for parameter " << localindex << " equal to zero. This alignable is used in the constraint"
0596                   << " '" << it.constraintName_
0597                   << "'. The id is: alignable->geomDetId().rawId() = " << ali->geomDetId().rawId() << ".";
0598             }
0599             output.push_back(std::make_pair(paramLabel, coeff));
0600           }
0601         }
0602       }
0603 
0604       if (std::find(createdConstraints.begin(), createdConstraints.end(), output) != createdConstraints.end()) {
0605         // check if linearly dependent constraint exists already:
0606         auto outFile = getFile(it, iHLS.first);
0607         if (outFile == nullptr) {
0608           throw cms::Exception("FileFindError") << "[PedeSteererWeakModeConstraints] Cannot find output file.";
0609         } else {
0610           *outFile << "! The constraint for this IOV/momentum range" << std::endl
0611                    << "! has been removed because the used parameters" << std::endl
0612                    << "! are not IOV or momentum-range dependent." << std::endl;
0613         }
0614         continue;
0615       }
0616       this->writeOutput(output, it, iHLS.first, sum_xi_x0);
0617       createdConstraints.push_back(output);
0618     }
0619   }
0620   this->closeOutputfiles();
0621 
0622   return nConstraints;
0623 }
0624 
0625 //_________________________________________________________________________
0626 bool PedeSteererWeakModeConstraints::checkMother(const Alignable* const lowleveldet, const Alignable* const HLS) const {
0627   if (lowleveldet->id() == HLS->id() && lowleveldet->alignableObjectId() == HLS->alignableObjectId()) {
0628     return true;
0629   } else {
0630     if (lowleveldet->mother() == nullptr)
0631       return false;
0632     else
0633       return this->checkMother(lowleveldet->mother(), HLS);
0634   }
0635 }
0636 
0637 //_________________________________________________________________________
0638 void PedeSteererWeakModeConstraints::verifyParameterNames(const edm::ParameterSet& pset, unsigned int psetnr) const {
0639   const auto parameterNames = pset.getParameterNames();
0640   for (const auto& name : parameterNames) {
0641     if (name != "coefficients" && name != "deadmodules" && name != "constraint" && name != "steerFilePrefix" &&
0642         name != "levels" && name != "excludedAlignables" && name != "ignoredInstances" && name != "downToLowestLevel") {
0643       throw cms::Exception("BadConfig") << "@SUB=PedeSteererWeakModeConstraints::verifyParameterNames:"
0644                                         << " Unknown parameter name '" << name << "' in PSet number " << psetnr
0645                                         << ". Maybe a typo?";
0646     }
0647   }
0648 }
0649 
0650 //_________________________________________________________________________
0651 const std::vector<std::pair<Alignable*, std::string> > PedeSteererWeakModeConstraints::makeLevelsFilenames(
0652     std::set<std::string>& steerFilePrefixContainer,
0653     const align::Alignables& alis,
0654     const std::string& steerFilePrefix) const {
0655   //check whether the prefix is unique
0656   if (steerFilePrefixContainer.find(steerFilePrefix) != steerFilePrefixContainer.end()) {
0657     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints] Steering file"
0658                                       << " prefix '" << steerFilePrefix << "' already exists. Specify unique names!";
0659   } else {
0660     steerFilePrefixContainer.insert(steerFilePrefix);
0661   }
0662 
0663   std::vector<std::pair<Alignable*, std::string> > levelsFilenames;
0664   for (const auto& ali : alis) {
0665     std::stringstream n;
0666     n << steerFile_ << "_" << steerFilePrefix  //<< "_" << name
0667       << "_" << alignableObjectId_.idToString(ali->alignableObjectId()) << "_" << ali->id() << "_"
0668       << ali->alignableObjectId() << ".txt";
0669 
0670     levelsFilenames.push_back(std::make_pair(ali, n.str()));
0671   }
0672   return levelsFilenames;
0673 }
0674 
0675 //_________________________________________________________________________
0676 int PedeSteererWeakModeConstraints::verifyDeformationName(const std::string& name,
0677                                                           const std::vector<double>& coefficients) const {
0678   int sysdeformation = SystematicDeformations::kUnknown;
0679 
0680   if (name == "twist") {
0681     sysdeformation = SystematicDeformations::kTwist;
0682   } else if (name == "zexpansion") {
0683     sysdeformation = SystematicDeformations::kZexpansion;
0684   } else if (name == "sagitta") {
0685     sysdeformation = SystematicDeformations::kSagitta;
0686   } else if (name == "radial") {
0687     sysdeformation = SystematicDeformations::kRadial;
0688   } else if (name == "telescope") {
0689     sysdeformation = SystematicDeformations::kTelescope;
0690   } else if (name == "layerrotation") {
0691     sysdeformation = SystematicDeformations::kLayerRotation;
0692   } else if (name == "bowing") {
0693     sysdeformation = SystematicDeformations::kBowing;
0694   } else if (name == "skew") {
0695     sysdeformation = SystematicDeformations::kSkew;
0696   } else if (name == "elliptical") {
0697     sysdeformation = SystematicDeformations::kElliptical;
0698   }
0699 
0700   if (sysdeformation == SystematicDeformations::kUnknown) {
0701     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0702                                       << " specified configuration option '" << name << "' not known.";
0703   }
0704   if ((sysdeformation == SystematicDeformations::kSagitta || sysdeformation == SystematicDeformations::kElliptical ||
0705        sysdeformation == SystematicDeformations::kSkew) &&
0706       coefficients.size() != 2) {
0707     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0708                                       << " Excactly two parameters using the coefficient"
0709                                       << " variable have to be provided for the " << name << " constraint.";
0710   }
0711   if ((sysdeformation == SystematicDeformations::kTwist || sysdeformation == SystematicDeformations::kZexpansion ||
0712        sysdeformation == SystematicDeformations::kTelescope ||
0713        sysdeformation == SystematicDeformations::kLayerRotation || sysdeformation == SystematicDeformations::kRadial ||
0714        sysdeformation == SystematicDeformations::kBowing) &&
0715       coefficients.size() != 1) {
0716     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0717                                       << " Excactly ONE parameter using the coefficient"
0718                                       << " variable have to be provided for the " << name << " constraint.";
0719   }
0720 
0721   if (coefficients.empty()) {
0722     throw cms::Exception("BadConfig") << "[PedeSteererWeakModeConstraints]"
0723                                       << " At least one coefficient has to be specified.";
0724   }
0725   return sysdeformation;
0726 }
0727 
0728 //_________________________________________________________________________
0729 double PedeSteererWeakModeConstraints::getPhase(const std::vector<double>& coefficients) const {
0730   return coefficients.size() == 2 ? coefficients.at(1) : 0.0;  //treat second parameter as phase otherwise return 0
0731 }
0732 
0733 //_________________________________________________________________________
0734 PedeSteererWeakModeConstraints::~PedeSteererWeakModeConstraints() = default;