QA/Tracking/TpcClusterQA.cc

0001 #include "TpcClusterQA.h"
0002
0003 #include <g4detectors/PHG4TpcCylinderGeom.h>
0004 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0005
0006 #include <trackbase/ActsGeometry.h>
0007 #include <trackbase/TpcDefs.h>
0008 #include <trackbase/TrackFitUtils.h>
0009 #include <trackbase/TrkrCluster.h>
0010 #include <trackbase/TrkrClusterContainer.h>
0011 #include <trackbase/TrkrClusterHitAssoc.h>
0012 #include <trackbase/TrkrDefs.h>
0013 #include <trackbase/TrkrHit.h>
0014 #include <trackbase/TrkrHitSet.h>
0015 #include <trackbase/TrkrHitSetContainer.h>
0016
0017 #include <tpc/TpcDistortionCorrectionContainer.h>
0018 #include <tpc/TpcGlobalPositionWrapper.h>
0019
0020 #include <qautils/QAHistManagerDef.h>
0021 #include <qautils/QAUtil.h>
0022
0023 #include <fun4all/Fun4AllHistoManager.h>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <fun4all/SubsysReco.h>
0026
0027 #include <phool/PHCompositeNode.h>
0028 #include <phool/getClass.h>
0029
0030 #include <TH1.h>
0031 #include <TH2.h>
0032
0033 #include <format>
0034
0035 //____________________________________________________________________________..
0036 TpcClusterQA::TpcClusterQA(const std::string &name)
0037   : SubsysReco(name)
0038 {
0039 }
0040
0041 //____________________________________________________________________________..
0042 int TpcClusterQA::InitRun(PHCompositeNode *topNode)
0043 {
0044   // find tpc geometry
0045   auto *geomContainer =
0046       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0047   if (!geomContainer)
0048   {
0049     std::cout << PHWHERE << " unable to find DST node CYLINDERCELLGEOM_SVTX" << std::endl;
0050     return Fun4AllReturnCodes::ABORTRUN;
0051   }
0052   // TPC has 3 regions, inner, mid and outer
0053   std::vector<int> region_layer_low = {7, 23, 39};
0054   std::vector<int> region_layer_high = {22, 38, 54};
0055
0056   // make a layer to region multimap
0057   const auto range = geomContainer->get_begin_end();
0058   for (auto iter = range.first; iter != range.second; ++iter)
0059   {
0060     m_layers.insert(iter->first);
0061
0062     for (int region = 0; region < 3; ++region)
0063     {
0064       if (iter->first >= region_layer_low[region] && iter->first <= region_layer_high[region])
0065       {
0066         m_layerRegionMap.insert(std::make_pair(iter->first, region));
0067       }
0068     }
0069   }
0070
0071   createHistos();
0072
0073   return Fun4AllReturnCodes::EVENT_OK;
0074 }
0075
0076 //____________________________________________________________________________..
0077 int TpcClusterQA::process_event(PHCompositeNode *topNode)
0078 {
0079   auto *clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0080   if (!clusterContainer)
0081   {
0082     return Fun4AllReturnCodes::ABORTEVENT;
0083   }
0084   auto *geomContainer =
0085       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0086   auto *hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0087   if (!hitmap)
0088   {
0089     std::cout << PHWHERE << "No hitmap found, bailing" << std::endl;
0090     return Fun4AllReturnCodes::ABORTEVENT;
0091   }
0092   auto *tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0093   if (!tGeometry)
0094   {
0095     std::cout << PHWHERE << "No acts geometry on node tree, bailing" << std::endl;
0096     return Fun4AllReturnCodes::ABORTEVENT;
0097   }
0098
0099   struct HistoList
0100   {
0101     TH1 *crphisize_side0 = nullptr;
0102     TH1 *crphisize_side1 = nullptr;
0103     TH1 *czsize = nullptr;
0104     TH1 *crphierr = nullptr;
0105     TH1 *czerr = nullptr;
0106     TH1 *cedge = nullptr;
0107     TH1 *coverlap = nullptr;
0108     TH1 *cxposition_side0 = nullptr;
0109     TH1 *cxposition_side1 = nullptr;
0110     TH1 *cyposition_side0 = nullptr;
0111     TH1 *cyposition_side1 = nullptr;
0112     TH1 *czposition_side0 = nullptr;
0113     TH1 *czposition_side1 = nullptr;
0114   };
0115
0116   int hitsetkeynum = 0;
0117   using HistoMap = std::map<int, HistoList>;
0118   HistoMap histos;
0119
0120   for (const auto &region : {0, 1, 2})
0121   {
0122     HistoList hist;
0123     hist.crphisize_side0 = h_phisize_side0[region];
0124     hist.crphisize_side1 = h_phisize_side1[region];
0125     hist.czsize = h_zsize[region];
0126     hist.crphierr = h_rphierror[region];
0127     hist.czerr = h_zerror[region];
0128     hist.cedge = h_clusedge[region];
0129     hist.coverlap = h_clusoverlap[region];
0130
0131     hist.cxposition_side0 = h_clusxposition_side0[region];
0132     hist.cxposition_side1 = h_clusxposition_side1[region];
0133
0134     hist.cyposition_side0 = h_clusyposition_side0[region];
0135     hist.cyposition_side1 = h_clusyposition_side1[region];
0136
0137     hist.czposition_side0 = h_cluszposition_side0[region];
0138     hist.czposition_side1 = h_cluszposition_side1[region];
0139
0140     histos.insert(std::make_pair(region, hist));
0141   }
0142   auto fill = [](TH1 *h, float val)
0143   { if (h) { h->Fill(val);
0144 } };
0145
0146   TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets(TrkrDefs::TrkrId::tpcId);
0147   for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
0148        hitsetiter != all_hitsets.second;
0149        ++hitsetiter)
0150   {
0151     auto hitsetkey = hitsetiter->first;
0152     TrkrHitSet *hitset = hitsetiter->second;
0153     if (TrkrDefs::getTrkrId(hitsetkey) != TrkrDefs::TrkrId::tpcId)
0154     {
0155       continue;
0156     }
0157     int hitlayer = TrkrDefs::getLayer(hitsetkey);
0158     // auto sector = TpcDefs::getSectorId(hitsetkey);
0159     auto m_side = TpcDefs::getSide(hitsetkey);
0160     auto hitrangei = hitset->getHits();
0161     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0162          hitr != hitrangei.second;
0163          ++hitr)
0164     {
0165       auto hitkey = hitr->first;
0166       // auto hit = hitr->second;
0167       // auto adc = hit->getAdc();
0168       auto hitpad = TpcDefs::getPad(hitkey);
0169       auto m_hittbin = TpcDefs::getTBin(hitkey);
0170       // Check TrackResiduals.cc
0171       auto *geoLayer = geomContainer->GetLayerCellGeom(hitlayer);
0172       auto phi = geoLayer->get_phicenter(hitpad, m_side);
0173       auto radius = geoLayer->get_radius();
0174       auto m_hitgx = radius * std::cos(phi);
0175       auto m_hitgy = radius * std::sin(phi);
0176       float AdcClockPeriod = geoLayer->get_zstep();
0177       float m_zdriftlength = m_hittbin * tGeometry->get_drift_velocity() * AdcClockPeriod;
0178       double NZBinsSide = 249;  // physical z bins per TPC side
0179       double tdriftmax = AdcClockPeriod * NZBinsSide;
0180       auto m_hitgz = (tdriftmax * tGeometry->get_drift_velocity()) - m_zdriftlength;
0181       if (m_side == 0)
0182       {
0183         m_hitgz *= -1;
0184       }
0185       // geoLayer->identify(std::cout);
0186       h_hitpositions->Fill(m_hitgx, m_hitgy);
0187       if (m_side == 0)
0188       {
0189         h_hitzpositions_side0->Fill(m_hitgz);
0190       }
0191       if (m_side == 1)
0192       {
0193         h_hitzpositions_side1->Fill(m_hitgz);
0194       }
0195     }
0196   }
0197   float nclusperevent[24] = {0};
0198   for (auto &hsk : clusterContainer->getHitSetKeys(TrkrDefs::TrkrId::tpcId))
0199   {
0200     int numclusters = 0;
0201     auto range = clusterContainer->getClusters(hsk);
0202     int sector = TpcDefs::getSectorId(hsk);
0203     int side = TpcDefs::getSide(hsk);
0204     if (side > 0)
0205     {
0206       sector += 12;
0207     }
0208     for (auto iter = range.first; iter != range.second; ++iter)
0209     {
0210       const auto cluskey = iter->first;
0211       auto *const cluster = iter->second;  // auto cluster = clusters->findCluster(key);
0212       auto glob = tGeometry->getGlobalPosition(cluskey, cluster);
0213       auto sclusgx = glob.x();
0214       auto sclusgy = glob.y();
0215       auto sclusgz = glob.z();
0216
0217       const auto it = m_layerRegionMap.find(TrkrDefs::getLayer(cluskey));
0218       int region = it->second;
0219       const auto hiter = histos.find(region);
0220       if (hiter == histos.end())
0221       {
0222         continue;
0223       }
0224       fill(hiter->second.czsize, cluster->getZSize());
0225       fill(hiter->second.crphierr, cluster->getRPhiError());
0226       fill(hiter->second.czerr, cluster->getZError());
0227       fill(hiter->second.cedge, cluster->getEdge());
0228       fill(hiter->second.coverlap, cluster->getOverlap());
0229
0230       if (side == 0)
0231       {
0232         fill(hiter->second.crphisize_side0, cluster->getPhiSize());
0233         fill(hiter->second.cxposition_side0, sclusgx);
0234         fill(hiter->second.cyposition_side0, sclusgy);
0235         fill(hiter->second.czposition_side0, sclusgz);
0236       }
0237       if (side == 1)
0238       {
0239         fill(hiter->second.crphisize_side1, cluster->getPhiSize());
0240         fill(hiter->second.cxposition_side1, sclusgx);
0241         fill(hiter->second.cyposition_side1, sclusgy);
0242         fill(hiter->second.czposition_side1, sclusgz);
0243       }
0244
0245       numclusters++;
0246     }
0247
0248     nclusperevent[sector] += numclusters;
0249     h_totalclusters->Fill(hitsetkeynum, numclusters);
0250     m_totalClusters += numclusters;
0251     hitsetkeynum++;
0252   }
0253   for (int i = 0; i < 24; i++)
0254   {
0255     h_clusterssector->Fill(i, nclusperevent[i]);
0256     m_clustersPerSector[i] += nclusperevent[i];
0257   }
0258
0259   m_event++;
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263 int TpcClusterQA::EndRun(const int /*runnumber*/)
0264 {
0265   return Fun4AllReturnCodes::EVENT_OK;
0266 }
0267
0268 std::string TpcClusterQA::getHistoPrefix() const
0269 {
0270   return std::string("h_") + Name() + std::string("_");
0271 }
0272 void TpcClusterQA::createHistos()
0273 {
0274   auto *hm = QAHistManagerDef::getHistoManager();
0275   assert(hm);
0276
0277   {
0278     h_clusterssector = new TH2F(std::string(getHistoPrefix() + "ncluspersector").c_str(),
0279                                 "TPC Clusters per event per sector", 24, 0, 24, 5000, 0, 5000);
0280     h_clusterssector->GetXaxis()->SetTitle("Sector number");
0281     h_clusterssector->GetYaxis()->SetTitle("Clusters per event");
0282     hm->registerHisto(h_clusterssector);
0283   }
0284   for (const auto &region : {0, 1, 2})
0285   {
0286     {
0287       h_phisize_side0[region] = new TH1F(std::format("{}phisize_side0_{}", getHistoPrefix(), region).c_str(),
0288                                          std::format("TPC (side 0) cluster #phi size region_{}", region).c_str(), 10, 0, 10);
0289       h_phisize_side0[region]->GetXaxis()->SetTitle("Cluster #phi_{size}");
0290       hm->registerHisto(h_phisize_side0[region]);
0291     }
0292     {
0293       h_phisize_side1[region] = new TH1F(std::format("{}phisize_side1_{}", getHistoPrefix(), region).c_str(),
0294                                          std::format("TPC (side 1) cluster #phi size region_{}", region).c_str(), 10, 0, 10);
0295       h_phisize_side1[region]->GetXaxis()->SetTitle("Cluster #phi_{size}");
0296       hm->registerHisto(h_phisize_side1[region]);
0297     }
0298     {
0299       h_zsize[region] = new TH1F(std::format("{}zsize_{}", getHistoPrefix(), region).c_str(),
0300                                  std::format("TPC cluster z size region_{}", region).c_str(), 10, 0, 10);
0301       h_zsize[region]->GetXaxis()->SetTitle("Cluster z_{size}");
0302       hm->registerHisto(h_zsize[region]);
0303     }
0304     {
0305       h_rphierror[region] = new TH1F(std::format("{}rphi_error_{}", getHistoPrefix(), region).c_str(),
0306                                      std::format("TPC r#Delta#phi error region_{}", region).c_str(), 100, 0, 0.075);
0307       h_rphierror[region]->GetXaxis()->SetTitle("r#Delta#phi error [cm]");
0308       hm->registerHisto(h_rphierror[region]);
0309     }
0310     {
0311       h_zerror[region] = new TH1F(std::format("{}z_error_{}", getHistoPrefix(), region).c_str(),
0312                                   std::format("TPC z error region_{}", region).c_str(), 100, 0, 0.18);
0313       h_zerror[region]->GetXaxis()->SetTitle("z error [cm]");
0314       hm->registerHisto(h_zerror[region]);
0315     }
0316     {
0317       h_clusedge[region] = new TH1F(std::format("{}clusedge_{}", getHistoPrefix(), region).c_str(),
0318                                     std::format("TPC hits on edge region_{}", region).c_str(), 30, 0, 30);
0319       h_clusedge[region]->GetXaxis()->SetTitle("Cluster edge");
0320       hm->registerHisto(h_clusedge[region]);
0321     }
0322     {
0323       h_clusoverlap[region] = new TH1F(std::format("{}clusoverlap_{}", getHistoPrefix(), region).c_str(),
0324                                        std::format("TPC clus overlap region_{}", region).c_str(), 30, 0, 30);
0325       h_clusoverlap[region]->GetXaxis()->SetTitle("Cluster overlap");
0326       hm->registerHisto(h_clusoverlap[region]);
0327     }
0328     {
0329       h_clusxposition_side0[region] = new TH1F(std::format("{}clusxposition_side0_{}", getHistoPrefix(), region).c_str(),
0330                                                std::format("TPC cluster x position side 0 region_{}", region).c_str(), 210 * 2, -105, 105);
0331       h_clusxposition_side0[region]->GetXaxis()->SetTitle("x (cm)");
0332       hm->registerHisto(h_clusxposition_side0[region]);
0333     }
0334     {
0335       h_clusxposition_side1[region] = new TH1F(std::format("{}clusxposition_side1_{}", getHistoPrefix(), region).c_str(),
0336                                                std::format("TPC cluster x position side 1 region_{}", region).c_str(), 210 * 2, -105, 105);
0337       h_clusxposition_side1[region]->GetXaxis()->SetTitle("x (cm)");
0338       hm->registerHisto(h_clusxposition_side1[region]);
0339     }
0340     {
0341       h_clusyposition_side0[region] = new TH1F(std::format("{}clusyposition_side0_{}", getHistoPrefix(), region).c_str(),
0342                                                std::format("TPC cluster y position side 0 region_{}", region).c_str(), 210 * 2, -105, 105);
0343       h_clusyposition_side0[region]->GetXaxis()->SetTitle("y (cm)");
0344       hm->registerHisto(h_clusyposition_side0[region]);
0345     }
0346     {
0347       h_clusyposition_side1[region] = new TH1F(std::format("{}clusyposition_side1_{}", getHistoPrefix(), region).c_str(),
0348                                                std::format("TPC cluster y position side 1 region_{}", region).c_str(), 210 * 2, -105, 105);
0349       h_clusyposition_side1[region]->GetXaxis()->SetTitle("y (cm)");
0350       hm->registerHisto(h_clusyposition_side1[region]);
0351     }
0352     {
0353       h_cluszposition_side0[region] = new TH1F(std::format("{}cluszposition_side0_{}", getHistoPrefix(), region).c_str(),
0354                                                std::format("TPC cluster z position side 0 region_{}", region).c_str(), 210 * 2, -105, 105);
0355       h_cluszposition_side0[region]->GetXaxis()->SetTitle("z (cm)");
0356       hm->registerHisto(h_cluszposition_side0[region]);
0357     }
0358     {
0359       h_cluszposition_side1[region] = new TH1F(std::format("{}cluszposition_side1_{}", getHistoPrefix(), region).c_str(),
0360                                                std::format("TPC cluster z position side 1 region_{}", region).c_str(), 210 * 2, -105, 105);
0361       h_cluszposition_side1[region]->GetXaxis()->SetTitle("z (cm)");
0362       hm->registerHisto(h_cluszposition_side1[region]);
0363     }
0364   }
0365
0366   {
0367     h_totalclusters = new TH2F(std::string(getHistoPrefix() + "stotal_clusters").c_str(),
0368                                "TPC clusters per hitsetkey", 1152, 0, 1152, 10000, 0, 10000);
0369     h_totalclusters->GetXaxis()->SetTitle("Hitsetkey number");
0370     h_totalclusters->GetYaxis()->SetTitle("Number of clusters");
0371     hm->registerHisto(h_totalclusters);
0372   }
0373
0374   {
0375     h_hitpositions = new TH2F(std::string(getHistoPrefix() + "hit_positions").c_str(),
0376                               "Histogram of hit x y positions", 160, 0, 80, 160, 0, 80);
0377     h_hitpositions->GetXaxis()->SetTitle("x (cm)");
0378     h_hitpositions->GetYaxis()->SetTitle("y (cm)");
0379     hm->registerHisto(h_hitpositions);
0380   }
0381   {
0382     h_hitzpositions_side0 = new TH1F(std::string(getHistoPrefix() + "hitz_positions_side0").c_str(),
0383                                      "Histogram of hit z positions side=0", 105 * 4, -105, 105);
0384     h_hitzpositions_side0->GetXaxis()->SetTitle("z (cm)");
0385     hm->registerHisto(h_hitzpositions_side0);
0386   }
0387   {
0388     h_hitzpositions_side1 = new TH1F(std::string(getHistoPrefix() + "hitz_positions_side1").c_str(),
0389                                      "Histogram of hit z positions side=1", 105 * 4, -105, 105);
0390     h_hitzpositions_side1->GetXaxis()->SetTitle("z (cm)");
0391     hm->registerHisto(h_hitzpositions_side1);
0392   }
0393   return;
0394 }