AlignmentRelative/src/AlignmentTask.cc

0001 /****************************************************************************
0002 * Authors:
0003 *  Jan Kašpar (jan.kaspar@gmail.com)
0004 ****************************************************************************/
0005
0006 #include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h"
0007 #include "CalibPPS/AlignmentRelative/interface/AlignmentConstraint.h"
0008
0009 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0010 #include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"
0011
0012 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0013
0014 #include <algorithm>
0015
0016 //----------------------------------------------------------------------------------------------------
0017
0018 using namespace std;
0019 using namespace edm;
0020
0021 //----------------------------------------------------------------------------------------------------
0022
0023 AlignmentTask::AlignmentTask(const ParameterSet &ps)
0024     : resolveShR(ps.getParameter<bool>("resolveShR")),
0025       resolveShZ(ps.getParameter<bool>("resolveShZ")),
0026       resolveRotZ(ps.getParameter<bool>("resolveRotZ")),
0027
0028       oneRotZPerPot(ps.getParameter<bool>("oneRotZPerPot")),
0029       useEqualMeanUMeanVRotZConstraints(ps.getParameter<bool>("useEqualMeanUMeanVRotZConstraints")),
0030
0031       fixedDetectorsConstraints(ps.getParameterSet("fixedDetectorsConstraints")),
0032       standardConstraints(ps.getParameterSet("standardConstraints")) {
0033   if (resolveShR) {
0034     quantityClasses.push_back(qcShR1);
0035     quantityClasses.push_back(qcShR2);
0036   }
0037
0038   if (resolveShZ) {
0039     quantityClasses.push_back(qcShZ);
0040   }
0041
0042   if (resolveRotZ) {
0043     quantityClasses.push_back(qcRotZ);
0044   }
0045 }
0046
0047 //----------------------------------------------------------------------------------------------------
0048
0049 void AlignmentTask::buildGeometry(const vector<unsigned int> &rpDecIds,
0050                                   const vector<unsigned int> &excludedSensors,
0051                                   const CTPPSGeometry *input,
0052                                   double z0,
0053                                   AlignmentGeometry &geometry) {
0054   geometry.z0 = z0;
0055
0056   // traverse full known geometry
0057   for (auto it = input->beginSensor(); it != input->endSensor(); ++it) {
0058     // skip excluded sensors
0059     if (find(excludedSensors.begin(), excludedSensors.end(), it->first) != excludedSensors.end())
0060       continue;
0061
0062     // is RP selected?
0063     const CTPPSDetId detId(it->first);
0064     const unsigned int rpDecId = 100 * detId.arm() + 10 * detId.station() + detId.rp();
0065     if (find(rpDecIds.begin(), rpDecIds.end(), rpDecId) == rpDecIds.end())
0066       continue;
0067
0068     // extract geometry data
0069     CTPPSGeometry::Vector c = input->localToGlobal(detId, CTPPSGeometry::Vector(0., 0., 0.));
0070     CTPPSGeometry::Vector d1 = input->localToGlobal(detId, CTPPSGeometry::Vector(1., 0., 0.)) - c;
0071     CTPPSGeometry::Vector d2 = input->localToGlobal(detId, CTPPSGeometry::Vector(0., 1., 0.)) - c;
0072
0073     // for strips: is it U plane?
0074     bool isU = false;
0075     if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0076       TotemRPDetId stripDetId(it->first);
0077       unsigned int rpNum = stripDetId.rp();
0078       unsigned int plNum = stripDetId.plane();
0079       isU = (plNum % 2 != 0);
0080       if (rpNum == 2 || rpNum == 3)
0081         isU = !isU;
0082     }
0083
0084     DetGeometry dg(c.z() - z0, c.x(), c.y(), isU);
0085     dg.setDirection(1, d1.x(), d1.y(), d1.z());
0086     dg.setDirection(2, d2.x(), d2.y(), d2.z());
0087     geometry.insert(it->first, dg);
0088   }
0089 }
0090
0091 //----------------------------------------------------------------------------------------------------
0092
0093 void AlignmentTask::buildIndexMaps() {
0094   // remove old mapping
0095   mapMeasurementIndeces.clear();
0096   mapQuantityIndeces.clear();
0097
0098   // loop over all classes
0099   for (const auto &qcl : quantityClasses) {
0100     // create entry for this class
0101     mapMeasurementIndeces[qcl];
0102
0103     // loop over all sensors
0104     unsigned int idxMeas = 0;
0105     unsigned int idxQuan = 0;
0106     for (const auto &git : geometry.getSensorMap()) {
0107       const unsigned int detId = git.first;
0108       const unsigned int subdetId = CTPPSDetId(git.first).subdetId();
0109
0110       // update measurement map
0111       if (qcl == qcShR1) {
0112         if (subdetId == CTPPSDetId::sdTimingDiamond)
0113           mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++;
0114         if (subdetId == CTPPSDetId::sdTrackingPixel)
0115           mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++;
0116       }
0117
0118       if (qcl == qcShR2) {
0119         if (subdetId == CTPPSDetId::sdTrackingStrip)
0120           mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++;
0121         if (subdetId == CTPPSDetId::sdTrackingPixel)
0122           mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++;
0123       }
0124
0125       if (qcl == qcShZ) {
0126         if (subdetId == CTPPSDetId::sdTrackingStrip)
0127           mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++;
0128         if (subdetId == CTPPSDetId::sdTrackingPixel)
0129           mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++;
0130         if (subdetId == CTPPSDetId::sdTrackingPixel)
0131           mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++;
0132         if (subdetId == CTPPSDetId::sdTimingDiamond)
0133           mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++;
0134       }
0135
0136       if (qcl == qcRotZ) {
0137         if (subdetId == CTPPSDetId::sdTrackingStrip)
0138           mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++;
0139         if (subdetId == CTPPSDetId::sdTrackingPixel)
0140           mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++;
0141         if (subdetId == CTPPSDetId::sdTrackingPixel)
0142           mapMeasurementIndeces[qcl][{detId, 2}] = idxMeas++;
0143         if (subdetId == CTPPSDetId::sdTimingDiamond)
0144           mapMeasurementIndeces[qcl][{detId, 1}] = idxMeas++;
0145       }
0146
0147       // update quantity map
0148       if (qcl == qcShR1) {
0149         if (subdetId == CTPPSDetId::sdTimingDiamond)
0150           mapQuantityIndeces[qcl][detId] = idxQuan++;
0151         if (subdetId == CTPPSDetId::sdTrackingPixel)
0152           mapQuantityIndeces[qcl][detId] = idxQuan++;
0153       }
0154
0155       if (qcl == qcShR2) {
0156         if (subdetId == CTPPSDetId::sdTrackingStrip)
0157           mapQuantityIndeces[qcl][detId] = idxQuan++;
0158         if (subdetId == CTPPSDetId::sdTrackingPixel)
0159           mapQuantityIndeces[qcl][detId] = idxQuan++;
0160       }
0161
0162       if (qcl == qcShZ) {
0163         if (subdetId == CTPPSDetId::sdTrackingStrip)
0164           mapQuantityIndeces[qcl][detId] = idxQuan++;
0165         if (subdetId == CTPPSDetId::sdTimingDiamond)
0166           mapQuantityIndeces[qcl][detId] = idxQuan++;
0167         if (subdetId == CTPPSDetId::sdTrackingPixel)
0168           mapQuantityIndeces[qcl][detId] = idxQuan++;
0169       }
0170
0171       if (qcl == qcRotZ) {
0172         if (subdetId == CTPPSDetId::sdTrackingStrip)
0173           mapQuantityIndeces[qcl][detId] = idxQuan++;
0174         if (subdetId == CTPPSDetId::sdTimingDiamond)
0175           mapQuantityIndeces[qcl][detId] = idxQuan++;
0176         if (subdetId == CTPPSDetId::sdTrackingPixel)
0177           mapQuantityIndeces[qcl][detId] = idxQuan++;
0178       }
0179     }
0180   }
0181 }
0182
0183 //----------------------------------------------------------------------------------------------------
0184
0185 signed int AlignmentTask::getMeasurementIndex(QuantityClass cl, unsigned int detId, unsigned int dirIdx) const {
0186   auto clit = mapMeasurementIndeces.find(cl);
0187   if (clit == mapMeasurementIndeces.end())
0188     return -1;
0189
0190   auto it = clit->second.find({detId, dirIdx});
0191   if (it == clit->second.end())
0192     return -1;
0193
0194   return it->second;
0195 }
0196
0197 //----------------------------------------------------------------------------------------------------
0198
0199 signed int AlignmentTask::getQuantityIndex(QuantityClass cl, unsigned int detId) const {
0200   auto clit = mapQuantityIndeces.find(cl);
0201   if (clit == mapQuantityIndeces.end())
0202     return -1;
0203
0204   auto it = clit->second.find(detId);
0205   if (it == clit->second.end())
0206     return -1;
0207
0208   return it->second;
0209 }
0210
0211 //----------------------------------------------------------------------------------------------------
0212
0213 string AlignmentTask::quantityClassTag(QuantityClass qc) const {
0214   switch (qc) {
0215     case qcShR1:
0216       return "ShR1";
0217     case qcShR2:
0218       return "ShR2";
0219     case qcShZ:
0220       return "ShZ";
0221     case qcRotZ:
0222       return "RotZ";
0223   }
0224
0225   throw cms::Exception("PPS") << "Unknown quantity class " << qc << ".";
0226 }
0227
0228 //----------------------------------------------------------------------------------------------------
0229
0230 unsigned int AlignmentTask::measurementsOfClass(QuantityClass qc) const {
0231   auto it = mapMeasurementIndeces.find(qc);
0232   if (it == mapMeasurementIndeces.end())
0233     return 0;
0234   else
0235     return it->second.size();
0236 }
0237
0238 //----------------------------------------------------------------------------------------------------
0239
0240 unsigned int AlignmentTask::quantitiesOfClass(QuantityClass qc) const {
0241   auto it = mapQuantityIndeces.find(qc);
0242   if (it == mapQuantityIndeces.end())
0243     return 0;
0244   else
0245     return it->second.size();
0246 }
0247
0248 //----------------------------------------------------------------------------------------------------
0249
0250 void AlignmentTask::buildFixedDetectorsConstraints(vector<AlignmentConstraint> &constraints) const {
0251   for (auto &quantityClass : quantityClasses) {
0252     // get input
0253     const string &tag = quantityClassTag(quantityClass);
0254
0255     const ParameterSet &classSettings = fixedDetectorsConstraints.getParameterSet(tag.c_str());
0256     vector<unsigned int> ids(classSettings.getParameter<vector<unsigned int>>("ids"));
0257     vector<double> values(classSettings.getParameter<vector<double>>("values"));
0258
0259     if (ids.size() != values.size())
0260       throw cms::Exception("PPS") << "Different number of constraint ids and values for " << tag << ".";
0261
0262     // determine number of constraints
0263     unsigned int size = ids.size();
0264
0265     // just one basic constraint
0266     if (oneRotZPerPot && quantityClass == qcRotZ) {
0267       if (size > 1)
0268         size = 1;
0269     }
0270
0271     // build constraints
0272     for (unsigned int j = 0; j < size; j++) {
0273       // prepare empty constraint
0274       AlignmentConstraint ac;
0275
0276       for (auto &qcit : quantityClasses) {
0277         ac.coef[qcit].ResizeTo(quantitiesOfClass(qcit));
0278         ac.coef[qcit].Zero();
0279       }
0280
0281       // set constraint name
0282       char buf[40];
0283       sprintf(buf, "%s: fixed plane %4u", tag.c_str(), ids[j]);
0284       ac.name = buf;
0285
0286       // get quantity index
0287       signed int qIndex = getQuantityIndex(quantityClass, ids[j]);
0288       if (qIndex < 0)
0289         throw cms::Exception("AlignmentTask::BuildFixedDetectorsConstraints")
0290             << "Quantity index for class " << quantityClass << " and id " << ids[j] << " is " << qIndex;
0291
0292       // set constraint coefficient and value
0293       ac.coef[quantityClass][qIndex] = 1.;
0294       ac.val = values[j] * 1E-3;
0295
0296       // save constraint
0297       constraints.push_back(ac);
0298     }
0299
0300     if (oneRotZPerPot && quantityClass == qcRotZ)
0301       buildOneRotZPerPotConstraints(constraints);
0302   }
0303 }
0304
0305 //----------------------------------------------------------------------------------------------------
0306
0307 void AlignmentTask::buildStandardConstraints(vector<AlignmentConstraint> &constraints) const {
0308   const vector<unsigned int> &decUnitIds = standardConstraints.getParameter<vector<unsigned int>>("units");
0309
0310   // count planes in RPs
0311   map<unsigned int, unsigned int> planesPerPot;
0312   for (const auto &it : geometry.getSensorMap()) {
0313     CTPPSDetId detId(it.first);
0314     planesPerPot[detId.rpId()]++;
0315   }
0316
0317   // ShR constraints
0318   if (resolveShR) {
0319     for (const auto &decUnitId : decUnitIds) {
0320       // prepare empty constraints
0321       AlignmentConstraint ac_X;
0322       for (auto &qcit : quantityClasses) {
0323         ac_X.coef[qcit].ResizeTo(quantitiesOfClass(qcit));
0324         ac_X.coef[qcit].Zero();
0325       }
0326       ac_X.val = 0;
0327
0328       AlignmentConstraint ac_Y(ac_X);
0329
0330       // set constraint names
0331       char buf[50];
0332       sprintf(buf, "ShR: unit %u, MeanX=0", decUnitId);
0333       ac_X.name = buf;
0334       sprintf(buf, "ShR: unit %u, MeanY=0", decUnitId);
0335       ac_Y.name = buf;
0336
0337       // traverse geometry
0338       for (const auto &git : geometry.getSensorMap()) {
0339         // stop is sensor not in the selected arm
0340         CTPPSDetId senId(git.first);
0341         unsigned int senDecUnit = senId.arm() * 100 + senId.station() * 10;
0342         if (senId.rp() > 2)
0343           senDecUnit += 1;
0344
0345         if (senDecUnit != decUnitId)
0346           continue;
0347
0348         // fill constraint for strip sensors
0349         if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0350           signed int qIndex = getQuantityIndex(qcShR2, git.first);
0351           if (qIndex < 0)
0352             throw cms::Exception("AlignmentTask::BuildStandardConstraints")
0353                 << "Cannot get quantity index for class " << qcShR2 << " and sensor id " << git.first << ".";
0354
0355           // determine weight
0356           const double weight = 1. / planesPerPot[senId.rpId()];
0357
0358           // set constraint coefficients
0359           ac_X.coef[qcShR2][qIndex] = git.second.getDirectionData(2).dx * weight;
0360           ac_Y.coef[qcShR2][qIndex] = git.second.getDirectionData(2).dy * weight;
0361         }
0362
0363         // fill constraint for pixel sensors
0364         if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0365           // get quantity indeces
0366           const signed int qIndex1 = getQuantityIndex(qcShR1, git.first);
0367           if (qIndex1 < 0)
0368             throw cms::Exception("AlignmentTask::BuildStandardConstraints")
0369                 << "Cannot get quantity index for class " << qcShR1 << " and sensor id " << git.first << ".";
0370
0371           const signed int qIndex2 = getQuantityIndex(qcShR2, git.first);
0372           if (qIndex2 < 0)
0373             throw cms::Exception("AlignmentTask::BuildStandardConstraints")
0374                 << "Cannot get quantity index for class " << qcShR2 << " and sensor id " << git.first << ".";
0375
0376           // determine weight (two constraints per plane)
0377           const double weight = 0.5 / planesPerPot[senId.rpId()];
0378
0379           // get geometry
0380           const double d1x = git.second.getDirectionData(1).dx;
0381           const double d1y = git.second.getDirectionData(1).dy;
0382           const double d2x = git.second.getDirectionData(2).dx;
0383           const double d2y = git.second.getDirectionData(2).dy;
0384
0385           // calculate coefficients, by inversion of this matrix relation
0386           //  [ s1 ] = [ d1x  d1y ] * [ de x ]
0387           //  [ s2 ]   [ d2x  d2y ]   [ de y ]
0388           const double D = d1x * d2y - d1y * d2x;
0389           const double coef_x_s1 = +d2y / D;
0390           const double coef_y_s1 = -d2x / D;
0391           const double coef_x_s2 = -d1y / D;
0392           const double coef_y_s2 = +d1x / D;
0393
0394           // set constraint coefficients
0395           ac_X.coef[qcShR1][qIndex1] = coef_x_s1 * weight;
0396           ac_Y.coef[qcShR1][qIndex1] = coef_y_s1 * weight;
0397           ac_X.coef[qcShR2][qIndex2] = coef_x_s2 * weight;
0398           ac_Y.coef[qcShR2][qIndex2] = coef_y_s2 * weight;
0399         }
0400       }
0401
0402       // add constraints
0403       constraints.push_back(ac_X);
0404       constraints.push_back(ac_Y);
0405     }
0406   }
0407
0408   // RotZ constraints
0409   if (resolveRotZ) {
0410     for (const auto &decUnitId : decUnitIds) {
0411       // prepare empty constraints
0412       AlignmentConstraint ac;
0413       for (unsigned int i = 0; i < quantityClasses.size(); i++) {
0414         ac.coef[quantityClasses[i]].ResizeTo(quantitiesOfClass(quantityClasses[i]));
0415         ac.coef[quantityClasses[i]].Zero();
0416       }
0417       ac.val = 0;
0418
0419       char buf[50];
0420       sprintf(buf, "RotZ: unit %u, Mean=0", decUnitId);
0421       ac.name = buf;
0422
0423       // traverse geometry
0424       for (const auto &git : geometry.getSensorMap()) {
0425         // stop is sensor not in the selected arm
0426         CTPPSDetId senId(git.first);
0427         unsigned int senDecUnit = senId.arm() * 100 + senId.station() * 10;
0428         if (senId.rp() > 2)
0429           senDecUnit += 1;
0430
0431         if (senDecUnit != decUnitId)
0432           continue;
0433
0434         // determine weight
0435         const double weight = 1. / planesPerPot[senId.rpId()];
0436
0437         // set coefficient
0438         signed int qIndex = getQuantityIndex(qcRotZ, git.first);
0439         ac.coef[qcRotZ][qIndex] = weight;
0440       }
0441
0442       constraints.push_back(ac);
0443     }
0444   }
0445
0446   if (resolveRotZ && oneRotZPerPot)
0447     buildOneRotZPerPotConstraints(constraints);
0448
0449   if (resolveRotZ && useEqualMeanUMeanVRotZConstraints)
0450     buildEqualMeanUMeanVRotZConstraints(constraints);
0451 }
0452
0453 //----------------------------------------------------------------------------------------------------
0454
0455 void AlignmentTask::buildOneRotZPerPotConstraints(std::vector<AlignmentConstraint> &constraints) const {
0456   // build map rp id --> sensor ids
0457   map<unsigned int, vector<unsigned int>> m;
0458   for (const auto &p : geometry.getSensorMap()) {
0459     CTPPSDetId detId(p.first);
0460     CTPPSDetId rpId = detId.rpId();
0461     unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0462     m[decRPId].push_back(p.first);
0463   }
0464
0465   // traverse all RPs
0466   for (const auto &p : m) {
0467     // build n_planes-1 constraints
0468     unsigned int prev_detId = 0;
0469     for (const auto &detId : p.second) {
0470       if (prev_detId != 0) {
0471         AlignmentConstraint ac;
0472
0473         char buf[100];
0474         sprintf(buf, "RotZ: RP %u, plane %u = plane %u", p.first, prev_detId, detId);
0475         ac.name = buf;
0476
0477         ac.val = 0;
0478
0479         for (auto &qcit : quantityClasses) {
0480           ac.coef[qcit].ResizeTo(quantitiesOfClass(qcit));
0481           ac.coef[qcit].Zero();
0482         }
0483
0484         signed int qIdx1 = getQuantityIndex(qcRotZ, prev_detId);
0485         signed int qIdx2 = getQuantityIndex(qcRotZ, detId);
0486
0487         ac.coef[qcRotZ][qIdx1] = +1.;
0488         ac.coef[qcRotZ][qIdx2] = -1.;
0489
0490         constraints.push_back(ac);
0491       }
0492
0493       prev_detId = detId;
0494     }
0495   }
0496 }
0497
0498 //----------------------------------------------------------------------------------------------------
0499
0500 void AlignmentTask::buildEqualMeanUMeanVRotZConstraints(vector<AlignmentConstraint> &constraints) const {
0501   // build map strip rp id --> pair( vector of U planes, vector of V planes )
0502   map<unsigned int, pair<vector<unsigned int>, vector<unsigned int>>> m;
0503   for (const auto &p : geometry.getSensorMap()) {
0504     CTPPSDetId detId(p.first);
0505
0506     if (detId.subdetId() != CTPPSDetId::sdTrackingStrip)
0507       continue;
0508
0509     CTPPSDetId rpId = detId.rpId();
0510     unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0511
0512     if (p.second.isU)
0513       m[decRPId].first.push_back(p.first);
0514     else
0515       m[decRPId].second.push_back(p.first);
0516   }
0517
0518   // loop over RPs
0519   for (const auto &p : m) {
0520     AlignmentConstraint ac;
0521
0522     char buf[100];
0523     sprintf(buf, "RotZ: RP %u, MeanU = MeanV", p.first);
0524     ac.name = buf;
0525
0526     ac.val = 0;
0527
0528     for (auto &qcit : quantityClasses) {
0529       ac.coef[qcit].ResizeTo(quantitiesOfClass(qcit));
0530       ac.coef[qcit].Zero();
0531     }
0532
0533     for (auto &&proj : {"U", "V"}) {
0534       const auto &planes = (proj == string("U")) ? p.second.first : p.second.second;
0535       const double c = ((proj == string("U")) ? -1. : +1.) / planes.size();
0536
0537       for (const auto &plane : planes) {
0538         signed int qIdx = getQuantityIndex(qcRotZ, plane);
0539         ac.coef[qcRotZ][qIdx] = c;
0540
0541         TotemRPDetId plId(plane);
0542       }
0543     }
0544
0545     constraints.push_back(ac);
0546   }
0547 }