packages/tpccalib/MicromegasDriftEvaluator.cc

0001 #include "MicromegasDriftEvaluator.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <g4detectors/PHG4CylinderGeomContainer.h>
0005 #include <micromegas/CylinderGeomMicromegas.h>
0006 #include <micromegas/MicromegasDefs.h>
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/PHNodeIterator.h>
0009 #include <phool/getClass.h>
0010 #include <tpc/TpcGlobalPositionWrapper.h>
0011 #include <trackbase/ActsGeometry.h>
0012 #include <trackbase/TrackFitUtils.h>
0013 #include <trackbase/TrkrCluster.h>
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase/TrkrDefs.h>
0016 #include <trackbase_historic/SvtxTrack.h>
0017 #include <trackbase_historic/SvtxTrackMap.h>
0018
0019 #include <TCanvas.h>
0020 #include <TDirectory.h>
0021 #include <TF1.h>
0022 #include <TF2.h>
0023 #include <TFile.h>
0024 #include <TH1.h>
0025 #include <TH2.h>
0026 #include <TH3.h>
0027 #include <TLegend.h>
0028 #include <TParameter.h>
0029 #include <TVector3.h>
0030
0031 #include <cassert>
0032 #include <climits>
0033 #include <cmath>
0034 #include <format>
0035 #include <iostream>
0036 #include <memory>
0037
0038 namespace
0039 {
0040
0041   template <class T>
0042   class range_adaptor
0043   {
0044    public:
0045     explicit range_adaptor(const T& range)
0046       : m_range(range)
0047     {
0048     }
0049     const typename T::first_type& begin() { return m_range.first; }
0050     const typename T::second_type& end() { return m_range.second; }
0051
0052    private:
0053     T m_range;
0054   };
0055
0056   template <class T>
0057   constexpr T square(T x)
0058   {
0059     return x * x;
0060   }
0061   template <class T>
0062   inline T get_r(T x, T y)
0063   {
0064     return std::sqrt(square(x) + square(y));
0065   }
0066
0067   double normalize_angle(double phi)
0068   {
0069     while (phi < 0)
0070     {
0071       phi += 2 * M_PI;
0072     }
0073     while (phi >= 2 * M_PI)
0074     {
0075       phi -= 2 * M_PI;
0076     }
0077     return phi;
0078   }
0079
0080   bool phi_in_range(double phi, double min, double max)
0081   {
0082     phi = normalize_angle(phi);
0083     min = normalize_angle(min);
0084     max = normalize_angle(max);
0085     return (min < max) ? (phi >= min && phi <= max)
0086                        : (phi >= min || phi <= max);
0087   }
0088
0089   //  This function is identical to the version in MicromegasTrackEvaluator_hp.cc
0090
0091   bool helix_plane_intersection(
0092       double t_min,
0093       double t_max,
0094       double zmin,
0095       double zmax,
0096       double R,
0097       double X0,
0098       double Y0,
0099       double intersect_rz,
0100       double slope_rz,
0101       const TVector3& ptile,
0102       const TVector3& ntile,
0103       TVector3& intersect)
0104   {
0105     // Number of iterations and tolerance for Newton Raphson method
0106     const int max_iter = 10;
0107     const double tol = 1e-6;
0108
0109     // Define C
0110     double C = ntile.X() * (X0 - ptile.X()) + ntile.Y() * (Y0 - ptile.Y()) + ntile.Z() * (intersect_rz - ptile.Z());
0111
0112     // Defines the function and the corresponding derivative to be used in the Newton Raphson method
0113     auto f = [&](double t)
0114     {
0115       double xt = X0 + R * std::cos(t);
0116       double yt = Y0 + R * std::sin(t);
0117       double Rt = std::sqrt(xt * xt + yt * yt);
0118       return ntile.X() * R * std::cos(t) + ntile.Y() * R * std::sin(t) + ntile.Z() * slope_rz * Rt + C;
0119     };
0120
0121     auto df = [&](double t)
0122     {
0123       double xt = X0 + R * std::cos(t);
0124       double yt = Y0 + R * std::sin(t);
0125       double Rt = std::sqrt(xt * xt + yt * yt);
0126       return -ntile.X() * R * std::sin(t) + ntile.Y() * R * std::cos(t) + ntile.Z() * R * slope_rz * (Y0 * std::cos(t) - X0 * std::sin(t)) / Rt;
0127     };
0128
0129     auto solve_from = [&](double t_seed, TVector3& result) -> bool
0130     {
0131       double t = t_seed;
0132       for (int i = 0; i < max_iter; ++i)
0133       {
0134         double ft = f(t);
0135         double dft = df(t);
0136         if (std::abs(dft) < 1e-8)
0137         {
0138           return false;
0139         }
0140         double t_new = t - ft / dft;
0141
0142         double x = X0 + R * std::cos(t_new);
0143         double y = Y0 + R * std::sin(t_new);
0144         double Rt_n = std::sqrt(x * x + y * y);
0145         double z = slope_rz * Rt_n + intersect_rz;
0146         double phi = std::atan2(y, x);
0147
0148         TVector3 cand(x, y, z);
0149         bool phi_ok = phi_in_range(phi, t_min - 1e-4, t_max + 1e-4);
0150         bool z_ok = (z >= zmin - 1e-4 && z <= zmax + 1e-4);
0151         bool proj_ok = (std::abs(ntile.Dot(cand - ptile)) <= 0.05);
0152
0153         if (std::abs(t_new - t) < tol && proj_ok && phi_ok && z_ok)
0154         {
0155           result = cand;
0156           return true;
0157         }
0158         t = t_new;
0159       }
0160       return false;
0161     };
0162
0163     auto wrap = [&](double t)
0164     {
0165       while (t > t_max)
0166       {
0167         t -= 2 * M_PI;
0168       }
0169       while (t < t_min)
0170       {
0171         t += 2 * M_PI;
0172       }
0173       return t;
0174     };
0175
0176     std::vector<double> t_seeds;
0177     double t_center = 0.5 * (t_min + t_max);
0178     double delta = 2.0 * M_PI / 3.0;
0179
0180     // Wrap the angle
0181     for (int i = 0; i < 3; ++i)
0182     {
0183       double t = wrap(t_center + i * delta);
0184       t_seeds.push_back(t);
0185     }
0186
0187     // Looks for the solution within the tile acceptance in three different phi seeds in the Newton-Raphson (helix_plane could have more than one solution)
0188     for (double t_seed : t_seeds)
0189     {
0190       if (solve_from(t_seed, intersect))
0191       {
0192         return true;
0193       }
0194     }
0195     return false;
0196   }
0197
0198   // this is a piecewise fit function for the drift velocity plot
0199   double fit_function_2d(double* x, double* par)
0200   {
0201     const int itile = static_cast<int>(std::floor(x[0]));
0202     const double z = x[1];
0203     if (itile < 0 || itile >= 8)
0204     {
0205       TF2::RejectPoint();
0206       return 0.;
0207     }
0208     return par[itile + 1] + par[0] * z;
0209   }
0210
0211 // root fitting does not like const parameters suggested by clang-tidy
0212 // using NOLINT to suppress this warning
0213   double linear_function(double* x, double* par) // NOLINT(readability-non-const-parameter)
0214   {
0215     return par[0] * x[0] + par[1];
0216   }
0217
0218   const std::array<const char*, 8> k_tile_names =
0219       {"SCOZ", "SCIZ", "NCIZ", "NCOZ", "SEZ", "NEZ", "SWZ", "NWZ"};
0220
0221 }  // namespace
0222
0223 MicromegasDriftEvaluator::MicromegasDriftEvaluator(const std::string& name)
0224   : SubsysReco(name)
0225 {
0226 }
0227
0228 // ---------------------------------------------------------------------------
0229 int MicromegasDriftEvaluator::Init(PHCompositeNode* topNode)
0230 {
0231   std::cout << Name() << "::Init"
0232             << " drift_velocity=" << m_drift_velocity << " cm/ns"
0233             << " min_tpc_layer=" << m_min_tpc_layer
0234             << " max_tpc_layer=" << m_max_tpc_layer
0235             << std::endl;
0236
0237   PHNodeIterator iter(topNode);
0238   auto* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0239   if (!dstNode)
0240   {
0241     std::cerr << Name() << "::Init - DST node missing" << std::endl;
0242     return Fun4AllReturnCodes::ABORTEVENT;
0243   }
0244
0245   iter = PHNodeIterator(dstNode);
0246   auto* evalNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "EVAL"));
0247   if (!evalNode)
0248   {
0249     evalNode = new PHCompositeNode("EVAL");
0250     dstNode->addNode(evalNode);
0251   }
0252
0253   auto* newNode = new PHIODataNode<PHObject>(new Container, "MicromegasDriftEvaluator::Container", "PHObject");
0254   newNode->SplitLevel(99);
0255   evalNode->addNode(newNode);
0256
0257   m_hist3D = new TH3F("MicromegasDriftEval_hist3D", ";tile;z_{track} (cm);#Deltaz (track#minuscluster) (cm)", 8, 0, 8, 220, -110, 110, 100, -10, 10);
0258   m_hist3D->SetDirectory(nullptr);
0259
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263 // ---------------------------------------------------------------------------
0264 int MicromegasDriftEvaluator::InitRun(PHCompositeNode* topNode)
0265 {
0266   return load_nodes(topNode);
0267 }
0268
0269 // ---------------------------------------------------------------------------
0270 int MicromegasDriftEvaluator::process_event(PHCompositeNode* topNode)
0271 {
0272   const auto res = load_nodes(topNode);
0273   if (res != Fun4AllReturnCodes::EVENT_OK)
0274   {
0275     return res;
0276   }
0277
0278   if (m_container)
0279   {
0280     m_container->Reset();
0281   }
0282   evaluate_tracks();
0283
0284   return Fun4AllReturnCodes::EVENT_OK;
0285 }
0286
0287 // ---------------------------------------------------------------------------
0288 int MicromegasDriftEvaluator::End(PHCompositeNode* /*topNode*/)
0289 {
0290   if (!m_hist3D)
0291   {
0292     std::cerr << Name() << "::End - histogram not found, skipping fit." << std::endl;
0293     return Fun4AllReturnCodes::EVENT_OK;
0294   }
0295
0296   const int nEntries = static_cast<int>(m_hist3D->GetEntries());
0297   std::cout << Name() << "::End - fitting " << nEntries << " entries" << std::endl;
0298
0299   auto* h_fit = new TH2F("h_fit", "", 8, 0, 8, 220, -110, 110);
0300   h_fit->SetDirectory(nullptr);
0301
0302   for (int j = 0; j < 8; ++j)
0303   {
0304     m_hist3D->GetXaxis()->SetRange(j + 1, j + 1);
0305     auto* h2d = static_cast<TH2F*>(m_hist3D->Project3D("zy"));
0306     h2d->SetName(std::format("h_{}", k_tile_names[j]).c_str());
0307     h2d->SetDirectory(nullptr);
0308
0309     // Fit vertical slices; require a minimum of 10 entries per slice
0310     h2d->FitSlicesY(nullptr, 0, -1, 10);
0311     auto* h_mean = static_cast<TH1F*>(gDirectory->Get(std::format("h_{}_1", k_tile_names[j]).c_str()));
0312
0313     if (!h_mean)
0314     {
0315       delete h2d;
0316       continue;
0317     }
0318
0319     for (int i = 0; i < h_mean->GetNbinsX(); ++i)
0320     {
0321       const double entries = h2d->Integral(i + 1, i + 1, 1, m_hist3D->GetNbinsZ());
0322       if (entries > 0)
0323       {
0324         h_fit->SetBinContent(j + 1, i + 1, h_mean->GetBinContent(i + 1));
0325       }
0326     }
0327     delete h2d;
0328   }
0329
0330   m_hist3D->GetXaxis()->SetRange(0, 0);
0331
0332   // This part fits the 8 micromegas modules one at a time but constrains the slopes to be identical. This eliminates the need for perfect translational TPOT alignment
0333   auto* fit2d = new TF2("fit2d", fit_function_2d, 0, 8, -110, 110, 9);
0334   for (int i = 0; i < 9; ++i)
0335   {
0336     fit2d->SetParameter(i, 0.0);
0337   }
0338
0339   h_fit->Fit(fit2d, "0R");
0340
0341   const double slope = fit2d->GetParameter(0);
0342   const double slope_err = fit2d->GetParError(0);
0343   const double new_drift = m_drift_velocity / (1.0 + slope);
0344   const double drift_err = m_drift_velocity / std::pow(1.0 + slope, 2) * slope_err;
0345
0346   std::cout << Name() << "::End" << " slope=" << slope << " input_drift=" << m_drift_velocity << " cm/ns" << " new_drift=" << new_drift << " cm/ns +/- " << drift_err << " cm/ns" << std::endl;
0347
0348   // Plot the whole thing
0349   auto* canvas = new TCanvas("drift_calib", "Drift velocity calibration", 2000, 1000);
0350   canvas->Divide(4, 2);
0351
0352   for (int j = 0; j < 8; ++j)
0353   {
0354     canvas->cd(j + 1);
0355
0356     m_hist3D->GetXaxis()->SetRange(j + 1, j + 1);
0357     auto* h2d = static_cast<TH2F*>(m_hist3D->Project3D("zy"));
0358     h2d->SetName(std::format("hplot_{}", k_tile_names[j]).c_str());
0359     h2d->SetTitle(std::format("{};z_{{track}} (cm);#Deltaz (track#minuscluster) (cm)", k_tile_names[j]).c_str());
0360     h2d->SetStats(false);
0361     h2d->Draw("COLZ");
0362
0363     auto* h_fit_proj = h_fit->ProjectionY(std::format("h_fit_proj_{}", j).c_str(), j + 1, j + 1);  // These give you the Gaussian means for each slice
0364     h_fit_proj->SetMarkerStyle(20);
0365     h_fit_proj->SetMarkerColor(kRed);
0366     h_fit_proj->SetLineColor(kBlack);
0367
0368     auto* f1d = new TF1(std::format("f1d_{}", j).c_str(), linear_function, -110, 110, 2);
0369     f1d->SetParameter(0, slope);
0370     f1d->SetParameter(1, fit2d->GetParameter(j + 1));
0371     f1d->SetLineColor(kGreen + 2);
0372     f1d->SetLineWidth(2);
0373     f1d->Draw("same");
0374
0375     auto* leg = new TLegend(0.35, 0.75, 0.92, 0.92);
0376     leg->SetHeader(std::format("{} entries, v_{{in}}={:.2f} m/ms", nEntries, m_drift_velocity * 1e4).c_str(), "C");
0377     leg->AddEntry(h_fit_proj, "Gaussian slice mean", "p");
0378     leg->AddEntry(f1d, std::format("slope={:.4f}  v_{{new}}={:.3f}#pm{:.3f} m/ms", slope, new_drift * 1e4, drift_err * 1e4).c_str(), "l");
0379     leg->Draw();
0380   }
0381
0382   m_hist3D->GetXaxis()->SetRange(0, 0);
0383
0384   canvas->SaveAs(m_plot_filename.c_str());
0385   std::cout << Name() << "::End - QA canvas saved to " << m_plot_filename << std::endl;
0386
0387   // write histograms, fit and results to a ROOT file
0388   if (!m_root_filename.empty())
0389   {
0390     std::unique_ptr<TFile> outfile(TFile::Open(m_root_filename.c_str(), "RECREATE"));
0391     if (outfile && !outfile->IsZombie())
0392     {
0393       outfile->cd();
0394       m_hist3D->Write();
0395       h_fit->Write("h_fit_micromegas");
0396       fit2d->Write();
0397       canvas->Write();
0398       TParameter<double>("slope", slope).Write();
0399       TParameter<double>("drift_velocity_in", m_drift_velocity).Write();
0400       TParameter<double>("drift_velocity_new", new_drift).Write();
0401       TParameter<double>("drift_velocity_err", drift_err).Write();
0402       outfile->Close();
0403       std::cout << Name() << "::End - histograms and fit results saved to " << m_root_filename << std::endl;
0404     }
0405     else
0406     {
0407       std::cerr << Name() << "::End - could not open " << m_root_filename << " for writing." << std::endl;
0408     }
0409   }
0410
0411   delete canvas;
0412   delete fit2d;
0413   delete h_fit;
0414
0415   return Fun4AllReturnCodes::EVENT_OK;
0416 }
0417
0418 // ---------------------------------------------------------------------------
0419 int MicromegasDriftEvaluator::load_nodes(PHCompositeNode* topNode)
0420 {
0421   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0422   assert(m_tGeometry);
0423
0424   m_micromegas_geomcontainer = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0425   assert(m_micromegas_geomcontainer);
0426
0427   m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_trackmapname);
0428
0429   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0430   assert(m_cluster_map);
0431
0432   m_container = findNode::getClass<Container>(topNode, "MicromegasDriftEvaluator::Container");
0433   assert(m_container);
0434
0435   m_globalPositionWrapper.loadNodes(topNode);
0436
0437   return Fun4AllReturnCodes::EVENT_OK;
0438 }
0439
0440 // ---------------------------------------------------------------------------
0441 void MicromegasDriftEvaluator::evaluate_tracks()
0442 {
0443   if (!(m_tGeometry && m_micromegas_geomcontainer && m_track_map && m_cluster_map && m_container && m_hist3D))
0444   {
0445     return;
0446   }
0447
0448   m_container->clear_tracks();
0449
0450   for (const auto& [track_id, track] : *m_track_map)
0451   {
0452     // valid crossing
0453     const auto crossing = track->get_crossing();
0454     if (crossing == SHRT_MAX)
0455     {
0456       continue;
0457     }
0458
0459     std::vector<Acts::Vector3> tpc_positions;
0460
0461     // Also count clusters per subsystem for the cuts
0462     unsigned int n_tpc = 0;
0463     unsigned int n_mvtx = 0;
0464     unsigned int n_intt = 0;
0465     unsigned int n_mm = 0;
0466
0467     for (const auto* seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0468     {
0469       if (!seed)
0470       {
0471         continue;
0472       }
0473       for (auto it = seed->begin_cluster_keys(); it != seed->end_cluster_keys(); ++it)
0474       {
0475         const auto ckey = *it;
0476         const auto detid = TrkrDefs::getTrkrId(ckey);
0477         const auto layer = TrkrDefs::getLayer(ckey);
0478
0479         switch (detid)
0480         {
0481         case TrkrDefs::tpcId:
0482           ++n_tpc;
0483           if (layer >= m_min_tpc_layer && layer < m_max_tpc_layer)
0484           {
0485             auto* const cl = m_cluster_map->findCluster(ckey);
0486             if (cl)
0487             {
0488               tpc_positions.push_back(
0489                   m_globalPositionWrapper.getGlobalPositionDistortionCorrected(
0490                       ckey, cl, crossing));
0491             }
0492           }
0493           break;
0494         case TrkrDefs::mvtxId:
0495           ++n_mvtx;
0496           break;
0497         case TrkrDefs::inttId:
0498           ++n_intt;
0499           break;
0500         case TrkrDefs::micromegasId:
0501           ++n_mm;
0502           break;
0503         default:
0504           break;
0505         }
0506       }
0507     }
0508
0509     // need at least 3 TPC clusters in range
0510     if (tpc_positions.size() < 3)
0511     {
0512       continue;
0513     }
0514
0515     const auto [slope_rz, intersect_rz] = TrackFitUtils::line_fit(tpc_positions);
0516     const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(tpc_positions);
0517
0518     // reject badly reconstructed / low-pT tracks
0519     if (R < 40.0)
0520     {
0521       continue;
0522     }
0523
0524     const auto mm_range = m_micromegas_geomcontainer->get_begin_end();
0525     for (const auto& [mm_layer, base_layergeom] : range_adaptor(mm_range))
0526     {
0527       const auto* layergeom = static_cast<const CylinderGeomMicromegas*>(base_layergeom);
0528       assert(layergeom);
0529
0530       // skip the phi layer. Only the z-view layer matters here
0531       if (layergeom->get_segmentation_type() !=
0532           MicromegasDefs::SegmentationType::SEGMENTATION_Z)
0533       {
0534         continue;
0535       }
0536
0537       const double layer_radius = layergeom->get_radius();
0538       auto [xplus, yplus, xminus, yminus] =
0539           TrackFitUtils::circle_circle_intersection(layer_radius, R, X0, Y0);
0540
0541       if (!std::isfinite(xplus))
0542       {
0543         continue;
0544       }
0545
0546       // pick the solution closest in phi to the last TPC cluster
0547       const double last_phi = std::atan2(tpc_positions.back().y(), tpc_positions.back().x());
0548       const double phi_plus = std::atan2(yplus, xplus);
0549       const double phi_minus = std::atan2(yminus, xminus);
0550       const double phi = (std::abs(last_phi - phi_plus) < std::abs(last_phi - phi_minus)) ? phi_plus : phi_minus;
0551
0552       const double r_cyl = layer_radius;
0553       const double z_cyl = intersect_rz + slope_rz * r_cyl;
0554       const TVector3 world_cyl(r_cyl * std::cos(phi), r_cyl * std::sin(phi), z_cyl);
0555
0556       const int tileid = layergeom->find_tile_cylindrical(world_cyl);
0557       if (tileid < 0)
0558       {
0559         continue;
0560       }
0561
0562       const auto tile_center = layergeom->get_world_from_local_coords(tileid, m_tGeometry, {0, 0});
0563       const TVector3 ptile(tile_center.x(), tile_center.y(), tile_center.z());
0564
0565       const auto tile_norm = layergeom->get_world_from_local_vect(tileid, m_tGeometry, {0, 0, 1});
0566       const TVector3 ntile(tile_norm.x(), tile_norm.y(), tile_norm.z());
0567
0568       const auto phi_range = layergeom->get_phi_range(tileid, m_tGeometry);
0569       const double zmin = layergeom->get_zmin();
0570       const double zmax = layergeom->get_zmax();
0571
0572       TVector3 intersection;
0573       if (!helix_plane_intersection(phi_range.first, phi_range.second, zmin, zmax, R, X0, Y0, intersect_rz, slope_rz, ptile, ntile, intersection))
0574       {
0575         continue;
0576       }
0577
0578       const auto local_intersection = layergeom->get_local_from_world_coords(tileid, m_tGeometry, {intersection.x(), intersection.y(), intersection.z()});
0579       const double y_local = local_intersection.y();
0580
0581       if (std::abs(y_local) > m_y_local_cut)
0582       {
0583         continue;
0584       }
0585
0586       // find the nearest TPOT cluster
0587       const auto hitsetkey = MicromegasDefs::genHitSetKey(mm_layer, MicromegasDefs::SegmentationType::SEGMENTATION_Z, tileid);
0588       const auto clusrange = m_cluster_map->getClusters(hitsetkey);
0589
0590       double dmin = -1;
0591       ClusterStruct best_cluster;
0592
0593       for (const auto& [ckey, cl] : range_adaptor(clusrange))
0594       {
0595         const double cl_y_local = cl->getLocalY();
0596         const double d = std::abs(y_local - cl_y_local);
0597         if (dmin < 0 || d < dmin)
0598         {
0599           dmin = d;
0600           const auto gpos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cl, crossing);
0601           best_cluster._layer = mm_layer;
0602           best_cluster._tile = tileid;
0603           best_cluster._z = gpos.z();
0604         }
0605       }
0606
0607       // require cluster within the z search window
0608       if (dmin < 0 || dmin > m_z_search_win)
0609       {
0610         continue;
0611       }
0612
0613       // fill track struct and histogram
0614       TrackStruct track_struct;
0615       track_struct._chisquare = track->get_chisq();
0616       track_struct._ndf = track->get_ndf();
0617       track_struct._nclusters_tpc = n_tpc;
0618       track_struct._nclusters_mvtx = n_mvtx;
0619       track_struct._nclusters_intt = n_intt;
0620       track_struct._nclusters_micromegas = n_mm;
0621
0622       track_struct._trk_state_z._layer = mm_layer;
0623       track_struct._trk_state_z._tile = tileid;
0624       track_struct._trk_state_z._z = intersection.z();
0625       track_struct._trk_state_z._y_local = y_local;
0626
0627       track_struct._found_cluster_z = best_cluster;
0628
0629       const double z_track = track_struct._trk_state_z._z;
0630       const double z_cluster = track_struct._found_cluster_z._z;
0631       m_hist3D->Fill(tileid + 0.5, z_track, z_track - z_cluster);
0632
0633       m_container->add_track(track_struct);
0634       break;
0635     }
0636   }
0637 }