packages/tpc/TpcRawDataTree.cc

0001
0002 #include "TpcRawDataTree.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005
0006 #include <phool/getClass.h>
0007
0008 #include <Event/Event.h>
0009 #include <Event/packet.h>
0010
0011 #include <TFile.h>
0012 #include <TH1.h>
0013 #include <TH2.h>
0014 #include <TTree.h>
0015
0016 #include <cassert>
0017 #include <iostream>
0018 #include <memory>
0019
0020 //____________________________________________________________________________..
0021 TpcRawDataTree::TpcRawDataTree(const std::string &name)
0022   : SubsysReco("TpcRawDataTree")
0023   , m_fname(name)
0024 {
0025   // reserve memory for max ADC samples
0026   m_adcSamples.resize(1024, 0);
0027   M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
0028 }
0029
0030 //____________________________________________________________________________..
0031 int TpcRawDataTree::InitRun(PHCompositeNode * /*unused*/)
0032 {
0033   sectorNum = m_fname;
0034   size_t pos = sectorNum.find("TPC_ebdc");
0035   sectorNum.erase(sectorNum.begin(), sectorNum.begin() + pos + 8);
0036   sectorNum.erase(sectorNum.begin() + 2, sectorNum.end());
0037   if (sectorNum.at(0) == '0')
0038   {
0039     sectorNum.erase(sectorNum.begin(), sectorNum.begin() + 1);
0040   }
0041   if (stoi(sectorNum) > 11)
0042   {
0043     side = 1;
0044   }
0045
0046   m_file = TFile::Open(m_fname.c_str(), "recreate");
0047   assert(m_file->IsOpen());
0048
0049   m_PacketTree = new TTree("PacketTree", "Each entry is one packet");
0050
0051   m_PacketTree->Branch("packet", &m_packet, "packet/I");
0052   m_PacketTree->Branch("frame", &m_frame, "frame/I");
0053   m_PacketTree->Branch("nWaveormInFrame", &m_nWaveormInFrame, "nWaveormInFrame/I");
0054   m_PacketTree->Branch("nTaggerInFrame", &m_nTaggerInFrame, "nTaggerInFrame/I");
0055   m_PacketTree->Branch("maxFEECount", &m_maxFEECount, "maxFEECount/I");
0056
0057   m_SampleTree = new TTree("SampleTree", "Each entry is one waveform");
0058
0059   m_SampleTree->Branch("packet", &m_packet, "packet/I");
0060   m_SampleTree->Branch("frame", &m_frame, "frame/I");
0061   m_SampleTree->Branch("nWaveormInFrame", &m_nWaveormInFrame, "nWaveormInFrame/I");
0062   m_SampleTree->Branch("maxFEECount", &m_maxFEECount, "maxFEECount/I");
0063   m_SampleTree->Branch("nSamples", &m_nSamples, "nSamples/I");
0064   m_SampleTree->Branch("adcSamples", &m_adcSamples[0], "adcSamples[nSamples]/s");
0065   m_SampleTree->Branch("fee", &m_fee, "fee/I");
0066   m_SampleTree->Branch("sampaAddress", &m_sampaAddress, "sampaAddress/I");
0067   m_SampleTree->Branch("sampaChannel", &m_sampaChannel, "sampaChannel/I");
0068   m_SampleTree->Branch("Channel", &m_Channel, "Channel/I");
0069   m_SampleTree->Branch("BCO", &m_BCO, "BCO/I");
0070   m_SampleTree->Branch("checksum", &m_checksum, "checksum/I");
0071   m_SampleTree->Branch("checksumError", &m_checksumError, "checksumError/I");
0072
0073   m_TaggerTree = new TTree("TaggerTree", "Each entry is one tagger for level 1 trigger or endat tag");
0074
0075   m_TaggerTree->Branch("packet", &m_packet, "packet/I");
0076   m_TaggerTree->Branch("frame", &m_frame, "frame/I");
0077   m_TaggerTree->Branch("tagger_type", &m_tagger_type, "tagger_type/s");
0078   m_TaggerTree->Branch("is_endat", &m_is_endat, "is_endat/b");
0079   m_TaggerTree->Branch("is_lvl1", &m_is_lvl1, "is_lvl1/b");
0080   m_TaggerTree->Branch("bco", &m_bco, "bco/l");
0081   m_TaggerTree->Branch("lvl1_count", &m_lvl1_count, "lvl1_count/i");
0082   m_TaggerTree->Branch("endat_count", &m_endat_count, "endat_count/i");
0083   m_TaggerTree->Branch("last_bco", &m_last_bco, "last_bco/l");
0084   m_TaggerTree->Branch("modebits", &m_modebits, "modebits/b");
0085
0086   R1_hist = new TH1F("R1_hist", "R1_hist", 1024, -0.5, 1023.5);
0087   R2_hist = new TH1F("R2_hist", "R2_hist", 1024, -0.5, 1023.5);
0088   R3_hist = new TH1F("R3_hist", "R3_hist", 1024, -0.5, 1023.5);
0089
0090   R1_time = new TH2F("R1_time", "R1_time", 360, -0.5, 359.5, 1024, -0.5, 1023.5);
0091   R2_time = new TH2F("R2_time", "R2_time", 360, -0.5, 359.5, 1024, -0.5, 1023.5);
0092   R3_time = new TH2F("R3_time", "R3_time", 360, -0.5, 359.5, 1024, -0.5, 1023.5);
0093
0094   TotalFEE = new TH1F("TotalFEE", "Total FEE", 26, -0.5, 25.5);
0095   TotalFEEsampa = new TH1F("TotalFEEsampa", "Total FEE + sampa", 26 * 8, -0.5, 25 * 8 - .5);
0096   TotalFRAME = new TH1F("TotalFRAME", "Total FRAME", 21, -0.5, 20.5);
0097
0098   checksumError_fee = new TH1F("FEEWithError", "FEE with Error", 26, -0.5, 25.5);
0099   checksumError_feesampa = new TH1F("FEEsampaWithError", "FEE*8+sampa with Error", 26 * 8, -0.5, 25 * 8 - .5);
0100   checksumError_frame = new TH1F("FRAMEWithError", "FRAME with Error", 21, -0.5, 20.5);
0101
0102   if (m_includeXYPos)
0103   {
0104     m_SampleTree->Branch("xPos", &m_xPos, "xPos/d");
0105     m_SampleTree->Branch("yPos", &m_yPos, "yPos/d");
0106   }
0107
0108   return Fun4AllReturnCodes::EVENT_OK;
0109 }
0110
0111 //____________________________________________________________________________..
0112 int TpcRawDataTree::process_event(PHCompositeNode *topNode)
0113 {
0114   Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0115   if (_event == nullptr)
0116   {
0117     std::cout << "TpcRawDataTree::Process_Event Event not found" << std::endl;
0118     return -1;
0119   }
0120   if (_event->getEvtType() >= 8)  /// special events
0121   {
0122     return Fun4AllReturnCodes::DISCARDEVENT;
0123   }
0124
0125   m_frame = _event->getEvtSequence();
0126
0127   for (int packet : m_packets)
0128   {
0129     if (Verbosity())
0130     {
0131       std::cout << __PRETTY_FUNCTION__ << " : decoding packet " << packet << std::endl;
0132     }
0133
0134     m_packet = packet;
0135
0136     std::unique_ptr<Packet> p(_event->getPacket(m_packet));
0137     if (!p)
0138     {
0139       if (Verbosity())
0140       {
0141         std::cout << __PRETTY_FUNCTION__ << " : missing packet " << packet << std::endl;
0142       }
0143
0144       continue;
0145     }
0146
0147     m_nWaveormInFrame = p->iValue(0, "NR_WF");
0148     m_nTaggerInFrame = p->lValue(0, "N_TAGGER");
0149     m_maxFEECount = p->iValue(0, "MAX_FEECOUNT");
0150
0151     for (int t = 0; t < m_nTaggerInFrame; t++)
0152     {
0153       /*uint16_t*/ m_tagger_type = (uint16_t) (p->lValue(t, "TAGGER_TYPE"));
0154       /*uint8_t*/ m_is_endat = (uint8_t) (p->lValue(t, "IS_ENDAT"));
0155       /*uint8_t*/ m_is_lvl1 = (uint8_t) (p->lValue(t, "IS_LEVEL1_TRIGGER"));
0156       /*uint64_t*/ m_bco = (uint64_t) (p->lValue(t, "BCO"));
0157       /*uint32_t*/ m_lvl1_count = (uint32_t) (p->lValue(t, "LEVEL1_COUNT"));
0158       /*uint32_t*/ m_endat_count = (uint32_t) (p->lValue(t, "ENDAT_COUNT"));
0159       /*uint64_t*/ m_last_bco = (uint64_t) (p->lValue(t, "LAST_BCO"));
0160       /*uint8_t*/ m_modebits = (uint8_t) (p->lValue(t, "MODEBITS"));
0161
0162       m_TaggerTree->Fill();
0163     }
0164
0165     for (int wf = 0; wf < m_nWaveormInFrame; wf++)
0166     {
0167       m_BCO = p->iValue(wf, "BCO");
0168       m_nSamples = p->iValue(wf, "SAMPLES");
0169       m_fee = p->iValue(wf, "FEE");
0170
0171       m_sampaAddress = p->iValue(wf, "SAMPAADDRESS");
0172       m_sampaChannel = p->iValue(wf, "SAMPACHANNEL");
0173       m_Channel = p->iValue(wf, "CHANNEL");
0174       m_checksum = p->iValue(wf, "CHECKSUM");
0175       m_checksumError = p->iValue(wf, "CHECKSUMERROR");
0176
0177       TH1 *fillHist;
0178       TH2 *fillHist2D;
0179
0180       if (m_fee == 2 ||
0181           m_fee == 4 ||
0182           m_fee == 3 ||
0183           m_fee == 13 ||
0184           m_fee == 17 ||
0185           m_fee == 16)
0186       {
0187         fillHist = R1_hist;
0188         fillHist2D = R1_time;
0189       }
0190       else if (m_fee == 11 ||
0191                m_fee == 12 ||
0192                m_fee == 19 ||
0193                m_fee == 18 ||
0194                m_fee == 01 ||
0195                m_fee == 00 ||
0196                m_fee == 16 ||
0197                m_fee == 15)
0198       {
0199         fillHist = R2_hist;
0200         fillHist2D = R2_time;
0201       }
0202       else
0203       {
0204         fillHist = R3_hist;
0205         fillHist2D = R3_time;
0206       }
0207
0208       assert(m_nSamples < (int) m_adcSamples.size());  // no need for movements in memory allocation
0209       for (int s = 0; s < m_nSamples; s++)
0210       {
0211         m_adcSamples[s] = p->iValue(wf, s);
0212         if (m_checksumError == 0)
0213         {
0214           fillHist->Fill(m_adcSamples[s]);
0215           fillHist2D->Fill(s, m_adcSamples[s]);
0216         }
0217         else
0218         {
0219           checksumError_fee->Fill(m_fee);
0220           checksumError_feesampa->Fill((m_fee * 8. + m_sampaAddress));
0221           checksumError_frame->Fill(m_frame);
0222         }
0223         TotalFEE->Fill(m_fee);
0224         TotalFEEsampa->Fill((m_fee * 8. + m_sampaAddress));
0225         TotalFRAME->Fill(m_frame);
0226       }
0227       if (m_includeXYPos)
0228       {
0229         int feeM = FEE_map[m_fee];
0230         if (FEE_R[m_fee] == 2)
0231         {
0232           feeM += 6;
0233         }
0234         if (FEE_R[m_fee] == 3)
0235         {
0236           feeM += 14;
0237         }
0238         int layer = M.getLayer(feeM, m_Channel);
0239         if (layer != 0)
0240         {
0241           double R = M.getR(feeM, m_Channel);
0242           double phi = M.getPhi(feeM, m_Channel) + (stod(sectorNum) - side * 12.) * M_PI / 6.;
0243           R /= 10.;  // convert mm to cm
0244
0245           m_xPos = R * cos(phi);
0246           m_yPos = R * sin(phi);
0247         }
0248         else
0249         {
0250           m_xPos = 0.;
0251           m_yPos = 0.;
0252         }
0253       }
0254       m_SampleTree->Fill();
0255     }
0256
0257     m_PacketTree->Fill();
0258   }  //   for (int packet : m_packets)
0259
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263 //____________________________________________________________________________..
0264 int TpcRawDataTree::End(PHCompositeNode * /*topNode*/)
0265 {
0266   checksumError_fee->Divide(TotalFEE);
0267   checksumError_feesampa->Divide(TotalFEEsampa);
0268   checksumError_frame->Divide(TotalFRAME);
0269
0270   TotalFEE->SetDirectory(nullptr);
0271   TotalFEEsampa->SetDirectory(nullptr);
0272   TotalFRAME->SetDirectory(nullptr);
0273
0274   m_file->Write();
0275
0276   std::cout << __PRETTY_FUNCTION__ << " : completed saving to " << m_file->GetName() << std::endl;
0277   m_file->ls();
0278
0279   m_file->Close();
0280
0281   return Fun4AllReturnCodes::EVENT_OK;
0282 }