File indexing completed on 2025-08-06 08:18:49
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 <stddef.h> // for size_t
0022 #include <stdint.h> // for uint16_t, int32_t
0023 #include <algorithm> // for max, sort
0024 #include <cassert>
0025 #include <cmath>
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 std::cout << __PRETTY_FUNCTION__ << " : Found TpcRawHitContainer Node "
0067 << thisNode->getName() << std::endl;
0068
0069 TpcRawHitContainer *rawhitcont = (TpcRawHitContainer *) thisNode->getData();
0070 if (rawhitcont)
0071 {
0072 rawhitcont_vec.push_back(rawhitcont);
0073 }
0074 }
0075 }
0076 }
0077
0078 auto hm = QAHistManagerDef::getHistoManager();
0079 assert(hm);
0080
0081 for (int s = 0; s < 24; s++)
0082 {
0083 h_nhits_sectors[s] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s)).c_str()));
0084 h_nhits_sectors_fees[s] = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees").c_str()));
0085 h_nhits_sectors_laser[s] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_laser").c_str()));
0086 h_nhits_sectors_fees_laser[s] = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees_laser").c_str()));
0087 for (int r = 0; r < 3; r++)
0088 {
0089 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)).c_str()));
0090 h_adc[s][r] = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "adc_sec" + std::to_string(s) + "_R" + std::to_string(r)).c_str()));
0091 }
0092 }
0093
0094 h_xy_N = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "xyPos_North").c_str()));
0095 h_xy_S = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "xyPos_South").c_str()));
0096
0097 return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100
0101 int TpcRawHitQA::process_event(PHCompositeNode * )
0102 {
0103
0104 std::vector<int> sectors;
0105 std::vector<uint16_t> fees;
0106
0107 float nhit_sectors[24] = {0};
0108 float nhit_sectors_fees[24][26] = {{0}};
0109 float nhit_sectors_laser[24] = {0};
0110 float nhit_sectors_fees_laser[24][26] = {{0}};
0111
0112 unsigned int raw_hit_num = 0;
0113 for (TpcRawHitContainer *&rawhitcont : rawhitcont_vec)
0114 {
0115 raw_hit_num = rawhitcont->get_nhits();
0116 for (unsigned int i = 0; i < raw_hit_num; i++)
0117 {
0118 auto hit = rawhitcont->get_hit(i);
0119 uint16_t sam = hit->get_samples();
0120 int32_t packet_id = hit->get_packetid();
0121 int ep = (packet_id - 4000) % 10;
0122 int sector = (packet_id - 4000 - ep) / 10;
0123 uint16_t fee = hit->get_fee();
0124 int channel = hit->get_channel();
0125 int region = FEE_R[fee];
0126 int feeM = FEE_map[fee];
0127 if (FEE_R[fee] == 2)
0128 {
0129 feeM += 6;
0130 }
0131 if (FEE_R[fee] == 3)
0132 {
0133 feeM += 14;
0134 }
0135 double R = M.getR(feeM, channel);
0136 double phi = 0;
0137 if (sector < 12)
0138 {
0139 phi = M.getPhi(feeM, channel) + (sector) *M_PI / 6;
0140 }
0141 else if (sector >= 12)
0142 {
0143 phi = M.getPhi(feeM, channel) + (18 - sector) * M_PI / 6;
0144 }
0145
0146 float median = 60;
0147 float stdDev = 3;
0148 if (longPresample)
0149 {
0150 std::vector<int> values;
0151 values.reserve(sam);
0152
0153
0154 for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(hit->CreateAdcIterator());
0155 !adc_iterator->IsDone();
0156 adc_iterator->Next())
0157 {
0158 values.push_back((int) adc_iterator->CurrentAdc());
0159 }
0160 std::sort(values.begin(), values.end());
0161 size_t size = values.size();
0162 if (size % 2 == 0)
0163 {
0164 median = (values[size / 2 - 1] + values[size / 2]) / 2.0;
0165 }
0166 else
0167 {
0168 median = values[size / 2];
0169 }
0170 std::vector<int> selectedValues;
0171 for (int value : values)
0172 {
0173 if (value >= median - 40 && value <= median + 40)
0174 {
0175 selectedValues.push_back(value);
0176 }
0177 }
0178
0179 if (selectedValues.size() > 0)
0180 {
0181 float mean = 0.0;
0182 for (int value : selectedValues)
0183 {
0184 mean += value;
0185 }
0186 mean /= selectedValues.size();
0187 float sumSquares = 0.0;
0188 for (int value : selectedValues)
0189 {
0190 float diff = value - mean;
0191 sumSquares += std::pow(diff, 2);
0192 }
0193 float variance = sumSquares / selectedValues.size();
0194 stdDev = std::sqrt(variance);
0195 }
0196 }
0197
0198
0199
0200
0201 for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(hit->CreateAdcIterator());
0202 !adc_iterator->IsDone();
0203 adc_iterator->Next())
0204 {
0205 const uint16_t sampleN = adc_iterator->CurrentTimeBin();
0206 const uint16_t adc = adc_iterator->CurrentAdc();
0207 if (adc - median <= (std::max(5 * stdDev, (float) 20.)))
0208 {
0209 continue;
0210 }
0211
0212 if (sector < 12)
0213 {
0214 h_xy_N->Fill(R * cos(phi), R * sin(phi));
0215 }
0216 if (sector >= 12)
0217 {
0218 h_xy_S->Fill(R * cos(phi), R * sin(phi));
0219 }
0220 h_nhits_sam[sector][region - 1]->Fill(sampleN);
0221 h_adc[sector][region - 1]->Fill(adc - median);
0222 if (sampleN >= 400 && sampleN <= 430)
0223 {
0224 nhit_sectors_laser[sector]++;
0225 nhit_sectors_fees_laser[sector][fee]++;
0226 continue;
0227 }
0228 nhit_sectors[sector]++;
0229 nhit_sectors_fees[sector][fee]++;
0230 }
0231 }
0232 }
0233
0234
0235 if (raw_hit_num == 0)
0236 {
0237 return Fun4AllReturnCodes::EVENT_OK;
0238 }
0239
0240 for (int s = 0; s < 24; s++)
0241 {
0242 h_nhits_sectors[s]->Fill(nhit_sectors[s]);
0243 h_nhits_sectors_laser[s]->Fill(nhit_sectors_laser[s]);
0244 for (int f = 0; f < 26; f++)
0245 {
0246 h_nhits_sectors_fees[s]->Fill(f, nhit_sectors_fees[s][f]);
0247 h_nhits_sectors_fees_laser[s]->Fill(f, nhit_sectors_fees_laser[s][f]);
0248 }
0249 }
0250
0251 return Fun4AllReturnCodes::EVENT_OK;
0252 }
0253
0254
0255 int TpcRawHitQA::EndRun(const int )
0256 {
0257 auto hm = QAHistManagerDef::getHistoManager();
0258 assert(hm);
0259
0260 return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263
0264
0265 std::string TpcRawHitQA::getHistoPrefix() const { return std::string("h_") + Name() + std::string("_"); }
0266
0267 void TpcRawHitQA::createHistos()
0268 {
0269 auto hm = QAHistManagerDef::getHistoManager();
0270 assert(hm);
0271
0272 for (int s = 0; s < 24; s++)
0273 {
0274 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s)).c_str(),
0275 std::string("Number of Hits in Sector " + std::to_string(s) + ";Number of Hits/Event;Entries").c_str(), 100, 0, 30000));
0276 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees").c_str(),
0277 std::string("Sector " + std::to_string(s) + " Fee Hit Distribution;FEE;Number of Hits/Event").c_str(), 26, -0.5, 25.5, 100, 0, 3000));
0278 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_laser").c_str(),
0279 std::string("Laser Hits in Sector " + std::to_string(s) + ";Number of Hits/Event;Entries").c_str(), 100, 0, 1000));
0280 hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "nhits_sec" + std::to_string(s) + "_fees_laser").c_str(),
0281 std::string("Sector " + std::to_string(s) + " Fee Laser Hits;FEE;Number of Hits/Event").c_str(), 26, -0.5, 25.5, 100, 0, 500));
0282 for (int r = 0; r < 3; r++)
0283 {
0284 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nhits_sample_sec" + std::to_string(s) + "_R" + std::to_string(r)).c_str(),
0285 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));
0286 hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "adc_sec" + std::to_string(s) + "_R" + std::to_string(r)).c_str(),
0287 std::string("Sector " + std::to_string(s) + " ADC Distribution;ADC-pedestal [ADU];Entries").c_str(), 281, -100, 1024));
0288 }
0289 }
0290
0291 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));
0292 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));
0293 }
