packages/micromegas/MicromegasCombinedDataDecoder.cc

0001 /*!
0002  * \file MicromegasCombinedDataDecoder.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "MicromegasCombinedDataDecoder.h"
0007 #include "MicromegasDefs.h"
0008
0009 #include <trackbase/TrkrHitSet.h>
0010 #include <trackbase/TrkrHitSetContainerv1.h>
0011 #include <trackbase/TrkrHitv2.h>
0012
0013 #include <ffarawobjects/MicromegasRawHit.h>
0014 #include <ffarawobjects/MicromegasRawHitContainer.h>
0015
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <fun4all/Fun4AllServer.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/getClass.h>
0022
0023 #include <algorithm>
0024 #include <cassert>
0025 #include <memory>
0026
0027 namespace
0028 {
0029   /*
0030    * returns true if a given channel for a given FEE is permanently masked
0031    * for now all channels from 128 to 255, for FEE 8 (SCOZ) are masked
0032    */
0033   bool channel_is_permanently_masked( int fee_id, int channel )
0034   { return fee_id==8 && channel>=128; }
0035
0036 }
0037
0038 //_________________________________________________________
0039 MicromegasCombinedDataDecoder::MicromegasCombinedDataDecoder(const std::string& name)
0040   : SubsysReco(name)
0041 {
0042 }
0043
0044 //_____________________________________________________________________
0045 int MicromegasCombinedDataDecoder::Init(PHCompositeNode* /*topNode*/)
0046 {
0047   // print configuration
0048   std::cout
0049       << "MicromegasCombinedDataDecoder::Init -"
0050       << " m_calibration_filename: "
0051       << (m_calibration_filename.empty() ? "unspecified" : m_calibration_filename)
0052       << std::endl;
0053
0054   std::cout
0055       << "MicromegasCombinedDataDecoder::Init -"
0056       << " m_hot_channel_map_filename: "
0057       << (m_hot_channel_map_filename.empty() ? "unspecified" : m_hot_channel_map_filename)
0058       << std::endl;
0059
0060   std::cout << "MicromegasCombinedDataDecoder::Init - m_n_sigma: " << m_n_sigma << std::endl;
0061   std::cout << "MicromegasCombinedDataDecoder::Init - m_min_adc: " << m_min_adc << std::endl;
0062   std::cout << "MicromegasCombinedDataDecoder::Init - m_sample_min: " << m_sample_min << std::endl;
0063   std::cout << "MicromegasCombinedDataDecoder::Init - m_sample_max: " << m_sample_max << std::endl;
0064
0065   // read calibrations
0066   m_calibration_data.read(m_calibration_filename);
0067
0068   // read hot channels
0069   if (!m_hot_channel_map_filename.empty())
0070   {
0071     m_hot_channels.read(m_hot_channel_map_filename);
0072   }
0073
0074   return Fun4AllReturnCodes::EVENT_OK;
0075 }
0076
0077 //____________________________________________________________________________..
0078 int MicromegasCombinedDataDecoder::InitRun(PHCompositeNode* topNode)
0079 {
0080   // get dst node
0081   PHNodeIterator iter(topNode);
0082   auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0083   if (!dstNode)
0084   {
0085     std::cout << "MicromegasCombinedDataDecoder::InitRun - DST Node missing, doing nothing." << std::endl;
0086     exit(1);
0087   }
0088
0089   // create hitset container if needed
0090   auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0091   if (!hitsetcontainer)
0092   {
0093     std::cout << "MicromegasCombinedDataDecoder::InitRun - creating TRKR_HITSET." << std::endl;
0094
0095     // find or create TRKR node
0096     PHNodeIterator dstiter(dstNode);
0097     auto trkrnode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0098     if (!trkrnode)
0099     {
0100       trkrnode = new PHCompositeNode("TRKR");
0101       dstNode->addNode(trkrnode);
0102     }
0103
0104     // create container and add to the tree
0105     hitsetcontainer = new TrkrHitSetContainerv1;
0106     auto newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
0107     trkrnode->addNode(newNode);
0108   }
0109   auto rawhitcontainer = findNode::getClass<MicromegasRawHitContainer>(topNode, m_rawhitnodename);
0110   if(!rawhitcontainer)
0111   {
0112     Fun4AllServer* se = Fun4AllServer::instance();
0113     se->unregisterSubsystem(this);
0114     std::cout << PHWHERE << "Removing TPOT unpacker, no raw hit container" << std::endl;
0115   }
0116
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 //___________________________________________________________________________
0121 int MicromegasCombinedDataDecoder::process_event(PHCompositeNode* topNode)
0122 {
0123   // load relevant nodes
0124   // Get the TrkrHitSetContainer node
0125   auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0126   assert(trkrhitsetcontainer);
0127
0128   // load raw hits container
0129   auto rawhitcontainer = findNode::getClass<MicromegasRawHitContainer>(topNode, m_rawhitnodename);
0130   assert(rawhitcontainer);
0131
0132   // loop over raw hits
0133   if (Verbosity())
0134   {
0135     std::cout << "MicromegasCombinedDataDecoder::process_event - hits: " << rawhitcontainer->get_nhits() << std::endl;
0136   }
0137
0138   int n_signal_hits = 0;
0139
0140   bool first = true;
0141   uint64_t first_lvl1_bco = 0;
0142
0143   for (unsigned int ihit = 0; ihit < rawhitcontainer->get_nhits(); ++ihit)
0144   {
0145     const auto rawhit = rawhitcontainer->get_hit(ihit);
0146     const auto packet_id = rawhit->get_packetid();
0147
0148     if (first)
0149     {
0150       first = false;
0151       first_lvl1_bco = rawhit->get_gtm_bco();
0152     }
0153
0154     // make sure packet is valid
0155     if (std::find(std::begin(MicromegasDefs::m_packet_ids), std::end(MicromegasDefs::m_packet_ids), packet_id) == std::end(MicromegasDefs::m_packet_ids))
0156     {
0157       std::cout << "MicromegasCombinedDataDecoder::process_event - invalid packet: " << packet_id << std::endl;
0158       continue;
0159     }
0160
0161     // get fee id, apply mapping to current fiber set, for backward compatibility
0162     const int fee = m_mapping.get_new_fee_id(rawhit->get_fee());
0163     const auto channel = rawhit->get_channel();
0164
0165     // check if channel is permanently masked
0166     if( channel_is_permanently_masked(fee, channel ))
0167     {
0168       continue;
0169     }
0170
0171     // map fee and channel to physical hitsetid and physical strip
0172     // get hitset key matching this fee
0173     const TrkrDefs::hitsetkey hitsetkey = m_mapping.get_hitsetkey(fee);
0174     if (!hitsetkey)
0175     {
0176       continue;
0177     }
0178
0179     // get matching layer, tile, physical strip
0180     int layer = int(TrkrDefs::getLayer(hitsetkey));
0181     int tile = int(MicromegasDefs::getTileId(hitsetkey));
0182     int strip = m_mapping.get_physical_strip(fee, channel);
0183     if (strip < 0)
0184     {
0185       continue;
0186     }
0187
0188     // check agains hot channels
0189     if (m_hot_channels.is_hot_channel(layer, tile, strip))
0190     {
0191       continue;
0192     }
0193
0194     // get channel rms and pedestal from calibration data
0195     const double pedestal = m_calibration_data.get_pedestal(fee, channel);
0196     const double rms = m_calibration_data.get_rms(fee, channel);
0197
0198     // a rms of zero means the calibration has failed. the data is unusable
0199     if (rms <= 0)
0200     {
0201       continue;
0202     }
0203
0204     // loop over sample_range find maximum
0205     const auto sample_range = std::make_pair(rawhit->get_sample_begin(), rawhit->get_sample_end());
0206     std::vector<uint16_t> adc_list;
0207     for (auto is = std::max(m_sample_min, sample_range.first); is < std::min(m_sample_max, sample_range.second); ++is)
0208     {
0209       const uint16_t adc = rawhit->get_adc(is);
0210       if (adc != MicromegasDefs::m_adc_invalid)
0211       {
0212         adc_list.push_back(adc);
0213       }
0214     }
0215
0216     if (adc_list.empty())
0217     {
0218       continue;
0219     }
0220
0221     // get max adc value in range
0222     /* TODO: use more advanced signal processing */
0223     auto max_adc = *std::max_element(adc_list.begin(), adc_list.end());
0224
0225     // compare to hard min_adc value
0226     if (max_adc < m_min_adc)
0227     {
0228       continue;
0229     }
0230
0231     // compare to threshold
0232     if (max_adc < pedestal + m_n_sigma * rms)
0233     {
0234       continue;
0235     }
0236
0237     if (Verbosity())
0238     {
0239       const auto bco = rawhit->get_gtm_bco();
0240       std::cout << "MicromegasCombinedDataDecoder::process_event -"
0241                 << " bco: " << bco
0242                 << " layer: " << layer
0243                 << " tile: " << tile
0244                 << " channel: " << channel
0245                 << " strip: " << strip
0246                 << " adc: " << max_adc
0247                 << std::endl;
0248     }
0249
0250     // get matching hitset
0251     const auto hitset_it = trkrhitsetcontainer->findOrAddHitSet(hitsetkey);
0252
0253     // generate hit key
0254     const TrkrDefs::hitkey hitkey = MicromegasDefs::genHitKey(strip);
0255
0256     // find existing hit, or create
0257     auto hit = hitset_it->second->getHit(hitkey);
0258     if (hit)
0259     {
0260       // std::cout << "MicromegasCombinedDataDecoder::process_event - duplicated hit, hitsetkey: " << hitsetkey << " strip: " << strip << std::endl;
0261       continue;
0262     }
0263
0264     // create hit, assign adc and insert in hitset
0265     hit = new TrkrHitv2;
0266     hit->setAdc(max_adc);
0267     hitset_it->second->addHitSpecificKey(hitkey, hit);
0268
0269     // increment counter
0270     ++m_hitcounts[hitsetkey];
0271     ++n_signal_hits;
0272   }
0273
0274   if (Verbosity())
0275   {
0276     std::cout << "MicromegasCombinedDataDecoder::process_event - BCO: " << first_lvl1_bco << " n_signal_hits: " << n_signal_hits << std::endl;
0277   }
0278
0279   return Fun4AllReturnCodes::EVENT_OK;
0280 }
0281
0282 //_____________________________________________________________________
0283 int MicromegasCombinedDataDecoder::End(PHCompositeNode* /*topNode*/)
0284 {
0285   // if( Verbosity() )
0286   {
0287     for (const auto& [hitsetkey, count] : m_hitcounts)
0288     {
0289       std::cout << "MicromegasCombinedDataDecoder::End - hitsetkey: " << hitsetkey << ", hit count: " << count << std::endl;
0290     }
0291   }
0292
0293   return Fun4AllReturnCodes::EVENT_OK;
0294 }