QA/Tpc/TpcNoiseQA.cc

0001 // Includes
0002 #include "TpcNoiseQA.h"
0003
0004 #include <fun4all/Fun4AllHistoManager.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009
0010 #include <Event/Event.h>
0011 #include <Event/packet.h>
0012
0013 #include <TH2.h>
0014
0015 #include <qautils/QAHistManagerDef.h>
0016
0017 #include <cassert>
0018 #include <cstdint>
0019 #include <format>
0020 #include <iostream>
0021 #include <string>
0022 //
0023
0024 //____________________________________________________________________________..
0025 TpcNoiseQA::TpcNoiseQA(const std::string &name)
0026   : SubsysReco(name)
0027 {
0028   // reserves memory for max ADC samples
0029   m_adcSamples.resize(1024, 0);
0030
0031   M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
0032 }
0033
0034 //____________________________________________________________________________..
0035 int TpcNoiseQA::InitRun(PHCompositeNode * /*unused*/)
0036 {
0037   // // Creates data file and checks whether it was successfully opened
0038   // m_file = TFile::Open(m_fname.c_str(), "recreate");
0039   // assert(m_file->IsOpen());
0040
0041   // double r_bins_new[r_bins_N + 1];
0042   for (int i = 0; i < r_bins_N + 1; i++)
0043   {
0044     r_bins_new[i] = -1 * r_bins[r_bins_N - i];
0045     r_bins_new[i + r_bins_N + nphi + 1] = r_bins[i];
0046   }
0047   for (int i = 0; i < nphi + 1; i++)
0048   {
0049     r_bins_new[i + r_bins_N] = phi_bins[i];
0050   }
0051
0052   createHistos();
0053
0054   return Fun4AllReturnCodes::EVENT_OK;
0055 }
0056
0057 //____________________________________________________________________________..
0058 int TpcNoiseQA::process_event(PHCompositeNode *topNode)
0059 {
0060   // Defines object from class Event which calls getClass function from
0061   // findNode class
0062   Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0063
0064   // Checks if event exists and returns error if not
0065   if (_event == nullptr)
0066   {
0067     std::cout << "TPCRawDataTree::Process_Event - Event not found" << std::endl;
0068     return -1;
0069   }
0070
0071   // Checks if event is "special" and discards it if so
0072   if (_event->getEvtType() >= 8)  /// special events
0073   {
0074     return Fun4AllReturnCodes::DISCARDEVENT;
0075   }
0076
0077   for (int fee_no = 0; fee_no < 26; fee_no++)
0078   {
0079     for (int channel_no = 0; channel_no < 256; channel_no++)
0080     {
0081       ave_adc_fee_channel[fee_no][channel_no] = 0.0;
0082       std_adc_fee_channel[fee_no][channel_no] = 0.0;
0083       counts_adc_fee_channel[fee_no][channel_no] = 0.0;
0084     }
0085   }
0086
0087   // Call HistoManager
0088   auto hm = QAHistManagerDef::getHistoManager();
0089   assert(hm);
0090
0091   // Reference histograms initialized in header file to histos in HistoManager
0092   h_NPol_Ped_Mean = dynamic_cast<TH2F *>(hm->getHisto(std::format("{}NPol_Ped_Mean", getHistoPrefix())));
0093   h_NPol_Ped_RMS = dynamic_cast<TH2F *>(hm->getHisto(std::format("{}NPol_Ped_RMS", getHistoPrefix())));
0094   h_SPol_Ped_Mean = dynamic_cast<TH2F *>(hm->getHisto(std::format("{}SPol_Ped_Mean", getHistoPrefix())));
0095   h_SPol_Ped_RMS = dynamic_cast<TH2F *>(hm->getHisto(std::format("{}SPol_Ped_RMS", getHistoPrefix())));
0096   //
0097
0098   int sector = 0;
0099
0100   std::vector<Packet *> pktvec = _event->getPacketVector();
0101
0102   // Loop over packets in event
0103   for (auto packet : pktvec)
0104   {
0105     if (!packet)
0106     {
0107       if (Verbosity())
0108       {
0109         std::cout << __PRETTY_FUNCTION__ << " : missing packet " << std::endl;
0110       }
0111       continue;
0112     }
0113     int32_t packet_id = packet->getIdentifier();
0114
0115     if (Verbosity())
0116     {
0117       std::cout << __PRETTY_FUNCTION__ << " : decoding packet " << packet_id << std::endl;
0118     }
0119
0120     int ep = (packet_id - 4000) % 10;
0121     sector = (packet_id - 4000 - ep) / 10;
0122     if (sector > 11)
0123     {
0124       side = 1;
0125     }
0126     else
0127     {
0128       side = 0;
0129     }
0130
0131     // pull number of waveforms
0132     m_nWaveformInFrame = packet->iValue(0, "NR_WF");
0133
0134     for (int wf = 0; wf < m_nWaveformInFrame; wf++)
0135     {
0136       m_FEE = packet->iValue(wf, "FEE");
0137       m_Channel = packet->iValue(wf, "CHANNEL");
0138       m_nSamples = packet->iValue(wf, "SAMPLES");
0139
0140       // Checks if sample number and number of ADC values agrees
0141       // assert(m_nSamples < (int) m_adcSamples.size());
0142       if (m_nSamples > (int) m_adcSamples.size())
0143       {
0144         continue;
0145       }
0146
0147       dead = false;
0148       // Loop over samples in waveform
0149       for (int s = 0; s < m_nSamples; s++)
0150       {
0151         // Assign ADC value of sample s in waveform wf to adcSamples[s]
0152         m_adcSamples[s] = packet->iValue(wf, s);
0153
0154         if (m_adcSamples[s] == 0 || std::isnan(float(m_adcSamples[s])))
0155         {
0156           dead = true;
0157           break;
0158         }
0159       }
0160
0161       if (dead)
0162       {
0163         continue;
0164       }
0165
0166       for (int adc_sam_no = 0; adc_sam_no < m_nSamples; adc_sam_no++)
0167       {
0168         if (m_adcSamples[adc_sam_no] < 1024)
0169         {
0170           ave_adc_fee_channel[m_FEE][m_Channel] += m_adcSamples[adc_sam_no];
0171           std_adc_fee_channel[m_FEE][m_Channel] += pow(m_adcSamples[adc_sam_no], 2);
0172           counts_adc_fee_channel[m_FEE][m_Channel] += 1.0;
0173         }
0174       }
0175     }
0176   }
0177
0178   for (int fee_no = 0; fee_no < 26; fee_no++)
0179   {
0180     for (int channel_no = 0; channel_no < 256; channel_no++)
0181     {
0182       if (counts_adc_fee_channel[fee_no][channel_no] != 0.0)
0183       {
0184         temp1 = ave_adc_fee_channel[fee_no][channel_no] / counts_adc_fee_channel[fee_no][channel_no];
0185         temp2 = std_adc_fee_channel[fee_no][channel_no] / counts_adc_fee_channel[fee_no][channel_no];
0186         ave_adc_fee_channel[fee_no][channel_no] = temp1;
0187         std_adc_fee_channel[fee_no][channel_no] = temp2;
0188       }
0189       else
0190       {
0191         ave_adc_fee_channel[fee_no][channel_no] = 0.0;
0192         std_adc_fee_channel[fee_no][channel_no] = 0.0;
0193       }
0194
0195       // setting the mapp of the FEE
0196       feeM = FEE_map[fee_no];
0197       if (mod_arr[fee_no] == 2)
0198       {
0199         feeM += 6;
0200       }
0201       if (mod_arr[fee_no] == 3)
0202       {
0203         feeM += 14;
0204       }
0205       key = 256 * (feeM) + channel_no;
0206       R = M.getR(feeM, channel_no);
0207       phi = pow(-1, side) * M.getPhi(feeM, channel_no) + (sector - side * 12.0) * M_PI / 6.0;
0208
0209       if (phi < 0.0)
0210       {
0211         phi = phi + 2.0 * M_PI;
0212       }
0213
0214       pedMean = ave_adc_fee_channel[fee_no][channel_no];
0215       pedStdi = sqrt(std_adc_fee_channel[fee_no][channel_no] - pow(ave_adc_fee_channel[fee_no][channel_no], 2));
0216
0217       if (side == 0)
0218       {
0219         h_NPol_Ped_Mean->Fill(phi, R, pedMean);
0220         h_NPol_Ped_RMS->Fill(phi, R, pedStdi);
0221       }
0222
0223       else
0224       {
0225         h_SPol_Ped_Mean->Fill(phi, R, pedMean);
0226         h_SPol_Ped_RMS->Fill(phi, R, pedStdi);
0227       }
0228     }
0229   }
0230
0231   return Fun4AllReturnCodes::EVENT_OK;
0232 }
0233
0234 //____________________________________________________________________________..
0235 int TpcNoiseQA::End(PHCompositeNode * /*unused*/)
0236 {
0237   // // Set histogram directory to 0 so data is saved after closing file
0238   // h_NPol_Ped_Mean->SetDirectory(nullptr);
0239   // h_NPol_Ped_RMS->SetDirectory(nullptr);
0240   // h_SPol_Ped_Mean->SetDirectory(nullptr);
0241   // h_SPol_Ped_RMS->SetDirectory(nullptr);
0242
0243   // // Write histograms to file
0244   // m_file->cd();
0245   // h_NPol_Ped_Mean->Write();
0246   // h_NPol_Ped_RMS->Write();
0247   // h_SPol_Ped_Mean->Write();
0248   // h_SPol_Ped_RMS->Write();
0249
0250   // std::cout << __PRETTY_FUNCTION__ << " : completed saving to " << m_file->GetName() << std::endl;
0251
0252   // m_file->ls();
0253
0254   // // Close the file
0255   // m_file->Close();
0256
0257   return Fun4AllReturnCodes::EVENT_OK;
0258 }
0259
0260 //____________________________________________________________________________..
0261 std::string TpcNoiseQA::getHistoPrefix() const { return std::string("h_") + Name() + std::string("_"); }  // Define prefix to all histos in HistoManager
0262
0263 //____________________________________________________________________________..
0264 void TpcNoiseQA::createHistos()
0265 {
0266   // Initialize HistoManager
0267   auto hm = QAHistManagerDef::getHistoManager();
0268   assert(hm);
0269
0270   // Create and register histos in HistoManager
0271   {
0272     auto h = new TH2F(std::format("{}NPol_Ped_Mean", getHistoPrefix()).c_str(), ";x;y", (2 * r_bins_N + nphi + 1), r_bins_new, (2 * r_bins_N + nphi + 1), r_bins_new);
0273     hm->registerHisto(h);
0274   }
0275
0276   {
0277     auto h = new TH2F(std::format("{}NPol_Ped_RMS", getHistoPrefix()).c_str(), ";x;y", (2 * r_bins_N + nphi + 1), r_bins_new, (2 * r_bins_N + nphi + 1), r_bins_new);
0278     hm->registerHisto(h);
0279   }
0280
0281   {
0282     auto h = new TH2F(std::format("{}SPol_Ped_Mean", getHistoPrefix()).c_str(), ";x;y", (2 * r_bins_N + nphi + 1), r_bins_new, (2 * r_bins_N + nphi + 1), r_bins_new);
0283     hm->registerHisto(h);
0284   }
0285
0286   {
0287     auto h = new TH2F(std::format("{}SPol_Ped_RMS", getHistoPrefix()).c_str(), ";x;y", (2 * r_bins_N + nphi + 1), r_bins_new, (2 * r_bins_N + nphi + 1), r_bins_new);
0288     hm->registerHisto(h);
0289   }
0290 }