CommonAlignmentAlgorithm/src/AlignmentParameterStore.cc

0001 /**
0002  * \file AlignmentParameterStore.cc
0003  *
0004  *  $Revision: 1.31 $
0005  *  $Date: 2011/05/23 20:50:32 $
0006  *  (last update by $Author: mussgill $)
0007  */
0008
0009 // This class's header should be first
0010 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "Alignment/CommonAlignment/interface/Alignable.h"
0017 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
0018 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
0019
0020 #include "Alignment/CommonAlignmentParametrization/interface/RigidBodyAlignmentParameters.h"
0021 #include "Alignment/CommonAlignmentParametrization/interface/ParametersToParametersDerivatives.h"
0022 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentExtendedCorrelationsStore.h"
0023 #include "DataFormats/GeometryCommonDetAlgo/interface/AlignmentPositionError.h"
0024
0025 #include "Geometry/CommonTopologies/interface/SurfaceDeformation.h"
0026 #include "Geometry/CommonTopologies/interface/SurfaceDeformationFactory.h"
0027
0028 //__________________________________________________________________________________________________
0029 AlignmentParameterStore::AlignmentParameterStore(const align::Alignables& alis, const edm::ParameterSet& config)
0030     : theAlignables(alis) {
0031   if (config.getUntrackedParameter<bool>("UseExtendedCorrelations")) {
0032     theCorrelationsStore =
0033         new AlignmentExtendedCorrelationsStore(config.getParameter<edm::ParameterSet>("ExtendedCorrelationsConfig"));
0034   } else {
0035     theCorrelationsStore = new AlignmentCorrelationsStore();
0036   }
0037
0038   edm::LogInfo("Alignment") << "@SUB=AlignmentParameterStore"
0039                             << "Created with " << theAlignables.size() << " alignables.";
0040
0041   // set hierarchy vs averaging constraints
0042   theTypeOfConstraints = NONE;
0043   const std::string cfgStrTypeOfConstraints(config.getParameter<std::string>("TypeOfConstraints"));
0044   if (cfgStrTypeOfConstraints == "hierarchy") {
0045     theTypeOfConstraints = HIERARCHY_CONSTRAINTS;
0046   } else if (cfgStrTypeOfConstraints == "approximate_averaging") {
0047     theTypeOfConstraints = APPROX_AVERAGING_CONSTRAINTS;
0048     edm::LogWarning("Alignment")
0049         << "@SUB=AlignmentParameterStore"
0050         << "\n\n\n******* WARNING ******************************************\n"
0051         << "Using approximate implementation of averaging constraints."
0052         << "This is not recommended."
0053         << "Consider to use 'hierarchy' constraints:"
0054         << "  AlignmentProducer.ParameterStore.TypeOfConstraints = cms.string('hierarchy')\n\n\n";
0055   } else {
0056     edm::LogError("BadArgument") << "@SUB=AlignmentParameterStore"
0057                                  << "Unknown type of hierarchy constraints '" << cfgStrTypeOfConstraints << "'";
0058   }
0059 }
0060
0061 //__________________________________________________________________________________________________
0062 AlignmentParameterStore::~AlignmentParameterStore() { delete theCorrelationsStore; }
0063
0064 //__________________________________________________________________________________________________
0065 CompositeAlignmentParameters AlignmentParameterStore::selectParameters(
0066     const std::vector<AlignableDet*>& alignabledets) const {
0067   std::vector<AlignableDetOrUnitPtr> detOrUnits;
0068   detOrUnits.reserve(alignabledets.size());
0069
0070   std::vector<AlignableDet*>::const_iterator it, iEnd;
0071   for (it = alignabledets.begin(), iEnd = alignabledets.end(); it != iEnd; ++it)
0072     detOrUnits.push_back(AlignableDetOrUnitPtr(*it));
0073
0074   return this->selectParameters(detOrUnits);
0075 }
0076
0077 //__________________________________________________________________________________________________
0078 CompositeAlignmentParameters AlignmentParameterStore::selectParameters(
0079     const std::vector<AlignableDetOrUnitPtr>& alignabledets) const {
0080   align::Alignables alignables;
0081   std::map<AlignableDetOrUnitPtr, Alignable*> alidettoalimap;
0082   std::map<Alignable*, int> aliposmap;
0083   std::map<Alignable*, int> alilenmap;
0084   int nparam = 0;
0085
0086   // iterate over AlignableDet's
0087   for (std::vector<AlignableDetOrUnitPtr>::const_iterator iad = alignabledets.begin(); iad != alignabledets.end();
0088        ++iad) {
0089     Alignable* ali = alignableFromAlignableDet(*iad);
0090     if (ali) {
0091       alidettoalimap[*iad] = ali;  // Add to map
0092       // Check if Alignable already there, insert into vector if not
0093       if (find(alignables.begin(), alignables.end(), ali) == alignables.end()) {
0094         alignables.push_back(ali);
0095         AlignmentParameters* ap = ali->alignmentParameters();
0096         nparam += ap->numSelected();
0097       }
0098     }
0099   }
0100
0101   AlgebraicVector* selpar = new AlgebraicVector(nparam, 0);
0102   AlgebraicSymMatrix* selcov = new AlgebraicSymMatrix(nparam, 0);
0103
0104   // Fill in parameters and corresponding covariance matricess
0105   int ipos = 1;  // NOTE: .sub indices start from 1
0106   align::Alignables::const_iterator it1;
0107   for (it1 = alignables.begin(); it1 != alignables.end(); ++it1) {
0108     AlignmentParameters* ap = (*it1)->alignmentParameters();
0109     selpar->sub(ipos, ap->selectedParameters());
0110     selcov->sub(ipos, ap->selectedCovariance());
0111     int npar = ap->numSelected();
0112     aliposmap[*it1] = ipos;
0113     alilenmap[*it1] = npar;
0114     ipos += npar;
0115   }
0116
0117   // Fill in the correlations. Has to be an extra loop, because the
0118   // AlignmentExtendedCorrelationsStore (if used) needs the
0119   // alignables' covariance matrices already present.
0120   ipos = 1;
0121   for (it1 = alignables.begin(); it1 != alignables.end(); ++it1) {
0122     int jpos = 1;
0123
0124     // Look for correlations between alignables
0125     align::Alignables::const_iterator it2;
0126     for (it2 = alignables.begin(); it2 != it1; ++it2) {
0127       theCorrelationsStore->correlations(*it1, *it2, *selcov, ipos - 1, jpos - 1);
0128       jpos += (*it2)->alignmentParameters()->numSelected();
0129     }
0130
0131     ipos += (*it1)->alignmentParameters()->numSelected();
0132   }
0133
0134   AlignmentParametersData::DataContainer data(new AlignmentParametersData(selpar, selcov));
0135   CompositeAlignmentParameters aap(data, alignables, alidettoalimap, aliposmap, alilenmap);
0136
0137   return aap;
0138 }
0139
0140 //__________________________________________________________________________________________________
0141 CompositeAlignmentParameters AlignmentParameterStore::selectParameters(const align::Alignables& alignables) const {
0142   align::Alignables selectedAlignables;
0143   std::map<AlignableDetOrUnitPtr, Alignable*> alidettoalimap;  // This map won't be filled!!!
0144   std::map<Alignable*, int> aliposmap;
0145   std::map<Alignable*, int> alilenmap;
0146   int nparam = 0;
0147
0148   // iterate over Alignable's
0149   align::Alignables::const_iterator ita;
0150   for (ita = alignables.begin(); ita != alignables.end(); ++ita) {
0151     // Check if Alignable already there, insert into vector if not
0152     if (find(selectedAlignables.begin(), selectedAlignables.end(), *ita) == selectedAlignables.end()) {
0153       selectedAlignables.push_back(*ita);
0154       AlignmentParameters* ap = (*ita)->alignmentParameters();
0155       nparam += ap->numSelected();
0156     }
0157   }
0158
0159   AlgebraicVector* selpar = new AlgebraicVector(nparam, 0);
0160   AlgebraicSymMatrix* selcov = new AlgebraicSymMatrix(nparam, 0);
0161
0162   // Fill in parameters and corresponding covariance matrices
0163   int ipos = 1;  // NOTE: .sub indices start from 1
0164   align::Alignables::const_iterator it1;
0165   for (it1 = selectedAlignables.begin(); it1 != selectedAlignables.end(); ++it1) {
0166     AlignmentParameters* ap = (*it1)->alignmentParameters();
0167     selpar->sub(ipos, ap->selectedParameters());
0168     selcov->sub(ipos, ap->selectedCovariance());
0169     int npar = ap->numSelected();
0170     aliposmap[*it1] = ipos;
0171     alilenmap[*it1] = npar;
0172     ipos += npar;
0173   }
0174
0175   // Fill in the correlations. Has to be an extra loop, because the
0176   // AlignmentExtendedCorrelationsStore (if used) needs the
0177   // alignables' covariance matrices already present.
0178   ipos = 1;
0179   for (it1 = selectedAlignables.begin(); it1 != selectedAlignables.end(); ++it1) {
0180     int jpos = 1;
0181
0182     // Look for correlations between alignables
0183     align::Alignables::const_iterator it2;
0184     for (it2 = selectedAlignables.begin(); it2 != it1; ++it2) {
0185       theCorrelationsStore->correlations(*it1, *it2, *selcov, ipos - 1, jpos - 1);
0186       jpos += (*it2)->alignmentParameters()->numSelected();
0187     }
0188
0189     ipos += (*it1)->alignmentParameters()->numSelected();
0190   }
0191
0192   AlignmentParametersData::DataContainer data(new AlignmentParametersData(selpar, selcov));
0193   CompositeAlignmentParameters aap(data, selectedAlignables, alidettoalimap, aliposmap, alilenmap);
0194
0195   return aap;
0196 }
0197
0198 //__________________________________________________________________________________________________
0199 void AlignmentParameterStore::updateParameters(const CompositeAlignmentParameters& aap, bool updateCorrelations) {
0200   align::Alignables alignables = aap.components();
0201   const AlgebraicVector& parameters = aap.parameters();
0202   const AlgebraicSymMatrix& covariance = aap.covariance();
0203
0204   int ipos = 1;  // NOTE: .sub indices start from 1
0205
0206   // Loop over alignables
0207   for (align::Alignables::const_iterator it = alignables.begin(); it != alignables.end(); ++it) {
0208     // Update parameters and local covariance
0209     AlignmentParameters* ap = (*it)->alignmentParameters();
0210     int nsel = ap->numSelected();
0211     AlgebraicVector subvec = parameters.sub(ipos, ipos + nsel - 1);
0212     AlgebraicSymMatrix subcov = covariance.sub(ipos, ipos + nsel - 1);
0213     AlignmentParameters* apnew = ap->cloneFromSelected(subvec, subcov);
0214     (*it)->setAlignmentParameters(apnew);
0215
0216     // Now update correlations between detectors
0217     if (updateCorrelations) {
0218       int jpos = 1;
0219       for (align::Alignables::const_iterator it2 = alignables.begin(); it2 != it; ++it2) {
0220         theCorrelationsStore->setCorrelations(*it, *it2, covariance, ipos - 1, jpos - 1);
0221         jpos += (*it2)->alignmentParameters()->numSelected();
0222       }
0223     }
0224
0225     ipos += nsel;
0226   }
0227 }
0228
0229 //__________________________________________________________________________________________________
0230 align::Alignables AlignmentParameterStore::validAlignables(void) const {
0231   align::Alignables result;
0232   for (align::Alignables::const_iterator iali = theAlignables.begin(); iali != theAlignables.end(); ++iali)
0233     if ((*iali)->alignmentParameters()->isValid())
0234       result.push_back(*iali);
0235
0236   LogDebug("Alignment") << "@SUB=AlignmentParameterStore::validAlignables"
0237                         << "Valid alignables: " << result.size() << "out of " << theAlignables.size();
0238   return result;
0239 }
0240
0241 //__________________________________________________________________________________________________
0242 Alignable* AlignmentParameterStore::alignableFromAlignableDet(const AlignableDetOrUnitPtr& _alignableDet) const {
0243   AlignableDetOrUnitPtr alignableDet = _alignableDet;
0244   Alignable* mother = alignableDet;
0245   while (mother) {
0246     if (mother->alignmentParameters())
0247       return mother;
0248     mother = mother->mother();
0249   }
0250
0251   return nullptr;
0252 }
0253
0254 //__________________________________________________________________________________________________
0255 void AlignmentParameterStore::applyParameters(void) {
0256   align::Alignables::const_iterator iali;
0257   for (iali = theAlignables.begin(); iali != theAlignables.end(); ++iali)
0258     applyParameters(*iali);
0259 }
0260
0261 //__________________________________________________________________________________________________
0262 void AlignmentParameterStore::applyParameters(Alignable* alignable) {
0263   AlignmentParameters* pars = (alignable ? alignable->alignmentParameters() : nullptr);
0264   if (!pars) {
0265     throw cms::Exception("BadAlignable") << "applyParameters: provided alignable does not have alignment parameters";
0266   }
0267   pars->apply();
0268 }
0269
0270 //__________________________________________________________________________________________________
0271 void AlignmentParameterStore::resetParameters(void) {
0272   // Erase contents of correlation map
0273   theCorrelationsStore->resetCorrelations();
0274
0275   // Iterate over alignables in the store and reset parameters
0276   align::Alignables::const_iterator iali;
0277   for (iali = theAlignables.begin(); iali != theAlignables.end(); ++iali)
0278     resetParameters(*iali);
0279 }
0280
0281 //__________________________________________________________________________________________________
0282 void AlignmentParameterStore::resetParameters(Alignable* ali) {
0283   if (ali) {
0284     // Get alignment parameters for this alignable
0285     AlignmentParameters* ap = ali->alignmentParameters();
0286     if (ap) {
0287       int npar = ap->numSelected();
0288
0289       AlgebraicVector par(npar, 0);
0290       AlgebraicSymMatrix cov(npar, 0);
0291       AlignmentParameters* apnew = ap->cloneFromSelected(par, cov);
0292       ali->setAlignmentParameters(apnew);
0293       apnew->setValid(false);
0294     } else
0295       edm::LogError("BadArgument") << "@SUB=AlignmentParameterStore::resetParameters"
0296                                    << "alignable has no alignment parameter";
0297   } else
0298     edm::LogError("BadArgument") << "@SUB=AlignmentParameterStore::resetParameters"
0299                                  << "argument is NULL";
0300 }
0301
0302 //__________________________________________________________________________________________________
0303 void AlignmentParameterStore::cacheTransformations(void) {
0304   align::Alignables::const_iterator iali;
0305   for (iali = theAlignables.begin(); iali != theAlignables.end(); ++iali)
0306     (*iali)->cacheTransformation();
0307 }
0308
0309 //__________________________________________________________________________________________________
0310 void AlignmentParameterStore::cacheTransformations(const align::RunNumber& run) {
0311   for (const auto& iali : theAlignables)
0312     iali->cacheTransformation(run);
0313 }
0314
0315 //__________________________________________________________________________________________________
0316 void AlignmentParameterStore::restoreCachedTransformations(void) {
0317   align::Alignables::const_iterator iali;
0318   for (iali = theAlignables.begin(); iali != theAlignables.end(); ++iali)
0319     (*iali)->restoreCachedTransformation();
0320 }
0321
0322 //__________________________________________________________________________________________________
0323 void AlignmentParameterStore::restoreCachedTransformations(const align::RunNumber& run) {
0324   for (const auto& iali : theAlignables)
0325     iali->restoreCachedTransformation(run);
0326 }
0327
0328 //__________________________________________________________________________________________________
0329 void AlignmentParameterStore::acquireRelativeParameters(void) {
0330   unsigned int nAlignables = theAlignables.size();
0331
0332   for (unsigned int i = 0; i < nAlignables; ++i) {
0333     Alignable* ali = theAlignables[i];
0334
0335     RigidBodyAlignmentParameters* ap = dynamic_cast<RigidBodyAlignmentParameters*>(ali->alignmentParameters());
0336
0337     if (!ap)
0338       throw cms::Exception("BadAlignable") << "acquireRelativeParameters: "
0339                                            << "provided alignable does not have rigid body alignment parameters";
0340
0341     AlgebraicVector par(ap->size(), 0);
0342     AlgebraicSymMatrix cov(ap->size(), 0);
0343
0344     // Get displacement and transform global->local
0345     align::LocalVector dloc = ali->surface().toLocal(ali->displacement());
0346     par[0] = dloc.x();
0347     par[1] = dloc.y();
0348     par[2] = dloc.z();
0349
0350     // Transform to local euler angles
0351     align::EulerAngles euloc = align::toAngles(ali->surface().toLocal(ali->rotation()));
0352     par[3] = euloc(1);
0353     par[4] = euloc(2);
0354     par[5] = euloc(3);
0355
0356     // Clone parameters
0357     RigidBodyAlignmentParameters* apnew = ap->clone(par, cov);
0358
0359     ali->setAlignmentParameters(apnew);
0360   }
0361 }
0362
0363 //__________________________________________________________________________________________________
0364 // Get type/layer from Alignable
0365 // type: -6   -5   -4   -3   -2    -1     1     2    3    4    5    6
0366 //      TEC- TOB- TID- TIB- PxEC- PxBR- PxBr+ PxEC+ TIB+ TID+ TOB+ TEC+
0367 // Layers start from zero
0368 std::pair<int, int> AlignmentParameterStore::typeAndLayer(const Alignable* ali, const TrackerTopology* tTopo) const {
0369   return TrackerAlignableId().typeAndLayerFromDetId(ali->id(), tTopo);
0370 }
0371
0372 //__________________________________________________________________________________________________
0373 void AlignmentParameterStore::applyAlignableAbsolutePositions(const align::Alignables& alivec,
0374                                                               const AlignablePositions& newpos,
0375                                                               int& ierr) {
0376   unsigned int nappl = 0;
0377   ierr = 0;
0378
0379   // Iterate over list of alignables
0380   for (align::Alignables::const_iterator iali = alivec.begin(); iali != alivec.end(); ++iali) {
0381     Alignable* ali = *iali;
0382     align::ID id = ali->id();
0383     align::StructureType typeId = ali->alignableObjectId();
0384
0385     // Find corresponding entry in AlignablePositions
0386     bool found = false;
0387     for (AlignablePositions::const_iterator ipos = newpos.begin(); ipos != newpos.end(); ++ipos) {
0388       if (id == ipos->id() && typeId == ipos->objId()) {
0389         if (found) {
0390           edm::LogError("DuplicatePosition") << "New positions for alignable found more than once!";
0391         } else {
0392           // New position/rotation
0393           const align::PositionType& pnew = ipos->pos();
0394           const align::RotationType& rnew = ipos->rot();
0395           const std::vector<double>& dnew = ipos->deformationParameters();
0396           // Current position / rotation
0397           const align::PositionType& pold = ali->globalPosition();
0398           const align::RotationType& rold = ali->globalRotation();
0399           // Current surf. deformation
0400           std::vector<std::pair<int, SurfaceDeformation*> > dold_id_pairs;
0401           SurfaceDeformation* dold_obj = nullptr;
0402           SurfaceDeformationFactory::Type dtype = SurfaceDeformationFactory::kNoDeformations;
0403           std::vector<double> dold;
0404           if (1 == ali->surfaceDeformationIdPairs(dold_id_pairs)) {  // might not have any...
0405             dold_obj = dold_id_pairs[0].second;
0406             dold = dold_obj->parameters();
0407             dtype = (SurfaceDeformationFactory::Type)dold_obj->type();
0408           }
0409
0410           // shift needed to move from current to new position
0411           align::GlobalVector posDiff = pnew - pold;
0412           align::RotationType rotDiff = rold.multiplyInverse(rnew);
0413           align::rectify(rotDiff);  // correct for rounding errors
0414           ali->move(posDiff);
0415           ali->rotateInGlobalFrame(rotDiff);
0416           LogDebug("NewPosition") << "moving by:" << posDiff;
0417           LogDebug("NewRotation") << "rotating by:\n" << rotDiff;
0418
0419           // add the surface deformations
0420           // If an old surface deformation record exists, ensure that the added deformation has the same type and size.
0421           if (!dold.empty() && dtype != SurfaceDeformationFactory::kNoDeformations && dnew.size() == dold.size()) {
0422             std::vector<double> defDiff;
0423             defDiff.reserve(dold.size());
0424             for (unsigned int i = 0; i < dold.size(); i++)
0425               defDiff.push_back(dnew[i] - dold[i]);
0426             auto deform = SurfaceDeformationFactory::create(dtype, defDiff);
0427             edm::LogInfo("Alignment") << "@SUB=AlignmentParameterStore::applyAlignableAbsolutePositions"
0428                                       << "Adding surface deformation of type "
0429                                       << SurfaceDeformationFactory::surfaceDeformationTypeName(
0430                                              (SurfaceDeformationFactory::Type)deform->type())
0431                                       << ", size " << defDiff.size() << " and first element " << defDiff.at(0)
0432                                       << " to alignable with id / type: " << id << " / " << typeId;
0433             ali->addSurfaceDeformation(deform, true);
0434             delete deform;
0435           }
0436           // In case no old surface deformation record exists, only ensure that the new surface deformation record has size>0. Size check is done elsewhere.
0437           else if (!dnew.empty()) {
0438             auto deform = SurfaceDeformationFactory::create(dnew);
0439             edm::LogInfo("Alignment") << "@SUB=AlignmentParameterStore::applyAlignableAbsolutePositions"
0440                                       << "Setting surface deformation of type "
0441                                       << SurfaceDeformationFactory::surfaceDeformationTypeName(
0442                                              (SurfaceDeformationFactory::Type)deform->type())
0443                                       << ", size " << dnew.size() << " and first element " << dnew.at(0)
0444                                       << " to alignable with id / type: " << id << " / " << typeId;
0445             ali->addSurfaceDeformation(deform, true);  // Equivalent to setSurfaceDeformation in this case
0446             delete deform;
0447           }
0448           // If there is no new surface deformation record, do nothing.
0449
0450           // add position error
0451           // AlignmentPositionError ape(shift.x(),shift.y(),shift.z());
0452           // (*iali)->addAlignmentPositionError(ape);
0453           // (*iali)->addAlignmentPositionErrorFromRotation(rot);
0454
0455           found = true;
0456           ++nappl;
0457         }
0458       }
0459     }
0460   }
0461
0462   if (nappl < newpos.size())
0463     edm::LogError("Mismatch") << "Applied only " << nappl << " new positions"
0464                               << " out of " << newpos.size();
0465
0466   LogDebug("NewPositions") << "Applied new positions for " << nappl << " out of " << alivec.size() << " alignables.";
0467 }
0468
0469 //__________________________________________________________________________________________________
0470 void AlignmentParameterStore::applyAlignableRelativePositions(const align::Alignables& alivec,
0471                                                               const AlignableShifts& shifts,
0472                                                               int& ierr) {
0473   ierr = 0;
0474   unsigned int nappl = 0;
0475   unsigned int nAlignables = alivec.size();
0476
0477   for (unsigned int i = 0; i < nAlignables; ++i) {
0478     Alignable* ali = alivec[i];
0479
0480     align::ID id = ali->id();
0481     align::StructureType typeId = ali->alignableObjectId();
0482
0483     // Find corresponding entry in AlignableShifts
0484     bool found = false;
0485     for (AlignableShifts::const_iterator ipos = shifts.begin(); ipos != shifts.end(); ++ipos) {
0486       if (id == ipos->id() && typeId == ipos->objId()) {
0487         if (found) {
0488           edm::LogError("DuplicatePosition") << "New positions for alignable found more than once!";
0489         } else {
0490           // Current surf. deformation
0491           std::vector<std::pair<int, SurfaceDeformation*> > dold_id_pairs;
0492           SurfaceDeformation* dold_obj = nullptr;
0493           SurfaceDeformationFactory::Type dtype = SurfaceDeformationFactory::kNoDeformations;
0494           std::vector<double> dold;
0495           if (1 == ali->surfaceDeformationIdPairs(dold_id_pairs)) {  // might not have any...
0496             dold_obj = dold_id_pairs[0].second;
0497             dold = dold_obj->parameters();
0498             dtype = (SurfaceDeformationFactory::Type)dold_obj->type();
0499           }
0500
0501           ali->move(ipos->pos());
0502           ali->rotateInGlobalFrame(ipos->rot());
0503
0504           const std::vector<double>& defDiff = ipos->deformationParameters();
0505           // If an old surface deformation record exists, ensure that the added deformation has the same type and size.
0506           if (!dold.empty() && dtype != SurfaceDeformationFactory::kNoDeformations && defDiff.size() == dold.size()) {
0507             auto deform = SurfaceDeformationFactory::create(dtype, defDiff);
0508             edm::LogInfo("Alignment") << "@SUB=AlignmentParameterStore::applyAlignableRelativePositions"
0509                                       << "Adding surface deformation of type "
0510                                       << SurfaceDeformationFactory::surfaceDeformationTypeName(
0511                                              (SurfaceDeformationFactory::Type)deform->type())
0512                                       << ", size " << defDiff.size() << " and first element " << defDiff.at(0)
0513                                       << " to alignable with id / type: " << id << " / " << typeId;
0514             ali->addSurfaceDeformation(deform, true);
0515             delete deform;
0516           }
0517           // In case no old surface deformation record exists, only ensure that the new surface deformation record has size>0. Size check is done elsewhere.
0518           else if (!defDiff.empty()) {
0519             auto deform = SurfaceDeformationFactory::create(defDiff);
0520             edm::LogInfo("Alignment") << "@SUB=AlignmentParameterStore::applyAlignableRelativePositions"
0521                                       << "Setting surface deformation of type "
0522                                       << SurfaceDeformationFactory::surfaceDeformationTypeName(
0523                                              (SurfaceDeformationFactory::Type)deform->type())
0524                                       << ", size " << defDiff.size() << " and first element " << defDiff.at(0)
0525                                       << " to alignable with id / type: " << id << " / " << typeId;
0526             ali->addSurfaceDeformation(deform, true);  // Equivalent to setSurfaceDeformation in this case
0527             delete deform;
0528           }
0529           // If there is no new surface deformation record, do nothing.
0530
0531           // Add position error
0532           //AlignmentPositionError ape(pnew.x(),pnew.y(),pnew.z());
0533           //ali->addAlignmentPositionError(ape);
0534           //ali->addAlignmentPositionErrorFromRotation(rnew);
0535
0536           found = true;
0537           ++nappl;
0538         }
0539       }
0540     }
0541   }
0542
0543   if (nappl < shifts.size())
0544     edm::LogError("Mismatch") << "Applied only " << nappl << " new positions"
0545                               << " out of " << shifts.size();
0546
0547   LogDebug("NewPositions") << "Applied new positions for " << nappl << " alignables.";
0548 }
0549
0550 //__________________________________________________________________________________________________
0551 void AlignmentParameterStore::attachAlignmentParameters(const Parameters& parvec, int& ierr) {
0552   attachAlignmentParameters(theAlignables, parvec, ierr);
0553 }
0554
0555 //__________________________________________________________________________________________________
0556 void AlignmentParameterStore::attachAlignmentParameters(const align::Alignables& alivec,
0557                                                         const Parameters& parvec,
0558                                                         int& ierr) {
0559   int ipass = 0;
0560   int ifail = 0;
0561   ierr = 0;
0562
0563   // Iterate over alignables
0564   for (align::Alignables::const_iterator iali = alivec.begin(); iali != alivec.end(); ++iali) {
0565     // Iterate over Parameters
0566     bool found = false;
0567     for (Parameters::const_iterator ipar = parvec.begin(); ipar != parvec.end(); ++ipar) {
0568       // Get new alignment parameters
0569       AlignmentParameters* ap = *ipar;
0570
0571       // Check if parameters belong to alignable
0572       if (ap->alignable() == (*iali)) {
0573         if (!found) {
0574           (*iali)->setAlignmentParameters(ap);
0575           ++ipass;
0576           found = true;
0577         } else
0578           edm::LogError("Alignment") << "@SUB=AlignmentParameterStore::attachAlignmentParameters"
0579                                      << "More than one parameters for Alignable.";
0580       }
0581     }
0582     if (!found)
0583       ++ifail;
0584   }
0585   if (ifail > 0)
0586     ierr = -1;
0587
0588   LogDebug("attachAlignmentParameters") << " Parameters, Alignables: " << parvec.size() << "," << alivec.size()
0589                                         << "\n pass,fail: " << ipass << "," << ifail;
0590 }
0591
0592 //__________________________________________________________________________________________________
0593 void AlignmentParameterStore::attachCorrelations(const Correlations& cormap, bool overwrite, int& ierr) {
0594   attachCorrelations(theAlignables, cormap, overwrite, ierr);
0595 }
0596
0597 //__________________________________________________________________________________________________
0598 void AlignmentParameterStore::attachCorrelations(const align::Alignables& alivec,
0599                                                  const Correlations& cormap,
0600                                                  bool overwrite,
0601                                                  int& ierr) {
0602   ierr = 0;
0603   int icount = 0;
0604
0605   // Iterate over correlations
0606   for (Correlations::const_iterator icor = cormap.begin(); icor != cormap.end(); ++icor) {
0607     AlgebraicMatrix mat = (*icor).second;
0608     Alignable* ali1 = (*icor).first.first;
0609     Alignable* ali2 = (*icor).first.second;
0610
0611     // Check if alignables exist
0612     if (find(alivec.begin(), alivec.end(), ali1) != alivec.end() &&
0613         find(alivec.begin(), alivec.end(), ali2) != alivec.end()) {
0614       // Check if correlations already existing between these alignables
0615       if (!theCorrelationsStore->correlationsAvailable(ali1, ali2) || (overwrite)) {
0616         theCorrelationsStore->setCorrelations(ali1, ali2, mat);
0617         ++icount;
0618       } else
0619         edm::LogInfo("AlreadyExists") << "Correlation existing and not overwritten";
0620     } else
0621       edm::LogInfo("IgnoreCorrelation") << "Ignoring correlation with no alignables!";
0622   }
0623
0624   LogDebug("attachCorrelations") << " Alignables,Correlations: " << alivec.size() << "," << cormap.size()
0625                                  << "\n applied: " << icount;
0626 }
0627
0628 //__________________________________________________________________________________________________
0629 void AlignmentParameterStore::attachUserVariables(const align::Alignables& alivec,
0630                                                   const std::vector<AlignmentUserVariables*>& uvarvec,
0631                                                   int& ierr) {
0632   ierr = 0;
0633
0634   LogDebug("DumpArguments") << "size of alivec:   " << alivec.size() << "\nsize of uvarvec: " << uvarvec.size();
0635
0636   std::vector<AlignmentUserVariables*>::const_iterator iuvar = uvarvec.begin();
0637
0638   for (align::Alignables::const_iterator iali = alivec.begin(); iali != alivec.end(); ++iali, ++iuvar) {
0639     AlignmentParameters* ap = (*iali)->alignmentParameters();
0640     AlignmentUserVariables* uvarnew = (*iuvar);
0641     ap->setUserVariables(uvarnew);
0642   }
0643 }
0644
0645 //__________________________________________________________________________________________________
0646 void AlignmentParameterStore::setAlignmentPositionError(const align::Alignables& alivec,
0647                                                         double valshift,
0648                                                         double valrot) {
0649   unsigned int nAlignables = alivec.size();
0650
0651   for (unsigned int i = 0; i < nAlignables; ++i) {
0652     Alignable* ali = alivec[i];
0653
0654     // First reset APE
0655     AlignmentPositionError nulApe(0, 0, 0);
0656     ali->setAlignmentPositionError(nulApe, true);
0657
0658     // Set APE from displacement
0659     AlignmentPositionError ape(valshift, valshift, valshift);
0660     if (valshift > 0.)
0661       ali->addAlignmentPositionError(ape, true);
0662     else
0663       ali->setAlignmentPositionError(ape, true);
0664     // GF: Resetting and setting as above does not really make sense to me,
0665     //     and adding to zero or setting is the same! I'd just do
0666     //ali->setAlignmentPositionError(AlignmentPositionError ape(valshift,valshift,valshift),true);
0667
0668     // Set APE from rotation
0669     align::EulerAngles r(3);
0670     r(1) = valrot;
0671     r(2) = valrot;
0672     r(3) = valrot;
0673     ali->addAlignmentPositionErrorFromRotation(align::toMatrix(r), true);
0674   }
0675
0676   LogDebug("StoreAPE") << "Store APE from shift: " << valshift << "\nStore APE from rotation: " << valrot;
0677 }
0678
0679 //__________________________________________________________________________________________________
0680 bool AlignmentParameterStore ::hierarchyConstraints(const Alignable* ali,
0681                                                     const align::Alignables& aliComps,
0682                                                     std::vector<std::vector<ParameterId> >& paramIdsVecOut,
0683                                                     std::vector<std::vector<double> >& factorsVecOut,
0684                                                     bool all,
0685                                                     double epsilon) const {
0686   // Weak point if all = false:
0687   // Ignores constraints between non-subsequent levels in case the parameter is not considered in
0688   // the intermediate level, e.g. global z for dets and layers is aligned, but not for rods!
0689   if (!ali || !ali->alignmentParameters())
0690     return false;
0691
0692   const std::vector<bool>& aliSel = ali->alignmentParameters()->selector();
0693   paramIdsVecOut.clear();
0694   factorsVecOut.clear();
0695
0696   bool firstComp = true;
0697   for (align::Alignables::const_iterator iComp = aliComps.begin(), iCompE = aliComps.end(); iComp != iCompE; ++iComp) {
0698     const ParametersToParametersDerivatives p2pDerivs(**iComp, *ali);
0699     if (!p2pDerivs.isOK()) {
0700       // std::cerr << (*iComp)->alignmentParameters()->type() << " "
0701       //        << ali->alignmentParameters()->type() << std::endl;
0702       throw cms::Exception("BadConfig") << "AlignmentParameterStore::hierarchyConstraints"
0703                                         << " Bad match of types of AlignmentParameters classes.\n";
0704       return false;
0705     }
0706     const std::vector<bool>& aliCompSel = (*iComp)->alignmentParameters()->selector();
0707     for (unsigned int iParMast = 0, iParMastUsed = 0; iParMast < aliSel.size(); ++iParMast) {
0708       if (!all && !aliSel[iParMast])
0709         continue;       // no higher level parameter & constraint deselected
0710       if (firstComp) {  // fill output with empty arrays
0711         paramIdsVecOut.push_back(std::vector<ParameterId>());
0712         factorsVecOut.push_back(std::vector<double>());
0713       }
0714       for (unsigned int iParComp = 0; iParComp < aliCompSel.size(); ++iParComp) {
0715         if (aliCompSel[iParComp]) {
0716           double factor = 0.;
0717           if (theTypeOfConstraints == HIERARCHY_CONSTRAINTS) {
0718             // hierachy constraints
0719             factor = p2pDerivs(iParMast, iParComp);
0720           } else if (theTypeOfConstraints == APPROX_AVERAGING_CONSTRAINTS) {
0721             // CHK poor mans averaging constraints
0722             factor = p2pDerivs(iParMast, iParComp);
0723             if (iParMast < 3 && (iParComp % 9) >= 3)
0724               factor = 0.;
0725           }
0726           if (fabs(factor) > epsilon) {
0727             paramIdsVecOut[iParMastUsed].push_back(ParameterId(*iComp, iParComp));
0728             factorsVecOut[iParMastUsed].push_back(factor);
0729           }
0730         }
0731       }
0732       ++iParMastUsed;
0733     }
0734     firstComp = false;
0735   }  // end loop on components
0736
0737   return true;
0738 }