SNBackgroundRejection_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \brief    Remove SN background hits from a file
// \author   Justin Vasel <justin.vasel@gmail.com>
// \date     December 2019
////////////////////////////////////////////////////////////////////////

// STL includes
#include <iostream>

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#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_root_io/TFileService.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"

// NOvASoft includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "NovaTimingUtilities/TimingUtilities.h"
#include "RawData/DAQHeader.h"
#include "RawData/RawTrigger.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Supernova/BackgroundRejection/ColdHotChannelsRemover.h"
#include "Supernova/BackgroundRejection/DaqChannelMask.h"
#include "Supernova/BackgroundRejection/FEBFlasherRemover.h"
#include "Supernova/BackgroundRejection/HitVetoMap.h"
#include "Supernova/BackgroundRejection/LowHighADCRemover.h"
#include "Supernova/BackgroundRejection/SliceRemover.h"
#include "Supernova/BackgroundRejection/TrackActivityRemover.h"
#include "Supernova/BackgroundRejection/HighEnergyRemover.h"
#include "Utilities/HashFunctions.h"
#include "DAQDataFormats/RawDataBlock.h"
#include "DAQDataFormats/RawEvent.h"
#include "RawData/FlatDAQData.h"

// ROOT includes
#include "TTree.h"
#include "TVector3.h"


namespace sn {
  class SNBackgroundRejection;
}

class sn::SNBackgroundRejection : public art::EDProducer {
public:
  explicit SNBackgroundRejection(fhicl::ParameterSet const & p);

  SNBackgroundRejection(SNBackgroundRejection const &) = delete;
  SNBackgroundRejection(SNBackgroundRejection &&) = delete;
  SNBackgroundRejection & operator = (SNBackgroundRejection const &) = delete;
  SNBackgroundRejection & operator = (SNBackgroundRejection &&) = delete;

  float FOM(std::vector<rb::CellHit> hits);
  std::pair<uint64_t,uint64_t> MicrosliceTimeMarkers(rawdata::FlatDAQData& flat);
  void produce(art::Event& e);

  void endJob();


private:
  fhicl::ParameterSet const fParams;

  std::string fRawDataLabel;         ///< ART module label for raw data
  std::string fCellHitLabel;         ///< ART module label for cell hits
  std::string fClusterLabel;         ///< ART module label for clusters
  std::string fCleanupClusterLabel;  ///< ART module label for cleanup clusters
  std::string fTrackLabel;           ///< ART module label for tracks

  bool fVerbose;                     ///< Enable/disable verbose output

  int fNSigFinal;                    ///< Number of signal (MC) hits seen by job
  int fNBkgFinal;                    ///< Number of background hits seen by job

  int fTriggerTimeTotal;
  int fVetoTimeTotal;

  bool fRemoveFirstMicroslice;       ///< Enable/disable removing hits from first microslice

  bool fRemoveLowHighADC;            ///< Enable/disable removing these hits
  int fLowADCThresh;                 ///< Min. acceptable ADC
  int fHighADCThresh;                ///< Max. acceptable ADC

  bool fRemoveColdHotChannels;       ///< Enable/disable removing these hits
  float fColdRateThreshold;          ///< Min. acceptable channel hit rate
  float fHotRateThreshold;           ///< Max. acceptable channel hit rate
  float fRateDuration;               ///< Total time (s) to accumulate statistics

  bool fRemoveFEBFlashers;           ///< Enable/disable removing these hits
  unsigned int fInstigatorQ;         ///< Min. accumulated charge to be a "flash instigator"
  unsigned int fInstigatorT;         ///< Time over which fInstigatorQ charge must be accumulated
  unsigned int fVetoQ;               ///< Min. charge for a hit to be a "flash hit"
  unsigned int fVetoT;               ///< Max. time after the

  bool fRemoveTrackActivity;         ///< Enable/disable removing these hits
  int fMuonEndTimeCut;               ///< Min. time difference between hit and track end to reject [ns]
  int fMuonBodyTimeCut;              ///< Min. time difference between hit and track to reject [ns]
  int fMuonEndSphereRadius;          ///< Min. distance of hit from track end to reject [cm]
  int fMuonBodyCylinderRadius;       ///< Min. distance of hit from track to reject [cm]

  bool fRemovePhysicsSlices;         ///< Enable/disable removing these hits

  bool fRemoveHighEnergyActivity;    ///< Enable/disable removing hits from time windows with HE activity
  float fPeakThreshold;              ///< Number standard deviations above baseline for a peak that should be vetoed
  float fVetoWindowPadding;          ///< Fraction of peak width that should added to either side of TDC veto window

