QA/Tracking/InttClusterQA.cc

0001 #include "InttClusterQA.h"
0002
0003 #include <intt/CylinderGeomIntt.h>
0004
0005 #include <g4detectors/PHG4CylinderGeomContainer.h>
0006
0007 #include <trackbase/ActsGeometry.h>
0008 #include <trackbase/InttDefs.h>
0009 #include <trackbase/TrackFitUtils.h>
0010 #include <trackbase/TrkrCluster.h>
0011 #include <trackbase/TrkrClusterContainer.h>
0012 #include <trackbase/TrkrClusterHitAssoc.h>
0013 #include <trackbase/TrkrDefs.h>
0014 #include <trackbase/TrkrHitSet.h>
0015 #include <trackbase/TrkrHitSetContainerv1.h>
0016
0017 #include <qautils/QAHistManagerDef.h>
0018 #include <qautils/QAUtil.h>
0019
0020 #include <fun4all/Fun4AllHistoManager.h>
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/SubsysReco.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/getClass.h>
0026
0027 #include <TH1.h>
0028 #include <TH2.h>
0029
0030 #include <format>
0031
0032 //____________________________________________________________________________..
0033 InttClusterQA::InttClusterQA(const std::string &name)
0034   : SubsysReco(name)
0035 {
0036 }
0037
0038 //____________________________________________________________________________..
0039 int InttClusterQA::InitRun(PHCompositeNode * /*unused*/)
0040 {
0041   for (const auto &layer : {0, 1, 2, 3})
0042   {
0043     if (layer < 2)
0044     {
0045       m_layerLadderMap.insert(std::make_pair(layer, 12));
0046     }
0047     else
0048     {
0049       m_layerLadderMap.insert(std::make_pair(layer, 16));
0050     }
0051   }
0052
0053   createHistos();
0054   return Fun4AllReturnCodes::EVENT_OK;
0055 }
0056
0057 //____________________________________________________________________________..
0058 int InttClusterQA::process_event(PHCompositeNode *topNode)
0059 {
0060   auto *clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0061   if (!clusterContainer)
0062   {
0063     std::cout << PHWHERE << "No cluster container, bailing" << std::endl;
0064     return Fun4AllReturnCodes::ABORTEVENT;
0065   }
0066
0067   auto *trkrHitSetContainer = findNode::getClass<TrkrHitSetContainerv1>(topNode, "TRKR_HITSET");
0068   if (!trkrHitSetContainer)
0069   {
0070     std::cout << PHWHERE << "No trkrhitset container, bailing" << std::endl;
0071     return Fun4AllReturnCodes::ABORTEVENT;
0072   }
0073
0074   auto *tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0075   if (!tGeometry)
0076   {
0077     std::cout << PHWHERE << "No acts geometry on node tree, bailing" << std::endl;
0078     return Fun4AllReturnCodes::ABORTEVENT;
0079   }
0080
0081   for (auto &hsk : clusterContainer->getHitSetKeys(TrkrDefs::TrkrId::inttId))
0082   {
0083     int numclusters = 0;
0084     auto range = clusterContainer->getClusters(hsk);
0085     auto layer = TrkrDefs::getLayer(hsk);
0086     auto ladderphiid = InttDefs::getLadderPhiId(hsk);
0087     auto sensor = InttDefs::getLadderZId(hsk);
0088
0089     if (m_sensorInfo)
0090     {
0091       for (auto iter = range.first; iter != range.second; ++iter)
0092       {
0093         const auto cluskey = iter->first;
0094         auto *const cluster = iter->second;
0095         auto globalpos = tGeometry->getGlobalPosition(cluskey, cluster);
0096         auto phi = atan2(globalpos(1), globalpos(0));
0097         auto clayer = TrkrDefs::getLayer(cluskey);
0098         h_cluspersensor[(int) (layer) -3][(int) ladderphiid][(int) sensor]->Fill(cluster->getLocalY(), cluster->getLocalX());
0099         h_clusPhi_incl->Fill(phi);
0100         if (clayer == 3 || clayer == 4)
0101         {
0102           h_clusPhi_l34->Fill(phi);
0103           h_clusZ_clusPhi_l34->Fill(globalpos(2), phi);
0104         }
0105         else if (clayer == 5 || clayer == 6)
0106         {
0107           h_clusPhi_l56->Fill(phi);
0108           h_clusZ_clusPhi_l56->Fill(globalpos(2), phi);
0109         }
0110
0111         m_totalClusters++;
0112         numclusters++;
0113       }
0114       m_nclustersPerSensor[((int) layer) - 3][(int) ladderphiid][(int) sensor] += numclusters;
0115     }
0116     else
0117     {
0118       for (auto iter = range.first; iter != range.second; ++iter)
0119       {
0120         const auto cluskey = iter->first;
0121         auto *const cluster = iter->second;
0122         auto globalpos = tGeometry->getGlobalPosition(cluskey, cluster);
0123         auto phi = atan2(globalpos(1), globalpos(0));
0124         auto clayer = TrkrDefs::getLayer(cluskey);
0125         h_clusSize->Fill(cluster->getSize());
0126         h_clusPhi_incl->Fill(phi);
0127         if (clayer == 3 || clayer == 4)
0128         {
0129           h_clusPhi_l34->Fill(phi);
0130           h_clusZ_clusPhi_l34->Fill(globalpos(2), phi);
0131         }
0132         else if (clayer == 5 || clayer == 6)
0133         {
0134           h_clusPhi_l56->Fill(phi);
0135           h_clusZ_clusPhi_l56->Fill(globalpos(2), phi);
0136         }
0137       }
0138     }
0139   }
0140
0141   TrkrHitSetContainer::ConstRange hitsetrange = trkrHitSetContainer->getHitSets(TrkrDefs::TrkrId::inttId);
0142
0143   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first; hitsetitr != hitsetrange.second; ++hitsetitr)
0144   {
0145     int sensor_hits = hitsetitr->second->size();
0146     float sensor_occupancy = (float) sensor_hits / (128. * 26.);
0147     h_occupancy->Fill(100. * sensor_occupancy);
0148   }
0149
0150   m_event++;
0151
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154 int InttClusterQA::EndRun(const int /*runnumber*/)
0155 {
0156   return Fun4AllReturnCodes::EVENT_OK;
0157 }
0158
0159 std::string InttClusterQA::getHistoPrefix() const
0160 {
0161   return std::string("h_") + Name() + std::string("_");
0162 }
0163
0164 void InttClusterQA::createHistos()
0165 {
0166   auto *hm = QAHistManagerDef::getHistoManager();
0167   assert(hm);
0168
0169   h_occupancy = new TH1F(std::format("{}sensorOccupancy", getHistoPrefix()).c_str(), "INTT Sensor Occupancy", 100, 0, 5);
0170   h_occupancy->GetXaxis()->SetTitle("Sensor Occupancy [%]");
0171   h_occupancy->GetYaxis()->SetTitle("Entries");
0172   hm->registerHisto(h_occupancy);
0173   h_clusSize = new TH1F(std::format("{}clusterSize", getHistoPrefix()).c_str(), "INTT Cluster Size", 20, -0.5, 19.5);
0174   h_clusSize->GetXaxis()->SetTitle("Cluster Size");
0175   h_clusSize->GetYaxis()->SetTitle("Entries");
0176   hm->registerHisto(h_clusSize);
0177   h_clusPhi_incl = new TH1F(std::format("{}clusterPhi_incl", getHistoPrefix()).c_str(), "INTT Cluster Phi", 320, -3.2, 3.2);
0178   h_clusPhi_incl->GetXaxis()->SetTitle("Cluster (inner+outer) #phi [rad]");
0179   h_clusPhi_incl->GetYaxis()->SetTitle("Entries");
0180   hm->registerHisto(h_clusPhi_incl);
0181   h_clusPhi_l34 = new TH1F(std::format("{}clusterPhi_l34", getHistoPrefix()).c_str(), "INTT Cluster Phi", 320, -3.2, 3.2);
0182   h_clusPhi_l34->GetXaxis()->SetTitle("Cluster (inner) #phi [rad]");
0183   h_clusPhi_l34->GetYaxis()->SetTitle("Entries");
0184   hm->registerHisto(h_clusPhi_l34);
0185   h_clusPhi_l56 = new TH1F(std::format("{}clusterPhi_l56", getHistoPrefix()).c_str(), "INTT Cluster Phi", 320, -3.2, 3.2);
0186   h_clusPhi_l56->GetXaxis()->SetTitle("Cluster (outer) #phi [rad]");
0187   h_clusPhi_l56->GetYaxis()->SetTitle("Entries");
0188   hm->registerHisto(h_clusPhi_l56);
0189   h_clusZ_clusPhi_l34 = new TH2F(std::format("{}clusterZ_clusPhi_l34", getHistoPrefix()).c_str(), "INTT Cluster Z vs Cluster Phi", 55, cluszbin, 350, -3.5, 3.5);
0190   h_clusZ_clusPhi_l34->GetXaxis()->SetTitle("Cluster (inner) Z [cm]");
0191   h_clusZ_clusPhi_l34->GetYaxis()->SetTitle("Cluster (inner) #phi [rad]");
0192   hm->registerHisto(h_clusZ_clusPhi_l34);
0193   h_clusZ_clusPhi_l56 = new TH2F(std::format("{}clusterZ_clusPhi_l56", getHistoPrefix()).c_str(), "INTT Cluster Z vs Cluster Phi", 55, cluszbin, 350, -3.5, 3.5);
0194   h_clusZ_clusPhi_l56->GetXaxis()->SetTitle("Cluster (outer) Z [cm]");
0195   h_clusZ_clusPhi_l56->GetYaxis()->SetTitle("Cluster (outer) #phi [rad]");
0196   hm->registerHisto(h_clusZ_clusPhi_l56);
0197
0198   if (m_sensorInfo)
0199   {
0200     for (const auto &[layer, ladders] : m_layerLadderMap)
0201     {
0202       for (int ladder = 0; ladder < ladders; ladder++)
0203       {
0204         //! 4 sensor on each ladder
0205         for (int sensor = 0; sensor < 4; sensor++)
0206         {
0207           h_cluspersensor[layer][ladder][sensor] = new TH2F(std::format("{}ncluspersensor{}_{}_{}", getHistoPrefix(), layer, ladder, sensor).c_str(),
0208                                                             "INTT clusters per sensor", 100, -5, 5, 1000, -1, 1);
0209           h_cluspersensor[layer][ladder][sensor]->GetXaxis()->SetTitle("Local z [cm]");
0210           h_cluspersensor[layer][ladder][sensor]->GetYaxis()->SetTitle("Local rphi [cm]");
0211           hm->registerHisto(h_cluspersensor[layer][ladder][sensor]);
0212         }
0213       }
0214     }
0215   }
0216
0217   return;
0218 }