framework/fun4allraw/SingleMicromegasPoolInput_v2.cc

0001 #include "SingleMicromegasPoolInput_v2.h"
0002
0003 #include "Fun4AllStreamingInputManager.h"
0004 #include "InputManagerType.h"
0005
0006 #include <ffarawobjects/MicromegasRawHitContainerv3.h>
0007 #include <ffarawobjects/MicromegasRawHitv3.h>
0008
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <qautils/QAHistManagerDef.h>
0011 #include <qautils/QAUtil.h>
0012
0013 #include <frog/FROG.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0017 #include <phool/getClass.h>
0018 #include <phool/phool.h>
0019
0020 #include <Event/Event.h>
0021 #include <Event/EventTypes.h>
0022 #include <Event/Eventiterator.h>
0023 #include <Event/fileEventiterator.h>
0024 #include <Event/oncsSubConstants.h>
0025
0026 #include <TFile.h>
0027 #include <TH1.h>
0028 #include <TTree.h>
0029
0030 #include <algorithm>
0031 #include <bitset>
0032
0033 namespace
0034 {
0035   // maximum number of packets
0036   static constexpr int m_npackets_active = 2;
0037
0038   //! number of sampa chips per fee board
0039   static constexpr int m_nsampa_fee = 8;
0040
0041   //! number of fee boards
0042   static constexpr int m_nfee_max = 26;
0043
0044   /* see: https://git.racf.bnl.gov/gitea/Instrumentation/sampa_data/src/branch/fmtv2/README.md */
0045   enum SampaDataType
0046   {
0047     HEARTBEAT_T = 0b000,
0048     TRUNCATED_DATA_T = 0b001,
0049     TRUNCATED_TRIG_EARLY_DATA_T = 0b011,
0050     NORMAL_DATA_T = 0b100,
0051     LARGE_DATA_T = 0b101,
0052     TRIG_EARLY_DATA_T = 0b110,
0053     TRIG_EARLY_LARGE_DATA_T = 0b111,
0054   };
0055
0056   /* see: https://git.racf.bnl.gov/gitea/Instrumentation/sampa_data/src/branch/fmtv2/README.md */
0057   enum ModeBitType
0058   {
0059     BX_COUNTER_SYNC_T = 0,
0060     ELINK_HEARTBEAT_T = 1,
0061     SAMPA_EVENT_TRIGGER_T = 2,
0062     CLEAR_LV1_LAST_T = 6,
0063     CLEAR_LV1_ENDAT_T = 7
0064   };
0065
0066   static constexpr uint16_t FEE_PACKET_MAGIC_KEY_1 = 0xfe;
0067   static constexpr uint16_t FEE_PACKET_MAGIC_KEY_2 = 0xed;
0068
0069   static constexpr uint16_t FEE_MAGIC_KEY = 0xba00;
0070   static constexpr uint16_t GTM_MAGIC_KEY = 0xbb00;
0071   static constexpr uint16_t GTM_LVL1_ACCEPT_MAGIC_KEY = 0xbbf0;
0072   static constexpr uint16_t GTM_ENDAT_MAGIC_KEY = 0xbbf1;
0073   static constexpr uint16_t GTM_MODEBIT_MAGIC_KEY = 0xbbf2;
0074
0075   // header length
0076   static constexpr uint16_t HEADER_LENGTH = 7;
0077
0078   // packet length
0079   static constexpr uint16_t MAX_PACKET_LENGTH = 1025;
0080
0081   //_____________________________________________________________
0082   [[maybe_unused]] uint16_t reverseBits(const uint16_t& x)
0083   {
0084     uint16_t n = x;
0085     n = ((uint16_t)(n >> 1U) & 0x55555555U) | ((uint16_t)(n << 1U) & 0xaaaaaaaaU);
0086     n = ((uint16_t)(n >> 2U) & 0x33333333U) | ((uint16_t)(n << 2U) & 0xccccccccU);
0087     n = ((uint16_t)(n >> 4U) & 0x0f0f0f0fU) | ((uint16_t)(n << 4U) & 0xf0f0f0f0U);
0088     n = ((uint16_t)(n >> 8U) & 0x00ff00ffU) | ((uint16_t)(n << 8U) & 0xff00ff00U);
0089     // n = (n >> 16U) & 0x0000ffffU | (n << 16U) & 0xffff0000U;
0090     return n;
0091   }
0092
0093   //_____________________________________________________________
0094   [[maybe_unused]] uint16_t crc16( const std::deque<uint16_t>& data, const unsigned int index, const int l)
0095   {
0096     uint16_t crc = 0xffffU;
0097
0098     for (int i = 0; i < l; i++)
0099     {
0100       uint16_t x = data[index + i];
0101       crc ^= reverseBits(x);
0102       for (uint16_t k = 0; k < 16; k++)
0103       {
0104         crc = crc & 1U ? (uint16_t)(crc >> 1U) ^ 0xa001U : crc >> 1U;
0105       }
0106     }
0107     crc = reverseBits(crc);
0108     return crc;
0109   }
0110
0111   /// used to make sure buffer is properly cleared when leaving current scope
0112   class buffer_cleaner_t
0113   {
0114     public:
0115
0116     buffer_cleaner_t( std::deque<uint16_t>& buffer, uint16_t pkt_length ):
0117       m_buffer( buffer ),
0118       m_pkt_length( pkt_length )
0119     {}
0120
0121     ~buffer_cleaner_t()
0122     { m_buffer.erase(m_buffer.begin(), m_buffer.begin() + m_pkt_length + 1); }
0123
0124     //! deleted copy constructor
0125     buffer_cleaner_t( const buffer_cleaner_t& ) = delete;
0126
0127     //! deleted copy constructor
0128     buffer_cleaner_t( buffer_cleaner_t&& ) = delete;
0129
0130     //! assignment operator
0131     buffer_cleaner_t& operator = ( const buffer_cleaner_t& ) = delete;
0132
0133     //! assignment operator
0134     buffer_cleaner_t& operator = ( buffer_cleaner_t&& ) = delete;
0135
0136     private:
0137
0138     std::deque<uint16_t>& m_buffer;
0139     uint16_t m_pkt_length;
0140   };
0141
0142 }  // namespace
0143
0144
0145 //______________________________________________________________
0146 SingleMicromegasPoolInput_v2::SingleMicromegasPoolInput_v2(const std::string& name)
0147   : SingleStreamingInput(name)
0148 {
0149   SubsystemEnum(InputManagerType::MICROMEGAS);
0150   m_rawHitContainerName = "MICROMEGASRAWHIT";
0151 }
0152
0153 //______________________________________________________________
0154 SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2()
0155 {
0156   std::cout << "SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2 - runnumber: " << RunNumber() << std::endl;
0157
0158   if (m_evaluation_file && m_evaluation_tree)
0159   {
0160     m_evaluation_file->cd();
0161     m_evaluation_tree->Write();
0162     m_evaluation_file->Close();
0163   }
0164
0165   // timer statistics
0166   m_timer.print_stat();
0167
0168   // dropped waveforms
0169   for( const auto& [packet,counter]:m_waveform_counters )
0170   {
0171     std::cout << "SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2 -"
0172       << " packet: " << packet
0173       << " wf_total: " << counter.total
0174       << " wf_dropped_bco: " << counter.dropped_bco
0175       << " wf_dropped_pool: " << counter.dropped_pool
0176       << " ratio_bco: " << counter.dropped_fraction_bco()
0177       << " ratio_pool: " << counter.dropped_fraction_pool()
0178       << std::endl;
0179   }
0180   std::cout << std::endl;
0181
0182   // dropped heartbeats
0183   for( const auto& [packet,counter]:m_heartbeat_counters )
0184   {
0185     std::cout << "SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2 -"
0186       << " packet: " << packet
0187       << " hb_total: " << counter.total
0188       << " hb_dropped_bco: " << counter.dropped_bco
0189       << " ratio_bco: " << counter.dropped_fraction_bco()
0190       << std::endl;
0191   }
0192   std::cout << std::endl;
0193
0194   // drop per fee statistics
0195   for( const auto& [fee,counter]:m_fee_waveform_counters )
0196   {
0197     std::cout << "SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2 -"
0198       << " fee: " << fee
0199       << " wf_total: " << counter.total
0200       << " wf_dropped_bco: " << counter.dropped_bco
0201       << " ratio_bco: " << counter.dropped_fraction_bco()
0202       << " ratio_pool: " << counter.dropped_fraction_pool()
0203       << std::endl;
0204   }
0205   std::cout << std::endl;
0206
0207   // drop per fee statistics
0208   for( const auto& [fee,counter]:m_fee_heartbeat_counters )
0209   {
0210     std::cout << "SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2 -"
0211       << " fee: " << fee
0212       << " hb_total: " << counter.total
0213       << " hb_dropped_bco: " << counter.dropped_bco
0214       << " ratio_bco: " << counter.dropped_fraction_bco()
0215       << std::endl;
0216   }
0217   std::cout << std::endl;
0218
0219   // also printout adjusted multipliers
0220   for( const auto& [packet_id, bco_matching]:m_bco_matching_information_map )
0221   {
0222     const auto old_precision = std::cout.precision();
0223     std::cout << "SingleMicromegasPoolInput_v2::~SingleMicromegasPoolInput_v2 -"
0224       << " packet: " << packet_id
0225       << " multiplier: " <<  std::setprecision(9) << bco_matching.get_adjusted_multiplier() << std::setprecision(old_precision)
0226       << std::endl;
0227   }
0228
0229 }
0230
0231 //______________________________________________________________
0232 void SingleMicromegasPoolInput_v2::FillPool(const unsigned int /*nbclks*/)
0233 {
0234   if (AllDone())  // no more files and all events read
0235   {
0236     return;
0237   }
0238
0239   while (!GetEventiterator())  // at startup this is a null pointer
0240   {
0241     if (!OpenNextFile())
0242     {
0243       AllDone(1);
0244       return;
0245     }
0246   }
0247
0248   while (GetSomeMoreEvents())
0249   {
0250     // std::cout << "SingleMicromegasPoolInput_v2::FillPool" << std::endl;
0251     std::unique_ptr<Event> evt(GetEventiterator()->getNextEvent());
0252     while (!evt)
0253     {
0254       fileclose();
0255       if (!OpenNextFile())
0256       {
0257         AllDone(1);
0258         return;
0259       }
0260
0261       // get next event
0262       evt.reset(GetEventiterator()->getNextEvent());
0263     }
0264
0265     if (Verbosity() > 2)
0266     {
0267       std::cout << "Fetching next Event" << evt->getEvtSequence() << std::endl;
0268     }
0269
0270     if (evt->getEvtType() != DATAEVENT)
0271     {
0272       m_NumSpecialEvents++;
0273       continue;
0274     }
0275
0276     RunNumber(evt->getRunNumber());
0277     if (Verbosity() > 1)
0278     {
0279       evt->identify();
0280     }
0281
0282     const int npackets = evt->getPacketList(&plist[0], 10);
0283
0284     if (npackets == 10)
0285     {
0286       std::cout << "SingleMicromegasPoolInput_v2::FillPool - too many packets" << std::endl;
0287       exit(1);
0288     }
0289
0290     m_timer.restart();
0291     for (int i = 0; i < npackets; i++)
0292     {
0293       // keep pointer to local packet
0294       std::unique_ptr<Packet> packet(plist[i]);
0295
0296       // process
0297       process_packet( packet.get() );
0298     }
0299
0300     if( Verbosity()>2)
0301     {
0302       for( size_t i=0; i<m_feeData.size(); ++i )
0303       {
0304         std::cout << " SingleMicromegasPoolInput_v2::FillPool -"
0305           << " fee: " << i
0306           << " buffer size: " << m_feeData[i].size()
0307           << std::endl;
0308       }
0309     }
0310
0311     m_timer.stop();
0312   }
0313 }
0314
0315 //______________________________________________________________
0316 void SingleMicromegasPoolInput_v2::Print(const std::string& what) const
0317 {
0318   if (what == "ALL" || what == "FEE")
0319   {
0320     for (const auto& bcliter : m_BeamClockFEE)
0321     {
0322       std::cout << "Beam clock 0x" << std::hex << bcliter.first << std::dec << std::endl;
0323       for (auto feeiter : bcliter.second)
0324       {
0325         std::cout << "FEM: " << feeiter << std::endl;
0326       }
0327     }
0328   }
0329
0330   if (what == "ALL" || what == "FEEBCLK")
0331   {
0332     for (auto bcliter : m_FEEBclkMap)
0333     {
0334       std::cout << "FEE" << bcliter.first << " bclk: 0x"
0335                 << std::hex << bcliter.second << std::dec << std::endl;
0336     }
0337   }
0338
0339   if (what == "ALL" || what == "STORAGE")
0340   {
0341     for (const auto& bcliter : m_MicromegasRawHitMap)
0342     {
0343       std::cout << "Beam clock 0x" << std::hex << bcliter.first << std::dec << std::endl;
0344       for (auto feeiter : bcliter.second)
0345       {
0346         std::cout
0347             << "fee: " << feeiter->get_fee()
0348             << " at " << std::hex << feeiter << std::dec
0349             << std::endl;
0350       }
0351     }
0352   }
0353
0354   if (what == "ALL" || what == "STACK")
0355   {
0356     for (const auto& iter : m_BclkStack)
0357     {
0358       std::cout << "stacked bclk: 0x" << std::hex << iter << std::dec << std::endl;
0359     }
0360   }
0361 }
0362
0363 //____________________________________________________________________________
0364 void SingleMicromegasPoolInput_v2::CleanupUsedPackets(const uint64_t bclk, bool dropped)
0365 {
0366
0367   // delete all raw hits associated to bco smaller than reference, and remove from map
0368   for(auto iter = m_MicromegasRawHitMap.begin(); iter != m_MicromegasRawHitMap.end() && (iter->first <= bclk); iter = m_MicromegasRawHitMap.erase(iter))
0369   {
0370     for (const auto& rawhit : iter->second)
0371     {
0372       if( dropped )
0373       {
0374         // increment dropped waveform counter and histogram
0375         ++m_waveform_counters[rawhit->get_packetid()].dropped_pool;
0376         ++m_fee_waveform_counters[rawhit->get_fee()].dropped_pool;
0377         h_waveform_count_dropped_pool->Fill( std::to_string(rawhit->get_packetid()).c_str(), 1 );
0378         h_fee_waveform_count_dropped_pool->Fill( rawhit->get_fee(), 1 );
0379       }
0380       delete rawhit;
0381     }
0382   }
0383
0384   // cleanup bco stacks
0385   /* it erases all elements for which the bco is no greater than the provided one */
0386   m_BclkStack.erase(m_BclkStack.begin(), m_BclkStack.upper_bound(bclk));
0387   m_BeamClockFEE.erase(m_BeamClockFEE.begin(), m_BeamClockFEE.upper_bound(bclk));
0388   m_BeamClockPacket.erase(m_BeamClockPacket.begin(), m_BeamClockPacket.upper_bound(bclk));
0389
0390   // cleanup matching information
0391   for( auto&& bco_matching:m_bco_matching_information_map )
0392   { bco_matching.second.cleanup(bclk); }
0393 }
0394
0395 //_______________________________________________________
0396 void SingleMicromegasPoolInput_v2::ClearCurrentEvent()
0397 {
0398   std::cout << "SingleMicromegasPoolInput_v2::ClearCurrentEvent." << std::endl;
0399   uint64_t currentbclk = *m_BclkStack.begin();
0400   CleanupUsedPackets(currentbclk);
0401   return;
0402 }
0403
0404 //_______________________________________________________
0405 bool SingleMicromegasPoolInput_v2::GetSomeMoreEvents()
0406 {
0407   if (AllDone())
0408   {
0409     return false;
0410   }
0411
0412   // check minimum pool size
0413   if (m_MicromegasRawHitMap.size() < m_BcoPoolSize)
0414   {
0415     return true;
0416   }
0417
0418   // make sure that the latest BCO received by each FEEs is past the current BCO
0419   std::set<int> toerase;
0420   uint64_t lowest_bclk = m_MicromegasRawHitMap.begin()->first + m_BcoRange;
0421   for (auto bcliter : m_FEEBclkMap)
0422   {
0423     if (bcliter.second <= lowest_bclk)
0424     {
0425       uint64_t highest_bclk = m_MicromegasRawHitMap.rbegin()->first;
0426       if ((highest_bclk - m_MicromegasRawHitMap.begin()->first) < MaxBclkDiff())
0427       {
0428         return true;
0429       }
0430       else
0431       {
0432         std::cout << PHWHERE << Name() << ": erasing FEE " << bcliter.first
0433                   << " with stuck bclk: " << std::hex << bcliter.second
0434                   << " current bco range: 0x" << m_MicromegasRawHitMap.begin()->first
0435                   << ", to: 0x" << highest_bclk << ", delta: " << std::dec
0436                   << (highest_bclk - m_MicromegasRawHitMap.begin()->first)
0437                   << std::dec << std::endl;
0438         toerase.insert(bcliter.first);
0439       }
0440     }
0441   }
0442   for(auto& fee: toerase)
0443   {
0444     m_FEEBclkMap.erase(fee);
0445   }
0446
0447   return false;
0448 }
0449
0450 //_______________________________________________________
0451 void SingleMicromegasPoolInput_v2::CreateDSTNode(PHCompositeNode* topNode)
0452 {
0453   PHNodeIterator iter(topNode);
0454   auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0455   if (!dstNode)
0456   {
0457     dstNode = new PHCompositeNode("DST");
0458     topNode->addNode(dstNode);
0459   }
0460
0461   PHNodeIterator iterDst(dstNode);
0462   auto detNode = dynamic_cast<PHCompositeNode*>(iterDst.findFirst("PHCompositeNode", "MICROMEGAS"));
0463   if (!detNode)
0464   {
0465     detNode = new PHCompositeNode("MICROMEGAS");
0466     dstNode->addNode(detNode);
0467   }
0468
0469   auto container = findNode::getClass<MicromegasRawHitContainer>(detNode, m_rawHitContainerName);
0470   if (!container)
0471   {
0472     container = new MicromegasRawHitContainerv3;
0473     auto newNode = new PHIODataNode<PHObject>(container, m_rawHitContainerName, "PHObject");
0474     detNode->addNode(newNode);
0475   }
0476 }
0477
0478 //_______________________________________________________
0479 void SingleMicromegasPoolInput_v2::ConfigureStreamingInputManager()
0480 {
0481   if (StreamingInputManager())
0482   {
0483     StreamingInputManager()->SetMicromegasBcoRange(m_BcoRange);
0484     StreamingInputManager()->SetMicromegasNegativeBco(m_NegativeBco);
0485   }
0486 }
0487
0488 //_______________________________________________________
0489 void SingleMicromegasPoolInput_v2::FillBcoQA(uint64_t gtm_bco)
0490 {
0491   auto hm = QAHistManagerDef::getHistoManager();
0492   assert(hm);
0493
0494   // set of packets for which the gtm BCO is found at least once
0495   unsigned int n_waveforms = 0;
0496   std::set<int> found_packets;
0497   for (uint64_t gtm_bco_loc = gtm_bco - m_NegativeBco; gtm_bco_loc < gtm_bco + m_BcoRange - m_NegativeBco; ++gtm_bco_loc)
0498   {
0499     // packet ids
0500     const auto packet_iter = m_BeamClockPacket.find(gtm_bco_loc);
0501     if( packet_iter != m_BeamClockPacket.end() )
0502     { found_packets.insert( packet_iter->second.begin(), packet_iter->second.end() ); }
0503
0504     // waveforms
0505     const auto wf_iter = m_MicromegasRawHitMap.find(gtm_bco_loc);
0506     if (wf_iter != m_MicromegasRawHitMap.end())
0507     { n_waveforms += wf_iter->second.size(); }
0508   }
0509
0510   // per packet statistics
0511   h_packet_stat->Fill( "Reference", 1 );
0512
0513   for( const auto& packet_id:found_packets )
0514   { h_packet_stat->Fill( std::to_string(packet_id).c_str(), 1 ); }
0515   h_packet_stat->Fill( "All", found_packets.size()>= m_npackets_active );
0516
0517   // how many packet_id found for this BCO
0518   h_packet->Fill(found_packets.size());
0519
0520   // how many waveforms found for this BCO
0521   h_waveform->Fill(n_waveforms);
0522 }
0523 //_______________________________________________________
0524 void SingleMicromegasPoolInput_v2::createQAHistos()
0525 {
0526   auto hm = QAHistManagerDef::getHistoManager();
0527   assert(hm);
0528
0529   // number of packets found with BCO from felix matching reference BCO
0530   h_packet = new TH1D( "h_MicromegasBCOQA_npacket_bco", "TPOT Packet Count per GTM BCO; Matching BCO tagger count; GL1 trigger count", 10, 0, 10 );
0531
0532   // number of waveforms found with BCO from felix matching reference BCO
0533   h_waveform = new TH1D( "h_MicromegasBCOQA_nwaveform_bco", "TPOT Waveform Count per GTM BCO; Matching Waveform count; GL1 trigger count", 4100, 0, 4100 );
0534
0535   /*
0536    * first bin is the number of requested GL1 BCO, for reference
0537    * next two bins is the number of times the GL1 BCO is found in the taggers list for a given packet_id
0538    * last bin is the sum
0539    */
0540   h_packet_stat = new TH1D( "h_MicromegasBCOQA_packet_stat", "Matching Tagger count per packet; packet id; GL1 trigger count", m_npackets_active+2, 0, m_npackets_active+2 );
0541   h_packet_stat->GetXaxis()->SetBinLabel(1, "Reference" );
0542   h_packet_stat->GetXaxis()->SetBinLabel(2, "5001" );
0543   h_packet_stat->GetXaxis()->SetBinLabel(3, "5002" );
0544   h_packet_stat->GetXaxis()->SetBinLabel(4, "All" );
0545
0546   /*
0547    * first two bins is the number of heartbeat MODEBITS are received from the GTM, per packet
0548    * next bins are how many heartbeat packets are recieved from the FEEs on a per sampa basis
0549    */
0550   h_heartbeat_stat = new TH1D(
0551     "h_MicromegasBCOQA_heartbeat_stat", "Heartbeat count per SAMPA; sampa id; heartbeat count",
0552     m_npackets_active+m_nfee_max*m_nsampa_fee, 0, m_npackets_active+m_nfee_max*m_nsampa_fee );
0553   h_heartbeat_stat->GetXaxis()->SetBinLabel(1, "5001" );
0554   h_heartbeat_stat->GetXaxis()->SetBinLabel(2, "5002" );
0555
0556   // total number of waveform per packet
0557   h_waveform_count_total = new TH1D( "h_MicromegasBCOQA_waveform_count_total", "Total number of waveforms per packet", m_npackets_active, 0, m_npackets_active );
0558
0559   // number of dropped waveform per packet due to bco mismatch
0560   h_waveform_count_dropped_bco = new TH1D( "h_MicromegasBCOQA_waveform_count_dropped_bco", "Number of dropped waveforms per packet (bco)", m_npackets_active, 0, m_npackets_active );
0561
0562   // number of dropped waveform per packet due to fun4all pool mismatch
0563   h_waveform_count_dropped_pool = new TH1D( "h_MicromegasBCOQA_waveform_count_dropped_pool", "Number of dropped waveforms per packet (pool)", m_npackets_active, 0, m_npackets_active );
0564
0565   // total number of waveform per packet
0566   h_fee_waveform_count_total = new TH1D( "h_MicromegasBCOQA_fee_waveform_count_total", "Total number of waveforms per fee", m_nfee_max, 0, m_nfee_max );
0567
0568   // number of dropped waveform per fee due to bco mismatch
0569   h_fee_waveform_count_dropped_bco = new TH1D( "h_MicromegasBCOQA_fee_waveform_count_dropped_bco", "Number of dropped waveforms per fee (bco)", m_nfee_max, 0, m_nfee_max );
0570
0571   // number of dropped waveform per fee due to fun4all pool mismatch
0572   h_fee_waveform_count_dropped_pool = new TH1D( "h_MicromegasBCOQA_fee_waveform_count_dropped_pool", "Number of dropped waveforms per fee (pool)", m_nfee_max, 0, m_nfee_max );
0573
0574   // define axis
0575   for( const auto& h:std::initializer_list<TH1*>{h_waveform_count_total, h_waveform_count_dropped_bco, h_waveform_count_dropped_pool} )
0576   {
0577     h->GetXaxis()->SetBinLabel(1, "5001" );
0578     h->GetXaxis()->SetBinLabel(2, "5002" );
0579     h->GetYaxis()->SetTitle( "waveform count" );
0580     h->GetXaxis()->SetTitle( "packet id" );
0581   }
0582
0583   for( const auto& h:std::initializer_list<TH1*>{h_fee_waveform_count_total, h_fee_waveform_count_dropped_bco, h_fee_waveform_count_dropped_pool})
0584   {
0585     h->GetYaxis()->SetTitle( "waveform count" );
0586     h->GetXaxis()->SetTitle( "FEE id" );
0587   }
0588
0589   // register all histograms to histogram manager
0590   for( const auto& h:std::initializer_list<TH1*>{
0591     h_packet, h_waveform, h_packet_stat, h_heartbeat_stat,
0592     h_waveform_count_total, h_waveform_count_dropped_bco, h_waveform_count_dropped_pool,
0593     h_fee_waveform_count_total, h_fee_waveform_count_dropped_bco, h_fee_waveform_count_dropped_pool} )
0594   {
0595     h->SetFillStyle(1001);
0596     h->SetFillColor(kYellow);
0597     hm->registerHisto(h);
0598   }
0599
0600   // also create evaluation trees
0601   if( m_do_evaluation )
0602   {
0603     m_evaluation_file.reset(new TFile(m_evaluation_filename.c_str(), "RECREATE"));
0604     m_evaluation_tree = new TTree("T", "T");
0605     m_evaluation_tree->Branch("is_heartbeat", &m_waveform.is_heartbeat);
0606     m_evaluation_tree->Branch("packet_id", &m_waveform.packet_id);
0607     m_evaluation_tree->Branch("fee_id", &m_waveform.fee_id);
0608     m_evaluation_tree->Branch("channel", &m_waveform.channel);
0609
0610     m_evaluation_tree->Branch("gtm_bco_first", &m_waveform.gtm_bco_first);
0611     m_evaluation_tree->Branch("gtm_bco", &m_waveform.gtm_bco);
0612     m_evaluation_tree->Branch("gtm_bco_matched", &m_waveform.gtm_bco_matched);
0613
0614     m_evaluation_tree->Branch("fee_bco_first", &m_waveform.fee_bco_first);
0615     m_evaluation_tree->Branch("fee_bco", &m_waveform.fee_bco);
0616     m_evaluation_tree->Branch("fee_bco_predicted", &m_waveform.fee_bco_predicted);
0617     m_evaluation_tree->Branch("fee_bco_predicted_matched", &m_waveform.fee_bco_predicted_matched);
0618   }
0619
0620 }
0621
0622 //__________________________________________________________________________________
0623 void SingleMicromegasPoolInput_v2::process_packet(Packet* packet )
0624 {
0625   // check hit format
0626   const auto hitformat = packet->getHitFormat();
0627   if (hitformat != IDTPCFEEV4 && hitformat != IDTPCFEEV5)
0628   { return; }
0629
0630   // set clock multipliers if not alsready
0631   if( !MicromegasBcoMatchingInformation_v2::gtm_clock_multiplier_is_set() )
0632   {
0633     switch( hitformat )
0634     {
0635       case IDTPCFEEV4:
0636       std::cout << "SingleMicromegasPoolInput_v2::process_packet - hitformat: IDTPCFEEV4" << std::endl;
0637       MicromegasBcoMatchingInformation_v2::set_gtm_clock_multiplier(4.262916255);
0638       break;
0639
0640       case IDTPCFEEV5:
0641       std::cout << "SingleMicromegasPoolInput_v2::process_packet - hitformat: IDTPCFEEV5" << std::endl;
0642       MicromegasBcoMatchingInformation_v2::set_enable_multiplier_adjustment(false);
0643       MicromegasBcoMatchingInformation_v2::set_gtm_clock_multiplier(30./8);
0644       break;
0645
0646       default: break;
0647     }
0648   }
0649
0650   // get packet id
0651   const int packet_id = packet->getIdentifier();
0652
0653   // get relevant bco matching information and initialize
0654   auto& bco_matching_information = m_bco_matching_information_map[packet_id];
0655   bco_matching_information.set_verbosity(Verbosity());
0656
0657   // decode
0658   const int data_length = packet->getDataLength();  // 32bit length
0659   const int data_padding = packet->getPadding();  // 32bit padding
0660
0661   // maximum number of dma words
0662   const size_t dma_words_buffer = static_cast<unsigned long>(data_length) * 2 / DAM_DMA_WORD_LENGTH + 1;
0663
0664   // dma words
0665   std::vector<dma_word> buffer(dma_words_buffer);
0666
0667   int l2 = 0;
0668   packet->fillIntArray(reinterpret_cast<int*>(buffer.data()), data_length + DAM_DMA_WORD_LENGTH / 2, &l2, "DATA");
0669
0670   // sanity checks
0671   assert(l2 <= data_length);
0672
0673   // sanity checks
0674   l2 -= data_padding;
0675   assert(l2 >= 0);
0676
0677   // actual number of dma words
0678   const size_t dma_words =  static_cast<unsigned long>(l2) * 2 / DAM_DMA_WORD_LENGTH;
0679   assert(dma_words <= buffer.size());
0680
0681   // demultiplexer
0682   for (size_t index = 0; index < dma_words; ++index)
0683   {
0684     const auto& dma_word_data = buffer[index];
0685
0686     if ((dma_word_data.dma_header & 0xFF00U) == GTM_MAGIC_KEY)
0687     {
0688
0689       // decode gtm data
0690       decode_gtm_data(packet_id, dma_word_data);
0691
0692     } else if ((dma_word_data.dma_header & 0xFF00U) == FEE_MAGIC_KEY) {
0693
0694       // decode fee data
0695       // get fee id
0696       const unsigned int fee_id = dma_word_data.dma_header & 0xffU;
0697
0698       // populate fee buffer
0699       if (fee_id < MAX_FEECOUNT)
0700       {
0701         // NOLINTNEXTLINE(modernize-loop-convert)
0702         for (unsigned int i = 0; i < DAM_DMA_WORD_LENGTH - 1; i++)
0703         { m_feeData[fee_id].push_back(dma_word_data.data[i]); }
0704
0705         // immediate fee buffer processing to reduce memory consuption
0706         process_fee_data(packet_id, fee_id);
0707       }
0708     }
0709   }
0710
0711 }
0712
0713 //____________________________________________________________________
0714 void SingleMicromegasPoolInput_v2::decode_gtm_data( int packet_id, const SingleMicromegasPoolInput_v2::dma_word& gtm_word )
0715 {
0716   const unsigned char* gtm = reinterpret_cast<const unsigned char*>(&gtm_word);
0717   MicromegasBcoMatchingInformation_v2::gtm_payload payload;
0718
0719   payload.pkt_type = gtm[0]|(unsigned short)(gtm[1] << 8U);
0720
0721   // check packet type
0722   if (payload.pkt_type != GTM_LVL1_ACCEPT_MAGIC_KEY &&
0723     payload.pkt_type != GTM_ENDAT_MAGIC_KEY &&
0724     payload.pkt_type != GTM_MODEBIT_MAGIC_KEY)
0725   { return; }
0726
0727   payload.is_lvl1 = payload.pkt_type == GTM_LVL1_ACCEPT_MAGIC_KEY;
0728   payload.is_endat = payload.pkt_type == GTM_ENDAT_MAGIC_KEY;
0729   payload.is_modebit = payload.pkt_type == GTM_MODEBIT_MAGIC_KEY;
0730
0731   payload.bco = ((unsigned long long) gtm[2] << 0U) | ((unsigned long long) gtm[3] << 8U) | ((unsigned long long) gtm[4] << 16U) | ((unsigned long long) gtm[5] << 24U) | ((unsigned long long) gtm[6] << 32U) | (((unsigned long long) gtm[7]) << 40U);
0732   payload.lvl1_count = ((unsigned int) gtm[8] << 0U) | ((unsigned int) gtm[9] << 8U) | ((unsigned int) gtm[10] << 16U) | ((unsigned int) gtm[11] << 24U);
0733   payload.endat_count = ((unsigned int) gtm[12] << 0U) | ((unsigned int) gtm[13] << 8U) | ((unsigned int) gtm[14] << 16U) | ((unsigned int) gtm[15] << 24U);
0734   payload.last_bco = ((unsigned long long) gtm[16] << 0U) | ((unsigned long long) gtm[17] << 8U) | ((unsigned long long) gtm[18] << 16U) | ((unsigned long long) gtm[19] << 24U) | ((unsigned long long) gtm[20] << 32U) | (((unsigned long long) gtm[21]) << 40U);
0735   payload.modebits = gtm[22];
0736   payload.userbits = gtm[23];
0737
0738   // fill per packet heartbeat statistics
0739   const bool is_heartbeat=payload.modebits&(1U<<ELINK_HEARTBEAT_T);
0740   if( is_heartbeat )
0741   { h_heartbeat_stat->Fill(std::to_string(packet_id).c_str(),1); }
0742
0743   // save bco information
0744   auto& bco_matching_information = m_bco_matching_information_map[packet_id];
0745   bco_matching_information.save_gtm_bco_information(packet_id, payload);
0746
0747   if(payload.is_lvl1||payload.is_endat)
0748   {
0749     const auto& gtm_bco = payload.bco;
0750     m_BeamClockPacket[gtm_bco].insert(packet_id);
0751     m_BclkStack.insert(gtm_bco);
0752   }
0753
0754   // find reference from modebits, using BX_COUNTER_SYNC_T
0755   /*
0756   * This needs to be done even if the bco matching information is already verified
0757   * because any BX_COUNTER_SYNC_T event will break past references
0758   */
0759   bco_matching_information.find_reference_from_modebits(payload);
0760
0761   // store in running waveform
0762   if( m_do_evaluation )
0763   {
0764     m_waveform.packet_id = packet_id;
0765     m_waveform.gtm_bco_first = bco_matching_information.get_gtm_bco_first();
0766     m_waveform.gtm_bco = bco_matching_information.get_gtm_bco_last();
0767
0768     {
0769       const auto predicted = bco_matching_information.get_predicted_fee_bco(m_waveform.gtm_bco);;
0770       if( predicted ) m_waveform.fee_bco_predicted = predicted.value();
0771     }
0772
0773     m_waveform.fee_bco_first = bco_matching_information.get_fee_bco_first();
0774   }
0775 }
0776
0777 //____________________________________________________________________
0778 void SingleMicromegasPoolInput_v2::process_fee_data( int packet_id, unsigned int fee_id )
0779 {
0780   // get bco information
0781   auto& bco_matching_information = m_bco_matching_information_map[packet_id];
0782
0783   // get fee buffer
0784   assert(fee_id < m_feeData.size());
0785   auto& data_buffer = m_feeData[fee_id];
0786
0787   // process header
0788   while (HEADER_LENGTH <= data_buffer.size())
0789   {
0790
0791     // packet loop
0792     /* question: why check index [1] ? */
0793     if (data_buffer[1] != FEE_PACKET_MAGIC_KEY_1)
0794     {
0795       if (Verbosity())
0796       {
0797         std::cout << "SingleMicromegasPoolInput_v2::process_fee_data - Invalid FEE magic key at position 1 0x" << std::hex << data_buffer[1] << std::dec << std::endl;
0798       }
0799       data_buffer.pop_front();
0800       continue;
0801     }
0802
0803     if (data_buffer[2] != FEE_PACKET_MAGIC_KEY_2)
0804     {
0805       if (Verbosity())
0806       {
0807         std::cout << "SingleMicromegasPoolInput_v2::process_fee_data - Invalid FEE magic key at position 2 0x" << std::hex << data_buffer[2] << std::dec << std::endl;
0808       }
0809       data_buffer.pop_front();
0810       continue;
0811     }
0812
0813     // packet length
0814     // number of 10-bit words + 5 in this packet
0815     const uint16_t& pkt_length = data_buffer[0];
0816     if (pkt_length > MAX_PACKET_LENGTH)
0817     {
0818       if (Verbosity())
0819       {
0820         std::cout << "SingleMicromegasPoolInput_v2::process_fee_data - Invalid FEE pkt_length " << pkt_length << std::endl;
0821       }
0822       data_buffer.pop_front();
0823       continue;
0824     }
0825
0826     // compare to data size
0827     if (pkt_length + 1U > data_buffer.size())
0828     {
0829       // not enough data. will wait for more
0830       break;
0831     }
0832
0833     // create payload
0834     MicromegasBcoMatchingInformation_v2::fee_payload payload;
0835
0836     // number of 10-bit words in this packet
0837     payload.adc_length = data_buffer[0] - HEADER_LENGTH;
0838     payload.data_parity = data_buffer[4] >> 9U;
0839     payload.sampa_address = (uint16_t)(data_buffer[4] >> 5U) & 0xfU;
0840     payload.sampa_channel = data_buffer[4] & 0x1fU;
0841     payload.channel = data_buffer[4] & 0x1ffU;
0842     payload.type = (uint16_t)(data_buffer[3] >> 7U) & 0x7U;
0843     payload.user_word = data_buffer[3] & 0x7fU;
0844     payload.bx_timestamp = (uint32_t)((data_buffer[6] & 0x3ffU) << 10U)|(data_buffer[5] & 0x3ffU);
0845
0846     // fill heartbeat statistics
0847     const bool is_heartbeat(payload.type==HEARTBEAT_T);
0848     if( is_heartbeat )
0849     {
0850       const int sampa_id = payload.sampa_address + fee_id*m_nsampa_fee;
0851       h_heartbeat_stat->Fill(sampa_id + 2);
0852     }
0853
0854     // increment number of waveforms
0855     ++m_waveform_counters[packet_id].total;
0856     ++m_fee_waveform_counters[fee_id].total;
0857     h_waveform_count_total->Fill(std::to_string(packet_id).c_str(),1);
0858     h_fee_waveform_count_total->Fill(fee_id,1);
0859
0860     if( is_heartbeat )
0861     {
0862       ++m_heartbeat_counters[packet_id].total;
0863       ++m_fee_heartbeat_counters[fee_id].total;
0864     }
0865
0866     // crc
0867     payload.data_crc = data_buffer[pkt_length];
0868     // payload.calc_crc = crc16(data_buffer, 0, pkt_length);
0869
0870     // make sure buffer is cleaned as soon as we exit current scope
0871     /*
0872      * the buffer is cleared when buffer_cleaner is deleted, that is, as soon as it becomes out of scope
0873      * this allows to use the various 'continue' statements below, without having to care about the buffer being properly cleared
0874      */
0875     [[maybe_unused]] buffer_cleaner_t buffer_cleaner( data_buffer, pkt_length );
0876
0877     // check bco matching information
0878     if (!bco_matching_information.is_verified())
0879     {
0880       bco_matching_information.find_reference_from_data(payload);
0881     }
0882
0883     // if bco matching information is still not verified, drop the packet
0884     if (!bco_matching_information.is_verified())
0885     {
0886       ++m_waveform_counters[packet_id].dropped_bco;
0887       ++m_fee_waveform_counters[fee_id].dropped_bco;
0888       h_waveform_count_dropped_bco->Fill(std::to_string(packet_id).c_str(), 1);
0889       h_fee_waveform_count_dropped_bco->Fill(fee_id, 1);
0890
0891       if( is_heartbeat )
0892       {
0893         ++m_heartbeat_counters[packet_id].dropped_bco;
0894         ++m_fee_heartbeat_counters[fee_id].dropped_bco;
0895       }
0896
0897       continue;
0898     }
0899
0900     // try get gtm bco matching fee
0901     const auto& fee_bco = payload.bx_timestamp;
0902
0903     // find matching gtm bco
0904     uint64_t gtm_bco = 0;
0905     const auto result = bco_matching_information.find_gtm_bco(packet_id, fee_id, fee_bco);
0906     if (result)
0907     {
0908       // assign gtm bco
0909       gtm_bco = result.value();
0910     } else {
0911
0912       // increment counter and histogram
0913       ++m_waveform_counters[packet_id].dropped_bco;
0914       ++m_fee_waveform_counters[fee_id].dropped_bco;
0915       h_waveform_count_dropped_bco->Fill( std::to_string(packet_id).c_str(), 1 );
0916       h_fee_waveform_count_dropped_bco->Fill(fee_id, 1);
0917
0918       if( is_heartbeat )
0919       {
0920         ++m_heartbeat_counters[packet_id].dropped_bco;
0921         ++m_fee_heartbeat_counters[fee_id].dropped_bco;
0922       }
0923
0924       // skip the waverform
0925       continue;
0926     }
0927
0928     if( m_do_evaluation )
0929     {
0930       m_waveform.is_heartbeat = (payload.type == HEARTBEAT_T);
0931       m_waveform.fee_id = fee_id;
0932       m_waveform.channel = payload.channel;
0933       m_waveform.fee_bco = fee_bco;
0934
0935       m_waveform.gtm_bco_matched = gtm_bco;
0936       {
0937         const auto predicted = bco_matching_information.get_predicted_fee_bco(gtm_bco);;
0938         if( predicted ) m_waveform.fee_bco_predicted_matched = predicted.value();
0939       }
0940
0941       m_evaluation_tree->Fill();
0942     }
0943
0944     // ignore heartbeat waveforms
0945     if( is_heartbeat ) { continue; }
0946
0947     // store data from string
0948     // Format is (N sample) (start time), (1st sample)... (Nth sample)
0949     size_t pos = HEADER_LENGTH;
0950     while (pos < pkt_length )
0951     {
0952       const uint16_t& samples = data_buffer[pos++];
0953       const uint16_t& start_t = data_buffer[pos++];
0954       if (pos + samples > size_t(pkt_length) )
0955       {
0956         if (Verbosity())
0957         {
0958           std::cout << "SingleMicromegasPoolInput_v2::process_fee_data -"
0959             << " samples: " << samples
0960             << " pos: " << pos
0961             << " pkt_length: " << pkt_length
0962             << " format error"
0963             << std::endl;
0964         }
0965         break;
0966       }
0967
0968       std::vector<uint16_t> adc(samples);
0969       for (int i = 0; i < samples; ++i)
0970       { adc[i] = data_buffer[pos++]; }
0971
0972       // add
0973       payload.waveforms.emplace_back(start_t,std::move(adc));
0974     }
0975
0976     // create new hit
0977     auto newhit = std::make_unique<MicromegasRawHitv3>();
0978     newhit->set_bco(fee_bco);
0979     newhit->set_gtm_bco(gtm_bco);
0980
0981     // packet id, fee id, channel, etc.
0982     newhit->set_packetid(packet_id);
0983     newhit->set_fee(fee_id);
0984     newhit->set_channel(payload.channel);
0985     newhit->set_sampaaddress(payload.sampa_address);
0986     newhit->set_sampachannel(payload.sampa_channel);
0987
0988     // adc values
0989     for( auto&& [start_t, adc]:payload.waveforms )
0990     { newhit->move_adc_waveform(start_t,std::move(adc)); }
0991
0992     m_BeamClockFEE[gtm_bco].insert(fee_id);
0993     m_FEEBclkMap[fee_id] = gtm_bco;
0994
0995     if (StreamingInputManager())
0996     {
0997       StreamingInputManager()->AddMicromegasRawHit(gtm_bco, newhit.get());
0998     }
0999
1000     m_MicromegasRawHitMap[gtm_bco].push_back(newhit.release());
1001   }
1002 }