FilterAnalyzer_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       FilterAnalyzer
// Module Type: analyser
// File:        FilterAnalyzer_module.cc
// Author:      Karl Warburton (k.warburton@sheffield.ac.uk), April 2016
// 
// A quick analysis module for looking at the filter.
//
// Runs over an artdaq-formatted file and produces a tree or histogram,
// filled for each channel combination and the respective correlation.
////////////////////////////////////////////////////////////////////////

// framework
#include "art/Framework/Core/EDAnalyzer.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 "art_root_io/TFileService.h"
#include "art_root_io/TFileDirectory.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// lbne-artdaq
#include "lbne-raw-data/Overlays/TpcMilliSliceFragment.hh"
#include "artdaq-core/Data/Fragment.hh"

// larsoft
#include "lardataobj/RawData/RawDigit.h"
#include "lardataobj/RawData/raw.h"
#include "larcore/Geometry/Geometry.h"
#include "tpcFragmentToRawDigits.h"
#include "utilities/UnpackFragment.h"
#include "larevt/CalibrationDBI/Interface/DetPedestalService.h"
#include "larevt/CalibrationDBI/Interface/DetPedestalProvider.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
#include "dune/RunHistory/DetPedestalDUNE.h"
#include "cetlib/getenv.h"

// c++
#include <memory>
#include <array>
#include <iostream>
#include <fstream>
#include <sstream>

// ROOT
#include "TMath.h"
#include "TTree.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TVirtualFFT.h"
#include "TStyle.h"

namespace DAQToOffline {
  class FilterAnalyzer;
}

class DAQToOffline::FilterAnalyzer : public art::EDAnalyzer {
public:

  explicit FilterAnalyzer(fhicl::ParameterSet const & pset);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

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

  void analyze(art::Event const& evt) override;
  void reconfigure(const fhicl::ParameterSet &pset);

private:

  std::string fDigitModuleLabel;
  std::string fDigitModuleInstance;
  bool        fMakeADCPlots;

  TF1* fColFilterFunc;  ///< Parameterized collection filter function.
  TF1* fIndUFilterFunc; ///< Parameterized induction filter function.
  TF1* fIndVFilterFunc; ///< Parameterized induction filter function.

  const lariov::DetPedestalProvider& fPedestalRetrievalAlg = *(lar::providerFrom<lariov::DetPedestalService>());

  TH1F* RawDigitHistos[2048];
  TH1F* RawDigitFFT[2048];
  TH1F* RawDigitFFTCorrect[2048];
  TH1F* RawDigitFFTChannel[2048];
  TH1F* RawDigitFFTChannelFilter[2048];
  TH1F* RawDigitInvFFT[2048];

  TH2F* RawFFT_100KHz;
  TH2F* FixFFT_100KHz;
  TH2F* RawFFT_1000KHz;
  TH2F* FixFFT_1000KHz;  
};

