MillePedeAlignmentAlgorithm/src/PedeSteerer.cc

0001 /**
0002  * \file PedeSteerer.cc
0003  *
0004  *  \author    : Gero Flucke
0005  *  date       : October 2006
0006  *  $Revision: 1.42 $
0007  *  $Date: 2013/06/18 15:47:55 $
0008  *  (last update by $Author: jbehr $)
0009  */
0010
0011 #include "PedeSteerer.h"
0012 #include "PedeSteererWeakModeConstraints.h"
0013
0014 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
0015
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "Alignment/CommonAlignment/interface/Alignable.h"
0019 #include "Alignment/CommonAlignment/interface/Utilities.h"
0020 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0021 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterSelector.h"
0022 #include "Alignment/CommonAlignmentAlgorithm/interface/SelectionUserVariables.h"
0023 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
0024 #include "Alignment/CommonAlignmentParametrization/interface/RigidBodyAlignmentParameters.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 <atomic>
0042 #include <fstream>
0043 #include <sstream>
0044 #include <algorithm>
0045
0046 // from ROOT
0047 #include <TSystem.h>
0048 #include <TMath.h>
0049
0050 #include <iostream>
0051
0052 //_________________________________________________________________________
0053 PedeSteerer::PedeSteerer(AlignableTracker *aliTracker,
0054                          AlignableMuon *aliMuon,
0055                          AlignableExtras *aliExtras,
0056                          AlignmentParameterStore *store,
0057                          const PedeLabelerBase *labels,
0058                          const edm::ParameterSet &config,
0059                          const std::string &defaultDir,
0060                          bool noSteerFiles)
0061     : myParameterStore(store),
0062       myLabels(labels),
0063       alignableObjectId_{AlignableObjectId::commonObjectIdProvider(aliTracker, aliMuon)},
0064       myConfig(config),
0065       myDirectory(myConfig.getUntrackedParameter<std::string>("fileDir")),
0066       myRunDirectory(myConfig.getUntrackedParameter<std::string>("runDir")),
0067       myNoSteerFiles(noSteerFiles),
0068       myIsSteerFileDebug(myConfig.getUntrackedParameter<bool>("steerFileDebug")),
0069       myParameterSign(myConfig.getUntrackedParameter<int>("parameterSign")),
0070       theMinHieraConstrCoeff(myConfig.getParameter<double>("minHieraConstrCoeff")),
0071       theMinHieraParPerConstr(myConfig.getParameter<unsigned int>("minHieraParPerConstr")),
0072       theConstrPrecision(myConfig.getParameter<unsigned int>("constrPrecision")),
0073       theCoordMaster(nullptr) {
0074   if (myParameterSign != 1 && myParameterSign != -1) {
0075     cms::Exception("BadConfig") << "Expect PedeSteerer.parameterSign = +/-1, "
0076                                 << "found " << myParameterSign << ".";
0077   }
0078
0079   // Correct directory, needed before asking for fileName(..):
0080   if (myDirectory.empty())
0081     myDirectory = defaultDir;
0082   if (!myDirectory.empty() && myDirectory.find_last_of('/') != myDirectory.size() - 1) {
0083     myDirectory += '/';  // directory may need '/'
0084   }
0085
0086   if (myRunDirectory.empty())
0087     myRunDirectory = defaultDir;
0088   if (!myRunDirectory.empty() && myRunDirectory.find_last_of('/') != myRunDirectory.size() - 1) {
0089     myRunDirectory += '/';  // directory may need '/'
0090     myDirectory = myRunDirectory;
0091   }
0092
0093   const auto &alis = myParameterStore->alignables();
0094   if (!this->checkParameterChoices(alis)) {
0095   }  // anyway thrown exception
0096
0097   // Coordinate system selection and correction before everything
0098   theCoordDefiners = this->selectCoordinateAlis(alis);
0099   if (!theCoordDefiners.empty()) {  // Create steering with constraints to define coordinate system:
0100     // OK, some hacks:
0101     // - we want a composite with global coordinates where tracker and muon are components
0102     //   (to call RigidBodyAl.Param.->globalParameters() in correctToReferenceSystem(..))
0103     // - so we create a AlignableComposite and add tracker and muon
0104     // - but the addComponent(..) method is so cute that it calculates position from
0105     //   daughters' deepComponents()
0106     // - so we want to move it back to (0,0,0), but ali->move(..) would move daughters as well
0107     //   => move only the surface back
0108     // - this master object does not have a label for its parameters
0109     //   => some warnings if debug output selected in pedeSteer files
0110     // - we must not delete our new master (little mem. leak...) since that would delete
0111     //   the daughters as well!
0112     if (aliTracker) {
0113       theCoordMaster = new AlignableComposite(aliTracker->id(), align::invalid);
0114       theCoordMaster->addComponent(aliTracker);
0115     } else if (aliMuon) {
0116       theCoordMaster = new AlignableComposite(aliMuon->id(), align::invalid);
0117     } else {
0118       throw cms::Exception("BadConfig") << "[PedeSteerer]"
0119                                         << "Cannot define global coordinate system "
0120                                         << "with neither tracker nor muon!";
0121     }
0122     if (aliMuon)
0123       theCoordMaster->addComponent(aliMuon);  // tracker is already added if existing
0124     if (aliExtras) {                          // tracker and/or muon are already added if existing
0125       align::Alignables allExtras = aliExtras->components();
0126       for (const auto &it : allExtras)
0127         theCoordMaster->addComponent(it);
0128     }
0129
0130     theCoordMaster->recenterSurface();
0131
0132     if (this->isCorrectToRefSystem(theCoordDefiners)) {  // defined by 's' (MC): 'correct' misalignment
0133       this->correctToReferenceSystem();                  // really before 'defineCoordinates'?
0134     }
0135   }
0136 }
0137
0138 //___________________________________________________________________________
0139 PedeSteerer::~PedeSteerer() {
0140   // delete theCoordMaster; NO, see above
0141 }
0142
0143 //_________________________________________________________________________
0144 bool PedeSteerer::isNoHiera(const Alignable *ali) const {
0145   return (myNoHieraCollection.find(ali) != myNoHieraCollection.end());
0146 }
0147
0148 //_________________________________________________________________________
0149 double PedeSteerer::cmsToPedeFactor(unsigned int parNum) const {
0150   return 1.;  // mmh, otherwise would need to FIXME hierarchyConstraint...
0151
0152   switch (parNum) {
0153     case RigidBodyAlignmentParameters::dx:
0154     case RigidBodyAlignmentParameters::dy:
0155       return 1000.;  // cm to mum *1/10 to get smaller values
0156     case RigidBodyAlignmentParameters::dz:
0157       return 2500.;  // cm to mum *1/4
0158     case RigidBodyAlignmentParameters::dalpha:
0159     case RigidBodyAlignmentParameters::dbeta:
0160       return 1000.;  // rad to mrad (no first guess for sensitivity yet)
0161     case RigidBodyAlignmentParameters::dgamma:
0162       return 10000.;  // rad to mrad *10 to get larger values
0163     default:
0164       return 1.;
0165   }
0166 }
0167
0168 //_________________________________________________________________________
0169 unsigned int PedeSteerer::buildNoHierarchyCollection(const align::Alignables &alis) {
0170   myNoHieraCollection.clear();  // just in case of re-use...
0171
0172   for (const auto &iAli : alis) {
0173     AlignmentParameters *params = iAli->alignmentParameters();
0174     SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(params->userVariables());
0175     if (!selVar)
0176       continue;
0177     // Now check whether taking out of hierarchy is selected - must be consistent!
0178     unsigned int numNoHieraPar = 0;
0179     unsigned int numHieraPar = 0;
0180     for (unsigned int iParam = 0; static_cast<int>(iParam) < params->size(); ++iParam) {
0181       const char selector = selVar->fullSelection()[iParam];
0182       if (selector == 'C' || selector == 'F' || selector == 'H') {
0183         ++numNoHieraPar;
0184       } else if (selector == 'c' || selector == 'f' || selector == '1' || selector == 'r' || selector == 's') {
0185         ++numHieraPar;
0186       }  // else ... accept '0' as undetermined
0187     }
0188     if (numNoHieraPar) {  // Selected to be taken out.
0189       if (numHieraPar) {  // Inconsistent: Some parameters still in hierarchy ==> exception!
0190         throw cms::Exception("BadConfig")
0191             << "[PedeSteerer::buildNoHierarchyCollection] All active parameters of alignables to be "
0192             << " taken out of the hierarchy must be marked with capital letters 'C', 'F' or 'H'!";
0193       }
0194       bool isInHiera = false;  // Check whether Alignable is really part of hierarchy:
0195       auto mother = iAli;
0196       while ((mother = mother->mother())) {
0197         if (mother->alignmentParameters())
0198           isInHiera = true;  // could 'break;', but loop is short
0199       }
0200       // Complain, but keep collection short if not in hierarchy:
0201       if (isInHiera)
0202         myNoHieraCollection.insert(iAli);
0203       else
0204         edm::LogWarning("Alignment") << "@SUB=PedeSteerer::buildNoHierarchyCollection"
0205                                      << "Alignable not in hierarchy, no need to remove it!";
0206     }
0207   }  // end loop on alignables
0208
0209   return myNoHieraCollection.size();
0210 }
0211
0212 //_________________________________________________________________________
0213 bool PedeSteerer::checkParameterChoices(const align::Alignables &alis) const {
0214   for (const auto &iAli : alis) {
0215     AlignmentParameters *paras = iAli->alignmentParameters();
0216     SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(paras->userVariables());
0217     if (!selVar)
0218       continue;
0219     for (unsigned int iParam = 0; static_cast<int>(iParam) < paras->size(); ++iParam) {
0220       const char sel = selVar->fullSelection()[iParam];
0221       if (sel != 'f' && sel != 'F' && sel != 'c' && sel != 'C' && sel != '0' && sel != '1' && sel != 'H' &&
0222           sel != 'r' && sel != 's') {
0223         throw cms::Exception("BadConfig")
0224             << "[PedeSteerer::unknownParameterChoices] "
0225             << "Unexpected parameter selector '" << sel
0226             << "', use \n'f/F' (fix),\n'c/C' (fix at correct pos.),\n'1/H' (free),\n"
0227             << "'r/s' (free, but defining reference system, trying to correct misalignment if 's')"
0228             << " or \n'0' (ignore).\n"
0229             << "Capital letters mean that the Alignable is taken out of a possible hierarchy,\n"
0230             << "but must be used consistently for all its parameters.";
0231         return false;  // unreached
0232       }
0233     }
0234   }
0235
0236   return true;
0237 }
0238
0239 //_________________________________________________________________________
0240 std::pair<unsigned int, unsigned int> PedeSteerer::fixParameters(const align::Alignables &alis,
0241                                                                  const std::string &fileName) {
0242   // return number of parameters fixed at 0. and fixed at original position
0243   std::pair<unsigned int, unsigned int> numFixNumFixCor(0, 0);
0244
0245   std::ofstream *filePtr = nullptr;
0246
0247   for (const auto &iAli : alis) {
0248     AlignmentParameters *params = iAli->alignmentParameters();
0249     SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(params->userVariables());
0250     if (!selVar)
0251       continue;
0252
0253     for (unsigned int iParam = 0; static_cast<int>(iParam) < params->size(); ++iParam) {
0254       const unsigned int nInstances = myLabels->numberOfParameterInstances(iAli, iParam);
0255       for (unsigned int iInstance = 0; iInstance < nInstances; ++iInstance) {
0256         int whichFix = this->fixParameter(iAli, iInstance, iParam, selVar->fullSelection()[iParam], filePtr, fileName);
0257         if (whichFix == 1) {
0258           ++(numFixNumFixCor.first);
0259         } else if (whichFix == -1) {
0260           ++(numFixNumFixCor.second);
0261         }
0262       }
0263     }
0264   }
0265
0266   delete filePtr;  // automatically flushes, no problem if NULL ptr.
0267
0268   return numFixNumFixCor;
0269 }
0270
0271 //_________________________________________________________________________
0272 int PedeSteerer::fixParameter(Alignable *ali,
0273                               unsigned int iInstance,
0274                               unsigned int iParam,
0275                               char selector,
0276                               std::ofstream *&filePtr,
0277                               const std::string &fileName) {
0278   int result = 0;
0279   float fixAt = 0.;
0280   if (selector == 'c' || selector == 'C') {
0281     if (ali->alignmentParameters()->type() != AlignmentParametersFactory::kRigidBody) {
0282       throw cms::Exception("BadConfig")
0283           << "PedeSteerer::fixParameter: correction (c/C) possible only for RigidBodyParameters";
0284     }
0285     fixAt = -this->parameterSign() * RigidBodyAlignmentParameters(ali, true).parameters()[iParam];
0286     result = -1;
0287   } else if (selector == 'f' || selector == 'F') {
0288     result = 1;
0289   }
0290
0291   if (result) {
0292     if (!filePtr) {
0293       filePtr = this->createSteerFile(fileName, true);
0294       (*filePtr) << "Parameter\n";
0295     }
0296     std::ofstream &file = *filePtr;
0297
0298     const unsigned int aliLabel = myLabels->alignableLabelFromParamAndInstance(ali, iParam, iInstance);
0299     file << myLabels->parameterLabel(aliLabel, iParam) << "  " << fixAt * this->cmsToPedeFactor(iParam) << " -1.0";
0300     if (myIsSteerFileDebug) {  // debug
0301       const GlobalPoint position(ali->globalPosition());
0302       file << "  * id " << ali->id() << ", eta " << position.eta() << ", z " << position.z() << ", r "
0303            << position.perp() << ", phi " << position.phi();
0304     }
0305     file << "\n";
0306   }
0307
0308   return result;
0309 }
0310
0311 //_________________________________________________________________________
0312 align::Alignables PedeSteerer::selectCoordinateAlis(const align::Alignables &alis) const {
0313   align::Alignables coordAlis;
0314
0315   for (const auto &iAli : alis) {
0316     AlignmentParameters *params = iAli->alignmentParameters();
0317     SelectionUserVariables *selVar = dynamic_cast<SelectionUserVariables *>(params->userVariables());
0318     if (!selVar)
0319       continue;
0320     unsigned int refParam = 0;
0321     unsigned int nonRefParam = 0;
0322     for (unsigned int iParam = 0; static_cast<int>(iParam) < params->size(); ++iParam) {
0323       const char selector = selVar->fullSelection()[iParam];
0324       if (selector == 'r' || selector == 's') {
0325         ++refParam;
0326       } else if (selector != '0' && selector != 'f') {  // allow also 'c'?
0327         ++nonRefParam;
0328       }
0329     }
0330     // Check whether some 'r/s' selection string. If yes and selection makes sense, add to result:
0331     if (refParam) {
0332       if (nonRefParam) {
0333         throw cms::Exception("BadConfig")
0334             << "[PedeSteerer::selectCoordinateAlis] All active parameters of alignables defining "
0335             << "the coordinate system must be marked with 'r/s' (or fixed, 'f')!";
0336       } else {
0337         auto mother = iAli;
0338         while ((mother = mother->mother())) {
0339           if (mother->alignmentParameters()) {
0340             throw cms::Exception("BadConfig") << "[PedeSteerer::selectCoordinateAlis] "
0341                                               << "Alignables defining the coordinate system must "
0342                                               << "be highest level!";
0343           }
0344         }
0345         coordAlis.push_back(iAli);
0346       }
0347     }
0348   }  // end loop on alignables
0349
0350   return coordAlis;
0351 }
0352
0353 //_________________________________________________________________________
0354 void PedeSteerer::defineCoordinates(const align::Alignables &alis, Alignable *aliMaster, const std::string &fileName) {
0355   std::ofstream *filePtr = this->createSteerFile(fileName, true);
0356   (*filePtr) << "* Constraints to define coordinate system:\n";
0357   if (!aliMaster || aliMaster->alignmentParameters()) {
0358     throw cms::Exception("BadConfig") << "[PedeSteerer::defineCoordinates] "
0359                                       << "No master alignable or it has parameters!";
0360   }
0361   if (myIsSteerFileDebug) {  // See constructor comments about hack:
0362     edm::LogError("Alignment") << "@SUB=PedeSteerer::defineCoordinates"
0363                                << "Ignore following LogicErrors from PedeLabeler.";
0364   }
0365   AlignmentParameters *par = new RigidBodyAlignmentParameters(aliMaster, false);
0366   aliMaster->setAlignmentParameters(par);  // hierarchyConstraint needs parameters
0367   this->hierarchyConstraint(aliMaster, alis, *filePtr);
0368   aliMaster->setAlignmentParameters(nullptr);  // erase dummy parameters
0369
0370   delete filePtr;  // automatically flushes, no problem if NULL ptr.
0371 }
0372
0373 //_________________________________________________________________________
0374 bool PedeSteerer::isCorrectToRefSystem(const align::Alignables &coordDefiners) const {
0375   bool doCorrect = false;
0376   bool doNotCorrect = false;
0377   for (const auto &it : coordDefiners) {
0378     SelectionUserVariables *selVar =
0379         (it->alignmentParameters() ? dynamic_cast<SelectionUserVariables *>(it->alignmentParameters()->userVariables())
0380                                    : nullptr);
0381     if (!selVar)
0382       continue;  // is an error!?
0383
0384     for (unsigned int i = 0; i < selVar->fullSelection().size(); ++i) {
0385       if (selVar->fullSelection()[i] == 'r')
0386         doNotCorrect = true;
0387       else if (selVar->fullSelection()[i] == 's')
0388         doCorrect = true;
0389     }
0390   }
0391
0392   if (doCorrect && doNotCorrect) {
0393     throw cms::Exception("BadConfig")
0394         << "[PedeSteerer::doCorrectToRefSystem]: Parameter selection 's' and 'r' must not coexist!";
0395   }
0396
0397   return doCorrect;
0398 }
0399
0400 //_________________________________________________________________________
0401 void PedeSteerer::correctToReferenceSystem() {
0402   typedef RigidBodyAlignmentParameters RbPars;
0403   if (!theCoordMaster || theCoordDefiners.empty())
0404     return;  // nothing was defined
0405
0406   align::Alignables definerDets;  // or ...DetUnits
0407   for (const auto &it :
0408        theCoordDefiners) {  // find lowest level objects of alignables that define the coordinate system
0409     const auto &comp = it->deepComponents();
0410     definerDets.insert(definerDets.end(), comp.begin(), comp.end());
0411   }
0412
0413   for (unsigned int iLoop = 0;; ++iLoop) {  // iterate: shifts and rotations are not independent
0414     AlgebraicVector meanPars(RbPars::N_PARAM);
0415     for (const auto &it : definerDets) {                // sum up mean displacements/misrotations:
0416       meanPars += RbPars(it, true).globalParameters();  // requires theCoordMaster has global frame
0417     }
0418     meanPars /= definerDets.size();
0419     const align::Scalar squareSum = meanPars.normsq();
0420
0421     if (squareSum < 1.e-20)
0422       break;  // sqrt(1.e-20)=1.e-10: close enough to stop iterating
0423     if (iLoop == 0) {
0424       edm::LogInfo("Alignment") << "@SUB=PedeSteerer::correctToReferenceSystem"
0425                                 << "Loop " << iLoop << " "
0426                                 << "Mean misalignment of dets of defined coordinate system"
0427                                 << (squareSum < 1.e-20 ? ":" : " (will be iteratively corrected to < 1.e-10):")
0428                                 << meanPars;
0429     }
0430     if (iLoop >= 5) {  // 3 iterations should be safe, use 5 for 'more' safety...
0431       edm::LogError("Alignment") << "@SUB=PedeSteerer::correctToReferenceSystem"
0432                                  << "No convergence in " << iLoop << " iterations, "
0433                                  << "remaining misalignment: " << meanPars;
0434       break;
0435     }
0436
0437     const GlobalVector globalShift(meanPars[RbPars::dx], meanPars[RbPars::dy], meanPars[RbPars::dz]);
0438     theCoordMaster->move(-globalShift);  // sign to revert
0439     align::EulerAngles globalAngles(3);
0440     globalAngles[0] = meanPars[RbPars::dalpha];
0441     globalAngles[1] = meanPars[RbPars::dbeta];
0442     globalAngles[2] = meanPars[RbPars::dgamma];
0443     theCoordMaster->rotateInGlobalFrame(align::toMatrix(-globalAngles));  // sign to revert
0444   }
0445 }
0446
0447 //_________________________________________________________________________
0448 unsigned int PedeSteerer::hierarchyConstraints(const align::Alignables &alis, const std::string &fileName) {
0449   std::ofstream *filePtr = nullptr;
0450
0451   unsigned int nConstraints = 0;
0452   align::Alignables aliDaughts;
0453   for (const auto &iA : alis) {
0454     aliDaughts.clear();
0455     if (!(iA->firstCompsWithParams(aliDaughts))) {
0456       edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraints"
0457                                    << "Some but not all daughters of "
0458                                    << alignableObjectId_.idToString(iA->alignableObjectId()) << " with params!";
0459     }
0460     //     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::hierarchyConstraints"
0461     //                << aliDaughts.size() << " ali param components";
0462     if (aliDaughts.empty())
0463       continue;
0464     //     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::hierarchyConstraints"
0465     //                << aliDaughts.size() << " alignable components ("
0466     //                << iA->size() << " in total) for "
0467     //                << aliId.alignableTypeName(iA)
0468     //                << ", layer " << aliId.typeAndLayerFromAlignable(iA).second
0469     //                << ", position " << (iA)->globalPosition()
0470     //                << ", r = " << (iA)->globalPosition().perp();
0471     if (!filePtr)
0472       filePtr = this->createSteerFile(fileName, true);
0473     ++nConstraints;
0474     this->hierarchyConstraint(iA, aliDaughts, *filePtr);
0475   }
0476
0477   delete filePtr;  // automatically flushes, no problem if NULL ptr.
0478
0479   return nConstraints;
0480 }
0481 //_________________________________________________________________________
0482 void PedeSteerer::hierarchyConstraint(const Alignable *ali,
0483                                       const align::Alignables &components,
0484                                       std::ofstream &file) const {
0485   typedef AlignmentParameterStore::ParameterId ParameterId;
0486
0487   std::vector<std::vector<ParameterId> > paramIdsVec;
0488   std::vector<std::vector<double> > factorsVec;
0489   const bool allConstr = false;  // true; // make configurable?
0490   static std::atomic<bool> allConstrWarning{false};
0491   bool expected{false};
0492   if (allConstr && allConstrWarning.compare_exchange_strong(expected, true)) {
0493     edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraint"
0494                                  << "changed to use all 6 constraints";
0495   }
0496   if (!myParameterStore->hierarchyConstraints(
0497           ali, components, paramIdsVec, factorsVec, allConstr, theMinHieraConstrCoeff)) {
0498     edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraint"
0499                                  << "Problems from store.";
0500   }
0501
0502   for (unsigned int iConstr = 0; iConstr < paramIdsVec.size(); ++iConstr) {
0503     std::ostringstream aConstr;
0504
0505     const std::vector<ParameterId> &parIds = paramIdsVec[iConstr];
0506     const std::vector<double> &factors = factorsVec[iConstr];
0507     unsigned int nParPerConstr = 0;  // keep track of used factor/parId pair
0508     // parIds.size() == factors.size() granted by myParameterStore->hierarchyConstraints
0509     for (unsigned int iParam = 0; iParam < parIds.size(); ++iParam) {
0510       Alignable *aliSubComp = parIds[iParam].first;
0511       const unsigned int compParNum = parIds[iParam].second;
0512       if (this->isNoHiera(aliSubComp)) {
0513         if (myIsSteerFileDebug)
0514           aConstr << "* Taken out of hierarchy: ";
0515         continue;
0516       }
0517       const unsigned int aliLabel = myLabels->alignableLabel(aliSubComp);
0518       const unsigned int paramLabel = myLabels->parameterLabel(aliLabel, compParNum);
0519       // FIXME: multiply by cmsToPedeFactor(subcomponent)/cmsToPedeFactor(mother) (or vice a versa?)
0520       if (theConstrPrecision > 0)
0521         aConstr << paramLabel << "    " << std::setprecision(theConstrPrecision) << factors[iParam];
0522       else
0523         aConstr << paramLabel << "    " << factors[iParam];
0524       if (myIsSteerFileDebug) {  // debug
0525         aConstr << "   ! for param " << compParNum << " of a "
0526                 << alignableObjectId_.idToString(aliSubComp->alignableObjectId()) << " at "
0527                 << aliSubComp->globalPosition() << ", r=" << aliSubComp->globalPosition().perp();
0528       }
0529       aConstr << "\n";
0530       ++nParPerConstr;  // OK, we used one.
0531     }                   // end loop on params
0532
0533     //
0534     if (nParPerConstr && nParPerConstr >= theMinHieraParPerConstr) {  // Enough to make sense?
0535       if (myIsSteerFileDebug) {                                       //debug
0536         file << "\n* Nr. " << iConstr << " of a '" << alignableObjectId_.idToString(ali->alignableObjectId())
0537              << "' (label " << myLabels->alignableLabel(const_cast<Alignable *>(ali))  // ugly cast: FIXME!
0538              << "), position " << ali->globalPosition() << ", r = " << ali->globalPosition().perp();
0539       }
0540       file << "\nConstraint   0.\n" << aConstr.str();  // in future 'Wconstraint'?
0541     } else if (nParPerConstr > 0) {                    // no warning for trivial case...
0542       edm::LogWarning("Alignment") << "@SUB=PedeSteerer::hierarchyConstraint"
0543                                    << "Skip constraint on " << nParPerConstr << " parameter(s):\n"
0544                                    << aConstr.str();
0545     }
0546   }  // end loop on constraints
0547 }
0548
0549 //_________________________________________________________________________
0550 unsigned int PedeSteerer::presigmas(const std::vector<edm::ParameterSet> &cffPresi,
0551                                     const std::string &fileName,
0552                                     const align::Alignables &alis,
0553                                     AlignableTracker *aliTracker,
0554                                     AlignableMuon *aliMuon,
0555                                     AlignableExtras *aliExtras) {
0556   // We loop on given PSet's, each containing a parameter selection and the presigma value
0557   // The resulting presigmas are stored in a map with Alignable* as key.
0558   // This map, 'fileName' and 'alis' are passed further to create the steering file.
0559
0560   AlignmentParameterSelector selector(aliTracker, aliMuon, aliExtras);
0561   AlignablePresigmasMap aliPresiMap;  // map to store alis with presigmas of their parameters
0562   for (std::vector<edm::ParameterSet>::const_iterator iSet = cffPresi.begin(), iE = cffPresi.end(); iSet != iE;
0563        ++iSet) {  // loop on individual PSets defining ali-params with their presigma
0564     selector.clear();
0565     selector.addSelections((*iSet).getParameter<edm::ParameterSet>("Selector"));
0566     const auto &alis = selector.selectedAlignables();
0567     const std::vector<std::vector<char> > &sels = selector.selectedParameters();
0568     const float presigma = (*iSet).getParameter<double>("presigma");
0569     if (presigma <= 0.) {  // given presigma > 0., 0. later used if not (yet) chosen for parameter
0570       throw cms::Exception("BadConfig") << "[PedeSteerer::presigmas]: Pre-sigma must be > 0., but is " << presigma
0571                                         << ".";
0572     }
0573     // now loop on alis of present selection
0574     for (unsigned int iAli = 0; iAli < alis.size(); ++iAli) {
0575       std::vector<float> &presigmas = aliPresiMap[alis[iAli]];  // existing or empty, so ensure length:
0576       if (presigmas.size() < sels[iAli].size())
0577         presigmas.resize(sels[iAli].size(), 0.);
0578       for (unsigned int iParam = 0; iParam < sels[iAli].size(); ++iParam) {  // loop on parameters
0579         if (sels[iAli][iParam] != '0') {  // all but '0' means to apply the chosen presigma
0580           if (presigmas[iParam] != 0.) {  // reset forbidden (would make it order dependent!)
0581             throw cms::Exception("BadConfig")
0582                 << "[PedeSteerer::presigmas]: Try to set pre-sigma " << presigma << ", but already "
0583                 << "set " << presigmas[iParam] << " (for a "
0584                 << alignableObjectId_.idToString(alis[iAli]->alignableObjectId()) << ").";
0585           }
0586           presigmas[iParam] = presigma;
0587         }  // end if selected for presigma
0588       }    // end loop on params
0589     }      // end loop on alignables for given selection and presigma
0590   }        // end loop on PSets
0591
0592   if (aliPresiMap.empty())
0593     return 0;
0594   else
0595     return this->presigmasFile(fileName, alis, aliPresiMap);
0596 }
0597
0598 //_________________________________________________________________________
0599 unsigned int PedeSteerer::presigmasFile(const std::string &fileName,
0600                                         const align::Alignables &alis,
0601                                         const AlignablePresigmasMap &aliPresiMap) {
0602   // Check if 'alis' are in aliPresiMap,
0603   // if yes apply presigma - but NOT if parameter is fixed!
0604   std::ofstream *filePtr = nullptr;
0605
0606   unsigned int nPresiParam = 0;
0607   for (const auto &iAli : alis) {
0608     // Any presigma chosen for alignable?
0609     AlignablePresigmasMap::const_iterator presigmasIt = aliPresiMap.find(iAli);
0610     if (presigmasIt == aliPresiMap.end())
0611       continue;  // no presigma chosen for alignable
0612
0613     // Why does the following not work? It does with CMSSW_1_3_X on SLC3...
0614     // const AlignablePresigmasMap::data_type &presigmas = presigmasIt->second;
0615     const std::vector<float> &presigmas = presigmasIt->second;  // I want to hide float or double...
0616     for (unsigned int iParam = 0; iParam < presigmas.size(); ++iParam) {
0617       // Now check whether a presigma value > 0. chosen:
0618       if (presigmas[iParam] <= 0.)
0619         continue;  // must be positive, '<' checked above
0620       // Do not apply presigma to inactive or fixed values.
0621       if (!(iAli->alignmentParameters()->selector()[iParam]))
0622         continue;
0623       SelectionUserVariables *selVar =
0624           dynamic_cast<SelectionUserVariables *>(iAli->alignmentParameters()->userVariables());
0625       const char selChar = (selVar ? selVar->fullSelection()[iParam] : '1');
0626       if (selChar == 'f' || selChar == 'F' || selChar == 'c' || selChar == 'C')
0627         continue;
0628       // Finally create and write steering file:
0629       if (!filePtr) {
0630         filePtr = this->createSteerFile(fileName, true);
0631         (*filePtr) << "* Presigma values for active parameters: \nParameter\n";
0632       }
0633       const unsigned int aliLabel = myLabels->alignableLabel(iAli);
0634       (*filePtr) << myLabels->parameterLabel(aliLabel, iParam) << "   0.   "
0635                  << presigmas[iParam] * fabs(this->cmsToPedeFactor(iParam));
0636       if (myIsSteerFileDebug) {
0637         (*filePtr) << "  ! for a " << alignableObjectId_.idToString((iAli)->alignableObjectId());
0638       }
0639       (*filePtr) << '\n';
0640
0641       ++nPresiParam;
0642     }  // end loop on parameters for alignables with chosen presigmas
0643   }    // end loop on alignables
0644
0645   delete filePtr;  // close properly file
0646   return nPresiParam;
0647 }
0648
0649 //_________________________________________________________________________
0650 std::ofstream *PedeSteerer::createSteerFile(const std::string &name, bool addToList) {
0651   const std::string realName(myNoSteerFiles ? "/dev/null" : name.c_str());
0652
0653   std::ofstream *result = new std::ofstream(realName.c_str(), std::ios::out);
0654   if (!result || !result->is_open()) {
0655     delete result;  // needed before exception in case just open failed
0656     throw cms::Exception("FileOpenProblem") << "[PedeSteerer::createSteerFile]"
0657                                             << "Could not open " << realName << " as output file.";
0658   } else if (addToList) {
0659     mySteeringFiles.push_back(realName);  // keep track
0660   }
0661
0662   return result;
0663 }
0664
0665 //_________________________________________________________________________
0666 std::string PedeSteerer::fileName(const std::string &addendum) const {
0667   std::string name(myDirectory);
0668   name += myConfig.getParameter<std::string>("steerFile");
0669   name += addendum;
0670   name += ".txt";
0671
0672   return name;
0673 }
0674
0675 //___________________________________________________________________________
0676 void PedeSteerer::buildSubSteer(AlignableTracker *aliTracker, AlignableMuon *aliMuon, AlignableExtras *aliExtras) {
0677   const auto &alis = myParameterStore->alignables();
0678
0679   if (theCoordMaster && !theCoordDefiners.empty()) {
0680     const std::string nameCoordFile(this->fileName("Coord"));
0681     this->defineCoordinates(theCoordDefiners, theCoordMaster, nameCoordFile);
0682     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer" << theCoordDefiners.size()
0683                               << " highest level objects define the "
0684                               << "coordinate system, steering file " << nameCoordFile << ".";
0685   }
0686
0687   const std::string nameFixFile(this->fileName("FixPara"));
0688   const std::pair<unsigned int, unsigned int> nFixFixCor(this->fixParameters(alis, nameFixFile));
0689   if (nFixFixCor.first != 0 || nFixFixCor.second != 0) {
0690     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer" << nFixFixCor.first << " parameters fixed at 0. and "
0691                               << nFixFixCor.second << " at 'original' position, "
0692                               << "steering file " << nameFixFile << ".";
0693   }
0694
0695   if (this->buildNoHierarchyCollection(alis)) {  // before hierarchyConstraints(..)
0696     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer" << myNoHieraCollection.size()
0697                               << " alignables taken out of hierarchy.";
0698   }
0699
0700   const std::string nameHierarchyFile(this->fileName("Hierarchy"));
0701   unsigned int nConstraint = this->hierarchyConstraints(alis, nameHierarchyFile);
0702   if (nConstraint) {
0703     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer"
0704                               << "Hierarchy constraints for " << nConstraint << " alignables, "
0705                               << "steering file " << nameHierarchyFile << ".";
0706   }
0707
0708   //construct the systematic geometry deformations
0709   if (!(myConfig.getParameter<std::vector<edm::ParameterSet> >("constraints")).empty()) {
0710     PedeSteererWeakModeConstraints GeometryConstraints(
0711         aliTracker,
0712         myLabels,
0713         myConfig.getParameter<std::vector<edm::ParameterSet> >("constraints"),
0714         myConfig.getParameter<std::string>("steerFile"));
0715
0716     //prepare the output files
0717     //Get the data structure in which the configuration data are stored.
0718     //The relation between the ostream* and the corresponding file name needs to be filled
0719     auto &ConstraintsConfigContainer = GeometryConstraints.getConfigData();
0720
0721     //loop over all configured constraints
0722     for (auto &it : ConstraintsConfigContainer) {
0723       //each level has its own constraint which means the output is stored in a separate file
0724       for (const auto &ilevelsFilename : it.levelsFilenames_) {
0725         it.mapFileName_.insert(
0726             std::make_pair(ilevelsFilename.second, this->createSteerFile(ilevelsFilename.second, true)));
0727       }
0728     }
0729
0730     unsigned int nGeometryConstraint = GeometryConstraints.constructConstraints(alis);
0731     if (nGeometryConstraint) {
0732       edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer"
0733                                 << "Geometry constraints for " << nGeometryConstraint << " alignables.";
0734     }
0735   }
0736
0737   const std::string namePresigmaFile(this->fileName("Presigma"));
0738   unsigned int nPresigma = this->presigmas(myConfig.getParameter<std::vector<edm::ParameterSet> >("Presigmas"),
0739                                            namePresigmaFile,
0740                                            alis,
0741                                            aliTracker,
0742                                            aliMuon,
0743                                            aliExtras);
0744   if (nPresigma) {
0745     edm::LogInfo("Alignment") << "@SUB=PedeSteerer::buildSubSteer"
0746                               << "Presigma values set for " << nPresigma << " parameters, "
0747                               << "steering file " << namePresigmaFile << ".";
0748   }
0749
0750   // Delete all SelectionUserVariables now? They will anyway be overwritten by MillePedeVariables...
0751 }
0752
0753 //_________________________________________________________________________
0754 std::string PedeSteerer::buildMasterSteer(const std::vector<std::string> &binaryFiles) {
0755   const std::string nameMasterSteer(this->fileName("Master"));
0756   std::ofstream *mainSteerPtr = this->createSteerFile(nameMasterSteer, false);
0757   if (!mainSteerPtr)
0758     return "";
0759
0760   // add external steering files, if any
0761   std::vector<std::string> addfiles = myConfig.getParameter<std::vector<std::string> >("additionalSteerFiles");
0762   mySteeringFiles.insert(mySteeringFiles.end(), addfiles.begin(), addfiles.end());
0763
0764   // add steering files to master steering file
0765   std::ofstream &mainSteerRef = *mainSteerPtr;
0766   for (unsigned int iFile = 0; iFile < mySteeringFiles.size(); ++iFile) {
0767     mainSteerRef << mySteeringFiles[iFile] << "\n";
0768   }
0769
0770   // add binary files to master steering file
0771   mainSteerRef << "\nCfiles\n";
0772   for (unsigned int iFile = 0; iFile < binaryFiles.size(); ++iFile) {
0773     if (myRunDirectory.empty())
0774       mainSteerRef << binaryFiles[iFile] << "\n";
0775     else
0776       mainSteerRef << "../" + binaryFiles[iFile] << "\n";
0777   }
0778
0779   // add method
0780   mainSteerRef << "\nmethod  " << myConfig.getParameter<std::string>("method") << "\n";
0781
0782   // add further options
0783   const std::vector<std::string> opt(myConfig.getParameter<std::vector<std::string> >("options"));
0784   mainSteerRef << "\n* Outlier treatment and other options \n";
0785   for (unsigned int i = 0; i < opt.size(); ++i) {
0786     mainSteerRef << opt[i] << "\n";
0787   }
0788
0789   delete mainSteerPtr;  // close (and flush) again
0790
0791   return nameMasterSteer;
0792 }
0793
0794 //_________________________________________________________________________
0795 int PedeSteerer::runPede(const std::string &masterSteer) const {
0796   if (masterSteer.empty()) {
0797     edm::LogError("Alignment") << "@SUB=PedeSteerer::runPede"
0798                                << "Empty master steer file, stop";
0799     return 0;  //false;
0800   }
0801
0802   // Change pede command if running in different directory
0803   std::string command;
0804   if (myRunDirectory.empty()) {
0805     command = myConfig.getUntrackedParameter<std::string>("pedeCommand");
0806     (command += " ") += masterSteer;
0807   } else {
0808     command = "(cd " + myRunDirectory + " && ";
0809     command += myConfig.getUntrackedParameter<std::string>("pedeCommand");
0810     (command += " ") += "../" + masterSteer;
0811     command += ")";
0812   }
0813
0814   const std::string dump(myConfig.getUntrackedParameter<std::string>("pedeDump"));
0815   if (!dump.empty()) {
0816     command += " > ";
0817     (command += myDirectory) += dump;
0818   }
0819
0820   edm::LogInfo("Alignment") << "@SUB=PedeSteerer::runPede"
0821                             << "Start running " << command;
0822   // FIXME: Recommended interface to system commands?
0823   int shellReturn = gSystem->Exec(command.c_str());
0824   edm::LogInfo("Alignment") << "@SUB=PedeSteerer::runPede"
0825                             << "Command returns " << shellReturn;
0826
0827   return shellReturn;
0828 }