Newer
Older
TB_Chris / TbAlignment / src / .svn / text-base / TbAlignmentMinuit0.cpp.svn-base
  1. #include "TbAlignmentMinuit0.h"
  2.  
  3. DECLARE_TOOL_FACTORY(TbAlignmentMinuit0)
  4.  
  5. //=============================================================================
  6. // Standard constructor, initializes variables
  7. //=============================================================================
  8. TbAlignmentMinuit0::TbAlignmentMinuit0(const std::string& type,
  9. const std::string& name,
  10. const IInterface* parent)
  11. : TbAlignmentMinuitBase(type, name, parent) {
  12.  
  13. declareProperty("ReferencePlane", m_referencePlane = 999);
  14.  
  15. m_dofsDefault = {true, true, false, true, true, true};
  16. }
  17.  
  18. //=============================================================================
  19. // Destructor
  20. //=============================================================================
  21. TbAlignmentMinuit0::~TbAlignmentMinuit0() {}
  22.  
  23. //=============================================================================
  24. // Calculate the chi2.
  25. //=============================================================================
  26. void TbAlignmentMinuit0::chi2(double& f, double* par, double* /*g*/) {
  27.  
  28. // Assign new aligment constants
  29. for (auto im = m_modules.begin(), end = m_modules.end(); im != end; ++im) {
  30. const unsigned int i = im - m_modules.begin();
  31. (*im)->setAlignment(par[6 * i + 0], par[6 * i + 1], par[6 * i + 2],
  32. par[6 * i + 3], par[6 * i + 4], par[6 * i + 5]);
  33. }
  34.  
  35. // Loop over tracks
  36. f = 0.;
  37. for (auto it = m_tracks.begin(), end = m_tracks.end(); it != end; ++it) {
  38. // Update global coordinates of the clusters
  39. SmartRefVector<LHCb::TbCluster> clusters = (*it)->track()->clusters();
  40. for (auto ic = clusters.begin(), end = clusters.end(); ic != end; ++ic) {
  41. Gaudi::XYZPoint pLocal((*ic)->xloc(), (*ic)->yloc(), 0.);
  42. const unsigned int plane = (*ic)->plane();
  43. Gaudi::XYZPoint pGlobal = geomSvc()->localToGlobal(pLocal, plane);
  44. (*ic)->setX(pGlobal.x());
  45. (*ic)->setY(pGlobal.y());
  46. (*ic)->setZ(pGlobal.z());
  47. }
  48. // Refit track with new cluster positions
  49. m_trackFit->fit((*it)->track());
  50. // Add the new track chi2 to the overall chi2
  51. f += (*it)->track()->chi2();
  52. }
  53. }
  54. //=============================================================================
  55. // Main alignment function.
  56. //=============================================================================
  57. void TbAlignmentMinuit0::align(
  58. std::vector<TbAlignmentTrack*>& alignmentTracks) {
  59.  
  60. TbAlignmentMinuitBase::align(alignmentTracks);
  61. info() << "Minuit technique 0" << endmsg;
  62. double arglist[2];
  63. arglist[0] = 10000;
  64. arglist[1] = 1.e-2;
  65.  
  66. for (unsigned int iteration = 0; iteration < m_nIterations; ++iteration) {
  67. info() << "Iteration " << iteration + 1 << "/" << m_nIterations << endmsg;
  68. // Align detector modules one at a time
  69. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  70. // Skip reference plane and masked planes.
  71. if (i == m_referencePlane || masked(i)) continue;
  72. // Wobble this plane and fix the others.
  73. for (unsigned int j = 0; j < m_nPlanes; ++j) {
  74. if (i != j) {
  75. m_fitter->FixParameter(6 * j + 0); // x
  76. m_fitter->FixParameter(6 * j + 1); // y
  77. m_fitter->FixParameter(6 * j + 2); // z
  78. m_fitter->FixParameter(6 * j + 3); // Rx
  79. m_fitter->FixParameter(6 * j + 4); // Ry
  80. m_fitter->FixParameter(6 * j + 5); // Rz
  81. } else {
  82. info() << "*** Wobbling detector " << j << endmsg;
  83. for (unsigned int k = 0; k < 6; ++k) {
  84. if (m_dofs[k]) {
  85. m_fitter->ReleaseParameter(6 * j + k);
  86. } else {
  87. m_fitter->FixParameter(6 * j + k);
  88. }
  89. }
  90. }
  91. }
  92. // Execute minimization and calculate proper error matrix
  93. m_fitter->ExecuteCommand("MIGRAD", arglist, 2);
  94. m_fitter->ExecuteCommand("HESSE", arglist, 1);
  95. if (msgLevel(MSG::INFO))
  96. m_fitter->ExecuteCommand("SHOW PARAMETERS", 0, 0);
  97. }
  98. }
  99. updateGeometry();
  100. }