File indexing completed on 2025-12-17 09:21:25
0001 #include "TpcRawHitQA.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004
0005 #include <fun4all/Fun4AllHistoManager.h>
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/Fun4AllServer.h>
0008
0009 #include <ffarawobjects/TpcRawHit.h>
0010 #include <ffarawobjects/TpcRawHitContainer.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h> // for PHIODataNode
0014 #include <phool/PHNode.h> // for PHNode
0015 #include <phool/PHNodeIterator.h> // for PHNodeIterator
0016 #include <phool/PHPointerListIterator.h>
0017
0018 #include <TH1.h>
0019 #include <TH2.h>
0020
0021 #include <algorithm> // for max, sort
0022 #include <cassert>
0023 #include <cmath>
0024 #include <cstddef> // for size_t
0025 #include <cstdint> // for uint16_t, int32_t
0026 #include <iostream> // for basic_ostream, operator<<
0027 #include <memory> // for unique_ptr
0028
0029
0030 TpcRawHitQA::TpcRawHitQA(const std::string &name)
0031 : SubsysReco(name)
0032 {
0033 M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
0034 }
0035
0036
0037 int TpcRawHitQA::InitRun(PHCompositeNode *topNode)
0038 {
0039 createHistos();
0040
0041 PHNodeIterator trkr_itr(topNode);
0042 PHCompositeNode *tpc_node = dynamic_cast<PHCompositeNode *>(trkr_itr.findFirst("PHCompositeNode", "TPC"));
0043 if (!tpc_node)
0044 {
0045 std::cout << __PRETTY_FUNCTION__ << " : ERROR : "
0046 << "No TPC node found, unregistering module" << std::endl;
0047 Fun4AllServer::instance()->unregisterSubsystem(this);
0048 return Fun4AllReturnCodes::EVENT_OK;
0049 }
0050
0051 PHNodeIterator tpc_itr(tpc_node);
0052 {
0053 PHPointerListIterator<PHNode> iter(tpc_itr.ls());
0054
0055 PHNode *thisNode_raw;
0056 while ((thisNode_raw = iter()))
0057 {
0058 if (thisNode_raw->getType() != "PHIODataNode")
0059 {
0060 continue;
0061 }
0062
0063 PHIODataNode<TpcRawHitContainer> *thisNode = static_cast<PHIODataNode<TpcRawHitContainer> *>(thisNode_raw);
0064 if (thisNode)
0065 {
0066 if (thisNode->getName() == "HEADER" || thisNode->getName().find("G4") != std::string::npos)
0067 {
0068 continue;
0069 }
0070 std::cout << __PRETTY_FUNCTION__ << " : Found TpcRawHitContainer Node "
0071 << thisNode->getName() << std::endl;
0072
0073 TpcRawHitContainer *rawhitcont = (TpcRawHitContainer *) thisNode->getData();
0074 if (rawhitcont)
0075 {
0076 rawhitcont_vec.push_back(rawhitcont);
0077 }
0078 }
0079 }
0080 }
0081
0082 auto *hm = QAHistManagerDef::getHistoManager();
0083 assert(hm);
0084
0085 for (int s = 0; s < 24; s++)
0086 {
0087 h_nhits_sectors[s] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s))));
0088 h_nhits_sectors_fees[s] = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees")));
0089 h_nhits_sectors_laser[s] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_laser")));
0090 h_nhits_sectors_fees_laser[s] = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees_laser")));
0091 for (int f = 0; f < 26; f++)
0092 {
0093 h_nhits_sectors_fees_sampas[s][f] = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees" + std::to_string(f) + "_sampas")));
0094 }
0095 for (int r = 0; r < 3; r++)
0096 {
0097 h_nhits_sam[s][r] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sample_sec" + std::to_string(s) + "_R" + std::to_string(r))));
0098 h_adc[s][r] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "adc_sec" + std::to_string(s) + "_R" + std::to_string(r))));
0099 }
0100 }
0101
0102 h_xy_N = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "xyPos_North")));
0103 h_xy_S = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "xyPos_South")));
0104
0105 return Fun4AllReturnCodes::EVENT_OK;
0106 }
0107
0108
0109 int TpcRawHitQA::process_event(PHCompositeNode * )
0110 {
0111
0112 std::vector<int> sectors;
0113 std::vector<uint16_t> fees;
0114 std::vector<uint16_t> sampas;
0115
0116 float nhit_sectors[24] = {0};
0117 float nhit_sectors_fees[24][26] = {{0}};
0118 float nhit_sectors_laser[24] = {0};
0119 float nhit_sectors_fees_laser[24][26] = {{0}};
0120 float nhit_sectors_fees_sampas[24][26][8] = {{{0}}};
0121
0122 unsigned int raw_hit_num = 0;
0123 for (TpcRawHitContainer *&rawhitcont : rawhitcont_vec)
0124 {
0125 raw_hit_num = rawhitcont->get_nhits();
0126 for (unsigned int i = 0; i < raw_hit_num; i++)
0127 {
0128 auto *hit = rawhitcont->get_hit(i);
0129 uint16_t sam = hit->get_samples();
0130 int32_t packet_id = hit->get_packetid();
0131 int ep = (packet_id - 4000) % 10;
0132 int sector = (packet_id - 4000 - ep) / 10;
0133 uint16_t fee = hit->get_fee();
0134 uint16_t sampa = hit->get_sampaaddress();
0135 int channel = hit->get_channel();
0136 int region = FEE_R[fee];
0137 int feeM = FEE_map[fee];
0138 if (FEE_R[fee] == 2)
0139 {
0140 feeM += 6;
0141 }
0142 if (FEE_R[fee] == 3)
0143 {
0144 feeM += 14;
0145 }
0146 double R = M.getR(feeM, channel);
0147 double phi = 0;
0148 if (sector < 12)
0149 {
0150 phi = M.getPhi(feeM, channel) + (sector) *M_PI / 6;
0151 }
0152 else if (sector >= 12)
0153 {
0154 phi = M.getPhi(feeM, channel) + (18 - sector) * M_PI / 6;
0155 }
0156
0157 float median = 60;
0158 float stdDev = 3;
0159 if (longPresample)
0160 {
0161 std::vector<int> values;
0162 values.reserve(sam);
0163
0164
0165 for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(hit->CreateAdcIterator());
0166 !adc_iterator->IsDone();
0167 adc_iterator->Next())
0168 {
0169 values.push_back((int) adc_iterator->CurrentAdc());
0170 }
0171 std::sort(values.begin(), values.end());
0172 size_t size = values.size();
0173 if (size % 2 == 0)
0174 {
0175 median = (values[size / 2 - 1] + values[size / 2]) / 2.0;
0176 }
0177 else
0178 {
0179 median = values[size / 2];
0180 }
0181 std::vector<int> selectedValues;
0182 for (int value : values)
0183 {
0184 if (value >= median - 40 && value <= median + 40)
0185 {
0186 selectedValues.push_back(value);
0187 }
0188 }
0189
0190 if (!selectedValues.empty())
0191 {
0192 float mean = 0.0;
0193 for (int value : selectedValues)
0194 {
0195 mean += value;
0196 }
0197 mean /= selectedValues.size();
0198 float sumSquares = 0.0;
0199 for (int value : selectedValues)
0200 {
0201 float diff = value - mean;
0202 sumSquares += std::pow(diff, 2);
0203 }
0204 float variance = sumSquares / selectedValues.size();
0205 stdDev = std::sqrt(variance);
0206 }
0207 }
0208
0209
0210
0211
0212 for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(hit->CreateAdcIterator());
0213 !adc_iterator->IsDone();
0214 adc_iterator->Next())
0215 {
0216 const uint16_t sampleN = adc_iterator->CurrentTimeBin();
0217 const uint16_t adc = adc_iterator->CurrentAdc();
0218 if (adc - median <= (std::max(5 * stdDev, (float) 20.)))
0219 {
0220 continue;
0221 }
0222
0223 if (sector < 12)
0224 {
0225 h_xy_N->Fill(R * cos(phi), R * sin(phi));
0226 }
0227 if (sector >= 12)
0228 {
0229 h_xy_S->Fill(R * cos(phi), R * sin(phi));
0230 }
0231 h_nhits_sam[sector][region - 1]->Fill(sampleN);
0232 h_adc[sector][region - 1]->Fill(adc - median);
0233 if (sampleN >= 400 && sampleN <= 430)
0234 {
0235 nhit_sectors_laser[sector]++;
0236 nhit_sectors_fees_laser[sector][fee]++;
0237 continue;
0238 }
0239 nhit_sectors[sector]++;
0240 nhit_sectors_fees[sector][fee]++;
0241 nhit_sectors_fees_sampas[sector][fee][sampa]++;
0242 }
0243 }
0244 }
0245
0246
0247 if (raw_hit_num == 0)
0248 {
0249 return Fun4AllReturnCodes::EVENT_OK;
0250 }
0251
0252 for (int s = 0; s < 24; s++)
0253 {
0254 h_nhits_sectors[s]->Fill(nhit_sectors[s]);
0255 h_nhits_sectors_laser[s]->Fill(nhit_sectors_laser[s]);
0256 for (int f = 0; f < 26; f++)
0257 {
0258 h_nhits_sectors_fees[s]->Fill(f, nhit_sectors_fees[s][f]);
0259 h_nhits_sectors_fees_laser[s]->Fill(f, nhit_sectors_fees_laser[s][f]);
0260 for (int sm = 0; sm < 8; sm++)
0261 {
0262 h_nhits_sectors_fees_sampas[s][f]->Fill(sm, nhit_sectors_fees_sampas[s][f][sm]);
0263 }
0264 }
0265 }
0266
0267 return Fun4AllReturnCodes::EVENT_OK;
0268 }
0269
0270
0271 int TpcRawHitQA::EndRun(const int )
0272 {
0273 auto *hm = QAHistManagerDef::getHistoManager();
0274 assert(hm);
0275
0276 return Fun4AllReturnCodes::EVENT_OK;
0277 }
0278
0279
0280
0281 std::string TpcRawHitQA::getHistoPrefix() const { return std::string("h_") + Name() + std::string("_"); }
0282
0283 void TpcRawHitQA::createHistos()
0284 {
0285 auto *hm = QAHistManagerDef::getHistoManager();
0286 assert(hm);
0287
0288 for (int s = 0; s < 24; s++)
0289 {
0290 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s)).c_str(),
0291 std::string("Number of Hits in Sector " + std::to_string(s) + ";Number of Hits/Event;Entries").c_str(), 100, 0, 30000));
0292 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees").c_str(),
0293 std::string("Sector " + std::to_string(s) + " Fee Hit Distribution;FEE;Number of Hits/Event").c_str(), 26, -0.5, 25.5, 101, -15, 3015));
0294 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_laser").c_str(),
0295 std::string("Laser Hits in Sector " + std::to_string(s) + ";Number of Hits/Event;Entries").c_str(), 100, 0, 1000));
0296 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees_laser").c_str(),
0297 std::string("Sector " + std::to_string(s) + " Fee Laser Hits;FEE;Number of Hits/Event").c_str(), 26, -0.5, 25.5, 101, -2.5, 502.5));
0298 for (int f = 0; f < 26; f++)
0299 {
0300 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees" + std::to_string(f) + "_sampas").c_str(),
0301 std::string("Sector " + std::to_string(s) + " Fee " + std::to_string(f) + " Sampa Hit Distribution;Sampa;Number of Hits/Event").c_str(), 8, -0.5, 7.5, 101, -15, 3015));
0302 }
0303 for (int r = 0; r < 3; r++)
0304 {
0305 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nhits_sample_sec" + std::to_string(s) + "_R" + std::to_string(r)).c_str(),
0306 std::string("Sector " + std::to_string(s) + " Sample Time Distribution;Time Bin [1/20 MHz];Number of Total Hits").c_str(), 1051, -0.5, 1050.5));
0307 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "adc_sec" + std::to_string(s) + "_R" + std::to_string(r)).c_str(),
0308 std::string("Sector " + std::to_string(s) + " ADC Distribution;ADC-pedestal [ADU];Entries").c_str(), 281, -100, 1024));
0309 }
0310 }
0311
0312 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "xyPos_North").c_str(), "Hit XY distribution (North);X [mm];Y [mm]", 400, -800, 800, 400, -800, 800));
0313 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "xyPos_South").c_str(), "Hit XY distribution (South);X [mm];Y [mm]", 400, -800, 800, 400, -800, 800));
0314 }
