00001
00011 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeAlignmentAlgorithm.h"
00012
00013
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015
00016 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00017
00018
00019
00020 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeMonitor.h"
00021 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeVariables.h"
00022 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeVariablesIORoot.h"
00023 #include "Mille.h"
00024 #include "PedeSteerer.h"
00025 #include "PedeReader.h"
00026 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
00027 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerPluginFactory.h"
00028
00029 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryBase.h"
00030 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryPlugin.h"
00031
00032 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
00033 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentIORoot.h"
00034
00035 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
00036 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
00037 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
00038
00039
00040 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
00041 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
00042 #include "Alignment/CommonAlignment/interface/AlignableExtras.h"
00043
00044 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00045 #include "DataFormats/TrackReco/interface/Track.h"
00046 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
00047 #include "DataFormats/Alignment/interface/TkFittedLasBeam.h"
00048 #include "Alignment/LaserAlignment/interface/TsosVectorCollection.h"
00049
00050 #include <Geometry/CommonDetUnit/interface/GeomDetUnit.h>
00051 #include <Geometry/CommonDetUnit/interface/GeomDetType.h>
00052
00053 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00054
00055 #include <fstream>
00056 #include <sstream>
00057 #include <algorithm>
00058
00059 #include <TMath.h>
00060 #include <TMatrixDSymEigen.h>
00061 typedef TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer;
00062 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
00063 typedef TrajectoryFactoryBase::ReferenceTrajectoryCollection RefTrajColl;
00064
00065
00066 #include <iostream>
00067 #include "Alignment/SurveyAnalysis/interface/SurveyPxbImage.h"
00068 #include "Alignment/SurveyAnalysis/interface/SurveyPxbImageLocalFit.h"
00069 #include "Alignment/SurveyAnalysis/interface/SurveyPxbImageReader.h"
00070 #include "Alignment/SurveyAnalysis/interface/SurveyPxbDicer.h"
00071 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00072 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00073
00074 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
00075
00076
00077
00078 MillePedeAlignmentAlgorithm::MillePedeAlignmentAlgorithm(const edm::ParameterSet &cfg) :
00079 AlignmentAlgorithmBase(cfg),
00080 theConfig(cfg), theMode(this->decodeMode(theConfig.getUntrackedParameter<std::string>("mode"))),
00081 theDir(theConfig.getUntrackedParameter<std::string>("fileDir")),
00082 theAlignmentParameterStore(0), theAlignables(), theAlignableNavigator(0),
00083 theMonitor(0), theMille(0), thePedeLabels(0), thePedeSteer(0),
00084 theTrajectoryFactory(0),
00085 theMinNumHits(cfg.getParameter<unsigned int>("minNumHits")),
00086 theMaximalCor2D(cfg.getParameter<double>("max2Dcorrelation"))
00087 {
00088 if (!theDir.empty() && theDir.find_last_of('/') != theDir.size()-1) theDir += '/';
00089 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm" << "Start in mode '"
00090 << theConfig.getUntrackedParameter<std::string>("mode")
00091 << "' with output directory '" << theDir << "'.";
00092 }
00093
00094
00095
00096 MillePedeAlignmentAlgorithm::~MillePedeAlignmentAlgorithm()
00097 {
00098 delete theAlignableNavigator;
00099 theAlignableNavigator = 0;
00100 delete theMille;
00101 theMille = 0;
00102 delete theMonitor;
00103 theMonitor = 0;
00104 delete thePedeSteer;
00105 thePedeSteer = 0;
00106 delete thePedeLabels;
00107 thePedeLabels = 0;
00108 delete theTrajectoryFactory;
00109 theTrajectoryFactory = 0;
00110 }
00111
00112
00113
00114 void MillePedeAlignmentAlgorithm::initialize(const edm::EventSetup &setup,
00115 AlignableTracker *tracker, AlignableMuon *muon, AlignableExtras *extras,
00116 AlignmentParameterStore *store)
00117 {
00118 if (muon) {
00119 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
00120 << "Running with AlignabeMuon not yet tested.";
00121 }
00122
00123 theAlignableNavigator = new AlignableNavigator(extras, tracker, muon);
00124 theAlignmentParameterStore = store;
00125 theAlignables = theAlignmentParameterStore->alignables();
00126
00127 edm::ParameterSet pedeLabelerCfg(theConfig.getParameter<edm::ParameterSet>("pedeLabeler"));
00128 edm::VParameterSet RunRangeSelectionVPSet(theConfig.getUntrackedParameter<edm::VParameterSet>("RunRangeSelection"));
00129 pedeLabelerCfg.addUntrackedParameter<edm::VParameterSet>("RunRangeSelection",
00130 RunRangeSelectionVPSet);
00131
00132 std::string labelerPlugin = "PedeLabeler";
00133 if (RunRangeSelectionVPSet.size()>0) {
00134 labelerPlugin = "RunRangeDependentPedeLabeler";
00135 if (pedeLabelerCfg.exists("plugin")) {
00136 std::string labelerPluginCfg = pedeLabelerCfg.getParameter<std::string>("plugin");
00137 if ((labelerPluginCfg!="PedeLabeler" && labelerPluginCfg!="RunRangeDependentPedeLabeler") ||
00138 pedeLabelerCfg.getUntrackedParameter<edm::VParameterSet>("parameterInstances").size()>0) {
00139 throw cms::Exception("BadConfig")
00140 << "MillePedeAlignmentAlgorithm::initialize"
00141 << "both RunRangeSelection and generic labeler specified in config file. "
00142 << "Please get rid of either one of them.\n";
00143 }
00144 }
00145 } else {
00146 if (pedeLabelerCfg.exists("plugin")) {
00147 labelerPlugin = pedeLabelerCfg.getParameter<std::string>("plugin");
00148 }
00149 }
00150
00151 if (!pedeLabelerCfg.exists("plugin")) {
00152 pedeLabelerCfg.addUntrackedParameter<std::string>("plugin", labelerPlugin);
00153 }
00154
00155 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
00156 << "Using plugin '" << labelerPlugin << "' to generate labels.";
00157
00158 thePedeLabels = PedeLabelerPluginFactory::get()->create(labelerPlugin,
00159 PedeLabelerBase::TopLevelAlignables(tracker, muon, extras),
00160 pedeLabelerCfg);
00161
00162
00163 edm::ParameterSet pedeSteerCfg(theConfig.getParameter<edm::ParameterSet>("pedeSteerer"));
00164 thePedeSteer = new PedeSteerer(tracker, muon, extras,
00165 theAlignmentParameterStore, thePedeLabels,
00166 pedeSteerCfg, theDir, !this->isMode(myPedeSteerBit));
00167
00168
00169
00170 const std::vector<edm::ParameterSet> mprespset
00171 (theConfig.getParameter<std::vector<edm::ParameterSet> >("pedeReaderInputs"));
00172 if (!mprespset.empty()) {
00173 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::initialize"
00174 << "Apply " << mprespset.end() - mprespset.begin()
00175 << " previous MillePede constants from 'pedeReaderInputs'.";
00176 }
00177
00178 RunRange runrange(cond::timeTypeSpecs[cond::runnumber].beginValue,
00179 cond::timeTypeSpecs[cond::runnumber].endValue);
00180 for (std::vector<edm::ParameterSet>::const_iterator iSet = mprespset.begin(), iE = mprespset.end();
00181 iSet != iE; ++iSet) {
00182 if (!this->readFromPede((*iSet), false, runrange)) {
00183 throw cms::Exception("BadConfig")
00184 << "MillePedeAlignmentAlgorithm::initialize: Problems reading input constants of "
00185 << "pedeReaderInputs entry " << iSet - mprespset.begin() << '.';
00186 }
00187 theAlignmentParameterStore->applyParameters();
00188
00189 theAlignmentParameterStore->resetParameters();
00190 }
00191
00192
00193 thePedeSteer->buildSubSteer(tracker, muon, extras);
00194
00195
00196 this->buildUserVariables(theAlignables);
00197
00198 if (this->isMode(myMilleBit)) {
00199 if (!theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles").empty() ||
00200 !theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles").empty()) {
00201 throw cms::Exception("BadConfig")
00202 << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' must be empty for "
00203 << "modes running mille.";
00204 }
00205 theMille = new Mille((theDir + theConfig.getParameter<std::string>("binaryFile")).c_str());
00206 const std::string moniFile(theConfig.getUntrackedParameter<std::string>("monitorFile"));
00207 if (moniFile.size()) theMonitor = new MillePedeMonitor((theDir + moniFile).c_str());
00208
00209
00210 const edm::ParameterSet fctCfg(theConfig.getParameter<edm::ParameterSet>("TrajectoryFactory"));
00211 const std::string fctName(fctCfg.getParameter<std::string>("TrajectoryFactoryName"));
00212 theTrajectoryFactory = TrajectoryFactoryPlugin::get()->create(fctName, fctCfg);
00213 }
00214
00215
00216 if (this->isMode(myPedeSteerBit)
00217 || !theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles").empty()) {
00218 this->doIO(0);
00219
00220
00221 const edm::ParameterSet pxbSurveyCfg(theConfig.getParameter<edm::ParameterSet>("surveyPixelBarrel"));
00222 theDoSurveyPixelBarrel = pxbSurveyCfg.getParameter<bool>("doSurvey");
00223 if (theDoSurveyPixelBarrel) addPxbSurvey(pxbSurveyCfg);
00224 }
00225 }
00226
00227 void MillePedeAlignmentAlgorithm::setParametersForRunRange(const RunRange &runrange)
00228 {
00229 if (this->isMode(myPedeReadBit)) {
00230 if (!this->readFromPede(theConfig.getParameter<edm::ParameterSet>("pedeReader"), true, runrange)) {
00231 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::terminate"
00232 << "Problems reading pede result, but applying!";
00233 }
00234 theAlignmentParameterStore->applyParameters();
00235 }
00236 }
00237
00238
00239
00240 void MillePedeAlignmentAlgorithm::terminate()
00241 {
00242 delete theMille;
00243 theMille = 0;
00244
00245 std::vector<std::string> files;
00246 if (this->isMode(myMilleBit) || !theConfig.getParameter<std::string>("binaryFile").empty()) {
00247 files.push_back(theDir + theConfig.getParameter<std::string>("binaryFile"));
00248 } else {
00249 const std::vector<std::string> plainFiles
00250 (theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
00251 for (std::vector<std::string>::const_iterator i = plainFiles.begin(), iEnd = plainFiles.end();
00252 i != iEnd; ++i) {
00253 files.push_back(theDir + *i);
00254 }
00255 }
00256 const std::string masterSteer(thePedeSteer->buildMasterSteer(files));
00257 if (this->isMode(myPedeRunBit)) {
00258 thePedeSteer->runPede(masterSteer);
00259 }
00260
00261 if (this->isMode(myPedeReadBit)) {
00262 RunRange runrange(cond::timeTypeSpecs[cond::runnumber].beginValue,
00263 cond::timeTypeSpecs[cond::runnumber].endValue);
00264 if (!this->readFromPede(theConfig.getParameter<edm::ParameterSet>("pedeReader"), true, runrange)) {
00265 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::terminate"
00266 << "Problems reading pede result, but applying!";
00267 }
00268 theAlignmentParameterStore->applyParameters();
00269 }
00270
00271 if (this->isMode(myMilleBit)) {
00272 this->doIO(1);
00273 } else if (this->isMode(myPedeReadBit)) {
00274 this->doIO(2);
00275 }
00276 }
00277
00278
00279
00280 void MillePedeAlignmentAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo)
00281 {
00282 if (!this->isMode(myMilleBit)) return;
00283 const ConstTrajTrackPairCollection &tracks = eventInfo.trajTrackPairs_;
00284
00285 if (theMonitor) {
00286 for (ConstTrajTrackPairCollection::const_iterator iTrajTrack = tracks.begin();
00287 iTrajTrack != tracks.end(); ++iTrajTrack) {
00288 theMonitor->fillTrack((*iTrajTrack).second);
00289 }
00290 }
00291
00292 const RefTrajColl trajectories(theTrajectoryFactory->trajectories(setup, tracks, eventInfo.beamSpot_));
00293
00294
00295 unsigned int refTrajCount = 0;
00296 for (RefTrajColl::const_iterator iRefTraj = trajectories.begin(), iRefTrajE = trajectories.end();
00297 iRefTraj != iRefTrajE; ++iRefTraj, ++refTrajCount) {
00298
00299 RefTrajColl::value_type refTrajPtr = *iRefTraj;
00300 if (theMonitor) theMonitor->fillRefTrajectory(refTrajPtr);
00301
00302 const std::pair<unsigned int, unsigned int> nHitXy = this->addReferenceTrajectory(eventInfo, refTrajPtr);
00303
00304 if (theMonitor && (nHitXy.first || nHitXy.second)) {
00305
00306
00307 const reco::Track *trackPtr =
00308 (trajectories.size() == tracks.size() ? tracks[refTrajCount].second : 0);
00309 theMonitor->fillUsedTrack(trackPtr, nHitXy.first, nHitXy.second);
00310 }
00311
00312 }
00313 }
00314
00315
00316
00317
00318 std::pair<unsigned int, unsigned int>
00319 MillePedeAlignmentAlgorithm::addReferenceTrajectory(const EventInfo &eventInfo,
00320 const RefTrajColl::value_type &refTrajPtr)
00321 {
00322 std::pair<unsigned int, unsigned int> hitResultXy(0,0);
00323 if (refTrajPtr->isValid()) {
00324
00325
00326 std::vector<AlignmentParameters*> parVec(refTrajPtr->recHits().size());
00327
00328 std::vector<bool> validHitVecY(refTrajPtr->recHits().size(), false);
00329
00330 for (unsigned int iHit = 0; iHit < refTrajPtr->recHits().size(); ++iHit) {
00331 const int flagXY = this->addMeasurementData(eventInfo, refTrajPtr, iHit, parVec[iHit]);
00332
00333 if (flagXY < 0) {
00334 hitResultXy.first = 0;
00335 break;
00336 } else {
00337 if (flagXY >= 1) ++hitResultXy.first;
00338 validHitVecY[iHit] = (flagXY >= 2);
00339 }
00340 }
00341
00342
00343 for (unsigned int iMsMeas = 0; iMsMeas < refTrajPtr->numberOfMsMeas(); ++iMsMeas) {
00344 this->addMsMeas(refTrajPtr, iMsMeas);
00345 }
00346
00347
00348 if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits) {
00349 theMille->kill();
00350 hitResultXy.first = hitResultXy.second = 0;
00351 } else {
00352 theMille->end();
00353
00354 for (unsigned int iHit = 0; iHit < validHitVecY.size(); ++iHit) {
00355 if (!parVec[iHit]) continue;
00356 MillePedeVariables *mpVar = static_cast<MillePedeVariables*>(parVec[iHit]->userVariables());
00357 mpVar->increaseHitsX();
00358 if (validHitVecY[iHit]) {
00359 mpVar->increaseHitsY();
00360 ++hitResultXy.second;
00361 }
00362 }
00363 }
00364 }
00365
00366 return hitResultXy;
00367 }
00368
00369
00370 void MillePedeAlignmentAlgorithm::endRun(const EventInfo &eventInfo, const EndRunInfo &runInfo,
00371 const edm::EventSetup &setup)
00372 {
00373 if(runInfo.tkLasBeams_ && runInfo.tkLasBeamTsoses_){
00374
00375 this->addLaserData(eventInfo, *(runInfo.tkLasBeams_), *(runInfo.tkLasBeamTsoses_));
00376 }
00377 }
00378
00379
00380 int MillePedeAlignmentAlgorithm::addMeasurementData(const EventInfo &eventInfo,
00381 const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00382 unsigned int iHit,
00383 AlignmentParameters *¶ms)
00384 {
00385 params = 0;
00386 theFloatBufferX.clear();
00387 theFloatBufferY.clear();
00388 theIntBuffer.clear();
00389
00390 const TrajectoryStateOnSurface &tsos = refTrajPtr->trajectoryStates()[iHit];
00391 const ConstRecHitPointer &recHitPtr = refTrajPtr->recHits()[iHit];
00392
00393 if (!recHitPtr->isValid()) return 0;
00394
00395
00396 AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
00397
00398 if (!this->globalDerivativesHierarchy(eventInfo,
00399 tsos, alidet, alidet, theFloatBufferX,
00400 theFloatBufferY, theIntBuffer, params)) {
00401 return -1;
00402 } else if (theFloatBufferX.empty()) {
00403 return 0;
00404 } else {
00405 return this->callMille(refTrajPtr, iHit, theIntBuffer, theFloatBufferX, theFloatBufferY);
00406 }
00407 }
00408
00409
00410 bool MillePedeAlignmentAlgorithm
00411 ::globalDerivativesHierarchy(const EventInfo &eventInfo,
00412 const TrajectoryStateOnSurface &tsos,
00413 Alignable *ali, const AlignableDetOrUnitPtr &alidet,
00414 std::vector<float> &globalDerivativesX,
00415 std::vector<float> &globalDerivativesY,
00416 std::vector<int> &globalLabels,
00417 AlignmentParameters *&lowestParams) const
00418 {
00419
00420 if (!ali) return true;
00421
00422 if (false && theMonitor && alidet != ali) theMonitor->fillFrameToFrame(alidet, ali);
00423
00424 AlignmentParameters *params = ali->alignmentParameters();
00425
00426 if (params) {
00427 if (!lowestParams) lowestParams = params;
00428
00429 bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
00430 const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
00431
00432 if (0 == alignableLabel) {
00433 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::globalDerivativesHierarchy"
00434 << "Label not found, skip Alignable.";
00435 return false;
00436 }
00437
00438 const std::vector<bool> &selPars = params->selector();
00439 const AlgebraicMatrix derivs(params->derivatives(tsos, alidet));
00440
00441
00442 for (unsigned int iSel = 0; iSel < selPars.size(); ++iSel) {
00443 if (selPars[iSel]) {
00444 globalDerivativesX.push_back(derivs[iSel][kLocalX]
00445 /thePedeSteer->cmsToPedeFactor(iSel));
00446 if (hasSplitParameters==true) {
00447 globalLabels.push_back(thePedeLabels->parameterLabel(ali, iSel, eventInfo, tsos));
00448 } else {
00449 globalLabels.push_back(thePedeLabels->parameterLabel(alignableLabel, iSel));
00450 }
00451 globalDerivativesY.push_back(derivs[iSel][kLocalY]
00452 /thePedeSteer->cmsToPedeFactor(iSel));
00453 }
00454 }
00455
00456 if (thePedeSteer->isNoHiera(ali)) return true;
00457 }
00458
00459 return this->globalDerivativesHierarchy(eventInfo,
00460 tsos, ali->mother(), alidet,
00461 globalDerivativesX, globalDerivativesY,
00462 globalLabels, lowestParams);
00463 }
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500 bool MillePedeAlignmentAlgorithm::is2D(const ConstRecHitPointer &recHit) const
00501 {
00502
00503
00504 if (recHit->dimension() < 2) {
00505 return false;
00506 } else if (recHit->detUnit()) {
00507 return recHit->detUnit()->type().isTrackerPixel();
00508 } else {
00509 if (dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit->hit())) {
00510
00511 return false;
00512 } else {
00513 return true;
00514 }
00515 }
00516 }
00517
00518
00519 bool MillePedeAlignmentAlgorithm::readFromPede(const edm::ParameterSet &mprespset, bool setUserVars,
00520 const RunRange &runrange)
00521 {
00522 bool allEmpty = this->areEmptyParams(theAlignables);
00523
00524 PedeReader reader(mprespset, *thePedeSteer, *thePedeLabels, runrange);
00525 std::vector<Alignable*> alis;
00526 bool okRead = reader.read(alis, setUserVars);
00527 bool numMatch = true;
00528
00529 std::stringstream out("Read ");
00530 out << alis.size() << " alignables";
00531 if (alis.size() != theAlignables.size()) {
00532 out << " while " << theAlignables.size() << " in store";
00533 numMatch = false;
00534 }
00535 if (!okRead) out << ", but problems in reading";
00536 if (!allEmpty) out << ", possibly overwriting previous settings";
00537 out << ".";
00538
00539 if (okRead && allEmpty) {
00540 if (numMatch) {
00541 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00542 } else if (alis.size()) {
00543 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00544 } else {
00545 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00546 return false;
00547 }
00548 return true;
00549 }
00550
00551 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00552 return false;
00553 }
00554
00555
00556 bool MillePedeAlignmentAlgorithm::areEmptyParams(const std::vector<Alignable*> &alignables) const
00557 {
00558
00559 for (std::vector<Alignable*>::const_iterator iAli = alignables.begin();
00560 iAli != alignables.end(); ++iAli) {
00561 const AlignmentParameters *params = (*iAli)->alignmentParameters();
00562 if (params) {
00563 const AlgebraicVector &parVec(params->parameters());
00564 const AlgebraicMatrix &parCov(params->covariance());
00565 for (int i = 0; i < parVec.num_row(); ++i) {
00566 if (parVec[i] != 0.) return false;
00567 for (int j = i; j < parCov.num_col(); ++j) {
00568 if (parCov[i][j] != 0.) return false;
00569 }
00570 }
00571 }
00572 }
00573
00574 return true;
00575 }
00576
00577
00578 unsigned int MillePedeAlignmentAlgorithm::doIO(int loop) const
00579 {
00580 unsigned int result = 0;
00581
00582 const std::string outFilePlain(theConfig.getParameter<std::string>("treeFile"));
00583 if (outFilePlain.empty()) {
00584 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
00585 << "treeFile parameter empty => skip writing for 'loop' " << loop;
00586 return result;
00587 }
00588
00589 const std::string outFile(theDir + outFilePlain);
00590
00591 AlignmentIORoot aliIO;
00592 int ioerr = 0;
00593 if (loop == 0) {
00594 aliIO.writeAlignableOriginalPositions(theAlignables, outFile.c_str(), loop, false, ioerr);
00595 if (ioerr) {
00596 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
00597 << "Problem " << ioerr << " in writeAlignableOriginalPositions";
00598 ++result;
00599 }
00600 } else {
00601 if (loop > 1) {
00602 const std::vector<std::string> inFiles
00603 (theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles"));
00604 const std::vector<std::string> binFiles
00605 (theConfig.getParameter<std::vector<std::string> >("mergeBinaryFiles"));
00606 if (inFiles.size() != binFiles.size()) {
00607 edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
00608 << "'vstring mergeTreeFiles' and 'vstring mergeBinaryFiles' "
00609 << "differ in size";
00610 }
00611 this->addHitStatistics(loop - 1, outFile, inFiles);
00612 }
00613 MillePedeVariablesIORoot millePedeIO;
00614 millePedeIO.writeMillePedeVariables(theAlignables, outFile.c_str(), loop, false, ioerr);
00615 if (ioerr) {
00616 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
00617 << "Problem " << ioerr << " writing MillePedeVariables";
00618 ++result;
00619 }
00620
00621
00622
00623
00624
00625
00626
00627 }
00628
00629
00630 aliIO.writeOrigRigidBodyAlignmentParameters(theAlignables, outFile.c_str(), loop, false, ioerr);
00631 if (ioerr) {
00632 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO" << "Problem " << ioerr
00633 << " in writeOrigRigidBodyAlignmentParameters, " << loop;
00634 ++result;
00635 }
00636 aliIO.writeAlignableAbsolutePositions(theAlignables, outFile.c_str(), loop, false, ioerr);
00637 if (ioerr) {
00638 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO" << "Problem " << ioerr
00639 << " in writeAlignableAbsolutePositions, " << loop;
00640 ++result;
00641 }
00642 aliIO.writeAlignableRelativePositions(theAlignables, outFile.c_str(), loop, false, ioerr);
00643 if (ioerr) {
00644 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO" << "Problem " << ioerr
00645 << " in writeAlignableRelativePositions, " << loop;
00646 ++result;
00647 }
00648
00649 return result;
00650 }
00651
00652
00653 void MillePedeAlignmentAlgorithm::buildUserVariables(const std::vector<Alignable*> &alis) const
00654 {
00655 for (std::vector<Alignable*>::const_iterator iAli = alis.begin(); iAli != alis.end(); ++iAli) {
00656 AlignmentParameters *params = (*iAli)->alignmentParameters();
00657 if (!params) {
00658 throw cms::Exception("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::buildUserVariables"
00659 << "No parameters for alignable";
00660 }
00661 MillePedeVariables *userVars = new MillePedeVariables(params->size(), thePedeLabels->alignableLabel(*iAli));
00662 params->setUserVariables(userVars);
00663 }
00664 }
00665
00666
00667 unsigned int MillePedeAlignmentAlgorithm::decodeMode(const std::string &mode) const
00668 {
00669 if (mode == "full") {
00670 return myMilleBit + myPedeSteerBit + myPedeRunBit + myPedeReadBit;
00671 } else if (mode == "mille") {
00672 return myMilleBit;
00673 } else if (mode == "pede") {
00674 return myPedeSteerBit + myPedeRunBit + myPedeReadBit;
00675 } else if (mode == "pedeSteer") {
00676 return myPedeSteerBit;
00677 } else if (mode == "pedeRun") {
00678 return myPedeSteerBit + myPedeRunBit + myPedeReadBit;
00679 } else if (mode == "pedeRead") {
00680 return myPedeReadBit;
00681 }
00682
00683 throw cms::Exception("BadConfig")
00684 << "Unknown mode '" << mode
00685 << "', use 'full', 'mille', 'pede', 'pedeRun', 'pedeSteer' or 'pedeRead'.";
00686
00687 return 0;
00688 }
00689
00690
00691 bool MillePedeAlignmentAlgorithm::addHitStatistics(int fromLoop, const std::string &outFile,
00692 const std::vector<std::string> &inFiles) const
00693 {
00694 bool allOk = true;
00695 int ierr = 0;
00696 MillePedeVariablesIORoot millePedeIO;
00697 if (inFiles.empty()) {
00698 const std::vector<AlignmentUserVariables*> mpVars =
00699 millePedeIO.readMillePedeVariables(theAlignables, outFile.c_str(), fromLoop, ierr);
00700 if (ierr || !this->addHits(theAlignables, mpVars)) {
00701 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addHitStatistics"
00702 << "Error " << ierr << " reading from " << outFile
00703 << ", tree " << fromLoop << ", or problems in addHits";
00704 allOk = false;
00705 }
00706 for (std::vector<AlignmentUserVariables*>::const_iterator i = mpVars.begin();
00707 i != mpVars.end(); ++i){
00708 delete *i;
00709 }
00710 } else {
00711 for (std::vector<std::string>::const_iterator iFile = inFiles.begin();
00712 iFile != inFiles.end(); ++iFile) {
00713 const std::string inFile(theDir + *iFile);
00714 const std::vector<AlignmentUserVariables*> mpVars =
00715 millePedeIO.readMillePedeVariables(theAlignables, inFile.c_str(), fromLoop, ierr);
00716 if (ierr || !this->addHits(theAlignables, mpVars)) {
00717 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addHitStatistics"
00718 << "Error " << ierr << " reading from " << inFile
00719 << ", tree " << fromLoop << ", or problems in addHits";
00720 allOk = false;
00721 }
00722 for (std::vector<AlignmentUserVariables*>::const_iterator i = mpVars.begin();
00723 i != mpVars.end(); ++i) {
00724 delete *i;
00725 }
00726 }
00727 }
00728
00729 return allOk;
00730 }
00731
00732
00733 bool MillePedeAlignmentAlgorithm::addHits(const std::vector<Alignable*> &alis,
00734 const std::vector<AlignmentUserVariables*> &mpVars) const
00735 {
00736 bool allOk = (mpVars.size() == alis.size());
00737 std::vector<AlignmentUserVariables*>::const_iterator iUser = mpVars.begin();
00738 for (std::vector<Alignable*>::const_iterator iAli = alis.begin();
00739 iAli != alis.end() && iUser != mpVars.end(); ++iAli, ++iUser) {
00740 MillePedeVariables *mpVarNew = dynamic_cast<MillePedeVariables*>(*iUser);
00741 AlignmentParameters *ps = (*iAli)->alignmentParameters();
00742 MillePedeVariables *mpVarOld = (ps ? dynamic_cast<MillePedeVariables*>(ps->userVariables()) : 0);
00743 if (!mpVarNew || !mpVarOld || mpVarOld->size() != mpVarNew->size()) {
00744 allOk = false;
00745 continue;
00746 }
00747
00748 mpVarOld->increaseHitsX(mpVarNew->hitsX());
00749 mpVarOld->increaseHitsY(mpVarNew->hitsY());
00750 }
00751
00752 return allOk;
00753 }
00754
00755
00756 void MillePedeAlignmentAlgorithm::makeGlobDerivMatrix(const std::vector<float> &globalDerivativesx,
00757 const std::vector<float> &globalDerivativesy,
00758 TMatrixF &aGlobalDerivativesM)
00759 {
00760
00761 for (unsigned int i = 0; i < globalDerivativesx.size(); ++i) {
00762 aGlobalDerivativesM(0,i) = globalDerivativesx[i];
00763 aGlobalDerivativesM(1,i) = globalDerivativesy[i];
00764 }
00765 }
00766
00767
00768 void MillePedeAlignmentAlgorithm::diagonalize
00769 (TMatrixDSym &aHitCovarianceM, TMatrixF &aLocalDerivativesM, TMatrixF &aHitResidualsM,
00770 TMatrixF &aGlobalDerivativesM) const
00771 {
00772 TMatrixDSymEigen myDiag(aHitCovarianceM);
00773 TMatrixD aTranfoToDiagonalSystem = myDiag.GetEigenVectors();
00774 TMatrixD aTranfoToDiagonalSystemInv = myDiag.GetEigenVectors( );
00775 TMatrixF aTranfoToDiagonalSystemInvF = myDiag.GetEigenVectors( );
00776 TMatrixD aMatrix = aTranfoToDiagonalSystemInv.Invert() * aHitCovarianceM * aTranfoToDiagonalSystem;
00777
00778
00779
00780
00781
00782 aHitCovarianceM = TMatrixDSym(2, aMatrix.GetMatrixArray());
00783 aTranfoToDiagonalSystemInvF.Invert();
00784
00785 aLocalDerivativesM = aTranfoToDiagonalSystemInvF * aLocalDerivativesM;
00786
00787
00788 aHitResidualsM = aTranfoToDiagonalSystemInvF * aHitResidualsM;
00789 aGlobalDerivativesM = aTranfoToDiagonalSystemInvF * aGlobalDerivativesM;
00790 }
00791
00792
00793 void MillePedeAlignmentAlgorithm
00794 ::addRefTrackMsMeas1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00795 unsigned int iMsMeas, TMatrixDSym &aHitCovarianceM,
00796 TMatrixF &aHitResidualsM, TMatrixF &aLocalDerivativesM)
00797 {
00798
00799
00800 const unsigned int xIndex = iMsMeas + refTrajPtr->numberOfHitMeas();
00801
00802
00803
00804 aHitCovarianceM(0,0)=refTrajPtr->measurementErrors()[xIndex][xIndex];
00805
00806
00807 aHitResidualsM(0,0)= refTrajPtr->measurements()[xIndex];
00808
00809
00810 const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
00811
00812
00813
00814 for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
00815 aLocalDerivativesM(0,i) = locDerivMatrix[xIndex][i];
00816 }
00817 }
00818
00819
00820 void MillePedeAlignmentAlgorithm
00821 ::addRefTrackData2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00822 unsigned int iTrajHit, TMatrixDSym &aHitCovarianceM,
00823 TMatrixF &aHitResidualsM, TMatrixF &aLocalDerivativesM)
00824 {
00825
00826
00827 const unsigned int xIndex = iTrajHit*2;
00828 const unsigned int yIndex = iTrajHit*2+1;
00829
00830
00831
00832 aHitCovarianceM(0,0)=refTrajPtr->measurementErrors()[xIndex][xIndex];
00833 aHitCovarianceM(0,1)=refTrajPtr->measurementErrors()[xIndex][yIndex];
00834 aHitCovarianceM(1,0)=refTrajPtr->measurementErrors()[yIndex][xIndex];
00835 aHitCovarianceM(1,1)=refTrajPtr->measurementErrors()[yIndex][yIndex];
00836
00837
00838 aHitResidualsM(0,0)= refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
00839 aHitResidualsM(1,0)= refTrajPtr->measurements()[yIndex] - refTrajPtr->trajectoryPositions()[yIndex];
00840
00841
00842 const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
00843
00844
00845
00846 for (int i = 0; i < locDerivMatrix.num_col(); ++i) {
00847 aLocalDerivativesM(0,i) = locDerivMatrix[xIndex][i];
00848 aLocalDerivativesM(1,i) = locDerivMatrix[yIndex][i];
00849 }
00850 }
00851
00852
00853 int MillePedeAlignmentAlgorithm
00854 ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00855 unsigned int iTrajHit, const std::vector<int> &globalLabels,
00856 const std::vector<float> &globalDerivativesX,
00857 const std::vector<float> &globalDerivativesY)
00858 {
00859 const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
00860
00861 if((aRecHit)->dimension() == 1) {
00862 return this->callMille1D(refTrajPtr, iTrajHit, globalLabels, globalDerivativesX);
00863 } else {
00864 return this->callMille2D(refTrajPtr, iTrajHit, globalLabels,
00865 globalDerivativesX, globalDerivativesY);
00866 }
00867 }
00868
00869
00870
00871 int MillePedeAlignmentAlgorithm
00872 ::callMille1D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00873 unsigned int iTrajHit, const std::vector<int> &globalLabels,
00874 const std::vector<float> &globalDerivativesX)
00875 {
00876 const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
00877 const unsigned int xIndex = iTrajHit*2;
00878
00879
00880 const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
00881 const int nLocal = locDerivMatrix.num_col();
00882 std::vector<float> localDerivatives(nLocal);
00883 for (unsigned int i = 0; i < localDerivatives.size(); ++i) {
00884 localDerivatives[i] = locDerivMatrix[xIndex][i];
00885 }
00886
00887
00888 float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
00889 float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
00890
00891
00892 const int nGlobal = globalDerivativesX.size();
00893
00894
00895
00896 theMille->mille(nLocal, &(localDerivatives[0]), nGlobal, &(globalDerivativesX[0]),
00897 &(globalLabels[0]), residX, hitErrX);
00898
00899 if (theMonitor) {
00900 theMonitor->fillDerivatives(aRecHit, &(localDerivatives[0]), nLocal,
00901 &(globalDerivativesX[0]), nGlobal, &(globalLabels[0]));
00902 theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit],
00903 iTrajHit, residX, hitErrX, false);
00904 }
00905
00906 return 1;
00907 }
00908
00909
00910 int MillePedeAlignmentAlgorithm
00911 ::callMille2D(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00912 unsigned int iTrajHit, const std::vector<int> &globalLabels,
00913 const std::vector<float> &globalDerivativesx,
00914 const std::vector<float> &globalDerivativesy)
00915 {
00916 const ConstRecHitPointer aRecHit(refTrajPtr->recHits()[iTrajHit]);
00917
00918 if((aRecHit)->dimension() != 2) {
00919 edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::callMille2D"
00920 << "You try to call method for 2D hits for a "
00921 << (aRecHit)->dimension()
00922 << "D Hit. Hit gets ignored!";
00923 return -1;
00924 }
00925
00926 TMatrixDSym aHitCovarianceM(2);
00927 TMatrixF aHitResidualsM(2,1);
00928 TMatrixF aLocalDerivativesM(2, refTrajPtr->derivatives().num_col());
00929
00930 this->addRefTrackData2D(refTrajPtr, iTrajHit, aHitCovarianceM,aHitResidualsM,aLocalDerivativesM);
00931 TMatrixF aGlobalDerivativesM(2,globalDerivativesx.size());
00932 this->makeGlobDerivMatrix(globalDerivativesx, globalDerivativesy, aGlobalDerivativesM);
00933
00934
00935
00936
00937 const double corr = aHitCovarianceM(0,1) / sqrt(aHitCovarianceM(0,0) * aHitCovarianceM(1,1));
00938 if (theMonitor) theMonitor->fillCorrelations2D(corr, aRecHit);
00939 bool diag = false;
00940 switch(aRecHit->geographicalId().subdetId()) {
00941 case SiStripDetId::TID:
00942 case SiStripDetId::TEC:
00943 if (aRecHit->geographicalId().det() == DetId::Tracker && TMath::Abs(corr) > theMaximalCor2D) {
00944 this->diagonalize(aHitCovarianceM, aLocalDerivativesM, aHitResidualsM, aGlobalDerivativesM);
00945 diag = true;
00946 }
00947 break;
00948 default:;
00949 }
00950
00951 float newResidX = aHitResidualsM(0,0);
00952 float newResidY = aHitResidualsM(1,0);
00953 float newHitErrX = TMath::Sqrt(aHitCovarianceM(0,0));
00954 float newHitErrY = TMath::Sqrt(aHitCovarianceM(1,1));
00955 float *newLocalDerivsX = aLocalDerivativesM[0].GetPtr();
00956 float *newLocalDerivsY = aLocalDerivativesM[1].GetPtr();
00957 float *newGlobDerivsX = aGlobalDerivativesM[0].GetPtr();
00958 float *newGlobDerivsY = aGlobalDerivativesM[1].GetPtr();
00959 const int nLocal = aLocalDerivativesM.GetNcols();
00960 const int nGlobal = aGlobalDerivativesM.GetNcols();
00961
00962 if (diag && (newHitErrX > newHitErrY)) {
00963
00964 std::swap(newResidX, newResidY);
00965 std::swap(newHitErrX, newHitErrY);
00966 std::swap(newLocalDerivsX, newLocalDerivsY);
00967 std::swap(newGlobDerivsX, newGlobDerivsY);
00968 }
00969
00970
00971
00972 theMille->mille(nLocal, newLocalDerivsX, nGlobal, newGlobDerivsX,
00973 &(globalLabels[0]), newResidX, newHitErrX);
00974
00975 if (theMonitor) {
00976 theMonitor->fillDerivatives(aRecHit, newLocalDerivsX, nLocal, newGlobDerivsX, nGlobal,
00977 &(globalLabels[0]));
00978 theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit],
00979 iTrajHit, newResidX, newHitErrX, false);
00980 }
00981 const bool isReal2DHit = this->is2D(aRecHit);
00982 if (isReal2DHit) {
00983 theMille->mille(nLocal, newLocalDerivsY, nGlobal, newGlobDerivsY,
00984 &(globalLabels[0]), newResidY, newHitErrY);
00985 if (theMonitor) {
00986 theMonitor->fillDerivatives(aRecHit, newLocalDerivsY, nLocal, newGlobDerivsY, nGlobal,
00987 &(globalLabels[0]));
00988 theMonitor->fillResiduals(aRecHit, refTrajPtr->trajectoryStates()[iTrajHit],
00989 iTrajHit, newResidY, newHitErrY, true);
00990 }
00991 }
00992
00993 return (isReal2DHit ? 2 : 1);
00994 }
00995
00996
00997 void MillePedeAlignmentAlgorithm
00998 ::addMsMeas(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr, unsigned int iMsMeas)
00999 {
01000 TMatrixDSym aHitCovarianceM(1);
01001 TMatrixF aHitResidualsM(1,1);
01002 TMatrixF aLocalDerivativesM(1, refTrajPtr->derivatives().num_col());
01003
01004 this->addRefTrackMsMeas1D(refTrajPtr, iMsMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
01005
01006
01007 TMatrixF aGlobalDerivativesM(1,1);
01008 aGlobalDerivativesM(0,0) = 0;
01009
01010 float newResidX = aHitResidualsM(0,0);
01011 float newHitErrX = TMath::Sqrt(aHitCovarianceM(0,0));
01012 float *newLocalDerivsX = aLocalDerivativesM[0].GetPtr();
01013 float *newGlobDerivsX = aGlobalDerivativesM[0].GetPtr();
01014 const int nLocal = aLocalDerivativesM.GetNcols();
01015 const int nGlobal = 0;
01016
01017 theMille->mille(nLocal, newLocalDerivsX, nGlobal, newGlobDerivsX,
01018 &nGlobal, newResidX, newHitErrX);
01019 }
01020
01021
01022 void MillePedeAlignmentAlgorithm::addLaserData(const EventInfo &eventInfo,
01023 const TkFittedLasBeamCollection &lasBeams,
01024 const TsosVectorCollection &lasBeamTsoses)
01025 {
01026 TsosVectorCollection::const_iterator iTsoses = lasBeamTsoses.begin();
01027 for(TkFittedLasBeamCollection::const_iterator iBeam = lasBeams.begin(), iEnd = lasBeams.end();
01028 iBeam != iEnd; ++iBeam, ++iTsoses){
01029
01030 edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::addLaserData"
01031 << "Beam " << iBeam->getBeamId() << " with "
01032 << iBeam->parameters().size() << " parameters and "
01033 << iBeam->getData().size() << " hits.\n There are "
01034 << iTsoses->size() << " TSOSes.";
01035
01036 this->addLasBeam(eventInfo, *iBeam, *iTsoses);
01037 }
01038 }
01039
01040
01041 void MillePedeAlignmentAlgorithm::addLasBeam(const EventInfo &eventInfo,
01042 const TkFittedLasBeam &lasBeam,
01043 const std::vector<TrajectoryStateOnSurface> &tsoses)
01044 {
01045 AlignmentParameters *dummyPtr = 0;
01046 std::vector<float> lasLocalDerivsX;
01047 const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId());
01048
01049 for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
01050 if (!tsoses[iHit].isValid()) continue;
01051
01052 theFloatBufferX.clear();
01053 theFloatBufferY.clear();
01054 theIntBuffer.clear();
01055 lasLocalDerivsX.clear();
01056
01057 const SiStripLaserRecHit2D &hit = lasBeam.getData()[iHit];
01058 AlignableDetOrUnitPtr lasAli(theAlignableNavigator->alignableFromDetId(hit.getDetId()));
01059 this->globalDerivativesHierarchy(eventInfo,
01060 tsoses[iHit], lasAli, lasAli,
01061 theFloatBufferX, theFloatBufferY, theIntBuffer, dummyPtr);
01062
01063 for (unsigned int nFitParams = 0;
01064 nFitParams < static_cast<unsigned int>(lasBeam.parameters().size());
01065 ++nFitParams) {
01066 const float derivative = lasBeam.derivatives()[iHit][nFitParams];
01067 if (nFitParams < lasBeam.firstFixedParameter()) {
01068 lasLocalDerivsX.push_back(derivative);
01069 } else {
01070 const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
01071 theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
01072 theFloatBufferX.push_back(derivative);
01073 }
01074 }
01075
01076 const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
01077
01078 const float error = 0.003;
01079
01080 theMille->mille(lasLocalDerivsX.size(), &(lasLocalDerivsX[0]), theFloatBufferX.size(),
01081 &(theFloatBufferX[0]), &(theIntBuffer[0]), residual, error);
01082 }
01083
01084 theMille->end();
01085 }
01086
01087 void MillePedeAlignmentAlgorithm::addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg)
01088 {
01089
01090 const bool doOutputOnStdout(pxbSurveyCfg.getParameter<bool>("doOutputOnStdout"));
01091 if (doOutputOnStdout) std::cout << "# Output from addPxbSurvey follows below because doOutputOnStdout is set to True" << std::endl;
01092
01093
01094 SurveyPxbDicer dicer(pxbSurveyCfg.getParameter<std::vector<edm::ParameterSet> >("toySurveyParameters"), pxbSurveyCfg.getParameter<unsigned int>("toySurveySeed"));
01095 std::ofstream outfile(pxbSurveyCfg.getUntrackedParameter<std::string>("toySurveyFile").c_str());
01096
01097
01098 std::vector<SurveyPxbImageLocalFit> measurements;
01099 std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
01100 SurveyPxbImageReader<SurveyPxbImageLocalFit> reader(filename, measurements, 800);
01101
01102
01103 for(std::vector<SurveyPxbImageLocalFit>::size_type i=0; i!=measurements.size(); i++)
01104 {
01105 if (doOutputOnStdout) std::cout << "Module " << i << ": ";
01106
01107
01108 AlignableDetOrUnitPtr mod1(theAlignableNavigator->alignableFromDetId(measurements[i].getIdFirst()));
01109 AlignableDetOrUnitPtr mod2(theAlignableNavigator->alignableFromDetId(measurements[i].getIdSecond()));
01110 const AlignableSurface& surf1 = mod1->surface();
01111 const AlignableSurface& surf2 = mod2->surface();
01112
01113
01114 const LocalPoint fidpoint0(-0.91,+3.30);
01115 const LocalPoint fidpoint1(+0.91,+3.30);
01116 const LocalPoint fidpoint2(+0.91,-3.30);
01117 const LocalPoint fidpoint3(-0.91,-3.30);
01118
01119
01120
01121
01122 const GlobalPoint surf2point0(surf2.toGlobal(fidpoint0));
01123 const GlobalPoint surf2point1(surf2.toGlobal(fidpoint1));
01124 const LocalPoint fidpoint0inSurf1frame(surf1.toLocal(surf2point0));
01125 const LocalPoint fidpoint1inSurf1frame(surf1.toLocal(surf2point1));
01126
01127
01128 SurveyPxbImageLocalFit::fidpoint_t fidpointvec;
01129 fidpointvec.push_back(fidpoint0inSurf1frame);
01130 fidpointvec.push_back(fidpoint1inSurf1frame);
01131 fidpointvec.push_back(fidpoint2);
01132 fidpointvec.push_back(fidpoint3);
01133
01134
01135 if (pxbSurveyCfg.getParameter<bool>("doToySurvey"))
01136 {
01137 dicer.doDice(fidpointvec,measurements[i].getIdPair(), outfile);
01138 }
01139
01140
01141 measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
01142 SurveyPxbImageLocalFit::localpars_t a;
01143 a = measurements[i].getLocalParameters();
01144 const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
01145
01146
01147 if (doOutputOnStdout)
01148 {
01149 std::cout << "a: " << a[0] << ", " << a[1] << ", " << a[2] << ", " << a[3]
01150 << " S= " << sqrt(a[2]*a[2]+a[3]*a[3])
01151 << " phi= " << atan(a[3]/a[2])
01152 << " chi2= " << chi2 << std::endl;
01153 }
01154 if (theMonitor)
01155 {
01156 theMonitor->fillPxbSurveyHistsChi2(chi2);
01157 theMonitor->fillPxbSurveyHistsLocalPars(a[0],a[1],sqrt(a[2]*a[2]+a[3]*a[3]),atan(a[3]/a[2]));
01158 }
01159
01160
01161 for(SurveyPxbImageLocalFit::count_t j=0; j!=SurveyPxbImageLocalFit::nMsrmts; j++)
01162 {
01163 theMille->mille((int)measurements[i].getLocalDerivsSize(),
01164 measurements[i].getLocalDerivsPtr(j),
01165 (int)measurements[i].getGlobalDerivsSize(),
01166 measurements[i].getGlobalDerivsPtr(j),
01167 measurements[i].getGlobalDerivsLabelPtr(j),
01168 measurements[i].getResiduum(j),
01169 measurements[i].getSigma(j));
01170 }
01171 theMille->end();
01172 }
01173 outfile.close();
01174 }
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271