CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Alignment/MillePedeAlignmentAlgorithm/src/MillePedeAlignmentAlgorithm.cc

Go to the documentation of this file.
00001 
00011 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeAlignmentAlgorithm.h"
00012 //#include "MillePedeAlignmentAlgorithm.h"
00013 
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 
00016 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00017 // in header, too
00018 // end in header, too
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"       // 'unpublished' interface located in src
00024 #include "PedeSteerer.h" // dito
00025 #include "PedeReader.h" // dito
00026 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h" // dito
00027 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerPluginFactory.h" // dito
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 // includes to make known that they inherit from Alignable:
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 // Includes for PXB survey
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 // Constructor ----------------------------------------------------------------
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 += '/';// may need '/'
00089   edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm" << "Start in mode '"
00090                             << theConfig.getUntrackedParameter<std::string>("mode")
00091                             << "' with output directory '" << theDir << "'.";
00092 }
00093 
00094 // Destructor ----------------------------------------------------------------
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 // Call at beginning of job ---------------------------------------------------
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   // 1) Create PedeSteerer: correct alignable positions for coordinate system selection
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   // 2) If requested, directly read in and apply result of previous pede run,
00169   //    assuming that correction from 1) was also applied to create the result:
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)) { // false: do not erase SelectionUserVariables
00183       throw cms::Exception("BadConfig")
00184         << "MillePedeAlignmentAlgorithm::initialize: Problems reading input constants of "
00185         << "pedeReaderInputs entry " << iSet - mprespset.begin() << '.';
00186     }
00187     theAlignmentParameterStore->applyParameters();
00188     // Needed to shut up later warning from checkAliParams:
00189     theAlignmentParameterStore->resetParameters();
00190   }
00191   
00192   // 3) Now create steerings with 'final' start position:
00193   thePedeSteer->buildSubSteer(tracker, muon, extras);
00194 
00195   // After (!) 1-3 of PedeSteerer which uses the SelectionUserVariables attached to the parameters:
00196   this->buildUserVariables(theAlignables); // for hit statistics and/or pede result
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());// add ', false);' for text output);
00206     const std::string moniFile(theConfig.getUntrackedParameter<std::string>("monitorFile"));
00207     if (moniFile.size()) theMonitor = new MillePedeMonitor((theDir + moniFile).c_str());
00208 
00209     // Get trajectory factory. In case nothing found, FrameWork will throw...
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   // FIXME: for PlotMillePede hit statistics stuff we also might want doIO(0)... ?
00216   if (this->isMode(myPedeSteerBit) // for pedeRun and pedeRead we might want to merge
00217       || !theConfig.getParameter<std::vector<std::string> >("mergeTreeFiles").empty()) {
00218     this->doIO(0);
00219 
00220   // Get config for survey and set flag accordingly
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 // Call at end of job ---------------------------------------------------------
00239 //____________________________________________________
00240 void MillePedeAlignmentAlgorithm::terminate()
00241 {
00242   delete theMille;// delete to close binary before running pede below (flush would be enough...)
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));// do only if myPedeSteerBit?
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)) { // if mille was run, we store trees with suffix _1...
00272     this->doIO(1);
00273   } else if (this->isMode(myPedeReadBit)) {// if pede runs otherwise, we use _2 (=> probably merge)
00274     this->doIO(2);
00275   }
00276 }
00277 
00278 // Run the algorithm on trajectories and tracks -------------------------------
00279 //____________________________________________________
00280 void MillePedeAlignmentAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo)
00281 {
00282   if (!this->isMode(myMilleBit)) return; // no theMille created...
00283   const ConstTrajTrackPairCollection &tracks = eventInfo.trajTrackPairs_;
00284 
00285   if (theMonitor) { // monitor input tracks
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   // Now loop over ReferenceTrajectoryCollection
00295   unsigned int refTrajCount = 0; // counter for track monitoring if 1 track per trajectory
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       // if track used (i.e. some hits), fill monitoring
00306       // track NULL ptr if trajectories and tracks do not match
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   } // end of reference trajectory and track loop
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     // to add hits if all fine:
00326     std::vector<AlignmentParameters*> parVec(refTrajPtr->recHits().size());
00327     // collect hit statistics, assuming that there are no y-only hits
00328     std::vector<bool> validHitVecY(refTrajPtr->recHits().size(), false);
00329     // Use recHits from ReferenceTrajectory (since they have the right order!):
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) { // problem
00334         hitResultXy.first = 0;
00335         break;
00336       } else { // hit is fine, increase x/y statistics
00337         if (flagXY >= 1) ++hitResultXy.first;
00338         validHitVecY[iHit] = (flagXY >= 2);
00339       } 
00340     } // end loop on hits
00341 
00342 // CHK add 'Multiple Scattering Measurements' for break points, broken lines
00343     for (unsigned int iMsMeas = 0; iMsMeas < refTrajPtr->numberOfMsMeas(); ++iMsMeas) {
00344       this->addMsMeas(refTrajPtr, iMsMeas);
00345     }
00346              
00347     // kill or end 'track' for mille, depends on #hits criterion
00348     if (hitResultXy.first == 0 || hitResultXy.first < theMinNumHits) {
00349       theMille->kill();
00350       hitResultXy.first = hitResultXy.second = 0; //reset
00351     } else {
00352       theMille->end();
00353       // take care about hit statistics as well
00354       for (unsigned int iHit = 0; iHit < validHitVecY.size(); ++iHit) {
00355         if (!parVec[iHit]) continue; // in case a non-selected alignable was hit (flagXY == 0)
00356         MillePedeVariables *mpVar = static_cast<MillePedeVariables*>(parVec[iHit]->userVariables());
00357         mpVar->increaseHitsX(); // every hit has an x-measurement, cf. above...
00358         if (validHitVecY[iHit]) {
00359           mpVar->increaseHitsY();
00360           ++hitResultXy.second;
00361         }
00362       }  
00363     }
00364   } // end if valid trajectory
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     // LAS beam treatment
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 *&params)
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   // ignore invalid hits
00393   if (!recHitPtr->isValid()) return 0;
00394 
00395   // get AlignableDet/Unit for this hit
00396   AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(recHitPtr->geographicalId()));
00397   
00398   if (!this->globalDerivativesHierarchy(eventInfo,
00399                                         tsos, alidet, alidet, theFloatBufferX, // 2x alidet, sic!
00400                                         theFloatBufferY, theIntBuffer, params)) {
00401     return -1; // problem
00402   } else if (theFloatBufferX.empty()) {
00403     return 0; // empty for X: no alignable for hit
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   // derivatives and labels are recursively attached
00420   if (!ali) return true; // no mother might be OK
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; // set parameters of lowest level
00428 
00429     bool hasSplitParameters = thePedeLabels->hasSplitParameters(ali);
00430     const unsigned int alignableLabel = thePedeLabels->alignableLabel(ali);
00431     
00432     if (0 == alignableLabel) { // FIXME: what about regardAllHits in Markus' code?
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     // cols: 2, i.e. x&y, rows: parameters, usually RigidBodyAlignmentParameters::N_PARAM
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     // Exclude mothers if Alignable selected to be no part of a hierarchy:
00456     if (thePedeSteer->isNoHiera(ali)) return true;
00457   }
00458   // Call recursively for mother, will stop if mother == 0:
00459   return this->globalDerivativesHierarchy(eventInfo,
00460                                           tsos, ali->mother(), alidet,
00461                                           globalDerivativesX, globalDerivativesY,
00462                                           globalLabels, lowestParams);
00463 }
00464 
00465 // //____________________________________________________
00466 // void MillePedeAlignmentAlgorithm
00467 // ::callMille(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
00468 //             unsigned int iTrajHit, MeasurementDirection xOrY,
00469 //             const std::vector<float> &globalDerivatives, const std::vector<int> &globalLabels)
00470 // {
00471 //   const unsigned int xyIndex = iTrajHit*2 + xOrY;
00472 //   // FIXME: here for residuum and sigma we could use KALMAN-Filter results
00473 //   const float residuum =
00474 //     refTrajPtr->measurements()[xyIndex] - refTrajPtr->trajectoryPositions()[xyIndex];
00475 //   const float covariance = refTrajPtr->measurementErrors()[xyIndex][xyIndex];
00476 //   const float sigma = (covariance > 0. ? TMath::Sqrt(covariance) : 0.);
00477 
00478 //   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
00479 
00480 //   std::vector<float> localDerivs(locDerivMatrix.num_col());
00481 //   for (unsigned int i = 0; i < localDerivs.size(); ++i) {
00482 //     localDerivs[i] = locDerivMatrix[xyIndex][i];
00483 //   }
00484 
00485 //   // &(vector[0]) is valid - as long as vector is not empty 
00486 //   // cf. https://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
00487 //   theMille->mille(localDerivs.size(), &(localDerivs[0]),
00488 //                globalDerivatives.size(), &(globalDerivatives[0]), &(globalLabels[0]),
00489 //                residuum, sigma);
00490 //   if (theMonitor) {
00491 //     theMonitor->fillDerivatives(refTrajPtr->recHits()[iTrajHit],localDerivs, globalDerivatives,
00492 //                              (xOrY == kLocalY));
00493 //     theMonitor->fillResiduals(refTrajPtr->recHits()[iTrajHit],
00494 //                            refTrajPtr->trajectoryStates()[iTrajHit],
00495 //                            iTrajHit, residuum, sigma, (xOrY == kLocalY));
00496 //   }
00497 // }
00498 
00499 //____________________________________________________
00500 bool MillePedeAlignmentAlgorithm::is2D(const ConstRecHitPointer &recHit) const
00501 {
00502   // FIXME: Check whether this is a reliable and recommended way to find out...
00503 
00504   if (recHit->dimension() < 2) {
00505     return false; // some muon and TIB/TOB stuff really has RecHit1D
00506   } else if (recHit->detUnit()) { // detunit in strip is 1D, in pixel 2D 
00507     return recHit->detUnit()->type().isTrackerPixel();
00508   } else { // stereo strips  (FIXME: endcap trouble due to non-parallel strips (wedge sensors)?)
00509     if (dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit->hit())) { // check persistent hit
00510       // projected: 1D measurement on 'glued' module
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; // FIXME: Should we check one by one? Or transfer 'alis' to the store?
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) { // as many alignables with result as trying to align
00541       edm::LogInfo("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00542     } else if (alis.size()) { // dead module do not get hits and no pede result
00543       edm::LogWarning("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00544     } else { // serious problem: no result read - and not all modules can be dead...
00545       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::readFromPede" << out.str();
00546       return false;
00547     }
00548     return true;
00549   }
00550   // the rest is not OK:
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 // // problem with following writeOrigRigidBodyAlignmentParameters
00621 //     aliIO.writeAlignmentParameters(theAlignables, outFile.c_str(), loop, false, ioerr);
00622 //     if (ioerr) {
00623 //       edm::LogError("Alignment") << "@SUB=MillePedeAlignmentAlgorithm::doIO"
00624 //                                  << "Problem " << ioerr << " in writeAlignmentParameters, " << loop;
00625 //       ++result;
00626 //     }
00627   }
00628   
00629   //  aliIO.writeAlignmentParameters(theAlignables, ("tmp"+outFile).c_str(), loop, false, ioerr);
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; // + myPedeSteerBit; // sic! Including production of steerig file. NO!
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; // sic! Including steering and reading of result.
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; // clean created objects
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; // clean created objects
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; // FIXME error etc.?
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   // Tranformation of matrix M is done by A^T*M*A, not A^{-1}*M*A.
00778   // But here A^T == A^{-1}, so we would only save CPU by Transpose()...
00779   // FIXME this - I guess simply use T(), not Transpose()...
00780   // TMatrixD aMatrix = aTranfoToDiagonalSystemInv.Transpose() * aHitCovarianceM
00781   //    * aTranfoToDiagonalSystem;
00782   aHitCovarianceM = TMatrixDSym(2, aMatrix.GetMatrixArray());
00783   aTranfoToDiagonalSystemInvF.Invert();
00784   //edm::LogInfo("Alignment") << "NEW HIT loca in matrix"<<aLocalDerivativesM(0,0);
00785   aLocalDerivativesM = aTranfoToDiagonalSystemInvF * aLocalDerivativesM;
00786   
00787   //edm::LogInfo("Alignment") << "NEW HIT loca in matrix after diag:"<<aLocalDerivativesM(0,0);
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   // This Method is valid for 1D measurements only
00800   const unsigned int xIndex = iMsMeas + refTrajPtr->numberOfHitMeas();
00801   // Covariance into a TMatrixDSym
00802   
00803   //aHitCovarianceM = new TMatrixDSym(1);
00804   aHitCovarianceM(0,0)=refTrajPtr->measurementErrors()[xIndex][xIndex];
00805   
00806   //theHitResidualsM= new TMatrixF(1,1);
00807   aHitResidualsM(0,0)= refTrajPtr->measurements()[xIndex];
00808   
00809   // Local Derivatives into a TMatrixDSym (to use matrix operations)
00810   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
00811   //  theLocalDerivativeNumber = locDerivMatrix.num_col();
00812   
00813   //theLocalDerivativesM = new TMatrixF(1,locDerivMatrix.num_col());
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   // This Method is valid for 2D measurements only
00826   
00827   const unsigned int xIndex = iTrajHit*2;
00828   const unsigned int yIndex = iTrajHit*2+1;
00829   // Covariance into a TMatrixDSym
00830 
00831   //aHitCovarianceM = new TMatrixDSym(2);
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   //theHitResidualsM= new TMatrixF(2,1);
00838   aHitResidualsM(0,0)= refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
00839   aHitResidualsM(1,0)= refTrajPtr->measurements()[yIndex] - refTrajPtr->trajectoryPositions()[yIndex];
00840   
00841   // Local Derivatives into a TMatrixDSym (to use matrix operations)
00842   const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
00843   //  theLocalDerivativeNumber = locDerivMatrix.num_col();
00844   
00845   //theLocalDerivativesM = new TMatrixF(2,locDerivMatrix.num_col());
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; // the even ones are local x
00878 
00879   // local derivatives
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   // residuum and error
00888   float residX = refTrajPtr->measurements()[xIndex] - refTrajPtr->trajectoryPositions()[xIndex];
00889   float hitErrX = TMath::Sqrt(refTrajPtr->measurementErrors()[xIndex][xIndex]);
00890 
00891   // number of global derivatives
00892   const int nGlobal = globalDerivativesX.size();
00893 
00894   // &(localDerivatives[0]) etc. are valid - as long as vector is not empty
00895   // cf. https://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
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   // below method fills above 3 matrices
00930   this->addRefTrackData2D(refTrajPtr, iTrajHit, aHitCovarianceM,aHitResidualsM,aLocalDerivativesM);
00931   TMatrixF aGlobalDerivativesM(2,globalDerivativesx.size());
00932   this->makeGlobDerivMatrix(globalDerivativesx, globalDerivativesy, aGlobalDerivativesM);
00933  
00934   // calculates correlation between Hit measurements
00935   // FIXME: Should take correlation (and resulting transformation) from original hit, 
00936   //        not 2x2 matrix from ReferenceTrajectory: That can come from error propagation etc.!
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; // diagonalise only tracker TID, TEC
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)) { // also for 2D hits?
00963     // measurement with smaller error is x-measurement (for !is2D do not fill y-measurement):
00964     std::swap(newResidX, newResidY);
00965     std::swap(newHitErrX, newHitErrY);
00966     std::swap(newLocalDerivsX, newLocalDerivsY);
00967     std::swap(newGlobDerivsX, newGlobDerivsY);
00968   }
00969 
00970   // &(globalLabels[0]) is valid - as long as vector is not empty 
00971   // cf. https://www.parashift.com/c++-faq-lite/containers.html#faq-34.3
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); // strip is 1D (except matched hits)
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);// true: y
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   // below method fills above 3 'matrices'
01004   this->addRefTrackMsMeas1D(refTrajPtr, iMsMeas, aHitCovarianceM, aHitResidualsM, aLocalDerivativesM);
01005   
01006   // no global parameters (use dummy 0)
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){ // beam/tsoses parallel!
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; // for globalDerivativesHierarchy()
01046   std::vector<float> lasLocalDerivsX; // buffer for local derivatives
01047   const unsigned int beamLabel = thePedeLabels->lasBeamLabel(lasBeam.getBeamId());// for global par
01048   
01049   for (unsigned int iHit = 0; iHit < tsoses.size(); ++iHit) {
01050     if (!tsoses[iHit].isValid()) continue;
01051     // clear buffer
01052     theFloatBufferX.clear();
01053     theFloatBufferY.clear();
01054     theIntBuffer.clear();
01055     lasLocalDerivsX.clear();
01056     // get alignables and global parameters
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     // fill derivatives vector from derivatives matrix
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()) { // first local beam parameters
01068         lasLocalDerivsX.push_back(derivative);
01069       } else {                                          // now global ones
01070         const unsigned int numPar = nFitParams - lasBeam.firstFixedParameter();
01071         theIntBuffer.push_back(thePedeLabels->parameterLabel(beamLabel, numPar));
01072         theFloatBufferX.push_back(derivative);
01073       }
01074     } // end loop over parameters
01075 
01076     const float residual = hit.localPosition().x() - tsoses[iHit].localPosition().x();
01077     // error from file or assume 0.003
01078     const float error = 0.003; // hit.localPositionError().xx(); sqrt???
01079     
01080     theMille->mille(lasLocalDerivsX.size(), &(lasLocalDerivsX[0]), theFloatBufferX.size(),
01081                     &(theFloatBufferX[0]), &(theIntBuffer[0]), residual, error);
01082   } // end of loop over hits
01083   
01084   theMille->end();
01085 }
01086 
01087 void MillePedeAlignmentAlgorithm::addPxbSurvey(const edm::ParameterSet &pxbSurveyCfg)
01088 {
01089         // do some printing, if requested
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         // instantiate a dicer object
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         // read data from file
01098         std::vector<SurveyPxbImageLocalFit> measurements;
01099         std::string filename(pxbSurveyCfg.getParameter<edm::FileInPath>("infile").fullPath());
01100         SurveyPxbImageReader<SurveyPxbImageLocalFit> reader(filename, measurements, 800);
01101 
01102         // loop over photographs (=measurements) and perform the fit
01103         for(std::vector<SurveyPxbImageLocalFit>::size_type i=0; i!=measurements.size(); i++)
01104         {
01105                 if (doOutputOnStdout) std::cout << "Module " << i << ": ";
01106 
01107                 // get the Alignables and their surfaces
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                 // the position of the fiducial points in local frame of a PXB module
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                 // We choose the local frame of the first module as reference,
01120                 // so take the fidpoints of the second module and calculate their
01121                 // positions in the reference frame
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                 // Create the vector for the fit
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                 // if toy survey is requested, dice the values now
01135                 if (pxbSurveyCfg.getParameter<bool>("doToySurvey"))
01136                 {
01137                         dicer.doDice(fidpointvec,measurements[i].getIdPair(), outfile);
01138                 }
01139                 
01140                 // do the fit
01141                 measurements[i].doFit(fidpointvec, thePedeLabels->alignableLabel(mod1), thePedeLabels->alignableLabel(mod2));
01142             SurveyPxbImageLocalFit::localpars_t a; // local pars from fit
01143                 a = measurements[i].getLocalParameters();
01144                 const SurveyPxbImageLocalFit::value_t chi2 = measurements[i].getChi2();
01145 
01146                 // do some reporting, if requested
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                 // pass the results from the local fit to mille
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 void MillePedeAlignmentAlgorithm::addBeamSpotConstraint
01179 (const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr,
01180  const edm::Event & iEvent )
01181 {
01182   // --- get BS from Event:
01183   const  reco::BeamSpot* bSpot;
01184   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
01185   iEvent.getByType(recoBeamSpotHandle);
01186   bSpot = recoBeamSpotHandle.product();
01187   
01188   //GlobalPoint gPointBs(0.,0.,0.);
01189   GlobalPoint gPointBs(bSpot->x0(), bSpot->y0(), bSpot->z0());
01190   const TrajectoryStateOnSurface trackTsos = refTrajPtr->trajectoryStates()[0];
01191   // create a FTS from innermost TSOS:
01192   FreeTrajectoryState innerFts = *(trackTsos.freeTrajectoryState());
01193   //create a TrajectoryStateClosestToBeamLine: 
01194   TrajectoryStateClosestToPointBuilder *tsctpBuilder = new TSCPBuilderNoMaterial();
01195   TrajectoryStateClosestToPoint tsctp = tsctpBuilder->operator()(trackTsos,gPointBs);
01196   FreeTrajectoryState pcaFts = tsctp.theState();
01197   edm::LogInfo("CHK") << " beamspot TSCP " << tsctp.referencePoint() << tsctp.position()
01198    << tsctp.perigeeParameters().vector();   
01199   const AlgebraicVector5 perigeeVecPars =  tsctp.perigeeParameters().vector(); 
01200   
01201   PerigeeConversions perigeeConv;
01202   const AlgebraicMatrix55& curv2perigee = perigeeConv.jacobianCurvilinear2Perigee(pcaFts);
01203   edm::LogInfo("CHK") << " beamspot C2P " << curv2perigee;
01204   
01205   //propagation
01206   AnalyticalPropagator propagator(&(innerFts.parameters().magneticField()), anyDirection);
01207   std::pair< TrajectoryStateOnSurface, double > tsosWithPath = propagator.propagateWithPath(pcaFts,trackTsos.surface());
01208   edm::LogInfo("CHK") << " beamspot s0 " << tsosWithPath.second;
01209   edm::LogInfo("CHK") << " beamspot t2c " << refTrajPtr->trajectoryToCurv();
01210   if (!tsosWithPath.first.isValid())  return; 
01211   
01212   // jacobian in curvilinear frame for propagation from the end point (inner TSOS) to the starting point (PCA) 
01213   AnalyticalCurvilinearJacobian curvJac( pcaFts.parameters(),
01214                                          tsosWithPath.first.globalPosition(),
01215                                          tsosWithPath.first.globalMomentum(), 
01216                                          tsosWithPath.second );
01217   int ierr;
01218   const AlgebraicMatrix55& matCurvJac = curvJac.jacobian().Inverse(ierr);
01219   edm::LogInfo("CHK") << " beamspot CurvJac " << matCurvJac;
01220   // jacobion trajectory to curvilinear
01221   // const AlgebraicMatrix &locDerivMatrix = refTrajPtr->derivatives();
01222   const AlgebraicMatrix55& traj2curv = asSMatrix<5,5>(refTrajPtr->trajectoryToCurv());   
01223   //combine the transformation jacobian:
01224   AlgebraicMatrix55 newJacobian = curv2perigee * matCurvJac * traj2curv;
01225   edm::LogInfo("CHK") << " beamspot newJac " << newJacobian;
01226 
01227   //get the Beam Spot residual
01228   const float residuumIP = -perigeeVecPars[3];
01229   const float sigmaIP = bSpot->BeamWidth();
01230   edm::LogInfo("CHK2") << " beamspot-res " << residuumIP << " " << perigeeVecPars[0] << " " << perigeeVecPars[1] << " " << perigeeVecPars[2];
01231       
01232   std::vector<float> ipLocalDerivs(5);
01233   for (unsigned int i = 0; i < ipLocalDerivs.size(); ++i) {
01234     ipLocalDerivs[i] = newJacobian(3,i);
01235   }
01236 
01237   // to be fixed: null global derivatives is right but the size is detemined from SelPar.size! 
01238   std::vector<float>  ipGlobalDerivs(4);
01239   for (unsigned int i = 0; i < 4; ++i) {
01240     ipGlobalDerivs[i] = 0.;
01241   }
01242   std::vector<int> ipGlobalLabels(4);
01243   for (unsigned int i = 0; i < 4; ++i) {
01244     ipGlobalLabels[i] = 0;
01245   }
01246 
01247   if(theAliBeamspot){
01248   double phi = perigeeVecPars[2];
01249   double dz  = perigeeVecPars[4];
01250   ipGlobalDerivs[0] = sin(phi);
01251   ipGlobalDerivs[1] = -cos(phi);
01252   ipGlobalDerivs[2] = sin(phi)*dz;
01253   ipGlobalDerivs[3] = -cos(phi)*dz;
01254  
01255   ipGlobalLabels[0] = 250000;
01256   ipGlobalLabels[1] = 250001;
01257   ipGlobalLabels[2] = 250002;
01258   ipGlobalLabels[3] = 250003;
01259   }
01260 
01261 
01262   theMille->mille(ipLocalDerivs.size(), &(ipLocalDerivs[0]),
01263           ipGlobalDerivs.size(), &(ipGlobalDerivs[0]), &(ipGlobalLabels[0]),
01264           residuumIP, sigmaIP);
01265 
01266 
01267   // delete new objects:
01268     delete tsctpBuilder;
01269     tsctpBuilder = NULL;
01270      
01271 } */