QA/Tracking/VertexQA.cc

0001
0002 #include "VertexQA.h"
0003
0004 #include <fun4all/Fun4AllHistoManager.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <qautils/QAHistManagerDef.h>
0008 #include <qautils/QAUtil.h>
0009
0010 #include <globalvertex/SvtxVertex.h>
0011 #include <globalvertex/SvtxVertexMap.h>
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014
0015 #include <TH2.h>
0016 #include <cassert>
0017 //____________________________________________________________________________..
0018 VertexQA::VertexQA(const std::string &name)
0019   : SubsysReco(name)
0020 {
0021 }
0022
0023 //____________________________________________________________________________..
0024 int VertexQA::InitRun(PHCompositeNode * /*unused*/)
0025 {
0026   createHistos();
0027   return Fun4AllReturnCodes::EVENT_OK;
0028 }
0029
0030 //____________________________________________________________________________..
0031 int VertexQA::process_event(PHCompositeNode *topNode)
0032 {
0033   // auto trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0034   auto *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
0035   if (!vertexmap)
0036   {
0037     std::cout << PHWHERE << "Missing node(s), can't continue" << std::endl;
0038     return Fun4AllReturnCodes::ABORTEVENT;
0039   }
0040
0041   auto *hm = QAHistManagerDef::getHistoManager();
0042   assert(hm);
0043
0044   auto *h_nvertex = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecovertices")));
0045   auto *h_vx = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vx")));
0046   auto *h_vy = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vy")));
0047   auto *h_vz = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vz")));
0048   auto *h_vt = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vt")));
0049   auto *h_vcrossing = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vertexcrossing")));
0050   auto *h_vchi2 = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vertexchi2")));
0051   auto *h_vndof = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vertexndof")));
0052   auto *h_ntrackpervertex = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntrackspervertex")));
0053
0054   m_vertices += vertexmap->size();
0055   h_nvertex->Fill(vertexmap->size());
0056   for (const auto &[key, vertex] : *vertexmap)
0057   {
0058     if (!vertex)
0059     {
0060       continue;
0061     }
0062
0063     float vx = vertex->get_x();
0064     float vy = vertex->get_y();
0065     float vz = vertex->get_z();
0066     float vt = vertex->get_t0();
0067     float vchi2 = vertex->get_chisq();
0068     int vndof = vertex->get_ndof();
0069     int vcrossing = vertex->get_beam_crossing();
0070
0071     h_vx->Fill(vx);
0072     h_vy->Fill(vy);
0073     h_vz->Fill(vz);
0074     h_vt->Fill(vt);
0075     h_vchi2->Fill(vchi2);
0076     h_vndof->Fill(vndof);
0077     h_vcrossing->Fill(vcrossing);
0078
0079     h_ntrackpervertex->Fill(vertex->size_tracks());
0080     // loop over all tracks on vertex
0081     // for (Vertex::TrackIter iter = vertex->begin_tracks();
0082     //     iter != vertex->end_tracks();
0083     //     ++iter)
0084     //{
0085     //  SvtxTrack* track = trackmap->get(*iter);
0086
0087     //  if (!track)
0088     //  {
0089     //    continue;
0090     //  }
0091
0092     //  float px = track->get_px();
0093     //  float py = track->get_py();
0094     //  float pz = track->get_pz();
0095     //  float pt = std::sqrt(QAG4Util::square(px) + QAG4Util::square(py));
0096     //  float eta = std::atanh(pz / std::sqrt(QAG4Util::square(pt) + QAG4Util::square(pz)));
0097     //  float phi = std::atan2(py, px);
0098     //  int tcrossing = track->get_crossing();
0099     //}
0100   }
0101
0102   m_event++;
0103   return Fun4AllReturnCodes::EVENT_OK;
0104 }
0105
0106 //____________________________________________________________________________..
0107 int VertexQA::EndRun(const int runnumber)
0108 {
0109   auto *hm = QAHistManagerDef::getHistoManager();
0110   assert(hm);
0111   TH2 *h_verticesperevent = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nverticesperrun")));
0112   // NOLINTNEXTLINE(bugprone-integer-division)
0113   h_verticesperevent->Fill(runnumber, m_vertices / m_event);
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 //____________________________________________________________________________..
0118 int VertexQA::End(PHCompositeNode * /*unused*/)
0119 {
0120   return Fun4AllReturnCodes::EVENT_OK;
0121 }
0122
0123 std::string VertexQA::getHistoPrefix() const
0124 {
0125   return std::string("h_") + Name() + std::string("_");
0126 }
0127
0128 void VertexQA::createHistos()
0129 {
0130   auto *hm = QAHistManagerDef::getHistoManager();
0131   assert(hm);
0132   {
0133     auto *h = new TH1F(std::string(getHistoPrefix() + "nrecovertices").c_str(),
0134                        "Num of reco vertices per event", 200, 0, 200);
0135     h->GetXaxis()->SetTitle("nVertices");
0136     hm->registerHisto(h);
0137   }
0138   {
0139     auto *h = new TH1F(std::string(getHistoPrefix() + "vx").c_str(),
0140                        "Vertex x", 100, -0.1, 0.1);
0141     h->GetXaxis()->SetTitle("vx [cm]");
0142     hm->registerHisto(h);
0143   }
0144   {
0145     auto *h = new TH1F(std::string(getHistoPrefix() + "vy").c_str(),
0146                        "Vertex y", 100, -0.1, 0.1);
0147     h->GetXaxis()->SetTitle("vy [cm]");
0148     hm->registerHisto(h);
0149   }
0150   {
0151     auto *h = new TH1F(std::string(getHistoPrefix() + "vz").c_str(),
0152                        "Vertex z", 100, -15, 15);
0153     h->GetXaxis()->SetTitle("vz [cm]");
0154     hm->registerHisto(h);
0155   }
0156   {
0157     auto *h = new TH1F(std::string(getHistoPrefix() + "vt").c_str(),
0158                        "Vertex t", 100, -1000, 20000);
0159     h->GetXaxis()->SetTitle("vt [ns]");
0160     hm->registerHisto(h);
0161   }
0162   {
0163     auto *h = new TH1F(std::string(getHistoPrefix() + "vertexcrossing").c_str(),
0164                        "Vertex beam bunch crossing", 100, -100, 300);
0165     h->GetXaxis()->SetTitle("vertex crossing");
0166     hm->registerHisto(h);
0167   }
0168   {
0169     auto *h = new TH1F(std::string(getHistoPrefix() + "vertexchi2").c_str(),
0170                        "Vertex chi2", 100, 0, 10000);
0171     h->GetXaxis()->SetTitle("vertex #chi2");
0172     hm->registerHisto(h);
0173   }
0174   {
0175     auto *h = new TH1F(std::string(getHistoPrefix() + "vertexndof").c_str(),
0176                        "Vertex ndof", 50, 0, 50);
0177     h->GetXaxis()->SetTitle("vertex ndof");
0178     hm->registerHisto(h);
0179   }
0180   {
0181     auto *h = new TH1F(std::string(getHistoPrefix() + "ntrackspervertex").c_str(),
0182                        "Num of tracks per vertex", 20, 0, 20);
0183     h->GetXaxis()->SetTitle("nTracks per vertex");
0184     hm->registerHisto(h);
0185   }
0186
0187   {
0188     auto *h = new TH2F(std::string(getHistoPrefix() + "nverticesperrun").c_str(),
0189                        "Num reconstructed vertices per run", m_runbins, m_beginRun, m_endRun, 100, 0, 1);
0190     h->GetYaxis()->SetTitle("N_{vertives}/event");
0191     h->GetXaxis()->SetTitle("Run number");
0192     hm->registerHisto(h);
0193   }
0194 }