QA/Tpc/TpcLaserQA.cc

0001 #include "TpcLaserQA.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004
0005 #include <tpc/LaserEventInfo.h>
0006 #include <trackbase/LaserCluster.h>
0007 #include <trackbase/LaserClusterContainer.h>
0008 #include <trackbase/TpcDefs.h>
0009
0010 #include <fun4all/Fun4AllHistoManager.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>  // for PHWHERE
0016
0017 #include <TH1.h>
0018 #include <TH2.h>
0019
0020 #include <cassert>
0021 #include <cmath>
0022 #include <format>
0023 #include <iostream>  // for basic_ostream, operator<<
0024 #include <map>       // for _Rb_tree_const_iterator
0025 #include <string>
0026 #include <utility>  // for get
0027
0028 TpcLaserQA::TpcLaserQA(const std::string &name)
0029   : SubsysReco(name)
0030 {
0031 }
0032
0033 int TpcLaserQA::InitRun(PHCompositeNode * /*topNode*/)
0034 {
0035   createHistos();
0036
0037   auto *hm = QAHistManagerDef::getHistoManager();
0038   assert(hm);
0039
0040   m_nLaserEvents = dynamic_cast<TH1 *>(hm->getHisto(std::format("{}nLaserEvents", getHistoPrefix())));
0041
0042   for (int side = 0; side < 2; side++)
0043   {
0044     m_TPCWheel[side] = dynamic_cast<TH2 *>(hm->getHisto(std::format("{}TPCWheel_{}", getHistoPrefix(), (side == 1 ? "North" : "South"))));
0045
0046     m_nLaserClusters[side] = dynamic_cast<TH1 *>(hm->getHisto(std::format("{}nLaserClusters_{}", getHistoPrefix(), (side == 1 ? "North" : "South"))));
0047     m_saturation[side] = dynamic_cast<TH2 *>(hm->getHisto(std::format("{}saturation_{}", getHistoPrefix(), (side == 1 ? "North" : "South"))));
0048
0049     for (int sec = 0; sec < 12; sec++)
0050     {
0051       m_sample_R1[side][sec] = dynamic_cast<TH1 *>(hm->getHisto(std::format("{}sample_R1_{}_{}", getHistoPrefix(), (side == 1 ? "North" : "South"), sec)));
0052       m_sample_R2[side][sec] = dynamic_cast<TH1 *>(hm->getHisto(std::format("{}sample_R2_{}_{}", getHistoPrefix(), (side == 1 ? "North" : "South"), sec)));
0053       m_sample_R3[side][sec] = dynamic_cast<TH1 *>(hm->getHisto(std::format("{}sample_R3_{}_{}", getHistoPrefix(), (side == 1 ? "North" : "South"), sec)));
0054     }
0055   }
0056
0057   return Fun4AllReturnCodes::EVENT_OK;
0058 }
0059
0060 int TpcLaserQA::process_event(PHCompositeNode *topNode)
0061 {
0062   LaserEventInfo *lei = findNode::getClass<LaserEventInfo>(topNode, "LaserEventInfo");
0063   if (!lei)
0064   {
0065     std::cout << PHWHERE << "LaserEventInfo Node missing, abort." << std::endl;
0066     return Fun4AllReturnCodes::ABORTRUN;
0067   }
0068
0069   LaserClusterContainer *lcc = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0070   if (!lcc)
0071   {
0072     std::cout << PHWHERE << "LASER_CLUSTER Node missing, abort." << std::endl;
0073     return Fun4AllReturnCodes::ABORTRUN;
0074   }
0075
0076   if (!lei->isLaserEvent())
0077   {
0078     return Fun4AllReturnCodes::EVENT_OK;
0079   }
0080
0081   /*
0082   if(lcc->size() < 1000)
0083   {
0084     return Fun4AllReturnCodes::EVENT_OK;
0085   }
0086   */
0087
0088   m_nLaserEvents->Fill(1.0);
0089
0090   int nN = 0;
0091   int nS = 0;
0092
0093   auto clusrange = lcc->getClusters();
0094   for (auto cmitr = clusrange.first;
0095        cmitr != clusrange.second;
0096        ++cmitr)
0097   {
0098     const auto &[cmkey, cmclus_orig] = *cmitr;
0099     LaserCluster *cmclus = cmclus_orig;
0100     int side = TpcDefs::getSide(cmkey);
0101
0102     if (side)
0103     {
0104       nN++;
0105     }
0106     else
0107     {
0108       nS++;
0109     }
0110
0111     const unsigned int nhits = cmclus->getNhits();
0112     for (unsigned int i = 0; i < nhits; i++)
0113     {
0114       float layer = cmclus->getHitLayer(i);
0115       float hitAdc = cmclus->getHitAdc(i);
0116       float hitIT = cmclus->getHitIT(i);
0117
0118       double phi = std::atan2(cmclus->getHitY(i), cmclus->getHitX(i));
0119       if (phi < -M_PI / 12.) { phi += 2 * M_PI;
0120 }
0121
0122       int mod = -1;
0123       double RValue = -999;
0124       if (layer >= 7 && layer <= 22)
0125       {
0126         mod = 0;
0127         RValue = 0.4;
0128       }
0129       else if (layer >= 23 && layer <= 38)
0130       {
0131         mod = 1;
0132         RValue = 0.6;
0133       }
0134       else if (layer >= 39 && layer <= 54)
0135       {
0136         mod = 2;
0137         RValue = 0.8;
0138       }
0139
0140       int sector = -1;
0141       for (int sec = 0; sec < 12; sec++)
0142       {
0143         if (phi >= sec * M_PI / 6.0 - M_PI / 12 && phi < sec * M_PI / 6.0 + M_PI / 12)
0144         {
0145           sector = sec;
0146           break;
0147         }
0148       }
0149
0150       if (mod == -1) { continue;
0151 }
0152
0153       m_TPCWheel[side]->Fill(phi, RValue, hitAdc);
0154       if (hitAdc > 900) { m_saturation[side]->Fill(phi, RValue);
0155 }
0156
0157       if (hitAdc > 100)
0158       {
0159         if (mod == 0)
0160         {
0161           m_sample_R1[side][sector]->Fill(hitIT);
0162         }
0163         else if (mod == 1)
0164         {
0165           m_sample_R2[side][sector]->Fill(hitIT);
0166         }
0167         if (mod == 2)
0168         {
0169           m_sample_R3[side][sector]->Fill(hitIT);
0170         }
0171       }
0172     }
0173   }
0174
0175   m_nLaserClusters[0]->Fill(nS);
0176   m_nLaserClusters[1]->Fill(nN);
0177
0178   return Fun4AllReturnCodes::EVENT_OK;
0179 }
0180
0181 std::string TpcLaserQA::getHistoPrefix() const { return std::string("h_") + Name() + std::string("_"); }  // Define prefix to all histos in HistoManager
0182
0183 void TpcLaserQA::createHistos()
0184 {
0185   auto *hm = QAHistManagerDef::getHistoManager();
0186   assert(hm);
0187
0188   auto *h1 = new TH1F(std::format("{}nLaserEvents", getHistoPrefix()).c_str(), "Number of Laser Events", 1, 0.5, 1.5);
0189   hm->registerHisto(h1);
0190
0191   for (int side = 0; side < 2; side++)
0192   {
0193     auto *h = new TH2F(std::format("{}TPCWheel_{}", getHistoPrefix(), (side == 1 ? "North" : "South")).c_str(), std::format("Laser Hit ADC per event {}", (side == 1 ? "North" : "South")).c_str(), 12, -M_PI / 12., 23. * M_PI / 12., 4, rBinEdges);
0194     hm->registerHisto(h);
0195
0196     auto *h2 = new TH1F(std::format("{}nLaserClusters_{}", getHistoPrefix(), (side == 1 ? "North" : "South")).c_str(), std::format("Number of Laser Clusters per Event {}", (side == 1 ? "North" : "South")).c_str(), 81, -50, 8050);
0197     if (side == 1)
0198     {
0199       h2->SetLineColor(kRed);
0200     }
0201     else
0202     {
0203       h2->SetLineColor(kBlue);
0204     }
0205     hm->registerHisto(h2);
0206
0207     auto *h3 = new TH2F(std::format("{}saturation_{}", getHistoPrefix(), (side == 1 ? "North" : "South")).c_str(), std::format("Number of Saturated Hits per event {}", (side == 1 ? "North" : "South")).c_str(), 12, -M_PI / 12., 23. * M_PI / 12., 4, rBinEdges);
0208     hm->registerHisto(h3);
0209
0210     for (int sec = 0; sec < 12; sec++)
0211     {
0212       auto *h4 = new TH1F(std::format("{}sample_R1_{}_{}", getHistoPrefix(), (side == 1 ? "North" : "South"), sec).c_str(), std::format("Diffuse Laser Time Sample Sector {} R1 {}", sec, (side == 1 ? "North" : "South")).c_str(), 51, 305.5, 356.5);
0213       h4->SetLineColor(kRed);
0214       auto *h5 = new TH1F(std::format("{}sample_R2_{}_{}", getHistoPrefix(), (side == 1 ? "North" : "South"), sec).c_str(), std::format("Diffuse Laser Time Sample Sector {} R2 {}", sec, (side == 1 ? "North" : "South")).c_str(), 51, 305.5, 356.5);
0215       h5->SetLineColor(kBlue);
0216       auto *h6 = new TH1F(std::format("{}sample_R3_{}_{}", getHistoPrefix(), (side == 1 ? "North" : "South"), sec).c_str(), std::format("Diffuse Laser Time Sample Sector {} R3 {}", sec, (side == 1 ? "North" : "South")).c_str(), 51, 305.5, 356.5);
0217       h6->SetLineColor(kGreen + 2);
0218
0219       hm->registerHisto(h4);
0220       hm->registerHisto(h5);
0221       hm->registerHisto(h6);
0222     }
0223   }
0224 }