File indexing completed on 2025-08-06 08:18:49
0001
0002 #include "CosmicTrackQA.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 <phool/PHCompositeNode.h>
0011 #include <phool/getClass.h>
0012 #include <trackbase/TrackFitUtils.h>
0013 #include <trackbase/TrkrClusterContainer.h>
0014 #include <trackbase_historic/SvtxTrack.h>
0015 #include <trackbase_historic/SvtxTrackMap.h>
0016
0017 #include <TH2.h>
0018 #include <trackbase/ActsGeometry.h>
0019
0020 CosmicTrackQA::CosmicTrackQA(const std::string &name)
0021 : SubsysReco(name)
0022 {
0023 }
0024
0025
0026 int CosmicTrackQA::InitRun(PHCompositeNode * )
0027 {
0028 createHistos();
0029 return Fun4AllReturnCodes::EVENT_OK;
0030 }
0031
0032
0033 int CosmicTrackQA::process_event(PHCompositeNode *topNode)
0034 {
0035 auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0036 auto *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0037
0038 auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0039 if (!trackmap or !clustermap or !geometry)
0040 {
0041 std::cout << PHWHERE << "Missing node(s), can't continue" << std::endl;
0042 return Fun4AllReturnCodes::ABORTEVENT;
0043 }
0044
0045 auto *hm = QAHistManagerDef::getHistoManager();
0046 assert(hm);
0047
0048 auto *h_ntrack = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks")));
0049 auto *h_nmaps = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nmaps")));
0050 auto *h_nintt = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nintt")));
0051 auto *h_ntpc = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc")));
0052 auto *h_nmms = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpot")));
0053 auto *h_lxresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "lxlineresiduals")));
0054 auto *h_lzresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "lzlineresiduals")));
0055 auto *h_gzresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "gzlineresiduals")));
0056 auto *h_gyresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "gylineresiduals")));
0057 auto *h_gxresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "gxlineresiduals")));
0058
0059 auto *h_lxfitresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "lxfitresiduals")));
0060 auto *h_lzfitresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "lzfitresiduals")));
0061 auto *h_gzfitresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "gzfitresiduals")));
0062 auto *h_gyfitresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "gyfitresiduals")));
0063 auto *h_gxfitresid = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "gxfitresiduals")));
0064
0065 m_tracks += trackmap->size();
0066 for (const auto &[key, track] : *trackmap)
0067 {
0068 if (!track)
0069 {
0070 continue;
0071 }
0072
0073 auto ckeys = get_cluster_keys(track);
0074 std::vector<Acts::Vector3> cluspos;
0075 TrackFitUtils::getTrackletClusters(geometry, clustermap, cluspos, ckeys);
0076
0077 auto lineFitParams = lineFitClusters(cluspos);
0078
0079
0080 if (std::fabs(std::get<1>(lineFitParams)) > 20)
0081 {
0082 continue;
0083 }
0084 float px = track->get_px();
0085 float py = track->get_py();
0086 float pz = track->get_pz();
0087 float pt = std::sqrt(QAG4Util::square(px) + QAG4Util::square(py));
0088 float eta = std::atanh(pz / std::sqrt(QAG4Util::square(pt) + QAG4Util::square(pz)));
0089 float phi = std::atan2(py, px);
0090
0091 int nmaps = 0;
0092 int nintt = 0;
0093 int ntpc = 0;
0094 int nmms = 0;
0095
0096 for (auto &ckey : ckeys)
0097 {
0098 switch (TrkrDefs::getTrkrId(ckey))
0099 {
0100 case TrkrDefs::mvtxId:
0101 nmaps++;
0102 break;
0103 case TrkrDefs::inttId:
0104 nintt++;
0105 break;
0106 case TrkrDefs::tpcId:
0107 ntpc++;
0108 break;
0109 case TrkrDefs::micromegasId:
0110 nmms++;
0111 break;
0112 default:
0113 std::cout << "unknown cluster key" << std::endl;
0114 break;
0115 }
0116 }
0117 if (nmaps == 0)
0118 {
0119 continue;
0120 }
0121
0122 h_ntrack->Fill(phi, eta);
0123
0124 int i = 0;
0125 for (auto &ckey : ckeys)
0126 {
0127 auto &glob = cluspos[i];
0128 i++;
0129 auto *cluster = clustermap->findCluster(ckey);
0130
0131 auto intersection = TrackFitUtils::surface_3Dline_intersection(ckey, cluster, geometry,
0132 std::get<0>(lineFitParams), std::get<1>(lineFitParams), std::get<2>(lineFitParams), std::get<3>(lineFitParams));
0133
0134 auto surf = geometry->maps().getSurface(ckey, cluster);
0135
0136 Acts::Vector3 surfnorm = surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(0, 0, 0), Acts::Vector3(0, 0, 0));
0137 float statelx = std::numeric_limits<float>::quiet_NaN();
0138 float statelz = std::numeric_limits<float>::quiet_NaN();
0139 if (!std::isnan(intersection.x()))
0140 {
0141 auto locstateres = surf->globalToLocal(geometry->geometry().getGeoContext(),
0142 intersection * Acts::UnitConstants::cm,
0143 surfnorm);
0144 if (locstateres.ok())
0145 {
0146 Acts::Vector2 loc = locstateres.value() / Acts::UnitConstants::cm;
0147 statelx = loc(0);
0148 statelz = loc(1);
0149 }
0150 else
0151 {
0152 Acts::Vector3 loc = surf->transform(geometry->geometry().getGeoContext()).inverse() * (intersection * Acts::UnitConstants::cm);
0153 loc /= Acts::UnitConstants::cm;
0154 statelx = loc(0);
0155 statelz = loc(1);
0156 }
0157 }
0158 unsigned int layer = TrkrDefs::getLayer(ckey);
0159 auto localclus = geometry->getLocalCoords(ckey, cluster);
0160
0161 h_lxresid->Fill(layer, localclus.x() - statelx);
0162 h_lzresid->Fill(layer, localclus.y() - statelz);
0163 h_gxresid->Fill(layer, glob.x() - intersection.x());
0164 h_gyresid->Fill(layer, glob.y() - intersection.y());
0165 h_gzresid->Fill(layer, glob.z() - intersection.z());
0166 }
0167
0168 for (auto &&iter = track->begin_states(); iter != track->end_states(); ++iter)
0169 {
0170 const auto &[pathlength, state] = *iter;
0171 if (pathlength == 0)
0172 {
0173 continue;
0174 }
0175 Acts::Vector3 stateglob(state->get_x(), state->get_y(), state->get_z());
0176 auto cluskey = state->get_cluskey();
0177 auto *cluster = clustermap->findCluster(cluskey);
0178 auto clusglob = geometry->getGlobalPosition(cluskey, cluster);
0179 auto clusloc = geometry->getLocalCoords(cluskey, cluster);
0180 auto surf = geometry->maps().getSurface(cluskey, cluster);
0181 auto stateloc = surf->globalToLocal(geometry->geometry().getGeoContext(),
0182 stateglob * Acts::UnitConstants::cm,
0183 surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(0, 0, 0), Acts::Vector3(0, 0, 0)));
0184 float statelx = NAN;
0185 float statelz = NAN;
0186 if (stateloc.ok())
0187 {
0188 Acts::Vector2 loc = stateloc.value() / Acts::UnitConstants::cm;
0189 statelx = loc(0);
0190 statelz = loc(1);
0191 }
0192 else
0193 {
0194 Acts::Vector3 loc = surf->transform(geometry->geometry().getGeoContext()).inverse() * (stateglob * Acts::UnitConstants::cm);
0195 loc /= Acts::UnitConstants::cm;
0196 statelx = loc(0);
0197 statelz = loc(1);
0198 }
0199 h_lxfitresid->Fill(TrkrDefs::getLayer(cluskey), clusloc.x() - statelx);
0200 h_lzfitresid->Fill(TrkrDefs::getLayer(cluskey), clusloc.y() - statelz);
0201 h_gxfitresid->Fill(TrkrDefs::getLayer(cluskey), clusglob.x() - stateglob.x());
0202 h_gyfitresid->Fill(TrkrDefs::getLayer(cluskey), clusglob.y() - stateglob.y());
0203 h_gzfitresid->Fill(TrkrDefs::getLayer(cluskey), clusglob.z() - stateglob.z());
0204 }
0205
0206 h_nmaps->Fill(nmaps);
0207 h_nintt->Fill(nintt);
0208 h_ntpc->Fill(ntpc);
0209 h_nmms->Fill(nmms);
0210 }
0211
0212 m_event++;
0213 return Fun4AllReturnCodes::EVENT_OK;
0214 }
0215 std::vector<TrkrDefs::cluskey> CosmicTrackQA::get_cluster_keys(SvtxTrack *track)
0216 {
0217 std::vector<TrkrDefs::cluskey> out;
0218 for (const auto &seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0219 {
0220 if (seed)
0221 {
0222 std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0223 }
0224 }
0225 return out;
0226 }
0227
0228 std::tuple<float, float, float, float>
0229 CosmicTrackQA::lineFitClusters(std::vector<Acts::Vector3> &positions)
0230 {
0231 TrackFitUtils::position_vector_t xypoints;
0232 TrackFitUtils::position_vector_t yzpoints;
0233 for (auto &pos : positions)
0234 {
0235 float clusr = std::sqrt(QAG4Util::square(pos.x()) + QAG4Util::square(pos.y()));
0236
0237 if (std::fabs(clusr) > 80 || std::fabs(clusr) < 30)
0238 {
0239 continue;
0240 }
0241 yzpoints.emplace_back(pos.z(), pos.y());
0242 xypoints.emplace_back(pos.x(), pos.y());
0243 }
0244
0245 auto xyparams = TrackFitUtils::line_fit(xypoints);
0246 auto yzparams = TrackFitUtils::line_fit(yzpoints);
0247
0248 return std::make_tuple(std::get<0>(xyparams),
0249 std::get<1>(xyparams),
0250 std::get<0>(yzparams),
0251 std::get<1>(yzparams));
0252 }
0253
0254 int CosmicTrackQA::EndRun(const int runnumber)
0255 {
0256 auto *hm = QAHistManagerDef::getHistoManager();
0257 assert(hm);
0258 TH2 *h_tracksperevent = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "ntracksperrun")));
0259
0260 h_tracksperevent->Fill(runnumber, m_tracks / m_event);
0261 return Fun4AllReturnCodes::EVENT_OK;
0262 }
0263
0264
0265 int CosmicTrackQA::End(PHCompositeNode * )
0266 {
0267 return Fun4AllReturnCodes::EVENT_OK;
0268 }
0269
0270 std::string CosmicTrackQA::getHistoPrefix() const
0271 {
0272 return std::string("h_") + Name() + std::string("_");
0273 }
0274
0275 void CosmicTrackQA::createHistos()
0276 {
0277 auto *hm = QAHistManagerDef::getHistoManager();
0278 assert(hm);
0279 {
0280 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpc").c_str(),
0281 "TPC clusters per track", 150, 0, 150);
0282 h->GetXaxis()->SetTitle("nTPC");
0283 hm->registerHisto(h);
0284 }
0285 {
0286 auto *h = new TH1F(std::string(getHistoPrefix() + "nmaps").c_str(),
0287 "MVTX clusters per track", 20, 0, 20);
0288 h->GetXaxis()->SetTitle("nMVTX");
0289 hm->registerHisto(h);
0290 }
0291 {
0292 auto *h = new TH1F(std::string(getHistoPrefix() + "nintt").c_str(),
0293 "INTT clusters per track", 20, 0, 20);
0294 h->GetXaxis()->SetTitle("nINTT");
0295 hm->registerHisto(h);
0296 }
0297 {
0298 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpot").c_str(),
0299 "TPOT clusters per track", 6, 0, 6);
0300 h->GetXaxis()->SetTitle("nTPOT");
0301 hm->registerHisto(h);
0302 }
0303 {
0304 auto *h = new TH2F(std::string(getHistoPrefix() + "lxlineresiduals").c_str(),
0305 "Local x HF fit residuals", 57, 0, 57, 1000, -10, 10);
0306 h->GetXaxis()->SetTitle("Layer");
0307 h->GetYaxis()->SetTitle("l_{x}^{clus}-l_{x}^{state} [cm]");
0308 hm->registerHisto(h);
0309 }
0310 {
0311 auto *h = new TH2F(std::string(getHistoPrefix() + "lzlineresiduals").c_str(),
0312 "Local z HF fit residuals", 57, 0, 57, 1000, -10, 10);
0313 h->GetXaxis()->SetTitle("Layer");
0314 h->GetYaxis()->SetTitle("l_{z}^{clus}-l_{z}^{state} [cm]");
0315 hm->registerHisto(h);
0316 }
0317 {
0318 auto *h = new TH2F(std::string(getHistoPrefix() + "gzlineresiduals").c_str(),
0319 "Global z HF fit residuals", 57, 0, 57, 1000, -10, 10);
0320 h->GetXaxis()->SetTitle("Layer");
0321 h->GetYaxis()->SetTitle("g_{z}^{clus}-g_{z}^{state} [cm]");
0322 hm->registerHisto(h);
0323 }
0324 {
0325 auto *h = new TH2F(std::string(getHistoPrefix() + "gylineresiduals").c_str(),
0326 "Global y HF fit residuals", 57, 0, 57, 1000, -10, 10);
0327 h->GetXaxis()->SetTitle("Layer");
0328 h->GetYaxis()->SetTitle("g_{y}^{clus}-g_{y}^{state} [cm]");
0329 hm->registerHisto(h);
0330 }
0331 {
0332 auto *h = new TH2F(std::string(getHistoPrefix() + "gxlineresiduals").c_str(),
0333 "Global x HF fit residuals", 57, 0, 57, 1000, -10, 10);
0334 h->GetXaxis()->SetTitle("Layer");
0335 h->GetYaxis()->SetTitle("g_{x}^{clus}-g_{x}^{state} [cm]");
0336 hm->registerHisto(h);
0337 }
0338 {
0339 auto *h = new TH2F(std::string(getHistoPrefix() + "nrecotracks").c_str(),
0340 "Num reconstructed tracks", 300, -3.14159, 3.1459, 100, -1.1, 1.1);
0341 h->GetXaxis()->SetTitle("#phi [rad]");
0342 h->GetYaxis()->SetTitle("#eta");
0343 hm->registerHisto(h);
0344 }
0345
0346 {
0347 auto *h = new TH2F(std::string(getHistoPrefix() + "ntracksperrun").c_str(),
0348 "Num reconstructed tracks per run", m_runbins, m_beginRun, m_endRun, 100, 0, 1);
0349 h->GetYaxis()->SetTitle("N_{tracks}/event");
0350 h->GetXaxis()->SetTitle("Run number");
0351 hm->registerHisto(h);
0352 }
0353
0354 {
0355 auto *h = new TH2F(std::string(getHistoPrefix() + "lxfitresiduals").c_str(),
0356 "Local x KF fit residuals", 57, 0, 57, 1000, -10, 10);
0357 h->GetXaxis()->SetTitle("Layer");
0358 h->GetYaxis()->SetTitle("l_{x}^{clus}-l_{x}^{state} [cm]");
0359 hm->registerHisto(h);
0360 }
0361 {
0362 auto *h = new TH2F(std::string(getHistoPrefix() + "lzfitresiduals").c_str(),
0363 "Local z KF fit residuals", 57, 0, 57, 1000, -10, 10);
0364 h->GetXaxis()->SetTitle("Layer");
0365 h->GetYaxis()->SetTitle("l_{z}^{clus}-l_{z}^{state} [cm]");
0366 hm->registerHisto(h);
0367 }
0368 {
0369 auto *h = new TH2F(std::string(getHistoPrefix() + "gzfitresiduals").c_str(),
0370 "Global z KF fit residuals", 57, 0, 57, 1000, -10, 10);
0371 h->GetXaxis()->SetTitle("Layer");
0372 h->GetYaxis()->SetTitle("g_{z}^{clus}-g_{z}^{state} [cm]");
0373 hm->registerHisto(h);
0374 }
0375 {
0376 auto *h = new TH2F(std::string(getHistoPrefix() + "gyfitresiduals").c_str(),
0377 "Global y KF fit residuals", 57, 0, 57, 1000, -10, 10);
0378 h->GetXaxis()->SetTitle("Layer");
0379 h->GetYaxis()->SetTitle("g_{y}^{clus}-g_{y}^{state} [cm]");
0380 hm->registerHisto(h);
0381 }
0382 {
0383 auto *h = new TH2F(std::string(getHistoPrefix() + "gxfitresiduals").c_str(),
0384 "Global x KF fit residuals", 57, 0, 57, 1000, -10, 10);
0385 h->GetXaxis()->SetTitle("Layer");
0386 h->GetYaxis()->SetTitle("g_{x}^{clus}-g_{x}^{state} [cm]");
0387 hm->registerHisto(h);
0388 }
0389 }