DAQToOffline::FilterAnalyzer::FilterAnalyzer(fhicl::ParameterSet const & pset) : art::EDAnalyzer(pset) {

  this->reconfigure(pset);
  gStyle->SetOptStat(0);

  art::ServiceHandle<art::TFileService> tfs;
  RawFFT_100KHz  = tfs->make<TH2F>("RawFFT_100KHz" , "Raw FFT for all channels less than 100 KHz; Channel Number; Frequency (KHz)"             , 2048, 0, 2048, 100, 0, 100 );
  FixFFT_100KHz  = tfs->make<TH2F>("FixFFT_100KHz" , "FFT for all channels less than 100 KHz after filtering; Channel Number; Frequency (KHz)" , 2048, 0, 2048, 100, 0, 100 );
  RawFFT_1000KHz = tfs->make<TH2F>("RawFFT_1000KHz", "Raw FFT for all channels less than 1000 KHz; Channel Number; Frequency (KHz)"            , 2048, 0, 2048, 1000, 0, 1000);
  FixFFT_1000KHz = tfs->make<TH2F>("FixFFT_1000KHz", "FFT for all channels less than 1000 KHz after filtering; Channel Number; Frequency (KHz)", 2048, 0, 2048, 1000, 0, 1000);

  int NumBins = 15000;
  for (int HistoChan=0; HistoChan<2048; HistoChan++) {
    
    // Make the RawDigit Histograms....
    std::stringstream oss1, oss2;
    oss1 << "RawDigit_"<<HistoChan;
    oss2 << "Raw Digit for Channel "<<HistoChan<<"; Time (us); ADC";
    std::string Name1 = oss1.str();
    std::string Title1= oss2.str();
    if (fMakeADCPlots) RawDigitHistos[HistoChan] = tfs->make<TH1F>(Name1.c_str(), Title1.c_str(), NumBins, 0, NumBins/2);
    else RawDigitHistos[HistoChan] = new TH1F(Name1.c_str(), Title1.c_str(), NumBins, 0, NumBins/2);  
    
    // Make the FFT Histograms....
    std::stringstream oss3, oss3a, oss4;
    oss3 << "FFT_"<<HistoChan;
    oss3a << "FFT_"<<HistoChan<<"_correct";
    oss4 << "FFT of Raw Digit for Channel "<<HistoChan<< " in frequency domain";
    std::string Name2 = oss3.str();
    std::string Name2a= oss3a.str();
    std::string Title2= oss4.str();
    RawDigitFFT[HistoChan] = new TH1F(Name2.c_str(), Title2.c_str(), NumBins, 0, NumBins);
    RawDigitFFTCorrect[HistoChan] = new TH1F(Name2a.c_str(),Title2.c_str(), NumBins, 0, NumBins);

    // Make the FFT Channel Histograms....
    std::stringstream oss5, oss5a, oss6;
    oss5 << "ChanFFT_"<<HistoChan;
    oss6 << "FFT of Raw Digit for Channel "<<HistoChan<<"; Frequency (KHz); Number";
    std::string Name3 = oss5.str();
    std::string Title3= oss6.str();
    RawDigitFFTChannel[HistoChan] = new TH1F(Name3.c_str(), Title3.c_str(), NumBins, 0, 2000);

    // Make the FFT Channel after filter Histograms....
    std::stringstream oss9, oss10;
    oss9 << "ChanFFT_"<<HistoChan<<"_Filter";
    oss10 << "FFT of Raw Digit for Channel "<<HistoChan<<" after filter is applied; Frequency (KHz); Number";
    std::string Name5 = oss9.str();
    std::string Title5= oss10.str();
    RawDigitFFTChannelFilter[HistoChan] = new TH1F(Name5.c_str(), Title5.c_str(), NumBins, 0, 2000);
   
    // Make the Inverse FFT Histograms....
    std::stringstream oss7, oss8;
    oss5 << "InvFFT_"<<HistoChan;
    oss6 << "Inverse FFT of the FFT Channel "<<HistoChan<<" histogram; Time (us); ADC";
    std::string Name4 = oss5.str();
    std::string Title4= oss6.str();
    RawDigitInvFFT[HistoChan] = new TH1F(Name4.c_str(), Title4.c_str(), NumBins, 0, NumBins/2);
  }
}

void DAQToOffline::FilterAnalyzer::reconfigure(fhicl::ParameterSet const& pset) {
  fDigitModuleLabel    = pset.get<std::string>("DigitModuleLabel");
  fDigitModuleInstance = pset.get<std::string>("DigitModuleInstance");
  fMakeADCPlots        = pset.get<bool>       ("MakeADCPlots");
  // Make the filter functions.
  std::string colFilt               = pset.get<std::string>("ColFilter");
  std::vector<double> colFiltParams = pset.get<std::vector<double> >("ColFilterParams");
  fColFilterFunc = new TF1("colFilter", colFilt.c_str());
  for(unsigned int i=0; i<colFiltParams.size(); ++i)
    fColFilterFunc->SetParameter(i, colFiltParams[i]);
    
  std::string indUFilt               = pset.get<std::string>("IndUFilter");
  std::vector<double> indUFiltParams = pset.get<std::vector<double> >("IndUFilterParams");
  fIndUFilterFunc = new TF1("indUFilter", indUFilt.c_str());
  for(unsigned int i=0; i<indUFiltParams.size(); ++i)
    fIndUFilterFunc->SetParameter(i, indUFiltParams[i]);

  std::string indVFilt               = pset.get<std::string>("IndVFilter");
  std::vector<double> indVFiltParams = pset.get<std::vector<double> >("IndVFilterParams");
  fIndVFilterFunc = new TF1("indVFilter", indVFilt.c_str());
  for(unsigned int i=0; i<indVFiltParams.size(); ++i)
    fIndVFilterFunc->SetParameter(i, indVFiltParams[i]);
}

