QA/Tracking/MvtxClusterQA.cc

0001
0002 #include "MvtxClusterQA.h"
0003
0004 #include <mvtx/CylinderGeom_Mvtx.h>
0005
0006 #include <g4detectors/PHG4CylinderGeomContainer.h>
0007
0008 #include <trackbase/ActsGeometry.h>
0009 #include <trackbase/MvtxDefs.h>
0010 #include <trackbase/TrackFitUtils.h>
0011 #include <trackbase/TrkrCluster.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013 #include <trackbase/TrkrClusterHitAssoc.h>
0014 #include <trackbase/TrkrDefs.h>
0015 #include <trackbase/TrkrHitSet.h>
0016 #include <trackbase/TrkrHitSetContainerv1.h>
0017
0018 #include <qautils/QAHistManagerDef.h>
0019 #include <qautils/QAUtil.h>
0020
0021 #include <fun4all/Fun4AllHistoManager.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/SubsysReco.h>
0024
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/getClass.h>
0027
0028 #include <TH1.h>
0029 #include <TH2.h>
0030
0031 #include <format>
0032
0033 //____________________________________________________________________________..
0034 MvtxClusterQA::MvtxClusterQA(const std::string &name)
0035   : SubsysReco(name)
0036 {
0037 }
0038
0039 //____________________________________________________________________________..
0040 int MvtxClusterQA::InitRun(PHCompositeNode *topNode)
0041 {
0042   auto *geomContainer = findNode::getClass<
0043       PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0044   if (!geomContainer)
0045   {
0046     std::cout << PHWHERE
0047               << " CYLINDERGEOM_MVTX  node not found on node tree"
0048               << std::endl;
0049     return Fun4AllReturnCodes::ABORTEVENT;
0050   }
0051   for (const auto &layer : {0, 1, 2})
0052   {
0053     auto *layergeom = dynamic_cast<CylinderGeom_Mvtx *>(geomContainer->GetLayerGeom(layer));
0054     if (!layergeom)
0055     {
0056       std::cout << PHWHERE << "Did not get layergeom for layer "
0057                 << layer << std::endl;
0058       return Fun4AllReturnCodes::ABORTEVENT;
0059     }
0060     int nstaves = layergeom->get_N_staves();
0061     m_layerStaveMap.insert(std::make_pair(layer, nstaves));
0062   }
0063
0064   createHistos();
0065   return Fun4AllReturnCodes::EVENT_OK;
0066 }
0067
0068 //____________________________________________________________________________..
0069 int MvtxClusterQA::process_event(PHCompositeNode *topNode)
0070 {
0071   auto *clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0072   if (!clusterContainer)
0073   {
0074     std::cout << PHWHERE << "No cluster container, bailing" << std::endl;
0075     return Fun4AllReturnCodes::ABORTEVENT;
0076   }
0077
0078   auto *trkrHitSetContainer = findNode::getClass<TrkrHitSetContainerv1>(topNode, "TRKR_HITSET");
0079   if (!trkrHitSetContainer)
0080   {
0081     std::cout << PHWHERE << "No trkrhitset container, bailing" << std::endl;
0082     return Fun4AllReturnCodes::ABORTEVENT;
0083   }
0084
0085   auto *tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0086   if (!tGeometry)
0087   {
0088     std::cout << PHWHERE << "No acts geometry on node tree, bailing" << std::endl;
0089     return Fun4AllReturnCodes::ABORTEVENT;
0090   }
0091
0092   int numclusters = 0;
0093   for (auto &hsk : clusterContainer->getHitSetKeys(TrkrDefs::TrkrId::mvtxId))
0094   {
0095     auto range = clusterContainer->getClusters(hsk);
0096     for (auto iter = range.first; iter != range.second; ++iter)
0097     {
0098       numclusters++;
0099     }
0100   }
0101
0102   for (auto &hsk : clusterContainer->getHitSetKeys(TrkrDefs::TrkrId::mvtxId))
0103   {
0104     int numclusters_chip = 0;
0105     auto range = clusterContainer->getClusters(hsk);
0106     auto layer = TrkrDefs::getLayer(hsk);
0107     auto stave = MvtxDefs::getStaveId(hsk);
0108     auto chip = MvtxDefs::getChipId(hsk);
0109
0110     if (m_chipInfo)
0111     {
0112       for (auto iter = range.first; iter != range.second; ++iter)
0113       {
0114         const auto cluskey = iter->first;
0115         auto *const cluster = iter->second;
0116         auto globalpos = tGeometry->getGlobalPosition(cluskey, cluster);
0117         auto phi = atan2(globalpos(1), globalpos(0));
0118         auto clayer = TrkrDefs::getLayer(cluskey);
0119         h_clusperchip[(int) layer][(int) stave][(int) chip]->Fill(cluster->getLocalY(), cluster->getLocalX());
0120         h_clusSize->Fill(cluster->getSize());
0121         h_clusSize_nClus->Fill(numclusters, cluster->getSize());
0122         h_clusPhi_incl->Fill(phi);
0123         if (clayer == 0)
0124         {
0125           h_clusPhi_l0->Fill(phi);
0126           h_clusZ_clusPhi_l0->Fill(globalpos(2), phi);
0127         }
0128         else if (clayer == 1)
0129         {
0130           h_clusPhi_l1->Fill(phi);
0131           h_clusZ_clusPhi_l1->Fill(globalpos(2), phi);
0132         }
0133         else if (clayer == 2)
0134         {
0135           h_clusPhi_l2->Fill(phi);
0136           h_clusZ_clusPhi_l2->Fill(globalpos(2), phi);
0137         }
0138         m_totalClusters++;
0139         numclusters_chip++;
0140       }
0141
0142       m_nclustersPerChip[(int) layer][(int) stave][(int) chip] += numclusters_chip;
0143     }
0144     else
0145     {
0146       for (auto iter = range.first; iter != range.second; ++iter)
0147       {
0148         const auto cluskey = iter->first;
0149         auto *const cluster = iter->second;
0150         auto globalpos = tGeometry->getGlobalPosition(cluskey, cluster);
0151         auto phi = atan2(globalpos(1), globalpos(0));
0152         auto clayer = TrkrDefs::getLayer(cluskey);
0153         h_clusSize->Fill(cluster->getSize());
0154         h_clusSize_nClus->Fill(numclusters, cluster->getSize());
0155         h_clusPhi_incl->Fill(phi);
0156         if (clayer == 0)
0157         {
0158           h_clusPhi_l0->Fill(phi);
0159           h_clusZ_clusPhi_l0->Fill(globalpos(2), phi);
0160         }
0161         else if (clayer == 1)
0162         {
0163           h_clusPhi_l1->Fill(phi);
0164           h_clusZ_clusPhi_l1->Fill(globalpos(2), phi);
0165         }
0166         else if (clayer == 2)
0167         {
0168           h_clusPhi_l2->Fill(phi);
0169           h_clusZ_clusPhi_l2->Fill(globalpos(2), phi);
0170         }
0171       }
0172     }
0173   }
0174
0175   TrkrHitSetContainer::ConstRange hitsetrange = trkrHitSetContainer->getHitSets(TrkrDefs::TrkrId::mvtxId);
0176
0177   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first; hitsetitr != hitsetrange.second; ++hitsetitr)
0178   {
0179     int chip_hits = hitsetitr->second->size();
0180     float chip_occupancy = (float) chip_hits / (512 * 1024);
0181     chip_occupancy = 100 * chip_occupancy;
0182     h_occupancy->Fill(chip_occupancy);
0183     int strobe = MvtxDefs::getStrobeId(hitsetitr->first);
0184     for (int i = 0; i < chip_hits; i++)
0185     {
0186       h_strobe->Fill(strobe);
0187     }
0188   }
0189
0190   m_event++;
0191   return Fun4AllReturnCodes::EVENT_OK;
0192 }
0193 int MvtxClusterQA::EndRun(const int /*runnumber*/)
0194 {
0195   return Fun4AllReturnCodes::EVENT_OK;
0196 }
0197 //____________________________________________________________________________..
0198
0199 std::string MvtxClusterQA::getHistoPrefix() const
0200 {
0201   return std::string("h_") + Name() + std::string("_");
0202 }
0203
0204 void MvtxClusterQA::createHistos()
0205 {
0206   auto *hm = QAHistManagerDef::getHistoManager();
0207   assert(hm);
0208
0209   h_occupancy = new TH1F(std::format("{}chipOccupancy", getHistoPrefix()).c_str(), "MVTX Chip Occupancy", 60, 0, 0.6);
0210   h_occupancy->GetXaxis()->SetTitle("Chip Occupancy [%]");
0211   h_occupancy->GetYaxis()->SetTitle("Entries");
0212   hm->registerHisto(h_occupancy);
0213   h_clusSize = new TH1F(std::format("{}clusterSize", getHistoPrefix()).c_str(), "MVTX Cluster Size", 50, -0.5, 49.5);
0214   h_clusSize->GetXaxis()->SetTitle("Cluster Size");
0215   h_clusSize->GetYaxis()->SetTitle("Entries");
0216   hm->registerHisto(h_clusSize);
0217   h_clusPhi_incl = new TH1F(std::format("{}clusterPhi_incl", getHistoPrefix()).c_str(), "MVTX Cluster Phi", 320, -3.2, 3.2);
0218   h_clusPhi_incl->GetXaxis()->SetTitle("Cluster (layer 0+1+2) #phi [rad]");
0219   h_clusPhi_incl->GetYaxis()->SetTitle("Entries");
0220   hm->registerHisto(h_clusPhi_incl);
0221   h_clusPhi_l0 = new TH1F(std::format("{}clusterPhi_l0", getHistoPrefix()).c_str(), "MVTX Cluster Phi", 320, -3.2, 3.2);
0222   h_clusPhi_l0->GetXaxis()->SetTitle("Cluster (layer 0) #phi [rad]");
0223   h_clusPhi_l0->GetYaxis()->SetTitle("Entries");
0224   hm->registerHisto(h_clusPhi_l0);
0225   h_clusPhi_l1 = new TH1F(std::format("{}clusterPhi_l1", getHistoPrefix()).c_str(), "MVTX Cluster Phi", 320, -3.2, 3.2);
0226   h_clusPhi_l1->GetXaxis()->SetTitle("Cluster (layer 1) #phi [rad]");
0227   h_clusPhi_l1->GetYaxis()->SetTitle("Entries");
0228   hm->registerHisto(h_clusPhi_l1);
0229   h_clusPhi_l2 = new TH1F(std::format("{}clusterPhi_l2", getHistoPrefix()).c_str(), "MVTX Cluster Phi", 320, -3.2, 3.2);
0230   h_clusPhi_l2->GetXaxis()->SetTitle("Cluster (layer 2) #phi [rad]");
0231   h_clusPhi_l2->GetYaxis()->SetTitle("Entries");
0232   hm->registerHisto(h_clusPhi_l2);
0233   h_clusZ_clusPhi_l0 = new TH2F(std::format("{}clusterZ_clusPhi_l0", getHistoPrefix()).c_str(), "MVTX Cluster Z vs Phi", 300, -15, 15, 350, -3.5, 3.5);
0234   h_clusZ_clusPhi_l0->GetXaxis()->SetTitle("Cluster (layer 2) Z [cm]");
0235   h_clusZ_clusPhi_l0->GetYaxis()->SetTitle("Cluster (layer 0) #phi [rad]");
0236   hm->registerHisto(h_clusZ_clusPhi_l0);
0237   h_clusZ_clusPhi_l1 = new TH2F(std::format("{}clusterZ_clusPhi_l1", getHistoPrefix()).c_str(), "MVTX Cluster Z vs Phi", 300, -15, 15, 350, -3.5, 3.5);
0238   h_clusZ_clusPhi_l1->GetXaxis()->SetTitle("Cluster (layer 2) Z [cm]");
0239   h_clusZ_clusPhi_l1->GetYaxis()->SetTitle("Cluster (layer 1) #phi [rad]");
0240   hm->registerHisto(h_clusZ_clusPhi_l1);
0241   h_clusZ_clusPhi_l2 = new TH2F(std::format("{}clusterZ_clusPhi_l2", getHistoPrefix()).c_str(), "MVTX Cluster Z vs Phi", 300, -15, 15, 350, -3.5, 3.5);
0242   h_clusZ_clusPhi_l2->GetXaxis()->SetTitle("Cluster (layer 2) Z [cm]");
0243   h_clusZ_clusPhi_l2->GetYaxis()->SetTitle("Cluster (layer 2) #phi [rad]");
0244   hm->registerHisto(h_clusZ_clusPhi_l2);
0245   h_strobe = new TH1I(std::format("{}strobeTiming", getHistoPrefix()).c_str(), "MVTX Strobe Timing per Hit", 32, -15, 16);
0246   h_strobe->GetXaxis()->SetTitle("Strobe BCO - GL1 BCO");
0247   h_strobe->GetYaxis()->SetTitle("Entries");
0248   hm->registerHisto(h_strobe);
0249   h_clusSize_nClus = new TH2F(std::format("{}clusSize_nCLus", getHistoPrefix()).c_str(), "MVTX Cluster Size vs Number of Clusters", 800, -0.5, 799.5, 25, -0.5, 24.5);
0250   h_clusSize_nClus->GetXaxis()->SetTitle("Number of Clusters");
0251   h_clusSize_nClus->GetYaxis()->SetTitle("Cluster Size");
0252   hm->registerHisto(h_clusSize_nClus);
0253
0254   if (m_chipInfo)
0255   {
0256     for (const auto &[layer, nstave] : m_layerStaveMap)
0257     {
0258       for (int stave = 0; stave < nstave; stave++)
0259       {
0260         //! 9 chips on each stave
0261         for (int chip = 0; chip < 9; chip++)
0262         {
0263           h_clusperchip[layer][stave][chip] = new TH2F(std::format("{}nclusperchip{}_{}_{}", getHistoPrefix(), layer, stave, chip).c_str(),
0264                                                        "MVTX clusters per chip", 2000, -2, 2, 2000, -1, 1);
0265           h_clusperchip[layer][stave][chip]->GetXaxis()->SetTitle("Local z [cm]");
0266           h_clusperchip[layer][stave][chip]->GetYaxis()->SetTitle("Local rphi [cm]");
0267           hm->registerHisto(h_clusperchip[layer][stave][chip]);
0268         }
0269       }
0270     }
0271   }
0272
0273   return;
0274 }