CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
MillePedeAlignmentAlgorithm.cc
Go to the documentation of this file.
1 
12 
17 
19 
24 #include "Alignment/MillePedeAlignmentAlgorithm/src/Mille.h" // 'unpublished' interface located in src
29 
32 
36 
40 
41 // includes to make known that they inherit from Alignable:
45 
50 
56 
60 
62 
63 #include <fstream>
64 #include <sstream>
65 #include <algorithm>
66 #include <sys/stat.h>
67 
68 #include <TMath.h>
72 
73 // Includes for PXB survey
80 
82 
83 using namespace gbl;
84 
85 // Constructor ----------------------------------------------------------------
86 //____________________________________________________
89  topoToken_(iC.esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
90  aliThrToken_(iC.esConsumes<AlignPCLThresholdsHG, AlignPCLThresholdsHGRcd, edm::Transition::BeginRun>()),
91  siPixelQualityToken_(iC.esConsumes<SiPixelQuality, SiPixelQualityFromDbRcd, edm::Transition::BeginRun>()),
92  geomToken_(iC.esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
93  theConfig(cfg),
94  theMode(this->decodeMode(theConfig.getUntrackedParameter<std::string>("mode"))),
95  theDir(theConfig.getUntrackedParameter<std::string>("fileDir")),
96  theAlignmentParameterStore(nullptr),
97  theAlignables(),
98  theTrajectoryFactory(
99  TrajectoryFactoryPlugin::get()->create(theConfig.getParameter<edm::ParameterSet>("TrajectoryFactory")
100  .getParameter<std::string>("TrajectoryFactoryName"),
101  theConfig.getParameter<edm::ParameterSet>("TrajectoryFactory"),
102  iC)),
103  theMinNumHits(cfg.getParameter<unsigned int>("minNumHits")),
104  theMaximalCor2D(cfg.getParameter<double>("max2Dcorrelation")),
105  firstIOV_(cfg.getUntrackedParameter<AlignmentAlgorithmBase::RunNumber>("firstIOV")),
106  ignoreFirstIOVCheck_(cfg.getUntrackedParameter<bool>("ignoreFirstIOVCheck")),
107  enableAlignableUpdates_(cfg.getUntrackedParameter<bool>("enableAlignableUpdates")),
108  theLastWrittenIov(0),
109  theGblDoubleBinary(cfg.getParameter<bool>("doubleBinary")),
110  runAtPCL_(cfg.getParameter<bool>("runAtPCL")),
111  ignoreHitsWithoutGlobalDerivatives_(cfg.getParameter<bool>("ignoreHitsWithoutGlobalDerivatives")),
112  skipGlobalPositionRcdCheck_(cfg.getParameter<bool>("skipGlobalPositionRcdCheck")),
113  uniqueRunRanges_(align::makeUniqueRunRanges(cfg.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection"),
114  cond::timeTypeSpecs[cond::runnumber].beginValue)),
115  enforceSingleIOVInput_(!(enableAlignableUpdates_ && areIOVsSpecified())),
116  lastProcessedRun_(cond::timeTypeSpecs[cond::runnumber].beginValue) {
117  pixelQuality = nullptr;
118  if (!theDir.empty() && theDir.find_last_of('/') != theDir.size() - 1)
119  theDir += '/'; // may need '/'
120  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm"
121  << "Start in mode '" << theConfig.getUntrackedParameter<std::string>("mode")
122  << "' with output directory '" << theDir << "'.";
123  if (this->isMode(myMilleBit)) {
124  theMille = std::make_unique<Mille>(
125  (theDir + theConfig.getParameter<std::string>("binaryFile")).c_str()); // add ', false);' for text output);
126  // use same file for GBL
127  theBinary = std::make_unique<MilleBinary>((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str(),
129  }
130 }
131 
132 // Destructor ----------------------------------------------------------------
133 //____________________________________________________
135 
136 // Call at beginning of job ---------------------------------------------------
137 //____________________________________________________
141  AlignableExtras *extras,
142  AlignmentParameterStore *store) {
143  if (muon) {
144  edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
145  << "Running with AlignabeMuon not yet tested.";
146  }
147 
151  const auto &iov_alignments = setup.get<TrackerAlignmentRcd>().validityInterval();
152  const auto &iov_surfaces = setup.get<TrackerSurfaceDeformationRcd>().validityInterval();
153  const auto &iov_errors = setup.get<TrackerAlignmentErrorExtendedRcd>().validityInterval();
154 
155  std::ostringstream message;
156  bool throwException{false};
157  if (iov_alignments.first().eventID().run() != MIN_VAL || iov_alignments.last().eventID().run() != MAX_VAL) {
158  message << "\nTrying to apply " << setup.get<TrackerAlignmentRcd>().key().name()
159  << " with multiple IOVs in tag without specifying 'RunRangeSelection'.\n"
160  << "Validity range is " << iov_alignments.first().eventID().run() << " - "
161  << iov_alignments.last().eventID().run() << "\n";
162  throwException = true;
163  }
164  if (iov_surfaces.first().eventID().run() != MIN_VAL || iov_surfaces.last().eventID().run() != MAX_VAL) {
165  message << "\nTrying to apply " << setup.get<TrackerSurfaceDeformationRcd>().key().name()
166  << " with multiple IOVs in tag without specifying 'RunRangeSelection'.\n"
167  << "Validity range is " << iov_surfaces.first().eventID().run() << " - "
168  << iov_surfaces.last().eventID().run() << "\n";
169  throwException = true;
170  }
171  if (iov_errors.first().eventID().run() != MIN_VAL || iov_errors.last().eventID().run() != MAX_VAL) {
172  message << "\nTrying to apply " << setup.get<TrackerAlignmentErrorExtendedRcd>().key().name()
173  << " with multiple IOVs in tag without specifying 'RunRangeSelection'.\n"
174  << "Validity range is " << iov_errors.first().eventID().run() << " - "
175  << iov_errors.last().eventID().run() << "\n";
176  throwException = true;
177  }
178  if (throwException) {
179  throw cms::Exception("DatabaseError") << "@SUB=MillePedeAlignmentAlgorithm::initialize" << message.str();
180  }
181  }
182 
183  //Retrieve tracker topology from geometry
184  const TrackerTopology *const tTopo = &setup.getData(topoToken_);
185 
186  //Retrieve the thresolds cuts from DB for the PCL
187  if (runAtPCL_) {
188  const auto &th = &setup.getData(aliThrToken_);
189  theThresholds = std::make_shared<AlignPCLThresholdsHG>();
190  storeThresholds(th->getNrecords(), th->getThreshold_Map(), th->getFloatMap());
191 
192  // Retrieve the SiPixelQuality object from setup
193  const SiPixelQuality &qual = setup.getData(siPixelQualityToken_);
194  // Create a new SiPixelQuality object on the heap using the copy constructor
195  pixelQuality = std::make_shared<SiPixelQuality>(qual);
196 
197  //Retrieve tracker geometry
198  const TrackerGeometry *tGeom = &setup.getData(geomToken_);
199  //Retrieve PixelTopologyMap
200  pixelTopologyMap = std::make_shared<PixelTopologyMap>(tGeom, tTopo);
201  }
202 
203  theAlignableNavigator = std::make_unique<AlignableNavigator>(extras, tracker, muon);
206 
207  edm::ParameterSet pedeLabelerCfg(theConfig.getParameter<edm::ParameterSet>("pedeLabeler"));
208  edm::VParameterSet RunRangeSelectionVPSet(theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection"));
209  pedeLabelerCfg.addUntrackedParameter<edm::VParameterSet>("RunRangeSelection", RunRangeSelectionVPSet);
210 
211  std::string labelerPlugin = "PedeLabeler";
212  if (!RunRangeSelectionVPSet.empty()) {
213  labelerPlugin = "RunRangeDependentPedeLabeler";
214  if (pedeLabelerCfg.exists("plugin")) {
215  std::string labelerPluginCfg = pedeLabelerCfg.getParameter<std::string>("plugin");
216  if ((labelerPluginCfg != "PedeLabeler" && labelerPluginCfg != "RunRangeDependentPedeLabeler") ||
217  !pedeLabelerCfg.getUntrackedParameter<edm::VParameterSet>("parameterInstances").empty()) {
218  throw cms::Exception("BadConfig") << "MillePedeAlignmentAlgorithm::initialize"
219  << "both RunRangeSelection and generic labeler specified in config file. "
220  << "Please get rid of either one of them.\n";
221  }
222  }
223  } else {
224  if (pedeLabelerCfg.exists("plugin")) {
225  labelerPlugin = pedeLabelerCfg.getParameter<std::string>("plugin");
226  }
227  }
228 
229  if (!pedeLabelerCfg.exists("plugin")) {
230  pedeLabelerCfg.addUntrackedParameter<std::string>("plugin", labelerPlugin);
231  }
232 
233  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
234  << "Using plugin '" << labelerPlugin << "' to generate labels.";
235 
236  thePedeLabels = std::shared_ptr<PedeLabelerBase>(PedeLabelerPluginFactory::get()->create(
237  labelerPlugin, PedeLabelerBase::TopLevelAlignables(tracker, muon, extras), pedeLabelerCfg));
238 
239  // 1) Create PedeSteerer: correct alignable positions for coordinate system selection
240  edm::ParameterSet pedeSteerCfg(theConfig.getParameter<edm::ParameterSet>("pedeSteerer"));
241  thePedeSteer = std::make_unique<PedeSteerer>(tracker,
242  muon,
243  extras,
245  thePedeLabels.get(),
246  pedeSteerCfg,
247  theDir,
248  !this->isMode(myPedeSteerBit));
249 
250  // 2) If requested, directly read in and apply result of previous pede run,
251  // assuming that correction from 1) was also applied to create the result:
252  const std::vector<edm::ParameterSet> mprespset(
253  theConfig.getParameter<std::vector<edm::ParameterSet> >("pedeReaderInputs"));
254  if (!mprespset.empty()) {
255  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
256  << "Apply " << mprespset.end() - mprespset.begin()
257  << " previous MillePede constants from 'pedeReaderInputs'.";
258  }
259 
260  // FIXME: add selection of run range via 'pedeReaderInputs'
261  // Note: Results for parameters of IntegratedCalibration's cannot be treated...
263  for (std::vector<edm::ParameterSet>::const_iterator iSet = mprespset.begin(), iE = mprespset.end(); iSet != iE;
264  ++iSet) {
265  // This read will ignore calibrations as long as they are not yet passed to Millepede
266  // during/before initialize(..) - currently addCalibrations(..) is called later in AlignmentProducer
267  if (!this->readFromPede((*iSet), false, runrange)) { // false: do not erase SelectionUserVariables
268  throw cms::Exception("BadConfig")
269  << "MillePedeAlignmentAlgorithm::initialize: Problems reading input constants of "
270  << "pedeReaderInputs entry " << iSet - mprespset.begin() << '.';
271  }
273  // Needed to shut up later warning from checkAliParams:
275  }
276 
277  // 3) Now create steerings with 'final' start position:
278  thePedeSteer->buildSubSteer(tracker, muon, extras);
279 
280  // After (!) 1-3 of PedeSteerer which uses the SelectionUserVariables attached to the parameters:
281  this->buildUserVariables(theAlignables); // for hit statistics and/or pede result
282 
283  if (this->isMode(myMilleBit)) {
284  if (!theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles").empty() ||
285  !theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles").empty()) {
286  throw cms::Exception("BadConfig") << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' must be empty for "
287  << "modes running mille.";
288  }
289  const std::string moniFile(theConfig.getUntrackedParameter<std::string>("monitorFile"));
290  if (!moniFile.empty())
291  theMonitor = std::make_unique<MillePedeMonitor>(tTopo, (theDir + moniFile).c_str());
292 
293  // Get trajectory factory. In case nothing found, FrameWork will throw...
294  }
295 
296  if (this->isMode(myPedeSteerBit)) {
297  // Get config for survey and set flag accordingly
298  const edm::ParameterSet pxbSurveyCfg(theConfig.getParameter<edm::ParameterSet>("surveyPixelBarrel"));
299  theDoSurveyPixelBarrel = pxbSurveyCfg.getParameter<bool>("doSurvey");
301  this->addPxbSurvey(pxbSurveyCfg);
302  }
303 }
304 
305 //____________________________________________________
307 
308 //____________________________________________________
309 bool MillePedeAlignmentAlgorithm::addCalibrations(const std::vector<IntegratedCalibrationBase *> &iCals) {
310  theCalibrations.insert(theCalibrations.end(), iCals.begin(), iCals.end());
311  thePedeLabels->addCalibrations(iCals);
312  return true;
313 }
314 
315 //____________________________________________________
317  const AlignPCLThresholdsHG::threshold_map &thresholdMap,
318  const AlignPCLThresholdsHG::param_map &floatMap) {
319  theThresholds->setAlignPCLThresholds(nRecords, thresholdMap);
320  theThresholds->setFloatMap(floatMap);
321  return true;
322 }
323 
324 //_____________________________________________________________________________
326  if (isMode(myMilleBit)) {
327  return true;
328  } else {
329  return false;
330  }
331 }
332 
333 //_____________________________________________________________________________
335  if (isMode(myPedeReadBit)) {
336  if (runAtPCL_) {
337  MillePedeFileReader mpReader(theConfig.getParameter<edm::ParameterSet>("MillePedeFileReader"),
341  pixelQuality);
342  mpReader.read();
343  return mpReader.storeAlignments();
344  } else {
345  return true;
346  }
347  } else {
348  return false;
349  }
350 }
351 
352 //____________________________________________________
354  if (this->isMode(myPedeReadBit)) {
355  if (not theAlignmentParameterStore) {
356  return false;
357  }
358  // restore initial positions, rotations and deformations
361  } else {
363  }
364 
365  // Needed to shut up later warning from checkAliParams:
367  // To avoid that they keep values from previous IOV if no new one in pede result
369 
370  if (!this->readFromPede(theConfig.getParameter<edm::ParameterSet>("pedeReader"), true, runrange)) {
371  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::setParametersForRunRange"
372  << "Problems reading pede result, but applying!";
373  }
375 
376  this->doIO(++theLastWrittenIov); // pre-increment!
377  }
378 
379  return true;
380 }
381 
382 // Call at end of job ---------------------------------------------------------
383 //____________________________________________________
386  theMille.reset(); // delete to close binary before running pede below (flush would be enough...)
387  theBinary.reset();
388 
389  std::vector<std::string> files;
390  if (this->isMode(myMilleBit) || !theConfig.getParameter<std::string>("binaryFile").empty()) {
391  files.push_back(theDir + theConfig.getParameter<std::string>("binaryFile"));
392  } else {
393  const std::vector<std::string> plainFiles(theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
394  files = getExistingFormattedFiles(plainFiles, theDir);
395  // Do some logging:
396  std::string filesForLogOutput;
397  for (const auto &file : files)
398  filesForLogOutput += " " + file + ",";
399  if (filesForLogOutput.length() != 0)
400  filesForLogOutput.pop_back();
401  edm::LogInfo("Alignment") << "Based on the config parameter mergeBinaryFiles, using the following "
402  << "files as input (assigned weights are indicated by ' -- <weight>'):"
403  << filesForLogOutput;
404  }
405 
407  return;
408 
409  // cache all positions, rotations and deformations
412  if (lastProcessedRun_ < uniqueRunRanges_.back().first) {
413  throw cms::Exception("BadConfig") << "@SUB=MillePedeAlignmentAlgorithm::terminate\n"
414  << "Last IOV of 'RunRangeSelection' has not been processed. "
415  << "Please reconfigure your source to process the runs at least up to "
416  << uniqueRunRanges_.back().first << ".";
417  }
418  auto lastCachedRun = uniqueRunRanges_.front().first;
419  for (const auto &runRange : uniqueRunRanges_) {
420  const auto run = runRange.first;
421  if (std::find(cachedRuns_.begin(), cachedRuns_.end(), run) == cachedRuns_.end()) {
424  } else {
425  lastCachedRun = run;
426  }
427  }
429  }
430 
431  const std::string masterSteer(thePedeSteer->buildMasterSteer(files)); // do only if myPedeSteerBit?
432  if (this->isMode(myPedeRunBit)) {
433  thePedeSteer->runPede(masterSteer);
434  }
435 
436  // parameters from pede are not yet applied,
437  // so we can still write start positions (but with hit statistics in case of mille):
438  this->doIO(0);
439  theLastWrittenIov = 0;
440 }
441 
443  const std::vector<std::string> &plainFiles, const std::string &theDir) {
444  std::vector<std::string> files;
445  for (const auto &plainFile : plainFiles) {
446  const std::string &theInputFileName = plainFile;
447  int theNumber = 0;
448  while (true) {
449  // Create a formatted version of the filename, with growing numbers
450  // If the parameter doesn't contain a formatting directive, it just stays unchanged
451  char theNumberedInputFileName[200];
452  sprintf(theNumberedInputFileName, theInputFileName.c_str(), theNumber);
453  std::string theCompleteInputFileName = theDir + theNumberedInputFileName;
454  const auto endOfStrippedFileName = theCompleteInputFileName.rfind(" --");
455  const auto strippedInputFileName = theCompleteInputFileName.substr(0, endOfStrippedFileName);
456  // Check if the file exists
457  struct stat buffer;
458  if (stat(strippedInputFileName.c_str(), &buffer) == 0) {
459  // If the file exists, add it to the list
460  files.push_back(theCompleteInputFileName);
461  if (theNumberedInputFileName == theInputFileName) {
462  // If the filename didn't contain a formatting directive, no reason to look any further, break out of the loop
463  break;
464  } else {
465  // Otherwise look for the next number
466  theNumber++;
467  }
468  } else {
469  // The file doesn't exist, break out of the loop
470  break;
471  }
472  }
473  // warning if unformatted (-> theNumber stays at 0) does not exist
474  if (theNumber == 0 && (files.empty() || files.back() != plainFile)) {
475  edm::LogWarning("Alignment") << "The input file '" << plainFile << "' does not exist.";
476  }
477  }
478  return files;
479 }
480 
481 // Run the algorithm on trajectories and tracks -------------------------------
482 //____________________________________________________
484  if (!this->isMode(myMilleBit))
485  return; // no theMille created...
486  const auto &tracks = eventInfo.trajTrackPairs();
487 
488  if (theMonitor) { // monitor input tracks
489  for (const auto &iTrajTrack : tracks) {
490  theMonitor->fillTrack(iTrajTrack.second);
491  }
492  }
493 
494  const RefTrajColl trajectories(theTrajectoryFactory->trajectories(setup, tracks, eventInfo.beamSpot()));
495 
496  // Now loop over ReferenceTrajectoryCollection
497  unsigned int refTrajCount = 0; // counter for track monitoring
498  const auto tracksPerTraj = theTrajectoryFactory->tracksPerTrajectory();
499  for (auto iRefTraj = trajectories.cbegin(), iRefTrajE = trajectories.cend(); iRefTraj != iRefTrajE;
500  ++iRefTraj, ++refTrajCount) {
501  const RefTrajColl::value_type &refTrajPtr = *iRefTraj;
502  if (theMonitor)
503  theMonitor->fillRefTrajectory(refTrajPtr);
504 
505  const auto nHitXy = this->addReferenceTrajectory(setup, eventInfo, refTrajPtr);
506 
507  if (theMonitor && (nHitXy.first || nHitXy.second)) {
508  // if trajectory used (i.e. some hits), fill monitoring
509  const auto offset = tracksPerTraj * refTrajCount;
510  for (unsigned int iTrack = 0; iTrack < tracksPerTraj; ++iTrack) {
511  theMonitor->fillUsedTrack(tracks[offset + iTrack].second, nHitXy.first, nHitXy.second);
512  }
513  }
514 
515  } // end of reference trajectory and track loop
516 }
517 
518 //____________________________________________________
519 std::pair<unsigned int, unsigned int> MillePedeAlignmentAlgorithm::addReferenceTrajectory(
520  const edm::EventSetup &setup, const EventInfo &eventInfo, const RefTrajColl::value_type &refTrajPtr) {
521  std::pair<unsigned int, unsigned int> hitResultXy(0, 0);
522  if (refTrajPtr->isValid()) {
523  // GblTrajectory?
524  if (!refTrajPtr->gblInput().empty()) {
525  // by construction: number of GblPoints == number of recHits or == zero !!!
526  unsigned int iHit = 0;
527  unsigned int numPointsWithMeas = 0;
528  std::vector<GblPoint>::iterator itPoint;
529  auto theGblInput = refTrajPtr->gblInput();
530  for (unsigned int iTraj = 0; iTraj < refTrajPtr->gblInput().size(); ++iTraj) {
531  for (itPoint = refTrajPtr->gblInput()[iTraj].first.begin(); itPoint < refTrajPtr->gblInput()[iTraj].first.end();
532  ++itPoint) {
533  if (this->addGlobalData(setup, eventInfo, refTrajPtr, iHit++, *itPoint) < 0)
534  return hitResultXy;
535  if (itPoint->numMeasurements() >= 1)
536  ++numPointsWithMeas;
537  }
538  }
539  hitResultXy.first = numPointsWithMeas;
540  // check #hits criterion
541  if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits)
542  return hitResultXy;
543  // construct GBL trajectory
544  if (refTrajPtr->gblInput().size() == 1) {
545  // from single track
546  GblTrajectory aGblTrajectory(refTrajPtr->gblInput()[0].first, refTrajPtr->nominalField() != 0);
547  // GBL fit trajectory
548  /*double Chi2;
549  int Ndf;
550  double lostWeight;
551  aGblTrajectory.fit(Chi2, Ndf, lostWeight);
552  std::cout << " GblFit: " << Chi2 << ", " << Ndf << ", " << lostWeight << std::endl; */
553  // write to MP binary file
554  if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
555  aGblTrajectory.milleOut(*theBinary);
556  }
557  if (refTrajPtr->gblInput().size() == 2) {
558  // from TwoBodyDecay
559  GblTrajectory aGblTrajectory(refTrajPtr->gblInput(),
560  refTrajPtr->gblExtDerivatives(),
561  refTrajPtr->gblExtMeasurements(),
562  refTrajPtr->gblExtPrecisions());
563  // write to MP binary file
564  if (aGblTrajectory.isValid() && aGblTrajectory.getNumPoints() >= theMinNumHits)
565  aGblTrajectory.milleOut(*theBinary);
566  }
567  } else {
568  // to add hits if all fine:
569  std::vector<AlignmentParameters *> parVec(refTrajPtr->recHits().size());
570  // collect hit statistics, assuming that there are no y-only hits
571  std::vector<bool> validHitVecY(refTrajPtr->recHits().size(), false);
572  // Use recHits from ReferenceTrajectory (since they have the right order!):
573  for (unsigned int iHit = 0; iHit < refTrajPtr->recHits().size(); ++iHit) {
574  const int flagXY = this->addMeasurementData(setup, eventInfo, refTrajPtr, iHit, parVec[iHit]);
575 
576  if (flagXY < 0) { // problem
577  hitResultXy.first = 0;
578  break;
579  } else { // hit is fine, increase x/y statistics
580  if (flagXY >= 1)
581  ++hitResultXy.first;
582  validHitVecY[iHit] = (flagXY >= 2);
583  }
584  } // end loop on hits
585 
586  // add virtual measurements
587  for (unsigned int iVirtualMeas = 0; iVirtualMeas < refTrajPtr->numberOfVirtualMeas(); ++iVirtualMeas) {
588  this->addVirtualMeas(refTrajPtr, iVirtualMeas);
589  }
590 
591  // kill or end 'track' for mille, depends on #hits criterion
592  if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits) {
593  theMille->kill();
594  hitResultXy.first = hitResultXy.second = 0; //reset
595  } else {
596  theMille->end();
597  // add x/y hit count to MillePedeVariables of parVec,
598  // returning number of y-hits of the reference trajectory
599  hitResultXy.second = this->addHitCount(parVec, validHitVecY);
600  //
601  }
602  }
603 
604  } // end if valid trajectory
605 
606  return hitResultXy;
607 }
608 
609 //____________________________________________________
610 unsigned int MillePedeAlignmentAlgorithm::addHitCount(const std::vector<AlignmentParameters *> &parVec,
611  const std::vector<bool> &validHitVecY) const {
612  // Loop on all hit information in the input arrays and count valid y-hits:
613  unsigned int nHitY = 0;
614  for (unsigned int iHit = 0; iHit < validHitVecY.size(); ++iHit) {
615  Alignable *ali = (parVec[iHit] ? parVec[iHit]->alignable() : nullptr);
616  // Loop upwards on hierarchy of alignables to add hits to all levels
617  // that are currently aligned. If only a non-selected alignable was hit,
618  // (i.e. flagXY == 0 in addReferenceTrajectory(..)), there is no loop at all...
619  while (ali) {
621  if (pars) { // otherwise hierarchy level not selected
622  // cast ensured by previous checks:
623  MillePedeVariables *mpVar = static_cast<MillePedeVariables *>(pars->userVariables());
624  // every hit has an x-measurement, cf. addReferenceTrajectory(..):
625  mpVar->increaseHitsX();
626  if (validHitVecY[iHit]) {
627  mpVar->increaseHitsY();
628  if (pars == parVec[iHit])
629  ++nHitY; // do not count hits twice
630  }
631  }
632  ali = ali->mother();
633  }
634  }
635 
636  return nHitY;
637 }
638 
640  if (run.run() < firstIOV_ && !ignoreFirstIOVCheck_) {
641  throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::beginRun\n"
642  << "Using data (run = " << run.run() << ") prior to the first defined IOV ("
643  << firstIOV_ << ").";
644  }
645 
646  lastProcessedRun_ = run.run();
647 
648  if (changed && enableAlignableUpdates_) {
649  const auto runNumber = run.run();
651  for (auto runRange = uniqueRunRanges_.crbegin(); runRange != uniqueRunRanges_.crend(); ++runRange) {
652  if (runNumber >= runRange->first) {
653  firstRun = runRange->first;
654  break;
655  }
656  }
657  if (std::find(cachedRuns_.begin(), cachedRuns_.end(), firstRun) != cachedRuns_.end()) {
658  const auto &geometryRcd = setup.get<IdealGeometryRecord>();
659  const auto &globalPosRcd = setup.get<GlobalPositionRcd>();
660  const auto &alignmentRcd = setup.get<TrackerAlignmentRcd>();
661  const auto &surfaceRcd = setup.get<TrackerSurfaceDeformationRcd>();
662  const auto &errorRcd = setup.get<TrackerAlignmentErrorExtendedRcd>();
663 
664  std::ostringstream message;
665  bool throwException{false};
666  message << "Trying to cache tracker alignment payloads for a run (" << runNumber << ") in an IOV (" << firstRun
667  << ") that was already cached.\n"
668  << "The following records in your input database tag have an IOV "
669  << "boundary that does not match your IOV definition:\n";
670  if (geometryRcd.validityInterval().first().eventID().run() > firstRun) {
671  message << " - IdealGeometryRecord '" << geometryRcd.key().name() << "' (since "
672  << geometryRcd.validityInterval().first().eventID().run() << ")\n";
673  throwException = true;
674  }
675  if (globalPosRcd.validityInterval().first().eventID().run() > firstRun) {
676  message << " - GlobalPositionRecord '" << globalPosRcd.key().name() << "' (since "
677  << globalPosRcd.validityInterval().first().eventID().run() << ")";
679  message << " --> ignored\n";
680  } else {
681  message << "\n";
682  throwException = true;
683  }
684  }
685  if (alignmentRcd.validityInterval().first().eventID().run() > firstRun) {
686  message << " - TrackerAlignmentRcd '" << alignmentRcd.key().name() << "' (since "
687  << alignmentRcd.validityInterval().first().eventID().run() << ")\n";
688  throwException = true;
689  }
690  if (surfaceRcd.validityInterval().first().eventID().run() > firstRun) {
691  message << " - TrackerSurfaceDeformationRcd '" << surfaceRcd.key().name() << "' (since "
692  << surfaceRcd.validityInterval().first().eventID().run() << ")\n";
693  throwException = true;
694  }
695  if (errorRcd.validityInterval().first().eventID().run() > firstRun) {
696  message << " - TrackerAlignmentErrorExtendedRcd '" << errorRcd.key().name() << "' (since "
697  << errorRcd.validityInterval().first().eventID().run() << ")\n";
698  throwException = true;
699  }
700 
701  if (throwException) {
702  throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::beginRun\n" << message.str();
703  }
704  } else {
705  cachedRuns_.push_back(firstRun);
707  }
708  }
709 }
710 
711 //____________________________________________________
713  const EndRunInfo &runInfo,
714  const edm::EventSetup &setup) {
715  if (runInfo.tkLasBeams() && runInfo.tkLasBeamTsoses()) {
716  // LAS beam treatment
717  this->addLaserData(eventInfo, *(runInfo.tkLasBeams()), *(runInfo.tkLasBeamTsoses()));
718  }
719  if (this->isMode(myMilleBit))
720  theMille->flushOutputFile();
721 }
722 
723 // Implementation of endRun that DOES get called. (Because we need it.)
725  if (this->isMode(myMilleBit))
726  theMille->flushOutputFile();
727 }
728 
729 //____________________________________________________
731  if (!runAtPCL_)
732  return;
733  if (this->isMode(myMilleBit)) {
734  theMille->resetOutputFile();
735  theBinary.reset(); // GBL output has to be considered since same binary file is used
736  theBinary = std::make_unique<MilleBinary>((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str(),
738  }
739 }
740 
741 //____________________________________________________
743  if (!runAtPCL_)
744  return;
745  if (this->isMode(myMilleBit))
746  theMille->flushOutputFile();
747 }
748 
749 //____________________________________________________
751  const EventInfo &eventInfo,
753  unsigned int iHit,
755  params = nullptr;
756  theFloatBufferX.clear();
757  theFloatBufferY.clear();
758  theIntBuffer.clear();
759 
760  const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
761  const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
762  // ignore invalid hits
763  if (!recHitPtr->isValid())
764  return 0;
765 
766  // First add the derivatives from IntegratedCalibration's,
767  // should even be OK if problems for "usual" derivatives from Alignables
768  this->globalDerivativesCalibration(recHitPtr,
769  tsos,
770  setup,
771  eventInfo, // input
774  theIntBuffer); // output
775 
776  // get AlignableDet/Unit for this hit
777  AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
778 
779  if (!this->globalDerivativesHierarchy(eventInfo,
780  tsos,
781  alidet,
782  alidet,
783  theFloatBufferX, // 2x alidet, sic!
785  theIntBuffer,
786  params)) {
787  return -1; // problem
789  return 0; // empty for X: no alignable for hit, nor calibrations
790  } else {
791  // store measurement even if no alignable or calibrations
792  // -> measurement used for pede-internal track-fit
793  return this->callMille(refTrajPtr, iHit, theIntBuffer, theFloatBufferX, theFloatBufferY);
794  }
795 }
796 
797 //____________________________________________________
798 
800  const EventInfo &eventInfo,
802  unsigned int iHit,
803  GblPoint &gblPoint) {
804  AlignmentParameters *params = nullptr;
805  std::vector<double> theDoubleBufferX, theDoubleBufferY;
806  theDoubleBufferX.clear();
807  theDoubleBufferY.clear();
808  theIntBuffer.clear();
809  int iret = 0;
810 
811  const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
812  const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
813  // ignore invalid hits
814  if (!recHitPtr->isValid())
815  return 0;
816 
817  // get AlignableDet/Unit for this hit
818  AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
819 
820  if (!this->globalDerivativesHierarchy(eventInfo,
821  tsos,
822  alidet,
823  alidet,
824  theDoubleBufferX, // 2x alidet, sic!
825  theDoubleBufferY,
826  theIntBuffer,
827  params)) {
828  return -1; // problem
829  }
830  //calibration parameters
831  int globalLabel;
832  std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
833  for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
834  // get all derivatives of this calibration // const unsigned int num =
835  (*iCalib)->derivatives(derivs, *recHitPtr, tsos, setup, eventInfo);
836  for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
837  // transfer label and x/y derivatives
838  globalLabel = thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second);
839  if (globalLabel > 0 && globalLabel <= 2147483647) {
840  theIntBuffer.push_back(globalLabel);
841  theDoubleBufferX.push_back(iValuesInd->first.first);
842  theDoubleBufferY.push_back(iValuesInd->first.second);
843  } else {
844  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addGlobalData"
845  << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
846  }
847  }
848  }
849  unsigned int numGlobals = theIntBuffer.size();
850  if (numGlobals > 0) {
851  Eigen::Matrix<double, 2, Eigen::Dynamic> globalDer{2, numGlobals};
852  for (unsigned int i = 0; i < numGlobals; ++i) {
853  globalDer(0, i) = theDoubleBufferX[i];
854  globalDer(1, i) = theDoubleBufferY[i];
855  }
856  gblPoint.addGlobals(theIntBuffer, globalDer);
857  iret = 1;
858  }
859  return iret;
860 }
861 
862 //____________________________________________________
864  const TrajectoryStateOnSurface &tsos,
865  Alignable *ali,
866  const AlignableDetOrUnitPtr &alidet,
867  std::vector<float> &globalDerivativesX,
868  std::vector<float> &globalDerivativesY,
869  std::vector<int> &globalLabels,
870  AlignmentParameters *&lowestParams) const {
871  // derivatives and labels are recursively attached
872  if (!ali)
873  return true; // no mother might be OK
874 
875  if (false && theMonitor && alidet != ali)
876  theMonitor->fillFrameToFrame(alidet, ali);
877 
879 
880  if (params) {
881  if (!lowestParams)
882  lowestParams = params; // set parameters of lowest level
883 
884  bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
885  const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
886 
887  if (0 == alignableLabel) { // FIXME: what about regardAllHits in Markus' code?
888  edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
889  << "Label not found, skip Alignable.";
890  return false;
891  }
892 
893  const std::vector<bool> &selPars = params->selector();
894  const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
895 
896  // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
897  for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
898  if (selPars[iSel]) {
899  globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
900  if (hasSplitParameters == true) {
901  globalLabels.push_back(thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos));
902  } else {
903  globalLabels.push_back(thePedeLabels->parameterLabel(alignableLabel, iSel));
904  }
905  globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
906  }
907  }
908  // Exclude mothers if Alignable selected to be no part of a hierarchy:
909  if (thePedeSteer->isNoHiera(ali))
910  return true;
911  }
912  // Call recursively for mother, will stop if mother == 0:
913  return this->globalDerivativesHierarchy(
914  eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
915 }
916 
917 //____________________________________________________
919  const TrajectoryStateOnSurface &tsos,
920  Alignable *ali,
921  const AlignableDetOrUnitPtr &alidet,
922  std::vector<double> &globalDerivativesX,
923  std::vector<double> &globalDerivativesY,
924  std::vector<int> &globalLabels,
925  AlignmentParameters *&lowestParams) const {
926  // derivatives and labels are recursively attached
927  if (!ali)
928  return true; // no mother might be OK
929 
930  if (false && theMonitor && alidet != ali)
931  theMonitor->fillFrameToFrame(alidet, ali);
932 
934 
935  if (params) {
936  if (!lowestParams)
937  lowestParams = params; // set parameters of lowest level
938 
939  bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
940  const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
941 
942  if (0 == alignableLabel) { // FIXME: what about regardAllHits in Markus' code?
943  edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
944  << "Label not found, skip Alignable.";
945  return false;
946  }
947 
948  const std::vector<bool> &selPars = params->selector();
949  const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
950  int globalLabel;
951 
952  // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
953  for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
954  if (selPars[iSel]) {
955  if (hasSplitParameters == true) {
956  globalLabel = thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos);
957  } else {
958  globalLabel = thePedeLabels->parameterLabel(alignableLabel, iSel);
959  }
960  if (globalLabel > 0 && globalLabel <= 2147483647) {
961  globalLabels.push_back(globalLabel);
962  globalDerivativesX.push_back(derivs[iSel][kLocalX] / thePedeSteer->cmsToPedeFactor(iSel));
963  globalDerivativesY.push_back(derivs[iSel][kLocalY] / thePedeSteer->cmsToPedeFactor(iSel));
964  } else {
965  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
966  << "Invalid label " << globalLabel << " <= 0 or > 2147483647";
967  }
968  }
969  }
970  // Exclude mothers if Alignable selected to be no part of a hierarchy:
971  if (thePedeSteer->isNoHiera(ali))
972  return true;
973  }
974  // Call recursively for mother, will stop if mother == 0:
975  return this->globalDerivativesHierarchy(
976  eventInfo, tsos, ali->mother(), alidet, globalDerivativesX, globalDerivativesY, globalLabels, lowestParams);
977 }
978 
979 //____________________________________________________
981  const TrajectoryStateOnSurface &tsos,
982  const edm::EventSetup &setup,
983  const EventInfo &eventInfo,
984  std::vector<float> &globalDerivativesX,
985  std::vector<float> &globalDerivativesY,
986  std::vector<int> &globalLabels) const {
987  std::vector<IntegratedCalibrationBase::ValuesIndexPair> derivs;
988  for (auto iCalib = theCalibrations.begin(); iCalib != theCalibrations.end(); ++iCalib) {
989  // get all derivatives of this calibration // const unsigned int num =
990  (*iCalib)->derivatives(derivs, *recHit, tsos, setup, eventInfo);
991  for (auto iValuesInd = derivs.begin(); iValuesInd != derivs.end(); ++iValuesInd) {
992  // transfer label and x/y derivatives
993  globalLabels.push_back(thePedeLabels->calibrationLabel(*iCalib, iValuesInd->second));
994  globalDerivativesX.push_back(iValuesInd->first.first);
995  globalDerivativesY.push_back(iValuesInd->first.second);
996  }
997  }
998 }
999 
1000 // //____________________________________________________
1001 // void MillePedeAlignmentAlgorithm
1002 // ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
1003 // unsigned int iTrajHit, MeasurementDirection xOrY,
1004 // const std::vector<float> &globalDerivatives, const std::vector<int> &globalLabels)
1005 // {
1006 // const unsigned int xyIndex = iTrajHit*2 + xOrY;
1007 // // FIXME: here for residuum and sigma we could use KALMAN-Filter results
1008 // const float residuum =
1009 // refTrajPtr->measurements()[xyIndex] - refTrajPtr->trajectoryPositions()[xyIndex];
1010 // const float covariance = refTrajPtr->measurementErrors()[xyIndex][xyIndex];
1011 // const float sigma = (covariance > 0. ? TMath::Sqrt(covariance) : 0.);
1012 
1013 // const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1014 
1015 // std::vector<float> localDerivs(locDerivMatrix.num_col());
1016 // for (unsigned int i = 0; i < localDerivs.size(); ++i) {
1017 // localDerivs[i] = locDerivMatrix[xyIndex][i];
1018 // }
1019 
1020 // // &(vector[0]) is valid - as long as vector is not empty
1021 // // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1022 // theMille->mille(localDerivs.size(), &(localDerivs[0]),
1023 // globalDerivatives.size(), &(globalDerivatives[0]), &(globalLabels[0]),
1024 // residuum, sigma);
1025 // if (theMonitor) {
1026 // theMonitor->fillDerivatives(refTrajPtr->recHits()[iTrajHit],localDerivs, globalDerivatives,
1027 // (xOrY == kLocalY));
1028 // theMonitor->fillResiduals(refTrajPtr->recHits()[iTrajHit],
1029 // refTrajPtr->trajectoryStates()[iTrajHit],
1030 // iTrajHit, residuum, sigma, (xOrY == kLocalY));
1031 // }
1032 // }
1033 
1034 //____________________________________________________
1036  // FIXME: Check whether this is a reliable and recommended way to find out...
1037 
1038  if (recHit->dimension() < 2) {
1039  return false; // some muon and TIB/TOB stuff really has RecHit1D
1040  } else if (recHit->detUnit()) { // detunit in strip is 1D, in pixel 2D
1041  return recHit->detUnit()->type().isTrackerPixel();
1042  } else { // stereo strips (FIXME: endcap trouble due to non-parallel strips (wedge sensors)?)
1043  if (dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit->hit())) { // check persistent hit
1044  // projected: 1D measurement on 'glued' module
1045  return false;
1046  } else {
1047  return true;
1048  }
1049  }
1050 }
1051 
1052 //__________________________________________________________________________________________________
1054  bool setUserVars,
1055  const RunRange &runrange) {
1056  bool allEmpty = this->areEmptyParams(theAlignables);
1057 
1058  PedeReader reader(mprespset, *thePedeSteer, *thePedeLabels, runrange);
1059  align::Alignables alis;
1060  bool okRead = reader.read(alis, setUserVars); // also may set params of IntegratedCalibration's
1061  bool numMatch = true;
1062 
1063  std::stringstream out;
1064  out << "Read " << alis.size() << " alignables";
1065  if (alis.size() != theAlignables.size()) {
1066  out << " while " << theAlignables.size() << " in store";
1067  numMatch = false; // FIXME: Should we check one by one? Or transfer 'alis' to the store?
1068  }
1069  if (!okRead)
1070  out << ", but problems in reading";
1071  if (!allEmpty)
1072  out << ", possibly overwriting previous settings";
1073  out << ".";
1074 
1075  if (okRead && allEmpty) {
1076  if (numMatch) { // as many alignables with result as trying to align
1077  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1078  } else if (!alis.empty()) { // dead module do not get hits and no pede result
1079  edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1080  } else { // serious problem: no result read - and not all modules can be dead...
1081  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1082  return false;
1083  }
1084  return true;
1085  }
1086  // the rest is not OK:
1087  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
1088  return false;
1089 }
1090 
1091 //__________________________________________________________________________________________________
1093  for (const auto &iAli : alignables) {
1094  const AlignmentParameters *params = iAli->alignmentParameters();
1095  if (params) {
1096  const auto &parVec(params->parameters());
1097  const auto &parCov(params->covariance());
1098  for (int i = 0; i < parVec.num_row(); ++i) {
1099  if (parVec[i] != 0.)
1100  return false;
1101  for (int j = i; j < parCov.num_col(); ++j) {
1102  if (parCov[i][j] != 0.)
1103  return false;
1104  }
1105  }
1106  }
1107  }
1108 
1109  return true;
1110 }
1111 
1112 //__________________________________________________________________________________________________
1113 unsigned int MillePedeAlignmentAlgorithm::doIO(int loop) const {
1114  unsigned int result = 0;
1115 
1116  const std::string outFilePlain(theConfig.getParameter<std::string>("treeFile"));
1117  if (outFilePlain.empty()) {
1118  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1119  << "treeFile parameter empty => skip writing for 'loop' " << loop;
1120  return result;
1121  }
1122 
1123  const std::string outFile(theDir + outFilePlain);
1124 
1125  AlignmentIORoot aliIO;
1126  int ioerr = 0;
1127  if (loop == 0) {
1128  aliIO.writeAlignableOriginalPositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1129  if (ioerr) {
1130  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1131  << "Problem " << ioerr << " in writeAlignableOriginalPositions";
1132  ++result;
1133  }
1134  } else if (loop == 1) {
1135  // only for first iov add hit counts, else 2x, 3x,... number of hits in IOV 2, 3,...
1136  const std::vector<std::string> inFiles(theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles"));
1137  const std::vector<std::string> binFiles(theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
1138  if (inFiles.size() != binFiles.size()) {
1139  edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1140  << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' "
1141  << "differ in size";
1142  }
1143  this->addHitStatistics(0, outFile, inFiles); // add hit info from tree 0 in 'infiles'
1144  }
1145  MillePedeVariablesIORoot millePedeIO;
1146  millePedeIO.writeMillePedeVariables(theAlignables, outFile.c_str(), loop, false, ioerr);
1147  if (ioerr) {
1148  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1149  << "Problem " << ioerr << " writing MillePedeVariables";
1150  ++result;
1151  }
1152 
1153  aliIO.writeOrigRigidBodyAlignmentParameters(theAlignables, outFile.c_str(), loop, false, ioerr);
1154  if (ioerr) {
1155  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1156  << "Problem " << ioerr << " in writeOrigRigidBodyAlignmentParameters, " << loop;
1157  ++result;
1158  }
1159  aliIO.writeAlignableAbsolutePositions(theAlignables, outFile.c_str(), loop, false, ioerr);
1160  if (ioerr) {
1161  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
1162  << "Problem " << ioerr << " in writeAlignableAbsolutePositions, " << loop;
1163  ++result;
1164  }
1165 
1166  return result;
1167 }
1168 
1169 //__________________________________________________________________________________________________
1171  for (const auto &iAli : alis) {
1172  AlignmentParameters *params = iAli->alignmentParameters();
1173  if (!params) {
1174  throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::buildUserVariables"
1175  << "No parameters for alignable";
1176  }
1177  MillePedeVariables *userVars = dynamic_cast<MillePedeVariables *>(params->userVariables());
1178  if (userVars) { // Just re-use existing, keeping label and numHits:
1179  for (unsigned int iPar = 0; iPar < userVars->size(); ++iPar) {
1180  // if (params->hierarchyLevel() > 0) {
1181  //std::cout << params->hierarchyLevel() << "\nBefore: " << userVars->parameter()[iPar];
1182  // }
1183  userVars->setAllDefault(iPar);
1184  //std::cout << "\nAfter: " << userVars->parameter()[iPar] << std::endl;
1185  }
1186  } else { // Nothing yet or erase wrong type:
1187  userVars = new MillePedeVariables(
1188  params->size(),
1189  thePedeLabels->alignableLabel(iAli),
1190  thePedeLabels->alignableTracker()->objectIdProvider().typeToName(iAli->alignableObjectId()));
1191  params->setUserVariables(userVars);
1192  }
1193  }
1194 }
1195 
1196 //__________________________________________________________________________________________________
1198  if (mode == "full") {
1200  } else if (mode == "mille") {
1201  return myMilleBit; // + myPedeSteerBit; // sic! Including production of steerig file. NO!
1202  } else if (mode == "pede") {
1204  } else if (mode == "pedeSteer") {
1205  return myPedeSteerBit;
1206  } else if (mode == "pedeRun") {
1207  return myPedeSteerBit + myPedeRunBit + myPedeReadBit; // sic! Including steering and reading of result.
1208  } else if (mode == "pedeRead") {
1209  return myPedeReadBit;
1210  }
1211 
1212  throw cms::Exception("BadConfig") << "Unknown mode '" << mode
1213  << "', use 'full', 'mille', 'pede', 'pedeRun', 'pedeSteer' or 'pedeRead'.";
1214 
1215  return 0;
1216 }
1217 
1218 //__________________________________________________________________________________________________
1220  const std::string &outFile,
1221  const std::vector<std::string> &inFiles) const {
1222  bool allOk = true;
1223  int ierr = 0;
1224  MillePedeVariablesIORoot millePedeIO;
1225  for (std::vector<std::string>::const_iterator iFile = inFiles.begin(); iFile != inFiles.end(); ++iFile) {
1226  const std::string inFile(theDir + *iFile);
1227  const std::vector<AlignmentUserVariables *> mpVars =
1228  millePedeIO.readMillePedeVariables(theAlignables, inFile.c_str(), fromIov, ierr);
1229  if (ierr || !this->addHits(theAlignables, mpVars)) {
1230  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addHitStatistics"
1231  << "Error " << ierr << " reading from " << inFile << ", tree " << fromIov
1232  << ", or problems in addHits";
1233  allOk = false;
1234  }
1235  for (std::vector<AlignmentUserVariables *>::const_iterator i = mpVars.begin(); i != mpVars.end(); ++i) {
1236  delete *i; // clean created objects
1237  }
1238  }
1239 
1240  return allOk;
1241 }
1242 
1243 //__________________________________________________________________________________________________
1245  const std::vector<AlignmentUserVariables *> &mpVars) const {
1246  bool allOk = (mpVars.size() == alis.size());
1247  std::vector<AlignmentUserVariables *>::const_iterator iUser = mpVars.begin();
1248  for (auto iAli = alis.cbegin(); iAli != alis.cend() && iUser != mpVars.end(); ++iAli, ++iUser) {
1249  MillePedeVariables *mpVarNew = dynamic_cast<MillePedeVariables *>(*iUser);
1250  AlignmentParameters *ps = (*iAli)->alignmentParameters();
1251  MillePedeVariables *mpVarOld = (ps ? dynamic_cast<MillePedeVariables *>(ps->userVariables()) : nullptr);
1252  if (!mpVarNew || !mpVarOld || mpVarOld->size() != mpVarNew->size()) {
1253  allOk = false;
1254  continue; // FIXME error etc.?
1255  }
1256 
1257  mpVarOld->increaseHitsX(mpVarNew->hitsX());
1258  mpVarOld->increaseHitsY(mpVarNew->hitsY());
1259  }
1260 
1261  return allOk;
1262 }
1263 
1264 //__________________________________________________________________________________________________
1265 template <typename GlobalDerivativeMatrix>
1266 void MillePedeAlignmentAlgorithm::makeGlobDerivMatrix(const std::vector<float> &globalDerivativesx,
1267  const std::vector<float> &globalDerivativesy,
1268  Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) {
1269  static_assert(GlobalDerivativeMatrix::RowsAtCompileTime == 2, "global derivative matrix must have two rows");
1270 
1271  for (size_t i = 0; i < globalDerivativesx.size(); ++i) {
1272  aGlobalDerivativesM(0, i) = globalDerivativesx[i];
1273  aGlobalDerivativesM(1, i) = globalDerivativesy[i];
1274  }
1275 }
1276 
1277 //__________________________________________________________________________________________________
1278 template <typename CovarianceMatrix,
1279  typename LocalDerivativeMatrix,
1280  typename ResidualMatrix,
1281  typename GlobalDerivativeMatrix>
1282 void MillePedeAlignmentAlgorithm::diagonalize(Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1283  Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM,
1284  Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1285  Eigen::MatrixBase<GlobalDerivativeMatrix> &aGlobalDerivativesM) const {
1287  "'aLocalDerivativesM' and 'aHitResidualsM' must have the "
1288  "same underlying scalar type");
1290  "'aLocalDerivativesM' and 'aGlobalDerivativesM' must have the "
1291  "same underlying scalar type");
1292 
1293  Eigen::SelfAdjointEigenSolver<typename CovarianceMatrix::PlainObject> myDiag{aHitCovarianceM};
1294  // eigenvectors of real symmetric matrices are orthogonal, i.e. invert == transpose
1295  auto aTranfoToDiagonalSystemInv =
1296  myDiag.eigenvectors().transpose().template cast<typename LocalDerivativeMatrix::Scalar>();
1297 
1298  aHitCovarianceM = myDiag.eigenvalues().asDiagonal();
1299  aLocalDerivativesM = aTranfoToDiagonalSystemInv * aLocalDerivativesM;
1300  aHitResidualsM = aTranfoToDiagonalSystemInv * aHitResidualsM;
1301  if (aGlobalDerivativesM.size() > 0) {
1302  // diagonalize only if measurement depends on alignables or calibrations
1303  aGlobalDerivativesM = aTranfoToDiagonalSystemInv * aGlobalDerivativesM;
1304  }
1305 }
1306 
1307 //__________________________________________________________________________________________________
1308 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1311  unsigned int iVirtualMeas,
1312  Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1313  Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1314  Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1315  // This Method is valid for 1D measurements only
1316 
1317  const unsigned int xIndex = iVirtualMeas + refTrajPtr->numberOfHitMeas();
1318 
1319  aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1320  aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex];
1321 
1322  const auto &locDerivMatrix = refTrajPtr->derivatives();
1323  for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1324  aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1325  }
1326 }
1327 
1328 //__________________________________________________________________________________________________
1329 template <typename CovarianceMatrix, typename ResidualMatrix, typename LocalDerivativeMatrix>
1331  unsigned int iTrajHit,
1332  Eigen::MatrixBase<CovarianceMatrix> &aHitCovarianceM,
1333  Eigen::MatrixBase<ResidualMatrix> &aHitResidualsM,
1334  Eigen::MatrixBase<LocalDerivativeMatrix> &aLocalDerivativesM) {
1335  // This Method is valid for 2D measurements only
1336 
1337  const unsigned int xIndex = iTrajHit * 2;
1338  const unsigned int yIndex = iTrajHit * 2 + 1;
1339 
1340  aHitCovarianceM(0, 0) = refTrajPtr->measurementErrors()[xIndex][xIndex];
1341  aHitCovarianceM(0, 1) = refTrajPtr->measurementErrors()[xIndex][yIndex];
1342  aHitCovarianceM(1, 0) = refTrajPtr->measurementErrors()[yIndex][xIndex];
1343  aHitCovarianceM(1, 1) = refTrajPtr->measurementErrors()[yIndex][yIndex];
1344 
1345  aHitResidualsM(0, 0) = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1346  aHitResidualsM(1, 0) = refTrajPtr->measurements()[yIndex] - refTrajPtr->trajectoryPositions()[yIndex];
1347 
1348  const auto &locDerivMatrix = refTrajPtr->derivatives();
1349  for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
1350  aLocalDerivativesM(0, i) = locDerivMatrix[xIndex][i];
1351  aLocalDerivativesM(1, i) = locDerivMatrix[yIndex][i];
1352  }
1353 }
1354 
1355 //__________________________________________________________________________________________________
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 //__________________________________________________________________________________________________
1372  unsigned int iTrajHit,
1373  const std::vector<int> &globalLabels,
1374  const std::vector<float> &globalDerivativesX) {
1375  const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1376  const unsigned int xIndex = iTrajHit * 2; // the even ones are local x
1377 
1378  // local derivatives
1379  const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
1380  const int nLocal = locDerivMatrix.num_col();
1381  std::vector<float> localDerivatives(nLocal);
1382  for (unsigned int i = 0; i < localDerivatives.size(); ++i) {
1383  localDerivatives[i] = locDerivMatrix[xIndex][i];
1384  }
1385 
1386  // residuum and error
1387  float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
1388  float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
1389 
1390  // number of global derivatives
1391  const int nGlobal = globalDerivativesX.size();
1392 
1393  // &(localDerivatives[0]) etc. are valid - as long as vector is not empty
1394  // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1395  theMille->mille(
1396  nLocal, &(localDerivatives[0]), nGlobal, &(globalDerivativesX[0]), &(globalLabels[0]), residX, hitErrX);
1397 
1398  if (theMonitor) {
1399  theMonitor->fillDerivatives(
1400  aRecHit, &(localDerivatives[0]), nLocal, &(globalDerivativesX[0]), nGlobal, &(globalLabels[0]));
1401  theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, residX, hitErrX, false);
1402  }
1403 
1404  return 1;
1405 }
1406 
1407 //__________________________________________________________________________________________________
1409  unsigned int iTrajHit,
1410  const std::vector<int> &globalLabels,
1411  const std::vector<float> &globalDerivativesx,
1412  const std::vector<float> &globalDerivativesy) {
1413  const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
1414 
1415  if ((aRecHit)->dimension() != 2) {
1416  edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::callMille2D"
1417  << "You try to call method for 2D hits for a " << (aRecHit)->dimension()
1418  << "D Hit. Hit gets ignored!";
1419  return -1;
1420  }
1421 
1422  Eigen::Matrix<double, 2, 2> aHitCovarianceM;
1423  Eigen::Matrix<float, 2, 1> aHitResidualsM;
1424  Eigen::Matrix<float, 2, Eigen::Dynamic> aLocalDerivativesM{2, refTrajPtr->derivatives().num_col()};
1425  // below method fills above 3 matrices
1426  this->addRefTrackData2D(refTrajPtr, iTrajHit, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1427  Eigen::Matrix<float, 2, Eigen::Dynamic> aGlobalDerivativesM{2, globalDerivativesx.size()};
1428  this->makeGlobDerivMatrix(globalDerivativesx, globalDerivativesy, aGlobalDerivativesM);
1429 
1430  // calculates correlation between Hit measurements
1431  // FIXME: Should take correlation (and resulting transformation) from original hit,
1432  // not 2x2 matrix from ReferenceTrajectory: That can come from error propagation etc.!
1433  const double corr = aHitCovarianceM(0, 1) / sqrt(aHitCovarianceM(0, 0) * aHitCovarianceM(1, 1));
1434  if (theMonitor)
1435  theMonitor->fillCorrelations2D(corr, aRecHit);
1436  bool diag = false; // diagonalise only tracker TID, TEC
1437  switch (aRecHit->geographicalId().subdetId()) {
1438  case SiStripDetId::TID:
1439  case SiStripDetId::TEC:
1440  if (aRecHit->geographicalId().det() == DetId::Tracker && TMath::Abs(corr) > theMaximalCor2D) {
1441  this->diagonalize(aHitCovarianceM, aLocalDerivativesM, aHitResidualsM, aGlobalDerivativesM);
1442  diag = true;
1443  }
1444  break;
1445  default:;
1446  }
1447 
1448  float newResidX = aHitResidualsM(0, 0);
1449  float newResidY = aHitResidualsM(1, 0);
1450  float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1451  float newHitErrY = TMath::Sqrt(aHitCovarianceM(1, 1));
1452 
1453  // change from column major (Eigen default) to row major to have row entries
1454  // in continuous memory
1455  std::vector<float> newLocalDerivs(aLocalDerivativesM.size());
1456  Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1457  newLocalDerivs.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1458  float *newLocalDerivsX = &(newLocalDerivs[0]);
1459  float *newLocalDerivsY = &(newLocalDerivs[aLocalDerivativesM.cols()]);
1460 
1461  // change from column major (Eigen default) to row major to have row entries
1462  // in continuous memory
1463  std::vector<float> newGlobDerivs(aGlobalDerivativesM.size());
1464  Eigen::Map<Eigen::Matrix<float, 2, Eigen::Dynamic, Eigen::RowMajor> >(
1465  newGlobDerivs.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1466  float *newGlobDerivsX = &(newGlobDerivs[0]);
1467  float *newGlobDerivsY = &(newGlobDerivs[aGlobalDerivativesM.cols()]);
1468 
1469  const int nLocal = aLocalDerivativesM.cols();
1470  const int nGlobal = aGlobalDerivativesM.cols();
1471 
1472  if (diag && (newHitErrX > newHitErrY)) { // also for 2D hits?
1473  // measurement with smaller error is x-measurement (for !is2D do not fill y-measurement):
1474  std::swap(newResidX, newResidY);
1475  std::swap(newHitErrX, newHitErrY);
1476  std::swap(newLocalDerivsX, newLocalDerivsY);
1477  std::swap(newGlobDerivsX, newGlobDerivsY);
1478  }
1479 
1480  // &(globalLabels[0]) is valid - as long as vector is not empty
1481  // cf. http://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
1482  theMille->mille(nLocal, newLocalDerivsX, nGlobal, newGlobDerivsX, &(globalLabels[0]), newResidX, newHitErrX);
1483 
1484  if (theMonitor) {
1485  theMonitor->fillDerivatives(aRecHit, newLocalDerivsX, nLocal, newGlobDerivsX, nGlobal, &(globalLabels[0]));
1486  theMonitor->fillResiduals(
1487  aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidX, newHitErrX, false);
1488  }
1489  const bool isReal2DHit = this->is2D(aRecHit); // strip is 1D (except matched hits)
1490  if (isReal2DHit) {
1491  theMille->mille(nLocal, newLocalDerivsY, nGlobal, newGlobDerivsY, &(globalLabels[0]), newResidY, newHitErrY);
1492  if (theMonitor) {
1493  theMonitor->fillDerivatives(aRecHit, newLocalDerivsY, nLocal, newGlobDerivsY, nGlobal, &(globalLabels[0]));
1494  theMonitor->fillResiduals(
1495  aRecHit, refTrajPtr->trajectoryStates()[iTrajHit], iTrajHit, newResidY, newHitErrY, true); // true: y
1496  }
1497  }
1498 
1499  return (isReal2DHit ? 2 : 1);
1500 }
1501 
1502 //__________________________________________________________________________________________________
1504  unsigned int iVirtualMeas) {
1505  Eigen::Matrix<double, 1, 1> aHitCovarianceM;
1506  Eigen::Matrix<float, 1, 1> aHitResidualsM;
1507  Eigen::Matrix<float, 1, Eigen::Dynamic> aLocalDerivativesM{1, refTrajPtr->derivatives().num_col()};
1508  // below method fills above 3 'matrices'
1509  this->addRefTrackVirtualMeas1D(refTrajPtr, iVirtualMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
1510 
1511  // no global parameters (use dummy 0)
1512  auto aGlobalDerivativesM = Eigen::Matrix<float, 1, 1>::Zero();
1513 
1514  float newResidX = aHitResidualsM(0, 0);
1515  float newHitErrX = TMath::Sqrt(aHitCovarianceM(0, 0));
1516  std::vector<float> newLocalDerivsX(aLocalDerivativesM.size());
1517  Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1518  newLocalDerivsX.data(), aLocalDerivativesM.rows(), aLocalDerivativesM.cols()) = aLocalDerivativesM;
1519 
1520  std::vector<float> newGlobDerivsX(aGlobalDerivativesM.size());
1521  Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(
1522  newGlobDerivsX.data(), aGlobalDerivativesM.rows(), aGlobalDerivativesM.cols()) = aGlobalDerivativesM;
1523 
1524  const int nLocal = aLocalDerivativesM.cols();
1525  const int nGlobal = 0;
1526 
1527  theMille->mille(nLocal, newLocalDerivsX.data(), nGlobal, newGlobDerivsX.data(), &nGlobal, newResidX, newHitErrX);
1528 }
1529 
1530 //____________________________________________________
1532  const TkFittedLasBeamCollection &lasBeams,
1533  const TsosVectorCollection &lasBeamTsoses) {
1534  TsosVectorCollection::const_iterator iTsoses = lasBeamTsoses.begin();
1535  for (TkFittedLasBeamCollection::const_iterator iBeam = lasBeams.begin(), iEnd = lasBeams.end(); iBeam != iEnd;
1536  ++iBeam, ++iTsoses) { // beam/tsoses parallel!
1537 
1538  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addLaserData"
1539  << "Beam " << iBeam->getBeamId() << " with " << iBeam->parameters().size()
1540  << " parameters and " << iBeam->getData().size() << " hits.\n There are "
1541  << iTsoses->size() << " TSOSes.";
1542 
1543  this->addLasBeam(eventInfo, *iBeam, *iTsoses);
1544  }
1545 }
1546 
1547 //____________________________________________________
1549  const TkFittedLasBeam &lasBeam,
1550  const std::vector<TrajectoryStateOnSurface> &tsoses) {
1551  AlignmentParameters *dummyPtr = nullptr; // for globalDerivativesHierarchy()
1552  std::vector<float> lasLocalDerivsX; // buffer for local derivatives
1553  const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId()); // for global par
1554 
1555  for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
1556  if (!tsoses[iHit].isValid())
1557  continue;
1558  // clear buffer
1559  theFloatBufferX.clear();
1560  theFloatBufferY.clear();
1561  theIntBuffer.clear();
1562  lasLocalDerivsX.clear();
1563  // get alignables and global parameters
1564  const SiStripLaserRecHit2D &hit = lasBeam.getData()[iHit];
1565  AlignableDetOrUnitPtr lasAli(theAlignableNavigator->alignableFromDetId(hit.getDetId()));
1567  eventInfo, tsoses[iHit], lasAli, lasAli, theFloatBufferX, theFloatBufferY, theIntBuffer, dummyPtr);
1568  // fill derivatives vector from derivatives matrix
1569  for (unsigned int nFitParams = 0; nFitParams < static_cast<unsigned int>(lasBeam.parameters().size());
1570  ++nFitParams) {
1571  const float derivative = lasBeam.derivatives()[iHit][nFitParams];
1572  if (nFitParams < lasBeam.firstFixedParameter()) { // first local beam parameters
1573  lasLocalDerivsX.push_back(derivative);
1574  } else { // now global ones
1575  const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
1576  theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
1577  theFloatBufferX.push_back(derivative);
1578  }
1579  } // end loop over parameters
1580 
1581  const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
1582  // error from file or assume 0.003
1583  const float error = 0.003; // hit.localPositionError().xx(); sqrt???
1584 
1585  theMille->mille(lasLocalDerivsX.size(),
1586  &(lasLocalDerivsX[0]),
1587  theFloatBufferX.size(),
1588  &(theFloatBufferX[0]),
1589  &(theIntBuffer[0]),
1590  residual,
1591  error);
1592  } // end of loop over hits
1593 
1594  theMille->end();
1595 }
1596 
1598  // do some printing, if requested
1599  const bool doOutputOnStdout(pxbSurveyCfg.getParameter<bool>("doOutputOnStdout"));
1600  if (doOutputOnStdout) {
1601  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1602  << "# Output from addPxbSurvey follows below because "
1603  << "doOutputOnStdout is set to True";
1604  }
1605 
1606  // instantiate a dicer object
1607  SurveyPxbDicer dicer(pxbSurveyCfg.getParameter<std::vector<edm::ParameterSet> >("toySurveyParameters"),
1608  pxbSurveyCfg.getParameter<unsigned int>("toySurveySeed"));
1609  std::ofstream outfile(pxbSurveyCfg.getUntrackedParameter<std::string>("toySurveyFile").c_str());
1610 
1611  // read data from file
1612  std::vector<SurveyPxbImageLocalFit> measurements;
1613  std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
1615 
1616  // loop over photographs (=measurements) and perform the fit
1617  for (std::vector<SurveyPxbImageLocalFit>::size_type i = 0; i != measurements.size(); i++) {
1618  if (doOutputOnStdout) {
1619  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1620  << "Module " << i << ": ";
1621  }
1622 
1623  // get the Alignables and their surfaces
1624  AlignableDetOrUnitPtr mod1(theAlignableNavigator->alignableFromDetId(measurements[i].getIdFirst()));
1625  AlignableDetOrUnitPtr mod2(theAlignableNavigator->alignableFromDetId(measurements[i].getIdSecond()));
1626  const AlignableSurface &surf1 = mod1->surface();
1627  const AlignableSurface &surf2 = mod2->surface();
1628 
1629  // the position of the fiducial points in local frame of a PXB module
1630  const LocalPoint fidpoint0(-0.91, +3.30);
1631  const LocalPoint fidpoint1(+0.91, +3.30);
1632  const LocalPoint fidpoint2(+0.91, -3.30);
1633  const LocalPoint fidpoint3(-0.91, -3.30);
1634 
1635  // We choose the local frame of the first module as reference,
1636  // so take the fidpoints of the second module and calculate their
1637  // positions in the reference frame
1638  const GlobalPoint surf2point0(surf2.toGlobal(fidpoint0));
1639  const GlobalPoint surf2point1(surf2.toGlobal(fidpoint1));
1640  const LocalPoint fidpoint0inSurf1frame(surf1.toLocal(surf2point0));
1641  const LocalPoint fidpoint1inSurf1frame(surf1.toLocal(surf2point1));
1642 
1643  // Create the vector for the fit
1645  fidpointvec.push_back(fidpoint0inSurf1frame);
1646  fidpointvec.push_back(fidpoint1inSurf1frame);
1647  fidpointvec.push_back(fidpoint2);
1648  fidpointvec.push_back(fidpoint3);
1649 
1650  // if toy survey is requested, dice the values now
1651  if (pxbSurveyCfg.getParameter<bool>("doToySurvey")) {
1652  dicer.doDice(fidpointvec, measurements[i].getIdPair(), outfile);
1653  }
1654 
1655  // do the fit
1656  measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
1657  SurveyPxbImageLocalFit::localpars_t a; // local pars from fit
1658  a = measurements[i].getLocalParameters();
1659  const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
1660 
1661  // do some reporting, if requested
1662  if (doOutputOnStdout) {
1663  edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addPxbSurvey"
1664  << "a: " << a[0] << ", " << a[1] << ", " << a[2] << ", " << a[3]
1665  << " S= " << sqrt(a[2] * a[2] + a[3] * a[3]) << " phi= " << atan(a[3] / a[2])
1666  << " chi2= " << chi2 << std::endl;
1667  }
1668  if (theMonitor) {
1669  theMonitor->fillPxbSurveyHistsChi2(chi2);
1670  theMonitor->fillPxbSurveyHistsLocalPars(a[0], a[1], sqrt(a[2] * a[2] + a[3] * a[3]), atan(a[3] / a[2]));
1671  }
1672 
1673  // pass the results from the local fit to mille
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 
1689  const auto runRangeSelection = theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection");
1690 
1691  if (runRangeSelection.empty())
1692  return false;
1693 
1694  const auto runRanges =
1696 
1697  return !(runRanges.empty());
1698 }
bool processesEvents() override
Returns whether MP should process events in the current configuration.
align::GlobalPoints toGlobal(const align::LocalPoints &) const
Return in global coord given a set of local points.
const TimeTypeSpecs timeTypeSpecs[]
Definition: Time.cc:16
std::shared_ptr< AlignPCLThresholdsHG > theThresholds
const ConstTrajTrackPairCollection & trajTrackPairs() const
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
void resetParameters(void)
reset parameters, correlations, user variables
const IOVSyncValue & first() const
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:58
bool supportsCalibrations() override
Returns whether MP supports calibrations.
Alignable * mother() const
Return pointer to container alignable (if any)
Definition: Alignable.h:91
std::unordered_map< std::string, std::vector< float > > param_map
std::unique_ptr< MillePedeMonitor > theMonitor
void increaseHitsX(unsigned int add=1)
increase hits for x-measurement
bool setParametersForRunRange(const RunRange &runrange) override
void globalDerivativesCalibration(const TransientTrackingRecHit::ConstRecHitPointer &recHit, const TrajectoryStateOnSurface &tsos, const edm::EventSetup &setup, const EventInfo &eventInfo, std::vector< float > &globalDerivativesX, std::vector< float > &globalDerivativesY, std::vector< int > &globalLabels) const
adding derivatives from integrated calibrations
bool addCalibrations(const std::vector< IntegratedCalibrationBase *> &iCals) override
Pass integrated calibrations to Millepede (they are not owned by Millepede!)
Derivative< X, A >::type derivative(const A &_)
Definition: Derivative.h:18
virtual void endRun(const EventInfo &, const EndRunInfo &, const edm::EventSetup &)
Run on run products, e.g. TkLAS.
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:132
const bool isValid(const Frame &aFrame, const FrameQuality &aQuality, const uint16_t aExpectedPos)
std::vector< coord_t > fidpoint_t
def create(alignables, pedeDump, additionalData, outputFile, config)
std::map< std::string, AlignPCLThreshold > threshold_map
void increaseHitsY(unsigned int add=1)
increase hits for y-measurement
Time_t beginValue
Definition: Time.h:41
void endLuminosityBlock(const edm::EventSetup &) override
called at end of luminosity block
const edm::ESGetToken< SiPixelQuality, SiPixelQualityFromDbRcd > siPixelQualityToken_
static constexpr auto TID
Definition: SiStripDetId.h:39
TrajectoryFactoryBase::ReferenceTrajectoryCollection RefTrajColl
std::vector< ParameterSet > VParameterSet
Definition: ParameterSet.h:35
std::unique_ptr< TrajectoryFactoryBase > theTrajectoryFactory
bool areEmptyParams(const align::Alignables &alignables) const
std::unique_ptr< gbl::MilleBinary > theBinary
bool exists(std::string const &parameterName) const
checks if a parameter exists
const Time_t MIN_VAL(0)
int addGlobalData(const edm::EventSetup &setup, const EventInfo &eventInfo, const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iHit, gbl::GblPoint &gblPoint)
void restoreCachedTransformations(void)
restore the previously cached position, rotation and other parameters
std::shared_ptr< PedeLabelerBase > thePedeLabels
bool globalDerivativesHierarchy(const EventInfo &eventInfo, const TrajectoryStateOnSurface &tsos, Alignable *ali, const AlignableDetOrUnitPtr &alidet, std::vector< float > &globalDerivativesX, std::vector< float > &globalDerivativesY, std::vector< int > &globalLabels, AlignmentParameters *&lowestParams) const
recursively adding derivatives and labels, false if problems
const align::Alignables & alignables(void) const
get all alignables
void applyParameters(void)
Obsolete: Use AlignableNavigator::alignableDetFromDetId and alignableFromAlignableDet.
std::vector< std::string > getExistingFormattedFiles(const std::vector< std::string > &plainFiles, const std::string &theDir)
bool isMode(unsigned int testMode) const
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > geomToken_
int callMille2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iTrajHit, const std::vector< int > &globalLabels, const std::vector< float > &globalDerivativesx, const std::vector< float > &globalDerivativesy)
Log< level::Error, false > LogError
const std::vector< SiStripLaserRecHit2D > & getData(void) const
access the collection of hits
Definition: TkLasBeam.h:26
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
void writeMillePedeVariables(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
void writeAlignableOriginalPositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr) override
write Alignable original (before misalignment) absolute positions
uint16_t size_type
void addLasBeam(const EventInfo &eventInfo, const TkFittedLasBeam &lasBeam, const std::vector< TrajectoryStateOnSurface > &tsoses)
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:112
define event information passed to algorithms
T getUntrackedParameter(std::string const &, T const &) const
int callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iTrajHit, const std::vector< int > &globalLabels, const std::vector< float > &globalDerivativesX, const std::vector< float > &globalDerivativesY)
calls callMille1D or callMille2D
TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer
U second(std::pair< T, U > const &p)
void beginLuminosityBlock(const edm::EventSetup &) override
called at begin of luminosity block (resets Mille binary in mille mode)
std::vector< align::RunNumber > cachedRuns_
std::string doDice(const fidpoint_t &fidpointvec, const idPair_t &id, const bool rotate=false)
unsigned int hitsX() const
get number of hits for x-measurement
int callMille1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iTrajHit, const std::vector< int > &globalLabels, const std::vector< float > &globalDerivativesX)
calls Mille for 1D hits
void addLaserData(const EventInfo &eventInfo, const TkFittedLasBeamCollection &tkLasBeams, const TsosVectorCollection &tkLasBeamTsoses)
int addMeasurementData(const edm::EventSetup &setup, const EventInfo &eventInfo, const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iHit, AlignmentParameters *&params)
virtual bool storeThresholds(const int &nRecords, const AlignPCLThresholdsHG::threshold_map &thresholdMap, const AlignPCLThresholdsHG::param_map &floatMap)
void beginRun(const edm::Run &run, const edm::EventSetup &setup, bool changed) override
called at begin of run
dictionary corr
CLHEP::HepMatrix AlgebraicMatrix
RunRanges makeNonOverlappingRunRanges(const edm::VParameterSet &runRanges, const RunNumber &defaultRun)
Definition: Utilities.cc:202
bool is2D(const TransientTrackingRecHit::ConstRecHitPointer &recHit) const
true if hit belongs to 2D detector (currently tracker specific)
void makeGlobDerivMatrix(const std::vector< float > &globalDerivativesx, const std::vector< float > &globalDerivativesy, Eigen::MatrixBase< GlobalDerivativeMatrix > &aGlobalDerivativesM)
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
void writeOrigRigidBodyAlignmentParameters(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr) override
write RigidBodyAlignmentParameters as applied on top of original positions
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > topoToken_
T sqrt(T t)
Definition: SSEVec.h:23
RunRanges makeUniqueRunRanges(const edm::VParameterSet &runRanges, const RunNumber &defaultRun)
Definition: Utilities.cc:241
const AlgebraicMatrix & derivatives() const
matrix of local derivatives: columns are parameters, rows are hits
std::shared_ptr< PixelTopologyMap > pixelTopologyMap
const edm::ESGetToken< AlignPCLThresholdsHG, AlignPCLThresholdsHGRcd > aliThrToken_
AlignmentParameterStore * theAlignmentParameterStore
directory for all kind of files
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
bool addHitStatistics(int fromLoop, const std::string &outFile, const std::vector< std::string > &inFiles) const
void diagonalize(Eigen::MatrixBase< CovarianceMatrix > &aHitCovarianceM, Eigen::MatrixBase< LocalDerivativeMatrix > &aLocalDerivativesM, Eigen::MatrixBase< ResidualMatrix > &aHitResidualsM, Eigen::MatrixBase< GlobalDerivativeMatrix > &aGlobalDerivativesM) const
std::vector< AlignmentUserVariables * > readMillePedeVariables(const align::Alignables &alivec, const char *filename, int iter, int &ierr)
Transition
Definition: Transition.h:12
bool addHits(const align::Alignables &alis, const std::vector< AlignmentUserVariables *> &mpVars) const
void addRefTrackVirtualMeas1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iVirtualMeas, Eigen::MatrixBase< CovarianceMatrix > &aHitCovarianceM, Eigen::MatrixBase< ResidualMatrix > &aHitResidualsM, Eigen::MatrixBase< LocalDerivativeMatrix > &aLocalDerivativesM)
adds data for a specific virtual measurement from reference trajectory
key
prepare the HTCondor submission files and eventually submit them
void addRefTrackData2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iTrajHit, Eigen::MatrixBase< CovarianceMatrix > &aHitCovarianceM, Eigen::MatrixBase< ResidualMatrix > &aHitResidualsM, Eigen::MatrixBase< LocalDerivativeMatrix > &aLocalDerivativesM)
adds data from reference trajectory from a specific Hit
unsigned int decodeMode(const std::string &mode) const
unsigned int addHitCount(const std::vector< AlignmentParameters *> &parVec, const std::vector< bool > &validHitVecY) const
std::unique_ptr< PedeSteerer > thePedeSteer
bool storeAlignments() override
Returns whether MP produced results to be stored.
void terminate() override
Called at end of job.
unsigned int getBeamId(void) const
return the full beam identifier
Definition: TkLasBeam.h:23
void initialize(const edm::EventSetup &setup, AlignableTracker *tracker, AlignableMuon *muon, AlignableExtras *extras, AlignmentParameterStore *store) override
Called at beginning of job.
std::vector< ConstRecHitPointer > ConstRecHitContainer
unsigned int firstFixedParameter() const
Log< level::Info, false > LogInfo
void addVirtualMeas(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iVirtualMeas)
adds data for virtual measurements from reference trajectory
RunNumber_t run() const
Definition: EventID.h:38
bool readFromPede(const edm::ParameterSet &mprespset, bool setUserVars, const RunRange &runrange)
read pede input defined by &#39;psetName&#39;, flag to create/not create MillePedeVariables ...
std::pair< unsigned int, unsigned int > addReferenceTrajectory(const edm::EventSetup &setup, const EventInfo &eventInfo, const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr)
fill mille for a trajectory, returning number of x/y hits ([0,0] if &#39;bad&#39; trajectory) ...
void buildUserVariables(const align::Alignables &alignables) const
add MillePedeVariables for each AlignmentParameters (exception if no parameters...)
void addUntrackedParameter(std::string const &name, T const &value)
Definition: ParameterSet.h:193
std::vector< IntegratedCalibrationBase * > theCalibrations
Class to hold one picture of the BPix survey.
std::vector< value_t > localpars_t
void writeAlignableAbsolutePositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr) override
write Alignable current absolute positions
std::vector< Alignable * > Alignables
Definition: Utilities.h:31
const EventID & eventID() const
Definition: IOVSyncValue.h:40
void run(const edm::EventSetup &setup, const EventInfo &eventInfo) override
Run the algorithm on trajectories and tracks.
HLT enums.
~MillePedeAlignmentAlgorithm() override
Destructor.
MillePedeAlignmentAlgorithm(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC)
Constructor.
align::RotationType toLocal(const align::RotationType &) const
Return in local frame a rotation given in global frame.
std::vector< TkFittedLasBeam > TkFittedLasBeamCollection
double a
Definition: hdecay.h:121
Eigen::Matrix< float, EcalPulseShape::TEMPLATESAMPLES, EcalPulseShape::TEMPLATESAMPLES > CovarianceMatrix
void cacheTransformations(void)
cache the current position, rotation and other parameters
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
unsigned int size() const
number of parameters
const std::string & fullPath() const
Definition: FileInPath.cc:144
static const count_t nMsrmts
eventInfo
add run, event number and lumi section
void addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg)
add measurement data from PXB survey
#define get
std::unique_ptr< AlignableNavigator > theAlignableNavigator
Log< level::Warning, false > LogWarning
Time_t endValue
Definition: Time.h:42
const std::vector< Scalar > & parameters() const
parallel to derivatives()
std::vector< std::vector< TrajectoryStateOnSurface > > TsosVectorCollection
define run information passed to algorithms (in endRun)
unsigned int hitsY() const
get number of hits for y-measurement
Constructor of the full muon geometry.
Definition: AlignableMuon.h:38
uint32_t dimension(pat::CandKinResolution::Parametrization parametrization)
Returns the number of free parameters in a parametrization (3 or 4)
ValidityInterval validityInterval() const
const Time_t MAX_VAL(std::numeric_limits< Time_t >::max())
std::vector< ReferenceTrajectoryPtr > ReferenceTrajectoryCollection
void throwException(const std::string &message, const std::string &methodName)
Definition: Exception.cc:12
bool setAllDefault(unsigned int nParam)
set default values for all data concerning nParam (false if nParam out of range)
static constexpr auto TEC
Definition: SiStripDetId.h:41
std::shared_ptr< SiPixelQuality > pixelQuality
Definition: Run.h:45