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