RceRawDecoder_module.cc

Go to the documentation of this file.

// art includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"

// artdaq and dune-raw-data includes
#include "dune-raw-data/Overlays/RceFragment.hh"
#include "artdaq-core/Data/Fragment.hh"
#include "artdaq-core/Data/ContainerFragment.hh"
#include "dune-raw-data/Overlays/FragmentType.hh"
#include "dune-raw-data/Services/ChannelMap/PdspChannelMapService.h"

// larsoft includes
#include "lardataobj/RawData/RawDigit.h"

// ROOT includes
#include "TH1.h"
#include "TH2F.h"
#include "TProfile.h"

// C++ Includes
#include <memory>
#include <iostream>

namespace dune {
  class RceRawDecoder;
}

class dune::RceRawDecoder : public art::EDProducer {
public:
  explicit RceRawDecoder(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  RceRawDecoder(RceRawDecoder const &) = delete;
  RceRawDecoder(RceRawDecoder &&) = delete;
  RceRawDecoder & operator = (RceRawDecoder const &) = delete;
  RceRawDecoder & operator = (RceRawDecoder &&) = delete;

  // Required functions.
  void produce(art::Event & e) override;
  void reconfigure(const fhicl::ParameterSet &pset);
  void beginJob() override;

  void setRootObjects();

private:
  typedef std::vector<raw::RawDigit> RawDigits;

  bool _process(
          const artdaq::Fragment& frag, 
          RawDigits& raw_digits);

  // Declare member data here.
  std::string _input_label; 
  std::string _output_label;
  bool _expect_container_fragments;

  TH1D* _h_nticks;
  TH1I* _h_all_adc_values;
  TH1I* _h_crate_numbers;
  TH1I* _h_slot_IDs;
  TH1I* _h_fiber_IDs;
  TH1I* _h_online_channels;
  TH1I* _h_offline_channels;
  std::vector<TH2F*> _h_mean_slot_channels;
  std::vector<TH2F*> _h_rms_slot_channels;      
  std::vector<TProfile*> _h_mean_slot_channels_pfx;
  std::vector<TProfile*> _h_rms_slot_channels_pfx;      
  std::vector<TH2I*> _h_fiber_persistent_wav;
  std::vector<int16_t> _buffer;
        
  // define my rms and mean functions
  float rmsADC(raw::RawDigit::ADCvector_t & wav);
  float meanADC(raw::RawDigit::ADCvector_t & wav);
};


dune::RceRawDecoder::RceRawDecoder(fhicl::ParameterSet const & pset)
  : EDProducer{pset}
// Initialize member data here.
{
  art::ServiceHandle<art::TFileService> fs;
  //fs->registerFileSwitchCallback(this, &RceRawDecoder::setRootObjects);
  setRootObjects();

  reconfigure(pset);
  // Call appropriate produces<>() functions here.
  //produces< std::vector<raw::RawDigit> > ( _output_label );  
  produces<RawDigits>( _output_label );  
}

void dune::RceRawDecoder::reconfigure(fhicl::ParameterSet const& pset) {

  _input_label = pset.get<std::string>("RawDataLabel");
  _output_label = pset.get<std::string>("OutputDataLabel");
  _expect_container_fragments = pset.get<bool>("ExpectContainerFragments", true);
}

void dune::RceRawDecoder::setRootObjects(){
        //clear vectors
        _h_mean_slot_channels.clear();
        _h_rms_slot_channels.clear();
        _h_mean_slot_channels_pfx.clear();
        _h_rms_slot_channels_pfx.clear();
        _h_fiber_persistent_wav.clear();
        // place to define the histograms
  art::ServiceHandle<art::TFileService> file_srv;
  _h_nticks = file_srv->make<TH1D>("rce_NTicks","TPC: Number of ticks",  100, 0, 20000);
  _h_nticks->SetXTitle("NTicks");
  _h_all_adc_values = file_srv->make<TH1I>("rce_All_ADC_values","TPC: All_ADC_values",  5000, 0, 5000);
  _h_all_adc_values->SetXTitle("ADC");
  _h_crate_numbers = file_srv->make<TH1I>("rce_All_crate_numbers","TPC: Crate Numbers",  6, 0, 6);
  _h_crate_numbers->SetXTitle("Crate Number");
  _h_slot_IDs = file_srv->make<TH1I>("rce_All_slot_IDs","TPC: Slot IDs",  30, 0, 30);
  _h_slot_IDs->SetXTitle("Slot ID");
  _h_fiber_IDs = file_srv->make<TH1I>("rce_All_fiber_IDs","TPC: Fiber IDs",  120, 0, 120);
  _h_fiber_IDs->SetXTitle("Fiber ID");
  _h_online_channels = file_srv->make<TH1I>("rce_Online_Channels", "TPC: Online Channel ID", 5000, 0, 5000);
  _h_offline_channels = file_srv->make<TH1I>("rce_Offline_Channels", "TPC: Offline Channel ID", 5000, 0, 5000);
  for(int i=0;i<30;i++) {
        _h_mean_slot_channels.push_back(file_srv->make<TH2F>(Form("Slot%d_Mean", i), Form("Slot%d:Mean_vs_SlotChannel", i), 512, 0, 512, 5000, .0, 5000)); //hard-coded for now. Will be moved to the analyzer module. Feb 7, 2018
        _h_rms_slot_channels.push_back(file_srv->make<TH2F>(Form("Slot%d_RMS", i), Form("Slot%d:RMS_vs_SlotChannel", i), 512, 0, 512, 5000, .0, 5000)); //hard-coded for now. Will be moved to the analyzer module. Feb 7, 2018
        _h_mean_slot_channels_pfx.push_back(file_srv->make<TProfile>(Form("Slot%d_Mean_pfx", i), Form("Slot%d:Mean_vs_SlotChannel_pfx", i), 512, 0, 512)); //hard-coded for now. Will be moved to the analyzer module. Feb 7, 2018
        _h_rms_slot_channels_pfx.push_back(file_srv->make<TProfile>(Form("Slot%d_RMS_pfx", i), Form("Slot%d:RMS_vs_SlotChannel_pfx", i), 512, 0, 512)); //hard-coded for now. Will be moved to the analyzer module. Feb 7, 2018
        
        _h_mean_slot_channels[i]->GetXaxis()->SetTitle("Slot Channel"); _h_mean_slot_channels[i]->GetYaxis()->SetTitle("Raw Mean"); 
        _h_rms_slot_channels[i]->GetXaxis()->SetTitle("Slot Channel"); _h_rms_slot_channels[i]->GetYaxis()->SetTitle("Raw RMS"); 
        _h_mean_slot_channels_pfx[i]->GetXaxis()->SetTitle("Slot Channel"); _h_mean_slot_channels_pfx[i]->GetYaxis()->SetTitle("Profiled Mean"); 
        _h_rms_slot_channels_pfx[i]->GetXaxis()->SetTitle("Slot Channel"); _h_rms_slot_channels_pfx[i]->GetYaxis()->SetTitle("Profiled RMS"); 
  }
  for(int i=0;i<120;i++) {
    _h_fiber_persistent_wav.push_back(file_srv->make<TH2I>(Form("Persistent_Waveform_Fiber#%d", i), Form("Persistent_Waveform_Fiber#%d", i), 10000, 0, 10000, 500, 0, 5000));
  }
  
}
void dune::RceRawDecoder::beginJob(){
}

void dune::RceRawDecoder::produce(art::Event & evt){
    // TODO Use MF_LOG_DEBUG
    MF_LOG_INFO("RceRawDecoder")
      << "-------------------- RCE RawDecoder -------------------";

  RawDigits raw_digits;
  unsigned int n_rce_frags = 0;

  if (_expect_container_fragments) {
    art::InputTag itag1(_input_label, "ContainerTPC");
    auto cont_frags = evt.getHandle<artdaq::Fragments>(itag1);
    art::EventNumber_t eventNumber = evt.event();
    // Check if there is Timing data in this event
    // Don't crash code if not present, just don't save anything
    try { cont_frags->size(); }
    catch(std::exception const&) {
      std::cout << "WARNING: Container TPC/RCE data not found in event " << eventNumber << std::endl;
      std::vector<raw::RawDigit> digits;
      evt.put(std::make_unique<std::vector<raw::RawDigit>>(std::move(digits)), _output_label);
      return;
    }
    //Check that the data is valid
    if(!cont_frags){
      MF_LOG_ERROR("RceRawDecoder")
          << "Run: " << evt.run()
                  << ", SubRun: " << evt.subRun()
          << ", Event: " << evt.event()
                  << " Container Fragments is NOT VALID";
    }
    
    for (auto const& cont : *cont_frags)
    {
      artdaq::ContainerFragment cont_frag(cont);
      for (size_t ii = 0; ii < cont_frag.block_count(); ++ii)
        {
          //artdaq::Fragment frag;
          //size_t frag_size = cont_frag.fragSize(ii);
          //frag.resizeBytes(frag_size);
          //memcpy(frag.headerAddress(), cont_frag.at(ii), frag_size);
          if (_process(*cont_frag[ii], raw_digits)) ++n_rce_frags;
        }
    }
  }
  else
  {
    art::InputTag itag2(_input_label, "TPC");
    auto frags = evt.getHandle<artdaq::Fragments>(itag2);

    // Check if there is Timing data in this event
    // Don't crash code if not present, just don't save anything
    art::EventNumber_t eventNumber = evt.event();
    try { frags->size(); }
    catch(std::exception const&) {
      std::cout << "WARNING: Raw TPC/RCE data not found in event " << eventNumber << std::endl;
      std::vector<raw::RawDigit> digits;
      evt.put(std::make_unique<std::vector<raw::RawDigit>>(std::move(digits)), _output_label);
      return;
    }

    //Check that the data is valid
    if(!frags.isValid()){
      MF_LOG_ERROR("RceRawDecoder")
          << "Run: " << evt.run()
                  << ", SubRun: " << evt.subRun()
          << ", Event: " << evt.event()
                  << " Fragments is NOT VALID";
    }

    for(auto const& frag: *frags)
    {
      if (_process(frag, raw_digits)) ++n_rce_frags;
    }
  }

  MF_LOG_INFO("RceRawDecoder")
      << " Processed " << n_rce_frags
      << " RCE Fragments, "
      << raw_digits.size()
      << " RawDigits.";

  evt.put(std::make_unique<decltype(raw_digits)>(std::move(raw_digits)),
          _output_label);
}

bool dune::RceRawDecoder::_process(
        const artdaq::Fragment& frag, 
        RawDigits& raw_digits
        )
{
  // FIXME: Remove hard-coded fragment type
  if((unsigned)frag.type() != 2) return false;

  MF_LOG_INFO("RceRawDecoder")
      << "   SequenceID = " << frag.sequenceID()
      << "   fragmentID = " << frag.fragmentID()
      << "   fragmentType = " << (unsigned)frag.type()
      << "   Timestamp =  " << frag.timestamp();
  art::ServiceHandle<dune::PdspChannelMapService> channelMap;
  dune::RceFragment rce(frag);
  
  uint32_t ch_counter = 0;
  for (int i = 0; i < rce.size(); ++i)
  {
     auto const * rce_stream = rce.get_stream(i);
     int n_ch = rce_stream->getNChannels();
     int n_ticks = rce_stream->getNTicks();
     auto const identifier = rce_stream->getIdentifier();
     uint32_t crateNumber = identifier.getCrate();
     uint32_t slotNumber = identifier.getSlot();
     uint32_t fiberNumber = identifier.getFiber();
     std::cout<<"crate, slot, fiber = "<<crateNumber<<", "<<slotNumber<<", "<<fiberNumber<<std::endl;
     uint32_t fiberID = (crateNumber*5+slotNumber)*4 + fiberNumber;
     _h_fiber_IDs->Fill(fiberID);
     MF_LOG_INFO("RceRawDecoder")
         << "RceFragment timestamp: " << rce_stream->getTimeStamp()
         << ", NChannels: " << n_ch
         << ", NTicks: " << n_ticks;

     _h_nticks->Fill(n_ticks);

     size_t buffer_size = n_ch * n_ticks;
     if (_buffer.capacity() < buffer_size)
     {
         MF_LOG_INFO("RceRawDecoder")
             << "Increase buffer size from " << _buffer.capacity()
             << " to " << buffer_size;

         _buffer.reserve(buffer_size);
     }

     int16_t* adcs = _buffer.data();
     rce_stream->getMultiChannelData(adcs);

     raw::RawDigit::ADCvector_t v_adc;
     for (int i_ch = 0; i_ch < n_ch; i_ch++)
     {
        if(i==0 && i_ch ==0) std::cout<<" ADCs for the the 100 ticks in the 1st channel of the 1st RCE "<<std::endl;
        v_adc.clear();
        for (int i_tick = 0; i_tick < n_ticks; i_tick++)
        {
            //print ADCs for the the 100 ticks in the 1st channel of the 1st RCE
            if(i==0 && i_ch==0 && i_tick<100) {
              std::cout<<adcs[i_tick]<<"\t";
              if(i_tick==99) std::cout<<std::endl;
            }
            v_adc.push_back(adcs[i_tick]);
            _h_all_adc_values->Fill(adcs[i_tick]);
            //save the 1st waveform in a fiber
            _h_fiber_persistent_wav.at(fiberID)->Fill(i_tick, adcs[i_tick]);
        }
        adcs += n_ticks;

        ch_counter++;
        int offlineChannel = -1;
        offlineChannel = channelMap->GetOfflineNumberFromDetectorElements(crateNumber, slotNumber, fiberNumber, i_ch, dune::PdspChannelMapService::kRCE);
        _h_offline_channels->Fill(offlineChannel);
        raw::RawDigit raw_digit(offlineChannel, n_ticks, v_adc);
        raw_digits.push_back(raw_digit);
        
        //=========== fill Mean and RMS with slotchannel number=====================
        uint32_t slotID = crateNumber*5 + slotNumber; //0 - 29
        _h_slot_IDs->Fill(slotID);
        uint32_t slotchannel = 0; // 0 - 127
        slotchannel = fiberNumber*128 + i_ch; //hard-coded. will be moved to analyzer module. Feb 7, 2018
        float mean = meanADC(v_adc);
        float rms = rmsADC(v_adc);
        _h_mean_slot_channels.at(slotID)->Fill(slotchannel, mean);
        _h_rms_slot_channels.at(slotID)->Fill(slotchannel, rms);
        _h_mean_slot_channels_pfx.at(slotID)->Fill(slotchannel, mean, 1);
        _h_rms_slot_channels_pfx.at(slotID)->Fill(slotchannel, rms, 1);
        //==========================================================================
        
        
     }
  }

  return true;
}

//-----------------------------------------------------------------------   
  // define RMS
  float dune::RceRawDecoder::rmsADC(raw::RawDigit::ADCvector_t &wav)
  {
    int n = wav.size();
    float sum = 0.;
    for(int i = 0; i < n; i++){
      if(wav[i]!=0) sum += wav[i];
    }
    float mean = sum / n;
    sum = 0;
    for(int i = 0; i < n; i++)
      {
        if (wav[i]!=0)     sum += (wav[i]-mean)*(wav[i]-mean);
      }
    return sqrt(sum / n);
  }

  //-----------------------------------------------------------------------  
  //define Mean
  float dune::RceRawDecoder::meanADC(raw::RawDigit::ADCvector_t &wav)
  {
    int n = wav.size();
    float sum = 0.;
    for(int i = 0; i < n; i++)
      {
        if (wav[i]!=0) sum += abs(wav[i]);
      }
    return sum / n;
  }

DEFINE_ART_MODULE(dune::RceRawDecoder)