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