  art::ServiceHandle<geo::Geometry> fGeom;
  sn::DaqChannelMask* fChannelMask;
};


// ............................................................................
sn::SNBackgroundRejection::SNBackgroundRejection(fhicl::ParameterSet const & p):
    EDProducer(p),
    fParams(p),
fRawDataLabel(p.get<std::string>("RawDataLabel")),
fCellHitLabel(p.get<std::string>("CellHitLabel")),
fClusterLabel(p.get<std::string>("ClusterLabel")),
fCleanupClusterLabel(p.get<std::string>("CleanupClusterLabel")),
fTrackLabel(p.get<std::string>("TrackLabel")),
fNSigFinal(0),
fNBkgFinal(0),
fTriggerTimeTotal(0),
fVetoTimeTotal(0),
fRemoveFirstMicroslice(p.get<bool>("RemoveFirstMicroslice")),
fRemoveLowHighADC(p.get<bool>("RemoveLowHighADC")),
fLowADCThresh(p.get<unsigned int>("LowADCThresh")),
fHighADCThresh(p.get<unsigned int>("HighADCThresh")),
fRemoveColdHotChannels(p.get<bool>("RemoveColdHotChannels")),
fColdRateThreshold(p.get<float>("ColdRateThreshold")),
fHotRateThreshold(p.get<float>("HotRateThreshold")),
fRateDuration(p.get<float>("RateDuration")),
fRemoveFEBFlashers(p.get<bool>("RemoveFEBFlashers")),
fInstigatorQ(p.get<unsigned int>("InstigatorQ")),
fInstigatorT(p.get<unsigned int>("InstigatorT")),
fVetoQ(p.get<unsigned int>("VetoQ")),
fVetoT(p.get<unsigned int>("VetoT")),
fRemoveTrackActivity(p.get<bool>("RemoveTrackActivity")),
fMuonEndTimeCut(p.get<int>("MuonEndTimeCut")),
fMuonBodyTimeCut(p.get<int>("MuonBodyTimeCut")),
fMuonEndSphereRadius(p.get<int>("MuonEndSphereRadius")),
fMuonBodyCylinderRadius(p.get<int>("MuonBodyCylinderRadius")),
fRemovePhysicsSlices(p.get<bool>("RemovePhysicsSlices")),
fRemoveHighEnergyActivity(p.get<bool>("RemoveHighEnergyActivity")),
fPeakThreshold(p.get<float>("PeakThreshold")),
fVetoWindowPadding(p.get<float>("VetoWindowPadding")),
fChannelMask(new sn::DaqChannelMask(fColdRateThreshold, fHotRateThreshold))
{
  produces<std::vector<rb::CellHit>>();
  produces<int>();
}


// ............................................................................
float sn::SNBackgroundRejection::FOM(std::vector<rb::CellHit> hits)
{
  int s = 0;
  int b = 0;

  for (rb::CellHit hit : hits) {
    (hit.IsMC()) ? ++s : ++b;
  }

  return (float)s / std::sqrt(s+b);
}

// ............................................................................
std::pair<uint64_t,uint64_t> sn::SNBackgroundRejection::MicrosliceTimeMarkers(rawdata::FlatDAQData& flat)
{
  uint64_t maxTimingMarker = 0x0;
  uint64_t minTimingMarker = 0xFFFFFFFFFFFFFFFF;

  // Unpack raw event
  daqdataformats::RawEvent rawEvt;
  rawEvt.readData(flat.getRawBufferPointer());

  // Loop over data blocks
  size_t nDataBlocks = rawEvt.getDataBlockNumber();
  for (size_t idxDB=0; idxDB<nDataBlocks; ++idxDB) {
    rawEvt.setFloatingDataBlock(idxDB);
    daqdataformats::RawDataBlock dataBlock = (*rawEvt.getFloatingDataBlock());

    // Loop over microslices
    size_t nMicroBlocks = dataBlock.getNumMicroBlocks();
    for (size_t idxMB=0; idxMB<nMicroBlocks; ++idxMB) {
      dataBlock.setFloatingMicroBlock(idxMB);

      // Extract max/min timing marker across all microslices
      uint64_t tm = dataBlock.getFloatingMicroBlock()->getMicroSlice()->getTimingMarker()->getTime();
      if (tm > maxTimingMarker) {
        maxTimingMarker = tm;
      }
      if (tm < minTimingMarker) {
        minTimingMarker = tm;
      }
    }
  }

  return std::make_pair(minTimingMarker, maxTimingMarker);
}


