framework/fun4allraw/SingleMicromegasPoolInput_v1.cc

0001 #include "SingleMicromegasPoolInput_v1.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
0023 #include <TFile.h>
0024 #include <TH1.h>
0025
0026 #include <algorithm>
0027 #include <bitset>
0028
0029 namespace
0030 {
0031
0032   //! mark invalid ADC values
0033   static constexpr uint16_t m_adc_invalid = 65000;
0034
0035   // maximum number of packets
0036   static constexpr int m_npackets_active = 2;
0037
0038   // minimum number of requested samples
0039   static constexpr int m_min_req_samples = 5;
0040
0041   /* see: https://git.racf.bnl.gov/gitea/Instrumentation/sampa_data/src/branch/fmtv2/README.md */
0042   enum SampaDataType
0043   {
0044     HEARTBEAT_T = 0b000,
0045     TRUNCATED_DATA_T = 0b001,
0046     TRUNCATED_TRIG_EARLY_DATA_T = 0b011,
0047     NORMAL_DATA_T = 0b100,
0048     LARGE_DATA_T = 0b101,
0049     TRIG_EARLY_DATA_T = 0b110,
0050     TRIG_EARLY_LARGE_DATA_T = 0b111,
0051   };
0052
0053 }  // namespace
0054
0055 //______________________________________________________________
0056 SingleMicromegasPoolInput_v1::SingleMicromegasPoolInput_v1(const std::string& name)
0057   : SingleStreamingInput(name)
0058 {
0059   SubsystemEnum(InputManagerType::MICROMEGAS);
0060   m_rawHitContainerName = "MICROMEGASRAWHIT";
0061 }
0062
0063 //______________________________________________________________
0064 SingleMicromegasPoolInput_v1::~SingleMicromegasPoolInput_v1()
0065 {
0066
0067   std::cout << "SingleMicromegasPoolInput_v1::~SingleMicromegasPoolInput_v1 - runnumber: " << RunNumber() << std::endl;
0068
0069   // timer statistics
0070   m_timer.print_stat();
0071
0072   for( const auto& [packet,counts]:m_waveform_count_total )
0073   {
0074     const auto dropped_bco =  m_waveform_count_dropped_bco[packet];
0075     const auto dropped_pool =  m_waveform_count_dropped_pool[packet];
0076     std::cout << "SingleMicromegasPoolInput_v1::~SingleMicromegasPoolInput_v1 -"
0077       << " packet: " << packet
0078       << " wf_total: " << counts
0079       << " wf_dropped_bco: " << dropped_bco
0080       << " wf_dropped_pool: " << dropped_pool
0081       << " ratio_bco: " << double(dropped_bco)/counts
0082       << " ratio_pool: " << double(dropped_pool)/counts
0083       << std::endl;
0084   }
0085   std::cout << std::endl;
0086
0087   // drop per fee statistics
0088   for( const auto& [fee,counts]:m_fee_waveform_count_total )
0089   {
0090     const auto dropped_bco =  m_fee_waveform_count_dropped_bco[fee];
0091     std::cout << "SingleMicromegasPoolInput_v1::~SingleMicromegasPoolInput_v1 -"
0092       << " fee: " << fee
0093       << " wf_total: " << counts
0094       << " wf_dropped_bco: " << dropped_bco
0095       << " ratio_bco: " << double(dropped_bco)/counts
0096       << std::endl;
0097
0098   }
0099   std::cout << std::endl;
0100
0101   // also printout adjusted multipliers
0102   for( const auto& [packet_id, bco_matching]:m_bco_matching_information_map )
0103   {
0104     const auto old_precision = std::cout.precision();
0105     std::cout << "SingleMicromegasPoolInput_v1::~SingleMicromegasPoolInput_v1 -"
0106       << " packet: " << packet_id
0107       << " multiplier: " <<  std::setprecision(9) << bco_matching.get_adjusted_multiplier() << std::setprecision(old_precision)
0108       << std::endl;
0109   }
0110
0111 }
0112
0113 //______________________________________________________________
0114 void SingleMicromegasPoolInput_v1::FillPool(const unsigned int /*nbclks*/)
0115 {
0116   if (AllDone())  // no more files and all events read
0117   {
0118     return;
0119   }
0120
0121   while (!GetEventiterator())  // at startup this is a null pointer
0122   {
0123     if (!OpenNextFile())
0124     {
0125       AllDone(1);
0126       return;
0127     }
0128   }
0129
0130   while (GetSomeMoreEvents())
0131   {
0132     std::unique_ptr<Event> evt(GetEventiterator()->getNextEvent());
0133     while (!evt)
0134     {
0135       fileclose();
0136       if (!OpenNextFile())
0137       {
0138         AllDone(1);
0139         return;
0140       }
0141
0142       // get next event
0143       evt.reset(GetEventiterator()->getNextEvent());
0144     }
0145
0146     if (Verbosity() > 2)
0147     {
0148       std::cout << "Fetching next Event" << evt->getEvtSequence() << std::endl;
0149     }
0150
0151     if (evt->getEvtType() != DATAEVENT)
0152     {
0153       m_NumSpecialEvents++;
0154       continue;
0155     }
0156
0157     RunNumber(evt->getRunNumber());
0158     if (Verbosity() > 1)
0159     {
0160       evt->identify();
0161     }
0162
0163     const int EventSequence = evt->getEvtSequence();
0164     const int npackets = evt->getPacketList(&plist[0], 10);
0165
0166     if (npackets == 10)
0167     {
0168       std::cout << "SingleMicromegasPoolInput_v1::FillPool - too many packets" << std::endl;
0169       exit(1);
0170     }
0171
0172     for (int i = 0; i < npackets; i++)
0173     {
0174       // keep pointer to local packet
0175       std::unique_ptr<Packet> packet(plist[i]);
0176
0177       // get packet id
0178       const auto packet_id = packet->getIdentifier();
0179
0180       if (Verbosity() > 1)
0181       {
0182         packet->identify();
0183       }
0184
0185       // get relevant bco matching information
0186       auto& bco_matching_information = m_bco_matching_information_map[packet_id];
0187       bco_matching_information.set_verbosity(Verbosity());
0188
0189       // read gtm bco information
0190       bco_matching_information.save_gtm_bco_information(packet.get());
0191
0192       // save BCO from tagger internally
0193       const int n_tagger = packet->lValue(0, "N_TAGGER");
0194       for (int t = 0; t < n_tagger; ++t)
0195       {
0196         const bool is_lvl1 = static_cast<uint8_t>(packet->lValue(t, "IS_LEVEL1_TRIGGER"));
0197         const bool is_endat = static_cast<uint8_t>(packet->lValue(t, "IS_ENDAT"));
0198         if (is_lvl1||is_endat)
0199         {
0200           const uint64_t gtm_bco = static_cast<uint64_t>(packet->lValue(t, "BCO"));
0201           m_BeamClockPacket[gtm_bco].insert(packet_id);
0202           m_BclkStack.insert(gtm_bco);
0203         }
0204       }
0205
0206       // loop over waveforms
0207       const int nwf = packet->iValue(0, "NR_WF");
0208
0209       // increment counter and histogram
0210       m_waveform_count_total[packet_id] += nwf;
0211       h_waveform_count_total->Fill( std::to_string(packet_id).c_str(), nwf );
0212
0213       if (Verbosity())
0214       {
0215         std::cout
0216           << "SingleMicromegasPoolInput_v1::FillPool -"
0217           << " packet: " << packet_id
0218           << " n_waveform: " << nwf
0219           << std::endl;
0220         bco_matching_information.print_gtm_bco_information();
0221       }
0222
0223       // find reference from modebits, using BX_COUNTER_SYNC_T
0224       /*
0225        * This needs to be done even if the bco matching information is already verified
0226        * because any BX_COUNTER_SYNC_T event will break past references
0227        */
0228       bco_matching_information.find_reference_from_modebits(packet.get());
0229
0230       // if bco matching information is not verified, try find reference from data
0231       if (!bco_matching_information.is_verified())
0232       {
0233         bco_matching_information.find_reference_from_data(packet.get());
0234       }
0235
0236       // if bco matching information is still not verified, drop the packet
0237       if (!bco_matching_information.is_verified())
0238       {
0239         std::cout << "SingleMicromegasPoolInput_v1::FillPool - bco_matching not verified, dropping packet" << std::endl;
0240         m_waveform_count_dropped_bco[packet_id] += nwf;
0241         h_waveform_count_dropped_bco->Fill( std::to_string(packet_id).c_str(), nwf );
0242         continue;
0243       }
0244
0245       for (int wf = 0; wf < nwf; wf++)
0246       {
0247         // get fee id
0248         const int fee_id = packet->iValue(wf, "FEE");
0249         ++m_fee_waveform_count_total[fee_id];
0250
0251         // get checksum_error and check
0252         const auto checksum_error = packet->iValue(wf, "CHECKSUMERROR");
0253         if (checksum_error)
0254         {
0255           continue;
0256         }
0257
0258         // get fee bco
0259         const unsigned int fee_bco = packet->iValue(wf, "BCO");
0260
0261         // find matching gtm bco
0262         uint64_t gtm_bco = 0;
0263         const auto result = bco_matching_information.find_gtm_bco(fee_bco);
0264
0265         if (result)
0266         {
0267           // assign gtm bco
0268           gtm_bco = result.value();
0269         }
0270         else
0271         {
0272           // increment counter and histogram
0273           ++m_waveform_count_dropped_bco[packet_id];
0274           ++m_fee_waveform_count_dropped_bco[fee_id];
0275           h_waveform_count_dropped_bco->Fill( std::to_string(packet_id).c_str(), 1 );
0276
0277           // skip the waverform
0278           continue;
0279         }
0280
0281         // get type
0282         // ignore heartbeat waveforms
0283         if( packet->iValue(wf, "TYPE" ) == HEARTBEAT_T ) continue;
0284
0285         // get number of samples and check
0286         const uint16_t samples = packet->iValue(wf, "SAMPLES");
0287         if (samples < m_min_req_samples)
0288         {
0289           continue;
0290         }
0291
0292         // create new hit
0293         auto newhit = std::make_unique<MicromegasRawHitv3>();
0294         newhit->set_bco(fee_bco);
0295         newhit->set_gtm_bco(gtm_bco);
0296
0297         // packet id, fee id, channel, etc.
0298         newhit->set_packetid(packet_id);
0299         newhit->set_fee(fee_id);
0300         newhit->set_channel(packet->iValue(wf, "CHANNEL"));
0301         newhit->set_sampaaddress(packet->iValue(wf, "SAMPAADDRESS"));
0302         newhit->set_sampachannel(packet->iValue(wf, "CHANNEL"));
0303
0304         // adc values
0305         for (uint16_t is = 0; is < samples; ++is)
0306         {
0307           uint16_t adc = packet->iValue(wf, is);
0308           if( adc == m_adc_invalid)
0309           {
0310             continue;
0311           }
0312
0313           uint16_t first = is;
0314           MicromegasRawHitv3::adc_list_t values;
0315           for( ;is<samples && (adc = packet->iValue(wf, is)) != m_adc_invalid; ++is )
0316           { values.push_back(adc); }
0317           newhit->move_adc_waveform( first, std::move(values));
0318
0319         }
0320
0321         m_BeamClockFEE[gtm_bco].insert(fee_id);
0322         m_FEEBclkMap[fee_id] = gtm_bco;
0323         if (Verbosity() > 2)
0324         {
0325           std::cout << "evtno: " << EventSequence
0326                     << ", hits: " << wf
0327                     << ", num waveforms: " << nwf
0328                     << ", bco: 0x" << std::hex << gtm_bco << std::dec
0329                     << ", FEE: " << fee_id << std::endl;
0330         }
0331
0332         if (StreamingInputManager())
0333         {
0334           StreamingInputManager()->AddMicromegasRawHit(gtm_bco, newhit.get());
0335         }
0336
0337         m_MicromegasRawHitMap[gtm_bco].push_back(newhit.release());
0338       }
0339     }
0340   }
0341 }
0342
0343 //______________________________________________________________
0344 void SingleMicromegasPoolInput_v1::Print(const std::string& what) const
0345 {
0346   if (what == "ALL" || what == "FEE")
0347   {
0348     for (const auto& bcliter : m_BeamClockFEE)
0349     {
0350       std::cout << "Beam clock 0x" << std::hex << bcliter.first << std::dec << std::endl;
0351       for (auto feeiter : bcliter.second)
0352       {
0353         std::cout << "FEM: " << feeiter << std::endl;
0354       }
0355     }
0356   }
0357
0358   if (what == "ALL" || what == "FEEBCLK")
0359   {
0360     for (auto bcliter : m_FEEBclkMap)
0361     {
0362       std::cout << "FEE" << bcliter.first << " bclk: 0x"
0363                 << std::hex << bcliter.second << std::dec << std::endl;
0364     }
0365   }
0366
0367   if (what == "ALL" || what == "STORAGE")
0368   {
0369     for (const auto& bcliter : m_MicromegasRawHitMap)
0370     {
0371       std::cout << "Beam clock 0x" << std::hex << bcliter.first << std::dec << std::endl;
0372       for (auto feeiter : bcliter.second)
0373       {
0374         std::cout
0375             << "fee: " << feeiter->get_fee()
0376             << " at " << std::hex << feeiter << std::dec
0377             << std::endl;
0378       }
0379     }
0380   }
0381
0382   if (what == "ALL" || what == "STACK")
0383   {
0384     for (const auto& iter : m_BclkStack)
0385     {
0386       std::cout << "stacked bclk: 0x" << std::hex << iter << std::dec << std::endl;
0387     }
0388   }
0389 }
0390
0391 //____________________________________________________________________________
0392 void SingleMicromegasPoolInput_v1::CleanupUsedPackets(const uint64_t bclk, bool dropped)
0393 {
0394
0395   // delete all raw hits associated to bco smaller than reference, and remove from map
0396   for(auto iter = m_MicromegasRawHitMap.begin(); iter != m_MicromegasRawHitMap.end() && (iter->first <= bclk); iter = m_MicromegasRawHitMap.erase(iter))
0397   {
0398     for (const auto& rawhit : iter->second)
0399     {
0400       if( dropped )
0401       {
0402         // increment dropped waveform counter and histogram
0403         ++m_waveform_count_dropped_pool[rawhit->get_packetid()];
0404         h_waveform_count_dropped_pool->Fill( std::to_string(rawhit->get_packetid()).c_str(), 1 );
0405       }
0406       delete rawhit;
0407     }
0408   }
0409
0410   // cleanup bco stacks
0411   /* it erases all elements for which the bco is no greater than the provided one */
0412   m_BclkStack.erase(m_BclkStack.begin(), m_BclkStack.upper_bound(bclk));
0413   m_BeamClockFEE.erase(m_BeamClockFEE.begin(), m_BeamClockFEE.upper_bound(bclk));
0414   m_BeamClockPacket.erase(m_BeamClockPacket.begin(), m_BeamClockPacket.upper_bound(bclk));
0415
0416   // cleanup matching information
0417   for( auto&& bco_matching:m_bco_matching_information_map )
0418   { bco_matching.second.cleanup(bclk); }
0419
0420 }
0421
0422 //_______________________________________________________
0423 void SingleMicromegasPoolInput_v1::ClearCurrentEvent()
0424 {
0425   std::cout << "SingleMicromegasPoolInput_v1::ClearCurrentEvent." << std::endl;
0426   uint64_t currentbclk = *m_BclkStack.begin();
0427   CleanupUsedPackets(currentbclk);
0428   return;
0429 }
0430
0431 //_______________________________________________________
0432 bool SingleMicromegasPoolInput_v1::GetSomeMoreEvents()
0433 {
0434   if (AllDone())
0435   {
0436     return false;
0437   }
0438
0439   // check minimum pool size
0440   if (m_MicromegasRawHitMap.size() < m_BcoPoolSize)
0441   {
0442     return true;
0443   }
0444
0445   // make sure that the latest BCO received by each FEEs is past the current BCO
0446   std::set<int> toerase;
0447   uint64_t lowest_bclk = m_MicromegasRawHitMap.begin()->first + m_BcoRange;
0448   for (auto bcliter : m_FEEBclkMap)
0449   {
0450     if (bcliter.second <= lowest_bclk)
0451     {
0452       uint64_t highest_bclk = m_MicromegasRawHitMap.rbegin()->first;
0453       if ((highest_bclk - m_MicromegasRawHitMap.begin()->first) < MaxBclkDiff())
0454       {
0455         return true;
0456       }
0457       else
0458       {
0459         std::cout << PHWHERE << Name() << ": erasing FEE " << bcliter.first
0460                   << " with stuck bclk: " << std::hex << bcliter.second
0461                   << " current bco range: 0x" << m_MicromegasRawHitMap.begin()->first
0462                   << ", to: 0x" << highest_bclk << ", delta: " << std::dec
0463                   << (highest_bclk - m_MicromegasRawHitMap.begin()->first)
0464                   << std::dec << std::endl;
0465         toerase.insert(bcliter.first);
0466       }
0467     }
0468   }
0469   for(auto& fee: toerase)
0470   {
0471     m_FEEBclkMap.erase(fee);
0472   }
0473
0474   return false;
0475 }
0476
0477 //_______________________________________________________
0478 void SingleMicromegasPoolInput_v1::CreateDSTNode(PHCompositeNode* topNode)
0479 {
0480   PHNodeIterator iter(topNode);
0481   auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0482   if (!dstNode)
0483   {
0484     dstNode = new PHCompositeNode("DST");
0485     topNode->addNode(dstNode);
0486   }
0487
0488   PHNodeIterator iterDst(dstNode);
0489   auto detNode = dynamic_cast<PHCompositeNode*>(iterDst.findFirst("PHCompositeNode", "MICROMEGAS"));
0490   if (!detNode)
0491   {
0492     detNode = new PHCompositeNode("MICROMEGAS");
0493     dstNode->addNode(detNode);
0494   }
0495
0496   auto container = findNode::getClass<MicromegasRawHitContainer>(detNode, m_rawHitContainerName);
0497   if (!container)
0498   {
0499     container = new MicromegasRawHitContainerv3();
0500     auto newNode = new PHIODataNode<PHObject>(container, m_rawHitContainerName, "PHObject");
0501     detNode->addNode(newNode);
0502   }
0503 }
0504
0505 //_______________________________________________________
0506 void SingleMicromegasPoolInput_v1::ConfigureStreamingInputManager()
0507 {
0508   if (StreamingInputManager())
0509   {
0510     StreamingInputManager()->SetMicromegasBcoRange(m_BcoRange);
0511     StreamingInputManager()->SetMicromegasNegativeBco(m_NegativeBco);
0512   }
0513 }
0514
0515 //_______________________________________________________
0516 void SingleMicromegasPoolInput_v1::FillBcoQA(uint64_t gtm_bco)
0517 {
0518   auto hm = QAHistManagerDef::getHistoManager();
0519   assert(hm);
0520
0521   // set of packets for which the gtm BCO is found at least once
0522   unsigned int n_waveforms = 0;
0523   std::set<int> found_packets;
0524   for (uint64_t gtm_bco_loc = gtm_bco - m_NegativeBco; gtm_bco_loc < gtm_bco + m_BcoRange - m_NegativeBco; ++gtm_bco_loc)
0525   {
0526     // packet ids
0527     const auto packet_iter = m_BeamClockPacket.find(gtm_bco_loc);
0528     if( packet_iter != m_BeamClockPacket.end() )
0529     { found_packets.insert( packet_iter->second.begin(), packet_iter->second.end() ); }
0530
0531     // waveforms
0532     const auto wf_iter = m_MicromegasRawHitMap.find(gtm_bco_loc);
0533     if (wf_iter != m_MicromegasRawHitMap.end())
0534     { n_waveforms += wf_iter->second.size(); }
0535   }
0536
0537   // per packet statistics
0538   h_packet_stat->Fill( "Reference", 1 );
0539
0540   for( const auto& packet_id:found_packets )
0541   { h_packet_stat->Fill( std::to_string(packet_id).c_str(), 1 ); }
0542   h_packet_stat->Fill( "All", found_packets.size()>= m_npackets_active );
0543
0544   // how many packet_id found for this BCO
0545   h_packet->Fill(found_packets.size());
0546
0547   // how many waveforms found for this BCO
0548   h_waveform->Fill(n_waveforms);
0549 }
0550 //_______________________________________________________
0551 void SingleMicromegasPoolInput_v1::createQAHistos()
0552 {
0553   auto hm = QAHistManagerDef::getHistoManager();
0554   assert(hm);
0555
0556   // number of packets found with BCO from felix matching reference BCO
0557   h_packet = new TH1I( "h_MicromegasBCOQA_npacket_bco", "TPOT Packet Count per GTM BCO; Matching BCO tagger count; GL1 trigger count", 10, 0, 10 );
0558
0559   // number of waveforms found with BCO from felix matching reference BCO
0560   h_waveform = new TH1I( "h_MicromegasBCOQA_nwaveform_bco", "TPOT Waveform Count per GTM BCO; Matching Waveform count; GL1 trigger count", 4100, 0, 4100 );
0561
0562   /*
0563    * first bin is the number of requested GL1 BCO, for reference
0564    * next two bins is the number of times the GL1 BCO is found in the taggers list for a given packet_id
0565    * last bin is the sum
0566    */
0567   h_packet_stat = new TH1I( "h_MicromegasBCOQA_packet_stat", "Matching Tagger count per packet; packet id; GL1 trigger count", m_npackets_active+2, 0, m_npackets_active+2 );
0568   h_packet_stat->GetXaxis()->SetBinLabel(1, "Reference" );
0569   h_packet_stat->GetXaxis()->SetBinLabel(2, "5001" );
0570   h_packet_stat->GetXaxis()->SetBinLabel(3, "5002" );
0571   h_packet_stat->GetXaxis()->SetBinLabel(4, "All" );
0572   h_packet_stat->GetYaxis()->SetTitle( "trigger count" );
0573
0574   // total number of waveform per packet
0575   h_waveform_count_total = new TH1I( "h_MicromegasBCOQA_waveform_count_total", "Total number of waveforms per packet", m_npackets_active, 0, m_npackets_active );
0576
0577   // number of dropped waveform per packet due to bco mismatch
0578   h_waveform_count_dropped_bco = new TH1I( "h_MicromegasBCOQA_waveform_count_dropped_bco", "Number of dropped waveforms per packet (bco)", m_npackets_active, 0, m_npackets_active );
0579
0580   // number of dropped waveform per packet due to fun4all pool mismatch
0581   h_waveform_count_dropped_pool = new TH1I( "h_MicromegasBCOQA_waveform_count_dropped_pool", "Number of dropped waveforms per packet (pool)", m_npackets_active, 0, m_npackets_active );
0582
0583   // define axis
0584   for( const auto& h:std::initializer_list<TH1*>{h_waveform_count_total, h_waveform_count_dropped_bco, h_waveform_count_dropped_pool} )
0585   {
0586     h->GetXaxis()->SetBinLabel(1, "5001" );
0587     h->GetXaxis()->SetBinLabel(2, "5002" );
0588     h->GetYaxis()->SetTitle( "waveform count" );
0589     h->GetXaxis()->SetTitle( "packet id" );
0590   }
0591
0592   // register all histograms to histogram manager
0593   for( const auto& h:std::initializer_list<TH1*>{h_packet, h_waveform, h_packet_stat, h_waveform_count_total, h_waveform_count_dropped_bco, h_waveform_count_dropped_pool} )
0594   {
0595     h->SetFillStyle(1001);
0596     h->SetFillColor(kYellow);
0597     hm->registerHisto(h);
0598   }
0599 }