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 = atan2(cmclus->getHitY(i), cmclus->getHitX(i));
0119       if (phi < -M_PI / 12.) phi += 2 * M_PI;
0120
0121       int mod = -1;
0122       double RValue = -999;
0123       if (layer >= 7 && layer <= 22)
0124       {
0125         mod = 0;
0126         RValue = 0.4;
0127       }
0128       else if (layer >= 23 && layer <= 38)
0129       {
0130         mod = 1;
0131         RValue = 0.6;
0132       }
0133       else if (layer >= 39 && layer <= 54)
0134       {
0135         mod = 2;
0136         RValue = 0.8;
0137       }
0138
0139       int sector = -1;
0140       for (int sec = 0; sec < 12; sec++)
0141       {
0142         if (phi >= sec * M_PI / 6.0 - M_PI / 12 && phi < sec * M_PI / 6.0 + M_PI / 12)
0143         {
0144           sector = sec;
0145           break;
0146         }
0147       }
0148
0149       if (mod == -1) continue;
0150
0151       m_TPCWheel[side]->Fill(phi, RValue, hitAdc);
0152       if (hitAdc > 900) m_saturation[side]->Fill(phi, RValue);
0153
0154       if (hitAdc > 100)
0155       {
0156         if (mod == 0)
0157         {
0158           m_sample_R1[side][sector]->Fill(hitIT);
0159         }
0160         else if (mod == 1)
0161         {
0162           m_sample_R2[side][sector]->Fill(hitIT);
0163         }
0164         if (mod == 2)
0165         {
0166           m_sample_R3[side][sector]->Fill(hitIT);
0167         }
0168       }
0169     }
0170   }
0171
0172   m_nLaserClusters[0]->Fill(nS);
0173   m_nLaserClusters[1]->Fill(nN);
0174
0175   return Fun4AllReturnCodes::EVENT_OK;
0176 }
0177
0178 std::string TpcLaserQA::getHistoPrefix() const { return std::string("h_") + Name() + std::string("_"); }  // Define prefix to all histos in HistoManager
0179
0180 void TpcLaserQA::createHistos()
0181 {
0182   auto hm = QAHistManagerDef::getHistoManager();
0183   assert(hm);
0184
0185   auto h1 = new TH1F(std::format("{}nLaserEvents", getHistoPrefix()).c_str(), "Number of Laser Events", 1, 0.5, 1.5);
0186   hm->registerHisto(h1);
0187
0188   for (int side = 0; side < 2; side++)
0189   {
0190     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);
0191     hm->registerHisto(h);
0192
0193     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);
0194     if (side == 1)
0195     {
0196       h2->SetLineColor(kRed);
0197     }
0198     else
0199     {
0200       h2->SetLineColor(kBlue);
0201     }
0202     hm->registerHisto(h2);
0203
0204     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);
0205     hm->registerHisto(h3);
0206
0207     for (int sec = 0; sec < 12; sec++)
0208     {
0209       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);
0210       h4->SetLineColor(kRed);
0211       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);
0212       h5->SetLineColor(kBlue);
0213       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);
0214       h6->SetLineColor(kGreen + 2);
0215
0216       hm->registerHisto(h4);
0217       hm->registerHisto(h5);
0218       hm->registerHisto(h6);
0219     }
0220   }
0221 }