Newer
Older
TB_Chris / TbAlignment / src / .svn / text-base / TbAlignmentBase.cpp.svn-base
  1. // ROOT
  2. #include <TH1D.h>
  3. #include <TProfile.h>
  4. #include <TFitResult.h>
  5. #include <TFitResultPtr.h>
  6.  
  7. // AIDA
  8. #include "AIDA/IHistogram1D.h"
  9. #include "AIDA/IProfile1D.h"
  10.  
  11. // Gaudi
  12. #include "GaudiUtils/Aida2ROOT.h"
  13.  
  14. // Local
  15. #include "TbAlignmentBase.h"
  16.  
  17. //=============================================================================
  18. // Standard constructor, initializes variables
  19. //=============================================================================
  20. TbAlignmentBase::TbAlignmentBase(const std::string& type,
  21. const std::string& name,
  22. const IInterface* parent)
  23. : GaudiHistoTool(type, name, parent),
  24. m_trackFit(nullptr),
  25. m_geomSvc(nullptr) {
  26.  
  27. declareInterface<TbAlignmentBase>(this);
  28.  
  29. declareProperty("TrackLocation",
  30. m_trackLocation = LHCb::TbTrackLocation::Default);
  31. declareProperty("DOFs", m_dofs = {});
  32. declareProperty("ClearTracks", m_clearTracks = true);
  33. declareProperty("MaskedPlanes", m_maskedPlanes = {});
  34. declareProperty("Monitoring", m_monitoring = false);
  35. declareProperty("MaxChi2", m_maxChi2 = 100.);
  36. }
  37.  
  38. //=============================================================================
  39. // Initialization
  40. //=============================================================================
  41. StatusCode TbAlignmentBase::initialize() {
  42.  
  43. // Initialise the base class.
  44. StatusCode sc = GaudiHistoTool::initialize();
  45. if (sc.isFailure()) return sc;
  46. m_modules = geomSvc()->modules();
  47. m_nPlanes = m_modules.size();
  48. m_masked.resize(m_nPlanes, false);
  49. for (const unsigned int plane : m_maskedPlanes) m_masked[plane] = true;
  50.  
  51. // Set the degrees of freedom.
  52. if (m_dofs.size() != 6) {
  53. info() << "Using the default degrees of freedom:" << endmsg;
  54. if (m_dofsDefault.size() != 6) {
  55. // Default DoFs are not defined. Set them.
  56. m_dofsDefault = {true, true, false, true, true, true};
  57. }
  58. m_dofs = m_dofsDefault;
  59. } else {
  60. info() << "Using the following degrees of freedom:" << endmsg;
  61. }
  62. // Print the degrees of freedom.
  63. const std::vector<std::string> labels = {"Translation X", "Translation Y",
  64. "Translation Z", "Rotation X",
  65. "Rotation Y", "Rotation Z"};
  66. for (unsigned int i = 0; i < 6; ++i) {
  67. const std::string on = m_dofs[i] ? "ON" : "OFF";
  68. info() << format(" %-14s %-3s", labels[i].c_str(), on.c_str()) << endmsg;
  69. }
  70.  
  71. // Set up the track fit tool.
  72. m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
  73. if (!m_trackFit) {
  74. return Error("Failed to initialise track fit.");
  75. }
  76. return StatusCode::SUCCESS;
  77. }
  78.  
  79. void TbAlignmentBase::plotResiduals(std::vector<TbAlignmentTrack*>& tracks,
  80. const std::string& tag) {
  81.  
  82. if (!m_monitoring) return;
  83. // Book histograms.
  84. std::vector<AIDA::IHistogram1D*> hResGlobalX;
  85. std::vector<AIDA::IHistogram1D*> hResGlobalY;
  86. std::vector<AIDA::IProfile1D*> hResGlobalXvsGlobalX;
  87. std::vector<AIDA::IProfile1D*> hResGlobalYvsGlobalY;
  88. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  89. const std::string plane = std::to_string(i);
  90. const std::string title = "Plane " + plane;
  91. // Distribution of unbiased global x-residuals.
  92. std::string name = tag + "/UnbiasedResiduals/GlobalX/Plane" + plane;
  93. double low(-0.2);
  94. double high(0.2);
  95. unsigned int bins(200);
  96. hResGlobalX.push_back(book1D(name, title, low, high, bins));
  97. // Distribution of unbiased global y-residuals.
  98. name = tag + "/UnbiasedResiduals/GlobalY/Plane" + plane;
  99. hResGlobalY.push_back(book1D(name, title, low, high, bins));
  100. // Profile of unbiased global x-residuals as function of global x.
  101. low = -20.;
  102. high = 20.;
  103. bins = 200;
  104. name = "UnbiasedResiduals/GlobalResXvsGlobalX/Plane" + plane;
  105. hResGlobalXvsGlobalX.push_back(
  106. bookProfile1D(name, title, low, high, bins));
  107. // Profile of unbiased global y-residuals as function of global y.
  108. name = "UnbiasedResiduals/GlobalResYvsGlobalY/Plane" + plane;
  109. hResGlobalYvsGlobalY.push_back(
  110. bookProfile1D(name, title, low, high, bins));
  111. }
  112. std::vector<double> ty(m_nPlanes, 0.);
  113. std::vector<double> yty(m_nPlanes, 0.);
  114. std::vector<double> y(m_nPlanes, 0.);
  115. std::vector<double> yy(m_nPlanes, 0.);
  116. std::vector<double> tyty(m_nPlanes, 0.);
  117.  
  118. for (auto& at : tracks) {
  119. // Get the track object.
  120. LHCb::TbTrack* track = at->track();
  121. // Loop over the clusters on the track.
  122. auto clusters = track->clusters();
  123. for (auto ic = clusters.cbegin(), end = clusters.cend(); ic != end; ++ic) {
  124. const unsigned int plane = (*ic)->plane();
  125. const LHCb::TbCluster* cluster = *ic;
  126. // Refit the track without this cluster.
  127. m_trackFit->maskPlane(plane);
  128. m_trackFit->fit(track);
  129. // Calculate the global x and y residuals.
  130. const Gaudi::XYZPoint intercept = geomSvc()->intercept(track, plane);
  131. const auto pLocal = Gaudi::XYZPoint(cluster->xloc(), cluster->yloc(), 0);
  132. const auto pGlobal = geomSvc()->localToGlobal(pLocal, plane);
  133. const double dxug = pGlobal.x() - intercept.x();
  134. const double dyug = pGlobal.y() - intercept.y();
  135.  
  136. ty[plane] += track->firstState().ty();
  137. yty[plane] += track->firstState().ty() * pGlobal.y();
  138. y[plane] += pGlobal.y();
  139. yy[plane] += pGlobal.y() * pGlobal.y();
  140. tyty[plane] += track->firstState().ty() * track->firstState().ty();
  141.  
  142. hResGlobalX[plane]->fill(dxug);
  143. hResGlobalY[plane]->fill(dyug);
  144. hResGlobalXvsGlobalX[plane]->fill(pGlobal.x(), dxug);
  145. hResGlobalYvsGlobalY[plane]->fill(pGlobal.y(), dyug);
  146. m_trackFit->unmaskPlane(plane);
  147. }
  148. m_trackFit->fit(track);
  149. // Loop over the associated clusters.
  150. auto aclusters = track->associatedClusters();
  151. for (auto it = aclusters.cbegin(), end = aclusters.cend(); it != end; ++it) {
  152. const unsigned int plane = (*it)->plane();
  153. // Calculate the global x and y residuals.
  154. const Gaudi::XYZPoint intercept = geomSvc()->intercept(track, plane);
  155. const auto pLocal = Gaudi::XYZPoint((*it)->xloc(), (*it)->yloc(), 0);
  156. const auto pGlobal = geomSvc()->localToGlobal(pLocal, plane);
  157. const double dxug = pGlobal.x() - intercept.x();
  158. const double dyug = pGlobal.y() - intercept.y();
  159. hResGlobalX[plane]->fill(dxug);
  160. hResGlobalY[plane]->fill(dyug);
  161. hResGlobalXvsGlobalX[plane]->fill(pGlobal.x(), dxug);
  162. hResGlobalYvsGlobalY[plane]->fill(pGlobal.y(), dyug);
  163. }
  164. }
  165.  
  166. if (msgLevel(MSG::DEBUG)) {
  167. const unsigned int nTracks = tracks.size();
  168. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  169. y[i] /= nTracks;
  170. ty[i] /= nTracks;
  171. yty[i] /= nTracks;
  172. yy[i] /= nTracks;
  173. tyty[i] /= nTracks;
  174. info() << "Plane " << i << ": Pearson-coefficient = "
  175. << (yty[i] - y[i] * ty[i]) /
  176. (sqrt(yy[i] - y[i] * y[i]) *
  177. sqrt(tyty[i] - ty[i] * ty[i])) << endmsg;
  178. }
  179. }
  180. // Fit the residual distributions and print the fit results.
  181. const std::string line(85, '-');
  182. info() << line << endmsg;
  183. info() << "Plane Mean X [\u03bcm] Mean Y [\u03bcm] "
  184. << "Sigma X [\u03bcm] Sigma Y [\u03bcm]" << endmsg;
  185. info() << line << endmsg;
  186. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  187. auto hx = Gaudi::Utils::Aida2ROOT::aida2root(hResGlobalX[i]);
  188. auto hy = Gaudi::Utils::Aida2ROOT::aida2root(hResGlobalY[i]);
  189. if (hx->GetEntries() == 0 || hy->GetEntries() == 0) continue;
  190. TFitResultPtr rx = hx->Fit("gaus", "QS0");
  191. TFitResultPtr ry = hy->Fit("gaus", "QS0");
  192. if (!rx.Get() || !ry.Get()) continue;
  193. info() << std::setprecision(4);
  194. info() << format(" %3u % 5.3f +/- %4.3f % 5.3f +/- %4.3f ", i,
  195. 1.e3 * rx->Parameter(1), 1.e3 * rx->Error(1),
  196. 1.e3 * ry->Parameter(1), 1.e3 * ry->Error(1))
  197. << format(" %4.3f +/- %4.3f %4.3f +/- %4.3f",
  198. 1.e3 * rx->Parameter(2), 1.e3 * rx->Error(2),
  199. 1.e3 * ry->Parameter(2), 1.e3 * ry->Error(2)) << endmsg;
  200. }
  201. // Fit the profiles (x/y residuals vs. global x/y) and print the results.
  202. info() << line << endmsg;
  203. info() << "Plane Slope X [\u03bcm] Slope Y [\u03bcm]" << endmsg;
  204. info() << line << endmsg;
  205. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  206. auto xgx = Gaudi::Utils::Aida2ROOT::aida2root(hResGlobalXvsGlobalX[i]);
  207. auto ygy = Gaudi::Utils::Aida2ROOT::aida2root(hResGlobalYvsGlobalY[i]);
  208. if (xgx->GetEntries() == 0 || ygy->GetEntries() == 0) continue;
  209. TFitResultPtr rbx = xgx->Fit("pol1", "QS0");
  210. TFitResultPtr rby = ygy->Fit("pol1", "QS0");
  211. if (!rbx.Get() || !rby.Get()) continue;
  212. info() << std::setprecision(4);
  213. info() << format(" %3u % 5.3f +/- %4.3f % 5.3f +/- %4.3f ", i,
  214. 1.e3 * rbx->Parameter(1), 1.e3 * rbx->Error(1),
  215. 1.e3 * rby->Parameter(1), 1.e3 * rby->Error(1))
  216. << endmsg;
  217. }
  218. info() << line << endmsg;
  219.  
  220. }
  221.  
  222. //=============================================================================
  223. // Collect alignment tracks (called at each event).
  224. //=============================================================================
  225. StatusCode TbAlignmentBase::execute(
  226. std::vector<TbAlignmentTrack*>& alignmentTracks) {
  227.  
  228. // Get the tracks.
  229. LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  230. if (!tracks) {
  231. return Error("No tracks in " + m_trackLocation);
  232. }
  233. // Add them to the alignment track store.
  234. for (LHCb::TbTrack* track : *tracks) {
  235. if (track->chi2() > m_maxChi2 || isEdge( track ) ) continue ;
  236. TbAlignmentTrack* alignmentTrack = new TbAlignmentTrack(track);
  237. alignmentTracks.push_back(alignmentTrack);
  238. }
  239. return StatusCode::SUCCESS;
  240. }