From 7ce5f9613e2d36e9b9eaefe7ed1317af21148124 Mon Sep 17 00:00:00 2001
From: Markus Prim <markus.prim@cern.ch>
Date: Fri, 5 Nov 2021 18:15:49 +0100
Subject: [PATCH] First prototyp for tracker station alignment. The prototype
operates in the global coordinate space and is thus only be able to align
xyz.
---
.../CMakeLists.txt | 18 +
.../GlobalChi2EquationGenerator/ReadMe | 4 +
.../src/GlobalChi2EquationGenerator.cxx | 558 ++++++++++++++++++
.../src/GlobalChi2EquationGenerator.h | 171 ++++++
.../GlobalChi2EquationGenerator_entries.cxx | 4 +
5 files changed, 755 insertions(+)
create mode 100644 Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/CMakeLists.txt
create mode 100644 Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/ReadMe
create mode 100755 Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.cxx
create mode 100755 Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.h
create mode 100644 Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/components/GlobalChi2EquationGenerator_entries.cxx
diff --git a/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/CMakeLists.txt b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/CMakeLists.txt
new file mode 100644
index 000000000..3f468a9df
--- /dev/null
+++ b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/CMakeLists.txt
@@ -0,0 +1,18 @@
+################################################################################
+# Package: GlobalChi2EquationGenerator
+################################################################################
+
+# Declare the package name:
+atlas_subdir( GlobalChi2EquationGenerator )
+
+# External dependencies:
+find_package( Eigen )
+
+find_package( ROOT COMPONENTS Core Tree MathCore)
+
+atlas_add_component( GlobalChi2EquationGenerator
+ src/components/*.cxx src/*.cxx src/*.h
+ INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${EIGEN_LIBRARIES} ${ROOT_LIBRARIES} AthenaBaseComps AthContainers GeoPrimitives Identifier GaudiKernel TrackerReadoutGeometry TrackerPrepRawData TrackerSpacePoint AthenaMonitoringKernelLib FaserDetDescr xAODEventInfo TrkEventPrimitives TrackerIdentifier FaserSiSpacePointToolLib TrkSpacePoint )
+
+atlas_install_python_modules( python/*.py )
diff --git a/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/ReadMe b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/ReadMe
new file mode 100644
index 000000000..9caf9456f
--- /dev/null
+++ b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/ReadMe
@@ -0,0 +1,4 @@
+GlobalChi2EquationGenerator:
+ Generates the matrix and vector for the global chi2 alignment equation.
+
+For each station, select the spacepoint triplets first, then perform a 3D line fit and require the chi2 less than chi2max.
diff --git a/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.cxx b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.cxx
new file mode 100755
index 000000000..8a59cca70
--- /dev/null
+++ b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.cxx
@@ -0,0 +1,558 @@
+/*
+ Copyright (C) 2002-2021 ChRN for the benefit of the ATLAS collaboration
+ */
+
+/***************************************************************************
+ -------------------
+ ***************************************************************************/
+
+//<<<<<< INCLUDES >>>>>>
+
+
+#include "TrackerSPFit.h"
+
+// For processing clusters
+#include "TrackerReadoutGeometry/SiLocalPosition.h"
+#include "TrackerReadoutGeometry/SiDetectorDesign.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+
+// Space point Classes,
+#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+
+
+// general Atlas classes
+#include "FaserDetDescr/FaserDetectorID.h"
+#include "xAODEventInfo/EventInfo.h"
+#include "StoreGate/ReadCondHandle.h"
+
+#include "AthenaMonitoringKernel/Monitored.h"
+
+#include <fstream>
+
+namespace Tracker
+{
+ //------------------------------------------------------------------------
+ TrackerSPFit::TrackerSPFit(const std::string& name,
+ ISvcLocator* pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+ , m_thistSvc("THistSvc", name)
+ {
+ declareProperty("SaveAllSPInfo", m_saveallsp=false);
+ declareProperty("chi2Max", m_chi2max=10);
+ declareProperty("OutputDir", m_outputdir="./");
+
+ m_t.resize(144, 1);
+ m_t.setZero();
+ m_M.resize(144, 144);
+ m_M.setZero();
+
+ for (unsigned int i=0; i < 8; ++i) {
+ TrackerSPFit::AlignVector18D t;
+ t.setZero();
+ m_ts.push_back(t);
+ TrackerSPFit::AlignMatrix18x18 M;
+ M.setZero();
+ m_Ms.push_back(M);
+ }
+ }
+
+ //-----------------------------------------------------------------------
+ StatusCode TrackerSPFit::initialize()
+ {
+ ATH_MSG_DEBUG( "TrackerSPFit::initialize()" );
+
+ CHECK(m_thistSvc.retrieve());
+ // Check that clusters, space points and ids have names
+ if(m_Sct_clcontainerKey.key().empty()){
+ ATH_MSG_FATAL( "SCTs selected and no name set for SCT clusters");
+ return StatusCode::FAILURE;
+ }
+ ATH_CHECK( m_Sct_clcontainerKey.initialize() );
+
+ if ( m_Sct_spcontainerKey.key().empty()){
+ ATH_MSG_FATAL( "SCTs selected and no name set for SCT clusters");
+ return StatusCode::FAILURE;
+ }
+ ATH_CHECK( m_Sct_spcontainerKey.initialize() );
+
+ m_tree= new TTree("spfit","tree");
+ TrackerSPFit::initializeTree();
+
+ // create containers (requires the Identifier Helpers)
+ ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
+
+ // Initialize the key of input SiElementPropertiesTable and SiDetectorElementCollection for SCT
+ ATH_CHECK(m_SCTDetEleCollKey.initialize());
+
+
+ CHECK(m_thistSvc->regTree("/TrackerSPFit/spfit",m_tree));
+ ATH_MSG_INFO( "TrackerSPFit::initialized" );
+ return StatusCode::SUCCESS;
+ }
+
+ //-------------------------------------------------------------------------
+ void TrackerSPFit::initializeTree(){
+ m_tree->Branch("evtId",&m_eventNumber);
+ m_tree->Branch("initial_track_chi2", &m_initial_track_chi2);
+ m_tree->Branch("initial_track_cx", &m_initial_track_cx);
+ m_tree->Branch("initial_track_mx", &m_initial_track_mx);
+ m_tree->Branch("initial_track_cy", &m_initial_track_cy);
+ m_tree->Branch("initial_track_my", &m_initial_track_my);
+ m_tree->Branch("initial_track_residual_x0", &m_initial_track_residual_x0);
+ m_tree->Branch("initial_track_residual_x1", &m_initial_track_residual_x1);
+ m_tree->Branch("initial_track_residual_x2", &m_initial_track_residual_x2);
+ m_tree->Branch("initial_track_residual_y0", &m_initial_track_residual_y0);
+ m_tree->Branch("initial_track_residual_y1", &m_initial_track_residual_y1);
+ m_tree->Branch("initial_track_residual_y2", &m_initial_track_residual_y2);
+ m_tree->Branch("x0", &m_x0);
+ m_tree->Branch("x1", &m_x1);
+ m_tree->Branch("x2", &m_x2);
+ m_tree->Branch("y0", &m_y0);
+ m_tree->Branch("y1", &m_y1);
+ m_tree->Branch("y2", &m_y2);
+ m_tree->Branch("sigma_x0", &m_sigma_x0);
+ m_tree->Branch("sigma_x1", &m_sigma_x1);
+ m_tree->Branch("sigma_x2", &m_sigma_x2);
+ m_tree->Branch("sigma_y0", &m_sigma_y0);
+ m_tree->Branch("sigma_y1", &m_sigma_y1);
+ m_tree->Branch("sigma_y2", &m_sigma_y2);
+ m_tree->Branch("z0", &m_z0);
+ m_tree->Branch("z1", &m_z1);
+ m_tree->Branch("z2", &m_z2);
+ m_tree->Branch("module0", &m_module0);
+ m_tree->Branch("module1", &m_module1);
+ m_tree->Branch("module2", &m_module2);
+ m_tree->Branch("x0_local", &m_x0_local);
+ m_tree->Branch("x1_local", &m_x1_local);
+ m_tree->Branch("x2_local", &m_x2_local);
+ m_tree->Branch("y0_local", &m_y0_local);
+ m_tree->Branch("y1_local", &m_y1_local);
+ m_tree->Branch("y2_local", &m_y2_local);
+ // m_tree->Branch("z0_local", &m_z0_local);
+ // m_tree->Branch("z1_local", &m_z1_local);
+ // m_tree->Branch("z2_local", &m_z2_local);
+ if(m_saveallsp){
+ m_tree->Branch("sp_all_x", &m_sp_all_x);
+ m_tree->Branch("sp_all_y", &m_sp_all_y);
+ m_tree->Branch("sp_all_z", &m_sp_all_z);
+ m_tree->Branch("sp_all_station", &m_sp_all_station);
+ m_tree->Branch("sp_all_layer", &m_sp_all_layer);
+ m_tree->Branch("sp_all_module", &m_sp_all_module);
+ }
+ }
+ //-------------------------------------------------------------------------
+ void TrackerSPFit::clearVariables() const {
+ m_initial_track_chi2.clear();
+ m_initial_track_cx.clear();
+ m_initial_track_mx.clear();
+ m_initial_track_cy.clear();
+ m_initial_track_my.clear();
+ m_initial_track_residual_x0.clear();
+ m_initial_track_residual_x1.clear();
+ m_initial_track_residual_x2.clear();
+ m_initial_track_residual_y0.clear();
+ m_initial_track_residual_y1.clear();
+ m_initial_track_residual_y2.clear();
+ m_x0.clear();
+ m_x1.clear();
+ m_x2.clear();
+ m_y0.clear();
+ m_y1.clear();
+ m_y2.clear();
+ m_sigma_x0.clear();
+ m_sigma_x1.clear();
+ m_sigma_x2.clear();
+ m_sigma_y0.clear();
+ m_sigma_y1.clear();
+ m_sigma_y2.clear();
+ m_z0.clear();
+ m_z1.clear();
+ m_z2.clear();
+ m_module0.clear();
+ m_module1.clear();
+ m_module2.clear();
+ m_x0_local.clear();
+ m_x1_local.clear();
+ m_x2_local.clear();
+ m_y0_local.clear();
+ m_y1_local.clear();
+ m_y2_local.clear();
+ // m_z0_local.clear();
+ // m_z1_local.clear();
+ // m_z2_local.clear();
+
+ if(m_saveallsp){
+ m_sp_all_x.clear();
+ m_sp_all_y.clear();
+ m_sp_all_z.clear();
+ m_sp_all_station.clear();
+ m_sp_all_layer.clear();
+ m_sp_all_module.clear();
+ }
+
+
+ }
+ //-------------------------------------------------------------------------
+
+ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
+ {
+
+ ++m_numberOfEvents;
+ m_eventNumber=m_numberOfEvents;
+ clearVariables();
+ const TrackerDD::SiDetectorElementCollection* elements = nullptr;
+
+
+ SG::ReadCondHandle<TrackerDD::SiDetectorElementCollection> sctDetEle(m_SCTDetEleCollKey, ctx);
+ elements = sctDetEle.retrieve();
+ if (elements==nullptr) {
+ ATH_MSG_FATAL("Pointer of SiDetectorElementCollection (" << m_SCTDetEleCollKey.fullKey() << ") could not be retrieved");
+ return StatusCode::SUCCESS;
+ }
+
+ // retrieve SCT spacepoint container
+ SG::ReadHandle<FaserSCT_SpacePointContainer> sct_spcontainer( m_Sct_spcontainerKey, ctx );
+ if (!sct_spcontainer.isValid()){
+ msg(MSG:: FATAL) << "Could not find the data object "<< sct_spcontainer.name() << " !" << endmsg;
+ return StatusCode::RECOVERABLE;
+ }
+
+ // save all the sp information
+ if(m_saveallsp){
+ ATH_MSG_DEBUG( "SCT spacepoint container found: " << sct_spcontainer->size() << " collections" );
+
+ for (const auto& colNext : *sct_spcontainer) {
+ ++m_numberOfSPCol;
+ Identifier elementID = colNext->identify();
+ for (const auto& sp : *colNext) {
+ ++m_numberOfSP;
+ m_sp_all_x.push_back(sp->globalPosition().x());
+ m_sp_all_y.push_back(sp->globalPosition().y());
+ m_sp_all_z.push_back(sp->globalPosition().z());
+ m_sp_all_station.push_back(m_idHelper->station(elementID));
+ m_sp_all_layer.push_back(m_idHelper->layer(elementID));
+ auto ietamodule = m_idHelper->eta_module(elementID);
+ auto iphimodule = m_idHelper->phi_module(elementID);
+ if(ietamodule < 0) {
+ iphimodule += 4;
+ }
+ m_sp_all_module.push_back(iphimodule);
+ }
+ }
+ }
+
+ // register the IdentifiableContainer into StoreGate
+
+ ATH_MSG_DEBUG( "SCT spacepoint container found: " << sct_spcontainer->size() << " collections" );
+
+ // Collect all space points in a flat vector
+
+ std::vector<SP_Seed> space_points_layer0;
+ std::vector<SP_Seed> space_points_layer1;
+ std::vector<SP_Seed> space_points_layer2;
+ for (const auto& colNext : *sct_spcontainer) {
+ Identifier elementID = colNext->identify();
+ ATH_MSG_DEBUG("reading SpacePoints collection "<<elementID.get_compact()<<" , and hash = "<<colNext->identifyHash() );
+
+ auto istation = m_idHelper->station(elementID);
+ auto ilayer = m_idHelper->layer(elementID);
+ auto ietamodule = m_idHelper->eta_module(elementID);
+ auto iphimodule = m_idHelper->phi_module(elementID);
+ auto imodule = iphimodule;
+ if(ietamodule < 0) {
+ imodule += 4;
+ }
+ auto size = colNext->size();
+ if (size == 0) {
+ ATH_MSG_DEBUG( "TrackerSPFit algorithm found no space points" );
+ } else {
+ //In a MT environment the cache maybe filled by another thread in which case this will delete the duplicate
+ for (const auto& sp : *colNext) {
+ // sp->globalPosition();
+ // sp->globCovariance(),
+ struct SP_Seed tmp{
+ *sp,
+ ilayer,
+ imodule
+ };
+ if(istation == 0) { // We want to align station by station, so make this an argument and call it repeatedly.
+ if(ilayer==0) space_points_layer0.push_back(tmp);
+ if(ilayer==1) space_points_layer1.push_back(tmp);
+ if(ilayer==2) space_points_layer2.push_back(tmp);
+ }
+ }
+ ATH_MSG_DEBUG( size << " SpacePoints successfully added to Container !" );
+ }
+ }
+ ATH_MSG_DEBUG( "TrackerSPFit number of spacepoints in each layer: " << space_points_layer0.size() << " " << space_points_layer1.size() << " " << space_points_layer2.size() );
+
+
+ // First check to filter out freak events in the cosmic / test beam data
+ if ((space_points_layer0.size() > 8) or (space_points_layer1.size() > 8) or (space_points_layer0.size() > 8)) {
+ return StatusCode::SUCCESS;
+ }
+
+ // Arange all possible combinations of hit triplets
+ std::vector<HitTriplet> hit_triplets;
+ for (const auto& sp0 : space_points_layer0) {
+ for (const auto& sp1 : space_points_layer1) {
+ for (const auto& sp2 : space_points_layer2) {
+ hit_triplets.emplace_back(HitTriplet({sp0.sp, sp1.sp, sp2.sp, sp0.module, sp1.module, sp2.module}));
+ }
+ }
+ }
+ ATH_MSG_DEBUG( "TrackerSPFit number of hit triplets: " << hit_triplets.size() );
+
+ // Second check to filter out freak events in the cosmic / test beam data.
+ if (hit_triplets.size() > 8) {
+ return StatusCode::SUCCESS;
+ }
+
+ for (const auto& hit_triplet : hit_triplets) {
+ const AlignVector3D x{
+ hit_triplet.sp0.globalPosition()[0],
+ hit_triplet.sp1.globalPosition()[0],
+ hit_triplet.sp2.globalPosition()[0]
+ };
+ const AlignVector3D y{
+ hit_triplet.sp0.globalPosition()[1],
+ hit_triplet.sp1.globalPosition()[1],
+ hit_triplet.sp2.globalPosition()[1]
+ };
+ const AlignVector3D z{
+ hit_triplet.sp0.globalPosition()[2],
+ hit_triplet.sp1.globalPosition()[2],
+ hit_triplet.sp2.globalPosition()[2]
+ };
+ const auto sigma_x0 = hit_triplet.sp0.globCovariance()(0, 0);
+ const auto sigma_x1 = hit_triplet.sp1.globCovariance()(0, 0);
+ const auto sigma_x2 = hit_triplet.sp2.globCovariance()(0, 0);
+ const auto sigma_y0 = hit_triplet.sp0.globCovariance()(1, 1);
+ const auto sigma_y1 = hit_triplet.sp1.globCovariance()(1, 1);
+ const auto sigma_y2 = hit_triplet.sp2.globCovariance()(1, 1);
+
+ const auto W_ = W(
+ 1./sigma_x0,
+ 1./sigma_x1,
+ 1./sigma_x2,
+ 1./sigma_y0,
+ 1./sigma_y1,
+ 1./sigma_y2
+ );
+
+ const auto track_fit_result = SolveLocalTrackFitWeighted(A(z), b(x,y), W_);
+ const AlignVector3D x_local{
+ hit_triplet.sp0.localParameters().x(),
+ hit_triplet.sp1.localParameters().x(),
+ hit_triplet.sp2.localParameters().x()
+ };
+ const AlignVector3D y_local{
+ hit_triplet.sp0.localParameters().y(),
+ hit_triplet.sp1.localParameters().y(),
+ hit_triplet.sp2.localParameters().y()
+ };
+
+ AlignVector6D xy;
+ xy << x[0], x[1], x[2], y[0], y[1], y[2];
+ // std::cout << "xy: " << xy.transpose() << std::endl;
+ // std::cout << "z: " << z.transpose() << std::endl;
+
+ const auto prediction_layer0 = TrackPrediction(z[0], track_fit_result);
+ const auto prediction_layer1 = TrackPrediction(z[1], track_fit_result);
+ const auto prediction_layer2 = TrackPrediction(z[2], track_fit_result);
+ AlignVector6D xy_predicted;
+ xy_predicted << prediction_layer0[0], prediction_layer1[0], prediction_layer2[0], prediction_layer0[1], prediction_layer1[1], prediction_layer2[1];
+ // std::cout << "pr: " << xy_predicted.transpose() << std::endl;
+
+ const auto residual = xy_predicted - xy;
+ // std::cout << "r: " << residual.transpose() << std::endl;
+
+ const double chi2 = residual.transpose() * W_ * residual;
+ const auto& module0 = hit_triplet.module0;
+ const auto& module1 = hit_triplet.module1;
+ const auto& module2 = hit_triplet.module2;
+
+
+ m_initial_track_chi2.push_back(chi2);
+ m_initial_track_cx.push_back(track_fit_result[0]);
+ m_initial_track_mx.push_back(track_fit_result[1]);
+ m_initial_track_cy.push_back(track_fit_result[2]);
+ m_initial_track_my.push_back(track_fit_result[3]);
+ m_initial_track_residual_x0.push_back(residual[0]);
+ m_initial_track_residual_x1.push_back(residual[1]);
+ m_initial_track_residual_x2.push_back(residual[2]);
+ m_initial_track_residual_y0.push_back(residual[3]);
+ m_initial_track_residual_y1.push_back(residual[4]);
+ m_initial_track_residual_y2.push_back(residual[5]);
+ m_x0.push_back(x[0]);
+ m_x1.push_back(x[1]);
+ m_x2.push_back(x[2]);
+ m_y0.push_back(y[0]);
+ m_y1.push_back(y[1]);
+ m_y2.push_back(y[2]);
+ m_z0.push_back(z[0]);
+ m_z1.push_back(z[1]);
+ m_z2.push_back(z[2]);
+ m_sigma_x0.push_back(sigma_x0);
+ m_sigma_x1.push_back(sigma_x1);
+ m_sigma_x2.push_back(sigma_x2);
+ m_sigma_y0.push_back(sigma_y0);
+ m_sigma_y1.push_back(sigma_y1);
+ m_sigma_y2.push_back(sigma_y2);
+ m_module0.push_back(module0);
+ m_module1.push_back(module1);
+ m_module2.push_back(module2);
+ m_x0_local.push_back(x_local[0]);
+ m_x1_local.push_back(x_local[1]);
+ m_x2_local.push_back(x_local[2]);
+ m_y0_local.push_back(y_local[0]);
+ m_y1_local.push_back(y_local[1]);
+ m_y2_local.push_back(y_local[2]);
+ // m_z0_local.push_back(z[0]);
+ // m_z1_local.push_back(z[1]);
+ // m_z2_local.push_back(z[2]);
+ const double chi2max = 10;
+ if (chi2 < chi2max) {
+ m_numberOfAlignmentTracks += 1;
+ ATH_MSG_INFO( "TrackerSPFit first track fit result: chi2: " << chi2 << " pars: " << track_fit_result[0] << " " << track_fit_result[1] << " " << track_fit_result[2] << " " << track_fit_result[3]);
+ const double mx = track_fit_result[1];
+ const double my = track_fit_result[3];
+ const auto dr_da = AlignmentDerivatives(x, y, mx, my);
+ const auto t = dr_da.transpose() * What(W_, A(z)) * residual;
+ const auto M = dr_da.transpose() * What(W_, A(z)) * dr_da;
+
+ // Save the local alignment matrices (module wise)
+ if ((module0 == module1) and (module0 == module2)) {
+ m_ts[module0] += t;
+ m_Ms[module0] += M;
+ }
+ // std::cout << dr_da << std::endl;
+
+ // In the following assign the blocks to the global matrix
+ const auto index0 = GetIndexInSolutionMatrix(0, module0);
+ const auto index1 = GetIndexInSolutionMatrix(1, module1);
+ const auto index2 = GetIndexInSolutionMatrix(2, module2);
+
+ // Add t terms in correct position of global alignment matrix
+ m_t.block<6, 1>(index0, 0) += t.block<6, 1>( 0, 0);
+ m_t.block<6, 1>(index1, 0) += t.block<6, 1>( 6, 0);
+ m_t.block<6, 1>(index2, 0) += t.block<6, 1>(12, 0);
+
+ // Add t terms in correct position of global alignment matrix
+ m_M.block<6, 6>(index0, index0) += M.block<6, 6>( 0, 0);
+ m_M.block<6, 6>(index0, index1) += M.block<6, 6>( 0, 6);
+ m_M.block<6, 6>(index0, index2) += M.block<6, 6>( 0, 12);
+
+ m_M.block<6, 6>(index1, index0) += M.block<6, 6>( 6, 0);
+ m_M.block<6, 6>(index1, index1) += M.block<6, 6>( 6, 6);
+ m_M.block<6, 6>(index1, index2) += M.block<6, 6>( 6, 12);
+
+ m_M.block<6, 6>(index2, index0) += M.block<6, 6>(12, 0);
+ m_M.block<6, 6>(index2, index1) += M.block<6, 6>(12, 6);
+ m_M.block<6, 6>(index2, index2) += M.block<6, 6>(12, 12);
+
+ }
+ }
+
+ if (hit_triplets.size() > 0) m_tree->Fill();
+ return StatusCode::SUCCESS;
+ }
+
+ //---------------------------------------------------------------------------
+ StatusCode TrackerSPFit::finalize()
+ {
+ ATH_MSG_INFO( "TrackerSPFit::finalize()" );
+ ATH_MSG_INFO( m_numberOfEvents << " events processed" );
+ ATH_MSG_INFO( m_numberOfSPCol<< " sct SP collections processed" );
+ ATH_MSG_INFO( m_numberOfSP<< " sct SP processed" );
+ ATH_MSG_INFO( m_numberOfAlignmentTracks << " tracks are used for alignment." );
+
+ const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", "\n");
+ std::ofstream file_M("./" + m_outputdir + "/M.csv");
+ file_M << m_M.format(CSVFormat);
+ file_M << std::endl;
+ std::ofstream file_t("./" + m_outputdir + "/t.csv");
+ file_t << m_t.format(CSVFormat);
+ file_t << std::endl;
+
+ for (unsigned int i=0; i<8; ++i) {
+ std::ofstream file_M("./" + m_outputdir + "/M" + std::to_string(i) + ".csv");
+ file_M << m_Ms[i].format(CSVFormat);
+ file_M << std::endl;
+ std::ofstream file_t("./" + m_outputdir + "/t" + std::to_string(i) + ".csv");
+ file_t << m_ts[i].format(CSVFormat);
+ file_t << std::endl;
+ }
+
+ return StatusCode::SUCCESS;
+ }
+
+ //--------------------------------------------------------------------------
+
+ TrackerSPFit::AlignVector3D TrackerSPFit::TrackPrediction(double z, const TrackerSPFit::AlignVector4D& parameters) const {
+ AlignVector3D prediction;
+ prediction << parameters[0] + parameters[1] * z, parameters[2] + parameters[3] * z, z;
+ return prediction;
+ }
+
+ TrackerSPFit::AlignMatrix6x4 TrackerSPFit::A(const TrackerSPFit::AlignVector3D& z) const {
+ TrackerSPFit::AlignMatrix6x4 A;
+ A << 1, z[0], 0, 0,
+ 1, z[1], 0, 0,
+ 1, z[2], 0, 0,
+ 0, 0, 1, z[0],
+ 0, 0, 1, z[1],
+ 0, 0, 1, z[2];
+ return A;
+ }
+
+ TrackerSPFit::AlignVector6D TrackerSPFit::b(const TrackerSPFit::AlignVector3D& x, const TrackerSPFit::AlignVector3D& y) const {
+ TrackerSPFit::AlignVector6D b;
+ b << x[0], x[1], x[2], y[0], y[1], y[2];
+ return b;
+ }
+
+ TrackerSPFit::AlignVector4D TrackerSPFit::SolveLocalTrackFitWeighted(const TrackerSPFit::AlignMatrix6x4& A, const TrackerSPFit::AlignVector6D& b, const TrackerSPFit::AlignDiagonalMatrix6x6& W) const {
+ return (A.transpose() * W * A).ldlt().solve(A.transpose() * W * b);
+ }
+
+ TrackerSPFit::AlignMatrix4x4 TrackerSPFit::TrackParameterCovariance(const TrackerSPFit::AlignMatrix6x4& A, const TrackerSPFit::AlignDiagonalMatrix6x6& W) const {
+ return (A.transpose() * A).inverse() * A.transpose() * W.inverse() * A * (A.transpose() * A).inverse();
+ }
+
+ TrackerSPFit::AlignDiagonalMatrix6x6 TrackerSPFit::W(const double w_x0, const double w_x1, const double w_x2, const double w_y0, const double w_y1, const double w_y2) const {
+ TrackerSPFit::AlignDiagonalMatrix6x6 W;
+ W.diagonal() << w_x0, w_x1, w_x2, w_y0, w_y1, w_y2;
+ return W;
+ }
+
+ TrackerSPFit::AlignMatrix6x6 TrackerSPFit::What(const TrackerSPFit::AlignDiagonalMatrix6x6& W, const TrackerSPFit::AlignMatrix6x4& A) const {
+ auto I = TrackerSPFit::AlignMatrix6x6::Identity();
+ return (I - A * (A.transpose() * W * A).inverse() * A.transpose() * W).transpose() * W;
+ }
+
+ TrackerSPFit::AlignMatrix6x18 TrackerSPFit::AlignmentDerivatives(const TrackerSPFit::AlignVector3D& x, const TrackerSPFit::AlignVector3D& y, double mx, double my) const {
+ TrackerSPFit::AlignVector6D dx_da;
+ dx_da << -1, 0, mx, 0, 0, 0;
+ TrackerSPFit::AlignVector6D dy_da;
+ dy_da << 0, -1, my, 0, 0, 0;
+
+ TrackerSPFit::AlignMatrix6x18 dr_da; // FIXME: Last three entries are currently 0, these are reserved for rotations_xy, that is when we need x, y
+ dr_da.setZero();
+ dr_da.block<1, 6>( 0, 0) += dx_da.transpose(); // Layer 0
+ dr_da.block<1, 6>( 1, 6) += dx_da.transpose(); // Layer 1
+ dr_da.block<1, 6>( 2, 12) += dx_da.transpose(); // Layer 2
+ dr_da.block<1, 6>( 3, 0) += dy_da.transpose(); // Layer 0
+ dr_da.block<1, 6>( 4, 6) += dy_da.transpose(); // Layer 1
+ dr_da.block<1, 6>( 5, 12) += dy_da.transpose(); // Layer 2
+ return dr_da;
+ }
+
+ unsigned int TrackerSPFit::GetIndexInSolutionMatrix(const unsigned int layer, const unsigned int module) const {
+ const unsigned int parameters_per_layer = 6;
+ const unsigned int layers_per_module = 3;
+ return layer * parameters_per_layer + module * parameters_per_layer * layers_per_module;
+ }
+
+}
diff --git a/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.h b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.h
new file mode 100755
index 000000000..95bd2f5fc
--- /dev/null
+++ b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/GlobalChi2EquationGenerator.h
@@ -0,0 +1,171 @@
+// -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef TRACKERSPFIT_TRACKERSPACEPOINTMAKERALG_H
+#define TRACKERSPFIT_TRACKERSPACEPOINTMAKERALG_H
+
+#include "StoreGate/ReadCondHandleKey.h"
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "Identifier/Identifier.h"
+
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "TrackerPrepRawData/TrackerClusterContainer.h"
+#include "TrackerReadoutGeometry/SiDetectorElementCollection.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePoint.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerSpacePoint/FaserSCT_SpacePointOverlapCollection.h"
+#include "TrkEventPrimitives/LocalParameters.h"
+
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ITHistSvc.h"
+
+#include <string>
+#include <vector>
+#include "TTree.h"
+
+#include <Eigen/Sparse>
+
+#include <string>
+
+class FaserSCT_ID;
+namespace Tracker
+{
+
+ class TrackerSPFit:public AthReentrantAlgorithm {
+
+ // Define typedefs which shortcut eigen matrices and vectors for the dimensions we need in our use case.
+ typedef Eigen::Matrix<double, 4, 4> AlignMatrix4x4;
+ typedef Eigen::Matrix<double, 6, 4> AlignMatrix6x4;
+ typedef Eigen::Matrix<double, 6, 6> AlignMatrix6x6;
+ typedef Eigen::Matrix<double, 6, 18> AlignMatrix6x18;
+ typedef Eigen::Matrix<double, 18, 18> AlignMatrix18x18;
+ typedef Eigen::Matrix<double, 144, 144> AlignMatrix144x144;
+ typedef Eigen::DiagonalMatrix<double, 6> AlignDiagonalMatrix6x6;
+ typedef Eigen::Matrix<double, 3, 1> AlignVector3D;
+ typedef Eigen::Matrix<double, 4, 1> AlignVector4D;
+ typedef Eigen::Matrix<double, 6, 1> AlignVector6D;
+ typedef Eigen::Matrix<double, 18, 1> AlignVector18D;
+ typedef Eigen::Matrix<double, 144, 1> AlignVector144D;
+
+ public:
+
+ /**
+ * @name AthReentrantAlgorithm methods
+ */
+ //@{
+ TrackerSPFit(const std::string& name,
+ ISvcLocator* pSvcLocator);
+
+ virtual ~TrackerSPFit() = default;
+
+ virtual StatusCode initialize() override;
+
+ virtual StatusCode execute (const EventContext& ctx) const override;
+
+ virtual StatusCode finalize() override;
+
+ struct SP_Seed{FaserSCT_SpacePoint sp; int layer; int module;};
+ struct HitTriplet{
+ FaserSCT_SpacePoint sp0;
+ FaserSCT_SpacePoint sp1;
+ FaserSCT_SpacePoint sp2;
+ int module0;
+ int module1;
+ int module2;
+ };
+
+ private:
+ AlignVector3D TrackPrediction(double z, const AlignVector4D& parameters) const;
+ AlignMatrix6x4 A(const AlignVector3D& z) const;
+ AlignVector6D b(const AlignVector3D& x, const AlignVector3D& y) const;
+ AlignDiagonalMatrix6x6 W(const double w_x0, const double w_x1, const double w_x2, const double w_y0, const double w_y1, const double w_y2) const;
+ AlignMatrix6x6 What(const AlignDiagonalMatrix6x6& W, const AlignMatrix6x4& A) const;
+ AlignVector4D SolveLocalTrackFitWeighted(const AlignMatrix6x4& A, const AlignVector6D& b, const AlignDiagonalMatrix6x6& W) const;
+ AlignMatrix4x4 TrackParameterCovariance(const AlignMatrix6x4& A, const AlignDiagonalMatrix6x6& W) const;
+ AlignMatrix6x18 AlignmentDerivatives(const AlignVector3D& x, const AlignVector3D& y, double mx, double my) const;
+ unsigned int GetIndexInSolutionMatrix(const unsigned int layer, const unsigned int module) const;
+
+ void initializeTree();
+ void clearVariables() const;
+ mutable TTree* m_tree;
+ TrackerSPFit() = delete;
+ TrackerSPFit(const TrackerSPFit&) =delete;
+ TrackerSPFit &operator=(const TrackerSPFit&) = delete;
+
+ SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_Sct_spcontainerKey{this, "SpacePointsSCTName", "SCT_SpacePointContainer"};
+ SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_Sct_clcontainerKey{this, "SCT_ClustersName", "SCT clContainer"};
+ SG::ReadCondHandleKey<TrackerDD::SiDetectorElementCollection> m_SCTDetEleCollKey{this, "SCTDetEleCollKey", "SCT_DetectorElementCollection", "Key of SiDetectorElementCollection for SCT"};
+
+ const FaserSCT_ID* m_idHelper{nullptr};
+ mutable std::atomic<int> m_numberOfEvents{0};
+ mutable std::atomic<int> m_numberOfSCT{0};
+ mutable std::atomic<int> m_numberOfTrack{0};
+ mutable std::atomic<int> m_numberOfSP{0};
+ mutable std::atomic<int> m_numberOfSPCol{0};
+ mutable std::atomic<int> m_numberOfAlignmentTracks{0};
+ //@}
+ ServiceHandle<ITHistSvc> m_thistSvc;
+
+ mutable int m_eventNumber;
+ mutable std::vector<double> m_initial_track_chi2;
+ mutable std::vector<double> m_initial_track_cx;
+ mutable std::vector<double> m_initial_track_mx;
+ mutable std::vector<double> m_initial_track_cy;
+ mutable std::vector<double> m_initial_track_my;
+ mutable std::vector<double> m_initial_track_residual_x0;
+ mutable std::vector<double> m_initial_track_residual_x1;
+ mutable std::vector<double> m_initial_track_residual_x2;
+ mutable std::vector<double> m_initial_track_residual_y0;
+ mutable std::vector<double> m_initial_track_residual_y1;
+ mutable std::vector<double> m_initial_track_residual_y2;
+ mutable std::vector<double> m_x0;
+ mutable std::vector<double> m_x1;
+ mutable std::vector<double> m_x2;
+ mutable std::vector<double> m_y0;
+ mutable std::vector<double> m_y1;
+ mutable std::vector<double> m_y2;
+ mutable std::vector<double> m_sigma_x0;
+ mutable std::vector<double> m_sigma_x1;
+ mutable std::vector<double> m_sigma_x2;
+ mutable std::vector<double> m_sigma_y0;
+ mutable std::vector<double> m_sigma_y1;
+ mutable std::vector<double> m_sigma_y2;
+ mutable std::vector<double> m_z0;
+ mutable std::vector<double> m_z1;
+ mutable std::vector<double> m_z2;
+ mutable std::vector<double> m_module0;
+ mutable std::vector<double> m_module1;
+ mutable std::vector<double> m_module2;
+ mutable std::vector<double> m_x0_local;
+ mutable std::vector<double> m_x1_local;
+ mutable std::vector<double> m_x2_local;
+ mutable std::vector<double> m_y0_local;
+ mutable std::vector<double> m_y1_local;
+ mutable std::vector<double> m_y2_local;
+ // mutable std::vector<double> m_z0_local;
+ // mutable std::vector<double> m_z1_local;
+ // mutable std::vector<double> m_z2_local;
+ bool m_saveallsp;
+ bool m_chi2max;
+ std::string m_outputdir;
+ mutable std::vector<double> m_sp_all_x;
+ mutable std::vector<double> m_sp_all_y;
+ mutable std::vector<double> m_sp_all_z;
+ mutable std::vector<int> m_sp_all_station;
+ mutable std::vector<int> m_sp_all_layer;
+ mutable std::vector<int> m_sp_all_module;
+
+ mutable Eigen::MatrixXd m_t;
+ mutable Eigen::MatrixXd m_M;
+
+ mutable std::vector<AlignVector18D> m_ts;
+ mutable std::vector<AlignMatrix18x18> m_Ms;
+
+ };
+
+}
+#endif // TrackerSpacePointMakerAlg_TRACKERSPACEPOINTMAKERALG_H
diff --git a/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/components/GlobalChi2EquationGenerator_entries.cxx b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/components/GlobalChi2EquationGenerator_entries.cxx
new file mode 100644
index 000000000..dae3ad64e
--- /dev/null
+++ b/Tracker/TrackerAlignAlgs/GlobalChi2EquationGenerator/src/components/GlobalChi2EquationGenerator_entries.cxx
@@ -0,0 +1,4 @@
+#include "../TrackerSPFit.h"
+
+DECLARE_COMPONENT( Tracker::TrackerSPFit )
+
--
GitLab