packages/trackreco/WeightedFitter.cc

0001 #include "WeightedFitter.h"
0002
0003 #include <trackbase/ActsGeometry.h>
0004 #include <trackbase/ActsSurfaceMaps.h>
0005 #include <trackbase/ClusterErrorPara.h>
0006 #include <trackbase/InttDefs.h>
0007 #include <trackbase/MvtxDefs.h>
0008 #include <trackbase/TpcDefs.h>
0009 #include <trackbase/TrkrCluster.h>
0010 #include <trackbase/TrkrClusterContainerv4.h>
0011
0012 #include <micromegas/MicromegasDefs.h>
0013
0014 #include <trackbase_historic/SvtxAlignmentState_v1.h>
0015 #include <trackbase_historic/SvtxAlignmentStateMap_v1.h>
0016 #include <trackbase_historic/SvtxTrack_v4.h>
0017 #include <trackbase_historic/SvtxTrackMap_v2.h>
0018 #include <trackbase_historic/SvtxTrackState_v3.h>
0019 #include <trackbase_historic/SvtxTrackSeed_v1.h>
0020 #include <trackbase_historic/SvtxTrackSeed_v2.h>
0021 #include <trackbase_historic/TrackSeed.h>
0022 #include <trackbase_historic/TrackSeed_v1.h>
0023 #include <trackbase_historic/TrackSeed_v2.h>
0024 #include <trackbase_historic/TrackSeedContainer.h>
0025 #include <trackbase_historic/WeightedTrackZeroField.h>
0026 #include <trackbase_historic/WeightedTrackMap.h>
0027
0028 #include <globalvertex/SvtxVertexMap_v1.h>
0029 #include <globalvertex/SvtxVertex_v2.h>
0030
0031 #include <fun4all/Fun4AllReturnCodes.h>
0032 #include <fun4all/Fun4AllServer.h>
0033
0034 #include <phool/getClass.h>
0035 #include <phool/PHCompositeNode.h>
0036 #include <phool/PHIODataNode.h>
0037 #include <phool/phool.h> // PHWHERE
0038
0039 #include <Acts/Definitions/Algebra.hpp>
0040 #include <Acts/Definitions/Units.hpp>
0041
0042 #include <TFile.h>
0043 #include <TNtuple.h>
0044
0045 #include <iostream>
0046 #include <functional>
0047 #include <limits>
0048 #include <map>
0049 #include <sstream>
0050 #include <stdexcept>
0051
0052 namespace { // anonymous
0053     template <typename T>
0054     T sqr (T const& t) { return t*t; }
0055 }
0056
0057 WeightedFitter::WeightedFitter (
0058     std::string const& name
0059 ) : SubsysReco (name),
0060     m_minimizer(ROOT::Math::Factory::CreateMinimizer("Minuit2"))
0061 {
0062     // Setup Minimizer
0063     m_minimizer->SetMaxFunctionCalls(1E8); // For Minuit/Minuit2
0064     m_minimizer->SetMaxIterations(1E6); // For GSL Minimizers--probably moot to call
0065     m_minimizer->SetTolerance(1E-8);
0066     m_minimizer->SetPrintLevel(0);
0067
0068     set_track_type<WeightedTrackZeroField>();
0069 }
0070
0071 WeightedFitter::~WeightedFitter (
0072 ) {
0073     // delete m_weighted_track; // m_weighted_track_map
0074     delete m_minimizer;
0075 }
0076
0077 void
0078 WeightedFitter::get_nodes (
0079     PHCompositeNode* top_node
0080 ) {
0081     std::vector<std::string> missing_node_names;
0082
0083     m_geometry = findNode::getClass<ActsGeometry>(top_node, m_geometry_node_name);
0084     if (!m_geometry) { missing_node_names.push_back(m_geometry_node_name); }
0085
0086     m_trkr_cluster_container = findNode::getClass<TrkrClusterContainer>(top_node, m_trkr_cluster_container_node_name);
0087     if (!m_trkr_cluster_container) { missing_node_names.push_back(m_trkr_cluster_container_node_name); }
0088
0089     if (m_which_tracks == k_svtx_tracks) {
0090         m_svtx_track_seed_container = findNode::getClass<TrackSeedContainer>(top_node, m_svtx_track_seed_container_node_name);
0091         if (!m_svtx_track_seed_container) { missing_node_names.push_back(m_svtx_track_seed_container_node_name); }
0092     }
0093
0094     if (m_which_tracks == k_svtx_tracks || m_which_tracks == k_silicon_tracks) {
0095         m_silicon_track_seed_container = findNode::getClass<TrackSeedContainer>(top_node, m_silicon_track_seed_container_node_name);
0096         if (!m_silicon_track_seed_container) { missing_node_names.push_back(m_silicon_track_seed_container_node_name); }
0097     }
0098
0099     if (m_which_tracks == k_svtx_tracks || m_which_tracks == k_tpc_tracks) {
0100         m_tpc_track_seed_container = findNode::getClass<TrackSeedContainer>(top_node, m_tpc_track_seed_container_node_name);
0101         if (!m_tpc_track_seed_container) { missing_node_names.push_back(m_tpc_track_seed_container_node_name); }
0102     }
0103
0104     if (m_use_vertex) {
0105         m_vertex_map = findNode::getClass<SvtxVertexMap>(top_node, m_vertex_map_node_name);
0106         if (!m_vertex_map) { missing_node_names.push_back(m_vertex_map_node_name); }
0107     }
0108
0109     if (!missing_node_names.empty()) {
0110         std::stringstream what;
0111         what
0112             << PHWHERE
0113             << " Couldn't get node(s): ";
0114         for (auto const& name : missing_node_names) { what << "\n " << name; }
0115         throw std::runtime_error(what.str());
0116     }
0117 }
0118
0119 void
0120 WeightedFitter::make_nodes (
0121     PHCompositeNode* top_node
0122 ) {
0123     PHNodeIterator dst_itr(top_node);
0124     PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(dst_itr.findFirst("PHCompositeNode", "DST"));
0125     if (!dst_node) {
0126         std::stringstream what;
0127         what
0128             << PHWHERE
0129             << " Couldn't get node: DST";
0130         throw std::runtime_error(what.str());
0131     }
0132
0133     PHCompositeNode* svtx_node = dynamic_cast<PHCompositeNode*>(dst_itr.findFirst("PHCompositeNode", "SVTX"));
0134     if (!svtx_node) {
0135         svtx_node = new PHCompositeNode("SVTX");
0136         dst_node->addNode(svtx_node);
0137     }
0138
0139     m_track_map = findNode::getClass<SvtxTrackMap>(top_node, m_track_map_node_name);
0140     if (!m_track_map) {
0141         m_track_map = new SvtxTrackMap_v2;
0142         PHIODataNode<PHObject>* track_map_node = new PHIODataNode<PHObject>(m_track_map, m_track_map_node_name, "PHObject");
0143         svtx_node->addNode(track_map_node);
0144     }
0145
0146     m_alignment_map = findNode::getClass<SvtxAlignmentStateMap>(top_node, m_alignment_map_node_name);
0147     if (!m_alignment_map) {
0148         m_alignment_map = new SvtxAlignmentStateMap_v1;
0149         PHIODataNode<PHObject>* alignment_map_node = new PHIODataNode<PHObject>(m_alignment_map, m_alignment_map_node_name, "PHObject");
0150         svtx_node->addNode(alignment_map_node);
0151     }
0152
0153     m_weighted_track_map = findNode::getClass<WeightedTrackMap>(top_node, m_weighted_track_map_node_name);
0154     if (!m_weighted_track_map) {
0155         m_weighted_track_map = new WeightedTrackMap;
0156         PHIODataNode<PHObject>* weighted_track_map_node = new PHIODataNode<PHObject>(m_weighted_track_map, m_weighted_track_map_node_name, "PHObject");
0157         svtx_node->addNode(weighted_track_map_node);
0158     }
0159 }
0160
0161 void
0162 WeightedFitter::make_ntuple (
0163 ) {
0164     if (m_ntuple_file_name.empty()) {
0165         if (m_file) { m_file->Close(); }
0166         m_file = nullptr;
0167
0168         delete m_ntuple;
0169         m_ntuple = nullptr;
0170
0171         return;
0172     }
0173
0174     m_file = TFile::Open(m_ntuple_file_name.c_str(), "RECREATE");
0175     if (!m_file) {
0176         std::stringstream what;
0177         what
0178             << PHWHERE
0179             << " Couldn't create file"
0180             << m_ntuple_file_name;
0181         throw std::runtime_error(what.str());
0182     }
0183
0184     delete m_ntuple;
0185     m_ntuple = new TNtuple (
0186         "ntp", "ntp",
0187         "event:track:"
0188         "fitstatus:"
0189         "nmaps:nintt:ntpc:nmms:"
0190         "cluslayer:"
0191         "clusstave:cluschip:clusstrobe:"
0192         "clusladderz:clusladderphi:clustimebucket:"
0193         "clusside:clussector:"
0194         "clussegtype:clustileid:"
0195         "cluslx:cluslz:cluselx:cluselz:"
0196         "clusgx:clusgy:clusgz:"
0197         "statelx:statelz:stateelx:stateelz:"
0198         "stategx:stategy:stategz:"
0199         "dXdx0:dXdy0:dXdtheta:dXdphi:"
0200         "dYdx0:dYdy0:dYdtheta:dYdphi:"
0201         "dca_xy:vtxgx:vtxgy:vtxgz:"
0202     );
0203     m_ntuple->SetDirectory(m_file);
0204 }
0205
0206 int
0207 WeightedFitter::InitRun (
0208     PHCompositeNode* top_node
0209 ) {
0210     m_global_position_wrapper.loadNodes(top_node);
0211     m_global_position_wrapper.set_suppressCrossing(true); // Suppresses print states in the case the crossing is SHRT_MAX
0212
0213     try {
0214         make_nodes(top_node);
0215         make_ntuple();
0216         return Fun4AllReturnCodes::EVENT_OK;
0217     } catch (std::exception const& e) {
0218         std::cout
0219             << e.what()
0220             << std::endl;
0221         return Fun4AllReturnCodes::ABORTEVENT;
0222     }
0223 }
0224
0225 int
0226 WeightedFitter::End (
0227     PHCompositeNode* /*top_node*/
0228 ) {
0229     if (m_ntuple && m_file) {
0230         m_ntuple->Write();
0231         m_file->Write();
0232         m_file->Close();
0233     }
0234
0235     return Fun4AllReturnCodes::EVENT_OK;
0236 }
0237
0238 int
0239 WeightedFitter::process_event (
0240     PHCompositeNode* top_node
0241 ) {
0242     try {
0243         get_nodes(top_node);
0244     } catch (std::exception const& e) {
0245         std::cout
0246             << e.what()
0247             << std::endl;
0248         return Fun4AllReturnCodes::ABORTEVENT;
0249     }
0250
0251     m_track_id = 0;
0252     m_track_map->Reset();
0253     m_alignment_map->Reset();
0254     m_weighted_track_map->Reset();
0255
0256     std::string track_seed_container_node_name{};
0257     TrackSeedContainer* track_seed_container{};
0258     switch (m_which_tracks) {
0259     case k_silicon_tracks:
0260         track_seed_container_node_name = m_silicon_track_seed_container_node_name;
0261         track_seed_container = m_silicon_track_seed_container;
0262         break;
0263     case k_tpc_tracks:
0264         track_seed_container_node_name = m_tpc_track_seed_container_node_name;
0265         track_seed_container = m_tpc_track_seed_container;
0266         break;
0267     case k_svtx_tracks:
0268         track_seed_container_node_name = m_svtx_track_seed_container_node_name;
0269         track_seed_container = m_svtx_track_seed_container;
0270         break;
0271     }
0272
0273     if (!track_seed_container) {
0274         if (Verbosity()) {
0275             std::cout
0276                 << PHWHERE
0277                 << " TrackSeedContainer " << track_seed_container_node_name << " is null"
0278                 << std::endl;
0279         }
0280         return Fun4AllReturnCodes::ABORTEVENT;
0281     }
0282
0283     if (Verbosity()) {
0284         std::cout
0285             << PHWHERE
0286             << " TrackSeedContainer " << track_seed_container_node_name
0287             << " identify(): ";
0288         track_seed_container->identify();
0289         std::cout << std::endl;
0290     }
0291
0292     for (auto* track_seed_ptr : *track_seed_container) {
0293         if (!track_seed_ptr) { continue; }
0294
0295         if (get_cluster_keys(track_seed_ptr)) {
0296             if (1 < Verbosity()) { std::cout << PHWHERE << "continuing" << std::endl; }
0297             continue;
0298         }
0299         if (get_points()) {
0300             if (1 < Verbosity()) { std::cout << PHWHERE << "continuing" << std::endl; }
0301             continue;
0302         }
0303         if (do_fit()) {
0304             if (1 < Verbosity()) { std::cout << PHWHERE << "continuing" << std::endl; }
0305             continue;
0306         }
0307         if (m_use_vertex && refit_with_vertex()) {
0308             if (1 < Verbosity()) { std::cout << PHWHERE << "continuing" << std::endl; }
0309             continue;
0310         }
0311         if (add_track()) {
0312             if (1 < Verbosity()) { std::cout << PHWHERE << "continuing" << std::endl; }
0313             continue;
0314         }
0315         ++m_track_id;
0316     }
0317
0318     if (Verbosity()) {
0319         std::cout
0320             << PHWHERE
0321             << " processed " << m_track_id << " tracks"
0322             << std::endl;
0323     }
0324
0325     return Fun4AllReturnCodes::EVENT_OK;
0326 }
0327
0328 bool
0329 WeightedFitter::get_cluster_keys (
0330     TrackSeed* track_seed
0331 ) {
0332     if (1 < Verbosity()) {
0333         std::cout
0334             << PHWHERE
0335             << " track seed identify: ";
0336         track_seed->identify();
0337         std::cout << std::endl;
0338     }
0339
0340     switch (m_which_tracks) {
0341     case k_silicon_tracks:
0342         m_silicon_seed = track_seed;
0343         m_tpc_seed = nullptr;
0344         break;
0345     case k_tpc_tracks:
0346         m_silicon_seed = nullptr;
0347         m_tpc_seed = track_seed;
0348         break;
0349     case k_svtx_tracks:
0350         m_silicon_seed = m_silicon_track_seed_container->get(track_seed->get_silicon_seed_index());
0351         m_tpc_seed = m_tpc_track_seed_container->get(track_seed->get_tpc_seed_index());
0352     }
0353     m_crossing = m_silicon_seed ? m_silicon_seed->get_crossing() : SHRT_MAX;
0354
0355     m_cluster_keys.clear();
0356     for (auto const* seed : {m_silicon_seed, m_tpc_seed}) {
0357         if (!seed) { continue; }
0358         std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(m_cluster_keys));
0359     }
0360
0361     return false;
0362 }
0363
0364 bool
0365 WeightedFitter::get_points (
0366 ) {
0367     m_num_mvtx = 0;
0368     m_num_intt = 0;
0369     m_num_tpc = 0;
0370     m_num_tpot = 0;
0371
0372     // delete m_weighted_track; // the track map owns the object
0373     m_weighted_track = m_track_factory();
0374     m_output_cluster_fit_points.clear();
0375
0376     // Assign all clusters (e.g., Si clusters) the same side
0377     // Choose mode as the value to set all points to
0378     // Print variations with high enough verbosity
0379     std::map<int, int> sides;
0380
0381     for (auto const& cluster_key : m_cluster_keys) {
0382         TrkrCluster* cluster = m_trkr_cluster_container->findCluster(cluster_key);
0383         if (!cluster) { continue; }
0384
0385         Surface const surf = m_geometry->maps().getSurface(cluster_key, cluster);
0386         if (!surf) { continue; }
0387
0388         auto local_to_global_transform = surf->transform(m_geometry->geometry().getGeoContext()); // in mm
0389         local_to_global_transform.translation() /= Acts::UnitConstants::cm; // converted to cm
0390         Eigen::Vector3d local_pos = Eigen::Vector3d { cluster->getLocalX(), cluster->getLocalY(), 0.0 }; // in cm
0391
0392         ClusterFitPoint point;
0393         point.cluster_key = cluster_key;
0394         point.cluster_position = (local_to_global_transform * local_pos);
0395         point.cluster_errors = { cluster->getRPhiError(), cluster->getZError() };
0396         point.sensor_local_to_global_transform = local_to_global_transform;
0397
0398         // Modify position, error using helper classes
0399         // Note that these calls are effectively NOOPs for non-TPC points
0400         point.cluster_position = m_global_position_wrapper.getGlobalPositionDistortionCorrected(cluster_key, cluster, m_crossing);
0401
0402         // note the "clusterRadius" argument (supplied as 0.0) is unused
0403         auto rphi_z_error_pair = ClusterErrorPara::get_clusterv5_modified_error(cluster, 0.0, cluster_key);
0404         point.cluster_errors = Eigen::Vector2d {std::sqrt(rphi_z_error_pair.first), std::sqrt(rphi_z_error_pair.second) };
0405
0406         int layer = TrkrDefs::getLayer(cluster_key);
0407
0408         // Include the point based on configuration handles
0409         // Both (AND) conditions must be met:
0410         // * We haven't specified a specific set of layers to include (so by default take all), OR we have explicitly specified to include the current layer
0411         // * We haven't explicitly specified to exclude the current layer
0412         if ( (m_fit_included_layers.empty() || m_fit_included_layers.contains(layer)) && !m_fit_excluded_layers.contains(layer) ) { m_weighted_track->push_back(point); }
0413         if ( (m_output_included_layers.empty() || m_output_included_layers.contains(layer)) && !m_output_excluded_layers.contains(layer) ) { m_output_cluster_fit_points.push_back(point); }
0414
0415         TrkrDefs::TrkrId trkr_id = static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(cluster_key));
0416         switch (trkr_id) {
0417             case TrkrDefs::mvtxId:
0418                 ++m_num_mvtx;
0419                 break;
0420             case TrkrDefs::inttId:
0421                 ++m_num_intt;
0422                 break;
0423             case TrkrDefs::tpcId:
0424                 ++m_num_tpc;
0425                 ++sides[TpcDefs::getSide(cluster_key)];
0426                 break;
0427             case TrkrDefs::micromegasId:
0428                 ++m_num_tpot;
0429                 break;
0430             default:
0431                 continue;
0432         }
0433     }
0434
0435     // Enforce minimum detector cluster requirements
0436     if (m_num_mvtx < m_min_num_mvtx) { return true; }
0437     if (m_num_intt < m_min_num_intt) { return true; }
0438     if (m_num_tpc < m_min_num_tpc) { return true; }
0439     if (m_num_tpot < m_min_num_tpot) { return true; }
0440
0441     if (m_reassign_sides && !sides.empty()) {
0442         m_side = sides[0] < sides[1] ? 1 : 0;
0443     }
0444
0445     return false;
0446 }
0447
0448 bool
0449 WeightedFitter::do_fit (
0450 ) {
0451     try {
0452         m_weighted_track->configure_minimizer(*m_minimizer);
0453     } catch (std::exception const& e) {
0454         if (Verbosity()) {
0455             std::cout
0456                 << PHWHERE << "\n"
0457                 << e.what()
0458                 << std::endl;
0459         }
0460         return true;
0461     }
0462
0463     std::function<double(double const*)> lambda = [&](double const* params) { return m_weighted_track->get_squared_error(params); };
0464     ROOT::Math::Functor functor(lambda, m_weighted_track->get_n_parameters());
0465     m_minimizer->SetFunction(functor);
0466
0467     bool fit_succeeded = m_minimizer->Minimize();
0468     if (1 < Verbosity()) {
0469         std::cout
0470             << PHWHERE
0471             << " fit " << (fit_succeeded ? "succeeded" : "failed")
0472             << " status: " << m_minimizer->Status()
0473             << std::endl;
0474     }
0475
0476     double const* params = m_minimizer->X();
0477     m_weighted_track->set_parameters(params);
0478
0479     //double const* errors = m_minimizer->Errors();
0480     int nparams = m_weighted_track->get_n_parameters();
0481     std::vector<double> cov_vector(nparams * nparams);
0482     m_minimizer->GetCovMatrix(cov_vector.data());
0483
0484     param_cov = Eigen::Matrix4d::Zero();
0485     for (int i = 0; i < nparams; ++i) {
0486         for (int j = 0; j < nparams; ++j) {
0487             param_cov(i, j) = cov_vector[i * nparams + j];
0488         }
0489     }
0490
0491     return !fit_succeeded;
0492 }
0493
0494 bool
0495 WeightedFitter::refit_with_vertex (
0496 ) {
0497     if (!m_vertex_map) {
0498         std::stringstream what;
0499         what
0500             << PHWHERE
0501             << " Couldn't get node: "
0502             << m_vertex_map_node_name;
0503         throw std::runtime_error(what.str());
0504     }
0505
0506     if (m_vertex_map->empty()) {
0507         if (1 < Verbosity()) {
0508             std::cout
0509                 << PHWHERE
0510                 << " vertex map is empty"
0511                 << std::endl;
0512         }
0513         return true;
0514     }
0515
0516     SvtxVertex* closest_vertex{};
0517     double min_dca = std::numeric_limits<double>::max();
0518     for (auto const& [vertex_id, svtx_vertex] : *m_vertex_map) {
0519         if (!svtx_vertex) { continue; }
0520
0521         Eigen::Vector3d vertex_pos {
0522             svtx_vertex->get_x(),
0523             svtx_vertex->get_y(),
0524             svtx_vertex->get_z(),
0525         };
0526
0527         double dca = (m_weighted_track->get_pca(vertex_pos) - vertex_pos).norm();
0528         if (dca < min_dca) {
0529             min_dca = dca;
0530             closest_vertex = svtx_vertex;
0531         }
0532     }
0533
0534     if (!closest_vertex) {
0535         std::stringstream what;
0536         what
0537             << PHWHERE
0538             << "no vertex after loop";
0539         throw std::runtime_error(what.str());
0540     }
0541
0542     m_weighted_track->use_vertex();
0543     for (int i = 0; i < 3; ++i) {
0544         m_weighted_track->vertex_position(i) = closest_vertex->get_position(i);
0545         for (int j = 0; j < 3; ++j) {
0546             m_weighted_track->vertex_covariance(i, j) = closest_vertex->get_error(i, j);
0547         }
0548     }
0549
0550     if (1 < Verbosity()) {
0551         std::cout
0552             << PHWHERE
0553             << " vertex: " << m_weighted_track->get_vertex_position().transpose()
0554             << " track pca: " << m_weighted_track->get_pca(m_weighted_track->get_vertex_position()).transpose()
0555             << " dca: " << min_dca
0556             << std::endl;
0557     }
0558
0559     return do_fit();
0560 }
0561
0562 bool
0563 WeightedFitter::add_track (
0564 ) {
0565     double const* params = m_minimizer->X();
0566     m_weighted_track->set_parameters(params);
0567     Eigen::Vector3d slope = m_weighted_track->get_slope_at_path_length(0); // Straight line tracks have uniform slope, we can pass any value for path_length here
0568
0569     SvtxTrack_v4 fitted_track;
0570     fitted_track.set_id(m_track_id);
0571     fitted_track.set_silicon_seed(m_silicon_seed);
0572     fitted_track.set_tpc_seed(m_tpc_seed);
0573     fitted_track.set_crossing(m_crossing);
0574     fitted_track.set_charge(m_tpc_seed ? m_tpc_seed->get_charge() : 0);
0575     fitted_track.set_x(params[0]);
0576     fitted_track.set_y(params[1]);
0577     fitted_track.set_z(0.0);
0578     fitted_track.set_px(slope(0));
0579     fitted_track.set_py(slope(1));
0580     fitted_track.set_pz(slope(2));
0581
0582     SvtxAlignmentStateMap::StateVec alignment_states;
0583     for (auto const& point : m_output_cluster_fit_points) {
0584         Acts::Vector3 intersection = m_weighted_track->get_intersection(point.sensor_local_to_global_transform);
0585         double path_length = m_weighted_track->get_path_length_of_intersection(point.sensor_local_to_global_transform);
0586
0587         SvtxTrackState_v3 svtx_track_state(path_length);
0588         svtx_track_state.set_x(intersection(0));
0589         svtx_track_state.set_y(intersection(1));
0590         svtx_track_state.set_z(intersection(2));
0591         svtx_track_state.set_px(slope(0));
0592         svtx_track_state.set_py(slope(1));
0593         svtx_track_state.set_pz(slope(2));
0594         svtx_track_state.set_name(std::to_string(point.cluster_key));
0595         svtx_track_state.set_cluskey(point.cluster_key);
0596
0597         Eigen::Matrix<double, 3, 4> Jacobian_fitpars_globpos;
0598         for (int i = 0; i < 4; ++i) { Jacobian_fitpars_globpos.col(i) = m_weighted_track->get_partial_derivative(i, path_length); }
0599         Eigen::Matrix3d globpos_cov = Jacobian_fitpars_globpos * param_cov * Jacobian_fitpars_globpos.transpose();
0600         for (int i = 0; i < 6; ++i) {
0601             for (int j = 0; j < 6; ++j) {
0602                 // use the rotated cluster uncertainties for the spatial component
0603                 // (these seem to be indices 0, 1, 2 judging by a comment in ActsTransformations.cc:187
0604                 if (i < 3 && j < 3) {
0605                     svtx_track_state.set_error(i, j, globpos_cov(i, j));
0606                 } else {
0607                     svtx_track_state.set_error(i, j, 0);
0608                 }
0609             }
0610         }
0611         fitted_track.insert_state(&svtx_track_state);
0612
0613         auto alignment_state = std::make_unique<SvtxAlignmentState_v1>();
0614         SvtxAlignmentState::GlobalMatrix global_derivative_matrix = SvtxAlignmentState::GlobalMatrix::Zero();
0615         SvtxAlignmentState::LocalMatrix local_derivative_matrix = SvtxAlignmentState::LocalMatrix::Zero();
0616
0617         Eigen::Vector2d residual = point.get_residuals(intersection);
0618         Eigen::Matrix<double, 2, 3> projection = m_weighted_track->get_projection(point.sensor_local_to_global_transform);
0619
0620         // PARTIAL derivatives of global alignment parameters
0621         // (the projection vectors take care of the implicit intersection dependency)
0622         std::array<Eigen::Vector3d, 6> global_derivatives {
0623             Eigen::Vector3d {1.0, 0.0, 0.0}.cross(intersection - point.sensor_local_to_global_transform.translation()),
0624             Eigen::Vector3d {0.0, 1.0, 0.0}.cross(intersection - point.sensor_local_to_global_transform.translation()),
0625             Eigen::Vector3d {0.0, 0.0, 1.0}.cross(intersection - point.sensor_local_to_global_transform.translation()),
0626             Eigen::Vector3d {1.0, 0.0, 0.0},
0627             Eigen::Vector3d {0.0, 1.0, 0.0},
0628             Eigen::Vector3d {0.0, 0.0, 1.0},
0629         };
0630
0631         for (int i = 0; i < 4; ++i) {
0632             local_derivative_matrix.col(i) = m_weighted_track->get_residual_derivative(i, point.sensor_local_to_global_transform);
0633         }
0634
0635         for (int i = 0; i < 6; ++i) {
0636             global_derivative_matrix.col(i) = projection * global_derivatives[i];
0637         }
0638
0639         alignment_state->set_cluster_key(point.cluster_key);
0640         alignment_state->set_residual(residual);
0641         alignment_state->set_local_derivative_matrix(local_derivative_matrix);
0642         alignment_state->set_global_derivative_matrix(global_derivative_matrix);
0643         alignment_states.push_back(alignment_state.release());
0644
0645         if (m_ntuple) {
0646             if (1 < Verbosity()) {
0647                 std::cout
0648                     << " Filling ntuple"
0649                     << std::endl;
0650             }
0651
0652             float layer = TrkrDefs::getLayer(point.cluster_key);
0653             float stave{-1};
0654             float chip{-1};
0655             float strobe{-1};
0656             float ladder_z{-1};
0657             float ladder_phi{-1};
0658             float time_bucket{-1};
0659             float side{-1};
0660             float sector{-1};
0661             float segtype{-1};
0662             float tileid{-1};
0663             TrkrDefs::TrkrId trkr_id = static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(point.cluster_key));
0664             switch (trkr_id) {
0665                 case TrkrDefs::mvtxId:
0666                     stave = MvtxDefs::getStaveId(point.cluster_key);
0667                     chip = MvtxDefs::getChipId(point.cluster_key);
0668                     strobe = MvtxDefs::getStrobeId(point.cluster_key);
0669                     break;
0670                 case TrkrDefs::inttId:
0671                     ladder_z = InttDefs::getLadderZId(point.cluster_key);
0672                     ladder_phi = InttDefs::getLadderPhiId(point.cluster_key);
0673                     time_bucket = InttDefs::getTimeBucketId(point.cluster_key);
0674                     break;
0675                 case TrkrDefs::tpcId:
0676                     side = TpcDefs::getSide(point.cluster_key);
0677                     sector = TpcDefs::getSectorId(point.cluster_key);
0678                     break;
0679                 case TrkrDefs::micromegasId:
0680                     segtype = (float) MicromegasDefs::getSegmentationType(point.cluster_key);
0681                     tileid = MicromegasDefs::getTileId(point.cluster_key);
0682                     break;
0683                 default:
0684                     continue;
0685             }
0686             if (m_reassign_sides) { side = m_side; }
0687
0688             Eigen::Affine3d global_to_local_transform = point.sensor_local_to_global_transform.inverse();
0689             Eigen::Vector3d cluster_local_position =  global_to_local_transform * point.cluster_position;
0690             Eigen::Vector3d intersection_local_position = global_to_local_transform * intersection;
0691
0692             float dca_xy = std::numeric_limits<float>::quiet_NaN();
0693             Eigen::Vector3f vertex = {
0694                 std::numeric_limits<float>::quiet_NaN(),
0695                 std::numeric_limits<float>::quiet_NaN(),
0696                 std::numeric_limits<float>::quiet_NaN(),
0697             };
0698
0699             if (m_use_vertex) {
0700                 vertex(0) = (float)m_weighted_track->vertex_position(0);
0701                 vertex(1) = (float)m_weighted_track->vertex_position(1);
0702                 vertex(2) = (float)m_weighted_track->vertex_position(2);
0703
0704                 // displacement between the associated vertex and the point of closest approach on the track
0705                 Eigen::Vector3d delta = m_weighted_track->get_pca(m_weighted_track->get_vertex_position()) - m_weighted_track->get_vertex_position();
0706                 dca_xy = std::sqrt(sqr(delta(0)) + sqr(delta(1)));
0707             }
0708
0709             float ntp_data[] = {
0710                 (float)Fun4AllServer::instance()->EventNumber(), (float)m_track_id,
0711                 (float)m_minimizer->Status(),
0712                 (float)m_num_mvtx, (float)m_num_intt, (float)m_num_tpc, (float)m_num_tpot,
0713                 layer,
0714                 stave, chip, strobe,
0715                 ladder_z, ladder_phi, time_bucket,
0716                 side, sector,
0717                 segtype, tileid,
0718                 (float)cluster_local_position(0), (float)cluster_local_position(1), (float)point.cluster_errors(0), (float)point.cluster_errors(1),
0719                 (float)point.cluster_position(0), (float)point.cluster_position(1), (float)point.cluster_position(2),
0720                 (float)intersection_local_position(0), (float)intersection_local_position(1), (float)point.cluster_errors(0), (float)point.cluster_errors(1),
0721                 (float)intersection(0), (float)intersection(1), (float)intersection(2),
0722                 (float)local_derivative_matrix(0, 0), (float)local_derivative_matrix(0, 1), (float)local_derivative_matrix(0, 2), (float)local_derivative_matrix(0, 3),
0723                 (float)local_derivative_matrix(1, 0), (float)local_derivative_matrix(1, 1), (float)local_derivative_matrix(1, 2), (float)local_derivative_matrix(1, 3),
0724                 dca_xy, vertex(0), vertex(1), vertex(2),
0725             };
0726             m_ntuple->Fill(ntp_data);
0727         }
0728
0729         if (2 < Verbosity()) {
0730             std::cout
0731                 << "\tcluster_key:  " << point.cluster_key << "\n"
0732                 << "\tlayer:        " << (int)TrkrDefs::getLayer(point.cluster_key) << "\n"
0733                 << "\tintersection: " << intersection.transpose() << "\n"
0734                 << "\tpos:          " << point.cluster_position.transpose() << "\n"
0735                 << "\to:            " << point.sensor_local_to_global_transform.translation().transpose() << "\n"
0736                 << "\tx:            " << point.sensor_local_to_global_transform.rotation().col(0).transpose() << "\n"
0737                 << "\ty:            " << point.sensor_local_to_global_transform.rotation().col(1).transpose() << "\n"
0738                 << "\tz:            " << point.sensor_local_to_global_transform.rotation().col(2).transpose() << "\n"
0739                 << "\tproj_x:       " << projection.row(0) << "\n"
0740                 << "\tproj_y:       " << projection.row(1) << "\n"
0741                 << std::endl;
0742         }
0743     }
0744
0745     // cuts...
0746     // return true;
0747
0748     m_track_map->insertWithKey(&fitted_track, m_track_id);
0749     m_alignment_map->insertWithKey(m_track_id, alignment_states);
0750     m_weighted_track_map->insert({m_track_id, m_weighted_track});
0751
0752     return false;
0753 }
0754