framework/fun4allraw/SingleTpcPoolInput.cc

0001 #include "SingleTpcPoolInput.h"
0002
0003 #include "Fun4AllStreamingInputManager.h"
0004 #include "InputManagerType.h"
0005
0006 #include <ffarawobjects/TpcRawHitContainerv2.h>
0007 #include <ffarawobjects/TpcRawHitv2.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0011 #include <phool/getClass.h>
0012 #include <phool/phool.h>
0013
0014 #include <Event/Event.h>
0015 #include <Event/EventTypes.h>
0016 #include <Event/Eventiterator.h>
0017
0018 #include <memory>
0019 #include <set>
0020
0021 SingleTpcPoolInput::SingleTpcPoolInput(const std::string &name)
0022   : SingleStreamingInput(name)
0023 {
0024   SubsystemEnum(InputManagerType::TPC);
0025   m_rawHitContainerName = "TPCRAWHIT";
0026 }
0027
0028 void SingleTpcPoolInput::FillPool(const uint64_t minBCO)
0029 {
0030   if (AllDone())  // no more files and all events read
0031   {
0032     return;
0033   }
0034   while (GetEventiterator() == nullptr)  // at startup this is a null pointer
0035   {
0036     if (!OpenNextFile())
0037     {
0038       AllDone(1);
0039       return;
0040     }
0041   }
0042   //  std::set<uint64_t> saved_beamclocks;
0043   while (GetSomeMoreEvents(0))
0044   {
0045     std::unique_ptr<Event> evt(GetEventiterator()->getNextEvent());
0046     while (!evt)
0047     {
0048       fileclose();
0049       if (!OpenNextFile())
0050       {
0051         AllDone(1);
0052         return;
0053       }
0054       evt.reset(GetEventiterator()->getNextEvent());
0055     }
0056     if (Verbosity() > 2)
0057     {
0058       std::cout << "Fetching next Event" << evt->getEvtSequence() << std::endl;
0059     }
0060     RunNumber(evt->getRunNumber());
0061     if (GetVerbosity() > 1)
0062     {
0063       evt->identify();
0064     }
0065     if (evt->getEvtType() != DATAEVENT)
0066     {
0067       m_NumSpecialEvents++;
0068       if (evt->getEvtType() == ENDRUNEVENT)
0069       {
0070         std::cout << "End run flag for " << Name() << " found, remaining TPC data is corrupted" << std::endl;
0071         AllDone(1);
0072         return;
0073       }
0074       continue;
0075     }
0076     int EventSequence = evt->getEvtSequence();
0077     std::vector<Packet *> pktvec = evt->getPacketVector();
0078     if (m_skipEarlyEvents)
0079     {
0080       for (auto packet : pktvec)
0081       {
0082         int numBCOs = packet->lValue(0, "N_TAGGER");
0083         bool bco_set = false;
0084         for (int j = 0; j < numBCOs; j++)
0085         {
0086           const auto is_lvl1 = static_cast<uint8_t>(packet->lValue(j, "IS_LEVEL1_TRIGGER"));
0087           const auto is_endat = static_cast<uint8_t>(packet->lValue(j, "IS_ENDAT"));
0088           if ((is_lvl1 || is_endat) && !bco_set)
0089           {
0090             bco_set = true;
0091             uint64_t bco = packet->lValue(j, "BCO");
0092             if (bco < minBCO)
0093             {
0094               continue;
0095             }
0096             m_skipEarlyEvents = false;
0097           }
0098         }
0099       }
0100     }
0101     if (m_skipEarlyEvents)
0102     {
0103       for (auto packet : pktvec)
0104       {
0105         delete packet;
0106       }
0107       continue;
0108     }
0109     for (auto packet : pktvec)
0110     {
0111       // get packet id
0112       const auto packet_id = packet->getIdentifier();
0113
0114       if (Verbosity() > 1)
0115       {
0116         packet->identify();
0117       }
0118
0119       // by default use previous bco clock for gtm bco
0120       auto &previous_bco = m_packet_bco[packet_id];
0121       uint64_t gtm_bco = previous_bco;
0122
0123       uint64_t m_nTaggerInFrame = packet->lValue(0, "N_TAGGER");
0124       bool skipthis = true;
0125       uint64_t largest_bco = 0;
0126       bool bco_set = false;
0127       for (uint64_t t = 0; t < m_nTaggerInFrame; t++)
0128       {
0129         // only store gtm_bco for level1 type of taggers (not ENDDAT)
0130         const auto is_lvl1 = static_cast<uint8_t>(packet->lValue(t, "IS_LEVEL1_TRIGGER"));
0131
0132         const auto is_endat = static_cast<uint8_t>(packet->lValue(t, "IS_ENDAT"));
0133         if ((is_lvl1 || is_endat) && !bco_set)
0134         {
0135           bco_set = true;
0136           gtm_bco = packet->lValue(t, "BCO");
0137           if (largest_bco < gtm_bco)
0138           {
0139             largest_bco = gtm_bco;
0140           }
0141           if (gtm_bco < minBCO)
0142           {
0143             continue;
0144           }
0145           if (Verbosity() > 0)
0146           {
0147             std::cout << "bco: 0x" << std::hex << gtm_bco << std::dec << std::endl;
0148           }
0149           // store
0150           skipthis = false;
0151           previous_bco = gtm_bco;
0152           if (m_BclkStackPacketMap.find(packet_id) == m_BclkStackPacketMap.end())
0153           {
0154             m_BclkStackPacketMap.insert(std::make_pair(packet_id, std::set<uint64_t>()));
0155           }
0156           m_BclkStackPacketMap[packet_id].insert(gtm_bco);
0157         }
0158       }
0159       if (skipthis)
0160       {
0161         if (Verbosity() > 1)
0162         {
0163           std::cout << "Largest bco: 0x" << std::hex << largest_bco << ", minbco 0x"
0164                     << minBCO << std::dec << ", evtno: " << EventSequence << std::endl;
0165         }
0166       }
0167       else
0168       {
0169         int m_nWaveFormInFrame = packet->iValue(0, "NR_WF");
0170         static int once = 0;
0171         for (int wf = 0; wf < m_nWaveFormInFrame; wf++)
0172         {
0173           if (m_TpcRawHitMap[gtm_bco].size() > 20000)
0174           {
0175             if (!once)
0176             {
0177               std::cout << "too many hits" << std::endl;
0178             }
0179             once++;
0180             continue;
0181           }
0182           else
0183           {
0184             if (once)
0185             {
0186               std::cout << "many more hits: " << once << std::endl;
0187             }
0188             once = 0;
0189           }
0190           bool checksumerror = (packet->iValue(wf, "CHECKSUMERROR") > 0);
0191           if (checksumerror)
0192           {
0193             continue;
0194           }
0195           bool parityerror = (packet->iValue(wf, "DATAPARITYERROR") > 0);
0196           auto newhit = std::make_unique<TpcRawHitv2>();
0197           int FEE = packet->iValue(wf, "FEE");
0198           newhit->set_bco(packet->iValue(wf, "BCO"));
0199
0200           // store gtm bco in hit
0201           newhit->set_gtm_bco(gtm_bco);
0202
0203           newhit->set_packetid(packet->getIdentifier());
0204           newhit->set_fee(FEE);
0205           newhit->set_channel(packet->iValue(wf, "CHANNEL"));
0206           newhit->set_sampaaddress(packet->iValue(wf, "SAMPAADDRESS"));
0207           newhit->set_sampachannel(packet->iValue(wf, "SAMPACHANNEL"));
0208           newhit->set_type(packet->iValue(wf, "TYPE"));
0209           newhit->set_userword(packet->iValue(wf, "USERWORD"));
0210           newhit->set_parity(packet->iValue(wf, "DATAPARITY"));
0211           newhit->set_checksum(packet->iValue(wf, "CHECKSUM"));
0212           newhit->set_checksumerror(checksumerror);
0213           newhit->set_parityerror(parityerror);
0214           //         // checksum and checksum error
0215           //         newhit->set_checksum( packet->iValue(iwf, "CHECKSUM") );
0216           //         newhit->set_checksum_error( packet->iValue(iwf, "CHECKSUMERROR") );
0217
0218           // samples
0219           // const uint16_t samples = packet->iValue(wf, "SAMPLES");
0220
0221           // Temp remedy as we set the time window as 425 for now (extended from previous 360
0222           // due to including of diffused laser flush)
0223           const uint16_t samples = m_max_tpc_time_samples;
0224
0225           newhit->set_samples(samples);
0226
0227           // adc values
0228           for (uint16_t is = 0; is < samples; ++is)
0229           {
0230             uint16_t adval = packet->iValue(wf, is);
0231
0232             // This is temporary fix for decoder change. Will be changed again for real ZS data decoding.
0233             // if(adval >= 64000){ newhit->set_samples(is); break;}
0234
0235             // With this, the hit is unseen from clusterizer
0236             if (adval <= 64000)
0237             {
0238               newhit->set_adc(is, adval);
0239             }
0240           }
0241
0242           m_BeamClockFEE[gtm_bco].insert(FEE);
0243           m_FEEBclkMap[FEE] = gtm_bco;
0244           if (Verbosity() > 2)
0245           {
0246             std::cout << "evtno: " << EventSequence
0247                       << ", hits: " << wf
0248                       << ", num waveforms: " << m_nWaveFormInFrame
0249                       << ", bco: 0x" << std::hex << gtm_bco << std::dec
0250                       << ", FEE: " << FEE << std::endl;
0251           }
0252           //          packet->convert();
0253           // if (m_TpcRawHitMap[gtm_bco].size() < 50000)
0254           // {
0255           if (StreamingInputManager())
0256           {
0257             StreamingInputManager()->AddTpcRawHit(gtm_bco, newhit.get());
0258           }
0259           m_TpcRawHitMap[gtm_bco].push_back(newhit.release());
0260           m_BclkStack.insert(gtm_bco);
0261           //    }
0262         }
0263       }
0264       delete packet;
0265     }
0266   }
0267
0268   //    Print("HITS");
0269   //  } while (m_TpcRawHitMap.size() < 10 || CheckPoolDepth(m_TpcRawHitMap.begin()->first));
0270 }
0271
0272 void SingleTpcPoolInput::Print(const std::string &what) const
0273 {
0274   if (what == "ALL" || what == "FEE")
0275   {
0276     for (const auto &bcliter : m_BeamClockFEE)
0277     {
0278       std::cout << "Beam clock 0x" << std::hex << bcliter.first << std::dec << std::endl;
0279       for (auto feeiter : bcliter.second)
0280       {
0281         std::cout << "FEM: " << feeiter << std::endl;
0282       }
0283     }
0284   }
0285   if (what == "ALL" || what == "FEEBCLK")
0286   {
0287     for (auto bcliter : m_FEEBclkMap)
0288     {
0289       std::cout << "FEE" << bcliter.first << " bclk: 0x"
0290                 << std::hex << bcliter.second << std::dec << std::endl;
0291     }
0292   }
0293   if (what == "ALL" || what == "HITS")
0294   {
0295     const auto bcliter = m_TpcRawHitMap.begin();
0296     {
0297       std::cout << Name() << ": Beam clock 0x" << std::hex << bcliter->first
0298                 << std::dec << ", Number of hits: " << bcliter->second.size()
0299                 << std::endl;
0300     }
0301   }
0302   if (what == "ALL" || what == "STORAGE")
0303   {
0304     for (const auto &bcliter : m_TpcRawHitMap)
0305     {
0306       std::cout << "Beam clock 0x" << std::hex << bcliter.first << std::dec << std::endl;
0307       for (auto feeiter : bcliter.second)
0308       {
0309         std::cout << "fee: " << feeiter->get_fee()
0310                   << " at " << std::hex << feeiter << std::dec << std::endl;
0311       }
0312     }
0313   }
0314   if (what == "ALL" || what == "STACK")
0315   {
0316     for (auto iter : m_BclkStack)
0317     {
0318       std::cout << "stacked bclk: 0x" << std::hex << iter << std::dec << std::endl;
0319     }
0320   }
0321 }
0322
0323 void SingleTpcPoolInput::CleanupUsedPackets(const uint64_t bclk)
0324 {
0325   if (Verbosity() > 2)
0326   {
0327     std::cout << "cleaning up bcos < 0x" << std::hex
0328               << bclk << std::dec << std::endl;
0329   }
0330   std::vector<uint64_t> toclearbclk;
0331   for (const auto &iter : m_TpcRawHitMap)
0332   {
0333     if (iter.first <= bclk)
0334     {
0335       for (auto pktiter : iter.second)
0336       {
0337         delete pktiter;
0338       }
0339       toclearbclk.push_back(iter.first);
0340     }
0341     else
0342     {
0343       break;
0344     }
0345   }
0346   // for (auto iter :  m_BeamClockFEE)
0347   // {
0348   //   iter.second.clear();
0349   // }
0350   for (auto iter : toclearbclk)
0351   {
0352     m_BclkStack.erase(iter);
0353     m_BeamClockFEE.erase(iter);
0354     m_TpcRawHitMap.erase(iter);
0355   }
0356 }
0357
0358 bool SingleTpcPoolInput::CheckPoolDepth(const uint64_t bclk)
0359 {
0360   // if (m_FEEBclkMap.size() < 10)
0361   // {
0362   //   std::cout << "not all FEEs in map: " << m_FEEBclkMap.size() << std::endl;
0363   //   return true;
0364   // }
0365   for (auto iter : m_FEEBclkMap)
0366   {
0367     if (Verbosity() > 2)
0368     {
0369       std::cout << "my bclk 0x" << std::hex << iter.second
0370                 << " req: 0x" << bclk << std::dec << std::endl;
0371     }
0372     if (iter.second < bclk)
0373     {
0374       if (Verbosity() > 1)
0375       {
0376         std::cout << "FEE " << iter.first << " beamclock 0x" << std::hex << iter.second
0377                   << " smaller than req bclk: 0x" << bclk << std::dec << std::endl;
0378       }
0379       return false;
0380     }
0381   }
0382   return true;
0383 }
0384
0385 void SingleTpcPoolInput::ClearCurrentEvent()
0386 {
0387   // called interactively, to get rid of the current event
0388   uint64_t currentbclk = *m_BclkStack.begin();
0389   //  std::cout << "clearing bclk 0x" << std::hex << currentbclk << std::dec << std::endl;
0390   CleanupUsedPackets(currentbclk);
0391   // m_BclkStack.erase(currentbclk);
0392   // m_BeamClockFEE.erase(currentbclk);
0393   return;
0394 }
0395
0396 bool SingleTpcPoolInput::GetSomeMoreEvents(const uint64_t ibclk)
0397 {
0398   if (AllDone())
0399   {
0400     return false;
0401   }
0402   if (m_TpcRawHitMap.empty())
0403   {
0404     return true;
0405   }
0406   uint64_t localbclk = ibclk;
0407   if (ibclk == 0)
0408   {
0409     if (m_TpcRawHitMap.empty())
0410     {
0411       return true;
0412     }
0413     localbclk = m_TpcRawHitMap.begin()->first;
0414   }
0415
0416   std::set<int> toerase;
0417   for (auto bcliter : m_FEEBclkMap)
0418   {
0419     if (bcliter.second <= localbclk)
0420     {
0421       uint64_t highest_bclk = m_TpcRawHitMap.rbegin()->first;
0422       if ((highest_bclk - m_TpcRawHitMap.begin()->first) < MaxBclkDiff())
0423       {
0424         // std::cout << "FEE " << bcliter.first << " bclk: "
0425         //      << std::hex << bcliter.second << ", req: " << lowest_bclk
0426         //       << " low: 0x" <<  m_TpcRawHitMap.begin()->first << ", high: " << highest_bclk << ", delta: " << std::dec << (highest_bclk-m_TpcRawHitMap.begin()->first)
0427         //      << std::dec << std::endl;
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_TpcRawHitMap.begin()->first
0435                   << ", to: 0x" << highest_bclk << ", delta: " << std::dec
0436                   << (highest_bclk - m_TpcRawHitMap.begin()->first)
0437                   << std::dec << std::endl;
0438         toerase.insert(bcliter.first);
0439       }
0440     }
0441   }
0442   for (auto iter : toerase)
0443   {
0444     m_FEEBclkMap.erase(iter);
0445   }
0446   return false;
0447   // if (CheckPoolDepth(m_TpcRawHitMap.begin()->first))
0448   // {
0449   //   if (m_TpcRawHitMap.size() >= 10)
0450   //   {
0451   //     return false;
0452   //   }
0453   // }
0454   // return true;
0455 }
0456
0457 void SingleTpcPoolInput::CreateDSTNode(PHCompositeNode *topNode)
0458 {
0459   PHNodeIterator iter(topNode);
0460   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0461   if (!dstNode)
0462   {
0463     dstNode = new PHCompositeNode("DST");
0464     topNode->addNode(dstNode);
0465   }
0466   PHNodeIterator iterDst(dstNode);
0467   PHCompositeNode *detNode = dynamic_cast<PHCompositeNode *>(iterDst.findFirst("PHCompositeNode", "TPC"));
0468   if (!detNode)
0469   {
0470     detNode = new PHCompositeNode("TPC");
0471     dstNode->addNode(detNode);
0472   }
0473   TpcRawHitContainer *tpchitcont = findNode::getClass<TpcRawHitContainer>(detNode, m_rawHitContainerName);
0474   if (!tpchitcont)
0475   {
0476     tpchitcont = new TpcRawHitContainerv2();
0477     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(tpchitcont, m_rawHitContainerName, "PHObject");
0478     detNode->addNode(newNode);
0479   }
0480 }
0481
0482 void SingleTpcPoolInput::ConfigureStreamingInputManager()
0483 {
0484   if (StreamingInputManager())
0485   {
0486     StreamingInputManager()->SetTpcBcoRange(m_BcoRange);
0487     StreamingInputManager()->SetTpcNegativeBco(m_NegativeBco);
0488   }
0489   return;
0490 }