EventWaveformDump_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
/// \file    EventWaveformDump_module.cc
/// \module  analyzer
/// \brief   Dump waveforms for a range of events
/// \author  Aidan Medcalf (University of Dallas) <amedcalf@fnal.gov>
/// \date    May 2019
////////////////////////////////////////////////////////////////////////////

// framework
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "art_root_io/TFileDirectory.h"
#include "canvas/Persistency/Common/FindManyP.h"

// nova
#include "RawData/RawBeamline.h"
#include "BeamlineRecoBase/BeamlineDigit.h"
#include "BeamlineRecoBase/BeamlinePID.h"
#include "BeamlineRecoBase/ToF.h"
#include "BeamlineUtils/BeamlineChannelMap.h"
#include "BeamlineUtils/BeamlineCalibration.h"

// ROOT
#include "TH1D.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TFile.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "THStack.h"
#include "TLine.h"
#include "TText.h"

// stl
#include <iostream>

// root

// -----------------------------------------------------------------------
namespace beamlinereco {

  class EventWaveformDump : public art::EDAnalyzer {

  public:

    EventWaveformDump(const fhicl::ParameterSet& pset);

    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob();
    void analyze(const art::Event& e);

  private:
    art::InputTag fRawToFDigitLabel;
    art::InputTag fRawCherenkovDigitLabel;
    art::InputTag fRecoToFDigitLabel;
    art::InputTag fRecoCherenkovDigitLabel;
    art::InputTag fRecoToFLabel;
    std::vector<std::pair<unsigned int, std::string>> fDumpChannels;
    bool fSaveWaveforms;
    bool fDrawMarkers;
    bool fDrawToFs;

    std::map<unsigned int, std::string> dumpChannelMap;

    // "Persistence plot" of all waveforms overlayed, per channel
    std::map<unsigned int, TH2D*> hPersistence;

    // services
    art::ServiceHandle<art::TFileService> tfs;
    art::ServiceHandle<beamlineutil::BeamlineChannelMap> fChannelMap;
    art::ServiceHandle<beamlineutil::BeamlineCalibration> fCalibration;

    unsigned int onlineChannel(ChannelID id);
  };

  DEFINE_ART_MODULE(EventWaveformDump)

}

unsigned int beamlinereco::EventWaveformDump::onlineChannel(ChannelID id) {
  if (id.System == 2) {
    switch (id.Detector) {
    case 0: return 8  + id.Channel;
    case 1: return 12 + id.Channel;
    case 2: return 4  + id.Channel;
    default: return 99999;
    }
  }
  return 99990;
}

// -----------------------------------------------------------------------
beamlinereco::EventWaveformDump::EventWaveformDump(const fhicl::ParameterSet& pset) : EDAnalyzer(pset) {
  this->reconfigure(pset);
}

// -----------------------------------------------------------------------
void beamlinereco::EventWaveformDump::reconfigure(const fhicl::ParameterSet& pset) {
  fRawToFDigitLabel        = pset.get<art::InputTag>("RawToFDigitLabel");
  fRawCherenkovDigitLabel  = pset.get<art::InputTag>("RawCherenkovDigitLabel");
  fRecoToFDigitLabel       = pset.get<art::InputTag>("RecoToFDigitLabel");
  fRecoCherenkovDigitLabel = pset.get<art::InputTag>("RecoCherenkovDigitLabel");
  fRecoToFLabel            = pset.get<art::InputTag>("RecoToFLabel");
  fDumpChannels            = pset.get<std::vector<std::pair<unsigned int, std::string>>>("DumpChannels");
  fSaveWaveforms           = pset.get<bool>("SaveWaveforms");
  fDrawMarkers             = pset.get<bool>("DrawMarkers");
  fDrawToFs                = pset.get<bool>("DrawToFs");

  dumpChannelMap.clear();
  for (auto ch : fDumpChannels) dumpChannelMap[ch.first] = ch.second;
}

// -----------------------------------------------------------------------
void beamlinereco::EventWaveformDump::beginJob() {
  art::TFileDirectory dir  = tfs->mkdir("Persistence");
  for (auto ch : fDumpChannels) {
    hPersistence[ch.first] = dir.make<TH2D>(Form("hPersistence_channel%d",ch.first),
                                            Form("Persistence (Channel %d);Time [ticks];ADC;",ch.first),
                                            1024, 0, 1024,
                                            4096, 0, 4096);
  }
}

