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