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     theBinary.reset();  // GBL output has to be considered since same binary file is used
0728     theBinary = std::make_unique<MilleBinary>((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str(),
0729                                               theGblDoubleBinary);
0730   }
0731 }
0732
0733 //____________________________________________________
0734 void MillePedeAlignmentAlgorithm::endLuminosityBlock(const edm::EventSetup &) {
0735   if (!runAtPCL_)
0736     return;
0737   if (this->isMode(myMilleBit))
0738     theMille->flushOutputFile();
0739 }
0740
0741 //____________________________________________________
0742 int MillePedeAlignmentAlgorithm::addMeasurementData(const edm::EventSetup &setup,
0743                                                     const EventInfo &eventInfo,
0744                                                     const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0745                                                     unsigned int iHit,
0746                                                     AlignmentParameters *&params) {
0747   params = nullptr;
0748   theFloatBufferX.clear();
0749   theFloatBufferY.clear();
0750   theIntBuffer.clear();
0751
0752   const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
0753   const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
0754   // ignore invalid hits
0755   if (!recHitPtr->isValid())
0756     return 0;
0757
0758   // First add the derivatives from IntegratedCalibration's,
0759   // should even be OK if problems for "usual" derivatives from Alignables
0760   this->globalDerivativesCalibration(recHitPtr,
0761                                      tsos,
0762                                      setup,
0763                                      eventInfo,  // input
0764                                      theFloatBufferX,
0765                                      theFloatBufferY,
0766                                      theIntBuffer);  // output
0767
0768   // get AlignableDet/Unit for this hit
0769   AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0770
0771   if (!this->globalDerivativesHierarchy(eventInfo,
0772                                         tsos,
0773                                         alidet,
0774                                         alidet,
0775                                         theFloatBufferX,  // 2x alidet, sic!
0776                                         theFloatBufferY,
0777                                         theIntBuffer,
0778                                         params)) {
0779     return -1;  // problem
0780   } else if (theFloatBufferX.empty() && ignoreHitsWithoutGlobalDerivatives_) {
0781     return 0;  // empty for X: no alignable for hit, nor calibrations
0782   } else {
0783     // store measurement even if no alignable or calibrations
0784     // -> measurement used for pede-internal track-fit
0785     return this->callMille(refTrajPtr, iHit, theIntBuffer, theFloatBufferX, theFloatBufferY);
0786   }
0787 }
0788
0789 //____________________________________________________
0790
0791 int MillePedeAlignmentAlgorithm::addGlobalData(const edm::EventSetup &setup,
0792                                                const EventInfo &eventInfo,
0793                                                const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0794                                                unsigned int iHit,
0795                                                GblPoint &gblPoint) {
0796   AlignmentParameters *params = nullptr;
0797   std::vector<double> theDoubleBufferX, theDoubleBufferY;
0798   theDoubleBufferX.clear();
0799   theDoubleBufferY.clear();
0800   theIntBuffer.clear();
0801   int iret = 0;
0802
0803   const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
0804   const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
0805   // ignore invalid hits
0806   if (!recHitPtr->isValid())
0807     return 0;
0808
0809   // get AlignableDet/Unit for this hit
0810   AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0811
0812   if (!this->globalDerivativesHierarchy(eventInfo,
0813                                         tsos,
0814                                         alidet,
0815                                         alidet,
0816                                         theDoubleBufferX,  // 2x alidet, sic!
0817                                         theDoubleBufferY,
0818                                         theIntBuffer,
0819                                         params)) {
0820     return -1;  // problem
0821   }
0822   //calibration parameters
0823   int globalLabel;
0824   std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0825   for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0826     // get all derivatives of this calibration // const unsigned int num =
0827     (*iCalib)->derivatives(derivs, *recHitPtr, tsos, setup, eventInfo);
0828     for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0829       // transfer label and x/y derivatives
0830       globalLabel = thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second);
0831       if (globalLabel > 0 && globalLabel <= 2147483647) {
0832         theIntBuffer.push_back(globalLabel);
0833         theDoubleBufferX.push_back(iValuesInd->first.first);
0834         theDoubleBufferY.push_back(iValuesInd->first.second);
0835       } else {
0836         edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addGlobalData"
0837                                    << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
0838       }
0839     }
0840   }
0841   unsigned int numGlobals = theIntBuffer.size();
0842   if (numGlobals > 0) {
0843     Eigen::Matrix<double, 2, Eigen::Dynamic> globalDer{2, numGlobals};
0844     for (unsigned int i = 0; i < numGlobals; ++i) {
0845       globalDer(0, i) = theDoubleBufferX[i];
0846       globalDer(1, i) = theDoubleBufferY[i];
0847     }
0848     gblPoint.addGlobals(theIntBuffer, globalDer);
0849     iret = 1;
0850   }
0851   return iret;
0852 }
0853
0854 //____________________________________________________
0855 bool MillePedeAlignmentAlgorithm ::globalDerivativesHierarchy(const EventInfo &eventInfo,
0856                                                               const TrajectoryStateOnSurface &tsos,
0857                                                               Alignable *ali,
0858                                                               const AlignableDetOrUnitPtr &alidet,
0859                                                               std::vector<float> &globalDerivativesX,
0860                                                               std::vector<float> &globalDerivativesY,
0861                                                               std::vector<int> &globalLabels,
0862                                                               AlignmentParameters *&lowestParams) const {
0863   // derivatives and labels are recursively attached
0864   if (!ali)
0865     return true;  // no mother might be OK
0866
0867   if (false && theMonitor && alidet != ali)
0868     theMonitor->fillFrameToFrame(alidet, ali);
0869
0870   AlignmentParameters *params = ali->alignmentParameters();
0871
0872   if (params) {
0873     if (!lowestParams)
0874       lowestParams = params;  // set parameters of lowest level
0875
0876     bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0877     const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0878
0879     if (0 == alignableLabel) {  // FIXME: what about regardAllHits in Markus' code?
0880       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0881                                    << "Label not found, skip Alignable.";
0882       return false;
0883     }
0884
0885     const std::vector<bool> &selPars = params->selector();
0886     const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
0887
0888     // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
0889     for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
0890       if (selPars[iSel]) {
0891         globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
0892         if (hasSplitParameters == true) {
0893           globalLabels.push_back(thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos));
0894         } else {
0895           globalLabels.push_back(thePedeLabels->parameterLabel(alignableLabel, iSel));
0896         }
0897         globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
0898       }
0899     }
0900     // Exclude mothers if Alignable selected to be no part of a hierarchy:
0901     if (thePedeSteer->isNoHiera(ali))
0902       return true;
0903   }
0904   // Call recursively for mother, will stop if mother == 0:
0905   return this->globalDerivativesHierarchy(
0906       eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
0907 }
0908
0909 //____________________________________________________
0910 bool MillePedeAlignmentAlgorithm ::globalDerivativesHierarchy(const EventInfo &eventInfo,
0911                                                               const TrajectoryStateOnSurface &tsos,
0912                                                               Alignable *ali,
0913                                                               const AlignableDetOrUnitPtr &alidet,
0914                                                               std::vector<double> &globalDerivativesX,
0915                                                               std::vector<double> &globalDerivativesY,
0916                                                               std::vector<int> &globalLabels,
0917                                                               AlignmentParameters *&lowestParams) const {
0918   // derivatives and labels are recursively attached
0919   if (!ali)
0920     return true;  // no mother might be OK
0921
0922   if (false && theMonitor && alidet != ali)
0923     theMonitor->fillFrameToFrame(alidet, ali);
0924
0925   AlignmentParameters *params = ali->alignmentParameters();
0926
0927   if (params) {
0928     if (!lowestParams)
0929       lowestParams = params;  // set parameters of lowest level
0930
0931     bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0932     const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0933
0934     if (0 == alignableLabel) {  // FIXME: what about regardAllHits in Markus' code?
0935       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0936                                    << "Label not found, skip Alignable.";
0937       return false;
0938     }
0939
0940     const std::vector<bool> &selPars = params->selector();
0941     const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
0942     int globalLabel;
0943
0944     // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
0945     for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
0946       if (selPars[iSel]) {
0947         if (hasSplitParameters == true) {
0948           globalLabel = thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos);
0949         } else {
0950           globalLabel = thePedeLabels->parameterLabel(alignableLabel, iSel);
0951         }
0952         if (globalLabel > 0 && globalLabel <= 2147483647) {
0953           globalLabels.push_back(globalLabel);
0954           globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
0955           globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
0956         } else {
0957           edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0958                                      << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
0959         }
0960       }
0961     }
0962     // Exclude mothers if Alignable selected to be no part of a hierarchy:
0963     if (thePedeSteer->isNoHiera(ali))
0964       return true;
0965   }
0966   // Call recursively for mother, will stop if mother == 0:
0967   return this->globalDerivativesHierarchy(
0968       eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
0969 }
0970
0971 //____________________________________________________
0972 void MillePedeAlignmentAlgorithm::globalDerivativesCalibration(const TransientTrackingRecHit::ConstRecHitPointer &recHit,
0973                                                                const TrajectoryStateOnSurface &tsos,
0974                                                                const edm::EventSetup &setup,
0975                                                                const EventInfo &eventInfo,
0976                                                                std::vector<float> &globalDerivativesX,
0977                                                                std::vector<float> &globalDerivativesY,
0978                                                                std::vector<int> &globalLabels) const {
0979   std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0980   for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0981     // get all derivatives of this calibration // const unsigned int num =
0982     (*iCalib)->derivatives(derivs, *recHit, tsos, setup, eventInfo);
0983     for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0984       // transfer label and x/y derivatives
0985       globalLabels.push_back(thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second));
0986       globalDerivativesX.push_back(iValuesInd->first.first);
0987       globalDerivativesY.push_back(iValuesInd->first.second);
0988     }
0989   }
0990 }
0991
0992 // //____________________________________________________
0993 // void MillePedeAlignmentAlgorithm
0994 // ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0995 //             unsigned int iTrajHit, MeasurementDirection xOrY,
0996 //             const std::vector<float> &globalDerivatives, const std::vector<int> &globalLabels)
0997 // {
0998 //   const unsigned int xyIndex = iTrajHit*2 + xOrY;
0999 //   // FIXME: here for residuum and sigma we could use KALMAN-Filter results
1000 //   const float residuum =
1001 //     refTrajPtr->measurements()[xyIndex] - refTrajPtr->trajectoryPositions()[xyIndex];
1002 //   const float covariance = refTrajPtr->measurementErrors()[xyIndex][xyIndex];
1003 //   const float sigma = (covariance > 0. ? TMath::Sqrt(covariance) : 0.);
1004
1005 //   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1006
1007 //   std::vector<float> localDerivs(locDerivMatrix.num_col());
1008 //   for (unsigned int i = 0; i < localDerivs.size(); ++i) {
1009 //     localDerivs[i] = locDerivMatrix[xyIndex][i];
1010 //   }
1011
1012 //   // &(vector[0]) is valid - as long as vector is not empty
1013 //   // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1014 //   theMille->mille(localDerivs.size(), &(localDerivs[0]),
1015 //                globalDerivatives.size(), &(globalDerivatives[0]), &(globalLabels[0]),
1016 //                residuum, sigma);
1017 //   if (theMonitor) {
1018 //     theMonitor->fillDerivatives(refTrajPtr->recHits()[iTrajHit],localDerivs, globalDerivatives,
1019 //                              (xOrY == kLocalY));
1020 //     theMonitor->fillResiduals(refTrajPtr->recHits()[iTrajHit],
1021 //                            refTrajPtr->trajectoryStates()[iTrajHit],
1022 //                            iTrajHit, residuum, sigma, (xOrY == kLocalY));
1023 //   }
1024 // }
1025
1026 //____________________________________________________
1027 bool MillePedeAlignmentAlgorithm::is2D(const ConstRecHitPointer &recHit) const {
1028   // FIXME: Check whether this is a reliable and recommended way to find out...
1029
1030   if (recHit->dimension() < 2) {
1031     return false;                  // some muon and TIB/TOB stuff really has RecHit1D
1032   } else if (recHit->detUnit()) {  // detunit in strip is 1D, in pixel 2D
1033     return recHit->detUnit()->type().isTrackerPixel();
1034   } else {  // stereo strips  (FIXME: endcap trouble due to non-parallel strips (wedge sensors)?)
1035     if (dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit->hit())) {  // check persistent hit
1036       // projected: 1D measurement on 'glued' module
1037       return false;
1038     } else {
1039       return true;
1040     }
1041   }
1042 }
1043
1044 //__________________________________________________________________________________________________
1045 bool MillePedeAlignmentAlgorithm::readFromPede(const edm::ParameterSet &mprespset,
1046                                                bool setUserVars,
1047                                                const RunRange &runrange) {
1048   bool allEmpty = this->areEmptyParams(theAlignables);
1049
1050   PedeReader reader(mprespset, *thePedeSteer, *thePedeLabels, runrange);
1051   align::Alignables alis;
1052   bool okRead = reader.read(alis, setUserVars);  // also may set params of IntegratedCalibration's
1053   bool numMatch = true;
1054
1055   std::stringstream out;
1056   out << "Read " << alis.size() << " alignables";
1057   if (alis.size() != theAlignables.size()) {
1058     out << " while " << theAlignables.size() << " in store";
1059     numMatch = false;  // FIXME: Should we check one by one? Or transfer 'alis' to the store?
1060   }
1061   if (!okRead)
1062     out << ", but problems in reading";
1063   if (!allEmpty)
1064     out << ", possibly overwriting previous settings";
1065   out << ".";
1066
1067   if (okRead && allEmpty) {
1068     if (numMatch) {  // as many alignables with result as trying to align
1069       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1070     } else if (!alis.empty()) {  // dead module do not get hits and no pede result
1071       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1072     } else {  // serious problem: no result read - and not all modules can be dead...
1073       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1074       return false;
1075     }
1076     return true;
1077   }
1078   // the rest is not OK:
1079   edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1080   return false;
1081 }
1082
1083 //__________________________________________________________________________________________________
1084 bool MillePedeAlignmentAlgorithm::areEmptyParams(const align::Alignables &alignables) const {
1085   for (const auto &iAli : alignables) {
1086     const AlignmentParameters *params = iAli->alignmentParameters();
1087     if (params) {
1088       const auto &parVec(params->parameters());
1089       const auto &parCov(params->covariance());
1090       for (int i = 0; i < parVec.num_row(); ++i) {
1091         if (parVec[i] != 0.)
1092           return false;
1093         for (int j = i; j < parCov.num_col(); ++j) {
1094           if (parCov[i][j] != 0.)
1095             return false;
1096         }
1097       }
1098     }
1099   }
1100
1101   return true;
1102 }
1103
1104 //__________________________________________________________________________________________________
1105 unsigned int MillePedeAlignmentAlgorithm::doIO(int loop) const {
1106   unsigned int result = 0;
1107
1108   const std::string outFilePlain(theConfig.getParameter<std::string>("treeFile"));
1109   if (outFilePlain.empty()) {
1110     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1111                               << "treeFile parameter empty => skip writing for 'loop' " << loop;
1112     return result;
1113   }
1114
1115   const std::string outFile(theDir + outFilePlain);
1116
1117   AlignmentIORoot aliIO;
1118   int ioerr = 0;
1119   if (loop == 0) {
1120     aliIO.writeAlignableOriginalPositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1121     if (ioerr) {
1122       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1123                                  << "Problem " << ioerr << " in writeAlignableOriginalPositions";
1124       ++result;
1125     }
1126   } else if (loop == 1) {
1127     // only for first iov add hit counts, else 2x, 3x,... number of hits in IOV 2, 3,...
1128     const std::vector<std::string> inFiles(theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles"));
1129     const std::vector<std::string> binFiles(theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
1130     if (inFiles.size() != binFiles.size()) {
1131       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1132                                    << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' "
1133                                    << "differ in size";
1134     }
1135     this->addHitStatistics(0, outFile, inFiles);  // add hit info from tree 0 in 'infiles'
1136   }
1137   MillePedeVariablesIORoot millePedeIO;
1138   millePedeIO.writeMillePedeVariables(theAlignables, outFile.c_str(), loop, false, ioerr);
1139   if (ioerr) {
1140     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1141                                << "Problem " << ioerr << " writing MillePedeVariables";
1142     ++result;
1143   }
1144
1145   aliIO.writeOrigRigidBodyAlignmentParameters(theAlignables, outFile.c_str(), loop, false, ioerr);
1146   if (ioerr) {
1147     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1148                                << "Problem " << ioerr << " in writeOrigRigidBodyAlignmentParameters, " << loop;
1149     ++result;
1150   }
1151   aliIO.writeAlignableAbsolutePositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1152   if (ioerr) {
1153     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1154                                << "Problem " << ioerr << " in writeAlignableAbsolutePositions, " << loop;
1155     ++result;
1156   }
1157
1158   return result;
1159 }
1160
1161 //__________________________________________________________________________________________________
1162 void MillePedeAlignmentAlgorithm::buildUserVariables(const align::Alignables &alis) const {
1163   for (const auto &iAli : alis) {
1164     AlignmentParameters *params = iAli->alignmentParameters();
1165     if (!params) {
1166       throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::buildUserVariables"
1167                                         << "No parameters for alignable";
1168     }
1169     MillePedeVariables *userVars = dynamic_cast<MillePedeVariables *>(params->userVariables());
1170     if (userVars) {  // Just re-use existing, keeping label and numHits:
1171       for (unsigned int iPar = 0; iPar < userVars->size(); ++iPar) {
1172         //      if (params->hierarchyLevel() > 0) {
1173         //std::cout << params->hierarchyLevel() << "\nBefore: " << userVars->parameter()[iPar];
1174         //      }
1175         userVars->setAllDefault(iPar);
1176         //std::cout << "\nAfter: " << userVars->parameter()[iPar] << std::endl;
1177       }
1178     } else {  // Nothing yet or erase wrong type:
1179       userVars = new MillePedeVariables(
1180           params->size(),
1181           thePedeLabels->alignableLabel(iAli),
1182           thePedeLabels->alignableTracker()->objectIdProvider().typeToName(iAli->alignableObjectId()));
1183       params->setUserVariables(userVars);
1184     }
1185   }
1186 }
1187
1188 //__________________________________________________________________________________________________
1189 unsigned int MillePedeAlignmentAlgorithm::decodeMode(const std::string &mode) const {
1190   if (mode == "full") {
1191     return myMilleBit + myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1192   } else if (mode == "mille") {
1193     return myMilleBit;  // + myPedeSteerBit; // sic! Including production of steerig file. NO!
1194   } else if (mode == "pede") {
1195     return myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1196   } else if (mode == "pedeSteer") {
1197     return myPedeSteerBit;
1198   } else if (mode == "pedeRun") {
1199     return myPedeSteerBit + myPedeRunBit + myPedeReadBit;  // sic! Including steering and reading of result.
1200   } else if (mode == "pedeRead") {
1201     return myPedeReadBit;
1202   }
1203
1204   throw cms::Exception("BadConfig") << "Unknown mode '" << mode
1205                                     << "', use 'full', 'mille', 'pede', 'pedeRun', 'pedeSteer' or 'pedeRead'.";
1206
1207   return 0;
1208 }
1209
1210 //__________________________________________________________________________________________________
1211 bool MillePedeAlignmentAlgorithm::addHitStatistics(int fromIov,
1212                                                    const std::string &outFile,
1213                                                    const std::vector<std::string> &inFiles) const {
1214   bool allOk = true;
1215   int ierr = 0;
1216   MillePedeVariablesIORoot millePedeIO;
1217   for (std::vector<std::string>::const_iterator iFile = inFiles.begin(); iFile != inFiles.end(); ++iFile) {
1218     const std::string inFile(theDir + *iFile);
1219     const std::vector<AlignmentUserVariables *> mpVars =
1220         millePedeIO.readMillePedeVariables(theAlignables, inFile.c_str(), fromIov, ierr);
1221     if (ierr || !this->addHits(theAlignables, mpVars)) {
1222       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addHitStatistics"
1223                                  << "Error " << ierr << " reading from " << inFile << ", tree " << fromIov
1224                                  << ", or problems in addHits";
1225       allOk = false;
1226     }
1227     for (std::vector<AlignmentUserVariables *>::const_iterator i = mpVars.begin(); i != mpVars.end(); ++i) {
1228       delete *i;  // clean created objects
1229     }
1230   }
1231
1232   return allOk;
1233 }
1234
1235 //__________________________________________________________________________________________________
1236 bool MillePedeAlignmentAlgorithm::addHits(const align::Alignables &alis,
1237                                           const std::vector<AlignmentUserVariables *> &mpVars) const {
1238   bool allOk = (mpVars.size() == alis.size());
1239   std::vector<AlignmentUserVariables *>::const_iterator iUser = mpVars.begin();
1240   for (auto iAli = alis.cbegin(); iAli != alis.cend() && iUser != mpVars.end(); ++iAli, ++iUser) {
1241     MillePedeVariables *mpVarNew = dynamic_cast<MillePedeVariables *>(*iUser);
1242     AlignmentParameters *ps = (*iAli)->alignmentParameters();
1243     MillePedeVariables *mpVarOld = (ps ? dynamic_cast<MillePedeVariables *>(ps->userVariables()) : nullptr);
1244     if (!mpVarNew || !mpVarOld || mpVarOld->size() != mpVarNew->size()) {
1245       allOk = false;
1246       continue;  // FIXME error etc.?
1247     }
1248
1249     mpVarOld->increaseHitsX(mpVarNew->hitsX());
1250     mpVarOld->increaseHitsY(mpVarNew->hitsY());
1251   }
1252
1253   return allOk;
1254 }
1255
1256 //__________________________________________________________________________________________________
1257 template <typename GlobalDerivativeMatrix>
1258 void MillePedeAlignmentAlgorithm::makeGlobDerivMatrix(const std::vector<float> &globalDerivativesx,
1259                                                       const std::vector<float> &globalDerivativesy,
1260                                                       Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) {
1261   static_assert(GlobalDerivativeMatrix::RowsAtCompileTime == 2, "global derivative matrix must have two rows");
1262
1263   for (size_t i = 0; i < globalDerivativesx.size(); ++i) {
1264     aGlobalDerivativesM(0, i) = globalDerivativesx[i];
1265     aGlobalDerivativesM(1, i) = globalDerivativesy[i];
1266   }
1267 }
1268
1269 //__________________________________________________________________________________________________
1270 template <typename CovarianceMatrix,
1271           typename LocalDerivativeMatrix,
1272           typename ResidualMatrix,
1273           typename GlobalDerivativeMatrix>
1274 void MillePedeAlignmentAlgorithm::diagonalize(Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1275                                               Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM,
1276                                               Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1277                                               Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) const {
1278   static_assert(std::is_same<typename LocalDerivativeMatrix::Scalar, typename ResidualMatrix::Scalar>::value,
1279                 "'aLocalDerivativesM' and 'aHitResidualsM' must have the "
1280                 "same underlying scalar type");
1281   static_assert(std::is_same<typename LocalDerivativeMatrix::Scalar, typename GlobalDerivativeMatrix::Scalar>::value,
1282                 "'aLocalDerivativesM' and 'aGlobalDerivativesM' must have the "
1283                 "same underlying scalar type");
1284
1285   Eigen::SelfAdjointEigenSolver<typename CovarianceMatrix::PlainObject> myDiag{aHitCovarianceM};
1286   // eigenvectors of real symmetric matrices are orthogonal, i.e. invert == transpose
1287   auto aTranfoToDiagonalSystemInv =
1288       myDiag.eigenvectors().transpose().template cast<typename LocalDerivativeMatrix::Scalar>();
1289
1290   aHitCovarianceM = myDiag.eigenvalues().asDiagonal();
1291   aLocalDerivativesM = aTranfoToDiagonalSystemInv * aLocalDerivativesM;
1292   aHitResidualsM = aTranfoToDiagonalSystemInv * aHitResidualsM;
1293   if (aGlobalDerivativesM.size() > 0) {
1294     // diagonalize only if measurement depends on alignables or calibrations
1295     aGlobalDerivativesM = aTranfoToDiagonalSystemInv * aGlobalDerivativesM;
1296   }
1297 }
1298
1299 //__________________________________________________________________________________________________
1300 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1301 void MillePedeAlignmentAlgorithm ::addRefTrackVirtualMeas1D(
1302     const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1303     unsigned int iVirtualMeas,
1304     Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1305     Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1306     Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1307   // This Method is valid for 1D measurements only
1308
1309   const unsigned int xIndex = iVirtualMeas + refTrajPtr->numberOfHitMeas();
1310
1311   aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1312   aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex];
1313
1314   const auto &locDerivMatrix = refTrajPtr->derivatives();
1315   for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1316     aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1317   }
1318 }
1319
1320 //__________________________________________________________________________________________________
1321 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1322 void MillePedeAlignmentAlgorithm ::addRefTrackData2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1323                                                      unsigned int iTrajHit,
1324                                                      Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1325                                                      Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1326                                                      Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1327   // This Method is valid for 2D measurements only
1328
1329   const unsigned int xIndex = iTrajHit * 2;
1330   const unsigned int yIndex = iTrajHit * 2 + 1;
1331
1332   aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1333   aHitCovarianceM(0, 1) = refTrajPtr->measurementErrors()[xIndex][yIndex];
1334   aHitCovarianceM(1, 0) = refTrajPtr->measurementErrors()[yIndex][xIndex];
1335   aHitCovarianceM(1, 1) = refTrajPtr->measurementErrors()[yIndex][yIndex];
1336
1337   aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1338   aHitResidualsM(1, 0) = refTrajPtr->measurements()[yIndex] - refTrajPtr->trajectoryPositions()[yIndex];
1339
1340   const auto &locDerivMatrix = refTrajPtr->derivatives();
1341   for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1342     aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1343     aLocalDerivativesM(1, i) = locDerivMatrix[yIndex][i];
1344   }
1345 }
1346
1347 //__________________________________________________________________________________________________
1348 int MillePedeAlignmentAlgorithm ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1349                                             unsigned int iTrajHit,
1350                                             const std::vector<int> &globalLabels,
1351                                             const std::vector<float> &globalDerivativesX,
1352                                             const std::vector<float> &globalDerivativesY) {
1353   const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1354
1355   if ((aRecHit)->dimension() == 1) {
1356     return this->callMille1D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX);
1357   } else {
1358     return this->callMille2D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX, globalDerivativesY);
1359   }
1360 }
1361
1362 //__________________________________________________________________________________________________
1363 int MillePedeAlignmentAlgorithm ::callMille1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1364                                               unsigned int iTrajHit,
1365                                               const std::vector<int> &globalLabels,
1366                                               const std::vector<float> &globalDerivativesX) {
1367   const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1368   const unsigned int xIndex = iTrajHit * 2;  // the even ones are local x
1369
1370   // local derivatives
1371   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1372   const int nLocal = locDerivMatrix.num_col();
1373   std::vector<float> localDerivatives(nLocal);
1374   for (unsigned int i = 0; i < localDerivatives.size(); ++i) {
1375     localDerivatives[i] = locDerivMatrix[xIndex][i];
1376   }
1377
1378   // residuum and error
1379   float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1380   float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
1381
1382   // number of global derivatives
1383   const int nGlobal = globalDerivativesX.size();
1384
1385   // &(localDerivatives[0]) etc. are valid - as long as vector is not empty
1386   // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1387   theMille->mille(
1388       nLocal, &(localDerivatives[0]), nGlobal, &(globalDerivativesX[0]), &(globalLabels[0]), residX, hitErrX);
1389
1390   if (theMonitor) {
1391     theMonitor->fillDerivatives(
1392         aRecHit, &(localDerivatives[0]), nLocal, &(globalDerivativesX[0]), nGlobal, &(globalLabels[0]));
1393     theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, residX, hitErrX, false);
1394   }
1395
1396   return 1;
1397 }
1398
1399 //__________________________________________________________________________________________________
1400 int MillePedeAlignmentAlgorithm ::callMille2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1401                                               unsigned int iTrajHit,
1402                                               const std::vector<int> &globalLabels,
1403                                               const std::vector<float> &globalDerivativesx,
1404                                               const std::vector<float> &globalDerivativesy) {
1405   const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1406
1407   if ((aRecHit)->dimension() != 2) {
1408     edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::callMille2D"
1409                                << "You try to call method for 2D hits for a " << (aRecHit)->dimension()
1410                                << "D Hit. Hit gets ignored!";
1411     return -1;
1412   }
1413
1414   Eigen::Matrix<double, 2, 2> aHitCovarianceM;
1415   Eigen::Matrix<float, 2, 1> aHitResidualsM;
1416   Eigen::Matrix<float, 2, Eigen::Dynamic> aLocalDerivativesM{2, refTrajPtr->derivatives().num_col()};
1417   // below method fills above 3 matrices
1418   this->addRefTrackData2D(refTrajPtr, iTrajHit, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1419   Eigen::Matrix<float, 2, Eigen::Dynamic> aGlobalDerivativesM{2, globalDerivativesx.size()};
1420   this->makeGlobDerivMatrix(globalDerivativesx, globalDerivativesy, aGlobalDerivativesM);
1421
1422   // calculates correlation between Hit measurements
1423   // FIXME: Should take correlation (and resulting transformation) from original hit,
1424   //        not 2x2 matrix from ReferenceTrajectory: That can come from error propagation etc.!
1425   const double corr = aHitCovarianceM(0, 1) / sqrt(aHitCovarianceM(0, 0) * aHitCovarianceM(1, 1));
1426   if (theMonitor)
1427     theMonitor->fillCorrelations2D(corr, aRecHit);
1428   bool diag = false;  // diagonalise only tracker TID, TEC
1429   switch (aRecHit->geographicalId().subdetId()) {
1430     case SiStripDetId::TID:
1431     case SiStripDetId::TEC:
1432       if (aRecHit->geographicalId().det() == DetId::Tracker && TMath::Abs(corr) > theMaximalCor2D) {
1433         this->diagonalize(aHitCovarianceM, aLocalDerivativesM, aHitResidualsM, aGlobalDerivativesM);
1434         diag = true;
1435       }
1436       break;
1437     default:;
1438   }
1439
1440   float newResidX = aHitResidualsM(0, 0);
1441   float newResidY = aHitResidualsM(1, 0);
1442   float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1443   float newHitErrY = TMath::Sqrt(aHitCovarianceM(1, 1));
1444
1445   // change from column major (Eigen default) to row major to have row entries
1446   // in continuous memory
1447   std::vector<float> newLocalDerivs(aLocalDerivativesM.size());
1448   Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1449       newLocalDerivs.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1450   float *newLocalDerivsX = &(newLocalDerivs[0]);
1451   float *newLocalDerivsY = &(newLocalDerivs[aLocalDerivativesM.cols()]);
1452
1453   // change from column major (Eigen default) to row major to have row entries
1454   // in continuous memory
1455   std::vector<float> newGlobDerivs(aGlobalDerivativesM.size());
1456   Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1457       newGlobDerivs.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1458   float *newGlobDerivsX = &(newGlobDerivs[0]);
1459   float *newGlobDerivsY = &(newGlobDerivs[aGlobalDerivativesM.cols()]);
1460
1461   const int nLocal = aLocalDerivativesM.cols();
1462   const int nGlobal = aGlobalDerivativesM.cols();
1463
1464   if (diag && (newHitErrX > newHitErrY)) {  // also for 2D hits?
1465     // measurement with smaller error is x-measurement (for !is2D do not fill y-measurement):
1466     std::swap(newResidX, newResidY);
1467     std::swap(newHitErrX, newHitErrY);
1468     std::swap(newLocalDerivsX, newLocalDerivsY);
1469     std::swap(newGlobDerivsX, newGlobDerivsY);
1470   }
1471
1472   // &(globalLabels[0]) is valid - as long as vector is not empty
1473   // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1474   theMille->mille(nLocal, newLocalDerivsX, nGlobal, newGlobDerivsX, &(globalLabels[0]), newResidX, newHitErrX);
1475
1476   if (theMonitor) {
1477     theMonitor->fillDerivatives(aRecHit, newLocalDerivsX, nLocal, newGlobDerivsX, nGlobal, &(globalLabels[0]));
1478     theMonitor->fillResiduals(
1479         aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidX, newHitErrX, false);
1480   }
1481   const bool isReal2DHit = this->is2D(aRecHit);  // strip is 1D (except matched hits)
1482   if (isReal2DHit) {
1483     theMille->mille(nLocal, newLocalDerivsY, nGlobal, newGlobDerivsY, &(globalLabels[0]), newResidY, newHitErrY);
1484     if (theMonitor) {
1485       theMonitor->fillDerivatives(aRecHit, newLocalDerivsY, nLocal, newGlobDerivsY, nGlobal, &(globalLabels[0]));
1486       theMonitor->fillResiduals(
1487           aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidY, newHitErrY, true);  // true: y
1488     }
1489   }
1490
1491   return (isReal2DHit ? 2 : 1);
1492 }
1493
1494 //__________________________________________________________________________________________________
1495 void MillePedeAlignmentAlgorithm ::addVirtualMeas(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1496                                                   unsigned int iVirtualMeas) {
1497   Eigen::Matrix<double, 1, 1> aHitCovarianceM;
1498   Eigen::Matrix<float, 1, 1> aHitResidualsM;
1499   Eigen::Matrix<float, 1, Eigen::Dynamic> aLocalDerivativesM{1, refTrajPtr->derivatives().num_col()};
1500   // below method fills above 3 'matrices'
1501   this->addRefTrackVirtualMeas1D(refTrajPtr, iVirtualMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1502
1503   // no global parameters (use dummy 0)
1504   auto aGlobalDerivativesM = Eigen::Matrix<float, 1, 1>::Zero();
1505
1506   float newResidX = aHitResidualsM(0, 0);
1507   float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1508   std::vector<float> newLocalDerivsX(aLocalDerivativesM.size());
1509   Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1510       newLocalDerivsX.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1511
1512   std::vector<float> newGlobDerivsX(aGlobalDerivativesM.size());
1513   Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1514       newGlobDerivsX.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1515
1516   const int nLocal = aLocalDerivativesM.cols();
1517   const int nGlobal = 0;
1518
1519   theMille->mille(nLocal, newLocalDerivsX.data(), nGlobal, newGlobDerivsX.data(), &nGlobal, newResidX, newHitErrX);
1520 }
1521
1522 //____________________________________________________
1523 void MillePedeAlignmentAlgorithm::addLaserData(const EventInfo &eventInfo,
1524                                                const TkFittedLasBeamCollection &lasBeams,
1525                                                const TsosVectorCollection &lasBeamTsoses) {
1526   TsosVectorCollection::const_iterator iTsoses = lasBeamTsoses.begin();
1527   for (TkFittedLasBeamCollection::const_iterator iBeam = lasBeams.begin(), iEnd = lasBeams.end(); iBeam != iEnd;
1528        ++iBeam, ++iTsoses) {  // beam/tsoses parallel!
1529
1530     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addLaserData"
1531                               << "Beam " << iBeam->getBeamId() << " with " << iBeam->parameters().size()
1532                               << " parameters and " << iBeam->getData().size() << " hits.\n There are "
1533                               << iTsoses->size() << " TSOSes.";
1534
1535     this->addLasBeam(eventInfo, *iBeam, *iTsoses);
1536   }
1537 }
1538
1539 //____________________________________________________
1540 void MillePedeAlignmentAlgorithm::addLasBeam(const EventInfo &eventInfo,
1541                                              const TkFittedLasBeam &lasBeam,
1542                                              const std::vector<TrajectoryStateOnSurface> &tsoses) {
1543   AlignmentParameters *dummyPtr = nullptr;                                          // for globalDerivativesHierarchy()
1544   std::vector<float> lasLocalDerivsX;                                               // buffer for local derivatives
1545   const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId());  // for global par
1546
1547   for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
1548     if (!tsoses[iHit].isValid())
1549       continue;
1550     // clear buffer
1551     theFloatBufferX.clear();
1552     theFloatBufferY.clear();
1553     theIntBuffer.clear();
1554     lasLocalDerivsX.clear();
1555     // get alignables and global parameters
1556     const SiStripLaserRecHit2D &hit = lasBeam.getData()[iHit];
1557     AlignableDetOrUnitPtr lasAli(theAlignableNavigator->alignableFromDetId(hit.getDetId()));
1558     this->globalDerivativesHierarchy(
1559         eventInfo, tsoses[iHit], lasAli, lasAli, theFloatBufferX, theFloatBufferY, theIntBuffer, dummyPtr);
1560     // fill derivatives vector from derivatives matrix
1561     for (unsigned int nFitParams = 0; nFitParams < static_cast<unsigned int>(lasBeam.parameters().size());
1562          ++nFitParams) {
1563       const float derivative = lasBeam.derivatives()[iHit][nFitParams];
1564       if (nFitParams < lasBeam.firstFixedParameter()) {  // first local beam parameters
1565         lasLocalDerivsX.push_back(derivative);
1566       } else {  // now global ones
1567         const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
1568         theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
1569         theFloatBufferX.push_back(derivative);
1570       }
1571     }  // end loop over parameters
1572
1573     const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
1574     // error from file or assume 0.003
1575     const float error = 0.003;  // hit.localPositionError().xx(); sqrt???
1576
1577     theMille->mille(lasLocalDerivsX.size(),
1578                     &(lasLocalDerivsX[0]),
1579                     theFloatBufferX.size(),
1580                     &(theFloatBufferX[0]),
1581                     &(theIntBuffer[0]),
1582                     residual,
1583                     error);
1584   }  // end of loop over hits
1585
1586   theMille->end();
1587 }
1588
1589 void MillePedeAlignmentAlgorithm::addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg) {
1590   // do some printing, if requested
1591   const bool doOutputOnStdout(pxbSurveyCfg.getParameter<bool>("doOutputOnStdout"));
1592   if (doOutputOnStdout) {
1593     edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1594                               << "# Output from addPxbSurvey follows below because "
1595                               << "doOutputOnStdout is set to True";
1596   }
1597
1598   // instantiate a dicer object
1599   SurveyPxbDicer dicer(pxbSurveyCfg.getParameter<std::vector<edm::ParameterSet> >("toySurveyParameters"),
1600                        pxbSurveyCfg.getParameter<unsigned int>("toySurveySeed"));
1601   std::ofstream outfile(pxbSurveyCfg.getUntrackedParameter<std::string>("toySurveyFile").c_str());
1602
1603   // read data from file
1604   std::vector<SurveyPxbImageLocalFit> measurements;
1605   std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
1606   SurveyPxbImageReader<SurveyPxbImageLocalFit> reader(filename, measurements, 800);
1607
1608   // loop over photographs (=measurements) and perform the fit
1609   for (std::vector<SurveyPxbImageLocalFit>::size_type i = 0; i != measurements.size(); i++) {
1610     if (doOutputOnStdout) {
1611       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1612                                 << "Module " << i << ": ";
1613     }
1614
1615     // get the Alignables and their surfaces
1616     AlignableDetOrUnitPtr mod1(theAlignableNavigator->alignableFromDetId(measurements[i].getIdFirst()));
1617     AlignableDetOrUnitPtr mod2(theAlignableNavigator->alignableFromDetId(measurements[i].getIdSecond()));
1618     const AlignableSurface &surf1 = mod1->surface();
1619     const AlignableSurface &surf2 = mod2->surface();
1620
1621     // the position of the fiducial points in local frame of a PXB module
1622     const LocalPoint fidpoint0(-0.91, +3.30);
1623     const LocalPoint fidpoint1(+0.91, +3.30);
1624     const LocalPoint fidpoint2(+0.91, -3.30);
1625     const LocalPoint fidpoint3(-0.91, -3.30);
1626
1627     // We choose the local frame of the first module as reference,
1628     // so take the fidpoints of the second module and calculate their
1629     // positions in the reference frame
1630     const GlobalPoint surf2point0(surf2.toGlobal(fidpoint0));
1631     const GlobalPoint surf2point1(surf2.toGlobal(fidpoint1));
1632     const LocalPoint fidpoint0inSurf1frame(surf1.toLocal(surf2point0));
1633     const LocalPoint fidpoint1inSurf1frame(surf1.toLocal(surf2point1));
1634
1635     // Create the vector for the fit
1636     SurveyPxbImageLocalFit::fidpoint_t fidpointvec;
1637     fidpointvec.push_back(fidpoint0inSurf1frame);
1638     fidpointvec.push_back(fidpoint1inSurf1frame);
1639     fidpointvec.push_back(fidpoint2);
1640     fidpointvec.push_back(fidpoint3);
1641
1642     // if toy survey is requested, dice the values now
1643     if (pxbSurveyCfg.getParameter<bool>("doToySurvey")) {
1644       dicer.doDice(fidpointvec, measurements[i].getIdPair(), outfile);
1645     }
1646
1647     // do the fit
1648     measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
1649     SurveyPxbImageLocalFit::localpars_t a;  // local pars from fit
1650     a = measurements[i].getLocalParameters();
1651     const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
1652
1653     // do some reporting, if requested
1654     if (doOutputOnStdout) {
1655       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1656                                 << "a: " << a[0] << ", " << a[1] << ", " << a[2] << ", " << a[3]
1657                                 << " S= " << sqrt(a[2] * a[2] + a[3] * a[3]) << " phi= " << atan(a[3] / a[2])
1658                                 << " chi2= " << chi2 << std::endl;
1659     }
1660     if (theMonitor) {
1661       theMonitor->fillPxbSurveyHistsChi2(chi2);
1662       theMonitor->fillPxbSurveyHistsLocalPars(a[0], a[1], sqrt(a[2] * a[2] + a[3] * a[3]), atan(a[3] / a[2]));
1663     }
1664
1665     // pass the results from the local fit to mille
1666     for (SurveyPxbImageLocalFit::count_t j = 0; j != SurveyPxbImageLocalFit::nMsrmts; j++) {
1667       theMille->mille((int)measurements[i].getLocalDerivsSize(),
1668                       measurements[i].getLocalDerivsPtr(j),
1669                       (int)measurements[i].getGlobalDerivsSize(),
1670                       measurements[i].getGlobalDerivsPtr(j),
1671                       measurements[i].getGlobalDerivsLabelPtr(j),
1672                       measurements[i].getResiduum(j),
1673                       measurements[i].getSigma(j));
1674     }
1675     theMille->end();
1676   }
1677   outfile.close();
1678 }
1679
1680 bool MillePedeAlignmentAlgorithm::areIOVsSpecified() const {
1681   const auto runRangeSelection = theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection");
1682
1683   if (runRangeSelection.empty())
1684     return false;
1685
1686   const auto runRanges =
1687       align::makeNonOverlappingRunRanges(runRangeSelection, cond::timeTypeSpecs[cond::runnumber].beginValue);
1688
1689   return !(runRanges.empty());
1690 }