QA/Tracking/TrackFittingQA.cc

0001 #include "TrackFittingQA.h"
0002
0003 #include <trackbase_historic/SvtxAlignmentState.h>
0004 #include <trackbase_historic/SvtxTrack.h>
0005 #include <trackbase_historic/SvtxTrackMap.h>
0006 #include <trackbase_historic/SvtxTrackState.h>
0007
0008 #include <qautils/QAHistManagerDef.h>
0009 #include <qautils/QAUtil.h>
0010
0011 #include <fun4all/Fun4AllHistoManager.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/getClass.h>
0017
0018 #include <Rtypes.h>
0019 #include <TH1F.h>
0020 #include <TH2F.h>
0021 #include <TProfile.h>
0022
0023 #include <format>
0024
0025 namespace
0026 {
0027   template <typename T>
0028   inline T square (T const& t) { return t * t; }
0029
0030   template <class T>
0031   class range_adaptor
0032   {
0033    public:
0034     explicit range_adaptor(
0035         T const& begin,
0036         T const& end)
0037       : m_begin(begin)
0038       , m_end(end)
0039     {
0040     }
0041     T const& begin() { return m_begin; }
0042     T const& end() { return m_end; }
0043
0044    private:
0045     T m_begin;
0046     T m_end;
0047   };
0048 }  // namespace
0049
0050 TrackFittingQA::TrackFittingQA(const std::string& name)
0051   : SubsysReco(name)
0052 {
0053   // ...
0054 }
0055
0056 int TrackFittingQA::Init(PHCompositeNode* /*unused*/)
0057 {
0058   auto* hm = QAHistManagerDef::getHistoManager();
0059   if (!hm)
0060   {
0061     std::cout
0062         << PHWHERE << "\n"
0063         << "\tCould not get QAHistManager\n"
0064         << "\tAborting\n"
0065         << std::endl;
0066     return Fun4AllReturnCodes::ABORTRUN;
0067   }
0068
0069   // Create histograms
0070
0071   for (int charge = 0; charge < 2; ++charge)
0072   {
0073     std::string charge_str = charge ? "positively" : "negatively";
0074     EColor color = charge ? kRed : kBlue;
0075
0076     delete m_quality_hist[charge];
0077     m_quality_hist[charge] = new TH1F(
0078         std::format("h_{}_quality_{}_charged_tracks", Name(), charge_str).c_str(),
0079         std::format("quality distribution ({} charged tracks);quality;Counts", charge_str).c_str(),
0080         50, m_quality_hist_xrange.first, m_quality_hist_xrange.second);
0081     hm->registerHisto(m_quality_hist[charge]);
0082     m_quality_hist[charge]->SetMarkerColor(color);
0083     m_quality_hist[charge]->SetLineColor(color);
0084     // ...
0085
0086     delete m_p_hist[charge];
0087     m_p_hist[charge] = new TH1F(
0088         std::format("h_{}_p_{}_charged_tracks", Name(), charge_str).c_str(),
0089         std::format("p distribution ({} charged tracks);p (GeV);Counts", charge_str).c_str(),
0090         50, m_p_hist_xrange.first, m_p_hist_xrange.second);
0091     hm->registerHisto(m_p_hist[charge]);
0092     m_p_hist[charge]->SetMarkerColor(color);
0093     m_p_hist[charge]->SetLineColor(color);
0094     // ...
0095
0096     delete m_pt_hist[charge];
0097     m_pt_hist[charge] = new TH1F(
0098         std::format("h_{}_pt_{}_charged_tracks", Name(), charge_str).c_str(),
0099         std::format("pt distribution ({} charged tracks);pt (GeV);Counts", charge_str).c_str(),
0100         50, m_pt_hist_xrange.first, m_pt_hist_xrange.second);
0101     hm->registerHisto(m_pt_hist[charge]);
0102     m_pt_hist[charge]->SetMarkerColor(color);
0103     m_pt_hist[charge]->SetLineColor(color);
0104     // ...
0105
0106     delete m_pt_err_hist[charge];
0107     m_pt_err_hist[charge] = new TH2F(
0108         std::format("h_{}_pt_err_{}_charged_tracks", Name(), charge_str).c_str(),
0109         std::format("relative pt err vs pt distribution ({} charged tracks);pt (GeV);pt err / pt (unitless)", charge_str).c_str(),
0110         50, m_pt_err_hist_xrange.first, m_pt_err_hist_xrange.second,
0111         50, m_pt_err_hist_yrange.first, m_pt_err_hist_yrange.second);
0112     hm->registerHisto(m_pt_err_hist[charge]);
0113     m_pt_err_hist[charge]->SetMarkerColor(color);
0114     m_pt_err_hist[charge]->SetLineColor(color);
0115     // ...
0116
0117     delete m_pt_err_profile[charge];
0118     m_pt_err_profile[charge] = new TProfile(
0119         std::format("h_{}_pt_err_{}_charged_tracks_profile", Name(), charge_str).c_str(),
0120         std::format("relative pt err vs pt distribution ({} charged tracks);pt (GeV);pt err / pt (unitless)", charge_str).c_str(),
0121         50, m_pt_err_hist_xrange.first, m_pt_err_hist_xrange.second,
0122         m_pt_err_hist_yrange.first, m_pt_err_hist_yrange.second);
0123     hm->registerHisto(m_pt_err_profile[charge]);
0124     m_pt_err_profile[charge]->SetMarkerColor(color);
0125     m_pt_err_profile[charge]->SetLineColor(color);
0126     // ...
0127
0128     delete m_eta_hist[charge];
0129     m_eta_hist[charge] = new TH1F(
0130         std::format("h_{}_eta_{}_charged_tracks", Name(), charge_str).c_str(),
0131         std::format("eta distribution ({} charged tracks);eta;Counts", charge_str).c_str(),
0132         50, -1.2, 1.2);
0133     hm->registerHisto(m_eta_hist[charge]);
0134     m_eta_hist[charge]->SetMarkerColor(color);
0135     m_eta_hist[charge]->SetLineColor(color);
0136     // ...
0137
0138     delete m_phi_eta_hist[charge];
0139     m_phi_eta_hist[charge] = new TH2F(
0140         std::format("h_{}_phi_eta_{}_charged_tracks", Name(), charge_str).c_str(),
0141         std::format("phi-eta distribution ({} charged tracks);phi (radians);eta", charge_str).c_str(),
0142         50, -3.2, 3.2,
0143         50, -1.2, 1.2);
0144     hm->registerHisto(m_phi_eta_hist[charge]);
0145     // ...
0146
0147     delete m_mvtx_states_hist[charge];
0148     m_mvtx_states_hist[charge] = new TH1F(
0149         std::format("h_{}_mvtx_states_{}_charged_tracks", Name(), charge_str).c_str(),
0150         std::format("mvtx states distribution ({} charged tracks);number of mvtx states;tracks with that number of mvtx states", charge_str).c_str(),
0151         6, -0.5, 5.5);
0152     hm->registerHisto(m_mvtx_states_hist[charge]);
0153     m_mvtx_states_hist[charge]->SetMarkerColor(color);
0154     m_mvtx_states_hist[charge]->SetLineColor(color);
0155     // ...
0156
0157     delete m_intt_states_hist[charge];
0158     m_intt_states_hist[charge] = new TH1F(
0159         std::format("h_{}_intt_states_{}_charged_tracks", Name(), charge_str).c_str(),
0160         std::format("intt states distribution ({} charged tracks);number of intt states;tracks with that number of intt states", charge_str).c_str(),
0161         5, -0.5, 4.5);
0162     hm->registerHisto(m_intt_states_hist[charge]);
0163     m_intt_states_hist[charge]->SetMarkerColor(color);
0164     m_intt_states_hist[charge]->SetLineColor(color);
0165     // ...
0166
0167     delete m_tpc_states_hist[charge];
0168     m_tpc_states_hist[charge] = new TH1F(
0169         std::format("h_{}_tpc_states_{}_charged_tracks", Name(), charge_str).c_str(),
0170         std::format("tpc states distribution ({} charged tracks);number of tpc states;tracks with that number of tpc states", charge_str).c_str(),
0171         61, -0.5, 60.5);
0172     hm->registerHisto(m_tpc_states_hist[charge]);
0173     m_tpc_states_hist[charge]->SetMarkerColor(color);
0174     m_tpc_states_hist[charge]->SetLineColor(color);
0175     // ...
0176
0177     delete m_tpot_states_hist[charge];
0178     m_tpot_states_hist[charge] = new TH1F(
0179         std::format("h_{}_tpot_states_{}_charged_tracks", Name(), charge_str).c_str(),
0180         std::format("tpot states distribution ({} charged tracks);number of tpot states;tracks with that number of tpot states", charge_str).c_str(),
0181         4, -0.5, 3.5);
0182     hm->registerHisto(m_tpot_states_hist[charge]);
0183     m_tpot_states_hist[charge]->SetMarkerColor(color);
0184     m_tpot_states_hist[charge]->SetLineColor(color);
0185     // ...
0186   }
0187
0188   return Fun4AllReturnCodes::EVENT_OK;
0189 }
0190
0191 int TrackFittingQA::InitRun(
0192     PHCompositeNode* top_node)
0193 {
0194   // F4A will not actually ABORTRUN unless that return code is issued here
0195   auto* track_map = findNode::getClass<SvtxTrackMap>(top_node, m_track_map_node_name);
0196   if (!track_map)
0197   {
0198     std::cout
0199         << PHWHERE << "\n"
0200         << "\tCould not get track map:\n"
0201         << "\t\"" << m_track_map_node_name << "\"\n"
0202         << "\tAborting\n"
0203         << std::endl;
0204     return Fun4AllReturnCodes::ABORTRUN;
0205   }
0206
0207   return Fun4AllReturnCodes::EVENT_OK;
0208 }
0209
0210 int TrackFittingQA::process_event(
0211     PHCompositeNode* top_node)
0212 {
0213   auto* track_map = findNode::getClass<SvtxTrackMap>(top_node, m_track_map_node_name);
0214   if (!track_map)
0215   {
0216     std::cout
0217         << PHWHERE << "\n"
0218         << "\tCould not get track map:\n"
0219         << "\t\"" << m_track_map_node_name << "\"\n"
0220         << "\tAborting\n"
0221         << std::endl;
0222     return Fun4AllReturnCodes::ABORTRUN;
0223   }
0224
0225   for (auto const& [idkey, track] : *track_map)
0226   {
0227     // count states
0228     std::map<TrkrDefs::TrkrId, int> counters = {
0229         {TrkrDefs::mvtxId, 0},
0230         {TrkrDefs::inttId, 0},
0231         {TrkrDefs::tpcId, 0},
0232         {TrkrDefs::micromegasId, 0},
0233     };
0234
0235     for (auto const& [path_length, state] : range_adaptor(track->begin_states(), track->end_states()))
0236     {
0237       // There is an additional state representing the vertex at the beginning of the map,
0238       // but getTrkrId will return 0 for its corresponding cluster
0239       // Identify it as having path_length identically equal to 0
0240       if (path_length == 0) { continue; }
0241
0242       auto trkr_id = static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(state->get_cluskey()));
0243       auto itr = counters.find(trkr_id);
0244       if (itr == counters.end()) { continue; }
0245       ++itr->second;
0246     }
0247
0248     // Cuts
0249     if ( track->get_quality() < m_min_quality ) { continue; }
0250     if ( m_max_quality < track->get_quality() ) { continue; }
0251     if ( track->get_p() < m_min_p ) { continue; }
0252     if ( track->get_pt() < m_min_pt ) { continue; }
0253     if ( m_max_abs_eta < std::abs(track->get_eta()) ) { continue; }
0254     if ( counters[TrkrDefs::inttId] < m_min_intt_states ) { continue; }
0255     if ( counters[TrkrDefs::mvtxId] < m_min_mvtx_states ) { continue; }
0256     if ( counters[TrkrDefs::tpcId] < m_min_tpc_states ) { continue; }
0257     if ( counters[TrkrDefs::micromegasId] < m_min_tpot_states ) { continue; }
0258     if ( track->get_crossing() < m_min_crossing ) { continue; }
0259     if ( m_max_crossing < track->get_crossing() ) { continue; }
0260
0261     // Fill histograms if passing
0262     int charge = (0 < track->get_charge()) ? 1 : 0;
0263
0264     // Immitated from TrackResiduals.cc:1896-1899 (as of commit b75c9ef557d096f7d6ec9bbf2082d4291eddc4a1)
0265     float px = track->get_px();
0266     float py = track->get_py();
0267     float CVxx = track->get_error(3, 3);
0268     float CVxy = track->get_error(3, 4);
0269     float CVyy = track->get_error(4, 4);
0270     float pt_err = std::sqrt((CVxx * square(px) + 2 * CVxy * px * py + CVyy * square(py)) / (square(px) + square(py)));
0271
0272     m_quality_hist[charge]->Fill(track->get_quality());
0273     m_p_hist[charge]->Fill(track->get_p());
0274     m_pt_hist[charge]->Fill(track->get_pt());
0275     m_pt_err_hist[charge]->Fill(track->get_pt(), pt_err);
0276     m_pt_err_profile[charge]->Fill(track->get_pt(), pt_err);
0277     m_eta_hist[charge]->Fill(track->get_eta());
0278     m_phi_eta_hist[charge]->Fill(track->get_phi(), track->get_eta());
0279
0280     m_mvtx_states_hist[charge]->Fill(counters[TrkrDefs::mvtxId]);
0281     m_intt_states_hist[charge]->Fill(counters[TrkrDefs::inttId]);
0282     m_tpc_states_hist[charge]->Fill(counters[TrkrDefs::tpcId]);
0283     m_tpot_states_hist[charge]->Fill(counters[TrkrDefs::micromegasId]);
0284   }
0285
0286   return Fun4AllReturnCodes::EVENT_OK;
0287 }
0288
0289 int TrackFittingQA::ResetEvent(
0290     PHCompositeNode* /*unused*/)
0291 {
0292   return Fun4AllReturnCodes::EVENT_OK;
0293 }
0294
0295 int TrackFittingQA::EndRun(
0296     int const /*unused*/)
0297 {
0298   return Fun4AllReturnCodes::EVENT_OK;
0299 }
0300
0301 int TrackFittingQA::End(
0302     PHCompositeNode* /*unused*/)
0303 {
0304   return Fun4AllReturnCodes::EVENT_OK;
0305 }
0306
0307 int TrackFittingQA::Reset(
0308     PHCompositeNode* /*unused*/)
0309 {
0310   return Fun4AllReturnCodes::EVENT_OK;
0311 }
0312
0313 void TrackFittingQA::Print(
0314     std::string const& /*unused*/) const
0315 {
0316   // ...
0317 }