File indexing completed on 2025-01-04 00:29:04
0001
0002
0003
0004
0005
0006
0007
0008
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
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
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
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 += '/';
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());
0126
0127 theBinary = std::make_unique<MilleBinary>((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str(),
0128 theGblDoubleBinary);
0129 }
0130 }
0131
0132
0133
0134 MillePedeAlignmentAlgorithm::~MillePedeAlignmentAlgorithm() {}
0135
0136
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
0184 const TrackerTopology *const tTopo = &setup.getData(topoToken_);
0185
0186
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
0193 const SiPixelQuality &qual = setup.getData(siPixelQualityToken_);
0194
0195 pixelQuality = std::make_shared<SiPixelQuality>(qual);
0196
0197
0198 const TrackerGeometry *tGeom = &setup.getData(geomToken_);
0199
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
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
0251
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
0261
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
0266
0267 if (!this->readFromPede((*iSet), false, runrange)) {
0268 throw cms::Exception("BadConfig")
0269 << "MillePedeAlignmentAlgorithm::initialize: Problems reading input constants of "
0270 << "pedeReaderInputs entry " << iSet - mprespset.begin() << '.';
0271 }
0272 theAlignmentParameterStore->applyParameters();
0273
0274 theAlignmentParameterStore->resetParameters();
0275 }
0276
0277
0278 thePedeSteer->buildSubSteer(tracker, muon, extras);
0279
0280
0281 this->buildUserVariables(theAlignables);
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
0294 }
0295
0296 if (this->isMode(myPedeSteerBit)) {
0297
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
0359 if (enableAlignableUpdates_) {
0360 theAlignmentParameterStore->restoreCachedTransformations(runrange.first);
0361 } else {
0362 theAlignmentParameterStore->restoreCachedTransformations();
0363 }
0364
0365
0366 theAlignmentParameterStore->resetParameters();
0367
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);
0377 }
0378
0379 return true;
0380 }
0381
0382
0383
0384 void MillePedeAlignmentAlgorithm::terminate(const edm::EventSetup &iSetup) { terminate(); }
0385 void MillePedeAlignmentAlgorithm::terminate() {
0386 theMille.reset();
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
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
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));
0432 if (this->isMode(myPedeRunBit)) {
0433 thePedeSteer->runPede(masterSteer);
0434 }
0435
0436
0437
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
0450
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
0457 struct stat buffer;
0458 if (stat(strippedInputFileName.c_str(), &buffer) == 0) {
0459
0460 files.push_back(theCompleteInputFileName);
0461 if (theNumberedInputFileName == theInputFileName) {
0462
0463 break;
0464 } else {
0465
0466 theNumber++;
0467 }
0468 } else {
0469
0470 break;
0471 }
0472 }
0473
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
0482
0483 void MillePedeAlignmentAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo) {
0484 if (!this->isMode(myMilleBit))
0485 return;
0486 const auto &tracks = eventInfo.trajTrackPairs();
0487
0488 if (theMonitor) {
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
0497 unsigned int refTrajCount = 0;
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
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 }
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
0524 if (!refTrajPtr->gblInput().empty()) {
0525
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
0541 if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits)
0542 return hitResultXy;
0543
0544 if (refTrajPtr->gblInput().size() == 1) {
0545
0546 GblTrajectory aGblTrajectory(refTrajPtr->gblInput()[0].first, refTrajPtr->nominalField() != 0);
0547
0548
0549
0550
0551
0552
0553
0554 if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
0555 aGblTrajectory.milleOut(*theBinary);
0556 }
0557 if (refTrajPtr->gblInput().size() == 2) {
0558
0559 GblTrajectory aGblTrajectory(refTrajPtr->gblInput(),
0560 refTrajPtr->gblExtDerivatives(),
0561 refTrajPtr->gblExtMeasurements(),
0562 refTrajPtr->gblExtPrecisions());
0563
0564 if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
0565 aGblTrajectory.milleOut(*theBinary);
0566 }
0567 } else {
0568
0569 std::vector<AlignmentParameters *> parVec(refTrajPtr->recHits().size());
0570
0571 std::vector<bool> validHitVecY(refTrajPtr->recHits().size(), false);
0572
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) {
0577 hitResultXy.first = 0;
0578 break;
0579 } else {
0580 if (flagXY >= 1)
0581 ++hitResultXy.first;
0582 validHitVecY[iHit] = (flagXY >= 2);
0583 }
0584 }
0585
0586
0587 for (unsigned int iVirtualMeas = 0; iVirtualMeas < refTrajPtr->numberOfVirtualMeas(); ++iVirtualMeas) {
0588 this->addVirtualMeas(refTrajPtr, iVirtualMeas);
0589 }
0590
0591
0592 if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits) {
0593 theMille->kill();
0594 hitResultXy.first = hitResultXy.second = 0;
0595 } else {
0596 theMille->end();
0597
0598
0599 hitResultXy.second = this->addHitCount(parVec, validHitVecY);
0600
0601 }
0602 }
0603
0604 }
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
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
0617
0618
0619 while (ali) {
0620 AlignmentParameters *pars = ali->alignmentParameters();
0621 if (pars) {
0622
0623 MillePedeVariables *mpVar = static_cast<MillePedeVariables *>(pars->userVariables());
0624
0625 mpVar->increaseHitsX();
0626 if (validHitVecY[iHit]) {
0627 mpVar->increaseHitsY();
0628 if (pars == parVec[iHit])
0629 ++nHitY;
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
0717 this->addLaserData(eventInfo, *(runInfo.tkLasBeams()), *(runInfo.tkLasBeamTsoses()));
0718 }
0719 if (this->isMode(myMilleBit))
0720 theMille->flushOutputFile();
0721 }
0722
0723
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();
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 *¶ms) {
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
0763 if (!recHitPtr->isValid())
0764 return 0;
0765
0766
0767
0768 this->globalDerivativesCalibration(recHitPtr,
0769 tsos,
0770 setup,
0771 eventInfo,
0772 theFloatBufferX,
0773 theFloatBufferY,
0774 theIntBuffer);
0775
0776
0777 AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0778
0779 if (!this->globalDerivativesHierarchy(eventInfo,
0780 tsos,
0781 alidet,
0782 alidet,
0783 theFloatBufferX,
0784 theFloatBufferY,
0785 theIntBuffer,
0786 params)) {
0787 return -1;
0788 } else if (theFloatBufferX.empty() && ignoreHitsWithoutGlobalDerivatives_) {
0789 return 0;
0790 } else {
0791
0792
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
0814 if (!recHitPtr->isValid())
0815 return 0;
0816
0817
0818 AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0819
0820 if (!this->globalDerivativesHierarchy(eventInfo,
0821 tsos,
0822 alidet,
0823 alidet,
0824 theDoubleBufferX,
0825 theDoubleBufferY,
0826 theIntBuffer,
0827 params)) {
0828 return -1;
0829 }
0830
0831 int globalLabel;
0832 std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0833 for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0834
0835 (*iCalib)->derivatives(derivs, *recHitPtr, tsos, setup, eventInfo);
0836 for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0837
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
0872 if (!ali)
0873 return true;
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;
0883
0884 bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0885 const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0886
0887 if (0 == alignableLabel) {
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
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
0909 if (thePedeSteer->isNoHiera(ali))
0910 return true;
0911 }
0912
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
0927 if (!ali)
0928 return true;
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;
0938
0939 bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0940 const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0941
0942 if (0 == alignableLabel) {
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
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
0971 if (thePedeSteer->isNoHiera(ali))
0972 return true;
0973 }
0974
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
0990 (*iCalib)->derivatives(derivs, *recHit, tsos, setup, eventInfo);
0991 for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0992
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
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035 bool MillePedeAlignmentAlgorithm::is2D(const ConstRecHitPointer &recHit) const {
1036
1037
1038 if (recHit->dimension() < 2) {
1039 return false;
1040 } else if (recHit->detUnit()) {
1041 return recHit->detUnit()->type().isTrackerPixel();
1042 } else {
1043 if (dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit->hit())) {
1044
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);
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;
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) {
1077 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1078 } else if (!alis.empty()) {
1079 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1080 } else {
1081 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1082 return false;
1083 }
1084 return true;
1085 }
1086
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
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);
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) {
1179 for (unsigned int iPar = 0; iPar < userVars->size(); ++iPar) {
1180
1181
1182
1183 userVars->setAllDefault(iPar);
1184
1185 }
1186 } else {
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;
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;
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;
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;
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
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
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
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
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;
1377
1378
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
1387 float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1388 float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
1389
1390
1391 const int nGlobal = globalDerivativesX.size();
1392
1393
1394
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
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
1431
1432
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;
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
1454
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
1462
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)) {
1473
1474 std::swap(newResidX, newResidY);
1475 std::swap(newHitErrX, newHitErrY);
1476 std::swap(newLocalDerivsX, newLocalDerivsY);
1477 std::swap(newGlobDerivsX, newGlobDerivsY);
1478 }
1479
1480
1481
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);
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);
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
1509 this->addRefTrackVirtualMeas1D(refTrajPtr, iVirtualMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1510
1511
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) {
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;
1552 std::vector<float> lasLocalDerivsX;
1553 const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId());
1554
1555 for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
1556 if (!tsoses[iHit].isValid())
1557 continue;
1558
1559 theFloatBufferX.clear();
1560 theFloatBufferY.clear();
1561 theIntBuffer.clear();
1562 lasLocalDerivsX.clear();
1563
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
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()) {
1573 lasLocalDerivsX.push_back(derivative);
1574 } else {
1575 const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
1576 theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
1577 theFloatBufferX.push_back(derivative);
1578 }
1579 }
1580
1581 const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
1582
1583 const float error = 0.003;
1584
1585 theMille->mille(lasLocalDerivsX.size(),
1586 &(lasLocalDerivsX[0]),
1587 theFloatBufferX.size(),
1588 &(theFloatBufferX[0]),
1589 &(theIntBuffer[0]),
1590 residual,
1591 error);
1592 }
1593
1594 theMille->end();
1595 }
1596
1597 void MillePedeAlignmentAlgorithm::addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg) {
1598
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
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
1612 std::vector<SurveyPxbImageLocalFit> measurements;
1613 std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
1614 SurveyPxbImageReader<SurveyPxbImageLocalFit> reader(filename, measurements, 800);
1615
1616
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
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
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
1636
1637
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
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
1651 if (pxbSurveyCfg.getParameter<bool>("doToySurvey")) {
1652 dicer.doDice(fidpointvec, measurements[i].getIdPair(), outfile);
1653 }
1654
1655
1656 measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
1657 SurveyPxbImageLocalFit::localpars_t a;
1658 a = measurements[i].getLocalParameters();
1659 const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
1660
1661
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
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 }