// -----------------------------------------------------------------------
void beamlinereco::EventWaveformDump::analyze(const art::Event& e) {
  //std::cout << "+++++++++++++++Run " << e.run() << " Subrun " << e.subRun() << " Event " << e.event() << std::endl;

  std::vector<art::Ptr<rawdata::RawBeamlineDigit>> digits;
  std::vector<art::Ptr<rawdata::RawBeamlineDigit>> tofDigits;
  std::vector<art::Ptr<rawdata::RawBeamlineDigit>> cherenkovDigits;

  std::vector<art::Ptr<brb::BeamlineDigit>> recoDigits;
  std::vector<art::Ptr<brb::BeamlineDigit>> recoToFDigits;
  std::vector<art::Ptr<brb::BeamlineDigit>> recoCherenkovDigits;

  // Get raw tof digits
  art::Handle<std::vector<rawdata::RawBeamlineDigit>> tofDigitHandle;
  if (e.getByLabel(fRawToFDigitLabel, tofDigitHandle))
    art::fill_ptr_vector(tofDigits, tofDigitHandle);
  // Get raw cherenkov digits
  art::Handle<std::vector<rawdata::RawBeamlineDigit>> cherenkovDigitHandle;
  if (e.getByLabel(fRawCherenkovDigitLabel, cherenkovDigitHandle))
    art::fill_ptr_vector(cherenkovDigits, cherenkovDigitHandle);
  // Combine raw vectors
  digits.reserve(tofDigits.size() + cherenkovDigits.size());
  digits.insert(digits.end(), tofDigits.begin(), tofDigits.end());
  digits.insert(digits.end(), cherenkovDigits.begin(), cherenkovDigits.end());

  // Get reco tof digits
  art::Handle<std::vector<brb::BeamlineDigit>> recoToFDigitHandle;
  if (e.getByLabel(fRecoToFDigitLabel, recoToFDigitHandle))
    art::fill_ptr_vector(recoToFDigits, recoToFDigitHandle);
  // Get reco cherenkov digits
  art::Handle<std::vector<brb::BeamlineDigit>> recoCherenkovDigitHandle;
  if (e.getByLabel(fRecoCherenkovDigitLabel, recoCherenkovDigitHandle))
    art::fill_ptr_vector(recoCherenkovDigits, recoCherenkovDigitHandle);
  // Combine reco vectors
  recoDigits.reserve(recoToFDigits.size() + recoCherenkovDigits.size());
  recoDigits.insert(recoDigits.end(), recoToFDigits.begin(), recoToFDigits.end());
  recoDigits.insert(recoDigits.end(), recoCherenkovDigits.begin(), recoCherenkovDigits.end());

  // Get reco tofs
  art::Handle<std::vector<brb::ToF>> tofHandle;
  std::vector<art::Ptr<brb::ToF>>    tofs;
  if (e.getByLabel(fRecoToFLabel, tofHandle))
    art::fill_ptr_vector(tofs, tofHandle);

  // If we have no digits, exit early
  if (digits.size() < 1) return;

  bool haveChannels = false;
  for (auto d : digits) {
    // if channel is in the dump list
    if (dumpChannelMap.count(onlineChannel(d->ChannelID()))) {
      haveChannels = true;
      break;
    }
  }
  if (!haveChannels) return;

  for (auto digit : digits) {
    unsigned int signalchannel = onlineChannel(digit->ChannelID());
    //std::cout << "Run " << e.run() << " Subrun " << e.subRun() << " Event " << e.event() << " OnlineDigitChannel " << signalchannel << "      " << digit->ChannelID() << std::endl;
    if (dumpChannelMap.count(signalchannel) == 0) continue;
    //std::cout << "Channel: " << signalchannel << " about to fill persistence" << std::endl;
    for (unsigned int i=0; i<digit->NADC(); i++) {
      hPersistence[signalchannel]->Fill(i, digit->ADC(i));
    }
    //std::cout << "Filled" << std::endl;
  }

  art::TFileDirectory dir  = tfs->mkdir(Form("Run%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()));
  art::TFileDirectory hdir = dir.mkdir("Histograms");

  TLegend *leg = new TLegend(0.70,0.20,0.88,0.60);

  std::vector<TH1D*> waveforms;
  std::vector<TGraph*> markers;

  //std::cout << "------------------------------"
            //<< "Run " << e.run() << " Subrun " << e.subRun() << " Event " << e.event() << std::endl;
  //for (auto it : dumpChannelMap) { std::cout << "\t\t[" << it.first << ":" << it.second << "]" << std::endl; }
  // For each digit
  for (auto digit : digits) {
    unsigned int signalchannel = onlineChannel(digit->ChannelID());
    //std::cout << signalchannel << " " << dumpChannelMap.count(signalchannel) << std::endl;
    if (dumpChannelMap.count(signalchannel) == 0) continue;

    // Create histogram for waveform
    TH1D *waveform = hdir.make<TH1D>(Form("hr%ds%de%d_ch%d",e.run(),e.subRun(),e.event(),signalchannel), ";Time [ticks];ADC;", digit->NADC(), 0, digit->NADC());
    //std::cout << "channel: " << signalchannel << "\t\t" << digit->ChannelID() << std::endl;
    //std::cout << "\t\t\tDigit has " << digit->NADC() << std::endl;
    //std::cout << signalchannel << "    [";
    // Fill
    for (unsigned int i=0; i<digit->NADC(); i++) {
      //std::cout << digit->ADC(i) << ", ";
      waveform->Fill(i, digit->ADC(i));
    }
    //std::cout << "]" << std::endl;
    
    //std::cout << "Filled waveform" << std::endl;

    // Draw
    waveform->GetYaxis()->SetRangeUser(0,5000);
    waveform->SetLineWidth(2);
    leg->AddEntry(waveform, Form("Channel %d %s", signalchannel, dumpChannelMap[signalchannel].c_str()), "l");
    waveforms.push_back(waveform);

    if (fDrawMarkers) {
      std::vector<art::Ptr<brb::BeamlineDigit>> rds;
      for (auto rd : recoDigits) { if (onlineChannel(rd->ChannelID()) == signalchannel) rds.push_back(rd); }
      if (rds.size() > 0) {
        TGraph *g = new TGraph(rds.size() * 2);
        for (size_t i = 0; i<rds.size(); i++) {
          double starttime = rds[i]->StartTimeInNanoSec() - fCalibration->TimeCorrection(rds[i]->ChannelID());
          double peaktime  = rds[i]->PeakTimeInNanoSec() - fCalibration->TimeCorrection(rds[i]->ChannelID());
          //double starttime = rds[i]->StartTimeInNanoSec();
          //double peaktime  = rds[i]->PeakTimeInNanoSec();
          unsigned int startidx = std::min((unsigned int)1023, (unsigned int)std::max(0.0, (starttime/0.4)));
          unsigned int peakidx  = std::min((unsigned int)1023, (unsigned int)std::max(0.0, (peaktime/0.4)));
          //unsigned int startidx = (unsigned int)std::max(0.0, (starttime/0.4));
          //unsigned int peakidx = (unsigned int)std::max(0.0, (peaktime/0.4));
          //std::cout << "st=" << starttime << " si=" << startidx << " pt=" << peaktime << " pi=" << peakidx << std::endl;
          double startADC = digit->ADC(startidx);
          double peakADC  = digit->ADC(peakidx);
          g->SetPoint(2*i+0, starttime/0.4, startADC);
          g->SetPoint(2*i+1, peaktime/0.4, peakADC);
        }
        g->SetMarkerColor(1);
        g->SetMarkerStyle(3);
        g->SetMarkerSize(2);
        markers.push_back(g);
      }
    }
  }

  // We manually write canvases to the TFileSystem file, which is bad form
  TCanvas *canvas = dir.make<TCanvas>(Form("cRun%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()),  "", 800, 600);
  canvas->cd();

  THStack *stack = new THStack(Form("hsRun%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()), Form("hsRun%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()));
  for (auto w : waveforms) stack->Add(w, "hist");
  stack->Draw("nostack hist PLC");

  if (fDrawMarkers) {
    for (auto g : markers) g->Draw("same p");
  }

  if (fDrawToFs) {
    for (size_t i=0; i<tofs.size(); i++) {
      auto tof = tofs[i];
      TLine *l = new TLine(tof->Timestamps().first/0.4, 4900 - 200*i, tof->Timestamps().second/0.4, 4900 - 200*i);
      l->SetLineColor(1);
      l->SetLineWidth(2);
      l->SetLineStyle(2);
      l->Draw("same");
      double tx = (tof->Timestamps().first + tof->Timestamps().second)/2.0/0.4;
      double ty = 4900 - 200*i + 50;
      TText *t = new TText(tx, ty, Form("tof = %f ns", tof->Time()));
      t->SetTextAlign(21);
      t->SetTextSize(0.02);
      t->Draw();
    }
  }

  leg->Draw("same");
  canvas->Write();

  delete leg;
}