QA/Tracking/MicromegasClusterQA.cc

0001 #include "MicromegasClusterQA.h"
0002
0003 #include <micromegas/CylinderGeomMicromegas.h>
0004 #include <micromegas/MicromegasDefs.h>
0005
0006 #include <g4detectors/PHG4CylinderGeomContainer.h>
0007
0008 #include <trackbase/ActsGeometry.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/TrkrHit.h>
0015 #include <trackbase/TrkrHitSet.h>
0016 #include <trackbase/TrkrHitSetContainer.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 namespace
0035 {
0036   //! range adaptor to be able to use range-based for loop
0037   template <class T>
0038   class range_adaptor
0039   {
0040    public:
0041     explicit range_adaptor(const T& range)
0042       : m_range(range)
0043     {
0044     }
0045     const typename T::first_type& begin() { return m_range.first; }
0046     const typename T::second_type& end() { return m_range.second; }
0047
0048    private:
0049     T m_range;
0050   };
0051
0052   constexpr int m_max_cluster_count = 10;
0053   constexpr int m_max_cluster_size = 15;
0054   constexpr double m_max_cluster_charge = 5e3;
0055
0056 }  // namespace
0057
0058 //____________________________________________________________________________..
0059 MicromegasClusterQA::MicromegasClusterQA(const std::string& name)
0060   : SubsysReco(name)
0061 {
0062 }
0063
0064 //____________________________________________________________________________..
0065 int MicromegasClusterQA::InitRun(PHCompositeNode* topNode)
0066 {
0067   // print configuration
0068   std::cout << "MicromegasClusterQA::InitRun - m_use_default_pedestal: " << m_use_default_pedestal << std::endl;
0069   std::cout << "MicromegasClusterQA::InitRun - m_default_pedestal: " << m_default_pedestal << std::endl;
0070   std::cout
0071       << "MicromegasClusterQA::InitRun -"
0072       << " m_calibration_filename: "
0073       << (m_calibration_filename.empty() ? "unspecified" : m_calibration_filename)
0074       << std::endl;
0075
0076   // read calibrations
0077   if (!m_calibration_filename.empty())
0078   {
0079     m_calibration_data.read(m_calibration_filename);
0080   }
0081
0082   // get geometry and keep track of tiles per layer
0083   auto* geomContainer = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0084   assert(geomContainer);
0085
0086   // get number of layers
0087   m_nlayers = geomContainer->get_NLayers();
0088
0089   // loop over layers
0090   bool first = true;
0091   const auto range = geomContainer->get_begin_end();
0092   for (const auto& [layer, layergeom] : range_adaptor(range))
0093   {
0094     auto* const layergeom_mm = dynamic_cast<CylinderGeomMicromegas*>(layergeom);
0095     const int ntiles = layergeom_mm->get_tiles_count();
0096     const auto segmentation = layergeom_mm->get_segmentation_type();
0097
0098     for (int tile = 0; tile < ntiles; ++tile)
0099     {
0100       // generate hitset key get detector name and save
0101       const auto hitsetkey = MicromegasDefs::genHitSetKey(layer, segmentation, tile);
0102       const auto detector_name = m_mapping.get_detname_sphenix_from_hitsetkey(hitsetkey);
0103       m_detector_names.push_back(detector_name);
0104     }
0105
0106     // keep track of first layer
0107     if (first)
0108     {
0109       first = false;
0110       m_firstlayer = layer;
0111     }
0112   }
0113
0114   create_histograms();
0115   return Fun4AllReturnCodes::EVENT_OK;
0116 }
0117
0118 //____________________________________________________________________________..
0119 int MicromegasClusterQA::process_event(PHCompositeNode* topNode)
0120 {
0121   // acts geometry
0122   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0123   assert(m_tGeometry);
0124
0125   // hitset container
0126   m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0127   assert(m_hitsetcontainer);
0128
0129   // cluster map
0130   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0131   assert(m_cluster_map);
0132
0133   // cluster hit association map
0134   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0135   assert(m_cluster_hit_map);
0136
0137   auto* hm = QAHistManagerDef::getHistoManager();
0138   assert(hm);
0139
0140   // keep track of how many good clusters per detector
0141   std::array<int, MicromegasDefs::m_nfee> cluster_count = {};
0142   std::array<int, MicromegasDefs::m_nfee> good_cluster_count = {};
0143
0144   // first loop over TPOT hitsets
0145   for (const auto& [hitsetkey, hitset] : range_adaptor(m_hitsetcontainer->getHitSets(TrkrDefs::TrkrId::micromegasId)))
0146   {
0147     // get detector name, layer and tile associated to this hitset key
0148     const int layer = TrkrDefs::getLayer(hitsetkey);
0149     const int tile = MicromegasDefs::getTileId(hitsetkey);
0150     //    const auto detector_name = m_mapping.get_detname_sphenix_from_hitsetkey(hitsetkey);
0151
0152     // detector id
0153     const int detid = tile + (MicromegasDefs::m_ntiles * (layer - m_firstlayer));
0154
0155     // get clusters
0156     const auto cluster_range = m_cluster_map->getClusters(hitsetkey);
0157     cluster_count[detid] = std::distance(cluster_range.first, cluster_range.second);
0158
0159     // loop over clusters
0160     for (const auto& [ckey, cluster] : range_adaptor(cluster_range))
0161     {
0162       // find associated hits
0163       const auto hit_range = m_cluster_hit_map->getHits(ckey);
0164
0165       // store cluster size and fill cluster size histogram
0166       const int cluster_size = std::distance(hit_range.first, hit_range.second);
0167       m_h_cluster_size->Fill(detid, cluster_size);
0168
0169       // calculate cluster charge
0170       double cluster_charge = 0;
0171       for (const auto& [ckey2, hitkey] : range_adaptor(hit_range))
0172       {
0173         // get strip
0174         const auto strip = MicromegasDefs::getStrip(hitkey);
0175
0176         // get associated hit
0177         auto* const hit = hitset->getHit(hitkey);
0178         assert(hit);
0179
0180         // get adc, remove pedestal, increment total charge
0181         const auto adc = hit->getAdc();
0182         const double pedestal = m_use_default_pedestal ? m_default_pedestal : m_calibration_data.get_pedestal_mapped(hitsetkey, strip);
0183         cluster_charge += (adc - pedestal);
0184       }
0185
0186       // fill cluster charge histogram
0187       m_h_cluster_charge->Fill(detid, cluster_charge);
0188
0189       // increment good clusters
0190       /*
0191        * we cut on cluster charge > 200 to define good clusters.
0192        * This is consistent with what has been done offline so far.
0193        * other cuts considered could include cluster size, cluster strip, etc.
0194        */
0195       if (cluster_charge > 200)
0196       {
0197         ++good_cluster_count[detid];
0198       }
0199     }
0200   }
0201
0202   // fill reference and found cluster histograms
0203   for (int layer = 0; layer < m_nlayers; ++layer)
0204   {
0205     for (int tile = 0; tile < MicromegasDefs::m_ntiles; ++tile)
0206     {
0207       // get detector id
0208       const int detid = tile + (MicromegasDefs::m_ntiles * layer);
0209
0210       // fill multiplicity histogram
0211       m_h_cluster_multiplicity->Fill(detid, cluster_count[detid]);
0212
0213       // get reference detector id. It corresponds to the same tile, but on the other layer
0214       const int detid_ref = detid >= MicromegasDefs::m_ntiles ? detid - MicromegasDefs::m_ntiles : detid + MicromegasDefs::m_ntiles;
0215
0216       // check if there is exactly one good cluster in the reference detector
0217       if (good_cluster_count[detid_ref] == 1 && cluster_count[detid_ref] == 1)
0218       {
0219         // fill curent detector ref count
0220         m_h_cluster_count_ref->Fill(detid);
0221
0222         // if there is one or more cluster in the current detector, also fill good count
0223         if (cluster_count[detid])
0224         {
0225           m_h_cluster_count_found->Fill(detid);
0226         }
0227       }
0228     }
0229   }
0230
0231   return Fun4AllReturnCodes::EVENT_OK;
0232 }
0233
0234 //____________________________________________________________________________..
0235 std::string MicromegasClusterQA::get_histogram_prefix() const
0236 {
0237   return std::string("h_") + Name() + std::string("_");
0238 }
0239
0240 //____________________________________________________________________________..
0241 void MicromegasClusterQA::create_histograms()
0242 {
0243   auto* hm = QAHistManagerDef::getHistoManager();
0244   assert(hm);
0245
0246   // cluster count histograms
0247   const int n_detectors = m_detector_names.size();
0248   m_h_cluster_count_ref = new TH1F(std::format("{}clustercount_ref", get_histogram_prefix()).c_str(), "reference cluster count", n_detectors, 0, n_detectors);
0249   m_h_cluster_count_found = new TH1F(std::format("{}clustercount_found", get_histogram_prefix()).c_str(), "found cluster count", n_detectors, 0, n_detectors);
0250
0251   // cluster multiplicity, size and charge distributions
0252   m_h_cluster_multiplicity = new TH2F(std::format("{}cluster_multiplicity", get_histogram_prefix()).c_str(), "cluster multiplicity", n_detectors, 0, n_detectors, m_max_cluster_count, 0, m_max_cluster_count);
0253   m_h_cluster_size = new TH2F(std::format("{}cluster_size", get_histogram_prefix()).c_str(), "cluster size", n_detectors, 0, n_detectors, m_max_cluster_size, 0, m_max_cluster_size);
0254   m_h_cluster_charge = new TH2F(std::format("{}cluster_charge", get_histogram_prefix()).c_str(), "cluster charge", n_detectors, 0, n_detectors, 100, 0, m_max_cluster_charge);
0255
0256   // assign bin labels
0257   for (int i = 0; i < n_detectors; ++i)
0258   {
0259     for (TH1* h : std::vector<TH1*>({m_h_cluster_count_ref, m_h_cluster_count_found, m_h_cluster_multiplicity, m_h_cluster_size, m_h_cluster_charge}))
0260     {
0261       h->GetXaxis()->SetBinLabel(i + 1, m_detector_names[i].c_str());
0262     }
0263   }
0264
0265   // register
0266   for (TH1* h : std::vector<TH1*>({m_h_cluster_count_ref, m_h_cluster_count_found, m_h_cluster_multiplicity, m_h_cluster_size, m_h_cluster_charge}))
0267   {
0268     hm->registerHisto(h);
0269   }
0270
0271   return;
0272 }