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