MillePedeAlignmentAlgorithm/plugins/MillePedeAlignmentAlgorithm.cc

0001 /**
0002  * \file MillePedeAlignmentAlgorithm.cc
0003  *
0004  *  \author    : Gero Flucke
0005  *  date       : October 2006
0006  *  $Revision: 1.80 $
0007  *  $Date: 2013/01/07 20:21:32 $
0008  *  (last update by $Author: wmtan $)
0009  */
0010
0011 #include "MillePedeAlignmentAlgorithm.h"
0012
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/Run.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0019
0020 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeFileReader.h"
0021 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeMonitor.h"
0022 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeVariables.h"
0023 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeVariablesIORoot.h"
0024 #include "Alignment/MillePedeAlignmentAlgorithm/src/Mille.h"        // 'unpublished' interface located in src
0025 #include "Alignment/MillePedeAlignmentAlgorithm/src/PedeSteerer.h"  // ditto
0026 #include "Alignment/MillePedeAlignmentAlgorithm/src/PedeReader.h"   // ditto
0027 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
0028 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerPluginFactory.h"
0029
0030 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryBase.h"
0031 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryPlugin.h"
0032
0033 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0034 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentIORoot.h"
0035 #include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
0036
0037 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
0038 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
0039 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
0040
0041 // includes to make known that they inherit from Alignable:
0042 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0043 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0044 #include "Alignment/CommonAlignment/interface/AlignableExtras.h"
0045
0046 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0047 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0048 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0049 #include "CondFormats/AlignmentRecord/interface/TrackerSurfaceDeformationRcd.h"
0050
0051 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0052 #include "DataFormats/TrackReco/interface/Track.h"
0053 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0054 #include "DataFormats/Alignment/interface/TkFittedLasBeam.h"
0055 #include "Alignment/LaserAlignment/interface/TsosVectorCollection.h"
0056
0057 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0058 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0059 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0060
0061 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0062
0063 #include <fstream>
0064 #include <sstream>
0065 #include <algorithm>
0066 #include <sys/stat.h>
0067
0068 #include <TMath.h>
0069 typedef TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer;
0070 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
0071 typedef TrajectoryFactoryBase::ReferenceTrajectoryCollection RefTrajColl;
0072
0073 // Includes for PXB survey
0074 #include "Alignment/SurveyAnalysis/interface/SurveyPxbImage.h"
0075 #include "Alignment/SurveyAnalysis/interface/SurveyPxbImageLocalFit.h"
0076 #include "Alignment/SurveyAnalysis/interface/SurveyPxbImageReader.h"
0077 #include "Alignment/SurveyAnalysis/interface/SurveyPxbDicer.h"
0078 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0079 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0080
0081 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
0082
0083 using namespace gbl;
0084
0085 // Constructor ----------------------------------------------------------------
0086 //____________________________________________________
0087 MillePedeAlignmentAlgorithm::MillePedeAlignmentAlgorithm(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC)
0088     : AlignmentAlgorithmBase(cfg, iC),
0089       topoToken_(iC.esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
0090       aliThrToken_(iC.esConsumes<AlignPCLThresholdsHG, AlignPCLThresholdsHGRcd, edm::Transition::BeginRun>()),
0091       geomToken_(iC.esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0092       theConfig(cfg),
0093       theMode(this->decodeMode(theConfig.getUntrackedParameter<std::string>("mode"))),
0094       theDir(theConfig.getUntrackedParameter<std::string>("fileDir")),
0095       theAlignmentParameterStore(nullptr),
0096       theAlignables(),
0097       theTrajectoryFactory(
0098           TrajectoryFactoryPlugin::get()->create(theConfig.getParameter<edm::ParameterSet>("TrajectoryFactory")
0099                                                      .getParameter<std::string>("TrajectoryFactoryName"),
0100                                                  theConfig.getParameter<edm::ParameterSet>("TrajectoryFactory"),
0101                                                  iC)),
0102       theMinNumHits(cfg.getParameter<unsigned int>("minNumHits")),
0103       theMaximalCor2D(cfg.getParameter<double>("max2Dcorrelation")),
0104       firstIOV_(cfg.getUntrackedParameter<AlignmentAlgorithmBase::RunNumber>("firstIOV")),
0105       ignoreFirstIOVCheck_(cfg.getUntrackedParameter<bool>("ignoreFirstIOVCheck")),
0106       enableAlignableUpdates_(cfg.getUntrackedParameter<bool>("enableAlignableUpdates")),
0107       theLastWrittenIov(0),
0108       theGblDoubleBinary(cfg.getParameter<bool>("doubleBinary")),
0109       runAtPCL_(cfg.getParameter<bool>("runAtPCL")),
0110       ignoreHitsWithoutGlobalDerivatives_(cfg.getParameter<bool>("ignoreHitsWithoutGlobalDerivatives")),
0111       skipGlobalPositionRcdCheck_(cfg.getParameter<bool>("skipGlobalPositionRcdCheck")),
0112       uniqueRunRanges_(align::makeUniqueRunRanges(cfg.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection"),
0113                                                   cond::timeTypeSpecs[cond::runnumber].beginValue)),
0114       enforceSingleIOVInput_(!(enableAlignableUpdates_ && areIOVsSpecified())),
0115       lastProcessedRun_(cond::timeTypeSpecs[cond::runnumber].beginValue) {
0116   if (!theDir.empty() && theDir.find_last_of('/') != theDir.size() - 1)
0117     theDir += '/';  // may need '/'
0118   edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm"
0119                             << "Start in mode '" << theConfig.getUntrackedParameter<std::string>("mode")
0120                             << "' with output directory '" << theDir << "'.";
0121   if (this->isMode(myMilleBit)) {
0122     theMille = std::make_unique<Mille>(
0123         (theDir + theConfig.getParameter<std::string>("binaryFile")).c_str());  // add ', false);' for text output);
0124     // use same file for GBL
0125     theBinary = std::make_unique<MilleBinary>((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str(),
0126                                               theGblDoubleBinary);
0127   }
0128 }
0129
0130 // Destructor ----------------------------------------------------------------
0131 //____________________________________________________
0132 MillePedeAlignmentAlgorithm::~MillePedeAlignmentAlgorithm() {}
0133
0134 // Call at beginning of job ---------------------------------------------------
0135 //____________________________________________________
0136 void MillePedeAlignmentAlgorithm::initialize(const edm::EventSetup &setup,
0137                                              AlignableTracker *tracker,
0138                                              AlignableMuon *muon,
0139                                              AlignableExtras *extras,
0140                                              AlignmentParameterStore *store) {
0141   if (muon) {
0142     edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
0143                                  << "Running with AlignabeMuon not yet tested.";
0144   }
0145
0146   if (!runAtPCL_ && enforceSingleIOVInput_) {
0147     const auto MIN_VAL = cond::timeTypeSpecs[cond::runnumber].beginValue;
0148     const auto MAX_VAL = cond::timeTypeSpecs[cond::runnumber].endValue;
0149     const auto &iov_alignments = setup.get<TrackerAlignmentRcd>().validityInterval();
0150     const auto &iov_surfaces = setup.get<TrackerSurfaceDeformationRcd>().validityInterval();
0151     const auto &iov_errors = setup.get<TrackerAlignmentErrorExtendedRcd>().validityInterval();
0152
0153     std::ostringstream message;
0154     bool throwException{false};
0155     if (iov_alignments.first().eventID().run() != MIN_VAL || iov_alignments.last().eventID().run() != MAX_VAL) {
0156       message << "\nTrying to apply " << setup.get<TrackerAlignmentRcd>().key().name()
0157               << " with multiple IOVs in tag without specifying 'RunRangeSelection'.\n"
0158               << "Validity range is " << iov_alignments.first().eventID().run() << " - "
0159               << iov_alignments.last().eventID().run() << "\n";
0160       throwException = true;
0161     }
0162     if (iov_surfaces.first().eventID().run() != MIN_VAL || iov_surfaces.last().eventID().run() != MAX_VAL) {
0163       message << "\nTrying to apply " << setup.get<TrackerSurfaceDeformationRcd>().key().name()
0164               << " with multiple IOVs in tag without specifying 'RunRangeSelection'.\n"
0165               << "Validity range is " << iov_surfaces.first().eventID().run() << " - "
0166               << iov_surfaces.last().eventID().run() << "\n";
0167       throwException = true;
0168     }
0169     if (iov_errors.first().eventID().run() != MIN_VAL || iov_errors.last().eventID().run() != MAX_VAL) {
0170       message << "\nTrying to apply " << setup.get<TrackerAlignmentErrorExtendedRcd>().key().name()
0171               << " with multiple IOVs in tag without specifying 'RunRangeSelection'.\n"
0172               << "Validity range is " << iov_errors.first().eventID().run() << " - "
0173               << iov_errors.last().eventID().run() << "\n";
0174       throwException = true;
0175     }
0176     if (throwException) {
0177       throw cms::Exception("DatabaseError") << "@SUB=MillePedeAlignmentAlgorithm::initialize" << message.str();
0178     }
0179   }
0180
0181   //Retrieve tracker topology from geometry
0182   const TrackerTopology *const tTopo = &setup.getData(topoToken_);
0183
0184   //Retrieve the thresolds cuts from DB for the PCL
0185   if (runAtPCL_) {
0186     const auto &th = &setup.getData(aliThrToken_);
0187     theThresholds = std::make_shared<AlignPCLThresholdsHG>();
0188     storeThresholds(th->getNrecords(), th->getThreshold_Map(), th->getFloatMap());
0189
0190     //Retrieve tracker geometry
0191     const TrackerGeometry *tGeom = &setup.getData(geomToken_);
0192     //Retrieve PixelTopologyMap
0193     pixelTopologyMap = std::make_shared<PixelTopologyMap>(tGeom, tTopo);
0194   }
0195
0196   theAlignableNavigator = std::make_unique<AlignableNavigator>(extras, tracker, muon);
0197   theAlignmentParameterStore = store;
0198   theAlignables = theAlignmentParameterStore->alignables();
0199
0200   edm::ParameterSet pedeLabelerCfg(theConfig.getParameter<edm::ParameterSet>("pedeLabeler"));
0201   edm::VParameterSet RunRangeSelectionVPSet(theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection"));
0202   pedeLabelerCfg.addUntrackedParameter<edm::VParameterSet>("RunRangeSelection", RunRangeSelectionVPSet);
0203
0204   std::string labelerPlugin = "PedeLabeler";
0205   if (!RunRangeSelectionVPSet.empty()) {
0206     labelerPlugin = "RunRangeDependentPedeLabeler";
0207     if (pedeLabelerCfg.exists("plugin")) {
0208       std::string labelerPluginCfg = pedeLabelerCfg.getParameter<std::string>("plugin");
0209       if ((labelerPluginCfg != "PedeLabeler" && labelerPluginCfg != "RunRangeDependentPedeLabeler") ||
0210           !pedeLabelerCfg.getUntrackedParameter<edm::VParameterSet>("parameterInstances").empty()) {
0211         throw cms::Exception("BadConfig") << "MillePedeAlignmentAlgorithm::initialize"
0212                                           << "both RunRangeSelection and generic labeler specified in config file. "
0213                                           << "Please get rid of either one of them.\n";
0214       }
0215     }
0216   } else {
0217     if (pedeLabelerCfg.exists("plugin")) {
0218       labelerPlugin = pedeLabelerCfg.getParameter<std::string>("plugin");
0219     }
0220   }
0221
0222   if (!pedeLabelerCfg.exists("plugin")) {
0223     pedeLabelerCfg.addUntrackedParameter<std::string>("plugin", labelerPlugin);
0224   }
0225
0226   edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
0227                             << "Using plugin '" << labelerPlugin << "' to generate labels.";
0228
0229   thePedeLabels = std::shared_ptr<PedeLabelerBase>(PedeLabelerPluginFactory::get()->create(
0230       labelerPlugin, PedeLabelerBase::TopLevelAlignables(tracker, muon, extras), pedeLabelerCfg));
0231
0232   // 1) Create PedeSteerer: correct alignable positions for coordinate system selection
0233   edm::ParameterSet pedeSteerCfg(theConfig.getParameter<edm::ParameterSet>("pedeSteerer"));
0234   thePedeSteer = std::make_unique<PedeSteerer>(tracker,
0235                                                muon,
0236                                                extras,
0237                                                theAlignmentParameterStore,
0238                                                thePedeLabels.get(),
0239                                                pedeSteerCfg,
0240                                                theDir,
0241                                                !this->isMode(myPedeSteerBit));
0242
0243   // 2) If requested, directly read in and apply result of previous pede run,
0244   //    assuming that correction from 1) was also applied to create the result:
0245   const std::vector<edm::ParameterSet> mprespset(
0246       theConfig.getParameter<std::vector<edm::ParameterSet> >("pedeReaderInputs"));
0247   if (!mprespset.empty()) {
0248     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
0249                               << "Apply " << mprespset.end() - mprespset.begin()
0250                               << " previous MillePede constants from 'pedeReaderInputs'.";
0251   }
0252
0253   // FIXME: add selection of run range via 'pedeReaderInputs'
0254   // Note: Results for parameters of IntegratedCalibration's cannot be treated...
0255   RunRange runrange(cond::timeTypeSpecs[cond::runnumber].beginValue, cond::timeTypeSpecs[cond::runnumber].endValue);
0256   for (std::vector<edm::ParameterSet>::const_iterator iSet = mprespset.begin(), iE = mprespset.end(); iSet != iE;
0257        ++iSet) {
0258     // This read will ignore calibrations as long as they are not yet passed to Millepede
0259     // during/before initialize(..) - currently addCalibrations(..) is called later in AlignmentProducer
0260     if (!this->readFromPede((*iSet), false, runrange)) {  // false: do not erase SelectionUserVariables
0261       throw cms::Exception("BadConfig")
0262           << "MillePedeAlignmentAlgorithm::initialize: Problems reading input constants of "
0263           << "pedeReaderInputs entry " << iSet - mprespset.begin() << '.';
0264     }
0265     theAlignmentParameterStore->applyParameters();
0266     // Needed to shut up later warning from checkAliParams:
0267     theAlignmentParameterStore->resetParameters();
0268   }
0269
0270   // 3) Now create steerings with 'final' start position:
0271   thePedeSteer->buildSubSteer(tracker, muon, extras);
0272
0273   // After (!) 1-3 of PedeSteerer which uses the SelectionUserVariables attached to the parameters:
0274   this->buildUserVariables(theAlignables);  // for hit statistics and/or pede result
0275
0276   if (this->isMode(myMilleBit)) {
0277     if (!theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles").empty() ||
0278         !theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles").empty()) {
0279       throw cms::Exception("BadConfig") << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' must be empty for "
0280                                         << "modes running mille.";
0281     }
0282     const std::string moniFile(theConfig.getUntrackedParameter<std::string>("monitorFile"));
0283     if (!moniFile.empty())
0284       theMonitor = std::make_unique<MillePedeMonitor>(tTopo, (theDir + moniFile).c_str());
0285
0286     // Get trajectory factory. In case nothing found, FrameWork will throw...
0287   }
0288
0289   if (this->isMode(myPedeSteerBit)) {
0290     // Get config for survey and set flag accordingly
0291     const edm::ParameterSet pxbSurveyCfg(theConfig.getParameter<edm::ParameterSet>("surveyPixelBarrel"));
0292     theDoSurveyPixelBarrel = pxbSurveyCfg.getParameter<bool>("doSurvey");
0293     if (theDoSurveyPixelBarrel)
0294       this->addPxbSurvey(pxbSurveyCfg);
0295   }
0296 }
0297
0298 //____________________________________________________
0299 bool MillePedeAlignmentAlgorithm::supportsCalibrations() { return true; }
0300
0301 //____________________________________________________
0302 bool MillePedeAlignmentAlgorithm::addCalibrations(const std::vector<IntegratedCalibrationBase *> &iCals) {
0303   theCalibrations.insert(theCalibrations.end(), iCals.begin(), iCals.end());
0304   thePedeLabels->addCalibrations(iCals);
0305   return true;
0306 }
0307
0308 //____________________________________________________
0309 bool MillePedeAlignmentAlgorithm::storeThresholds(const int &nRecords,
0310                                                   const AlignPCLThresholdsHG::threshold_map &thresholdMap,
0311                                                   const AlignPCLThresholdsHG::param_map &floatMap) {
0312   theThresholds->setAlignPCLThresholds(nRecords, thresholdMap);
0313   theThresholds->setFloatMap(floatMap);
0314   return true;
0315 }
0316
0317 //_____________________________________________________________________________
0318 bool MillePedeAlignmentAlgorithm::processesEvents() {
0319   if (isMode(myMilleBit)) {
0320     return true;
0321   } else {
0322     return false;
0323   }
0324 }
0325
0326 //_____________________________________________________________________________
0327 bool MillePedeAlignmentAlgorithm::storeAlignments() {
0328   if (isMode(myPedeReadBit)) {
0329     if (runAtPCL_) {
0330       MillePedeFileReader mpReader(theConfig.getParameter<edm::ParameterSet>("MillePedeFileReader"),
0331                                    thePedeLabels,
0332                                    theThresholds,
0333                                    pixelTopologyMap);
0334       mpReader.read();
0335       return mpReader.storeAlignments();
0336     } else {
0337       return true;
0338     }
0339   } else {
0340     return false;
0341   }
0342 }
0343
0344 //____________________________________________________
0345 bool MillePedeAlignmentAlgorithm::setParametersForRunRange(const RunRange &runrange) {
0346   if (this->isMode(myPedeReadBit)) {
0347     if (not theAlignmentParameterStore) {
0348       return false;
0349     }
0350     // restore initial positions, rotations and deformations
0351     if (enableAlignableUpdates_) {
0352       theAlignmentParameterStore->restoreCachedTransformations(runrange.first);
0353     } else {
0354       theAlignmentParameterStore->restoreCachedTransformations();
0355     }
0356
0357     // Needed to shut up later warning from checkAliParams:
0358     theAlignmentParameterStore->resetParameters();
0359     // To avoid that they keep values from previous IOV if no new one in pede result
0360     this->buildUserVariables(theAlignables);
0361
0362     if (!this->readFromPede(theConfig.getParameter<edm::ParameterSet>("pedeReader"), true, runrange)) {
0363       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::setParametersForRunRange"
0364                                  << "Problems reading pede result, but applying!";
0365     }
0366     theAlignmentParameterStore->applyParameters();
0367
0368     this->doIO(++theLastWrittenIov);  // pre-increment!
0369   }
0370
0371   return true;
0372 }
0373
0374 // Call at end of job ---------------------------------------------------------
0375 //____________________________________________________
0376 void MillePedeAlignmentAlgorithm::terminate(const edm::EventSetup &iSetup) { terminate(); }
0377 void MillePedeAlignmentAlgorithm::terminate() {
0378   theMille.reset();  // delete to close binary before running pede below (flush would be enough...)
0379   theBinary.reset();
0380
0381   std::vector<std::string> files;
0382   if (this->isMode(myMilleBit) || !theConfig.getParameter<std::string>("binaryFile").empty()) {
0383     files.push_back(theDir + theConfig.getParameter<std::string>("binaryFile"));
0384   } else {
0385     const std::vector<std::string> plainFiles(theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
0386     files = getExistingFormattedFiles(plainFiles, theDir);
0387     // Do some logging:
0388     std::string filesForLogOutput;
0389     for (const auto &file : files)
0390       filesForLogOutput += " " + file + ",";
0391     if (filesForLogOutput.length() != 0)
0392       filesForLogOutput.pop_back();
0393     edm::LogInfo("Alignment") << "Based on the config parameter mergeBinaryFiles, using the following "
0394                               << "files as input (assigned weights are indicated by ' -- <weight>'):"
0395                               << filesForLogOutput;
0396   }
0397
0398   if (not theAlignmentParameterStore)
0399     return;
0400
0401   // cache all positions, rotations and deformations
0402   theAlignmentParameterStore->cacheTransformations();
0403   if (this->isMode(myPedeReadBit) && enableAlignableUpdates_) {
0404     if (lastProcessedRun_ < uniqueRunRanges_.back().first) {
0405       throw cms::Exception("BadConfig") << "@SUB=MillePedeAlignmentAlgorithm::terminate\n"
0406                                         << "Last IOV of 'RunRangeSelection' has not been processed. "
0407                                         << "Please reconfigure your source to process the runs at least up to "
0408                                         << uniqueRunRanges_.back().first << ".";
0409     }
0410     auto lastCachedRun = uniqueRunRanges_.front().first;
0411     for (const auto &runRange : uniqueRunRanges_) {
0412       const auto run = runRange.first;
0413       if (std::find(cachedRuns_.begin(), cachedRuns_.end(), run) == cachedRuns_.end()) {
0414         theAlignmentParameterStore->restoreCachedTransformations(lastCachedRun);
0415         theAlignmentParameterStore->cacheTransformations(run);
0416       } else {
0417         lastCachedRun = run;
0418       }
0419     }
0420     theAlignmentParameterStore->restoreCachedTransformations();
0421   }
0422
0423   const std::string masterSteer(thePedeSteer->buildMasterSteer(files));  // do only if myPedeSteerBit?
0424   if (this->isMode(myPedeRunBit)) {
0425     thePedeSteer->runPede(masterSteer);
0426   }
0427
0428   // parameters from pede are not yet applied,
0429   // so we can still write start positions (but with hit statistics in case of mille):
0430   this->doIO(0);
0431   theLastWrittenIov = 0;
0432 }
0433
0434 std::vector<std::string> MillePedeAlignmentAlgorithm::getExistingFormattedFiles(
0435     const std::vector<std::string> &plainFiles, const std::string &theDir) {
0436   std::vector<std::string> files;
0437   for (const auto &plainFile : plainFiles) {
0438     const std::string &theInputFileName = plainFile;
0439     int theNumber = 0;
0440     while (true) {
0441       // Create a formatted version of the filename, with growing numbers
0442       // If the parameter doesn't contain a formatting directive, it just stays unchanged
0443       char theNumberedInputFileName[200];
0444       sprintf(theNumberedInputFileName, theInputFileName.c_str(), theNumber);
0445       std::string theCompleteInputFileName = theDir + theNumberedInputFileName;
0446       const auto endOfStrippedFileName = theCompleteInputFileName.rfind(" --");
0447       const auto strippedInputFileName = theCompleteInputFileName.substr(0, endOfStrippedFileName);
0448       // Check if the file exists
0449       struct stat buffer;
0450       if (stat(strippedInputFileName.c_str(), &buffer) == 0) {
0451         // If the file exists, add it to the list
0452         files.push_back(theCompleteInputFileName);
0453         if (theNumberedInputFileName == theInputFileName) {
0454           // If the filename didn't contain a formatting directive, no reason to look any further, break out of the loop
0455           break;
0456         } else {
0457           // Otherwise look for the next number
0458           theNumber++;
0459         }
0460       } else {
0461         // The file doesn't exist, break out of the loop
0462         break;
0463       }
0464     }
0465     // warning if unformatted (-> theNumber stays at 0) does not exist
0466     if (theNumber == 0 && (files.empty() || files.back() != plainFile)) {
0467       edm::LogWarning("Alignment") << "The input file '" << plainFile << "' does not exist.";
0468     }
0469   }
0470   return files;
0471 }
0472
0473 // Run the algorithm on trajectories and tracks -------------------------------
0474 //____________________________________________________
0475 void MillePedeAlignmentAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo) {
0476   if (!this->isMode(myMilleBit))
0477     return;  // no theMille created...
0478   const auto &tracks = eventInfo.trajTrackPairs();
0479
0480   if (theMonitor) {  // monitor input tracks
0481     for (const auto &iTrajTrack : tracks) {
0482       theMonitor->fillTrack(iTrajTrack.second);
0483     }
0484   }
0485
0486   const RefTrajColl trajectories(theTrajectoryFactory->trajectories(setup, tracks, eventInfo.beamSpot()));
0487
0488   // Now loop over ReferenceTrajectoryCollection
0489   unsigned int refTrajCount = 0;  // counter for track monitoring
0490   const auto tracksPerTraj = theTrajectoryFactory->tracksPerTrajectory();
0491   for (auto iRefTraj = trajectories.cbegin(), iRefTrajE = trajectories.cend(); iRefTraj != iRefTrajE;
0492        ++iRefTraj, ++refTrajCount) {
0493     const RefTrajColl::value_type &refTrajPtr = *iRefTraj;
0494     if (theMonitor)
0495       theMonitor->fillRefTrajectory(refTrajPtr);
0496
0497     const auto nHitXy = this->addReferenceTrajectory(setup, eventInfo, refTrajPtr);
0498
0499     if (theMonitor && (nHitXy.first || nHitXy.second)) {
0500       // if trajectory used (i.e. some hits), fill monitoring
0501       const auto offset = tracksPerTraj * refTrajCount;
0502       for (unsigned int iTrack = 0; iTrack < tracksPerTraj; ++iTrack) {
0503         theMonitor->fillUsedTrack(tracks[offset + iTrack].second, nHitXy.first, nHitXy.second);
0504       }
0505     }
0506
0507   }  // end of reference trajectory and track loop
0508 }
0509
0510 //____________________________________________________
0511 std::pair<unsigned int, unsigned int> MillePedeAlignmentAlgorithm::addReferenceTrajectory(
0512     const edm::EventSetup &setup, const EventInfo &eventInfo, const RefTrajColl::value_type &refTrajPtr) {
0513   std::pair<unsigned int, unsigned int> hitResultXy(0, 0);
0514   if (refTrajPtr->isValid()) {
0515     // GblTrajectory?
0516     if (!refTrajPtr->gblInput().empty()) {
0517       // by construction: number of GblPoints == number of recHits or == zero !!!
0518       unsigned int iHit = 0;
0519       unsigned int numPointsWithMeas = 0;
0520       std::vector<GblPoint>::iterator itPoint;
0521       auto theGblInput = refTrajPtr->gblInput();
0522       for (unsigned int iTraj = 0; iTraj < refTrajPtr->gblInput().size(); ++iTraj) {
0523         for (itPoint = refTrajPtr->gblInput()[iTraj].first.begin(); itPoint < refTrajPtr->gblInput()[iTraj].first.end();
0524              ++itPoint) {
0525           if (this->addGlobalData(setup, eventInfo, refTrajPtr, iHit++, *itPoint) < 0)
0526             return hitResultXy;
0527           if (itPoint->numMeasurements() >= 1)
0528             ++numPointsWithMeas;
0529         }
0530       }
0531       hitResultXy.first = numPointsWithMeas;
0532       // check #hits criterion
0533       if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits)
0534         return hitResultXy;
0535       // construct GBL trajectory
0536       if (refTrajPtr->gblInput().size() == 1) {
0537         // from single track
0538         GblTrajectory aGblTrajectory(refTrajPtr->gblInput()[0].first, refTrajPtr->nominalField() != 0);
0539         // GBL fit trajectory
0540         /*double Chi2;
0541         int Ndf;
0542         double lostWeight;
0543         aGblTrajectory.fit(Chi2, Ndf, lostWeight);
0544         std::cout << " GblFit: " << Chi2 << ", " << Ndf << ", " << lostWeight << std::endl; */
0545         // write to MP binary file
0546         if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
0547           aGblTrajectory.milleOut(*theBinary);
0548       }
0549       if (refTrajPtr->gblInput().size() == 2) {
0550         // from TwoBodyDecay
0551         GblTrajectory aGblTrajectory(refTrajPtr->gblInput(),
0552                                      refTrajPtr->gblExtDerivatives(),
0553                                      refTrajPtr->gblExtMeasurements(),
0554                                      refTrajPtr->gblExtPrecisions());
0555         // write to MP binary file
0556         if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
0557           aGblTrajectory.milleOut(*theBinary);
0558       }
0559     } else {
0560       // to add hits if all fine:
0561       std::vector<AlignmentParameters *> parVec(refTrajPtr->recHits().size());
0562       // collect hit statistics, assuming that there are no y-only hits
0563       std::vector<bool> validHitVecY(refTrajPtr->recHits().size(), false);
0564       // Use recHits from ReferenceTrajectory (since they have the right order!):
0565       for (unsigned int iHit = 0; iHit < refTrajPtr->recHits().size(); ++iHit) {
0566         const int flagXY = this->addMeasurementData(setup, eventInfo, refTrajPtr, iHit, parVec[iHit]);
0567
0568         if (flagXY < 0) {  // problem
0569           hitResultXy.first = 0;
0570           break;
0571         } else {  // hit is fine, increase x/y statistics
0572           if (flagXY >= 1)
0573             ++hitResultXy.first;
0574           validHitVecY[iHit] = (flagXY >= 2);
0575         }
0576       }  // end loop on hits
0577
0578       // add virtual measurements
0579       for (unsigned int iVirtualMeas = 0; iVirtualMeas < refTrajPtr->numberOfVirtualMeas(); ++iVirtualMeas) {
0580         this->addVirtualMeas(refTrajPtr, iVirtualMeas);
0581       }
0582
0583       // kill or end 'track' for mille, depends on #hits criterion
0584       if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits) {
0585         theMille->kill();
0586         hitResultXy.first = hitResultXy.second = 0;  //reset
0587       } else {
0588         theMille->end();
0589         // add x/y hit count to MillePedeVariables of parVec,
0590         // returning number of y-hits of the reference trajectory
0591         hitResultXy.second = this->addHitCount(parVec, validHitVecY);
0592         //
0593       }
0594     }
0595
0596   }  // end if valid trajectory
0597
0598   return hitResultXy;
0599 }
0600
0601 //____________________________________________________
0602 unsigned int MillePedeAlignmentAlgorithm::addHitCount(const std::vector<AlignmentParameters *> &parVec,
0603                                                       const std::vector<bool> &validHitVecY) const {
0604   // Loop on all hit information in the input arrays and count valid y-hits:
0605   unsigned int nHitY = 0;
0606   for (unsigned int iHit = 0; iHit < validHitVecY.size(); ++iHit) {
0607     Alignable *ali = (parVec[iHit] ? parVec[iHit]->alignable() : nullptr);
0608     // Loop upwards on hierarchy of alignables to add hits to all levels
0609     // that are currently aligned. If only a non-selected alignable was hit,
0610     // (i.e. flagXY == 0 in addReferenceTrajectory(..)), there is no loop at all...
0611     while (ali) {
0612       AlignmentParameters *pars = ali->alignmentParameters();
0613       if (pars) {  // otherwise hierarchy level not selected
0614         // cast ensured by previous checks:
0615         MillePedeVariables *mpVar = static_cast<MillePedeVariables *>(pars->userVariables());
0616         // every hit has an x-measurement, cf. addReferenceTrajectory(..):
0617         mpVar->increaseHitsX();
0618         if (validHitVecY[iHit]) {
0619           mpVar->increaseHitsY();
0620           if (pars == parVec[iHit])
0621             ++nHitY;  // do not count hits twice
0622         }
0623       }
0624       ali = ali->mother();
0625     }
0626   }
0627
0628   return nHitY;
0629 }
0630
0631 void MillePedeAlignmentAlgorithm::beginRun(const edm::Run &run, const edm::EventSetup &setup, bool changed) {
0632   if (run.run() < firstIOV_ && !ignoreFirstIOVCheck_) {
0633     throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::beginRun\n"
0634                                       << "Using data (run = " << run.run() << ") prior to the first defined IOV ("
0635                                       << firstIOV_ << ").";
0636   }
0637
0638   lastProcessedRun_ = run.run();
0639
0640   if (changed && enableAlignableUpdates_) {
0641     const auto runNumber = run.run();
0642     auto firstRun = cond::timeTypeSpecs[cond::runnumber].beginValue;
0643     for (auto runRange = uniqueRunRanges_.crbegin(); runRange != uniqueRunRanges_.crend(); ++runRange) {
0644       if (runNumber >= runRange->first) {
0645         firstRun = runRange->first;
0646         break;
0647       }
0648     }
0649     if (std::find(cachedRuns_.begin(), cachedRuns_.end(), firstRun) != cachedRuns_.end()) {
0650       const auto &geometryRcd = setup.get<IdealGeometryRecord>();
0651       const auto &globalPosRcd = setup.get<GlobalPositionRcd>();
0652       const auto &alignmentRcd = setup.get<TrackerAlignmentRcd>();
0653       const auto &surfaceRcd = setup.get<TrackerSurfaceDeformationRcd>();
0654       const auto &errorRcd = setup.get<TrackerAlignmentErrorExtendedRcd>();
0655
0656       std::ostringstream message;
0657       bool throwException{false};
0658       message << "Trying to cache tracker alignment payloads for a run (" << runNumber << ") in an IOV (" << firstRun
0659               << ") that was already cached.\n"
0660               << "The following records in your input database tag have an IOV "
0661               << "boundary that does not match your IOV definition:\n";
0662       if (geometryRcd.validityInterval().first().eventID().run() > firstRun) {
0663         message << " - IdealGeometryRecord '" << geometryRcd.key().name() << "' (since "
0664                 << geometryRcd.validityInterval().first().eventID().run() << ")\n";
0665         throwException = true;
0666       }
0667       if (globalPosRcd.validityInterval().first().eventID().run() > firstRun) {
0668         message << " - GlobalPositionRecord '" << globalPosRcd.key().name() << "' (since "
0669                 << globalPosRcd.validityInterval().first().eventID().run() << ")";
0670         if (skipGlobalPositionRcdCheck_) {
0671           message << " --> ignored\n";
0672         } else {
0673           message << "\n";
0674           throwException = true;
0675         }
0676       }
0677       if (alignmentRcd.validityInterval().first().eventID().run() > firstRun) {
0678         message << " - TrackerAlignmentRcd '" << alignmentRcd.key().name() << "' (since "
0679                 << alignmentRcd.validityInterval().first().eventID().run() << ")\n";
0680         throwException = true;
0681       }
0682       if (surfaceRcd.validityInterval().first().eventID().run() > firstRun) {
0683         message << " - TrackerSurfaceDeformationRcd '" << surfaceRcd.key().name() << "' (since "
0684                 << surfaceRcd.validityInterval().first().eventID().run() << ")\n";
0685         throwException = true;
0686       }
0687       if (errorRcd.validityInterval().first().eventID().run() > firstRun) {
0688         message << " - TrackerAlignmentErrorExtendedRcd '" << errorRcd.key().name() << "' (since "
0689                 << errorRcd.validityInterval().first().eventID().run() << ")\n";
0690         throwException = true;
0691       }
0692
0693       if (throwException) {
0694         throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::beginRun\n" << message.str();
0695       }
0696     } else {
0697       cachedRuns_.push_back(firstRun);
0698       theAlignmentParameterStore->cacheTransformations(firstRun);
0699     }
0700   }
0701 }
0702
0703 //____________________________________________________
0704 void MillePedeAlignmentAlgorithm::endRun(const EventInfo &eventInfo,
0705                                          const EndRunInfo &runInfo,
0706                                          const edm::EventSetup &setup) {
0707   if (runInfo.tkLasBeams() && runInfo.tkLasBeamTsoses()) {
0708     // LAS beam treatment
0709     this->addLaserData(eventInfo, *(runInfo.tkLasBeams()), *(runInfo.tkLasBeamTsoses()));
0710   }
0711   if (this->isMode(myMilleBit))
0712     theMille->flushOutputFile();
0713 }
0714
0715 // Implementation of endRun that DOES get called. (Because we need it.)
0716 void MillePedeAlignmentAlgorithm::endRun(const EndRunInfo &runInfo, const edm::EventSetup &setup) {
0717   if (this->isMode(myMilleBit))
0718     theMille->flushOutputFile();
0719 }
0720
0721 //____________________________________________________
0722 void MillePedeAlignmentAlgorithm::beginLuminosityBlock(const edm::EventSetup &) {
0723   if (!runAtPCL_)
0724     return;
0725   if (this->isMode(myMilleBit))
0726     theMille->resetOutputFile();
0727 }
0728
0729 //____________________________________________________
0730 void MillePedeAlignmentAlgorithm::endLuminosityBlock(const edm::EventSetup &) {
0731   if (!runAtPCL_)
0732     return;
0733   if (this->isMode(myMilleBit))
0734     theMille->flushOutputFile();
0735 }
0736
0737 //____________________________________________________
0738 int MillePedeAlignmentAlgorithm::addMeasurementData(const edm::EventSetup &setup,
0739                                                     const EventInfo &eventInfo,
0740                                                     const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0741                                                     unsigned int iHit,
0742                                                     AlignmentParameters *&params) {
0743   params = nullptr;
0744   theFloatBufferX.clear();
0745   theFloatBufferY.clear();
0746   theIntBuffer.clear();
0747
0748   const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
0749   const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
0750   // ignore invalid hits
0751   if (!recHitPtr->isValid())
0752     return 0;
0753
0754   // First add the derivatives from IntegratedCalibration's,
0755   // should even be OK if problems for "usual" derivatives from Alignables
0756   this->globalDerivativesCalibration(recHitPtr,
0757                                      tsos,
0758                                      setup,
0759                                      eventInfo,  // input
0760                                      theFloatBufferX,
0761                                      theFloatBufferY,
0762                                      theIntBuffer);  // output
0763
0764   // get AlignableDet/Unit for this hit
0765   AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0766
0767   if (!this->globalDerivativesHierarchy(eventInfo,
0768                                         tsos,
0769                                         alidet,
0770                                         alidet,
0771                                         theFloatBufferX,  // 2x alidet, sic!
0772                                         theFloatBufferY,
0773                                         theIntBuffer,
0774                                         params)) {
0775     return -1;  // problem
0776   } else if (theFloatBufferX.empty() && ignoreHitsWithoutGlobalDerivatives_) {
0777     return 0;  // empty for X: no alignable for hit, nor calibrations
0778   } else {
0779     // store measurement even if no alignable or calibrations
0780     // -> measurement used for pede-internal track-fit
0781     return this->callMille(refTrajPtr, iHit, theIntBuffer, theFloatBufferX, theFloatBufferY);
0782   }
0783 }
0784
0785 //____________________________________________________
0786
0787 int MillePedeAlignmentAlgorithm::addGlobalData(const edm::EventSetup &setup,
0788                                                const EventInfo &eventInfo,
0789                                                const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0790                                                unsigned int iHit,
0791                                                GblPoint &gblPoint) {
0792   AlignmentParameters *params = nullptr;
0793   std::vector<double> theDoubleBufferX, theDoubleBufferY;
0794   theDoubleBufferX.clear();
0795   theDoubleBufferY.clear();
0796   theIntBuffer.clear();
0797   int iret = 0;
0798
0799   const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
0800   const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
0801   // ignore invalid hits
0802   if (!recHitPtr->isValid())
0803     return 0;
0804
0805   // get AlignableDet/Unit for this hit
0806   AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0807
0808   if (!this->globalDerivativesHierarchy(eventInfo,
0809                                         tsos,
0810                                         alidet,
0811                                         alidet,
0812                                         theDoubleBufferX,  // 2x alidet, sic!
0813                                         theDoubleBufferY,
0814                                         theIntBuffer,
0815                                         params)) {
0816     return -1;  // problem
0817   }
0818   //calibration parameters
0819   int globalLabel;
0820   std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0821   for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0822     // get all derivatives of this calibration // const unsigned int num =
0823     (*iCalib)->derivatives(derivs, *recHitPtr, tsos, setup, eventInfo);
0824     for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0825       // transfer label and x/y derivatives
0826       globalLabel = thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second);
0827       if (globalLabel > 0 && globalLabel <= 2147483647) {
0828         theIntBuffer.push_back(globalLabel);
0829         theDoubleBufferX.push_back(iValuesInd->first.first);
0830         theDoubleBufferY.push_back(iValuesInd->first.second);
0831       } else {
0832         edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addGlobalData"
0833                                    << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
0834       }
0835     }
0836   }
0837   unsigned int numGlobals = theIntBuffer.size();
0838   if (numGlobals > 0) {
0839     Eigen::Matrix<double, 2, Eigen::Dynamic> globalDer{2, numGlobals};
0840     for (unsigned int i = 0; i < numGlobals; ++i) {
0841       globalDer(0, i) = theDoubleBufferX[i];
0842       globalDer(1, i) = theDoubleBufferY[i];
0843     }
0844     gblPoint.addGlobals(theIntBuffer, globalDer);
0845     iret = 1;
0846   }
0847   return iret;
0848 }
0849
0850 //____________________________________________________
0851 bool MillePedeAlignmentAlgorithm ::globalDerivativesHierarchy(const EventInfo &eventInfo,
0852                                                               const TrajectoryStateOnSurface &tsos,
0853                                                               Alignable *ali,
0854                                                               const AlignableDetOrUnitPtr &alidet,
0855                                                               std::vector<float> &globalDerivativesX,
0856                                                               std::vector<float> &globalDerivativesY,
0857                                                               std::vector<int> &globalLabels,
0858                                                               AlignmentParameters *&lowestParams) const {
0859   // derivatives and labels are recursively attached
0860   if (!ali)
0861     return true;  // no mother might be OK
0862
0863   if (false && theMonitor && alidet != ali)
0864     theMonitor->fillFrameToFrame(alidet, ali);
0865
0866   AlignmentParameters *params = ali->alignmentParameters();
0867
0868   if (params) {
0869     if (!lowestParams)
0870       lowestParams = params;  // set parameters of lowest level
0871
0872     bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0873     const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0874
0875     if (0 == alignableLabel) {  // FIXME: what about regardAllHits in Markus' code?
0876       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0877                                    << "Label not found, skip Alignable.";
0878       return false;
0879     }
0880
0881     const std::vector<bool> &selPars = params->selector();
0882     const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
0883
0884     // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
0885     for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
0886       if (selPars[iSel]) {
0887         globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
0888         if (hasSplitParameters == true) {
0889           globalLabels.push_back(thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos));
0890         } else {
0891           globalLabels.push_back(thePedeLabels->parameterLabel(alignableLabel, iSel));
0892         }
0893         globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
0894       }
0895     }
0896     // Exclude mothers if Alignable selected to be no part of a hierarchy:
0897     if (thePedeSteer->isNoHiera(ali))
0898       return true;
0899   }
0900   // Call recursively for mother, will stop if mother == 0:
0901   return this->globalDerivativesHierarchy(
0902       eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
0903 }
0904
0905 //____________________________________________________
0906 bool MillePedeAlignmentAlgorithm ::globalDerivativesHierarchy(const EventInfo &eventInfo,
0907                                                               const TrajectoryStateOnSurface &tsos,
0908                                                               Alignable *ali,
0909                                                               const AlignableDetOrUnitPtr &alidet,
0910                                                               std::vector<double> &globalDerivativesX,
0911                                                               std::vector<double> &globalDerivativesY,
0912                                                               std::vector<int> &globalLabels,
0913                                                               AlignmentParameters *&lowestParams) const {
0914   // derivatives and labels are recursively attached
0915   if (!ali)
0916     return true;  // no mother might be OK
0917
0918   if (false && theMonitor && alidet != ali)
0919     theMonitor->fillFrameToFrame(alidet, ali);
0920
0921   AlignmentParameters *params = ali->alignmentParameters();
0922
0923   if (params) {
0924     if (!lowestParams)
0925       lowestParams = params;  // set parameters of lowest level
0926
0927     bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0928     const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0929
0930     if (0 == alignableLabel) {  // FIXME: what about regardAllHits in Markus' code?
0931       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0932                                    << "Label not found, skip Alignable.";
0933       return false;
0934     }
0935
0936     const std::vector<bool> &selPars = params->selector();
0937     const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
0938     int globalLabel;
0939
0940     // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
0941     for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
0942       if (selPars[iSel]) {
0943         if (hasSplitParameters == true) {
0944           globalLabel = thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos);
0945         } else {
0946           globalLabel = thePedeLabels->parameterLabel(alignableLabel, iSel);
0947         }
0948         if (globalLabel > 0 && globalLabel <= 2147483647) {
0949           globalLabels.push_back(globalLabel);
0950           globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
0951           globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
0952         } else {
0953           edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0954                                      << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
0955         }
0956       }
0957     }
0958     // Exclude mothers if Alignable selected to be no part of a hierarchy:
0959     if (thePedeSteer->isNoHiera(ali))
0960       return true;
0961   }
0962   // Call recursively for mother, will stop if mother == 0:
0963   return this->globalDerivativesHierarchy(
0964       eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
0965 }
0966
0967 //____________________________________________________
0968 void MillePedeAlignmentAlgorithm::globalDerivativesCalibration(const TransientTrackingRecHit::ConstRecHitPointer &recHit,
0969                                                                const TrajectoryStateOnSurface &tsos,
0970                                                                const edm::EventSetup &setup,
0971                                                                const EventInfo &eventInfo,
0972                                                                std::vector<float> &globalDerivativesX,
0973                                                                std::vector<float> &globalDerivativesY,
0974                                                                std::vector<int> &globalLabels) const {
0975   std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0976   for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0977     // get all derivatives of this calibration // const unsigned int num =
0978     (*iCalib)->derivatives(derivs, *recHit, tsos, setup, eventInfo);
0979     for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0980       // transfer label and x/y derivatives
0981       globalLabels.push_back(thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second));
0982       globalDerivativesX.push_back(iValuesInd->first.first);
0983       globalDerivativesY.push_back(iValuesInd->first.second);
0984     }
0985   }
0986 }
0987
0988 // //____________________________________________________
0989 // void MillePedeAlignmentAlgorithm
0990 // ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0991 //             unsigned int iTrajHit, MeasurementDirection xOrY,
0992 //             const std::vector<float> &globalDerivatives, const std::vector<int> &globalLabels)
0993 // {
0994 //   const unsigned int xyIndex = iTrajHit*2 + xOrY;
0995 //   // FIXME: here for residuum and sigma we could use KALMAN-Filter results
0996 //   const float residuum =
0997 //     refTrajPtr->measurements()[xyIndex] - refTrajPtr->trajectoryPositions()[xyIndex];
0998 //   const float covariance = refTrajPtr->measurementErrors()[xyIndex][xyIndex];
0999 //   const float sigma = (covariance > 0. ? TMath::Sqrt(covariance) : 0.);
1000
1001 //   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1002
1003 //   std::vector<float> localDerivs(locDerivMatrix.num_col());
1004 //   for (unsigned int i = 0; i < localDerivs.size(); ++i) {
1005 //     localDerivs[i] = locDerivMatrix[xyIndex][i];
1006 //   }
1007
1008 //   // &(vector[0]) is valid - as long as vector is not empty
1009 //   // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1010 //   theMille->mille(localDerivs.size(), &(localDerivs[0]),
1011 //                globalDerivatives.size(), &(globalDerivatives[0]), &(globalLabels[0]),
1012 //                residuum, sigma);
1013 //   if (theMonitor) {
1014 //     theMonitor->fillDerivatives(refTrajPtr->recHits()[iTrajHit],localDerivs, globalDerivatives,
1015 //                              (xOrY == kLocalY));
1016 //     theMonitor->fillResiduals(refTrajPtr->recHits()[iTrajHit],
1017 //                            refTrajPtr->trajectoryStates()[iTrajHit],
1018 //                            iTrajHit, residuum, sigma, (xOrY == kLocalY));
1019 //   }
1020 // }
1021
1022 //____________________________________________________
1023 bool MillePedeAlignmentAlgorithm::is2D(const ConstRecHitPointer &recHit) const {
1024   // FIXME: Check whether this is a reliable and recommended way to find out...
1025
1026   if (recHit->dimension() < 2) {
1027     return false;                  // some muon and TIB/TOB stuff really has RecHit1D
1028   } else if (recHit->detUnit()) {  // detunit in strip is 1D, in pixel 2D
1029     return recHit->detUnit()->type().isTrackerPixel();
1030   } else {  // stereo strips  (FIXME: endcap trouble due to non-parallel strips (wedge sensors)?)
1031     if (dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit->hit())) {  // check persistent hit
1032       // projected: 1D measurement on 'glued' module
1033       return false;
1034     } else {
1035       return true;
1036     }
1037   }
1038 }
1039
1040 //__________________________________________________________________________________________________
1041 bool MillePedeAlignmentAlgorithm::readFromPede(const edm::ParameterSet &mprespset,
1042                                                bool setUserVars,
1043                                                const RunRange &runrange) {
1044   bool allEmpty = this->areEmptyParams(theAlignables);
1045
1046   PedeReader reader(mprespset, *thePedeSteer, *thePedeLabels, runrange);
1047   align::Alignables alis;
1048   bool okRead = reader.read(alis, setUserVars);  // also may set params of IntegratedCalibration's
1049   bool numMatch = true;
1050
1051   std::stringstream out;
1052   out << "Read " << alis.size() << " alignables";
1053   if (alis.size() != theAlignables.size()) {
1054     out << " while " << theAlignables.size() << " in store";
1055     numMatch = false;  // FIXME: Should we check one by one? Or transfer 'alis' to the store?
1056   }
1057   if (!okRead)
1058     out << ", but problems in reading";
1059   if (!allEmpty)
1060     out << ", possibly overwriting previous settings";
1061   out << ".";
1062
1063   if (okRead && allEmpty) {
1064     if (numMatch) {  // as many alignables with result as trying to align
1065       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1066     } else if (!alis.empty()) {  // dead module do not get hits and no pede result
1067       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1068     } else {  // serious problem: no result read - and not all modules can be dead...
1069       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1070       return false;
1071     }
1072     return true;
1073   }
1074   // the rest is not OK:
1075   edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1076   return false;
1077 }
1078
1079 //__________________________________________________________________________________________________
1080 bool MillePedeAlignmentAlgorithm::areEmptyParams(const align::Alignables &alignables) const {
1081   for (const auto &iAli : alignables) {
1082     const AlignmentParameters *params = iAli->alignmentParameters();
1083     if (params) {
1084       const auto &parVec(params->parameters());
1085       const auto &parCov(params->covariance());
1086       for (int i = 0; i < parVec.num_row(); ++i) {
1087         if (parVec[i] != 0.)
1088           return false;
1089         for (int j = i; j < parCov.num_col(); ++j) {
1090           if (parCov[i][j] != 0.)
1091             return false;
1092         }
1093       }
1094     }
1095   }
1096
1097   return true;
1098 }
1099
1100 //__________________________________________________________________________________________________
1101 unsigned int MillePedeAlignmentAlgorithm::doIO(int loop) const {
1102   unsigned int result = 0;
1103
1104   const std::string outFilePlain(theConfig.getParameter<std::string>("treeFile"));
1105   if (outFilePlain.empty()) {
1106     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1107                               << "treeFile parameter empty => skip writing for 'loop' " << loop;
1108     return result;
1109   }
1110
1111   const std::string outFile(theDir + outFilePlain);
1112
1113   AlignmentIORoot aliIO;
1114   int ioerr = 0;
1115   if (loop == 0) {
1116     aliIO.writeAlignableOriginalPositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1117     if (ioerr) {
1118       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1119                                  << "Problem " << ioerr << " in writeAlignableOriginalPositions";
1120       ++result;
1121     }
1122   } else if (loop == 1) {
1123     // only for first iov add hit counts, else 2x, 3x,... number of hits in IOV 2, 3,...
1124     const std::vector<std::string> inFiles(theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles"));
1125     const std::vector<std::string> binFiles(theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
1126     if (inFiles.size() != binFiles.size()) {
1127       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1128                                    << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' "
1129                                    << "differ in size";
1130     }
1131     this->addHitStatistics(0, outFile, inFiles);  // add hit info from tree 0 in 'infiles'
1132   }
1133   MillePedeVariablesIORoot millePedeIO;
1134   millePedeIO.writeMillePedeVariables(theAlignables, outFile.c_str(), loop, false, ioerr);
1135   if (ioerr) {
1136     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1137                                << "Problem " << ioerr << " writing MillePedeVariables";
1138     ++result;
1139   }
1140
1141   aliIO.writeOrigRigidBodyAlignmentParameters(theAlignables, outFile.c_str(), loop, false, ioerr);
1142   if (ioerr) {
1143     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1144                                << "Problem " << ioerr << " in writeOrigRigidBodyAlignmentParameters, " << loop;
1145     ++result;
1146   }
1147   aliIO.writeAlignableAbsolutePositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1148   if (ioerr) {
1149     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1150                                << "Problem " << ioerr << " in writeAlignableAbsolutePositions, " << loop;
1151     ++result;
1152   }
1153
1154   return result;
1155 }
1156
1157 //__________________________________________________________________________________________________
1158 void MillePedeAlignmentAlgorithm::buildUserVariables(const align::Alignables &alis) const {
1159   for (const auto &iAli : alis) {
1160     AlignmentParameters *params = iAli->alignmentParameters();
1161     if (!params) {
1162       throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::buildUserVariables"
1163                                         << "No parameters for alignable";
1164     }
1165     MillePedeVariables *userVars = dynamic_cast<MillePedeVariables *>(params->userVariables());
1166     if (userVars) {  // Just re-use existing, keeping label and numHits:
1167       for (unsigned int iPar = 0; iPar < userVars->size(); ++iPar) {
1168         //      if (params->hierarchyLevel() > 0) {
1169         //std::cout << params->hierarchyLevel() << "\nBefore: " << userVars->parameter()[iPar];
1170         //      }
1171         userVars->setAllDefault(iPar);
1172         //std::cout << "\nAfter: " << userVars->parameter()[iPar] << std::endl;
1173       }
1174     } else {  // Nothing yet or erase wrong type:
1175       userVars = new MillePedeVariables(
1176           params->size(),
1177           thePedeLabels->alignableLabel(iAli),
1178           thePedeLabels->alignableTracker()->objectIdProvider().typeToName(iAli->alignableObjectId()));
1179       params->setUserVariables(userVars);
1180     }
1181   }
1182 }
1183
1184 //__________________________________________________________________________________________________
1185 unsigned int MillePedeAlignmentAlgorithm::decodeMode(const std::string &mode) const {
1186   if (mode == "full") {
1187     return myMilleBit + myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1188   } else if (mode == "mille") {
1189     return myMilleBit;  // + myPedeSteerBit; // sic! Including production of steerig file. NO!
1190   } else if (mode == "pede") {
1191     return myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1192   } else if (mode == "pedeSteer") {
1193     return myPedeSteerBit;
1194   } else if (mode == "pedeRun") {
1195     return myPedeSteerBit + myPedeRunBit + myPedeReadBit;  // sic! Including steering and reading of result.
1196   } else if (mode == "pedeRead") {
1197     return myPedeReadBit;
1198   }
1199
1200   throw cms::Exception("BadConfig") << "Unknown mode '" << mode
1201                                     << "', use 'full', 'mille', 'pede', 'pedeRun', 'pedeSteer' or 'pedeRead'.";
1202
1203   return 0;
1204 }
1205
1206 //__________________________________________________________________________________________________
1207 bool MillePedeAlignmentAlgorithm::addHitStatistics(int fromIov,
1208                                                    const std::string &outFile,
1209                                                    const std::vector<std::string> &inFiles) const {
1210   bool allOk = true;
1211   int ierr = 0;
1212   MillePedeVariablesIORoot millePedeIO;
1213   for (std::vector<std::string>::const_iterator iFile = inFiles.begin(); iFile != inFiles.end(); ++iFile) {
1214     const std::string inFile(theDir + *iFile);
1215     const std::vector<AlignmentUserVariables *> mpVars =
1216         millePedeIO.readMillePedeVariables(theAlignables, inFile.c_str(), fromIov, ierr);
1217     if (ierr || !this->addHits(theAlignables, mpVars)) {
1218       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addHitStatistics"
1219                                  << "Error " << ierr << " reading from " << inFile << ", tree " << fromIov
1220                                  << ", or problems in addHits";
1221       allOk = false;
1222     }
1223     for (std::vector<AlignmentUserVariables *>::const_iterator i = mpVars.begin(); i != mpVars.end(); ++i) {
1224       delete *i;  // clean created objects
1225     }
1226   }
1227
1228   return allOk;
1229 }
1230
1231 //__________________________________________________________________________________________________
1232 bool MillePedeAlignmentAlgorithm::addHits(const align::Alignables &alis,
1233                                           const std::vector<AlignmentUserVariables *> &mpVars) const {
1234   bool allOk = (mpVars.size() == alis.size());
1235   std::vector<AlignmentUserVariables *>::const_iterator iUser = mpVars.begin();
1236   for (auto iAli = alis.cbegin(); iAli != alis.cend() && iUser != mpVars.end(); ++iAli, ++iUser) {
1237     MillePedeVariables *mpVarNew = dynamic_cast<MillePedeVariables *>(*iUser);
1238     AlignmentParameters *ps = (*iAli)->alignmentParameters();
1239     MillePedeVariables *mpVarOld = (ps ? dynamic_cast<MillePedeVariables *>(ps->userVariables()) : nullptr);
1240     if (!mpVarNew || !mpVarOld || mpVarOld->size() != mpVarNew->size()) {
1241       allOk = false;
1242       continue;  // FIXME error etc.?
1243     }
1244
1245     mpVarOld->increaseHitsX(mpVarNew->hitsX());
1246     mpVarOld->increaseHitsY(mpVarNew->hitsY());
1247   }
1248
1249   return allOk;
1250 }
1251
1252 //__________________________________________________________________________________________________
1253 template <typename GlobalDerivativeMatrix>
1254 void MillePedeAlignmentAlgorithm::makeGlobDerivMatrix(const std::vector<float> &globalDerivativesx,
1255                                                       const std::vector<float> &globalDerivativesy,
1256                                                       Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) {
1257   static_assert(GlobalDerivativeMatrix::RowsAtCompileTime == 2, "global derivative matrix must have two rows");
1258
1259   for (size_t i = 0; i < globalDerivativesx.size(); ++i) {
1260     aGlobalDerivativesM(0, i) = globalDerivativesx[i];
1261     aGlobalDerivativesM(1, i) = globalDerivativesy[i];
1262   }
1263 }
1264
1265 //__________________________________________________________________________________________________
1266 template <typename CovarianceMatrix,
1267           typename LocalDerivativeMatrix,
1268           typename ResidualMatrix,
1269           typename GlobalDerivativeMatrix>
1270 void MillePedeAlignmentAlgorithm::diagonalize(Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1271                                               Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM,
1272                                               Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1273                                               Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) const {
1274   static_assert(std::is_same<typename LocalDerivativeMatrix::Scalar, typename ResidualMatrix::Scalar>::value,
1275                 "'aLocalDerivativesM' and 'aHitResidualsM' must have the "
1276                 "same underlying scalar type");
1277   static_assert(std::is_same<typename LocalDerivativeMatrix::Scalar, typename GlobalDerivativeMatrix::Scalar>::value,
1278                 "'aLocalDerivativesM' and 'aGlobalDerivativesM' must have the "
1279                 "same underlying scalar type");
1280
1281   Eigen::SelfAdjointEigenSolver<typename CovarianceMatrix::PlainObject> myDiag{aHitCovarianceM};
1282   // eigenvectors of real symmetric matrices are orthogonal, i.e. invert == transpose
1283   auto aTranfoToDiagonalSystemInv =
1284       myDiag.eigenvectors().transpose().template cast<typename LocalDerivativeMatrix::Scalar>();
1285
1286   aHitCovarianceM = myDiag.eigenvalues().asDiagonal();
1287   aLocalDerivativesM = aTranfoToDiagonalSystemInv * aLocalDerivativesM;
1288   aHitResidualsM = aTranfoToDiagonalSystemInv * aHitResidualsM;
1289   if (aGlobalDerivativesM.size() > 0) {
1290     // diagonalize only if measurement depends on alignables or calibrations
1291     aGlobalDerivativesM = aTranfoToDiagonalSystemInv * aGlobalDerivativesM;
1292   }
1293 }
1294
1295 //__________________________________________________________________________________________________
1296 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1297 void MillePedeAlignmentAlgorithm ::addRefTrackVirtualMeas1D(
1298     const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1299     unsigned int iVirtualMeas,
1300     Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1301     Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1302     Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1303   // This Method is valid for 1D measurements only
1304
1305   const unsigned int xIndex = iVirtualMeas + refTrajPtr->numberOfHitMeas();
1306
1307   aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1308   aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex];
1309
1310   const auto &locDerivMatrix = refTrajPtr->derivatives();
1311   for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1312     aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1313   }
1314 }
1315
1316 //__________________________________________________________________________________________________
1317 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1318 void MillePedeAlignmentAlgorithm ::addRefTrackData2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1319                                                      unsigned int iTrajHit,
1320                                                      Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1321                                                      Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1322                                                      Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1323   // This Method is valid for 2D measurements only
1324
1325   const unsigned int xIndex = iTrajHit * 2;
1326   const unsigned int yIndex = iTrajHit * 2 + 1;
1327
1328   aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1329   aHitCovarianceM(0, 1) = refTrajPtr->measurementErrors()[xIndex][yIndex];
1330   aHitCovarianceM(1, 0) = refTrajPtr->measurementErrors()[yIndex][xIndex];
1331   aHitCovarianceM(1, 1) = refTrajPtr->measurementErrors()[yIndex][yIndex];
1332
1333   aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1334   aHitResidualsM(1, 0) = refTrajPtr->measurements()[yIndex] - refTrajPtr->trajectoryPositions()[yIndex];
1335
1336   const auto &locDerivMatrix = refTrajPtr->derivatives();
1337   for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1338     aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1339     aLocalDerivativesM(1, i) = locDerivMatrix[yIndex][i];
1340   }
1341 }
1342
1343 //__________________________________________________________________________________________________
1344 int MillePedeAlignmentAlgorithm ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1345                                             unsigned int iTrajHit,
1346                                             const std::vector<int> &globalLabels,
1347                                             const std::vector<float> &globalDerivativesX,
1348                                             const std::vector<float> &globalDerivativesY) {
1349   const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1350
1351   if ((aRecHit)->dimension() == 1) {
1352     return this->callMille1D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX);
1353   } else {
1354     return this->callMille2D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX, globalDerivativesY);
1355   }
1356 }
1357
1358 //__________________________________________________________________________________________________
1359 int MillePedeAlignmentAlgorithm ::callMille1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1360                                               unsigned int iTrajHit,
1361                                               const std::vector<int> &globalLabels,
1362                                               const std::vector<float> &globalDerivativesX) {
1363   const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1364   const unsigned int xIndex = iTrajHit * 2;  // the even ones are local x
1365
1366   // local derivatives
1367   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1368   const int nLocal = locDerivMatrix.num_col();
1369   std::vector<float> localDerivatives(nLocal);
1370   for (unsigned int i = 0; i < localDerivatives.size(); ++i) {
1371     localDerivatives[i] = locDerivMatrix[xIndex][i];
1372   }
1373
1374   // residuum and error
1375   float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1376   float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
1377
1378   // number of global derivatives
1379   const int nGlobal = globalDerivativesX.size();
1380
1381   // &(localDerivatives[0]) etc. are valid - as long as vector is not empty
1382   // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1383   theMille->mille(
1384       nLocal, &(localDerivatives[0]), nGlobal, &(globalDerivativesX[0]), &(globalLabels[0]), residX, hitErrX);
1385
1386   if (theMonitor) {
1387     theMonitor->fillDerivatives(
1388         aRecHit, &(localDerivatives[0]), nLocal, &(globalDerivativesX[0]), nGlobal, &(globalLabels[0]));
1389     theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, residX, hitErrX, false);
1390   }
1391
1392   return 1;
1393 }
1394
1395 //__________________________________________________________________________________________________
1396 int MillePedeAlignmentAlgorithm ::callMille2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1397                                               unsigned int iTrajHit,
1398                                               const std::vector<int> &globalLabels,
1399                                               const std::vector<float> &globalDerivativesx,
1400                                               const std::vector<float> &globalDerivativesy) {
1401   const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1402
1403   if ((aRecHit)->dimension() != 2) {
1404     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::callMille2D"
1405                                << "You try to call method for 2D hits for a " << (aRecHit)->dimension()
1406                                << "D Hit. Hit gets ignored!";
1407     return -1;
1408   }
1409
1410   Eigen::Matrix<double, 2, 2> aHitCovarianceM;
1411   Eigen::Matrix<float, 2, 1> aHitResidualsM;
1412   Eigen::Matrix<float, 2, Eigen::Dynamic> aLocalDerivativesM{2, refTrajPtr->derivatives().num_col()};
1413   // below method fills above 3 matrices
1414   this->addRefTrackData2D(refTrajPtr, iTrajHit, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1415   Eigen::Matrix<float, 2, Eigen::Dynamic> aGlobalDerivativesM{2, globalDerivativesx.size()};
1416   this->makeGlobDerivMatrix(globalDerivativesx, globalDerivativesy, aGlobalDerivativesM);
1417
1418   // calculates correlation between Hit measurements
1419   // FIXME: Should take correlation (and resulting transformation) from original hit,
1420   //        not 2x2 matrix from ReferenceTrajectory: That can come from error propagation etc.!
1421   const double corr = aHitCovarianceM(0, 1) / sqrt(aHitCovarianceM(0, 0) * aHitCovarianceM(1, 1));
1422   if (theMonitor)
1423     theMonitor->fillCorrelations2D(corr, aRecHit);
1424   bool diag = false;  // diagonalise only tracker TID, TEC
1425   switch (aRecHit->geographicalId().subdetId()) {
1426     case SiStripDetId::TID:
1427     case SiStripDetId::TEC:
1428       if (aRecHit->geographicalId().det() == DetId::Tracker && TMath::Abs(corr) > theMaximalCor2D) {
1429         this->diagonalize(aHitCovarianceM, aLocalDerivativesM, aHitResidualsM, aGlobalDerivativesM);
1430         diag = true;
1431       }
1432       break;
1433     default:;
1434   }
1435
1436   float newResidX = aHitResidualsM(0, 0);
1437   float newResidY = aHitResidualsM(1, 0);
1438   float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1439   float newHitErrY = TMath::Sqrt(aHitCovarianceM(1, 1));
1440
1441   // change from column major (Eigen default) to row major to have row entries
1442   // in continuous memory
1443   std::vector<float> newLocalDerivs(aLocalDerivativesM.size());
1444   Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1445       newLocalDerivs.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1446   float *newLocalDerivsX = &(newLocalDerivs[0]);
1447   float *newLocalDerivsY = &(newLocalDerivs[aLocalDerivativesM.cols()]);
1448
1449   // change from column major (Eigen default) to row major to have row entries
1450   // in continuous memory
1451   std::vector<float> newGlobDerivs(aGlobalDerivativesM.size());
1452   Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1453       newGlobDerivs.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1454   float *newGlobDerivsX = &(newGlobDerivs[0]);
1455   float *newGlobDerivsY = &(newGlobDerivs[aGlobalDerivativesM.cols()]);
1456
1457   const int nLocal = aLocalDerivativesM.cols();
1458   const int nGlobal = aGlobalDerivativesM.cols();
1459
1460   if (diag && (newHitErrX > newHitErrY)) {  // also for 2D hits?
1461     // measurement with smaller error is x-measurement (for !is2D do not fill y-measurement):
1462     std::swap(newResidX, newResidY);
1463     std::swap(newHitErrX, newHitErrY);
1464     std::swap(newLocalDerivsX, newLocalDerivsY);
1465     std::swap(newGlobDerivsX, newGlobDerivsY);
1466   }
1467
1468   // &(globalLabels[0]) is valid - as long as vector is not empty
1469   // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1470   theMille->mille(nLocal, newLocalDerivsX, nGlobal, newGlobDerivsX, &(globalLabels[0]), newResidX, newHitErrX);
1471
1472   if (theMonitor) {
1473     theMonitor->fillDerivatives(aRecHit, newLocalDerivsX, nLocal, newGlobDerivsX, nGlobal, &(globalLabels[0]));
1474     theMonitor->fillResiduals(
1475         aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidX, newHitErrX, false);
1476   }
1477   const bool isReal2DHit = this->is2D(aRecHit);  // strip is 1D (except matched hits)
1478   if (isReal2DHit) {
1479     theMille->mille(nLocal, newLocalDerivsY, nGlobal, newGlobDerivsY, &(globalLabels[0]), newResidY, newHitErrY);
1480     if (theMonitor) {
1481       theMonitor->fillDerivatives(aRecHit, newLocalDerivsY, nLocal, newGlobDerivsY, nGlobal, &(globalLabels[0]));
1482       theMonitor->fillResiduals(
1483           aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidY, newHitErrY, true);  // true: y
1484     }
1485   }
1486
1487   return (isReal2DHit ? 2 : 1);
1488 }
1489
1490 //__________________________________________________________________________________________________
1491 void MillePedeAlignmentAlgorithm ::addVirtualMeas(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1492                                                   unsigned int iVirtualMeas) {
1493   Eigen::Matrix<double, 1, 1> aHitCovarianceM;
1494   Eigen::Matrix<float, 1, 1> aHitResidualsM;
1495   Eigen::Matrix<float, 1, Eigen::Dynamic> aLocalDerivativesM{1, refTrajPtr->derivatives().num_col()};
1496   // below method fills above 3 'matrices'
1497   this->addRefTrackVirtualMeas1D(refTrajPtr, iVirtualMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1498
1499   // no global parameters (use dummy 0)
1500   auto aGlobalDerivativesM = Eigen::Matrix<float, 1, 1>::Zero();
1501
1502   float newResidX = aHitResidualsM(0, 0);
1503   float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1504   std::vector<float> newLocalDerivsX(aLocalDerivativesM.size());
1505   Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1506       newLocalDerivsX.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1507
1508   std::vector<float> newGlobDerivsX(aGlobalDerivativesM.size());
1509   Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1510       newGlobDerivsX.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1511
1512   const int nLocal = aLocalDerivativesM.cols();
1513   const int nGlobal = 0;
1514
1515   theMille->mille(nLocal, newLocalDerivsX.data(), nGlobal, newGlobDerivsX.data(), &nGlobal, newResidX, newHitErrX);
1516 }
1517
1518 //____________________________________________________
1519 void MillePedeAlignmentAlgorithm::addLaserData(const EventInfo &eventInfo,
1520                                                const TkFittedLasBeamCollection &lasBeams,
1521                                                const TsosVectorCollection &lasBeamTsoses) {
1522   TsosVectorCollection::const_iterator iTsoses = lasBeamTsoses.begin();
1523   for (TkFittedLasBeamCollection::const_iterator iBeam = lasBeams.begin(), iEnd = lasBeams.end(); iBeam != iEnd;
1524        ++iBeam, ++iTsoses) {  // beam/tsoses parallel!
1525
1526     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addLaserData"
1527                               << "Beam " << iBeam->getBeamId() << " with " << iBeam->parameters().size()
1528                               << " parameters and " << iBeam->getData().size() << " hits.\n There are "
1529                               << iTsoses->size() << " TSOSes.";
1530
1531     this->addLasBeam(eventInfo, *iBeam, *iTsoses);
1532   }
1533 }
1534
1535 //____________________________________________________
1536 void MillePedeAlignmentAlgorithm::addLasBeam(const EventInfo &eventInfo,
1537                                              const TkFittedLasBeam &lasBeam,
1538                                              const std::vector<TrajectoryStateOnSurface> &tsoses) {
1539   AlignmentParameters *dummyPtr = nullptr;                                          // for globalDerivativesHierarchy()
1540   std::vector<float> lasLocalDerivsX;                                               // buffer for local derivatives
1541   const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId());  // for global par
1542
1543   for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
1544     if (!tsoses[iHit].isValid())
1545       continue;
1546     // clear buffer
1547     theFloatBufferX.clear();
1548     theFloatBufferY.clear();
1549     theIntBuffer.clear();
1550     lasLocalDerivsX.clear();
1551     // get alignables and global parameters
1552     const SiStripLaserRecHit2D &hit = lasBeam.getData()[iHit];
1553     AlignableDetOrUnitPtr lasAli(theAlignableNavigator->alignableFromDetId(hit.getDetId()));
1554     this->globalDerivativesHierarchy(
1555         eventInfo, tsoses[iHit], lasAli, lasAli, theFloatBufferX, theFloatBufferY, theIntBuffer, dummyPtr);
1556     // fill derivatives vector from derivatives matrix
1557     for (unsigned int nFitParams = 0; nFitParams < static_cast<unsigned int>(lasBeam.parameters().size());
1558          ++nFitParams) {
1559       const float derivative = lasBeam.derivatives()[iHit][nFitParams];
1560       if (nFitParams < lasBeam.firstFixedParameter()) {  // first local beam parameters
1561         lasLocalDerivsX.push_back(derivative);
1562       } else {  // now global ones
1563         const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
1564         theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
1565         theFloatBufferX.push_back(derivative);
1566       }
1567     }  // end loop over parameters
1568
1569     const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
1570     // error from file or assume 0.003
1571     const float error = 0.003;  // hit.localPositionError().xx(); sqrt???
1572
1573     theMille->mille(lasLocalDerivsX.size(),
1574                     &(lasLocalDerivsX[0]),
1575                     theFloatBufferX.size(),
1576                     &(theFloatBufferX[0]),
1577                     &(theIntBuffer[0]),
1578                     residual,
1579                     error);
1580   }  // end of loop over hits
1581
1582   theMille->end();
1583 }
1584
1585 void MillePedeAlignmentAlgorithm::addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg) {
1586   // do some printing, if requested
1587   const bool doOutputOnStdout(pxbSurveyCfg.getParameter<bool>("doOutputOnStdout"));
1588   if (doOutputOnStdout) {
1589     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1590                               << "# Output from addPxbSurvey follows below because "
1591                               << "doOutputOnStdout is set to True";
1592   }
1593
1594   // instantiate a dicer object
1595   SurveyPxbDicer dicer(pxbSurveyCfg.getParameter<std::vector<edm::ParameterSet> >("toySurveyParameters"),
1596                        pxbSurveyCfg.getParameter<unsigned int>("toySurveySeed"));
1597   std::ofstream outfile(pxbSurveyCfg.getUntrackedParameter<std::string>("toySurveyFile").c_str());
1598
1599   // read data from file
1600   std::vector<SurveyPxbImageLocalFit> measurements;
1601   std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
1602   SurveyPxbImageReader<SurveyPxbImageLocalFit> reader(filename, measurements, 800);
1603
1604   // loop over photographs (=measurements) and perform the fit
1605   for (std::vector<SurveyPxbImageLocalFit>::size_type i = 0; i != measurements.size(); i++) {
1606     if (doOutputOnStdout) {
1607       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1608                                 << "Module " << i << ": ";
1609     }
1610
1611     // get the Alignables and their surfaces
1612     AlignableDetOrUnitPtr mod1(theAlignableNavigator->alignableFromDetId(measurements[i].getIdFirst()));
1613     AlignableDetOrUnitPtr mod2(theAlignableNavigator->alignableFromDetId(measurements[i].getIdSecond()));
1614     const AlignableSurface &surf1 = mod1->surface();
1615     const AlignableSurface &surf2 = mod2->surface();
1616
1617     // the position of the fiducial points in local frame of a PXB module
1618     const LocalPoint fidpoint0(-0.91, +3.30);
1619     const LocalPoint fidpoint1(+0.91, +3.30);
1620     const LocalPoint fidpoint2(+0.91, -3.30);
1621     const LocalPoint fidpoint3(-0.91, -3.30);
1622
1623     // We choose the local frame of the first module as reference,
1624     // so take the fidpoints of the second module and calculate their
1625     // positions in the reference frame
1626     const GlobalPoint surf2point0(surf2.toGlobal(fidpoint0));
1627     const GlobalPoint surf2point1(surf2.toGlobal(fidpoint1));
1628     const LocalPoint fidpoint0inSurf1frame(surf1.toLocal(surf2point0));
1629     const LocalPoint fidpoint1inSurf1frame(surf1.toLocal(surf2point1));
1630
1631     // Create the vector for the fit
1632     SurveyPxbImageLocalFit::fidpoint_t fidpointvec;
1633     fidpointvec.push_back(fidpoint0inSurf1frame);
1634     fidpointvec.push_back(fidpoint1inSurf1frame);
1635     fidpointvec.push_back(fidpoint2);
1636     fidpointvec.push_back(fidpoint3);
1637
1638     // if toy survey is requested, dice the values now
1639     if (pxbSurveyCfg.getParameter<bool>("doToySurvey")) {
1640       dicer.doDice(fidpointvec, measurements[i].getIdPair(), outfile);
1641     }
1642
1643     // do the fit
1644     measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
1645     SurveyPxbImageLocalFit::localpars_t a;  // local pars from fit
1646     a = measurements[i].getLocalParameters();
1647     const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
1648
1649     // do some reporting, if requested
1650     if (doOutputOnStdout) {
1651       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1652                                 << "a: " << a[0] << ", " << a[1] << ", " << a[2] << ", " << a[3]
1653                                 << " S= " << sqrt(a[2] * a[2] + a[3] * a[3]) << " phi= " << atan(a[3] / a[2])
1654                                 << " chi2= " << chi2 << std::endl;
1655     }
1656     if (theMonitor) {
1657       theMonitor->fillPxbSurveyHistsChi2(chi2);
1658       theMonitor->fillPxbSurveyHistsLocalPars(a[0], a[1], sqrt(a[2] * a[2] + a[3] * a[3]), atan(a[3] / a[2]));
1659     }
1660
1661     // pass the results from the local fit to mille
1662     for (SurveyPxbImageLocalFit::count_t j = 0; j != SurveyPxbImageLocalFit::nMsrmts; j++) {
1663       theMille->mille((int)measurements[i].getLocalDerivsSize(),
1664                       measurements[i].getLocalDerivsPtr(j),
1665                       (int)measurements[i].getGlobalDerivsSize(),
1666                       measurements[i].getGlobalDerivsPtr(j),
1667                       measurements[i].getGlobalDerivsLabelPtr(j),
1668                       measurements[i].getResiduum(j),
1669                       measurements[i].getSigma(j));
1670     }
1671     theMille->end();
1672   }
1673   outfile.close();
1674 }
1675
1676 bool MillePedeAlignmentAlgorithm::areIOVsSpecified() const {
1677   const auto runRangeSelection = theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection");
1678
1679   if (runRangeSelection.empty())
1680     return false;
1681
1682   const auto runRanges =
1683       align::makeNonOverlappingRunRanges(runRangeSelection, cond::timeTypeSpecs[cond::runnumber].beginValue);
1684
1685   return !(runRanges.empty());
1686 }