void DAQToOffline::FilterAnalyzer::analyze(art::Event const& evt) {

  art::ServiceHandle<geo::Geometry> geo;

  art::Handle<std::vector<raw::RawDigit> > rawDigitHandle;
  evt.getByLabel(fDigitModuleLabel, fDigitModuleInstance, rawDigitHandle);
  std::vector<raw::RawDigit> const& rawDigitVector(*rawDigitHandle);
  
  std::vector< std::pair<int,int> > ZeroFreq;
  //ZeroFreq.push_back( std::make_pair(276 , 285 ) );
  //ZeroFreq.push_back( std::make_pair(558 , 568 ) );
  //ZeroFreq.push_back( std::make_pair(837 , 849 ) );
  //ZeroFreq.push_back( std::make_pair(1116, 1127) );
  //ZeroFreq.push_back( std::make_pair(4340, 5205) );
  
  for (size_t DigLoop=0; DigLoop < rawDigitVector.size(); ++DigLoop) {
    
    int Channel     = rawDigitVector[DigLoop].Channel();
    size_t NADC     = rawDigitVector[DigLoop].NADC();
    double Pedestal = rawDigitVector[DigLoop].GetPedestal();
    const geo::View_t view = geo->View(Channel);
    //std::cout << "Looking at rawDigitVector["<<DigLoop<<"] it was on channel " << rawDigitVector[DigLoop].Channel() << "("<<Channel<<") it is in View " << view
    //        << ", NADC is " << rawDigitVector[DigLoop].NADC() << " ("<<NADC<<")"
    //        << ", pedestal is " << rawDigitVector[DigLoop].GetPedestal() << " ("<<Pedestal<<")"
    //        << std::endl;
    
    // Fill the RawDigit histogram for this histogram.
    for (size_t ADCs=0; ADCs < NADC; ++ADCs) {
      RawDigitHistos[Channel] -> SetBinContent( ADCs+1, rawDigitVector[DigLoop].ADC(ADCs)-Pedestal );
    }
    for (int ww=NADC; ww<16384; ++ww)
      RawDigitHistos[Channel] -> SetBinContent( ww, 0 );
    // Make the FFT for this channel.
    RawDigitHistos[Channel] -> FFT( RawDigitFFT[Channel] ,"MAG");
    for (int bin = 0; bin < RawDigitFFTCorrect[Channel]->GetNbinsX(); ++bin) {
      double BinVal = RawDigitFFT[Channel]->GetBinContent(bin+1);
      RawDigitFFTCorrect[Channel] -> SetBinContent(bin+1, BinVal );
      double freq = 2000. * bin / (double)RawDigitFFTCorrect[Channel]->GetNbinsX();
      if (freq < 1000 && BinVal < 1e5) {
        RawFFT_1000KHz -> Fill( Channel,freq, BinVal);
        if (freq < 100)  {
          RawFFT_100KHz-> Fill( Channel,freq, BinVal);
        }
      }
    }

    // I want to do an inverse FFT, so need to convert the tranformed FFT into an array....
    int NBins = RawDigitFFT[Channel]->GetNbinsX();
    //double Re[NADC], Im[NADC];
    std::unique_ptr<double[]> Re( new double[NADC]);
    std::unique_ptr<double[]> Im( new double[NADC]);
    TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
    fft->GetPointsComplex(Re.get(),Im.get());
    
    // Set the noisy frequency range bins to zero.
    for (size_t aa=0; aa<ZeroFreq.size(); ++aa) {
      for (int bb=ZeroFreq[aa].first; bb<ZeroFreq[aa].second; ++bb) {
        double ReMeanVal=0;
        double ImMeanVal=0;
        int Range = 50;
        for (int cc=0; cc<Range; ++cc) {
          ReMeanVal += Re[ZeroFreq[aa].first-cc] + Re[ZeroFreq[aa].second+cc];
          ImMeanVal += Im[ZeroFreq[aa].first-cc] + Im[ZeroFreq[aa].second+cc];
        }
        ReMeanVal = ReMeanVal / Range;
        Re[bb]    = Re[1500-bb] = ReMeanVal;
        ImMeanVal = ImMeanVal / Range;
        Im[bb]    = Im[1500-bb] = ImMeanVal;
        double MagVal = pow ( Re[bb]*Re[bb] + Im[bb]*Im[bb], 0.5);
        RawDigitFFTCorrect[Channel] -> SetBinContent(bb+1      , MagVal );
        RawDigitFFTCorrect[Channel] -> SetBinContent(15000-bb+1, MagVal );
      }
    }
        
    // Renormalise the X axis, so it is in frequency.
    for (int bin = 0; bin < RawDigitFFTCorrect[Channel]->GetNbinsX(); ++bin){
      RawDigitFFTChannel[Channel] -> SetBinContent( bin+1, RawDigitFFTCorrect[Channel]->GetBinContent(bin+1) );
    }
    // Apply the filter...    
    for (int bin = 0; bin < NBins; ++bin) {
      double freq = 2000. * bin / NBins;
      if (view == geo::kU) { // U plane 
        Re[bin] = Re[bin]*fIndUFilterFunc->Eval(freq);
        Im[bin] = Im[bin]*fIndUFilterFunc->Eval(freq);
      } else if ( view == geo::kV) { // V plane
        Re[bin] = Re[bin]*fIndVFilterFunc->Eval(freq);
        Im[bin] = Im[bin]*fIndVFilterFunc->Eval(freq);
      } else if ( view == geo::kZ) { // Collection plane
        Re[bin] = Re[bin]*fColFilterFunc->Eval(freq);
        Im[bin] = Im[bin]*fColFilterFunc->Eval(freq);
      }
      double MagVal =  pow ( Re[bin]*Re[bin] + Im[bin]*Im[bin], 0.5);
      if (TMath::IsNaN(MagVal)) MagVal = 0;
      RawDigitFFTChannelFilter[Channel] -> SetBinContent( bin+1, MagVal );
            
      // Now do the big histograms...
      if (freq < 1000 && RawDigitFFTChannelFilter[Channel]->GetBinContent(bin+1) < 1e5) {
        FixFFT_1000KHz -> Fill( Channel,freq, MagVal );
        if (freq < 100)  {
          FixFFT_100KHz-> Fill( Channel,freq, MagVal );
        }
      }
    }
    
    // I have applied the filter so now transform back....
    TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &NBins, "C2R");
    fft_back->SetPointsComplex(Re.get(),Im.get());
    fft_back->Transform();
    TH1 *hb=0;
    hb = TH1::TransformHisto(fft_back, hb, "Re");
    for (int BinNum=0; BinNum<NBins; ++BinNum) {
      //if (TMath::IsNaN(hb->GetBinContent(BinNum+1))) {
      ////std::cout << "A bin entry is NaN." << std::endl;
      //RawDigitInvFFT[Channel] -> Fill( BinNum, 0 );
      //} else {
      RawDigitInvFFT[Channel] -> Fill( BinNum, hb->GetBinContent(BinNum+1) / NBins );
      //}
    }
    RawDigitInvFFT[Channel] -> SetXTitle("Time (us)");
    RawDigitInvFFT[Channel] -> SetYTitle("ADC");
  }
  return;
}

DEFINE_ART_MODULE(DAQToOffline::FilterAnalyzer)