// ............................................................................
void sn::SNBackgroundRejection::produce(art::Event& e)
{
  std::unique_ptr<std::vector<rb::CellHit>> product_cellhits(new std::vector<rb::CellHit>);
  std::unique_ptr<int> product_livetime(new int);

  /* Get the FlatDAQData from the event */
  art::Handle<std::vector<rawdata::FlatDAQData>> hdlFlatDaqData;
  e.getByLabel(fRawDataLabel, hdlFlatDaqData);
  rawdata::FlatDAQData flatDaqData = hdlFlatDaqData->at(0);

  /* Get the trigger information from event */
  art::Handle<std::vector<rawdata::RawTrigger>> hdlTrigger;
  e.getByLabel(fRawDataLabel, hdlTrigger);
  const rawdata::RawTrigger rawTrigger = hdlTrigger->at(0);
  uint64_t triggerStartTime = rawTrigger.fTriggerTimingMarker_TimeStart;
  // NOTE: 1 unit of trigger length = 500 ns = 32 TDC
  uint32_t triggerLengthTDC = rawTrigger.fTriggerRange_TriggerLength * 32;

  /* Get RawDigits from event, put them into a vector, and sort them by time */
  art::Handle<std::vector<rawdata::RawDigit>> hdlDigits;
  e.getByLabel(fRawDataLabel, hdlDigits);
  std::vector<art::Ptr<rawdata::RawDigit>> digits;
  digits.reserve(hdlDigits->size());
  int nHitSigBefore = 0;
  int nHitBkgBefore = 0;
  for (unsigned int idx=0; idx<hdlDigits->size(); ++idx){
    art::Ptr<rawdata::RawDigit> digit(hdlDigits, idx);
    digits.push_back(digit);
    (digit->IsMC()) ? ++nHitSigBefore : ++nHitBkgBefore;
  }

  /* Get CellHits from event, put them into a vector, and sort them by time */
  art::Handle<std::vector<rb::CellHit>> hdlHits;
  e.getByLabel(fCellHitLabel, hdlHits);
  std::vector<art::Ptr<rb::CellHit>> hits;
  hits.reserve(hdlHits->size());
  for (unsigned int idx=0; idx<hdlHits->size(); ++idx){
    art::Ptr<rb::CellHit> hit(hdlHits, idx);

    // Add hit to channel map
    if (fChannelMask->CurrentDuration() < fRateDuration) {
      fChannelMask->AddHit(*hit);
    }

    hits.push_back(hit);
  }
  rb::SortByTime(hits);

  /* Accumulate channel hit stats or calculate rates  */
  if (fChannelMask->CurrentDuration() < fRateDuration) {
    fChannelMask->IncrementDuration((float)rawTrigger.fTriggerRange_TriggerLength * 500e-9);
  } else {
    if (!fChannelMask->RatesCalculated()) fChannelMask->CalculateRates();
    if (fVerbose) std::cout << "Masked " << fChannelMask->MaskSize() << " channels." << std::endl;
    fChannelMask->Print();
  }

  /* Get Clusters from event */
  art::Handle<std::vector<rb::Cluster>> hdlClusters;
  e.getByLabel(fClusterLabel, hdlClusters);

  /* Get Cleanup Clusters from event */
  art::Handle<std::vector<rb::Cluster>> hdlCleanupClusters;
  e.getByLabel(fCleanupClusterLabel, hdlCleanupClusters);
  std::vector<art::Ptr<rb::Cluster>> cleanupClusters;
  art::fill_ptr_vector(cleanupClusters, hdlCleanupClusters);

  /* Get Tracks from event, put them in a ptrVector, and sort them by time */
  art::Handle<std::vector<rb::Track>> hdlTracks;
  e.getByLabel(fTrackLabel, hdlTracks);
  art::FindOneP<rb::Cluster> clusterFromTrack(hdlTracks, e, fTrackLabel);
  std::vector<art::Ptr<rb::Track>> tracks;
  art::fill_ptr_vector(tracks, hdlTracks);
  rb::SortTracksByTime(tracks);


  /* ..... BEGIN BACKGROUND REJECTION ..... */

  // Run each algorithm, add to the hit veto map
  if (fVerbose) std::cout << "=========== SN BACKGROUND REMOVAL ===========" << std::endl << std::endl;

  // Initialize hit veto map & hit removal classes
  sn::HitVetoMap* vetoMapHits = new sn::HitVetoMap();

  /* Remove first microslice */
  if (fRemoveFirstMicroslice) {
    if (fVerbose) std::cout << "-- FIRST (DUPLICATED) MICROSLICE ------------" << std::endl;
    std::pair<uint64_t,uint64_t> microsliceMarkers = this->MicrosliceTimeMarkers(flatDaqData);
    std::cout << "• uslice start time = " << microsliceMarkers.first << std::endl;
    std::cout << "• Remove hits between times: [" << microsliceMarkers.first << ", " << microsliceMarkers.first+3200 << "]" << std::endl;
    for (art::Ptr<rb::CellHit> hit : hits) {
      // Add first 50us to time window veto map (duplicate microslices)
      // NOTE: 3200 64-MHz clock ticks = 50 us
      if (hit->TDC() + triggerStartTime > microsliceMarkers.first + 3200) break;
      vetoMapHits->AddHit(*hit);
    }

    if (fVerbose) std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
  }

  /* Hot/Cold Channel Remove */
  if (fRemoveColdHotChannels) {
    if (fVerbose) std::cout << "-- COLD/HOT CHANNEL HITS --------------------" << std::endl;
    sn::ColdHotChannelsRemover* coldHotChanRemover = new sn::ColdHotChannelsRemover(*vetoMapHits, *fChannelMask, fColdRateThreshold, fHotRateThreshold);
    coldHotChanRemover->remove(digits);
    if (fVerbose) {
      if (!fChannelMask->RatesCalculated()) {
        std::cout << "• Accumulating statistics for cold/hot channel map..." << std::endl;
      } else {
        std::cout << "• Found " << coldHotChanRemover->NumRemovedCold() << " cold channel hits." << std::endl;
        std::cout << "• Found " << coldHotChanRemover->NumRemovedHot()  << " hot channel hits." << std::endl;
      }
      std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
    }
    delete coldHotChanRemover;
  }

  /* FEB Flasher Remove */
  if (fRemoveFEBFlashers) {
    if (fVerbose) std::cout << "-- FEB FLASHERS -----------------------------" << std::endl;
    sn::FEBFlasherRemover* flashRemover = new sn::FEBFlasherRemover(*vetoMapHits, fInstigatorQ, fInstigatorT, fVetoQ, fVetoT);
    flashRemover->remove(digits);
    if (fVerbose) {
      std::cout << "• Found " << flashRemover->NInstigators() << " instigators of " << flashRemover->NFlasherHits() << " flasher hits." << std::endl;
      std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
    }
    delete flashRemover;
  }

  /* Low/High ADC Remove */
  if (fRemoveLowHighADC) {
    if (fVerbose) std::cout << "-- LOW/HIGH-ADC HITS ------------------------" << std::endl;
    sn::LowHighADCRemover* adcRemover = new sn::LowHighADCRemover(*vetoMapHits, fLowADCThresh, fHighADCThresh);
    adcRemover->remove(digits);
    if (fVerbose) {
      std::cout << "• Found " << adcRemover->NumRemovedLow() << " low-ADC hits." << std::endl;
      std::cout << "• Found " << adcRemover->NumRemovedHigh() << " high-ADC hits." << std::endl;
      std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
    }
    delete adcRemover;
  }

  /* Track & Michel Remove */
  if (fRemoveTrackActivity) {
    if (fVerbose) std::cout << "-- TRACK ACTIVITY HITS ----------------------" << std::endl;
    sn::TrackActivityRemover* trackActivityRemover = new sn::TrackActivityRemover(*vetoMapHits, fParams, fGeom, clusterFromTrack);
    trackActivityRemover->SetVetoRegion(fMuonEndTimeCut, fMuonBodyTimeCut, fMuonEndSphereRadius, fMuonBodyCylinderRadius);
    trackActivityRemover->remove(hits, hdlTracks);
    if (fVerbose) {
      std::cout << "• Found " << trackActivityRemover->NumRemovedMichel() << " hits associated with Michel electrons." << std::endl;
      std::cout << "• Found " << trackActivityRemover->NumRemovedBody()   << " hits near track body." << std::endl;
      std::cout << "• Found " << trackActivityRemover->NumRemovedTrack()  << " track hits." << std::endl;
      std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
    }
    delete trackActivityRemover;
  }

  /* Slice Hits Remove */
  if (fRemovePhysicsSlices) {
    if (fVerbose) std::cout << "-- SLICE HITS ----------------------" << std::endl;
    sn::SliceRemover* sliceRemover = new sn::SliceRemover(*vetoMapHits);
    sliceRemover->remove(hits, cleanupClusters);
    if (fVerbose) {
      std::cout << "• Found " << sliceRemover->NumRemoved()
                << " hits associated with " << cleanupClusters.size()
                << " physics clusters." << std::endl;
      std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
    }
    delete sliceRemover;
  }

  /* High-Energy Activity Remove */
  uint32_t TDCVetoed = 0;
  if (fRemoveHighEnergyActivity) {
    if (fVerbose) std::cout << "-- HIGH-ENERGY ACTIVITY ---------------------" << std::endl;
    sn::HighEnergyRemover* highEnergyRemover = new sn::HighEnergyRemover(*vetoMapHits);
    highEnergyRemover->SetPeakThreshold(fPeakThreshold);
    highEnergyRemover->SetVetoWindowPadding(fVetoWindowPadding);
    highEnergyRemover->remove(hits);
    TDCVetoed = highEnergyRemover->TotalTDCVetoed();
    if (fVerbose) {
      std::cout << "• Found " << highEnergyRemover->NPeaks() << " peaks" << std::endl;
      std::cout << "• Vetoed " << highEnergyRemover->TotalTDCVetoed()
                << " / " << triggerLengthTDC << " TDC ("
                << (float)highEnergyRemover->TotalTDCVetoed() / (float)triggerLengthTDC * 100.0
                << "%)" <<  std::endl;
      std::cout << "Veto map size: " << vetoMapHits->Size() << std::endl << std::endl;
    }
    delete highEnergyRemover;
  }

  *product_livetime = triggerLengthTDC - TDCVetoed;
  fTriggerTimeTotal += triggerLengthTDC;
  fVetoTimeTotal += TDCVetoed;

  int nHitSigAfter = 0;
  int nHitBkgAfter = 0;

  // Invert hit veto map
  for (art::Ptr<rb::CellHit> hit : hits) {
    if (!vetoMapHits->ContainsHit(*hit)) {
      product_cellhits->push_back(*hit);
      (hit->IsMC()) ? ++nHitSigAfter : ++nHitBkgAfter;
      (hit->IsMC()) ? ++fNSigFinal : ++fNBkgFinal;
    }
  }

  delete vetoMapHits;

  /* ..... END BACKGROUND REJECTION ..... */


  if (fVerbose) {
    std::cout << "----------------------------------------------" << std::endl;
    std::cout << "Total hits in event: " << hits.size() << std::endl;
    std::cout << "Hits filtered:       " << hits.size() - product_cellhits->size() << std::endl;
    std::cout << "Hits passed:         " << product_cellhits->size() << std::endl;
    std::cout << "Signal loss:         " << nHitSigBefore - nHitSigAfter
              << " (" << (float)(nHitSigBefore-nHitSigAfter)/(float)nHitSigBefore * 100.0 << "%)" << std::endl;
    std::cout << "Bkgrnd loss:         " << nHitBkgBefore - nHitBkgAfter
              << " (" << (float)(nHitBkgBefore-nHitBkgAfter)/(float)nHitBkgBefore * 100.0 << "%)" << std::endl;
    std::cout << "Rejection power:     " << ((float)(hits.size() - product_cellhits->size()) / (float)hits.size()) * 100.0
              << "%" << std::endl << std::endl;
  }

  // Calculate FOM, write value and config params to TTree.
  if (fVerbose) {
    float fom = this->FOM(*product_cellhits);
    std::cout << "FOM = " << fom << std::endl;
    std::cout << "==============================================" << std::endl << std::endl;
  }

  e.put(std::move(product_livetime));
  e.put(std::move(product_cellhits));
  return;
}

void sn::SNBackgroundRejection::endJob()
{
  if (fVerbose) {
    std::cout << std::endl;
    std::cout << "----------------------------------------------"  << std::endl;
    std::cout << "JOB-END SUMMARY"                                 << std::endl;
    std::cout << "----------------------------------------------"  << std::endl;
    std::cout << "Final number of sig hits: " << fNSigFinal        << std::endl;
    std::cout << "Final number of bkg hits: " << fNBkgFinal        << std::endl;
    std::cout << "Final FOM:                " << ((float)fNSigFinal)/std::sqrt(fNBkgFinal+fNSigFinal) << std::endl;
    std::cout << "Final trigger time (TDC): " << fTriggerTimeTotal << std::endl;
    std::cout << "Final veto time (TDC):    " << fVetoTimeTotal    << std::endl;
    std::cout << "Final livetime (ms):      " << (fTriggerTimeTotal-fVetoTimeTotal)*15.625*1e-6 << std::endl;
    std::cout << "Final livetime loss (%):  " << ((float)fVetoTimeTotal / (float)fTriggerTimeTotal)*100 << std::endl;
    std::cout << "----------------------------------------------"  << std::endl;
  }
}

DEFINE_ART_MODULE(sn::SNBackgroundRejection)