Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-02 08:13:41 

0001 #include "StateClusterResidualsQA.h"
0002 
0003 #include <trackbase/InttDefs.h>
0004 #include <trackbase/MvtxDefs.h>
0005 #include <trackbase/TpcDefs.h>
0006 #include <trackbase/TrkrCluster.h>
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/ActsGeometry.h>
0009 
0010 #include <trackbase_historic/SvtxAlignmentState.h>
0011 #include <trackbase_historic/SvtxTrack.h>
0012 #include <trackbase_historic/SvtxTrackMap.h>
0013 #include <trackbase_historic/SvtxTrackState.h>
0014 
0015 #include <qautils/QAHistManagerDef.h>
0016 #include <qautils/QAUtil.h>
0017 
0018 #include <fun4all/Fun4AllHistoManager.h>
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 #include <fun4all/SubsysReco.h>
0021 
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/getClass.h>
0024 
0025 #include <Rtypes.h>
0026 #include <TH1F.h>
0027 #include <TH2F.h>
0028 
0029 #include <format>
0030 
0031 namespace
0032 {
0033   template <typename T>
0034   inline T square (T const& t) { return t * t; }
0035 
0036   template <class T>
0037   class range_adaptor
0038   {
0039    public:
0040     explicit range_adaptor(
0041         T const& begin,
0042         T const& end)
0043       : m_begin(begin)
0044       , m_end(end)
0045     {
0046     }
0047     T const& begin() { return m_begin; }
0048     T const& end() { return m_end; }
0049 
0050    private:
0051     T m_begin;
0052     T m_end;
0053   };
0054 }  // namespace
0055 
0056 StateClusterResidualsQA::StateClusterResidualsQA(const std::string& name)
0057   : SubsysReco(name)
0058 {
0059 }
0060 
0061 int StateClusterResidualsQA::InitRun(
0062     PHCompositeNode* top_node)
0063 {
0064   createHistos();
0065 
0066   // F4A will not actually ABORTRUN unless that return code is issued here
0067   auto* track_map = findNode::getClass<SvtxTrackMap>(top_node, m_track_map_node_name);
0068   if (!track_map)
0069   {
0070     std::cout
0071         << PHWHERE << "\n"
0072         << "\tCould not get track map:\n"
0073         << "\t\"" << m_track_map_node_name << "\"\n"
0074         << "\tAborting\n"
0075         << std::endl;
0076     return Fun4AllReturnCodes::ABORTRUN;
0077   }
0078 
0079   auto* cluster_map = findNode::getClass<TrkrClusterContainer>(top_node, m_clusterContainerName);
0080   if (!cluster_map)
0081   {
0082     std::cout
0083         << PHWHERE << "\n"
0084         << "\tCould not get cluster map:\n"
0085         << "\t\"" << m_clusterContainerName << "\"\n"
0086         << "\tAborting\n"
0087         << std::endl;
0088     return Fun4AllReturnCodes::ABORTRUN;
0089   }
0090   
0091   auto *geometry = findNode::getClass<ActsGeometry>(top_node, "ActsGeometry");
0092   if (!geometry)
0093   {
0094     std::cout
0095         << PHWHERE << "\n"
0096         << "\tCould not get ActsGeometry:\n"
0097         << "\t\"" << "ActsGeometry" << "\"\n"
0098         << "\tAborting\n"
0099         << std::endl;
0100     return Fun4AllReturnCodes::ABORTRUN;
0101   }
0102 
0103   auto* hm = QAHistManagerDef::getHistoManager();
0104   if (!hm)
0105   {
0106     std::cout
0107         << PHWHERE << "\n"
0108         << "\tCould not get QAHistManager\n"
0109         << "\tAborting\n"
0110         << std::endl;
0111     return Fun4AllReturnCodes::ABORTRUN;
0112   }
0113 
0114   for (const auto& cfg : m_pending)
0115   {
0116     if (m_use_local_coords)
0117     {
0118       m_histograms_x.push_back(dynamic_cast<TH1 *>(hm->getHisto(std::string(cfg.name + "_local_rphi"))));
0119       m_histograms_y.push_back(dynamic_cast<TH1 *>(hm->getHisto(std::string(cfg.name + "_local_z"))));
0120       m_histograms_layer_x.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_local_layer_rphi"))));
0121       m_histograms_layer_y.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_local_layer_z"))));
0122       m_histograms_phi_x.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_local_phi_rphi"))));
0123       m_histograms_phi_y.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_local_phi_z"))));
0124       m_histograms_eta_x.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_local_eta_rphi"))));
0125       m_histograms_eta_y.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_local_eta_z"))));
0126     }
0127     else
0128     {
0129       m_histograms_x.push_back(dynamic_cast<TH1 *>(hm->getHisto(std::string(cfg.name + "_x"))));
0130       m_histograms_y.push_back(dynamic_cast<TH1 *>(hm->getHisto(std::string(cfg.name + "_y"))));
0131       m_histograms_z.push_back(dynamic_cast<TH1 *>(hm->getHisto(std::string(cfg.name + "_z"))));
0132       m_histograms_layer_x.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_layer_x"))));
0133       m_histograms_layer_y.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_layer_y"))));
0134       m_histograms_layer_z.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_layer_z"))));
0135       m_histograms_phi_x.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_phi_x"))));
0136       m_histograms_phi_y.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_phi_y"))));
0137       m_histograms_phi_z.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_phi_z"))));
0138       m_histograms_eta_x.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_eta_x"))));
0139       m_histograms_eta_y.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_eta_y"))));
0140       m_histograms_eta_z.push_back(dynamic_cast<TH2 *>(hm->getHisto(std::string(cfg.name + "_eta_z"))));
0141     }
0142   }
0143   
0144   return Fun4AllReturnCodes::EVENT_OK;
0145 }
0146 
0147 int StateClusterResidualsQA::process_event(PHCompositeNode* top_node)
0148 {
0149   auto* track_map = findNode::getClass<SvtxTrackMap>(top_node, m_track_map_node_name);
0150   auto *cluster_map = findNode::getClass<TrkrClusterContainer>(top_node, m_clusterContainerName);
0151   auto *geometry = findNode::getClass<ActsGeometry>(top_node, "ActsGeometry");
0152 
0153   for (auto const& [idkey, track] : *track_map)
0154   {
0155     if (!track)
0156     {
0157       continue;
0158     }
0159 
0160     // count states
0161     std::map<TrkrDefs::TrkrId, int> counters = {
0162         {TrkrDefs::mvtxId, 0},
0163         {TrkrDefs::inttId, 0},
0164         {TrkrDefs::tpcId, 0},
0165         {TrkrDefs::micromegasId, 0},
0166     };
0167 
0168     for (auto const& [path_length, state] : range_adaptor(track->begin_states(), track->end_states()))
0169     {
0170       // There is an additional state representing the vertex at the beginning of the map,
0171       // but getTrkrId will return 0 for its corresponding cluster
0172       // Identify it as having path_length identically equal to 0
0173       if (path_length == 0) { continue; }
0174 
0175       auto trkr_id = static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(state->get_cluskey()));
0176       auto itr = counters.find(trkr_id);
0177       if (itr == counters.end()) { continue; }
0178       ++itr->second;
0179     }
0180 
0181     float track_eta = track->get_eta();
0182     float track_phi = track->get_phi();
0183     float track_pt = track->get_pt();
0184     int h = 0;
0185     for (const auto& cfg : m_pending)
0186     {
0187       if (cfg.charge != 0)
0188       {
0189         if ((cfg.charge < 0) && track->get_positive_charge())
0190         {
0191           continue;
0192         }
0193         if ((cfg.charge > 0) && !(track->get_positive_charge()))
0194         {
0195           continue;
0196         }
0197       }
0198       if (cfg.min_mvtx_clusters <= counters[TrkrDefs::mvtxId] && cfg.max_mvtx_clusters >= counters[TrkrDefs::mvtxId]
0199           && cfg.min_intt_clusters <= counters[TrkrDefs::inttId] && cfg.max_intt_clusters >= counters[TrkrDefs::inttId]
0200           && cfg.min_tpc_clusters <= counters[TrkrDefs::tpcId] && cfg.max_tpc_clusters >= counters[TrkrDefs::tpcId]
0201           && cfg.phi_min <= track_phi && cfg.phi_max >= track_phi
0202           && cfg.eta_min <= track_eta && cfg.eta_max >= track_eta
0203           && cfg.pt_min <= track_pt && cfg.pt_max >= track_pt)
0204       {
0205         for (auto const& [path_length, state] : range_adaptor(track->begin_states(), track->end_states()))
0206         {
0207           if (path_length == 0) { continue; }
0208     
0209           auto *cluster = cluster_map->findCluster(state->get_cluskey());
0210           if (!cluster) 
0211           {
0212             continue;
0213           }          
0214 
0215           float state_x; 
0216           float state_y;
0217           float state_z;
0218           float cluster_x;
0219           float cluster_y;
0220           float cluster_z;
0221           if (m_use_local_coords == true)
0222           {
0223             state_x = state->get_localX();
0224             state_y = state->get_localY();
0225             Acts::Vector2 loc = geometry->getLocalCoords(state->get_cluskey(), cluster);
0226             cluster_x = loc.x();
0227             cluster_y = loc.y();
0228             m_histograms_x[h]->Fill(state_x - cluster_x);
0229             m_histograms_y[h]->Fill(state_y - cluster_y);
0230             m_histograms_layer_x[h]->Fill(TrkrDefs::getLayer(state->get_cluskey()), state_x - cluster_x);
0231             m_histograms_layer_y[h]->Fill(TrkrDefs::getLayer(state->get_cluskey()), state_y - cluster_y);
0232             m_histograms_phi_x[h]->Fill(state->get_phi(), state_x - cluster_x);
0233             m_histograms_phi_y[h]->Fill(state->get_phi(), state_y - cluster_y);
0234             m_histograms_eta_x[h]->Fill(state->get_eta(), state_x - cluster_x);
0235             m_histograms_eta_y[h]->Fill(state->get_eta(), state_y - cluster_y);
0236           }
0237           else
0238           {
0239             state_x = state->get_x();
0240             state_y = state->get_y();
0241             state_z = state->get_z();
0242             Acts::Vector3 glob = geometry->getGlobalPosition(state->get_cluskey(), cluster);
0243             cluster_x = glob.x();
0244             cluster_y = glob.y();
0245             cluster_z = glob.z();
0246             m_histograms_x[h]->Fill(state_x - cluster_x);
0247             m_histograms_y[h]->Fill(state_y - cluster_y);
0248             m_histograms_z[h]->Fill(state_z - cluster_z);
0249             m_histograms_layer_x[h]->Fill(TrkrDefs::getLayer(state->get_cluskey()), state_x - cluster_x);
0250             m_histograms_layer_y[h]->Fill(TrkrDefs::getLayer(state->get_cluskey()), state_y - cluster_y);
0251             m_histograms_layer_z[h]->Fill(TrkrDefs::getLayer(state->get_cluskey()), state_z - cluster_z);
0252             m_histograms_phi_x[h]->Fill(state->get_phi(), state_x - cluster_x);
0253             m_histograms_phi_y[h]->Fill(state->get_phi(), state_y - cluster_y);
0254             m_histograms_phi_z[h]->Fill(state->get_phi(), state_z - cluster_z);
0255             m_histograms_eta_x[h]->Fill(state->get_eta(), state_x - cluster_x);
0256             m_histograms_eta_y[h]->Fill(state->get_eta(), state_y - cluster_y);
0257             m_histograms_eta_z[h]->Fill(state->get_eta(), state_z - cluster_z);
0258           }
0259         }
0260       } 
0261       ++h;
0262     }
0263   }
0264 
0265   return Fun4AllReturnCodes::EVENT_OK;
0266 }
0267 
0268 void StateClusterResidualsQA::createHistos()
0269 {
0270   auto *hm = QAHistManagerDef::getHistoManager();
0271   assert(hm);
0272 
0273   for (const auto& cfg : m_pending)
0274   {
0275     if (m_use_local_coords)
0276     {
0277       TH1F* h_new_x = new TH1F(
0278         (cfg.name + "_local_rphi").c_str(),
0279         ";State-Cluster Local r#phi Residual [cm];Entries",
0280         m_nBins, cfg.rphi_local_lower, cfg.rphi_local_upper);
0281       h_new_x->SetMarkerColor(kBlue);
0282       h_new_x->SetLineColor(kBlue);
0283       hm->registerHisto(h_new_x);
0284       TH1F* h_new_y = new TH1F(
0285         (cfg.name + "_local_z").c_str(),
0286         ";State-Cluster Local Z Residual [cm];Entries",
0287         m_nBins, cfg.z_local_lower, cfg.z_local_upper);
0288       h_new_y->SetMarkerColor(kBlue);
0289       h_new_y->SetLineColor(kBlue);
0290       hm->registerHisto(h_new_y);
0291       TH2F* h_new_layer_x = new TH2F(
0292         (cfg.name + "_local_layer_rphi").c_str(),
0293         ";Layer Number;State-Cluster Local r#phi Residual [cm]",
0294         60, -0.5, 59.5, m_nBins, cfg.rphi_local_lower, cfg.rphi_local_upper);
0295       hm->registerHisto(h_new_layer_x);
0296       TH2F* h_new_layer_y = new TH2F(
0297         (cfg.name + "_local_layer_z").c_str(),
0298         ";Layer Number;State-Cluster Local Z Residual [cm]",
0299         60, -0.5, 59.5, m_nBins, cfg.z_local_lower, cfg.z_local_upper);
0300       hm->registerHisto(h_new_layer_y);
0301       TH2F* h_new_phi_x = new TH2F(
0302         (cfg.name + "_local_phi_rphi").c_str(),
0303         ";#phi [rad];State-Cluster Local r#phi Residual [cm]",
0304         50, -3.2, 3.2, m_nBins, cfg.rphi_local_lower, cfg.rphi_local_upper);
0305       hm->registerHisto(h_new_phi_x);
0306       TH2F* h_new_phi_y = new TH2F(
0307         (cfg.name + "_local_phi_z").c_str(),
0308         ";#phi [rad];State-Cluster Local Z Residual [cm]",
0309         50, -3.2, 3.2, m_nBins, cfg.z_local_lower, cfg.z_local_upper);
0310       hm->registerHisto(h_new_phi_y);
0311       TH2F* h_new_eta_x = new TH2F(
0312         (cfg.name + "_local_eta_rphi").c_str(),
0313         ";#eta;State-Cluster Local r#phi Residual [cm]",
0314         50, -1.1, 1.1, m_nBins, cfg.rphi_local_lower, cfg.rphi_local_upper);
0315       hm->registerHisto(h_new_eta_x);
0316       TH2F* h_new_eta_y = new TH2F(
0317         (cfg.name + "_local_eta_z").c_str(),
0318         ";#eta;State-Cluster Local Z Residual [cm]",
0319         50, -1.1, 1.1, m_nBins, cfg.z_local_lower, cfg.z_local_upper);
0320       hm->registerHisto(h_new_eta_y);
0321     }
0322     else
0323     {
0324       TH1F* h_new_x = new TH1F(
0325         (cfg.name + "_x").c_str(),
0326         ";State-Cluster X Residual [cm];Entries",
0327         m_nBins, cfg.x_lower, cfg.x_upper);
0328       h_new_x->SetMarkerColor(kBlue);
0329       h_new_x->SetLineColor(kBlue);
0330       hm->registerHisto(h_new_x);
0331       TH1F* h_new_y = new TH1F(
0332         (cfg.name + "_y").c_str(),
0333         ";State-Cluster Y Residual [cm];Entries",
0334         m_nBins, cfg.y_lower, cfg.y_upper);
0335       h_new_y->SetMarkerColor(kBlue);
0336       h_new_y->SetLineColor(kBlue);
0337       hm->registerHisto(h_new_y);
0338       TH1F* h_new_z = new TH1F(
0339         (cfg.name + "_z").c_str(),
0340         ";State-Cluster Z Residual [cm];Entries",
0341         m_nBins, cfg.z_lower, cfg.z_upper);
0342       h_new_z->SetMarkerColor(kBlue);
0343       h_new_z->SetLineColor(kBlue);
0344       hm->registerHisto(h_new_z);
0345       TH2F* h_new_layer_x = new TH2F(
0346         (cfg.name + "_layer_x").c_str(),
0347         ";Layer Number;State-Cluster Local X Residual [cm]",
0348         60, -0.5, 59.5, m_nBins, cfg.x_lower, cfg.x_upper);
0349       hm->registerHisto(h_new_layer_x);
0350       TH2F* h_new_layer_y = new TH2F(
0351         (cfg.name + "_layer_y").c_str(),
0352         ";Layer Number;State-Cluster Local Y Residual [cm]",
0353         60, -0.5, 59.5, m_nBins, cfg.y_lower, cfg.y_upper);
0354       hm->registerHisto(h_new_layer_y);
0355       TH2F* h_new_layer_z = new TH2F(
0356         (cfg.name + "_layer_z").c_str(),
0357         ";Layer Number;State-Cluster Local Z Residual [cm]",
0358         60, -0.5, 59.5, m_nBins, cfg.z_lower, cfg.z_upper);
0359       hm->registerHisto(h_new_layer_z);
0360       TH2F* h_new_phi_x = new TH2F(
0361         (cfg.name + "_phi_x").c_str(),
0362         ";#phi [rad];State-Cluster Local X Residual [cm]",
0363         50, -3.2, 3.2, m_nBins, cfg.x_lower, cfg.x_upper);
0364       hm->registerHisto(h_new_phi_x);
0365       TH2F* h_new_phi_y = new TH2F(
0366         (cfg.name + "_phi_y").c_str(),
0367         ";#phi [rad];State-Cluster Local Y Residual [cm]",
0368         50, -3.2, 3.2, m_nBins, cfg.y_lower, cfg.y_upper);
0369       hm->registerHisto(h_new_phi_y);
0370       TH2F* h_new_phi_z = new TH2F(
0371         (cfg.name + "_phi_z").c_str(),
0372         ";#phi [rad];State-Cluster Local Z Residual [cm]",
0373         50, -3.2, 3.2, m_nBins, cfg.z_lower, cfg.z_upper);
0374       hm->registerHisto(h_new_phi_z);
0375       TH2F* h_new_eta_x = new TH2F(
0376         (cfg.name + "_eta_x").c_str(),
0377         ";#eta;State-Cluster Local X Residual [cm]",
0378         50, -1.1, 1.1, m_nBins, cfg.x_lower, cfg.x_upper);
0379       hm->registerHisto(h_new_eta_x);
0380       TH2F* h_new_eta_y = new TH2F(
0381         (cfg.name + "_eta_y").c_str(),
0382         ";#eta;State-Cluster Local Y Residual [cm]",
0383         50, -1.1, 1.1, m_nBins, cfg.y_lower, cfg.y_upper);
0384       hm->registerHisto(h_new_eta_y);
0385       TH2F* h_new_eta_z = new TH2F(
0386         (cfg.name + "_eta_z").c_str(),
0387         ";#eta;State-Cluster Local Z Residual [cm]",
0388         50, -1.1, 1.1, m_nBins, cfg.z_lower, cfg.z_upper);
0389       hm->registerHisto(h_new_eta_z);
0390     }
0391   }
0392 }
0393 
0394 int StateClusterResidualsQA::EndRun(const int /*unused*/)
0395 {
0396   auto *hm = QAHistManagerDef::getHistoManager();
0397   assert(hm);
0398 
0399   return Fun4AllReturnCodes::EVENT_OK;
0400 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!