QA/Tracking/SiliconDriftQA.cc

0001 #include "SiliconDriftQA.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004
0005 #include <fun4all/Fun4AllHistoManager.h>
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/getClass.h>
0010 #include <phool/phool.h>  // for PHWHERE
0011
0012 #include <trackbase/TrkrDefs.h>
0013 #include <trackbase_historic/SvtxTrack.h>
0014 #include <trackbase_historic/SvtxTrackMap.h>
0015 #include <trackbase_historic/TrackSeedHelper.h>
0016
0017 #include <TDirectory.h>
0018 #include <TF1.h>
0019 #include <TF2.h>
0020 #include <TH1.h>
0021 #include <TH2.h>
0022
0023 #include <cassert>
0024 #include <climits>
0025 #include <cmath>
0026 #include <format>
0027 #include <iostream>
0028 #include <string>
0029
0030 namespace
0031 {
0032   //! pt
0033   template <class T>
0034   T get_pt(const T& px, const T& py)
0035   {
0036     return std::sqrt(px * px + py * py);
0037   }
0038
0039   //! piecewise fit function used for the drift velocity extraction
0040   //  par[0] = constrained slope
0041   //  par[1] = offset for eta < 0
0042   //  par[2] = offset for eta >= 0
0043   double fit_function_2d(double* x, double* par)
0044   {
0045     const int ieta = static_cast<int>(std::floor(x[0]));
0046     const double z = x[1];
0047     if (ieta < 0 || ieta > 1)
0048     {
0049       TF2::RejectPoint();
0050       return 0.;
0051     }
0052     return par[ieta + 1] + par[0] * z;
0053   }
0054
0055   //! suffixes used in histogram names for the two eta bins
0056   const char* k_eta_suffix[2] = {"negeta", "poseta"};
0057
0058   //! number of z bins of the dz vs z histograms
0059   constexpr int k_nzbins = 200;
0060
0061 }  // namespace
0062
0063 //____________________________________________________________________________..
0064 SiliconDriftQA::SiliconDriftQA(const std::string& name)
0065   : SubsysReco(name)
0066 {
0067 }
0068
0069 //____________________________________________________________________________..
0070 int SiliconDriftQA::InitRun(PHCompositeNode* /*topNode*/)
0071 {
0072   createHistos();
0073
0074   // reference histograms initialized in header file to histos in HistoManager
0075   auto* hm = QAHistManagerDef::getHistoManager();
0076   assert(hm);
0077
0078   for (int ieta = 0; ieta < 2; ieta++)
0079   {
0080     h_zsi_dz[ieta] = dynamic_cast<TH2*>(hm->getHisto(std::format("{}zsi_dz_{}", getHistoPrefix(), k_eta_suffix[ieta])));
0081   }
0082   h_dz = dynamic_cast<TH1*>(hm->getHisto(std::format("{}dz", getHistoPrefix())));
0083   h_ntracks = dynamic_cast<TH1*>(hm->getHisto(std::format("{}ntracks", getHistoPrefix())));
0084   h_driftSummary = dynamic_cast<TH1*>(hm->getHisto(std::format("{}driftSummary", getHistoPrefix())));
0085
0086   return Fun4AllReturnCodes::EVENT_OK;
0087 }
0088
0089 //____________________________________________________________________________..
0090 int SiliconDriftQA::process_event(PHCompositeNode* topNode)
0091 {
0092   auto* track_map = findNode::getClass<SvtxTrackMap>(topNode, m_trackmapname);
0093   if (!track_map)
0094   {
0095     std::cout << PHWHERE << " " << m_trackmapname << " node missing, abort." << std::endl;
0096     return Fun4AllReturnCodes::ABORTRUN;
0097   }
0098
0099   int naccepted = 0;
0100
0101   for (const auto& [track_id, track] : *track_map)
0102   {
0103     // require valid beam-crossing
0104     const auto crossing = track->get_crossing();
0105     if (crossing == SHRT_MAX)
0106     {
0107       if (Verbosity())
0108       {
0109         std::cout << PHWHERE << " invalid crossing, track ignored." << std::endl;
0110       }
0111       continue;
0112     }
0113
0114     // require both seeds
0115     const auto* si_seed = track->get_silicon_seed();
0116     const auto* tpc_seed = track->get_tpc_seed();
0117     if (!si_seed || !tpc_seed)
0118     {
0119       continue;
0120     }
0121
0122     // count clusters per subsystem
0123     unsigned int n_tpc = 0;
0124     unsigned int n_mvtx = 0;
0125     unsigned int n_intt = 0;
0126
0127     for (const auto* seed : {si_seed, tpc_seed})
0128     {
0129       for (auto it = seed->begin_cluster_keys(); it != seed->end_cluster_keys(); ++it)
0130       {
0131         switch (TrkrDefs::getTrkrId(*it))
0132         {
0133         case TrkrDefs::tpcId:
0134           ++n_tpc;
0135           break;
0136         case TrkrDefs::mvtxId:
0137           ++n_mvtx;
0138           break;
0139         case TrkrDefs::inttId:
0140           ++n_intt;
0141           break;
0142         default:
0143           break;
0144         }
0145       }
0146     }
0147
0148     // apply selection cuts
0149     if (n_tpc < m_min_nclusters_tpc)
0150     {
0151       continue;
0152     }
0153     if (n_mvtx < m_min_nclusters_mvtx)
0154     {
0155       continue;
0156     }
0157     if (n_intt < m_min_nclusters_intt)
0158     {
0159       continue;
0160     }
0161
0162     const float eta = tpc_seed->get_eta();
0163     if (std::abs(eta) > m_max_eta)
0164     {
0165       continue;
0166     }
0167
0168     const float pt = get_pt(track->get_px(), track->get_py());
0169     if (pt < m_min_pt)
0170     {
0171       continue;
0172     }
0173
0174     // get seed z positions at POCA
0175     const auto si_pos = TrackSeedHelper::get_xyz(si_seed);
0176     const auto tpc_pos = TrackSeedHelper::get_xyz(tpc_seed);
0177
0178     const float z_si = si_pos.z();
0179     const float z_tpc = tpc_pos.z();
0180
0181     // dz = (z_tpc + sign(eta)*crossing*crossing_interval*dv) - z_si
0182     const double sign_eta = (eta >= 0) ? 1.0 : -1.0;
0183     const float z_tpc_corr = z_tpc + sign_eta * crossing * m_crossing_interval * m_drift_velocity;
0184     const float dz = z_tpc_corr - z_si;
0185
0186     // fill histograms
0187     const int ieta = (eta >= 0) ? 1 : 0;
0188     h_zsi_dz[ieta]->Fill(z_si, dz);
0189     h_dz->Fill(dz);
0190
0191     ++naccepted;
0192   }
0193
0194   h_ntracks->Fill(naccepted);
0195
0196   return Fun4AllReturnCodes::EVENT_OK;
0197 }
0198
0199 //____________________________________________________________________________..
0200 int SiliconDriftQA::End(PHCompositeNode* /*topNode*/)
0201 {
0202   if (!(h_zsi_dz[0] && h_zsi_dz[1] && h_driftSummary))
0203   {
0204     std::cout << PHWHERE << " histograms not found, skipping drift velocity fit." << std::endl;
0205     return Fun4AllReturnCodes::EVENT_OK;
0206   }
0207
0208   const int nEntries = static_cast<int>(h_zsi_dz[0]->GetEntries() + h_zsi_dz[1]->GetEntries());
0209   if (Verbosity())
0210   {
0211     std::cout << Name() << "::End - fitting " << nEntries << " entries" << std::endl;
0212   }
0213
0214   // record input drift velocity even if the fit is skipped
0215   h_driftSummary->SetBinContent(3, m_drift_velocity);
0216
0217   if (nEntries < 2 * m_min_slice_entries)
0218   {
0219     std::cout << Name() << "::End - not enough entries (" << nEntries << "), skipping drift velocity fit." << std::endl;
0220     return Fun4AllReturnCodes::EVENT_OK;
0221   }
0222
0223   // build mean-dz TH2 via FitSlicesY, one eta bin at a time
0224   // x = eta bin [0,2), y = z_si (cm), content = mean dz (cm)
0225   auto* h_fit = new TH2F("h_fit_silicon", "", 2, 0, 2, k_nzbins, -m_max_z, m_max_z);
0226   h_fit->SetDirectory(nullptr);
0227
0228   for (int ieta = 0; ieta < 2; ++ieta)
0229   {
0230     auto* h2d = h_zsi_dz[ieta];
0231
0232     // fit vertical slices; require a minimum of m_min_slice_entries per slice
0233     TObjArray slices;
0234     slices.SetOwner(kTRUE);
0235     h2d->FitSlicesY(nullptr, 0, -1, m_min_slice_entries, "QNR", &slices);
0236     auto* h_mean = dynamic_cast<TH1*>(slices.At(1));
0237     if (!h_mean)
0238     {
0239       continue;
0240     }
0241
0242     for (int iz = 1; iz <= h_mean->GetNbinsX(); ++iz)
0243     {
0244       const double entries = h2d->Integral(iz, iz, 1, h2d->GetNbinsY());
0245       if (entries > 0)
0246       {
0247         h_fit->SetBinContent(ieta + 1, iz, h_mean->GetBinContent(iz));
0248       }
0249     }
0250   }
0251
0252   // 2D piecewise fit: shared slope + per-eta offset
0253   auto* fit2d = new TF2("fit2d_silicon", fit_function_2d, 0, 2, -m_max_z, m_max_z, 3);
0254   for (int i = 0; i < 3; ++i)
0255   {
0256     fit2d->SetParameter(i, 0.0);
0257   }
0258   h_fit->Fit(fit2d, "0RQ");
0259
0260   const double slope = fit2d->GetParameter(0);
0261   const double slope_err = fit2d->GetParError(0);
0262   const double off_neg = fit2d->GetParameter(1);  // ieta=0, eta<0
0263   const double off_pos = fit2d->GetParameter(2);  // ieta=1, eta>=0
0264
0265   const double dv_new = m_drift_velocity / (1.0 + slope);
0266   const double dv_err = m_drift_velocity / std::pow(1.0 + slope, 2) * slope_err;
0267   const double t0_new = (off_pos - off_neg) / (2.0 * dv_new);
0268
0269   std::cout << Name() << "::End"
0270             << " slope=" << slope
0271             << " dv_in=" << m_drift_velocity << " cm/ns"
0272             << " dv_new=" << dv_new << " +/- " << dv_err << " cm/ns"
0273             << " t0_new=" << t0_new << " ns"
0274             << std::endl;
0275
0276   // store fit results in the summary histogram
0277   h_driftSummary->SetBinContent(1, slope);
0278   h_driftSummary->SetBinContent(2, slope_err);
0279   h_driftSummary->SetBinContent(4, dv_new);
0280   h_driftSummary->SetBinContent(5, dv_err);
0281   h_driftSummary->SetBinContent(6, t0_new);
0282
0283   delete fit2d;
0284   delete h_fit;
0285
0286   return Fun4AllReturnCodes::EVENT_OK;
0287 }
0288
0289 //____________________________________________________________________________..
0290 std::string SiliconDriftQA::getHistoPrefix() const
0291 {
0292   // define prefix to all histos in HistoManager
0293   return std::string("h_") + Name() + std::string("_");
0294 }
0295
0296 //____________________________________________________________________________..
0297 void SiliconDriftQA::createHistos()
0298 {
0299   // initialize HistoManager
0300   auto* hm = QAHistManagerDef::getHistoManager();
0301   assert(hm);
0302
0303   // create and register histos in HistoManager
0304   for (int ieta = 0; ieta < 2; ieta++)
0305   {
0306     auto* h = new TH2F(std::format("{}zsi_dz_{}", getHistoPrefix(), k_eta_suffix[ieta]).c_str(),
0307                        std::format("{};z_{{silicon}} (cm);#Deltaz_{{TPC-silicon}} (cm)",
0308                                    (ieta == 0 ? "#eta_{TPC} < 0" : "#eta_{TPC} #geq 0"))
0309                            .c_str(),
0310                        k_nzbins, -m_max_z, m_max_z, 200, -m_max_dz, m_max_dz);
0311     hm->registerHisto(h);
0312   }
0313
0314   {
0315     auto* h = new TH1F(std::format("{}dz", getHistoPrefix()).c_str(),
0316                        ";#Deltaz_{TPC-silicon} (cm);tracks", 200, -m_max_dz, m_max_dz);
0317     hm->registerHisto(h);
0318   }
0319
0320   {
0321     auto* h = new TH1F(std::format("{}ntracks", getHistoPrefix()).c_str(),
0322                        ";accepted tracks per event;events", 50, -0.5, 49.5);
0323     hm->registerHisto(h);
0324   }
0325
0326   {
0327     // summary of the drift velocity fit performed in End()
0328     auto* h = new TH1F(std::format("{}driftSummary", getHistoPrefix()).c_str(),
0329                        "drift velocity fit summary", 6, 0.5, 6.5);
0330     h->GetXaxis()->SetBinLabel(1, "slope");
0331     h->GetXaxis()->SetBinLabel(2, "slope_err");
0332     h->GetXaxis()->SetBinLabel(3, "v_{in} (cm/ns)");
0333     h->GetXaxis()->SetBinLabel(4, "v_{new} (cm/ns)");
0334     h->GetXaxis()->SetBinLabel(5, "v_{new} err (cm/ns)");
0335     h->GetXaxis()->SetBinLabel(6, "t_{0} (ns)");
0336     hm->registerHisto(h);
0337   }
0338 }