File indexing completed on 2023-03-17 10:39:28
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 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 += '/';
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());
0124
0125 theBinary = std::make_unique<MilleBinary>((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str(),
0126 theGblDoubleBinary);
0127 }
0128 }
0129
0130
0131
0132 MillePedeAlignmentAlgorithm::~MillePedeAlignmentAlgorithm() {}
0133
0134
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
0182 const TrackerTopology *const tTopo = &setup.getData(topoToken_);
0183
0184
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
0191 const TrackerGeometry *tGeom = &setup.getData(geomToken_);
0192
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
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
0244
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
0254
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
0259
0260 if (!this->readFromPede((*iSet), false, runrange)) {
0261 throw cms::Exception("BadConfig")
0262 << "MillePedeAlignmentAlgorithm::initialize: Problems reading input constants of "
0263 << "pedeReaderInputs entry " << iSet - mprespset.begin() << '.';
0264 }
0265 theAlignmentParameterStore->applyParameters();
0266
0267 theAlignmentParameterStore->resetParameters();
0268 }
0269
0270
0271 thePedeSteer->buildSubSteer(tracker, muon, extras);
0272
0273
0274 this->buildUserVariables(theAlignables);
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
0287 }
0288
0289 if (this->isMode(myPedeSteerBit)) {
0290
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
0351 if (enableAlignableUpdates_) {
0352 theAlignmentParameterStore->restoreCachedTransformations(runrange.first);
0353 } else {
0354 theAlignmentParameterStore->restoreCachedTransformations();
0355 }
0356
0357
0358 theAlignmentParameterStore->resetParameters();
0359
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);
0369 }
0370
0371 return true;
0372 }
0373
0374
0375
0376 void MillePedeAlignmentAlgorithm::terminate(const edm::EventSetup &iSetup) { terminate(); }
0377 void MillePedeAlignmentAlgorithm::terminate() {
0378 theMille.reset();
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
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
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));
0424 if (this->isMode(myPedeRunBit)) {
0425 thePedeSteer->runPede(masterSteer);
0426 }
0427
0428
0429
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
0442
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
0449 struct stat buffer;
0450 if (stat(strippedInputFileName.c_str(), &buffer) == 0) {
0451
0452 files.push_back(theCompleteInputFileName);
0453 if (theNumberedInputFileName == theInputFileName) {
0454
0455 break;
0456 } else {
0457
0458 theNumber++;
0459 }
0460 } else {
0461
0462 break;
0463 }
0464 }
0465
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
0474
0475 void MillePedeAlignmentAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo) {
0476 if (!this->isMode(myMilleBit))
0477 return;
0478 const auto &tracks = eventInfo.trajTrackPairs();
0479
0480 if (theMonitor) {
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
0489 unsigned int refTrajCount = 0;
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
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 }
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
0516 if (!refTrajPtr->gblInput().empty()) {
0517
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
0533 if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits)
0534 return hitResultXy;
0535
0536 if (refTrajPtr->gblInput().size() == 1) {
0537
0538 GblTrajectory aGblTrajectory(refTrajPtr->gblInput()[0].first, refTrajPtr->nominalField() != 0);
0539
0540
0541
0542
0543
0544
0545
0546 if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
0547 aGblTrajectory.milleOut(*theBinary);
0548 }
0549 if (refTrajPtr->gblInput().size() == 2) {
0550
0551 GblTrajectory aGblTrajectory(refTrajPtr->gblInput(),
0552 refTrajPtr->gblExtDerivatives(),
0553 refTrajPtr->gblExtMeasurements(),
0554 refTrajPtr->gblExtPrecisions());
0555
0556 if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
0557 aGblTrajectory.milleOut(*theBinary);
0558 }
0559 } else {
0560
0561 std::vector<AlignmentParameters *> parVec(refTrajPtr->recHits().size());
0562
0563 std::vector<bool> validHitVecY(refTrajPtr->recHits().size(), false);
0564
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) {
0569 hitResultXy.first = 0;
0570 break;
0571 } else {
0572 if (flagXY >= 1)
0573 ++hitResultXy.first;
0574 validHitVecY[iHit] = (flagXY >= 2);
0575 }
0576 }
0577
0578
0579 for (unsigned int iVirtualMeas = 0; iVirtualMeas < refTrajPtr->numberOfVirtualMeas(); ++iVirtualMeas) {
0580 this->addVirtualMeas(refTrajPtr, iVirtualMeas);
0581 }
0582
0583
0584 if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits) {
0585 theMille->kill();
0586 hitResultXy.first = hitResultXy.second = 0;
0587 } else {
0588 theMille->end();
0589
0590
0591 hitResultXy.second = this->addHitCount(parVec, validHitVecY);
0592
0593 }
0594 }
0595
0596 }
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
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
0609
0610
0611 while (ali) {
0612 AlignmentParameters *pars = ali->alignmentParameters();
0613 if (pars) {
0614
0615 MillePedeVariables *mpVar = static_cast<MillePedeVariables *>(pars->userVariables());
0616
0617 mpVar->increaseHitsX();
0618 if (validHitVecY[iHit]) {
0619 mpVar->increaseHitsY();
0620 if (pars == parVec[iHit])
0621 ++nHitY;
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
0709 this->addLaserData(eventInfo, *(runInfo.tkLasBeams()), *(runInfo.tkLasBeamTsoses()));
0710 }
0711 if (this->isMode(myMilleBit))
0712 theMille->flushOutputFile();
0713 }
0714
0715
0716 void MillePedeAlignmentAlgorithm::endRun(const EndRunInfo &runInfo, const edm::EventSetup &setup) {
0717 if (this->isMode(myMilleBit))
0718 theMille->flushOutputFile();
0719 }
0720
0721
0722 void MillePedeAlignmentAlgorithm::beginLuminosityBlock(const edm::EventSetup &) {
0723 if (!runAtPCL_)
0724 return;
0725 if (this->isMode(myMilleBit))
0726 theMille->resetOutputFile();
0727 }
0728
0729
0730 void MillePedeAlignmentAlgorithm::endLuminosityBlock(const edm::EventSetup &) {
0731 if (!runAtPCL_)
0732 return;
0733 if (this->isMode(myMilleBit))
0734 theMille->flushOutputFile();
0735 }
0736
0737
0738 int MillePedeAlignmentAlgorithm::addMeasurementData(const edm::EventSetup &setup,
0739 const EventInfo &eventInfo,
0740 const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0741 unsigned int iHit,
0742 AlignmentParameters *¶ms) {
0743 params = nullptr;
0744 theFloatBufferX.clear();
0745 theFloatBufferY.clear();
0746 theIntBuffer.clear();
0747
0748 const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
0749 const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
0750
0751 if (!recHitPtr->isValid())
0752 return 0;
0753
0754
0755
0756 this->globalDerivativesCalibration(recHitPtr,
0757 tsos,
0758 setup,
0759 eventInfo,
0760 theFloatBufferX,
0761 theFloatBufferY,
0762 theIntBuffer);
0763
0764
0765 AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0766
0767 if (!this->globalDerivativesHierarchy(eventInfo,
0768 tsos,
0769 alidet,
0770 alidet,
0771 theFloatBufferX,
0772 theFloatBufferY,
0773 theIntBuffer,
0774 params)) {
0775 return -1;
0776 } else if (theFloatBufferX.empty() && ignoreHitsWithoutGlobalDerivatives_) {
0777 return 0;
0778 } else {
0779
0780
0781 return this->callMille(refTrajPtr, iHit, theIntBuffer, theFloatBufferX, theFloatBufferY);
0782 }
0783 }
0784
0785
0786
0787 int MillePedeAlignmentAlgorithm::addGlobalData(const edm::EventSetup &setup,
0788 const EventInfo &eventInfo,
0789 const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
0790 unsigned int iHit,
0791 GblPoint &gblPoint) {
0792 AlignmentParameters *params = nullptr;
0793 std::vector<double> theDoubleBufferX, theDoubleBufferY;
0794 theDoubleBufferX.clear();
0795 theDoubleBufferY.clear();
0796 theIntBuffer.clear();
0797 int iret = 0;
0798
0799 const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
0800 const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
0801
0802 if (!recHitPtr->isValid())
0803 return 0;
0804
0805
0806 AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
0807
0808 if (!this->globalDerivativesHierarchy(eventInfo,
0809 tsos,
0810 alidet,
0811 alidet,
0812 theDoubleBufferX,
0813 theDoubleBufferY,
0814 theIntBuffer,
0815 params)) {
0816 return -1;
0817 }
0818
0819 int globalLabel;
0820 std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0821 for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0822
0823 (*iCalib)->derivatives(derivs, *recHitPtr, tsos, setup, eventInfo);
0824 for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0825
0826 globalLabel = thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second);
0827 if (globalLabel > 0 && globalLabel <= 2147483647) {
0828 theIntBuffer.push_back(globalLabel);
0829 theDoubleBufferX.push_back(iValuesInd->first.first);
0830 theDoubleBufferY.push_back(iValuesInd->first.second);
0831 } else {
0832 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addGlobalData"
0833 << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
0834 }
0835 }
0836 }
0837 unsigned int numGlobals = theIntBuffer.size();
0838 if (numGlobals > 0) {
0839 Eigen::Matrix<double, 2, Eigen::Dynamic> globalDer{2, numGlobals};
0840 for (unsigned int i = 0; i < numGlobals; ++i) {
0841 globalDer(0, i) = theDoubleBufferX[i];
0842 globalDer(1, i) = theDoubleBufferY[i];
0843 }
0844 gblPoint.addGlobals(theIntBuffer, globalDer);
0845 iret = 1;
0846 }
0847 return iret;
0848 }
0849
0850
0851 bool MillePedeAlignmentAlgorithm ::globalDerivativesHierarchy(const EventInfo &eventInfo,
0852 const TrajectoryStateOnSurface &tsos,
0853 Alignable *ali,
0854 const AlignableDetOrUnitPtr &alidet,
0855 std::vector<float> &globalDerivativesX,
0856 std::vector<float> &globalDerivativesY,
0857 std::vector<int> &globalLabels,
0858 AlignmentParameters *&lowestParams) const {
0859
0860 if (!ali)
0861 return true;
0862
0863 if (false && theMonitor && alidet != ali)
0864 theMonitor->fillFrameToFrame(alidet, ali);
0865
0866 AlignmentParameters *params = ali->alignmentParameters();
0867
0868 if (params) {
0869 if (!lowestParams)
0870 lowestParams = params;
0871
0872 bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0873 const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0874
0875 if (0 == alignableLabel) {
0876 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0877 << "Label not found, skip Alignable.";
0878 return false;
0879 }
0880
0881 const std::vector<bool> &selPars = params->selector();
0882 const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
0883
0884
0885 for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
0886 if (selPars[iSel]) {
0887 globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
0888 if (hasSplitParameters == true) {
0889 globalLabels.push_back(thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos));
0890 } else {
0891 globalLabels.push_back(thePedeLabels->parameterLabel(alignableLabel, iSel));
0892 }
0893 globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
0894 }
0895 }
0896
0897 if (thePedeSteer->isNoHiera(ali))
0898 return true;
0899 }
0900
0901 return this->globalDerivativesHierarchy(
0902 eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
0903 }
0904
0905
0906 bool MillePedeAlignmentAlgorithm ::globalDerivativesHierarchy(const EventInfo &eventInfo,
0907 const TrajectoryStateOnSurface &tsos,
0908 Alignable *ali,
0909 const AlignableDetOrUnitPtr &alidet,
0910 std::vector<double> &globalDerivativesX,
0911 std::vector<double> &globalDerivativesY,
0912 std::vector<int> &globalLabels,
0913 AlignmentParameters *&lowestParams) const {
0914
0915 if (!ali)
0916 return true;
0917
0918 if (false && theMonitor && alidet != ali)
0919 theMonitor->fillFrameToFrame(alidet, ali);
0920
0921 AlignmentParameters *params = ali->alignmentParameters();
0922
0923 if (params) {
0924 if (!lowestParams)
0925 lowestParams = params;
0926
0927 bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
0928 const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
0929
0930 if (0 == alignableLabel) {
0931 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0932 << "Label not found, skip Alignable.";
0933 return false;
0934 }
0935
0936 const std::vector<bool> &selPars = params->selector();
0937 const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
0938 int globalLabel;
0939
0940
0941 for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
0942 if (selPars[iSel]) {
0943 if (hasSplitParameters == true) {
0944 globalLabel = thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos);
0945 } else {
0946 globalLabel = thePedeLabels->parameterLabel(alignableLabel, iSel);
0947 }
0948 if (globalLabel > 0 && globalLabel <= 2147483647) {
0949 globalLabels.push_back(globalLabel);
0950 globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
0951 globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
0952 } else {
0953 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
0954 << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
0955 }
0956 }
0957 }
0958
0959 if (thePedeSteer->isNoHiera(ali))
0960 return true;
0961 }
0962
0963 return this->globalDerivativesHierarchy(
0964 eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
0965 }
0966
0967
0968 void MillePedeAlignmentAlgorithm::globalDerivativesCalibration(const TransientTrackingRecHit::ConstRecHitPointer &recHit,
0969 const TrajectoryStateOnSurface &tsos,
0970 const edm::EventSetup &setup,
0971 const EventInfo &eventInfo,
0972 std::vector<float> &globalDerivativesX,
0973 std::vector<float> &globalDerivativesY,
0974 std::vector<int> &globalLabels) const {
0975 std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
0976 for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
0977
0978 (*iCalib)->derivatives(derivs, *recHit, tsos, setup, eventInfo);
0979 for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
0980
0981 globalLabels.push_back(thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second));
0982 globalDerivativesX.push_back(iValuesInd->first.first);
0983 globalDerivativesY.push_back(iValuesInd->first.second);
0984 }
0985 }
0986 }
0987
0988
0989
0990
0991
0992
0993
0994
0995
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 bool MillePedeAlignmentAlgorithm::is2D(const ConstRecHitPointer &recHit) const {
1024
1025
1026 if (recHit->dimension() < 2) {
1027 return false;
1028 } else if (recHit->detUnit()) {
1029 return recHit->detUnit()->type().isTrackerPixel();
1030 } else {
1031 if (dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit->hit())) {
1032
1033 return false;
1034 } else {
1035 return true;
1036 }
1037 }
1038 }
1039
1040
1041 bool MillePedeAlignmentAlgorithm::readFromPede(const edm::ParameterSet &mprespset,
1042 bool setUserVars,
1043 const RunRange &runrange) {
1044 bool allEmpty = this->areEmptyParams(theAlignables);
1045
1046 PedeReader reader(mprespset, *thePedeSteer, *thePedeLabels, runrange);
1047 align::Alignables alis;
1048 bool okRead = reader.read(alis, setUserVars);
1049 bool numMatch = true;
1050
1051 std::stringstream out;
1052 out << "Read " << alis.size() << " alignables";
1053 if (alis.size() != theAlignables.size()) {
1054 out << " while " << theAlignables.size() << " in store";
1055 numMatch = false;
1056 }
1057 if (!okRead)
1058 out << ", but problems in reading";
1059 if (!allEmpty)
1060 out << ", possibly overwriting previous settings";
1061 out << ".";
1062
1063 if (okRead && allEmpty) {
1064 if (numMatch) {
1065 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1066 } else if (!alis.empty()) {
1067 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1068 } else {
1069 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1070 return false;
1071 }
1072 return true;
1073 }
1074
1075 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1076 return false;
1077 }
1078
1079
1080 bool MillePedeAlignmentAlgorithm::areEmptyParams(const align::Alignables &alignables) const {
1081 for (const auto &iAli : alignables) {
1082 const AlignmentParameters *params = iAli->alignmentParameters();
1083 if (params) {
1084 const auto &parVec(params->parameters());
1085 const auto &parCov(params->covariance());
1086 for (int i = 0; i < parVec.num_row(); ++i) {
1087 if (parVec[i] != 0.)
1088 return false;
1089 for (int j = i; j < parCov.num_col(); ++j) {
1090 if (parCov[i][j] != 0.)
1091 return false;
1092 }
1093 }
1094 }
1095 }
1096
1097 return true;
1098 }
1099
1100
1101 unsigned int MillePedeAlignmentAlgorithm::doIO(int loop) const {
1102 unsigned int result = 0;
1103
1104 const std::string outFilePlain(theConfig.getParameter<std::string>("treeFile"));
1105 if (outFilePlain.empty()) {
1106 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1107 << "treeFile parameter empty => skip writing for 'loop' " << loop;
1108 return result;
1109 }
1110
1111 const std::string outFile(theDir + outFilePlain);
1112
1113 AlignmentIORoot aliIO;
1114 int ioerr = 0;
1115 if (loop == 0) {
1116 aliIO.writeAlignableOriginalPositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1117 if (ioerr) {
1118 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1119 << "Problem " << ioerr << " in writeAlignableOriginalPositions";
1120 ++result;
1121 }
1122 } else if (loop == 1) {
1123
1124 const std::vector<std::string> inFiles(theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles"));
1125 const std::vector<std::string> binFiles(theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
1126 if (inFiles.size() != binFiles.size()) {
1127 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1128 << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' "
1129 << "differ in size";
1130 }
1131 this->addHitStatistics(0, outFile, inFiles);
1132 }
1133 MillePedeVariablesIORoot millePedeIO;
1134 millePedeIO.writeMillePedeVariables(theAlignables, outFile.c_str(), loop, false, ioerr);
1135 if (ioerr) {
1136 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1137 << "Problem " << ioerr << " writing MillePedeVariables";
1138 ++result;
1139 }
1140
1141 aliIO.writeOrigRigidBodyAlignmentParameters(theAlignables, outFile.c_str(), loop, false, ioerr);
1142 if (ioerr) {
1143 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1144 << "Problem " << ioerr << " in writeOrigRigidBodyAlignmentParameters, " << loop;
1145 ++result;
1146 }
1147 aliIO.writeAlignableAbsolutePositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1148 if (ioerr) {
1149 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1150 << "Problem " << ioerr << " in writeAlignableAbsolutePositions, " << loop;
1151 ++result;
1152 }
1153
1154 return result;
1155 }
1156
1157
1158 void MillePedeAlignmentAlgorithm::buildUserVariables(const align::Alignables &alis) const {
1159 for (const auto &iAli : alis) {
1160 AlignmentParameters *params = iAli->alignmentParameters();
1161 if (!params) {
1162 throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::buildUserVariables"
1163 << "No parameters for alignable";
1164 }
1165 MillePedeVariables *userVars = dynamic_cast<MillePedeVariables *>(params->userVariables());
1166 if (userVars) {
1167 for (unsigned int iPar = 0; iPar < userVars->size(); ++iPar) {
1168
1169
1170
1171 userVars->setAllDefault(iPar);
1172
1173 }
1174 } else {
1175 userVars = new MillePedeVariables(
1176 params->size(),
1177 thePedeLabels->alignableLabel(iAli),
1178 thePedeLabels->alignableTracker()->objectIdProvider().typeToName(iAli->alignableObjectId()));
1179 params->setUserVariables(userVars);
1180 }
1181 }
1182 }
1183
1184
1185 unsigned int MillePedeAlignmentAlgorithm::decodeMode(const std::string &mode) const {
1186 if (mode == "full") {
1187 return myMilleBit + myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1188 } else if (mode == "mille") {
1189 return myMilleBit;
1190 } else if (mode == "pede") {
1191 return myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1192 } else if (mode == "pedeSteer") {
1193 return myPedeSteerBit;
1194 } else if (mode == "pedeRun") {
1195 return myPedeSteerBit + myPedeRunBit + myPedeReadBit;
1196 } else if (mode == "pedeRead") {
1197 return myPedeReadBit;
1198 }
1199
1200 throw cms::Exception("BadConfig") << "Unknown mode '" << mode
1201 << "', use 'full', 'mille', 'pede', 'pedeRun', 'pedeSteer' or 'pedeRead'.";
1202
1203 return 0;
1204 }
1205
1206
1207 bool MillePedeAlignmentAlgorithm::addHitStatistics(int fromIov,
1208 const std::string &outFile,
1209 const std::vector<std::string> &inFiles) const {
1210 bool allOk = true;
1211 int ierr = 0;
1212 MillePedeVariablesIORoot millePedeIO;
1213 for (std::vector<std::string>::const_iterator iFile = inFiles.begin(); iFile != inFiles.end(); ++iFile) {
1214 const std::string inFile(theDir + *iFile);
1215 const std::vector<AlignmentUserVariables *> mpVars =
1216 millePedeIO.readMillePedeVariables(theAlignables, inFile.c_str(), fromIov, ierr);
1217 if (ierr || !this->addHits(theAlignables, mpVars)) {
1218 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addHitStatistics"
1219 << "Error " << ierr << " reading from " << inFile << ", tree " << fromIov
1220 << ", or problems in addHits";
1221 allOk = false;
1222 }
1223 for (std::vector<AlignmentUserVariables *>::const_iterator i = mpVars.begin(); i != mpVars.end(); ++i) {
1224 delete *i;
1225 }
1226 }
1227
1228 return allOk;
1229 }
1230
1231
1232 bool MillePedeAlignmentAlgorithm::addHits(const align::Alignables &alis,
1233 const std::vector<AlignmentUserVariables *> &mpVars) const {
1234 bool allOk = (mpVars.size() == alis.size());
1235 std::vector<AlignmentUserVariables *>::const_iterator iUser = mpVars.begin();
1236 for (auto iAli = alis.cbegin(); iAli != alis.cend() && iUser != mpVars.end(); ++iAli, ++iUser) {
1237 MillePedeVariables *mpVarNew = dynamic_cast<MillePedeVariables *>(*iUser);
1238 AlignmentParameters *ps = (*iAli)->alignmentParameters();
1239 MillePedeVariables *mpVarOld = (ps ? dynamic_cast<MillePedeVariables *>(ps->userVariables()) : nullptr);
1240 if (!mpVarNew || !mpVarOld || mpVarOld->size() != mpVarNew->size()) {
1241 allOk = false;
1242 continue;
1243 }
1244
1245 mpVarOld->increaseHitsX(mpVarNew->hitsX());
1246 mpVarOld->increaseHitsY(mpVarNew->hitsY());
1247 }
1248
1249 return allOk;
1250 }
1251
1252
1253 template <typename GlobalDerivativeMatrix>
1254 void MillePedeAlignmentAlgorithm::makeGlobDerivMatrix(const std::vector<float> &globalDerivativesx,
1255 const std::vector<float> &globalDerivativesy,
1256 Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) {
1257 static_assert(GlobalDerivativeMatrix::RowsAtCompileTime == 2, "global derivative matrix must have two rows");
1258
1259 for (size_t i = 0; i < globalDerivativesx.size(); ++i) {
1260 aGlobalDerivativesM(0, i) = globalDerivativesx[i];
1261 aGlobalDerivativesM(1, i) = globalDerivativesy[i];
1262 }
1263 }
1264
1265
1266 template <typename CovarianceMatrix,
1267 typename LocalDerivativeMatrix,
1268 typename ResidualMatrix,
1269 typename GlobalDerivativeMatrix>
1270 void MillePedeAlignmentAlgorithm::diagonalize(Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1271 Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM,
1272 Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1273 Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) const {
1274 static_assert(std::is_same<typename LocalDerivativeMatrix::Scalar, typename ResidualMatrix::Scalar>::value,
1275 "'aLocalDerivativesM' and 'aHitResidualsM' must have the "
1276 "same underlying scalar type");
1277 static_assert(std::is_same<typename LocalDerivativeMatrix::Scalar, typename GlobalDerivativeMatrix::Scalar>::value,
1278 "'aLocalDerivativesM' and 'aGlobalDerivativesM' must have the "
1279 "same underlying scalar type");
1280
1281 Eigen::SelfAdjointEigenSolver<typename CovarianceMatrix::PlainObject> myDiag{aHitCovarianceM};
1282
1283 auto aTranfoToDiagonalSystemInv =
1284 myDiag.eigenvectors().transpose().template cast<typename LocalDerivativeMatrix::Scalar>();
1285
1286 aHitCovarianceM = myDiag.eigenvalues().asDiagonal();
1287 aLocalDerivativesM = aTranfoToDiagonalSystemInv * aLocalDerivativesM;
1288 aHitResidualsM = aTranfoToDiagonalSystemInv * aHitResidualsM;
1289 if (aGlobalDerivativesM.size() > 0) {
1290
1291 aGlobalDerivativesM = aTranfoToDiagonalSystemInv * aGlobalDerivativesM;
1292 }
1293 }
1294
1295
1296 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1297 void MillePedeAlignmentAlgorithm ::addRefTrackVirtualMeas1D(
1298 const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1299 unsigned int iVirtualMeas,
1300 Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1301 Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1302 Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1303
1304
1305 const unsigned int xIndex = iVirtualMeas + refTrajPtr->numberOfHitMeas();
1306
1307 aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1308 aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex];
1309
1310 const auto &locDerivMatrix = refTrajPtr->derivatives();
1311 for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1312 aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1313 }
1314 }
1315
1316
1317 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1318 void MillePedeAlignmentAlgorithm ::addRefTrackData2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1319 unsigned int iTrajHit,
1320 Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1321 Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1322 Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1323
1324
1325 const unsigned int xIndex = iTrajHit * 2;
1326 const unsigned int yIndex = iTrajHit * 2 + 1;
1327
1328 aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1329 aHitCovarianceM(0, 1) = refTrajPtr->measurementErrors()[xIndex][yIndex];
1330 aHitCovarianceM(1, 0) = refTrajPtr->measurementErrors()[yIndex][xIndex];
1331 aHitCovarianceM(1, 1) = refTrajPtr->measurementErrors()[yIndex][yIndex];
1332
1333 aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1334 aHitResidualsM(1, 0) = refTrajPtr->measurements()[yIndex] - refTrajPtr->trajectoryPositions()[yIndex];
1335
1336 const auto &locDerivMatrix = refTrajPtr->derivatives();
1337 for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1338 aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1339 aLocalDerivativesM(1, i) = locDerivMatrix[yIndex][i];
1340 }
1341 }
1342
1343
1344 int MillePedeAlignmentAlgorithm ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1345 unsigned int iTrajHit,
1346 const std::vector<int> &globalLabels,
1347 const std::vector<float> &globalDerivativesX,
1348 const std::vector<float> &globalDerivativesY) {
1349 const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1350
1351 if ((aRecHit)->dimension() == 1) {
1352 return this->callMille1D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX);
1353 } else {
1354 return this->callMille2D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX, globalDerivativesY);
1355 }
1356 }
1357
1358
1359 int MillePedeAlignmentAlgorithm ::callMille1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1360 unsigned int iTrajHit,
1361 const std::vector<int> &globalLabels,
1362 const std::vector<float> &globalDerivativesX) {
1363 const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1364 const unsigned int xIndex = iTrajHit * 2;
1365
1366
1367 const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1368 const int nLocal = locDerivMatrix.num_col();
1369 std::vector<float> localDerivatives(nLocal);
1370 for (unsigned int i = 0; i < localDerivatives.size(); ++i) {
1371 localDerivatives[i] = locDerivMatrix[xIndex][i];
1372 }
1373
1374
1375 float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1376 float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
1377
1378
1379 const int nGlobal = globalDerivativesX.size();
1380
1381
1382
1383 theMille->mille(
1384 nLocal, &(localDerivatives[0]), nGlobal, &(globalDerivativesX[0]), &(globalLabels[0]), residX, hitErrX);
1385
1386 if (theMonitor) {
1387 theMonitor->fillDerivatives(
1388 aRecHit, &(localDerivatives[0]), nLocal, &(globalDerivativesX[0]), nGlobal, &(globalLabels[0]));
1389 theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, residX, hitErrX, false);
1390 }
1391
1392 return 1;
1393 }
1394
1395
1396 int MillePedeAlignmentAlgorithm ::callMille2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1397 unsigned int iTrajHit,
1398 const std::vector<int> &globalLabels,
1399 const std::vector<float> &globalDerivativesx,
1400 const std::vector<float> &globalDerivativesy) {
1401 const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1402
1403 if ((aRecHit)->dimension() != 2) {
1404 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::callMille2D"
1405 << "You try to call method for 2D hits for a " << (aRecHit)->dimension()
1406 << "D Hit. Hit gets ignored!";
1407 return -1;
1408 }
1409
1410 Eigen::Matrix<double, 2, 2> aHitCovarianceM;
1411 Eigen::Matrix<float, 2, 1> aHitResidualsM;
1412 Eigen::Matrix<float, 2, Eigen::Dynamic> aLocalDerivativesM{2, refTrajPtr->derivatives().num_col()};
1413
1414 this->addRefTrackData2D(refTrajPtr, iTrajHit, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1415 Eigen::Matrix<float, 2, Eigen::Dynamic> aGlobalDerivativesM{2, globalDerivativesx.size()};
1416 this->makeGlobDerivMatrix(globalDerivativesx, globalDerivativesy, aGlobalDerivativesM);
1417
1418
1419
1420
1421 const double corr = aHitCovarianceM(0, 1) / sqrt(aHitCovarianceM(0, 0) * aHitCovarianceM(1, 1));
1422 if (theMonitor)
1423 theMonitor->fillCorrelations2D(corr, aRecHit);
1424 bool diag = false;
1425 switch (aRecHit->geographicalId().subdetId()) {
1426 case SiStripDetId::TID:
1427 case SiStripDetId::TEC:
1428 if (aRecHit->geographicalId().det() == DetId::Tracker && TMath::Abs(corr) > theMaximalCor2D) {
1429 this->diagonalize(aHitCovarianceM, aLocalDerivativesM, aHitResidualsM, aGlobalDerivativesM);
1430 diag = true;
1431 }
1432 break;
1433 default:;
1434 }
1435
1436 float newResidX = aHitResidualsM(0, 0);
1437 float newResidY = aHitResidualsM(1, 0);
1438 float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1439 float newHitErrY = TMath::Sqrt(aHitCovarianceM(1, 1));
1440
1441
1442
1443 std::vector<float> newLocalDerivs(aLocalDerivativesM.size());
1444 Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1445 newLocalDerivs.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1446 float *newLocalDerivsX = &(newLocalDerivs[0]);
1447 float *newLocalDerivsY = &(newLocalDerivs[aLocalDerivativesM.cols()]);
1448
1449
1450
1451 std::vector<float> newGlobDerivs(aGlobalDerivativesM.size());
1452 Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1453 newGlobDerivs.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1454 float *newGlobDerivsX = &(newGlobDerivs[0]);
1455 float *newGlobDerivsY = &(newGlobDerivs[aGlobalDerivativesM.cols()]);
1456
1457 const int nLocal = aLocalDerivativesM.cols();
1458 const int nGlobal = aGlobalDerivativesM.cols();
1459
1460 if (diag && (newHitErrX > newHitErrY)) {
1461
1462 std::swap(newResidX, newResidY);
1463 std::swap(newHitErrX, newHitErrY);
1464 std::swap(newLocalDerivsX, newLocalDerivsY);
1465 std::swap(newGlobDerivsX, newGlobDerivsY);
1466 }
1467
1468
1469
1470 theMille->mille(nLocal, newLocalDerivsX, nGlobal, newGlobDerivsX, &(globalLabels[0]), newResidX, newHitErrX);
1471
1472 if (theMonitor) {
1473 theMonitor->fillDerivatives(aRecHit, newLocalDerivsX, nLocal, newGlobDerivsX, nGlobal, &(globalLabels[0]));
1474 theMonitor->fillResiduals(
1475 aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidX, newHitErrX, false);
1476 }
1477 const bool isReal2DHit = this->is2D(aRecHit);
1478 if (isReal2DHit) {
1479 theMille->mille(nLocal, newLocalDerivsY, nGlobal, newGlobDerivsY, &(globalLabels[0]), newResidY, newHitErrY);
1480 if (theMonitor) {
1481 theMonitor->fillDerivatives(aRecHit, newLocalDerivsY, nLocal, newGlobDerivsY, nGlobal, &(globalLabels[0]));
1482 theMonitor->fillResiduals(
1483 aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidY, newHitErrY, true);
1484 }
1485 }
1486
1487 return (isReal2DHit ? 2 : 1);
1488 }
1489
1490
1491 void MillePedeAlignmentAlgorithm ::addVirtualMeas(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1492 unsigned int iVirtualMeas) {
1493 Eigen::Matrix<double, 1, 1> aHitCovarianceM;
1494 Eigen::Matrix<float, 1, 1> aHitResidualsM;
1495 Eigen::Matrix<float, 1, Eigen::Dynamic> aLocalDerivativesM{1, refTrajPtr->derivatives().num_col()};
1496
1497 this->addRefTrackVirtualMeas1D(refTrajPtr, iVirtualMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1498
1499
1500 auto aGlobalDerivativesM = Eigen::Matrix<float, 1, 1>::Zero();
1501
1502 float newResidX = aHitResidualsM(0, 0);
1503 float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1504 std::vector<float> newLocalDerivsX(aLocalDerivativesM.size());
1505 Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1506 newLocalDerivsX.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1507
1508 std::vector<float> newGlobDerivsX(aGlobalDerivativesM.size());
1509 Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1510 newGlobDerivsX.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1511
1512 const int nLocal = aLocalDerivativesM.cols();
1513 const int nGlobal = 0;
1514
1515 theMille->mille(nLocal, newLocalDerivsX.data(), nGlobal, newGlobDerivsX.data(), &nGlobal, newResidX, newHitErrX);
1516 }
1517
1518
1519 void MillePedeAlignmentAlgorithm::addLaserData(const EventInfo &eventInfo,
1520 const TkFittedLasBeamCollection &lasBeams,
1521 const TsosVectorCollection &lasBeamTsoses) {
1522 TsosVectorCollection::const_iterator iTsoses = lasBeamTsoses.begin();
1523 for (TkFittedLasBeamCollection::const_iterator iBeam = lasBeams.begin(), iEnd = lasBeams.end(); iBeam != iEnd;
1524 ++iBeam, ++iTsoses) {
1525
1526 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addLaserData"
1527 << "Beam " << iBeam->getBeamId() << " with " << iBeam->parameters().size()
1528 << " parameters and " << iBeam->getData().size() << " hits.\n There are "
1529 << iTsoses->size() << " TSOSes.";
1530
1531 this->addLasBeam(eventInfo, *iBeam, *iTsoses);
1532 }
1533 }
1534
1535
1536 void MillePedeAlignmentAlgorithm::addLasBeam(const EventInfo &eventInfo,
1537 const TkFittedLasBeam &lasBeam,
1538 const std::vector<TrajectoryStateOnSurface> &tsoses) {
1539 AlignmentParameters *dummyPtr = nullptr;
1540 std::vector<float> lasLocalDerivsX;
1541 const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId());
1542
1543 for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
1544 if (!tsoses[iHit].isValid())
1545 continue;
1546
1547 theFloatBufferX.clear();
1548 theFloatBufferY.clear();
1549 theIntBuffer.clear();
1550 lasLocalDerivsX.clear();
1551
1552 const SiStripLaserRecHit2D &hit = lasBeam.getData()[iHit];
1553 AlignableDetOrUnitPtr lasAli(theAlignableNavigator->alignableFromDetId(hit.getDetId()));
1554 this->globalDerivativesHierarchy(
1555 eventInfo, tsoses[iHit], lasAli, lasAli, theFloatBufferX, theFloatBufferY, theIntBuffer, dummyPtr);
1556
1557 for (unsigned int nFitParams = 0; nFitParams < static_cast<unsigned int>(lasBeam.parameters().size());
1558 ++nFitParams) {
1559 const float derivative = lasBeam.derivatives()[iHit][nFitParams];
1560 if (nFitParams < lasBeam.firstFixedParameter()) {
1561 lasLocalDerivsX.push_back(derivative);
1562 } else {
1563 const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
1564 theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
1565 theFloatBufferX.push_back(derivative);
1566 }
1567 }
1568
1569 const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
1570
1571 const float error = 0.003;
1572
1573 theMille->mille(lasLocalDerivsX.size(),
1574 &(lasLocalDerivsX[0]),
1575 theFloatBufferX.size(),
1576 &(theFloatBufferX[0]),
1577 &(theIntBuffer[0]),
1578 residual,
1579 error);
1580 }
1581
1582 theMille->end();
1583 }
1584
1585 void MillePedeAlignmentAlgorithm::addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg) {
1586
1587 const bool doOutputOnStdout(pxbSurveyCfg.getParameter<bool>("doOutputOnStdout"));
1588 if (doOutputOnStdout) {
1589 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1590 << "# Output from addPxbSurvey follows below because "
1591 << "doOutputOnStdout is set to True";
1592 }
1593
1594
1595 SurveyPxbDicer dicer(pxbSurveyCfg.getParameter<std::vector<edm::ParameterSet> >("toySurveyParameters"),
1596 pxbSurveyCfg.getParameter<unsigned int>("toySurveySeed"));
1597 std::ofstream outfile(pxbSurveyCfg.getUntrackedParameter<std::string>("toySurveyFile").c_str());
1598
1599
1600 std::vector<SurveyPxbImageLocalFit> measurements;
1601 std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
1602 SurveyPxbImageReader<SurveyPxbImageLocalFit> reader(filename, measurements, 800);
1603
1604
1605 for (std::vector<SurveyPxbImageLocalFit>::size_type i = 0; i != measurements.size(); i++) {
1606 if (doOutputOnStdout) {
1607 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1608 << "Module " << i << ": ";
1609 }
1610
1611
1612 AlignableDetOrUnitPtr mod1(theAlignableNavigator->alignableFromDetId(measurements[i].getIdFirst()));
1613 AlignableDetOrUnitPtr mod2(theAlignableNavigator->alignableFromDetId(measurements[i].getIdSecond()));
1614 const AlignableSurface &surf1 = mod1->surface();
1615 const AlignableSurface &surf2 = mod2->surface();
1616
1617
1618 const LocalPoint fidpoint0(-0.91, +3.30);
1619 const LocalPoint fidpoint1(+0.91, +3.30);
1620 const LocalPoint fidpoint2(+0.91, -3.30);
1621 const LocalPoint fidpoint3(-0.91, -3.30);
1622
1623
1624
1625
1626 const GlobalPoint surf2point0(surf2.toGlobal(fidpoint0));
1627 const GlobalPoint surf2point1(surf2.toGlobal(fidpoint1));
1628 const LocalPoint fidpoint0inSurf1frame(surf1.toLocal(surf2point0));
1629 const LocalPoint fidpoint1inSurf1frame(surf1.toLocal(surf2point1));
1630
1631
1632 SurveyPxbImageLocalFit::fidpoint_t fidpointvec;
1633 fidpointvec.push_back(fidpoint0inSurf1frame);
1634 fidpointvec.push_back(fidpoint1inSurf1frame);
1635 fidpointvec.push_back(fidpoint2);
1636 fidpointvec.push_back(fidpoint3);
1637
1638
1639 if (pxbSurveyCfg.getParameter<bool>("doToySurvey")) {
1640 dicer.doDice(fidpointvec, measurements[i].getIdPair(), outfile);
1641 }
1642
1643
1644 measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
1645 SurveyPxbImageLocalFit::localpars_t a;
1646 a = measurements[i].getLocalParameters();
1647 const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
1648
1649
1650 if (doOutputOnStdout) {
1651 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1652 << "a: " << a[0] << ", " << a[1] << ", " << a[2] << ", " << a[3]
1653 << " S= " << sqrt(a[2] * a[2] + a[3] * a[3]) << " phi= " << atan(a[3] / a[2])
1654 << " chi2= " << chi2 << std::endl;
1655 }
1656 if (theMonitor) {
1657 theMonitor->fillPxbSurveyHistsChi2(chi2);
1658 theMonitor->fillPxbSurveyHistsLocalPars(a[0], a[1], sqrt(a[2] * a[2] + a[3] * a[3]), atan(a[3] / a[2]));
1659 }
1660
1661
1662 for (SurveyPxbImageLocalFit::count_t j = 0; j != SurveyPxbImageLocalFit::nMsrmts; j++) {
1663 theMille->mille((int)measurements[i].getLocalDerivsSize(),
1664 measurements[i].getLocalDerivsPtr(j),
1665 (int)measurements[i].getGlobalDerivsSize(),
1666 measurements[i].getGlobalDerivsPtr(j),
1667 measurements[i].getGlobalDerivsLabelPtr(j),
1668 measurements[i].getResiduum(j),
1669 measurements[i].getSigma(j));
1670 }
1671 theMille->end();
1672 }
1673 outfile.close();
1674 }
1675
1676 bool MillePedeAlignmentAlgorithm::areIOVsSpecified() const {
1677 const auto runRangeSelection = theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection");
1678
1679 if (runRangeSelection.empty())
1680 return false;
1681
1682 const auto runRanges =
1683 align::makeNonOverlappingRunRanges(runRangeSelection, cond::timeTypeSpecs[cond::runnumber].beginValue);
1684
1685 return !(runRanges.empty());
1686 }